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-2024 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
SOURCE
21 #if defined HAVE_CONFIG_H 22 #include "config.h" 23 #endif 24 25 #include "abi_common.h" 26 27 MODULE m_xc_noncoll 28 29 use defs_basis 30 use m_abicore 31 use m_errors 32 33 implicit none 34 35 private 36 37 ! public procedures 38 public :: rotate_mag ! Rotate a non-collinear density wrt a magnetization 39 public :: rotate_back_mag ! Rotate back a collinear XC potential wrt a magnetization 40 public :: rotate_back_mag_dfpt ! Rotate back a collinear 1st-order XC potential wrt a magnetization 41 public :: test_rotations ! test whether methods in rotate_back_mag_dfpt give similar results 42 43 !Tolerance on magnetization norm 44 real(dp),parameter :: m_norm_min=tol8 45 46 !Default rotation method for DFPT 47 integer,parameter :: rotation_method_default=3 48 49 CONTAINS 50 51 !===========================================================
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
SOURCE
861 subroutine test_rotations(option,cplex) 862 863 !Arguments ------------------------------------ 864 integer , intent(in) :: option 865 integer , intent(in) :: cplex 866 867 !Local variables------------------------------- 868 real(dp) :: m0(1,3),vxc0(1,4),kxc(1,3) 869 real(dp) :: n1(cplex,4),vxc1_in(cplex,4),vxc1_out(cplex,4) 870 real(dp) :: delta_23(cplex,4) !,delta_12(cplex,4) 871 real(dp) :: m0_norm,dvdn,dvdz,err23 !,wrong_comp!,err12 872 real(dp) :: theta0,phi0,theta1,phi1,err,m1_norm 873 integer :: dir0,dir1 874 ! ************************************************************************* 875 876 DBG_ENTER("COLL") 877 878 ! if (option/=1 .and. option/=2 ) then 879 ! write(msg,'(3a,i0)')& 880 !& 'The argument option should be 1 or 2,',ch10,& 881 !& 'however, option=',option 882 ! ABI_BUG(msg) 883 ! end if 884 ! 885 ! if (sizein<1) then 886 ! write(msg,'(3a,i0)')& 887 !& ' The argument sizein should be a positive number,',ch10,& 888 !& ' however, sizein=',sizein 889 ! ABI_ERROR(msg) 890 ! end if 891 892 DBG_EXIT("COLL") 893 894 !write(*,*) 'VXC_NONCOLL TESTS================================================================' 895 if (cplex==1) then 896 !write(*,*) ' cplex=1------------------------------------------------------------------------' 897 898 !write(*,*) ' TEST: simple m* orietnations, bxc^(1) part' 899 dvdn=zero;dvdz=1.0! 900 err23=zero 901 do dir0=1,3 902 m0=zero; n1=zero 903 do dir1=2,4 904 m0(1,dir0)=0.1 905 m0_norm=sqrt(m0(1,1)**2+m0(1,2)**2+m0(1,3)**2) 906 n1(1,dir1)=0.8 ! any number would do here 907 908 vxc0=zero; ! no bxc^(0) part at all 909 910 vxc1_in(1,1)= dvdn+dvdz 911 vxc1_in(1,2)= dvdn-dvdz 912 913 call rotate_back_mag_dfpt(1,vxc1_in,vxc1_out,vxc0,kxc,n1,m0,1,1,rot_method=1) 914 call rotate_back_mag_dfpt(1,vxc1_in,vxc1_out,vxc0,kxc,n1,m0,1,1,rot_method=2) 915 delta_23=vxc1_out 916 call rotate_back_mag_dfpt(1,vxc1_in,vxc1_out,vxc0,kxc,n1,m0,1,1,rot_method=3) 917 delta_23=abs(delta_23-vxc1_out) 918 err=max(delta_23(1,1),delta_23(1,2),delta_23(1,3),delta_23(1,4)) 919 if (err23<err) err23=err; 920 921 enddo 922 enddo 923 !write(*,*) ' maximum mismatch between methods 2 and 3:',err23 924 925 !write(*,*) ' TEST: simple m* orietnations, bxc^(0) part' 926 927 err23=zero 928 dvdn=zero;dvdz=1.0 929 do dir0=1,3 930 m0=zero; n1=zero 931 do dir1=2,4 932 m0(1,dir0)=0.1 933 m0_norm=sqrt(m0(1,1)**2+m0(1,2)**2+m0(1,3)**2) 934 n1(1,dir1)=0.8 ! =m^1, any number would do here 935 936 vxc0(1,1) = dvdn+dvdz*m0(1,3)/m0_norm 937 vxc0(1,2) = dvdn-dvdz*m0(1,3)/m0_norm 938 vxc0(1,3) = dvdz*m0(1,1)/m0_norm 939 vxc0(1,4) =-dvdz*m0(1,2)/m0_norm 940 941 vxc1_in=zero !vxc^(1) collinear is zero 942 943 call rotate_back_mag_dfpt(1,vxc1_in,vxc1_out,vxc0,kxc,n1,m0,1,1,rot_method=1) 944 call rotate_back_mag_dfpt(1,vxc1_in,vxc1_out,vxc0,kxc,n1,m0,1,1,rot_method=2) 945 delta_23=vxc1_out 946 call rotate_back_mag_dfpt(1,vxc1_in,vxc1_out,vxc0,kxc,n1,m0,1,1,rot_method=3) 947 delta_23=abs(delta_23-vxc1_out) 948 err=maxval(abs(delta_23(1,:))) 949 if (err23<err) err23=err; 950 enddo 951 enddo 952 !write(*,*) ' maximum mismatch between methods 2 and 3:',err23 953 954 !write(*,*) ' TEST: general m0 orietnations, bxc^(0) part' 955 956 theta0=zero 957 err23=zero 958 m0_norm=0.3 959 do while(theta0<=pi) 960 phi0=zero 961 do while(phi0<=2*pi) 962 m0(1,1)=m0_norm*sin(theta0)*cos(phi0) 963 m0(1,2)=m0_norm*sin(theta0)*sin(phi0) 964 m0(1,3)=m0_norm*cos(theta0) 965 966 do dir1=2,4 967 n1=zero 968 n1(1,dir1)=0.8 ! =m^1, any number would do here 969 970 !vxc0=zero; ! 971 vxc0(1,1) = dvdn+dvdz*m0(1,3)/m0_norm 972 vxc0(1,2) = dvdn-dvdz*m0(1,3)/m0_norm 973 vxc0(1,3) = dvdz*m0(1,1)/m0_norm 974 vxc0(1,4) =-dvdz*m0(1,2)/m0_norm 975 976 vxc1_in=zero 977 978 !call rotate_back_mag_dfpt(vxc1_in,vxc1_out,vxc0,kxc,n1,m0,1,1,rot_method=1) 979 call rotate_back_mag_dfpt(1,vxc1_in,vxc1_out,vxc0,kxc,n1,m0,1,1,rot_method=2) 980 delta_23=vxc1_out 981 call rotate_back_mag_dfpt(1,vxc1_in,vxc1_out,vxc0,kxc,n1,m0,1,1,rot_method=3) 982 delta_23=abs(delta_23-vxc1_out) 983 err=maxval(abs(delta_23(1,:))) 984 if (err23<err) err23=err; 985 enddo 986 phi0=phi0+2*pi/100.0 987 enddo 988 theta0=theta0+pi/100.0 989 enddo 990 !write(*,*) ' maximum mismatch between methods 2 and 3:',err23 991 992 if(option==2) then 993 !write(*,*) ' TEST: general m* orietnations, bxc^(0) part' 994 dvdn=zero;dvdz=1.0 995 996 theta0=zero 997 err23=zero 998 m0_norm=0.3 999 m1_norm=10.5 1000 do while(theta0<=pi) !loops on orientation of m^(0) 1001 phi0=zero 1002 do while(phi0<=2*pi) 1003 m0(1,1)=m0_norm*sin(theta0)*cos(phi0) 1004 m0(1,2)=m0_norm*sin(theta0)*sin(phi0) 1005 m0(1,3)=m0_norm*cos(theta0) 1006 1007 vxc0(1,1) = dvdn+dvdz*m0(1,3)/m0_norm 1008 vxc0(1,2) = dvdn-dvdz*m0(1,3)/m0_norm 1009 vxc0(1,3) = dvdz*m0(1,1)/m0_norm 1010 vxc0(1,4) =-dvdz*m0(1,2)/m0_norm 1011 1012 vxc1_in=zero 1013 1014 theta1=zero 1015 do while(theta1<=pi) !loops on orientation of m^(1) 1016 phi1=zero 1017 do while(phi1<=2*pi) 1018 n1(1,1)=zero 1019 n1(1,2)=m1_norm*sin(theta1)*cos(phi1) 1020 n1(1,3)=m1_norm*sin(theta1)*sin(phi1) 1021 n1(1,4)=m1_norm*cos(theta1) 1022 1023 !vxc0=zero; ! 1024 !call rotate_back_mag_dfpt(vxc1_in,vxc1_out,vxc0,kxc,n1,m0,1,1,rot_method=1) 1025 call rotate_back_mag_dfpt(1,vxc1_in,vxc1_out,vxc0,kxc,n1,m0,1,1,rot_method=2) 1026 delta_23=vxc1_out 1027 call rotate_back_mag_dfpt(1,vxc1_in,vxc1_out,vxc0,kxc,n1,m0,1,1,rot_method=3) 1028 delta_23=abs(delta_23-vxc1_out) 1029 err=maxval(abs(delta_23(1,:))) 1030 if (err23<err) err23=err; 1031 phi1=phi1+2*pi/100.0 1032 enddo 1033 theta1=theta1+pi/100.0 1034 enddo 1035 1036 phi0=phi0+2*pi/100.0 1037 enddo 1038 theta0=theta0+pi/100.0 1039 enddo 1040 !write(*,*) ' maximum mismatch between methods 2 and 3:',err23 1041 endif 1042 1043 !else !cplex=2 1044 1045 endif 1046 1047 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
SOURCE
220 subroutine rotate_back_mag(vxc_in,vxc_out,mag,vectsize,& 221 & mag_norm_in) ! optional argument 222 223 !Arguments ------------------------------------ 224 !scalars 225 integer,intent(in) :: vectsize 226 !arrays 227 real(dp),intent(in) :: vxc_in(vectsize,2),mag(vectsize,3) 228 real(dp),intent(out) :: vxc_out(vectsize,4) 229 real(dp),intent(in),optional :: mag_norm_in(vectsize) 230 231 !Local variables------------------------------- 232 !scalars 233 integer :: ipt 234 logical :: has_mag_norm 235 real(dp) :: dvdn,dvdz,m_norm 236 !arrays 237 238 ! ************************************************************************* 239 240 !DBG_ENTER("COLL") 241 242 has_mag_norm=present(mag_norm_in) 243 244 do ipt=1,vectsize 245 if (has_mag_norm) then 246 m_norm=mag_norm_in(ipt) 247 else 248 m_norm=dsqrt(mag(ipt,1)**2+mag(ipt,2)**2+mag(ipt,3)**2) 249 end if 250 251 dvdn=half*(vxc_in(ipt,1)+vxc_in(ipt,2)) 252 253 if (m_norm>m_norm_min) then 254 dvdz=half*(vxc_in(ipt,1)-vxc_in(ipt,2))/m_norm 255 vxc_out(ipt,1)=dvdn+mag(ipt,3)*dvdz 256 vxc_out(ipt,2)=dvdn-mag(ipt,3)*dvdz 257 vxc_out(ipt,3)= mag(ipt,1)*dvdz 258 vxc_out(ipt,4)=-mag(ipt,2)*dvdz 259 else 260 vxc_out(ipt,1:2)=dvdn 261 vxc_out(ipt,3:4)=zero 262 end if 263 end do 264 265 !DBG_EXIT("COLL") 266 267 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
SOURCE
304 subroutine rotate_back_mag_dfpt(option,vxc1_in,vxc1_out,vxc,kxc,rho1,mag,vectsize,cplex,& 305 & mag_norm_in,rot_method) ! optional arguments 306 307 !Arguments ------------------------------------ 308 !scalars 309 integer,intent(in) :: vectsize 310 integer,intent(in) :: cplex 311 integer,intent(in) :: option 312 integer,intent(in),optional :: rot_method 313 !arrays 314 real(dp),intent(in) :: kxc(:,:),mag(vectsize,3),vxc(vectsize,4) 315 real(dp),intent(in) :: rho1(cplex*vectsize,4) 316 real(dp),intent(in) :: vxc1_in(cplex*vectsize,2) 317 real(dp),intent(in),optional :: mag_norm_in(vectsize) 318 real(dp),intent(out) :: vxc1_out(cplex*vectsize,4) 319 320 !Local variables------------------------------- 321 !scalars 322 integer :: ipt,rotation_method 323 logical :: has_mag_norm 324 logical :: small_angle 325 real(dp) :: bxc_over_m,d1,d2,d3,d4,dvdn,dvdz,fact,m_dot_m1,m_norm 326 real(dp) :: dvdn_re,dvdn_im,dvdz_re,dvdz_im 327 complex(dpc) :: rho_updn 328 real(dp) :: mdirx,mdiry,mdirz,mxy,mx1,my1,mz1,wx,wy,wx1,wy1 329 real(dp) :: theta0,theta1,theta1_re,theta1_im 330 real(dp) :: wx1_re,wx1_im 331 real(dp) :: wy1_re,wy1_im 332 real(dp) :: mx1_re,mx1_im,my1_re,my1_im,mz1_re,mz1_im 333 real(dp) :: m_dot_m1_re,m_dot_m1_im 334 real(dp) :: bxc 335 !arrays 336 real(dp) :: vxc_diag(2),v21tmp(2) 337 complex(dpc) :: r1tmp(2,2),u0(2,2),u0_1(2,2),u0_1r1(2,2),u0v1(2,2) 338 complex(dpc) :: rho1_updn(2,2),v1tmp(2,2),vxc1tmp(2,2) 339 complex(dpc) :: rho1_offdiag(2) 340 341 ! ************************************************************************* 342 343 !DBG_ENTER("COLL") 344 345 !Optional arguments 346 has_mag_norm=present(mag_norm_in) 347 rotation_method=rotation_method_default 348 if (present(rot_method)) rotation_method=rot_method 349 350 !Check Kxc 351 if (size(kxc)>3*vectsize) then 352 ABI_ERROR('Cannot use Kxc from GGA!') 353 end if 354 355 if((rotation_method==1.or.rotation_method==2).and.cplex==2) then 356 ABI_ERROR('rotation_method=1 and 2 are not available for cplex=2 case! use ixcrot=3') 357 endif 358 359 360 select case (rotation_method) 361 362 !---------------------------------------- 363 ! Taylor expansion of U rotation matrix 364 !---------------------------------------- 365 case (1) 366 367 do ipt=1,vectsize 368 369 if (has_mag_norm) then 370 m_norm=mag_norm_in(ipt) 371 else 372 m_norm=sqrt(mag(ipt,1)**2+mag(ipt,2)**2+mag(ipt,3)**2) 373 end if 374 375 if(m_norm>m_norm_min) then 376 377 ! Define the U^(0) transformation matrix 378 rho_updn=(mag(ipt,1)+(zero,one)*mag(ipt,2)) 379 d1=sqrt(( m_norm+mag(ipt,3))**2+abs(rho_updn)**2) 380 d2=sqrt((-m_norm+mag(ipt,3))**2+abs(rho_updn)**2) 381 d3=sqrt(( m_norm-mag(ipt,3))**2+abs(rho_updn)**2) 382 d4=sqrt(( m_norm+mag(ipt,3))**2-abs(rho_updn)**2) 383 u0(1,1)=( m_norm+mag(ipt,3))/d1 ! ( m + mz)/d1 384 u0(2,2)=rho_updn/d2 ! ( mx +imy)/d2 385 u0(1,2)=(-m_norm+mag(ipt,3))/d2 ! (-m + mz)/d2 386 u0(2,1)=rho_updn/d1 ! ( mx +imy)/d1 387 388 ! Define the inverse of U^(0): U^(0)^-1 389 if (abs(rho_updn) > m_norm_min) then 390 u0_1(1,1)= half*d1/m_norm 391 u0_1(2,2)= half*d2*(m_norm+mag(ipt,3))/(m_norm*rho_updn) 392 u0_1(1,2)= half*d1*(m_norm-mag(ipt,3))/(m_norm*rho_updn) 393 u0_1(2,1)=-half*d2/m_norm 394 else 395 u0 = zero 396 u0(1,1) = one 397 u0(2,2) = one 398 u0_1 = zero 399 u0_1(1,1) = one 400 u0_1(2,2) = one 401 end if 402 403 ! Diagonalize the GS Vxc^(0): U^(0)^-1 Vxc^(0) U^(0) 404 ! (Remember the abinit notation for vxc!) 405 vxc_diag(1)=half*(vxc(ipt,1)+vxc(ipt,2) & 406 & -sqrt((vxc(ipt,1)-vxc(ipt,2))**2 & 407 & +four*(vxc(ipt,3)**2+vxc(ipt,4)**2))) 408 vxc_diag(2)=half*(vxc(ipt,1)+vxc(ipt,2) & 409 & +sqrt((vxc(ipt,1)-vxc(ipt,2))**2 & 410 & +four*(vxc(ipt,3)**2+vxc(ipt,4)**2))) 411 v1tmp(1,1)=cmplx(real(vxc1_in(ipt,1),kind=dp),zero) 412 v1tmp(2,2)=cmplx(real(vxc1_in(ipt,2),kind=dp),zero) 413 414 !Tranforming the rhor1 with U0 415 rho1_updn(1,1)=half*(rho1(ipt,1)+rho1(ipt,4)) 416 rho1_updn(2,2)=half*(rho1(ipt,1)-rho1(ipt,4)) 417 rho1_updn(1,2)=half*(rho1(ipt,2)-(zero,one)*rho1(ipt,3)) 418 rho1_updn(2,1)=half*(rho1(ipt,2)+(zero,one)*rho1(ipt,3)) 419 u0_1r1=matmul(u0_1,rho1_updn) 420 r1tmp=matmul(u0_1r1,u0) 421 rho1_offdiag(1)=r1tmp(1,2) ; rho1_offdiag(2)=r1tmp(2,1) 422 423 if (option==0) then ! for xccc alone 424 v1tmp(1,2)=cmplx(zero,zero) 425 v1tmp(2,1)=cmplx(zero,zero) 426 else 427 v1tmp(1,2)=-(rho1_offdiag(1)/m_norm)*(vxc_diag(2)-vxc_diag(1)) 428 v1tmp(2,1)= (rho1_offdiag(2)/m_norm)*(vxc_diag(1)-vxc_diag(2)) 429 endif 430 431 !Rotate back the "diagonal" xc computing the term U^(0) Vxc1_^(1) U^(0)^-1 432 u0v1=matmul(u0,v1tmp) 433 vxc1tmp=matmul(u0v1,u0_1) 434 vxc1_out(ipt,1)=real(vxc1tmp(1,1),kind=dp) 435 vxc1_out(ipt,2)=real(vxc1tmp(2,2),kind=dp) 436 vxc1_out(ipt,3)=real( real(vxc1tmp(1,2)),kind=dp) 437 vxc1_out(ipt,4)=real(aimag(vxc1tmp(1,2)),kind=dp) 438 439 else ! Magnetization is zero 440 dvdn=(vxc1_in(ipt,1)+vxc1_in(ipt,2))*half 441 mx1=rho1(ipt,2) ; my1=rho1(ipt,3) ; mz1=rho1(ipt,4) 442 ! Compute Bxc/|m| from Kxc (zero limit) 443 bxc_over_m = half*(half*(kxc(ipt,1)+kxc(ipt,3))-kxc(ipt,2)) 444 vxc1_out(ipt,1)= dvdn + bxc_over_m*mz1 445 vxc1_out(ipt,2)= dvdn - bxc_over_m*mz1 446 vxc1_out(ipt,3)= bxc_over_m*mx1 447 vxc1_out(ipt,4)=-bxc_over_m*my1 448 end if 449 450 end do ! ipt 451 452 !---------------------------------------- 453 ! Analytical expression of U rotation matrix 454 !---------------------------------------- 455 case (2) 456 !Alternative method (explicitely calculated rotation matrices) 457 !Vxc^(1) = phixc^(1).Id + // <= change of "electrostatic" XC potential (phixc^(1) is denoted dvdn) 458 ! + bxc^(1)*( Udag^(0).sigma_z.U^(0) ) + // <= this part describes the change of XC magnetic field magnitude bxc^(1) 459 ! + bxc^(0)*( Udag^(1).sigma_z.U^(0) + Udag^(0).sigma_z.U^(1) ) // <= remaining terms describe the cost of magnetization rotation 460 461 select case(cplex) 462 case(1) 463 do ipt=1,vectsize 464 465 if (has_mag_norm) then 466 m_norm=mag_norm_in(ipt) 467 else 468 m_norm=dsqrt(mag(ipt,1)**2+mag(ipt,2)**2+mag(ipt,3)**2) 469 end if 470 471 472 mx1 =rho1(ipt,2); 473 my1 =rho1(ipt,3); 474 mz1 =rho1(ipt,4) 475 476 dvdn=(vxc1_in(ipt,1)+vxc1_in(ipt,2))*half !phixc^(1) 477 dvdz=(vxc1_in(ipt,1)-vxc1_in(ipt,2))*half !bxc^(1) 478 479 if (m_norm>m_norm_min) then 480 481 mxy = dsqrt(mag(ipt,1)**2+mag(ipt,2)**2) 482 small_angle=(mxy/m_norm<tol8) !condition for sin(x)~x to be valid 483 !even possible to set to tol6 484 mdirx=mag(ipt,1)/m_norm 485 mdiry=mag(ipt,2)/m_norm 486 mdirz=mag(ipt,3)/m_norm 487 488 ! dvdn is phixc^(1) (density only part) 489 ! dvdz is bxc^(1) (magnetization magnitude part) 490 491 !U^(0)*.Vxc1.U^(0) part 492 vxc1_out(ipt,1)= dvdn+dvdz*mdirz 493 vxc1_out(ipt,2)= dvdn-dvdz*mdirz 494 vxc1_out(ipt,3)= dvdz*mdirx ! Real part 495 vxc1_out(ipt,4)=-dvdz*mdiry ! Imaginary part, minus sign comes from sigma_y 496 497 !U^(1)*.Vxc0.U^(0) + U^(0)*.Vxc0.U^(1) part 498 499 !bxc = dsqrt(((vxc(ipt,1)-vxc(ipt,2))*half)**2+vxc(ipt,3)**2+vxc(ipt,4)**2) !bxc^(0) 500 bxc = (vxc(ipt,1)-vxc(ipt,2))*half/mag(ipt,3)*m_norm 501 if (.not.small_angle) then 502 wx = mag(ipt,2)/mxy 503 wy =-mag(ipt,1)/mxy 504 theta0 = dacos(mag(ipt,3)/m_norm) 505 506 theta1 = (mdirz*(mdirx*mx1+mdiry*my1))/mxy - mz1*mxy/m_norm**2 507 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 508 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 509 510 vxc1_out(ipt,1) = vxc1_out(ipt,1) - bxc*dsin(theta0)*theta1 511 vxc1_out(ipt,2) = vxc1_out(ipt,2) + bxc*dsin(theta0)*theta1 512 vxc1_out(ipt,3) = vxc1_out(ipt,3) - bxc*(wy1+dcos(theta0)*wy*theta1) 513 vxc1_out(ipt,4) = vxc1_out(ipt,4) - bxc*(wx1+dcos(theta0)*wx*theta1) 514 else 515 !zero order terms O(1) 516 vxc1_out(ipt,3) = vxc1_out(ipt,3) + bxc*mx1/abs(mag(ipt,3)) 517 vxc1_out(ipt,4) = vxc1_out(ipt,4) - bxc*my1/abs(mag(ipt,3)) 518 !first order terms O(theta) 519 fact = bxc/(mag(ipt,3)*abs(mag(ipt,3))) 520 vxc1_out(ipt,1) = vxc1_out(ipt,1) - (mag(ipt,1)*mx1+mag(ipt,2)*my1)*fact 521 vxc1_out(ipt,2) = vxc1_out(ipt,2) + (mag(ipt,1)*mx1+mag(ipt,2)*my1)*fact 522 vxc1_out(ipt,3) = vxc1_out(ipt,3) - mag(ipt,1)*mz1*fact 523 vxc1_out(ipt,4) = vxc1_out(ipt,4) + mag(ipt,2)*mz1*fact 524 endif 525 526 else ! Magnetization is zero 527 ! Compute Bxc/|m| from Kxc (zero limit) 528 bxc_over_m = half*(half*(kxc(ipt,1)+kxc(ipt,3))-kxc(ipt,2)) 529 vxc1_out(ipt,1)= dvdn + bxc_over_m*mz1 530 vxc1_out(ipt,2)= dvdn - bxc_over_m*mz1 531 vxc1_out(ipt,3)= bxc_over_m*mx1 532 vxc1_out(ipt,4)=-bxc_over_m*my1 533 end if 534 end do ! ipt 535 536 case(2) !cplex=2 537 538 do ipt=1,vectsize 539 540 if (has_mag_norm) then 541 m_norm=mag_norm_in(ipt) 542 else 543 m_norm=dsqrt(mag(ipt,1)**2+mag(ipt,2)**2+mag(ipt,3)**2) 544 end if 545 546 mx1_re= rho1(2*ipt-1,2); mx1_im= rho1(2*ipt,2) 547 my1_re= rho1(2*ipt-1,3); my1_im= rho1(2*ipt,3) 548 mz1_re= rho1(2*ipt-1,4); mz1_im= rho1(2*ipt,4) 549 550 dvdn_re=(vxc1_in(2*ipt-1,1)+vxc1_in(2*ipt-1,2))*half 551 dvdz_re=(vxc1_in(2*ipt-1,1)-vxc1_in(2*ipt-1,2))*half 552 dvdn_im=(vxc1_in(2*ipt ,1)+vxc1_in(2*ipt ,2))*half 553 dvdz_im=(vxc1_in(2*ipt ,1)-vxc1_in(2*ipt ,2))*half 554 555 if (m_norm>m_norm_min) then 556 557 mdirx=mag(ipt,1)/m_norm 558 mdiry=mag(ipt,2)/m_norm 559 mdirz=mag(ipt,3)/m_norm 560 561 mxy = dsqrt(mag(ipt,1)**2+mag(ipt,2)**2) 562 small_angle=(mxy/m_norm<tol8) !condition for sin(x)~x to be valid 563 ! 564 mdirx=mag(ipt,1)/m_norm 565 mdiry=mag(ipt,2)/m_norm 566 mdirz=mag(ipt,3)/m_norm 567 568 ! dvdn is phixc^(1) (density only part) 569 ! dvdz is bxc^(1) (magnetization magnitude part) 570 571 !U^(0)*.Vxc1.U^(0) part 572 vxc1_out(2*ipt-1,1)= dvdn_re+dvdz_re*mdirz 573 vxc1_out(2*ipt ,1)= dvdn_im+dvdz_im*mdirz 574 vxc1_out(2*ipt-1,2)= dvdn_re-dvdz_re*mdirz 575 vxc1_out(2*ipt ,2)= dvdn_im-dvdz_im*mdirz 576 !NOTE: change of definition of the potential matrix components 577 ! vxc1_out(:,3) = V_updn 578 ! vxc1_out(:,4) = i.V_updn 579 vxc1_out(2*ipt-1,3)= dvdz_re*mdirx + dvdz_im*mdiry !Re[ V^12] 580 vxc1_out(2*ipt ,3)= dvdz_im*mdirx - dvdz_re*mdiry !Im[ V^12] 581 582 !U^(1)*.Vxc0.U^(0) + U^(0)*.Vxc0.U^(1) part 583 !bxc = dsqrt(((vxc(ipt,1)-vxc(ipt,2))*half)**2+vxc(ipt,3)**2+vxc(ipt,4)**2) !bxc^(0) 584 bxc = (vxc(ipt,1)-vxc(ipt,2))*half/mag(ipt,3)*m_norm 585 if (.not.small_angle) then 586 wx = mag(ipt,2)/mxy 587 wy =-mag(ipt,1)/mxy 588 theta0 = dacos(mag(ipt,3)/m_norm) 589 590 591 theta1_re = (mdirz*(mdirx*mx1_re+mdiry*my1_re))/mxy - mz1_re*mxy/m_norm**2 592 theta1_im = (mdirz*(mdirx*mx1_im+mdiry*my1_im))/mxy - mz1_im*mxy/m_norm**2 593 594 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 595 wx1_im = (+mag(ipt,1)**2*my1_im - mag(ipt,1)*mag(ipt,2)*mx1_im)/mxy**2/m_norm 596 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 597 wy1_im = (-mag(ipt,2)**2*mx1_im + mag(ipt,1)*mag(ipt,2)*my1_im)/mxy**2/m_norm 598 599 !U^(1)*.Vxc0.U^(0) + U^(0)*.Vxc0.U^(1) 600 vxc1_out(2*ipt-1,1) = vxc1_out(2*ipt-1,1) - bxc*dsin(theta0)*theta1_re 601 vxc1_out(2*ipt ,1) = vxc1_out(2*ipt ,1) - bxc*dsin(theta0)*theta1_im 602 vxc1_out(2*ipt-1,2) = vxc1_out(2*ipt-1,2) + bxc*dsin(theta0)*theta1_re 603 vxc1_out(2*ipt ,2) = vxc1_out(2*ipt ,2) + bxc*dsin(theta0)*theta1_im 604 !cplex=1 part: 605 !v12 += -(bxc)*(wy1+dcos(theta0)*wy*theta1)- 606 ! -i.(bxc)*(wx1+dcos(theta0)*wx*theta1) 607 vxc1_out(2*ipt-1,3) = vxc1_out(2*ipt-1,3) - bxc*(wy1_re+dcos(theta0)*wy*theta1_re) 608 vxc1_out(2*ipt-1,3) = vxc1_out(2*ipt-1,3) + bxc*(wx1_im+dcos(theta0)*wx*theta1_im) 609 vxc1_out(2*ipt ,3) = vxc1_out(2*ipt ,3) - bxc*(wy1_im+dcos(theta0)*wy*theta1_im) 610 vxc1_out(2*ipt ,3) = vxc1_out(2*ipt ,3) - bxc*(wx1_re+dcos(theta0)*wx*theta1_re) 611 else 612 !small theta case: 613 !zero order terms O(1) 614 vxc1_out(2*ipt-1,3) = vxc1_out(2*ipt-1,3) + bxc*mx1_re/abs(mag(ipt,3)) 615 vxc1_out(2*ipt-1,3) = vxc1_out(2*ipt-1,3) + bxc*my1_im/abs(mag(ipt,3)) 616 vxc1_out(2*ipt ,3) = vxc1_out(2*ipt ,3) + bxc*mx1_im/abs(mag(ipt,3)) 617 vxc1_out(2*ipt ,3) = vxc1_out(2*ipt ,3) - bxc*my1_re/abs(mag(ipt,3)) 618 !first order terms: 619 fact = bxc/(mag(ipt,3)*abs(mag(ipt,3))) 620 vxc1_out(2*ipt-1,1) = vxc1_out(2*ipt-1,1) - (mag(ipt,1)*mx1_re+mag(ipt,2)*my1_re)*fact 621 vxc1_out(2*ipt-1,2) = vxc1_out(2*ipt-1,2) + (mag(ipt,1)*mx1_re+mag(ipt,2)*my1_re)*fact 622 vxc1_out(2*ipt ,1) = vxc1_out(2*ipt ,1) - (mag(ipt,1)*mx1_im+mag(ipt,2)*my1_im)*fact 623 vxc1_out(2*ipt ,2) = vxc1_out(2*ipt ,2) + (mag(ipt,1)*mx1_im+mag(ipt,2)*my1_im)*fact 624 625 vxc1_out(2*ipt-1,3) = vxc1_out(2*ipt-1,3) - mag(ipt,1)*mz1_re*fact 626 vxc1_out(2*ipt-1,3) = vxc1_out(2*ipt-1,3) - mag(ipt,2)*mz1_im*fact 627 vxc1_out(2*ipt ,3) = vxc1_out(2*ipt ,3) + mag(ipt,2)*mz1_re*fact 628 vxc1_out(2*ipt ,3) = vxc1_out(2*ipt ,3) - mag(ipt,1)*mz1_im*fact 629 endif 630 631 else ! Magnetization is practically zero 632 ! Compute Bxc/|m| from Kxc (zero limit) 633 bxc_over_m = half*(half*(kxc(ipt,1)+kxc(ipt,3))-kxc(ipt,2)) 634 vxc1_out(2*ipt-1,1)= dvdn_re + bxc_over_m*mz1_re 635 vxc1_out(2*ipt ,1)= dvdn_im + bxc_over_m*mz1_im 636 vxc1_out(2*ipt-1,2)= dvdn_re - bxc_over_m*mz1_re 637 vxc1_out(2*ipt ,2)= dvdn_im - bxc_over_m*mz1_im 638 vxc1_out(2*ipt-1,3)= bxc_over_m*( mx1_re+my1_im) 639 vxc1_out(2*ipt ,3)= bxc_over_m*(-my1_re+mx1_im) 640 end if 641 !finally reconstruct i.V^12 from V^12 642 vxc1_out(2*ipt-1,4) = vxc1_out(2*ipt ,3) ! Re[i.V^21] = Im[V^12] 643 vxc1_out(2*ipt ,4) = vxc1_out(2*ipt-1,3) ! Im[i.V^21] = Re[V^12] 644 645 end do ! ipt 646 647 end select 648 649 !---------------------------------------- 650 ! Explicit derivative of the rotated XC functional 651 !---------------------------------------- 652 case (3) 653 ! Brute-force derivative of Vxc 654 ! Explicit calculation of the rotated xc functional 655 ! (derivatives of the analytical expression) (Eq. A) 656 ! Vxc^(1) = phixc^(1).Id + // <= change of "electrostatic" XC potential (phixc^(1) is denoted dvdn) 657 ! + bxc^(1)*(sigma,m^(0))/|m^(0)| + // <= this term is equivalent to ( Udag^(0).sigma_z.U^(0) ) term in rotation_method=2 658 ! + 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) ) 659 ! - bxc^(0)*(sigma,m^(0))*(m^(1),m^(0))/|m^(0)|**3 660 select case(cplex) 661 case(1) 662 663 do ipt=1,vectsize 664 665 if (has_mag_norm) then 666 m_norm=mag_norm_in(ipt) 667 else 668 m_norm=sqrt(mag(ipt,1)**2+mag(ipt,2)**2+mag(ipt,3)**2) 669 end if 670 671 ! dvdn is phixc^(1) (density only part) 672 ! dvdz is bxc^(1) (magnetization magnitude part) 673 dvdn=(vxc1_in(ipt,1)+vxc1_in(ipt,2))*half 674 dvdz=(vxc1_in(ipt,1)-vxc1_in(ipt,2))*half 675 676 mx1=rho1(ipt,2) ; my1=rho1(ipt,3) ; mz1=rho1(ipt,4) 677 678 if(m_norm>m_norm_min) then 679 680 mdirx=mag(ipt,1)/m_norm; mdiry=mag(ipt,2)/m_norm; mdirz=mag(ipt,3)/m_norm 681 682 !This part describes the change of the magnitude of the xc magnetic field 683 !and the change of the scalar part of the xc electrostatic potential, 1st + 2nd term in Eq.A 684 !phixc^(1).Id + bxc^(1) (sigma,m^(0))/|m^(0)| 685 vxc1_out(ipt,1)= dvdn+dvdz*mdirz 686 vxc1_out(ipt,2)= dvdn-dvdz*mdirz 687 vxc1_out(ipt,3)= dvdz*mdirx ! Real part 688 vxc1_out(ipt,4)=-dvdz*mdiry ! Imaginary part, minus sign comes from sigma_y 689 690 if (option/=0) then 691 !Add remaining contributions comming from the change of magnetization direction 692 !projection of m^(1) on gs magnetization direction 693 m_dot_m1=(mdirx*rho1(ipt,2)+mdiry*rho1(ipt,3)+mdirz*rho1(ipt,4)) 694 695 bxc_over_m =-dsqrt(((vxc(ipt,1)-vxc(ipt,2))*half)**2+vxc(ipt,3)**2+vxc(ipt,4)**2) !this is bxc^(0) 696 bxc_over_m = bxc_over_m/m_norm 697 vxc1_out(ipt,1) = vxc1_out(ipt,1) + bxc_over_m*( mz1 - mdirz*m_dot_m1 ) ! 698 vxc1_out(ipt,2) = vxc1_out(ipt,2) + bxc_over_m*(-mz1 + mdirz*m_dot_m1 ) ! 699 vxc1_out(ipt,3) = vxc1_out(ipt,3) + bxc_over_m*( mx1 - mdirx*m_dot_m1 ) ! 700 vxc1_out(ipt,4) = vxc1_out(ipt,4) + bxc_over_m*(-my1 + mdiry*m_dot_m1 ) ! 701 endif 702 703 else 704 if (option/=0) then 705 !Compute bxc^(0)/|m| from kxc (|m^(0)| -> zero limit) 706 bxc_over_m = half*(half*(kxc(ipt,1)+kxc(ipt,3))-kxc(ipt,2)) 707 vxc1_out(ipt,1)= dvdn + bxc_over_m*mz1 708 vxc1_out(ipt,2)= dvdn - bxc_over_m*mz1 709 vxc1_out(ipt,3)= bxc_over_m*mx1 710 vxc1_out(ipt,4)=-bxc_over_m*my1 711 else 712 vxc1_out(ipt,1)= dvdn 713 vxc1_out(ipt,2)= dvdn 714 vxc1_out(ipt,3)= zero 715 vxc1_out(ipt,4)= zero 716 endif 717 end if 718 719 end do ! ipt 720 721 case(2) 722 !cplex=2 case 723 724 do ipt=1,vectsize 725 726 if (has_mag_norm) then 727 m_norm=mag_norm_in(ipt) 728 else 729 m_norm=sqrt(mag(ipt,1)**2+mag(ipt,2)**2+mag(ipt,3)**2) 730 end if 731 732 ! see cplex=1 case for details 733 dvdn_re=(vxc1_in(2*ipt-1,1)+vxc1_in(2*ipt-1,2))*half 734 dvdn_im=(vxc1_in(2*ipt ,1)+vxc1_in(2*ipt ,2))*half 735 dvdz_re=(vxc1_in(2*ipt-1,1)-vxc1_in(2*ipt-1,2))*half 736 dvdz_im=(vxc1_in(2*ipt ,1)-vxc1_in(2*ipt ,2))*half 737 738 mx1_re=rho1(2*ipt-1,2); mx1_im=rho1(2*ipt,2) 739 my1_re=rho1(2*ipt-1,3); my1_im=rho1(2*ipt,3) 740 mz1_re=rho1(2*ipt-1,4); mz1_im=rho1(2*ipt,4) 741 742 if(m_norm>m_norm_min) then 743 744 mdirx=mag(ipt,1)/m_norm; mdiry=mag(ipt,2)/m_norm; mdirz=mag(ipt,3)/m_norm 745 746 !first two terms: 747 vxc1_out(2*ipt-1,1)= dvdn_re+dvdz_re*mdirz 748 vxc1_out(2*ipt ,1)= dvdn_im+dvdz_im*mdirz 749 vxc1_out(2*ipt-1,2)= dvdn_re-dvdz_re*mdirz 750 vxc1_out(2*ipt ,2)= dvdn_im-dvdz_im*mdirz 751 !NOTE: change of definition of the potential matrix components 752 ! vxc1_out(:,3) = V_updn 753 ! vxc1_out(:,4) = i.V_dnup 754 755 ! 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) 756 ! 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|) 757 vxc1_out(2*ipt-1,3)= dvdz_re*mdirx + dvdz_im*mdiry !Re[V^12] 758 vxc1_out(2*ipt ,3)= dvdz_im*mdirx - dvdz_re*mdiry !Im[V^12] 759 vxc1_out(2*ipt-1,4)= dvdz_re*mdirx - dvdz_im*mdiry !Re[V^21] 760 vxc1_out(2*ipt ,4)= dvdz_im*mdirx + dvdz_re*mdiry !Im[V^21] 761 if (option/=0) then 762 763 !remaining contributions: 764 m_dot_m1_re= mdirx*mx1_re + mdiry*my1_re + mdirz*mz1_re 765 m_dot_m1_im= mdirx*mx1_im + mdiry*my1_im + mdirz*mz1_im 766 767 bxc_over_m =-dsqrt(((vxc(ipt,1)-vxc(ipt,2))*half)**2+vxc(ipt,3)**2+vxc(ipt,4)**2) !this is bxc^(0) 768 bxc_over_m = bxc_over_m/m_norm 769 !bxc_over_m = (vxc(ipt,1)-vxc(ipt,2))*half/mag(ipt,3) 770 771 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] 772 vxc1_out(2*ipt ,1) = vxc1_out(2*ipt ,1) + bxc_over_m*( mz1_im - mdirz*m_dot_m1_im ) ! Im[V^11] 773 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] 774 vxc1_out(2*ipt ,2) = vxc1_out(2*ipt ,2) + bxc_over_m*(-mz1_im + mdirz*m_dot_m1_im ) ! Im[V^22] 775 776 ! v12 += bxc_over_m*( (mx1 - mdirx*m_dot_m1 ) - i.( my1 - mdiry*m_dot_m1 ) ) <= see cplex=1 777 ! Re[v12] += bxc_over_m*( (mx1_re - mdirx*m_dot_m1_re) + ( my1_im - mdiry*m_dot_m1_im) ) 778 ! Im[v12] += bxc_over_m*( (mx1_im - mdirx*m_dot_m1_im) + (-my1_re + mdiry*m_dot_m1_re) ) 779 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] 780 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] 781 vxc1_out(2*ipt ,3) = vxc1_out(2*ipt ,3) + bxc_over_m*( mx1_im - mdirx*m_dot_m1_im ) ! Im[V^12] 782 vxc1_out(2*ipt ,3) = vxc1_out(2*ipt ,3) + bxc_over_m*(-my1_re + mdiry*m_dot_m1_re ) ! Im[V^12] 783 784 ! v21 += bxc_over_m*( (mx1 - mdirx*m_dot_m1 ) + i.( my1 - mdiry*m_dot_m1 ) ) 785 ! Re[v21] += bxc_over_m*( (mx1_re - mdirx*m_dot_m1_re) + (-my1_im + mdiry*m_dot_m1_im) ) 786 ! Im[v21] += bxc_over_m*( (mx1_im - mdirx*m_dot_m1_im) + ( my1_re - mdiry*m_dot_m1_re) ) 787 ! the 4th component is actually not v21, but rather i.v21, this will be adjusted later 788 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] 789 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] 790 vxc1_out(2*ipt ,4) = vxc1_out(2*ipt ,4) + bxc_over_m*( mx1_im - mdirx*m_dot_m1_im ) ! Im[V^21] 791 vxc1_out(2*ipt ,4) = vxc1_out(2*ipt ,4) + bxc_over_m*( my1_re - mdiry*m_dot_m1_re ) ! Im[V^21] 792 endif 793 else 794 if(option/=0) then 795 !Compute Bxc/|m| from Kxc (|m^(0)| -> zero limit) 796 bxc_over_m = half*(half*(kxc(ipt,1)+kxc(ipt,3))-kxc(ipt,2)) 797 vxc1_out(2*ipt-1,1)= dvdn_re + bxc_over_m*mz1_re 798 vxc1_out(2*ipt-1,2)= dvdn_re - bxc_over_m*mz1_re 799 vxc1_out(2*ipt ,1)= dvdn_im + bxc_over_m*mz1_im 800 vxc1_out(2*ipt ,2)= dvdn_im - bxc_over_m*mz1_im 801 802 vxc1_out(2*ipt-1,3)= bxc_over_m*(mx1_re+my1_im) 803 vxc1_out(2*ipt ,3)= bxc_over_m*(mx1_im-my1_re) 804 vxc1_out(2*ipt-1,4)= bxc_over_m*(mx1_re-my1_im) 805 vxc1_out(2*ipt ,4)= bxc_over_m*(mx1_im+my1_re) 806 else 807 vxc1_out(2*ipt-1,1)= dvdn_re 808 vxc1_out(2*ipt-1,2)= dvdn_re 809 vxc1_out(2*ipt ,1)= dvdn_im 810 vxc1_out(2*ipt ,2)= dvdn_im 811 812 vxc1_out(2*ipt-1,3)= zero 813 vxc1_out(2*ipt ,3)= zero 814 vxc1_out(2*ipt-1,4)= zero 815 vxc1_out(2*ipt ,4)= zero 816 endif 817 end if 818 819 !finally reconstruct i.V^21 from V^21 820 v21tmp(1) = vxc1_out(2*ipt-1,4) !Re[V^21] 821 v21tmp(2) = vxc1_out(2*ipt ,4) !Im[V^21] 822 823 vxc1_out(2*ipt-1,4) =-v21tmp(2) ! Re[i.V^21]=-Im[V^21] 824 vxc1_out(2*ipt ,4) = v21tmp(1) ! Im[i.V^21]= Re[V^21] 825 826 end do ! ipt 827 828 end select !cplex 829 830 end select ! rotation_method 831 832 !DBG_EXIT("COLL") 833 834 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|)
SOURCE
88 subroutine rotate_mag(rho_in,rho_out,mag,vectsize,cplex,& 89 & mag_norm_in,mag_norm_out,rho_out_format) ! optional arguments 90 91 !Arguments ------------------------------------ 92 !scalars 93 integer,intent(in) :: vectsize 94 integer,intent(in) :: cplex 95 integer,intent(in),optional :: rho_out_format 96 !arrays 97 real(dp),intent(in) :: rho_in(cplex*vectsize,4),mag(vectsize,3) 98 real(dp),intent(out) :: rho_out(cplex*vectsize,2) 99 real(dp),intent(in),optional :: mag_norm_in(vectsize) 100 real(dp),intent(out),optional :: mag_norm_out(vectsize) 101 102 !Local variables------------------------------- 103 !scalars 104 integer :: ipt 105 logical :: has_mag_norm,out_mag_norm 106 real(dp) :: m_norm,mm,rho_up,rhoin_dot_mag 107 real(dp) :: rhoin_dot_mag_re,rhoin_dot_mag_im 108 !arrays 109 110 ! ************************************************************************* 111 112 !DBG_ENTER("COLL") 113 114 has_mag_norm=present(mag_norm_in) 115 out_mag_norm=present(mag_norm_out) 116 117 if(cplex==1) then 118 do ipt=1,vectsize 119 if (has_mag_norm) then 120 m_norm=mag_norm_in(ipt) 121 else 122 m_norm=dsqrt(mag(ipt,1)**2+mag(ipt,2)**2+mag(ipt,3)**2) 123 end if 124 125 rhoin_dot_mag=rho_in(ipt,2)*mag(ipt,1)+rho_in(ipt,3)*mag(ipt,2) & 126 & +rho_in(ipt,4)*mag(ipt,3) 127 128 if(m_norm>m_norm_min)then 129 mm=rhoin_dot_mag/m_norm 130 rho_out(ipt,1)=half*(rho_in(ipt,1)+mm) 131 rho_out(ipt,2)=half*(rho_in(ipt,1)-mm) 132 else 133 rho_out(ipt,1)=half*rho_in(ipt,1) 134 rho_out(ipt,2)=half*rho_in(ipt,1) 135 end if 136 137 if (out_mag_norm) then 138 if (m_norm >m_norm_min) mag_norm_out(ipt)=m_norm 139 if (m_norm<=m_norm_min) mag_norm_out(ipt)=zero 140 end if 141 142 end do 143 144 else ! cplex==2 145 146 do ipt=1,vectsize 147 if (has_mag_norm) then 148 m_norm=mag_norm_in(ipt) 149 else 150 m_norm=dsqrt(mag(ipt,1)**2+mag(ipt,2)**2+mag(ipt,3)**2) 151 end if 152 153 ! real part of m.m^(1) 154 rhoin_dot_mag_re=rho_in(2*ipt-1,2)*mag(ipt,1)+rho_in(2*ipt-1,3)*mag(ipt,2) & 155 & +rho_in(2*ipt-1,4)*mag(ipt,3) 156 ! imaginary part of m.m^(1) 157 rhoin_dot_mag_im=rho_in(2*ipt ,2)*mag(ipt,1)+rho_in(2*ipt ,3)*mag(ipt,2) & 158 & +rho_in(2*ipt ,4)*mag(ipt,3) 159 if(m_norm>m_norm_min)then 160 mm=rhoin_dot_mag_re/m_norm 161 rho_out(2*ipt-1,1)=half*(rho_in(2*ipt-1,1)+mm) 162 rho_out(2*ipt-1,2)=half*(rho_in(2*ipt-1,1)-mm) 163 mm=rhoin_dot_mag_im/m_norm 164 rho_out(2*ipt ,1)=half*(rho_in(2*ipt ,1)+mm) 165 rho_out(2*ipt ,2)=half*(rho_in(2*ipt ,1)-mm) 166 else 167 rho_out(2*ipt-1,1)=half*rho_in(2*ipt-1,1) 168 rho_out(2*ipt-1,2)=half*rho_in(2*ipt-1,2) 169 rho_out(2*ipt ,1)=half*rho_in(2*ipt ,1) 170 rho_out(2*ipt ,2)=half*rho_in(2*ipt ,2) 171 end if 172 173 if (out_mag_norm) then 174 if (m_norm >m_norm_min) mag_norm_out(ipt)=m_norm 175 if (m_norm<=m_norm_min) mag_norm_out(ipt)=zero 176 end if 177 178 end do 179 180 end if 181 182 if (present(rho_out_format)) then 183 if (rho_out_format==2) then 184 do ipt=1,cplex*vectsize 185 rho_up=rho_out(ipt,1) 186 rho_out(ipt,1)=rho_up+rho_out(ipt,2) 187 rho_out(ipt,2)=rho_up 188 end do 189 end if 190 end if 191 192 !DBG_EXIT("COLL") 193 194 end subroutine rotate_mag