TABLE OF CONTENTS
- ABINIT/m_xc_noncoll
- ABINIT/m_xc_noncoll/test_rotations
- m_xc_noncoll/rotate_back_mag
- m_xc_noncoll/rotate_back_mag_dfpt
- m_xc_noncoll/rotate_mag
ABINIT/m_xc_noncoll [ Functions ]
NAME
m_xc_noncoll
FUNCTION
This module provides several routines used in non-collinear XC routines (rotation of the magnetization in order to align it)
COPYRIGHT
Copyright (C) 2001-2018 ABINIT group (EB, MT, FR, SPr) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
OUTPUT
PARENTS
CHILDREN
SOURCE
26 #if defined HAVE_CONFIG_H 27 #include "config.h" 28 #endif 29 30 #include "abi_common.h" 31 32 MODULE m_xc_noncoll 33 34 use defs_basis 35 use m_abicore 36 use m_errors 37 38 implicit none 39 40 private 41 42 ! public procedures 43 public :: rotate_mag ! Rotate a non-collinear density wrt a magnetization 44 public :: rotate_back_mag ! Rotate back a collinear XC potential wrt a magnetization 45 public :: rotate_back_mag_dfpt ! Rotate back a collinear 1st-order XC potential wrt a magnetization 46 public :: test_rotations ! test whether methods in rotate_back_mag_dfpt give similar results 47 48 !Tolerance on magnetization norm 49 real(dp),parameter :: m_norm_min=tol8 50 51 !Default rotation method for DFPT 52 integer,parameter :: rotation_method_default=3 53 54 CONTAINS 55 56 !===========================================================
ABINIT/m_xc_noncoll/test_rotations [ Functions ]
NAME
test_rotations
FUNCTION
Test three different methods in rotate_back_mag_dfpt
INPUTS
option= types of tests to perform
0=> only quick tests
1=> quick and slow tests
cplex = complex or real potential and first order magnetization
OUTPUT
SIDE EFFECTS
NOTES
For debug purposes
PARENTS
CHILDREN
rotate_back_mag_dfpt
SOURCE
907 subroutine test_rotations(option,cplex) 908 909 use defs_basis 910 911 !This section has been created automatically by the script Abilint (TD). 912 !Do not modify the following lines by hand. 913 #undef ABI_FUNC 914 #define ABI_FUNC 'test_rotations' 915 !End of the abilint section 916 917 implicit none 918 919 !Arguments ------------------------------------ 920 integer , intent(in) :: option 921 integer , intent(in) :: cplex 922 923 !Local variables------------------------------- 924 real(dp) :: m0(1,3),vxc0(1,4),kxc(1,3) 925 real(dp) :: n1(cplex,4),vxc1_in(cplex,4),vxc1_out(cplex,4) 926 real(dp) :: delta_23(cplex,4) !,delta_12(cplex,4) 927 real(dp) :: m0_norm,dvdn,dvdz,err23 !,wrong_comp!,err12 928 real(dp) :: theta0,phi0,theta1,phi1,err,m1_norm 929 integer :: dir0,dir1 930 ! ************************************************************************* 931 932 DBG_ENTER("COLL") 933 934 ! if (option/=1 .and. option/=2 ) then 935 ! write(msg,'(3a,i0)')& 936 !& 'The argument option should be 1 or 2,',ch10,& 937 !& 'however, option=',option 938 ! MSG_BUG(msg) 939 ! end if 940 ! 941 ! if (sizein<1) then 942 ! write(msg,'(3a,i0)')& 943 !& ' The argument sizein should be a positive number,',ch10,& 944 !& ' however, sizein=',sizein 945 ! MSG_ERROR(msg) 946 ! end if 947 948 DBG_EXIT("COLL") 949 950 !write(*,*) 'VXC_NONCOLL TESTS================================================================' 951 if (cplex==1) then 952 !write(*,*) ' cplex=1------------------------------------------------------------------------' 953 954 !write(*,*) ' TEST: simple m* orietnations, bxc^(1) part' 955 dvdn=zero;dvdz=1.0! 956 err23=zero 957 do dir0=1,3 958 m0=zero; n1=zero 959 do dir1=2,4 960 m0(1,dir0)=0.1 961 m0_norm=sqrt(m0(1,1)**2+m0(1,2)**2+m0(1,3)**2) 962 n1(1,dir1)=0.8 ! any number would do here 963 964 vxc0=zero; ! no bxc^(0) part at all 965 966 vxc1_in(1,1)= dvdn+dvdz 967 vxc1_in(1,2)= dvdn-dvdz 968 969 call rotate_back_mag_dfpt(1,vxc1_in,vxc1_out,vxc0,kxc,n1,m0,1,1,rot_method=1) 970 call rotate_back_mag_dfpt(1,vxc1_in,vxc1_out,vxc0,kxc,n1,m0,1,1,rot_method=2) 971 delta_23=vxc1_out 972 call rotate_back_mag_dfpt(1,vxc1_in,vxc1_out,vxc0,kxc,n1,m0,1,1,rot_method=3) 973 delta_23=abs(delta_23-vxc1_out) 974 err=max(delta_23(1,1),delta_23(1,2),delta_23(1,3),delta_23(1,4)) 975 if (err23<err) err23=err; 976 977 enddo 978 enddo 979 !write(*,*) ' maximum mismatch between methods 2 and 3:',err23 980 981 !write(*,*) ' TEST: simple m* orietnations, bxc^(0) part' 982 983 err23=zero 984 dvdn=zero;dvdz=1.0 985 do dir0=1,3 986 m0=zero; n1=zero 987 do dir1=2,4 988 m0(1,dir0)=0.1 989 m0_norm=sqrt(m0(1,1)**2+m0(1,2)**2+m0(1,3)**2) 990 n1(1,dir1)=0.8 ! =m^1, any number would do here 991 992 vxc0(1,1) = dvdn+dvdz*m0(1,3)/m0_norm 993 vxc0(1,2) = dvdn-dvdz*m0(1,3)/m0_norm 994 vxc0(1,3) = dvdz*m0(1,1)/m0_norm 995 vxc0(1,4) =-dvdz*m0(1,2)/m0_norm 996 997 vxc1_in=zero !vxc^(1) collinear is zero 998 999 call rotate_back_mag_dfpt(1,vxc1_in,vxc1_out,vxc0,kxc,n1,m0,1,1,rot_method=1) 1000 call rotate_back_mag_dfpt(1,vxc1_in,vxc1_out,vxc0,kxc,n1,m0,1,1,rot_method=2) 1001 delta_23=vxc1_out 1002 call rotate_back_mag_dfpt(1,vxc1_in,vxc1_out,vxc0,kxc,n1,m0,1,1,rot_method=3) 1003 delta_23=abs(delta_23-vxc1_out) 1004 err=maxval(abs(delta_23(1,:))) 1005 if (err23<err) err23=err; 1006 enddo 1007 enddo 1008 !write(*,*) ' maximum mismatch between methods 2 and 3:',err23 1009 1010 !write(*,*) ' TEST: general m0 orietnations, bxc^(0) part' 1011 1012 theta0=zero 1013 err23=zero 1014 m0_norm=0.3 1015 do while(theta0<=pi) 1016 phi0=zero 1017 do while(phi0<=2*pi) 1018 m0(1,1)=m0_norm*sin(theta0)*cos(phi0) 1019 m0(1,2)=m0_norm*sin(theta0)*sin(phi0) 1020 m0(1,3)=m0_norm*cos(theta0) 1021 1022 do dir1=2,4 1023 n1=zero 1024 n1(1,dir1)=0.8 ! =m^1, any number would do here 1025 1026 !vxc0=zero; ! 1027 vxc0(1,1) = dvdn+dvdz*m0(1,3)/m0_norm 1028 vxc0(1,2) = dvdn-dvdz*m0(1,3)/m0_norm 1029 vxc0(1,3) = dvdz*m0(1,1)/m0_norm 1030 vxc0(1,4) =-dvdz*m0(1,2)/m0_norm 1031 1032 vxc1_in=zero 1033 1034 !call rotate_back_mag_dfpt(vxc1_in,vxc1_out,vxc0,kxc,n1,m0,1,1,rot_method=1) 1035 call rotate_back_mag_dfpt(1,vxc1_in,vxc1_out,vxc0,kxc,n1,m0,1,1,rot_method=2) 1036 delta_23=vxc1_out 1037 call rotate_back_mag_dfpt(1,vxc1_in,vxc1_out,vxc0,kxc,n1,m0,1,1,rot_method=3) 1038 delta_23=abs(delta_23-vxc1_out) 1039 err=maxval(abs(delta_23(1,:))) 1040 if (err23<err) err23=err; 1041 enddo 1042 phi0=phi0+2*pi/100.0 1043 enddo 1044 theta0=theta0+pi/100.0 1045 enddo 1046 !write(*,*) ' maximum mismatch between methods 2 and 3:',err23 1047 1048 if(option==2) then 1049 !write(*,*) ' TEST: general m* orietnations, bxc^(0) part' 1050 dvdn=zero;dvdz=1.0 1051 1052 theta0=zero 1053 err23=zero 1054 m0_norm=0.3 1055 m1_norm=10.5 1056 do while(theta0<=pi) !loops on orientation of m^(0) 1057 phi0=zero 1058 do while(phi0<=2*pi) 1059 m0(1,1)=m0_norm*sin(theta0)*cos(phi0) 1060 m0(1,2)=m0_norm*sin(theta0)*sin(phi0) 1061 m0(1,3)=m0_norm*cos(theta0) 1062 1063 vxc0(1,1) = dvdn+dvdz*m0(1,3)/m0_norm 1064 vxc0(1,2) = dvdn-dvdz*m0(1,3)/m0_norm 1065 vxc0(1,3) = dvdz*m0(1,1)/m0_norm 1066 vxc0(1,4) =-dvdz*m0(1,2)/m0_norm 1067 1068 vxc1_in=zero 1069 1070 theta1=zero 1071 do while(theta1<=pi) !loops on orientation of m^(1) 1072 phi1=zero 1073 do while(phi1<=2*pi) 1074 n1(1,1)=zero 1075 n1(1,2)=m1_norm*sin(theta1)*cos(phi1) 1076 n1(1,3)=m1_norm*sin(theta1)*sin(phi1) 1077 n1(1,4)=m1_norm*cos(theta1) 1078 1079 !vxc0=zero; ! 1080 !call rotate_back_mag_dfpt(vxc1_in,vxc1_out,vxc0,kxc,n1,m0,1,1,rot_method=1) 1081 call rotate_back_mag_dfpt(1,vxc1_in,vxc1_out,vxc0,kxc,n1,m0,1,1,rot_method=2) 1082 delta_23=vxc1_out 1083 call rotate_back_mag_dfpt(1,vxc1_in,vxc1_out,vxc0,kxc,n1,m0,1,1,rot_method=3) 1084 delta_23=abs(delta_23-vxc1_out) 1085 err=maxval(abs(delta_23(1,:))) 1086 if (err23<err) err23=err; 1087 phi1=phi1+2*pi/100.0 1088 enddo 1089 theta1=theta1+pi/100.0 1090 enddo 1091 1092 phi0=phi0+2*pi/100.0 1093 enddo 1094 theta0=theta0+pi/100.0 1095 enddo 1096 !write(*,*) ' maximum mismatch between methods 2 and 3:',err23 1097 endif 1098 1099 !else !cplex=2 1100 1101 endif 1102 1103 end subroutine test_rotations
m_xc_noncoll/rotate_back_mag [ Types ]
[ Top ] [ m_xc_noncoll ] [ Types ]
NAME
rotate_back_mag
FUNCTION
Rotate back a collinear XC potential (stored as up+dn) with respect to a magnetization and give a non-collinear XC potential (stored as up_up, dn_dn, Re[up_dn], Im[up_dn]) Note: works only for cplex=1
INPUTS
vxc_in(vectsize,2)=input collinear XC potential mag(vectsize,3)=gs magnetization used for projection vectsize=size of vector fields [mag_norm_in(vectsize)]= --optional-- norm of mag(:) at each point of the grid
OUTPUT
vxc_out(vectsize,4)=output non-collinear XC potential
PARENTS
dfpt_mkvxc_noncoll,m_pawxc
CHILDREN
rotate_back_mag_dfpt
SOURCE
246 subroutine rotate_back_mag(vxc_in,vxc_out,mag,vectsize,& 247 & mag_norm_in) ! optional argument 248 249 250 !This section has been created automatically by the script Abilint (TD). 251 !Do not modify the following lines by hand. 252 #undef ABI_FUNC 253 #define ABI_FUNC 'rotate_back_mag' 254 !End of the abilint section 255 256 implicit none 257 258 !Arguments ------------------------------------ 259 !scalars 260 integer,intent(in) :: vectsize 261 !arrays 262 real(dp),intent(in) :: vxc_in(vectsize,2),mag(vectsize,3) 263 real(dp),intent(out) :: vxc_out(vectsize,4) 264 real(dp),intent(in),optional :: mag_norm_in(vectsize) 265 266 !Local variables------------------------------- 267 !scalars 268 integer :: ipt 269 logical :: has_mag_norm 270 real(dp) :: dvdn,dvdz,m_norm 271 !arrays 272 273 ! ************************************************************************* 274 275 !DBG_ENTER("COLL") 276 277 has_mag_norm=present(mag_norm_in) 278 279 do ipt=1,vectsize 280 if (has_mag_norm) then 281 m_norm=mag_norm_in(ipt) 282 else 283 m_norm=dsqrt(mag(ipt,1)**2+mag(ipt,2)**2+mag(ipt,3)**2) 284 end if 285 286 dvdn=half*(vxc_in(ipt,1)+vxc_in(ipt,2)) 287 288 if (m_norm>m_norm_min) then 289 dvdz=half*(vxc_in(ipt,1)-vxc_in(ipt,2))/m_norm 290 vxc_out(ipt,1)=dvdn+mag(ipt,3)*dvdz 291 vxc_out(ipt,2)=dvdn-mag(ipt,3)*dvdz 292 vxc_out(ipt,3)= mag(ipt,1)*dvdz 293 vxc_out(ipt,4)=-mag(ipt,2)*dvdz 294 else 295 vxc_out(ipt,1:2)=dvdn 296 vxc_out(ipt,3:4)=zero 297 end if 298 end do 299 300 !DBG_EXIT("COLL") 301 302 end subroutine rotate_back_mag
m_xc_noncoll/rotate_back_mag_dfpt [ Types ]
[ Top ] [ m_xc_noncoll ] [ Types ]
NAME
rotate_back_mag_dfpt
FUNCTION
Rotate back a 1st-order collinear XC potential (stored as up+dn) with respect to
a magnetization and give a 1st-order non-collinear XC potential
(stored as up_up, dn_dn, Re{up_dn}, Im{up_dn}).
INPUTS
mag(vectsize,3)=0-order magnetization used for projection
rho1(vectsize,4)=1st-order non-collinear density and magnetization
vxc(vectsize,4)=0-order non-collinear XC potential
kxc(vectsize,nkxc)=0-order XC kernel (associated to vxc)
vxc1_in(vectsize,2)=input 1st-order collinear XC potential
vectsize=size of vector fields
[mag_norm_in(vectsize)]= --optional-- norm of 0-order mag(:) at each point of the grid
[rot_method]=Select method used to compute rotation matrix (1, 2 or 3)
option=if 0, compute only the U^0 vxc^(1) U^0 part
if 1, full first order xc potential
NOTES
cplex=1:
V is stored as : V^11, V^22, Re[V^12], Im[V^12] (complex, hermitian)
N is stored as : n, m_x, m_y, m_z (real)
cplex=2:
V is stored as : V^11, V^22, V^12, i.V^21 (complex)
N is stored as : n, m_x, m_y, mZ (complex)
OUTPUT
vxc1_out(vectsize,4)=output 1st-order non-collinear XC potential
PARENTS
dfpt_mkvxc_noncoll,m_pawxc,m_xc_noncoll
CHILDREN
rotate_back_mag_dfpt
SOURCE
345 subroutine rotate_back_mag_dfpt(option,vxc1_in,vxc1_out,vxc,kxc,rho1,mag,vectsize,cplex,& 346 & mag_norm_in,rot_method) ! optional arguments 347 348 349 !This section has been created automatically by the script Abilint (TD). 350 !Do not modify the following lines by hand. 351 #undef ABI_FUNC 352 #define ABI_FUNC 'rotate_back_mag_dfpt' 353 !End of the abilint section 354 355 implicit none 356 357 !Arguments ------------------------------------ 358 !scalars 359 integer,intent(in) :: vectsize 360 integer,intent(in) :: cplex 361 integer,intent(in) :: option 362 integer,intent(in),optional :: rot_method 363 !arrays 364 real(dp),intent(in) :: kxc(:,:),mag(vectsize,3),vxc(vectsize,4) 365 real(dp),intent(in) :: rho1(cplex*vectsize,4) 366 real(dp),intent(in) :: vxc1_in(cplex*vectsize,2) 367 real(dp),intent(in),optional :: mag_norm_in(vectsize) 368 real(dp),intent(out) :: vxc1_out(cplex*vectsize,4) 369 370 !Local variables------------------------------- 371 !scalars 372 integer :: ipt,rotation_method 373 logical :: has_mag_norm 374 logical :: small_angle 375 real(dp) :: bxc_over_m,d1,d2,d3,d4,dvdn,dvdz,fact,m_dot_m1,m_norm 376 real(dp) :: dvdn_re,dvdn_im,dvdz_re,dvdz_im 377 complex(dpc) :: rho_updn 378 real(dp) :: mdirx,mdiry,mdirz,mxy,mx1,my1,mz1,wx,wy,wx1,wy1 379 real(dp) :: theta0,theta1,theta1_re,theta1_im 380 real(dp) :: wx1_re,wx1_im 381 real(dp) :: wy1_re,wy1_im 382 real(dp) :: mx1_re,mx1_im,my1_re,my1_im,mz1_re,mz1_im 383 real(dp) :: m_dot_m1_re,m_dot_m1_im 384 real(dp) :: bxc 385 !arrays 386 real(dp) :: vxc_diag(2),v21tmp(2) 387 complex(dpc) :: r1tmp(2,2),u0(2,2),u0_1(2,2),u0_1r1(2,2),u0v1(2,2) 388 complex(dpc) :: rho1_updn(2,2),v1tmp(2,2),vxc1tmp(2,2) 389 complex(dpc) :: rho1_offdiag(2) 390 391 ! ************************************************************************* 392 393 !DBG_ENTER("COLL") 394 395 !Optional arguments 396 has_mag_norm=present(mag_norm_in) 397 rotation_method=rotation_method_default 398 if (present(rot_method)) rotation_method=rot_method 399 400 !Check Kxc 401 if (size(kxc)>3*vectsize) then 402 MSG_ERROR('Cannot use Kxc from GGA!') 403 end if 404 405 if((rotation_method==1.or.rotation_method==2).and.cplex==2) then 406 MSG_ERROR('rotation_method=1 and 2 are not available for cplex=2 case! use ixcrot=3') 407 endif 408 409 410 select case (rotation_method) 411 412 !---------------------------------------- 413 ! Taylor expansion of U rotation matrix 414 !---------------------------------------- 415 case (1) 416 417 do ipt=1,vectsize 418 419 if (has_mag_norm) then 420 m_norm=mag_norm_in(ipt) 421 else 422 m_norm=sqrt(mag(ipt,1)**2+mag(ipt,2)**2+mag(ipt,3)**2) 423 end if 424 425 if(m_norm>m_norm_min) then 426 427 ! Define the U^(0) transformation matrix 428 rho_updn=(mag(ipt,1)+(zero,one)*mag(ipt,2)) 429 d1=sqrt(( m_norm+mag(ipt,3))**2+rho_updn**2) 430 d2=sqrt((-m_norm+mag(ipt,3))**2+rho_updn**2) 431 d3=sqrt(( m_norm-mag(ipt,3))**2+rho_updn**2) 432 d4=sqrt(( m_norm+mag(ipt,3))**2-rho_updn**2) 433 u0(1,1)=( m_norm+mag(ipt,3))/d1 ! ( m + mz)/d1 434 u0(2,2)=rho_updn/d2 ! ( mx +imy)/d2 435 u0(1,2)=(-m_norm+mag(ipt,3))/d2 ! (-m + mz)/d2 436 u0(2,1)=rho_updn/d1 ! ( mx +imy)/d1 437 438 ! Define the inverse of U^(0): U^(0)^-1 439 u0_1(1,1)= half*d1/m_norm 440 u0_1(2,2)= half*d2*(m_norm+mag(ipt,3))/(m_norm*rho_updn) 441 u0_1(1,2)= half*d1*(m_norm-mag(ipt,3))/(m_norm*rho_updn) 442 u0_1(2,1)=-half*d2/m_norm 443 444 ! Diagonalize the GS Vxc^(0): U^(0)^-1 Vxc^(0) U^(0) 445 ! (Remember the abinit notation for vxc!) 446 vxc_diag(1)=half*(vxc(ipt,1)+vxc(ipt,2) & 447 & -sqrt((vxc(ipt,1)-vxc(ipt,2))**2 & 448 & +four*(vxc(ipt,3)**2+vxc(ipt,4)**2))) 449 vxc_diag(2)=half*(vxc(ipt,1)+vxc(ipt,2) & 450 & +sqrt((vxc(ipt,1)-vxc(ipt,2))**2 & 451 & +four*(vxc(ipt,3)**2+vxc(ipt,4)**2))) 452 v1tmp(1,1)=cmplx(real(vxc1_in(ipt,1),kind=dp),zero) 453 v1tmp(2,2)=cmplx(real(vxc1_in(ipt,2),kind=dp),zero) 454 455 !Tranforming the rhor1 with U0 456 rho1_updn(1,1)=half*(rho1(ipt,1)+rho1(ipt,4)) 457 rho1_updn(2,2)=half*(rho1(ipt,1)-rho1(ipt,4)) 458 rho1_updn(1,2)=half*(rho1(ipt,2)-(zero,one)*rho1(ipt,3)) 459 rho1_updn(2,1)=half*(rho1(ipt,2)+(zero,one)*rho1(ipt,3)) 460 u0_1r1=matmul(u0_1,rho1_updn) 461 r1tmp=matmul(u0_1r1,u0) 462 rho1_offdiag(1)=r1tmp(1,2) ; rho1_offdiag(2)=r1tmp(2,1) 463 464 if (option==0) then ! for xccc alone 465 v1tmp(1,2)=cmplx(zero,zero) 466 v1tmp(2,1)=cmplx(zero,zero) 467 else 468 v1tmp(1,2)=-(rho1_offdiag(1)/m_norm)*(vxc_diag(2)-vxc_diag(1)) 469 v1tmp(2,1)= (rho1_offdiag(2)/m_norm)*(vxc_diag(1)-vxc_diag(2)) 470 endif 471 472 !Rotate back the "diagonal" xc computing the term U^(0) Vxc1_^(1) U^(0)^-1 473 u0v1=matmul(u0,v1tmp) 474 vxc1tmp=matmul(u0v1,u0_1) 475 vxc1_out(ipt,1)=real(vxc1tmp(1,1),kind=dp) 476 vxc1_out(ipt,2)=real(vxc1tmp(2,2),kind=dp) 477 vxc1_out(ipt,3)=real( real(vxc1tmp(1,2)),kind=dp) 478 vxc1_out(ipt,4)=real(aimag(vxc1tmp(1,2)),kind=dp) 479 480 else ! Magnetization is zero 481 dvdn=(vxc1_in(ipt,1)+vxc1_in(ipt,2))*half 482 mx1=rho1(ipt,2) ; my1=rho1(ipt,3) ; mz1=rho1(ipt,4) 483 ! Compute Bxc/|m| from Kxc (zero limit) 484 bxc_over_m = half*(half*(kxc(ipt,1)+kxc(ipt,3))-kxc(ipt,2)) 485 vxc1_out(ipt,1)= dvdn + bxc_over_m*mz1 486 vxc1_out(ipt,2)= dvdn - bxc_over_m*mz1 487 vxc1_out(ipt,3)= bxc_over_m*mx1 488 vxc1_out(ipt,4)=-bxc_over_m*my1 489 end if 490 491 end do ! ipt 492 493 !---------------------------------------- 494 ! Analytical expression of U rotation matrix 495 !---------------------------------------- 496 case (2) 497 !Alternative method (explicitely calculated rotation matrices) 498 !Vxc^(1) = phixc^(1).Id + // <= change of "electrostatic" XC potential (phixc^(1) is denoted dvdn) 499 ! + bxc^(1)*( Udag^(0).sigma_z.U^(0) ) + // <= this part describes the change of XC magnetic field magnitude bxc^(1) 500 ! + bxc^(0)*( Udag^(1).sigma_z.U^(0) + Udag^(0).sigma_z.U^(1) ) // <= remaining terms describe the cost of magnetization rotation 501 502 select case(cplex) 503 case(1) 504 do ipt=1,vectsize 505 506 if (has_mag_norm) then 507 m_norm=mag_norm_in(ipt) 508 else 509 m_norm=dsqrt(mag(ipt,1)**2+mag(ipt,2)**2+mag(ipt,3)**2) 510 end if 511 512 513 mx1 =rho1(ipt,2); 514 my1 =rho1(ipt,3); 515 mz1 =rho1(ipt,4) 516 517 dvdn=(vxc1_in(ipt,1)+vxc1_in(ipt,2))*half !phixc^(1) 518 dvdz=(vxc1_in(ipt,1)-vxc1_in(ipt,2))*half !bxc^(1) 519 520 if (m_norm>m_norm_min) then 521 522 mxy = dsqrt(mag(ipt,1)**2+mag(ipt,2)**2) 523 small_angle=(mxy/m_norm<tol8) !condition for sin(x)~x to be valid 524 !even possible to set to tol6 525 mdirx=mag(ipt,1)/m_norm 526 mdiry=mag(ipt,2)/m_norm 527 mdirz=mag(ipt,3)/m_norm 528 529 ! dvdn is phixc^(1) (density only part) 530 ! dvdz is bxc^(1) (magnetization magnitude part) 531 532 !U^(0)*.Vxc1.U^(0) part 533 vxc1_out(ipt,1)= dvdn+dvdz*mdirz 534 vxc1_out(ipt,2)= dvdn-dvdz*mdirz 535 vxc1_out(ipt,3)= dvdz*mdirx ! Real part 536 vxc1_out(ipt,4)=-dvdz*mdiry ! Imaginary part, minus sign comes from sigma_y 537 538 !U^(1)*.Vxc0.U^(0) + U^(0)*.Vxc0.U^(1) part 539 540 !bxc = dsqrt(((vxc(ipt,1)-vxc(ipt,2))*half)**2+vxc(ipt,3)**2+vxc(ipt,4)**2) !bxc^(0) 541 bxc = (vxc(ipt,1)-vxc(ipt,2))*half/mag(ipt,3)*m_norm 542 if (.not.small_angle) then 543 wx = mag(ipt,2)/mxy 544 wy =-mag(ipt,1)/mxy 545 theta0 = dacos(mag(ipt,3)/m_norm) 546 547 theta1 = (mdirz*(mdirx*mx1+mdiry*my1))/mxy - mz1*mxy/m_norm**2 548 wx1 = (+mag(ipt,1)**2*my1 - mag(ipt,1)*mag(ipt,2)*mx1)/mxy**2/m_norm ! wx1 multiplied by sin(theta)=mxy/m_norm 549 wy1 = (-mag(ipt,2)**2*mx1 + mag(ipt,1)*mag(ipt,2)*my1)/mxy**2/m_norm ! wx1 multiplied by sin(theta)=mxy/m_norm 550 551 vxc1_out(ipt,1) = vxc1_out(ipt,1) - bxc*dsin(theta0)*theta1 552 vxc1_out(ipt,2) = vxc1_out(ipt,2) + bxc*dsin(theta0)*theta1 553 vxc1_out(ipt,3) = vxc1_out(ipt,3) - bxc*(wy1+dcos(theta0)*wy*theta1) 554 vxc1_out(ipt,4) = vxc1_out(ipt,4) - bxc*(wx1+dcos(theta0)*wx*theta1) 555 else 556 !zero order terms O(1) 557 vxc1_out(ipt,3) = vxc1_out(ipt,3) + bxc*mx1/abs(mag(ipt,3)) 558 vxc1_out(ipt,4) = vxc1_out(ipt,4) - bxc*my1/abs(mag(ipt,3)) 559 !first order terms O(theta) 560 fact = bxc/(mag(ipt,3)*abs(mag(ipt,3))) 561 vxc1_out(ipt,1) = vxc1_out(ipt,1) - (mag(ipt,1)*mx1+mag(ipt,2)*my1)*fact 562 vxc1_out(ipt,2) = vxc1_out(ipt,2) + (mag(ipt,1)*mx1+mag(ipt,2)*my1)*fact 563 vxc1_out(ipt,3) = vxc1_out(ipt,3) - mag(ipt,1)*mz1*fact 564 vxc1_out(ipt,4) = vxc1_out(ipt,4) + mag(ipt,2)*mz1*fact 565 endif 566 567 else ! Magnetization is zero 568 ! Compute Bxc/|m| from Kxc (zero limit) 569 bxc_over_m = half*(half*(kxc(ipt,1)+kxc(ipt,3))-kxc(ipt,2)) 570 vxc1_out(ipt,1)= dvdn + bxc_over_m*mz1 571 vxc1_out(ipt,2)= dvdn - bxc_over_m*mz1 572 vxc1_out(ipt,3)= bxc_over_m*mx1 573 vxc1_out(ipt,4)=-bxc_over_m*my1 574 end if 575 end do ! ipt 576 577 case(2) !cplex=2 578 579 do ipt=1,vectsize 580 581 if (has_mag_norm) then 582 m_norm=mag_norm_in(ipt) 583 else 584 m_norm=dsqrt(mag(ipt,1)**2+mag(ipt,2)**2+mag(ipt,3)**2) 585 end if 586 587 mx1_re= rho1(2*ipt-1,2); mx1_im= rho1(2*ipt,2) 588 my1_re= rho1(2*ipt-1,3); my1_im= rho1(2*ipt,3) 589 mz1_re= rho1(2*ipt-1,4); mz1_im= rho1(2*ipt,4) 590 591 dvdn_re=(vxc1_in(2*ipt-1,1)+vxc1_in(2*ipt-1,2))*half 592 dvdz_re=(vxc1_in(2*ipt-1,1)-vxc1_in(2*ipt-1,2))*half 593 dvdn_im=(vxc1_in(2*ipt ,1)+vxc1_in(2*ipt ,2))*half 594 dvdz_im=(vxc1_in(2*ipt ,1)-vxc1_in(2*ipt ,2))*half 595 596 if (m_norm>m_norm_min) then 597 598 mdirx=mag(ipt,1)/m_norm 599 mdiry=mag(ipt,2)/m_norm 600 mdirz=mag(ipt,3)/m_norm 601 602 mxy = dsqrt(mag(ipt,1)**2+mag(ipt,2)**2) 603 small_angle=(mxy/m_norm<tol8) !condition for sin(x)~x to be valid 604 ! 605 mdirx=mag(ipt,1)/m_norm 606 mdiry=mag(ipt,2)/m_norm 607 mdirz=mag(ipt,3)/m_norm 608 609 ! dvdn is phixc^(1) (density only part) 610 ! dvdz is bxc^(1) (magnetization magnitude part) 611 612 !U^(0)*.Vxc1.U^(0) part 613 vxc1_out(2*ipt-1,1)= dvdn_re+dvdz_re*mdirz 614 vxc1_out(2*ipt ,1)= dvdn_im+dvdz_im*mdirz 615 vxc1_out(2*ipt-1,2)= dvdn_re-dvdz_re*mdirz 616 vxc1_out(2*ipt ,2)= dvdn_im-dvdz_im*mdirz 617 !NOTE: change of definition of the potential matrix components 618 ! vxc1_out(:,3) = V_updn 619 ! vxc1_out(:,4) = i.V_updn 620 vxc1_out(2*ipt-1,3)= dvdz_re*mdirx + dvdz_im*mdiry !Re[ V^12] 621 vxc1_out(2*ipt ,3)= dvdz_im*mdirx - dvdz_re*mdiry !Im[ V^12] 622 623 !U^(1)*.Vxc0.U^(0) + U^(0)*.Vxc0.U^(1) part 624 !bxc = dsqrt(((vxc(ipt,1)-vxc(ipt,2))*half)**2+vxc(ipt,3)**2+vxc(ipt,4)**2) !bxc^(0) 625 bxc = (vxc(ipt,1)-vxc(ipt,2))*half/mag(ipt,3)*m_norm 626 if (.not.small_angle) then 627 wx = mag(ipt,2)/mxy 628 wy =-mag(ipt,1)/mxy 629 theta0 = dacos(mag(ipt,3)/m_norm) 630 631 632 theta1_re = (mdirz*(mdirx*mx1_re+mdiry*my1_re))/mxy - mz1_re*mxy/m_norm**2 633 theta1_im = (mdirz*(mdirx*mx1_im+mdiry*my1_im))/mxy - mz1_im*mxy/m_norm**2 634 635 wx1_re = (+mag(ipt,1)**2*my1_re - mag(ipt,1)*mag(ipt,2)*mx1_re)/mxy**2/m_norm ! wx1 multiplied by sin(theta)=mxy/m_norm 636 wx1_im = (+mag(ipt,1)**2*my1_im - mag(ipt,1)*mag(ipt,2)*mx1_im)/mxy**2/m_norm 637 wy1_re = (-mag(ipt,2)**2*mx1_re + mag(ipt,1)*mag(ipt,2)*my1_re)/mxy**2/m_norm ! wy1 multiplied by sin(theta)=mxy/m_norm 638 wy1_im = (-mag(ipt,2)**2*mx1_im + mag(ipt,1)*mag(ipt,2)*my1_im)/mxy**2/m_norm 639 640 !U^(1)*.Vxc0.U^(0) + U^(0)*.Vxc0.U^(1) 641 vxc1_out(2*ipt-1,1) = vxc1_out(2*ipt-1,1) - bxc*dsin(theta0)*theta1_re 642 vxc1_out(2*ipt ,1) = vxc1_out(2*ipt ,1) - bxc*dsin(theta0)*theta1_im 643 vxc1_out(2*ipt-1,2) = vxc1_out(2*ipt-1,2) + bxc*dsin(theta0)*theta1_re 644 vxc1_out(2*ipt ,2) = vxc1_out(2*ipt ,2) + bxc*dsin(theta0)*theta1_im 645 !cplex=1 part: 646 !v12 += -(bxc)*(wy1+dcos(theta0)*wy*theta1)- 647 ! -i.(bxc)*(wx1+dcos(theta0)*wx*theta1) 648 vxc1_out(2*ipt-1,3) = vxc1_out(2*ipt-1,3) - bxc*(wy1_re+dcos(theta0)*wy*theta1_re) 649 vxc1_out(2*ipt-1,3) = vxc1_out(2*ipt-1,3) + bxc*(wx1_im+dcos(theta0)*wx*theta1_im) 650 vxc1_out(2*ipt ,3) = vxc1_out(2*ipt ,3) - bxc*(wy1_im+dcos(theta0)*wy*theta1_im) 651 vxc1_out(2*ipt ,3) = vxc1_out(2*ipt ,3) - bxc*(wx1_re+dcos(theta0)*wx*theta1_re) 652 else 653 !small theta case: 654 !zero order terms O(1) 655 vxc1_out(2*ipt-1,3) = vxc1_out(2*ipt-1,3) + bxc*mx1_re/abs(mag(ipt,3)) 656 vxc1_out(2*ipt-1,3) = vxc1_out(2*ipt-1,3) + bxc*my1_im/abs(mag(ipt,3)) 657 vxc1_out(2*ipt ,3) = vxc1_out(2*ipt ,3) + bxc*mx1_im/abs(mag(ipt,3)) 658 vxc1_out(2*ipt ,3) = vxc1_out(2*ipt ,3) - bxc*my1_re/abs(mag(ipt,3)) 659 !first order terms: 660 fact = bxc/(mag(ipt,3)*abs(mag(ipt,3))) 661 vxc1_out(2*ipt-1,1) = vxc1_out(2*ipt-1,1) - (mag(ipt,1)*mx1_re+mag(ipt,2)*my1_re)*fact 662 vxc1_out(2*ipt-1,2) = vxc1_out(2*ipt-1,2) + (mag(ipt,1)*mx1_re+mag(ipt,2)*my1_re)*fact 663 vxc1_out(2*ipt ,1) = vxc1_out(2*ipt ,1) - (mag(ipt,1)*mx1_im+mag(ipt,2)*my1_im)*fact 664 vxc1_out(2*ipt ,2) = vxc1_out(2*ipt ,2) + (mag(ipt,1)*mx1_im+mag(ipt,2)*my1_im)*fact 665 666 vxc1_out(2*ipt-1,3) = vxc1_out(2*ipt-1,3) - mag(ipt,1)*mz1_re*fact 667 vxc1_out(2*ipt-1,3) = vxc1_out(2*ipt-1,3) - mag(ipt,2)*mz1_im*fact 668 vxc1_out(2*ipt ,3) = vxc1_out(2*ipt ,3) + mag(ipt,2)*mz1_re*fact 669 vxc1_out(2*ipt ,3) = vxc1_out(2*ipt ,3) - mag(ipt,1)*mz1_im*fact 670 endif 671 672 else ! Magnetization is practically zero 673 ! Compute Bxc/|m| from Kxc (zero limit) 674 bxc_over_m = half*(half*(kxc(ipt,1)+kxc(ipt,3))-kxc(ipt,2)) 675 vxc1_out(2*ipt-1,1)= dvdn_re + bxc_over_m*mz1_re 676 vxc1_out(2*ipt ,1)= dvdn_im + bxc_over_m*mz1_im 677 vxc1_out(2*ipt-1,2)= dvdn_re - bxc_over_m*mz1_re 678 vxc1_out(2*ipt ,2)= dvdn_im - bxc_over_m*mz1_im 679 vxc1_out(2*ipt-1,3)= bxc_over_m*( mx1_re+my1_im) 680 vxc1_out(2*ipt ,3)= bxc_over_m*(-my1_re+mx1_im) 681 end if 682 !finally reconstruct i.V^12 from V^12 683 vxc1_out(2*ipt-1,4) = vxc1_out(2*ipt ,3) ! Re[i.V^21] = Im[V^12] 684 vxc1_out(2*ipt ,4) = vxc1_out(2*ipt-1,3) ! Im[i.V^21] = Re[V^12] 685 686 end do ! ipt 687 688 end select 689 690 !---------------------------------------- 691 ! Explicit derivative of the rotated XC functional 692 !---------------------------------------- 693 case (3) 694 ! Brute-force derivative of Vxc 695 ! Explicit calculation of the rotated xc functional 696 ! (derivatives of the analytical expression) (Eq. A) 697 ! Vxc^(1) = phixc^(1).Id + // <= change of "electrostatic" XC potential (phixc^(1) is denoted dvdn) 698 ! + bxc^(1)*(sigma,m^(0))/|m^(0)| + // <= this term is equivalent to ( Udag^(0).sigma_z.U^(0) ) term in rotation_method=2 699 ! + bxc^(0)*(sigma,m^(1)))/|m^(0)| - // <= the last terms are equivalent to ( Udag^(1).sigma_z.U^(0) + Udag^(0).sigma_z.U^(1) ) 700 ! - bxc^(0)*(sigma,m^(0))*(m^(1),m^(0))/|m^(0)|**3 701 select case(cplex) 702 case(1) 703 704 do ipt=1,vectsize 705 706 if (has_mag_norm) then 707 m_norm=mag_norm_in(ipt) 708 else 709 m_norm=sqrt(mag(ipt,1)**2+mag(ipt,2)**2+mag(ipt,3)**2) 710 end if 711 712 ! dvdn is phixc^(1) (density only part) 713 ! dvdz is bxc^(1) (magnetization magnitude part) 714 dvdn=(vxc1_in(ipt,1)+vxc1_in(ipt,2))*half 715 dvdz=(vxc1_in(ipt,1)-vxc1_in(ipt,2))*half 716 717 mx1=rho1(ipt,2) ; my1=rho1(ipt,3) ; mz1=rho1(ipt,4) 718 719 if(m_norm>m_norm_min) then 720 721 mdirx=mag(ipt,1)/m_norm; mdiry=mag(ipt,2)/m_norm; mdirz=mag(ipt,3)/m_norm 722 723 !This part describes the change of the magnitude of the xc magnetic field 724 !and the change of the scalar part of the xc electrostatic potential, 1st + 2nd term in Eq.A 725 !phixc^(1).Id + bxc^(1) (sigma,m^(0))/|m^(0)| 726 vxc1_out(ipt,1)= dvdn+dvdz*mdirz 727 vxc1_out(ipt,2)= dvdn-dvdz*mdirz 728 vxc1_out(ipt,3)= dvdz*mdirx ! Real part 729 vxc1_out(ipt,4)=-dvdz*mdiry ! Imaginary part, minus sign comes from sigma_y 730 731 if (option/=0) then 732 !Add remaining contributions comming from the change of magnetization direction 733 !projection of m^(1) on gs magnetization direction 734 m_dot_m1=(mdirx*rho1(ipt,2)+mdiry*rho1(ipt,3)+mdirz*rho1(ipt,4)) 735 736 bxc_over_m =-dsqrt(((vxc(ipt,1)-vxc(ipt,2))*half)**2+vxc(ipt,3)**2+vxc(ipt,4)**2) !this is bxc^(0) 737 bxc_over_m = bxc_over_m/m_norm 738 vxc1_out(ipt,1) = vxc1_out(ipt,1) + bxc_over_m*( mz1 - mdirz*m_dot_m1 ) ! 739 vxc1_out(ipt,2) = vxc1_out(ipt,2) + bxc_over_m*(-mz1 + mdirz*m_dot_m1 ) ! 740 vxc1_out(ipt,3) = vxc1_out(ipt,3) + bxc_over_m*( mx1 - mdirx*m_dot_m1 ) ! 741 vxc1_out(ipt,4) = vxc1_out(ipt,4) + bxc_over_m*(-my1 + mdiry*m_dot_m1 ) ! 742 endif 743 744 else 745 if (option/=0) then 746 !Compute bxc^(0)/|m| from kxc (|m^(0)| -> zero limit) 747 bxc_over_m = half*(half*(kxc(ipt,1)+kxc(ipt,3))-kxc(ipt,2)) 748 vxc1_out(ipt,1)= dvdn + bxc_over_m*mz1 749 vxc1_out(ipt,2)= dvdn - bxc_over_m*mz1 750 vxc1_out(ipt,3)= bxc_over_m*mx1 751 vxc1_out(ipt,4)=-bxc_over_m*my1 752 else 753 vxc1_out(ipt,1)= dvdn 754 vxc1_out(ipt,2)= dvdn 755 vxc1_out(ipt,3)= zero 756 vxc1_out(ipt,4)= zero 757 endif 758 end if 759 760 end do ! ipt 761 762 case(2) 763 !cplex=2 case 764 765 do ipt=1,vectsize 766 767 if (has_mag_norm) then 768 m_norm=mag_norm_in(ipt) 769 else 770 m_norm=sqrt(mag(ipt,1)**2+mag(ipt,2)**2+mag(ipt,3)**2) 771 end if 772 773 ! see cplex=1 case for details 774 dvdn_re=(vxc1_in(2*ipt-1,1)+vxc1_in(2*ipt-1,2))*half 775 dvdn_im=(vxc1_in(2*ipt ,1)+vxc1_in(2*ipt ,2))*half 776 dvdz_re=(vxc1_in(2*ipt-1,1)-vxc1_in(2*ipt-1,2))*half 777 dvdz_im=(vxc1_in(2*ipt ,1)-vxc1_in(2*ipt ,2))*half 778 779 mx1_re=rho1(2*ipt-1,2); mx1_im=rho1(2*ipt,2) 780 my1_re=rho1(2*ipt-1,3); my1_im=rho1(2*ipt,3) 781 mz1_re=rho1(2*ipt-1,4); mz1_im=rho1(2*ipt,4) 782 783 if(m_norm>m_norm_min) then 784 785 mdirx=mag(ipt,1)/m_norm; mdiry=mag(ipt,2)/m_norm; mdirz=mag(ipt,3)/m_norm 786 787 !first two terms: 788 vxc1_out(2*ipt-1,1)= dvdn_re+dvdz_re*mdirz 789 vxc1_out(2*ipt ,1)= dvdn_im+dvdz_im*mdirz 790 vxc1_out(2*ipt-1,2)= dvdn_re-dvdz_re*mdirz 791 vxc1_out(2*ipt ,2)= dvdn_im-dvdz_im*mdirz 792 !NOTE: change of definition of the potential matrix components 793 ! vxc1_out(:,3) = V_updn 794 ! vxc1_out(:,4) = i.V_dnup 795 796 ! V^12 = dvdz*mx/|m| - i.dvdz*my/|m| = (Re[dvdz]*mx/|m| + Im[dvdz]*my/|m|) + i.(Im[dvdz]*mx/|m| - Re[dvdz]*my/|m|) => vxc1(:,3) 797 ! V^21 = dvdz*mx/|m| + i.dvdz*my/|m| = (Re[dvdz]*mx/|m| - Im[dvdz]*my/|m|) + i.(Im[dvdz]*mx/|m| + Re[dvdz]*my/|m|) 798 vxc1_out(2*ipt-1,3)= dvdz_re*mdirx + dvdz_im*mdiry !Re[V^12] 799 vxc1_out(2*ipt ,3)= dvdz_im*mdirx - dvdz_re*mdiry !Im[V^12] 800 vxc1_out(2*ipt-1,4)= dvdz_re*mdirx - dvdz_im*mdiry !Re[V^21] 801 vxc1_out(2*ipt ,4)= dvdz_im*mdirx + dvdz_re*mdiry !Im[V^21] 802 if (option/=0) then 803 804 !remaining contributions: 805 m_dot_m1_re= mdirx*mx1_re + mdiry*my1_re + mdirz*mz1_re 806 m_dot_m1_im= mdirx*mx1_im + mdiry*my1_im + mdirz*mz1_im 807 808 bxc_over_m =-dsqrt(((vxc(ipt,1)-vxc(ipt,2))*half)**2+vxc(ipt,3)**2+vxc(ipt,4)**2) !this is bxc^(0) 809 bxc_over_m = bxc_over_m/m_norm 810 !bxc_over_m = (vxc(ipt,1)-vxc(ipt,2))*half/mag(ipt,3) 811 812 vxc1_out(2*ipt-1,1) = vxc1_out(2*ipt-1,1) + bxc_over_m*( mz1_re - mdirz*m_dot_m1_re ) ! Re[V^11] 813 vxc1_out(2*ipt ,1) = vxc1_out(2*ipt ,1) + bxc_over_m*( mz1_im - mdirz*m_dot_m1_im ) ! Im[V^11] 814 vxc1_out(2*ipt-1,2) = vxc1_out(2*ipt-1,2) + bxc_over_m*(-mz1_re + mdirz*m_dot_m1_re ) ! Re[V^22] 815 vxc1_out(2*ipt ,2) = vxc1_out(2*ipt ,2) + bxc_over_m*(-mz1_im + mdirz*m_dot_m1_im ) ! Im[V^22] 816 817 ! v12 += bxc_over_m*( (mx1 - mdirx*m_dot_m1 ) - i.( my1 - mdiry*m_dot_m1 ) ) <= see cplex=1 818 ! Re[v12] += bxc_over_m*( (mx1_re - mdirx*m_dot_m1_re) + ( my1_im - mdiry*m_dot_m1_im) ) 819 ! Im[v12] += bxc_over_m*( (mx1_im - mdirx*m_dot_m1_im) + (-my1_re + mdiry*m_dot_m1_re) ) 820 vxc1_out(2*ipt-1,3) = vxc1_out(2*ipt-1,3) + bxc_over_m*( mx1_re - mdirx*m_dot_m1_re ) ! Re[V^12] 821 vxc1_out(2*ipt-1,3) = vxc1_out(2*ipt-1,3) + bxc_over_m*( my1_im - mdiry*m_dot_m1_im ) ! Re[V^12] 822 vxc1_out(2*ipt ,3) = vxc1_out(2*ipt ,3) + bxc_over_m*( mx1_im - mdirx*m_dot_m1_im ) ! Im[V^12] 823 vxc1_out(2*ipt ,3) = vxc1_out(2*ipt ,3) + bxc_over_m*(-my1_re + mdiry*m_dot_m1_re ) ! Im[V^12] 824 825 ! v21 += bxc_over_m*( (mx1 - mdirx*m_dot_m1 ) + i.( my1 - mdiry*m_dot_m1 ) ) 826 ! Re[v21] += bxc_over_m*( (mx1_re - mdirx*m_dot_m1_re) + (-my1_im + mdiry*m_dot_m1_im) ) 827 ! Im[v21] += bxc_over_m*( (mx1_im - mdirx*m_dot_m1_im) + ( my1_re - mdiry*m_dot_m1_re) ) 828 ! the 4th component is actually not v21, but rather i.v21, this will be adjusted later 829 vxc1_out(2*ipt-1,4) = vxc1_out(2*ipt-1,4) + bxc_over_m*( mx1_re - mdirx*m_dot_m1_re ) ! Re[V^21] 830 vxc1_out(2*ipt-1,4) = vxc1_out(2*ipt-1,4) + bxc_over_m*(-my1_im + mdiry*m_dot_m1_im ) ! Re[V^21] 831 vxc1_out(2*ipt ,4) = vxc1_out(2*ipt ,4) + bxc_over_m*( mx1_im - mdirx*m_dot_m1_im ) ! Im[V^21] 832 vxc1_out(2*ipt ,4) = vxc1_out(2*ipt ,4) + bxc_over_m*( my1_re - mdiry*m_dot_m1_re ) ! Im[V^21] 833 endif 834 else 835 if(option/=0) then 836 !Compute Bxc/|m| from Kxc (|m^(0)| -> zero limit) 837 bxc_over_m = half*(half*(kxc(ipt,1)+kxc(ipt,3))-kxc(ipt,2)) 838 vxc1_out(2*ipt-1,1)= dvdn_re + bxc_over_m*mz1_re 839 vxc1_out(2*ipt-1,2)= dvdn_re - bxc_over_m*mz1_re 840 vxc1_out(2*ipt ,1)= dvdn_im + bxc_over_m*mz1_im 841 vxc1_out(2*ipt ,2)= dvdn_im - bxc_over_m*mz1_im 842 843 vxc1_out(2*ipt-1,3)= bxc_over_m*(mx1_re+my1_im) 844 vxc1_out(2*ipt ,3)= bxc_over_m*(mx1_im-my1_re) 845 vxc1_out(2*ipt-1,4)= bxc_over_m*(mx1_re-my1_im) 846 vxc1_out(2*ipt ,4)= bxc_over_m*(mx1_im+my1_re) 847 else 848 vxc1_out(2*ipt-1,1)= dvdn_re 849 vxc1_out(2*ipt-1,2)= dvdn_re 850 vxc1_out(2*ipt ,1)= dvdn_im 851 vxc1_out(2*ipt ,2)= dvdn_im 852 853 vxc1_out(2*ipt-1,3)= zero 854 vxc1_out(2*ipt ,3)= zero 855 vxc1_out(2*ipt-1,4)= zero 856 vxc1_out(2*ipt ,4)= zero 857 endif 858 end if 859 860 !finally reconstruct i.V^21 from V^21 861 v21tmp(1) = vxc1_out(2*ipt-1,4) !Re[V^21] 862 v21tmp(2) = vxc1_out(2*ipt ,4) !Im[V^21] 863 864 vxc1_out(2*ipt-1,4) =-v21tmp(2) ! Re[i.V^21]=-Im[V^21] 865 vxc1_out(2*ipt ,4) = v21tmp(1) ! Im[i.V^21]= Re[V^21] 866 867 end do ! ipt 868 869 end select !cplex 870 871 end select ! rotation_method 872 873 !DBG_EXIT("COLL") 874 875 end subroutine rotate_back_mag_dfpt
m_xc_noncoll/rotate_mag [ Types ]
[ Top ] [ m_xc_noncoll ] [ Types ]
NAME
rotate_mag
FUNCTION
Project (rotate) a non-collinear density (stored as density+magn.) on a magnetization and give a collinear density (stored as [up,dn] or [up+dn,up]). Align both z-axis.
INPUTS
rho_in(vectsize,4)=input non-collinear density and magnetization (1st or 0th order)
mag(vectsize,3)=gs magnetization used for projection (0th order magnetization)
vectsize=size of vector fields
[mag_norm_in(vectsize)]= --optional-- norm of mag(:) at each point of the grid
[rho_out_format]= 1=rho_out is stored as [up,dn]
2=rho_out is stored as [up+dn,up]
Default=1
rho_out(vectsize,2)=output (projected, collinear) (1st order if rho_in is 1st order NC density matrix)
[mag_norm_out(vectsize)]= --optional-- norm of mag(:) at each point of the grid
Explicit formulae:
rho_out_format=1
rho_out(1) = half*( rho_in(1) + (mag,rho_in(2:4))/|mag|) // where (*,*) is scalar product
rho_out(2) = half*( rho_in(1) - (mag,rho_in(2:4))/|mag|)
rho_out_format=2
rho_out(1) = rho_in(1)
rho_out(2) = half*( rho_in(1) + (mag,rho_in(2:4))/|mag|)
PARENTS
dfpt_mkvxc_noncoll,m_pawxc
CHILDREN
rotate_back_mag_dfpt
SOURCE
99 subroutine rotate_mag(rho_in,rho_out,mag,vectsize,cplex,& 100 & mag_norm_in,mag_norm_out,rho_out_format) ! optional arguments 101 102 103 !This section has been created automatically by the script Abilint (TD). 104 !Do not modify the following lines by hand. 105 #undef ABI_FUNC 106 #define ABI_FUNC 'rotate_mag' 107 !End of the abilint section 108 109 implicit none 110 111 !Arguments ------------------------------------ 112 !scalars 113 integer,intent(in) :: vectsize 114 integer,intent(in) :: cplex 115 integer,intent(in),optional :: rho_out_format 116 !arrays 117 real(dp),intent(in) :: rho_in(cplex*vectsize,4),mag(vectsize,3) 118 real(dp),intent(out) :: rho_out(cplex*vectsize,2) 119 real(dp),intent(in),optional :: mag_norm_in(vectsize) 120 real(dp),intent(out),optional :: mag_norm_out(vectsize) 121 122 !Local variables------------------------------- 123 !scalars 124 integer :: ipt 125 logical :: has_mag_norm,out_mag_norm 126 real(dp) :: m_norm,mm,rho_up,rhoin_dot_mag 127 real(dp) :: rhoin_dot_mag_re,rhoin_dot_mag_im 128 !arrays 129 130 ! ************************************************************************* 131 132 !DBG_ENTER("COLL") 133 134 has_mag_norm=present(mag_norm_in) 135 out_mag_norm=present(mag_norm_out) 136 137 if(cplex==1) then 138 do ipt=1,vectsize 139 if (has_mag_norm) then 140 m_norm=mag_norm_in(ipt) 141 else 142 m_norm=dsqrt(mag(ipt,1)**2+mag(ipt,2)**2+mag(ipt,3)**2) 143 end if 144 145 rhoin_dot_mag=rho_in(ipt,2)*mag(ipt,1)+rho_in(ipt,3)*mag(ipt,2) & 146 & +rho_in(ipt,4)*mag(ipt,3) 147 148 if(m_norm>m_norm_min)then 149 mm=rhoin_dot_mag/m_norm 150 rho_out(ipt,1)=half*(rho_in(ipt,1)+mm) 151 rho_out(ipt,2)=half*(rho_in(ipt,1)-mm) 152 else 153 rho_out(ipt,1)=half*rho_in(ipt,1) 154 rho_out(ipt,2)=half*rho_in(ipt,1) 155 end if 156 157 if (out_mag_norm) then 158 if (m_norm >m_norm_min) mag_norm_out(ipt)=m_norm 159 if (m_norm<=m_norm_min) mag_norm_out(ipt)=zero 160 end if 161 162 end do 163 164 else ! cplex==2 165 166 do ipt=1,vectsize 167 if (has_mag_norm) then 168 m_norm=mag_norm_in(ipt) 169 else 170 m_norm=dsqrt(mag(ipt,1)**2+mag(ipt,2)**2+mag(ipt,3)**2) 171 end if 172 173 ! real part of m.m^(1) 174 rhoin_dot_mag_re=rho_in(2*ipt-1,2)*mag(ipt,1)+rho_in(2*ipt-1,3)*mag(ipt,2) & 175 & +rho_in(2*ipt-1,4)*mag(ipt,3) 176 ! imaginary part of m.m^(1) 177 rhoin_dot_mag_im=rho_in(2*ipt ,2)*mag(ipt,1)+rho_in(2*ipt ,3)*mag(ipt,2) & 178 & +rho_in(2*ipt ,4)*mag(ipt,3) 179 if(m_norm>m_norm_min)then 180 mm=rhoin_dot_mag_re/m_norm 181 rho_out(2*ipt-1,1)=half*(rho_in(2*ipt-1,1)+mm) 182 rho_out(2*ipt-1,2)=half*(rho_in(2*ipt-1,1)-mm) 183 mm=rhoin_dot_mag_im/m_norm 184 rho_out(2*ipt ,1)=half*(rho_in(2*ipt ,1)+mm) 185 rho_out(2*ipt ,2)=half*(rho_in(2*ipt ,1)-mm) 186 else 187 rho_out(2*ipt-1,1)=half*rho_in(2*ipt-1,1) 188 rho_out(2*ipt-1,2)=half*rho_in(2*ipt-1,2) 189 rho_out(2*ipt ,1)=half*rho_in(2*ipt ,1) 190 rho_out(2*ipt ,2)=half*rho_in(2*ipt ,2) 191 end if 192 193 if (out_mag_norm) then 194 if (m_norm >m_norm_min) mag_norm_out(ipt)=m_norm 195 if (m_norm<=m_norm_min) mag_norm_out(ipt)=zero 196 end if 197 198 end do 199 200 end if 201 202 if (present(rho_out_format)) then 203 if (rho_out_format==2) then 204 do ipt=1,cplex*vectsize 205 rho_up=rho_out(ipt,1) 206 rho_out(ipt,1)=rho_up+rho_out(ipt,2) 207 rho_out(ipt,2)=rho_up 208 end do 209 end if 210 end if 211 212 !DBG_EXIT("COLL") 213 214 end subroutine rotate_mag