TABLE OF CONTENTS
ABINIT/m_spin_current [ Modules ]
NAME
m_spin_current
FUNCTION
COPYRIGHT
Copyright (C) 2005-2024 ABINIT group (Mver) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_spin_current 23 24 use m_errors 25 use m_abicore 26 use m_splines 27 use m_hdr 28 use m_dtset 29 use m_dtfil 30 31 use defs_datatypes, only : pseudopotential_type 32 use defs_abitypes, only : MPI_type 33 use m_io_tools, only : open_file 34 use m_pptools, only : printxsf 35 use m_geometry, only : xred2xcart 36 use m_fftcore, only : sphereboundary 37 use m_special_funcs, only : gamma_function 38 use m_fft, only : fourwf 39 40 implicit none 41 42 private
m_spin_current/spin_current [ Functions ]
[ Top ] [ m_spin_current ] [ Functions ]
NAME
spin_current
FUNCTION
INPUTS
atindx(natom)=index table for atoms (see gstate.f) atindx1(natom)=inverse of atindx cg(2,mcg)=wavefunctions (may be read from disk instead of input) dtfil <type(datafiles_type)>=variables related to files dtset <type(dataset_type)>=all input variables in this dataset eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree) gmet = reciprocal space metric gprimd = dimensionful reciprocal space vectors hdr <type(hdr_type)>=the header of wf, den and pot files kg(3,mpw*mkmem)=reduced (integer) coordinates of G vecs in basis sphere mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol mpi_enreg=information about MPI parallelization nattyp(dtset%ntypat)=number of atoms of each type nfftf = fft grid dimensions for fine grid ph1d = phase factors in 1 radial dimension psps <type(pseudopotential_type)>=variables related to pseudopotentials | mpsang= 1+maximum angular momentum rhog(2,nfftf)=Fourier transform of total electron density (including compensation density in PAW) rhor(nfftf,nspden)=total electron density (including compensation density in PAW) rmet = real space metric tensor symrec(3,3,nsym)=symmetries in reciprocal space, reduced coordinates ucvol = unit cell volume wffnow=unit number for current wf disk file ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point ylmgr(mpw*mkmem,3,mpsang*mpsang*useylm)= gradients of real spherical harmonics
OUTPUT
only output to file
SOURCE
89 subroutine spin_current(cg,dtfil,dtset,gprimd,hdr,kg,mcg,mpi_enreg,psps) 90 91 !Arguments ------------------------------------ 92 !scalars 93 !integer,intent(in) :: nfftf 94 !real(dp),intent(in) :: ucvol 95 integer,intent(in) :: mcg 96 type(MPI_type),intent(in) :: mpi_enreg 97 type(datafiles_type),intent(in) :: dtfil 98 type(dataset_type),intent(in) :: dtset 99 type(hdr_type),intent(inout) :: hdr 100 type(pseudopotential_type),intent(in) :: psps 101 !type(wffile_type),intent(in) :: wffnow 102 !arrays 103 !integer,intent(in) :: atindx(dtset%natom),atindx1(dtset%natom) 104 integer,intent(in) :: kg(3,dtset%mpw*dtset%mkmem) 105 !integer,intent(in) :: nattyp(dtset%ntypat) 106 !integer,intent(in) :: symrec(3,3,dtset%nsym) 107 real(dp),intent(in) :: cg(2,mcg) 108 !real(dp),intent(in) :: eigen(dtset%mband*dtset%nkpt*dtset%nsppol),gmet(3,3) 109 real(dp),intent(in) :: gprimd(3,3) 110 !real(dp),intent(in) :: rhog(2,nfftf),rhor(nfftf,dtset%nspden) 111 !real(dp),intent(in) :: rmet(3,3) 112 !real(dp),intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm) 113 !real(dp),intent(in) :: ylmgr(dtset%mpw*dtset%mkmem,3,psps%mpsang*psps%mpsang*psps%useylm) 114 !real(dp),intent(inout) :: ph1d(2,3*(2*dtset%mgfft+1)*dtset%natom) 115 116 !Local variables------------------------------- 117 !scalars 118 integer :: cplex,fft_option,i1 119 integer :: i2,i3,iband,icartdir,icg,ig 120 integer :: ikg,ikpt,iocc,irealsp,ispindir,ispinor,ispinorp 121 integer :: npw 122 integer :: icplex 123 integer :: realrecip 124 integer :: iatom,spcur_unit 125 real(dp) :: prefact_nk 126 real(dp) :: rescale_current 127 character(len=500) :: message 128 character(len=fnlen) :: filnam 129 !arrays 130 integer,allocatable :: gbound(:,:),kg_k(:,:) 131 real(dp),allocatable :: dpsidr(:,:,:,:,:,:) 132 real(dp),allocatable :: density(:,:,:,:) 133 real(dp),allocatable :: dummy_denpot(:,:,:) 134 real(dp),allocatable :: gpsi(:,:,:,:),kgcart(:,:) 135 real(dp),allocatable :: position_op(:,:,:,:) 136 real(dp),allocatable :: psi(:,:,:),psi_r(:,:,:,:,:) 137 real(dp),allocatable :: spincurrent(:,:,:,:,:) 138 real(dp),allocatable :: vso_realspace(:,:,:,:,:),datagrid(:) 139 real(dp) :: dummy_fofgout(0,0) 140 real(dp),allocatable :: xcart(:,:) 141 character :: spin_symbol(3) 142 character :: spinor_sym(2) 143 character(len=2) :: realimag(2) 144 !no_abirules 145 !real(dp),allocatable :: density_matrix(:,:,:,:,:) 146 !real(dp),allocatable :: vso_realspace_nl(:,:,:,:,:) 147 148 ! ************************************************************************* 149 150 !write(std_out,*) ' Entering subroutine spin_current ' 151 !write(std_out,*) ' dtset%ngfft = ', dtset%ngfft 152 !write(std_out,*) ' hdr%istwfk = ', hdr%istwfk 153 154 !===================== init and checks ============================ 155 !check if nspinor is 2 156 if (dtset%nspinor /= 2) then 157 write(message, '(a,i0)' )' nspinor must be 2, but it is ',dtset%nspinor 158 ABI_ERROR(message) 159 end if 160 161 if (dtset%nsppol /= 1) then 162 write(message, '(a,i0)' )' spin_current: nsppol must be 1 but it is ',dtset%nsppol 163 ABI_ERROR(message) 164 end if 165 166 if (dtset%mkmem /= dtset%nkpt) then 167 write(message, '(a,i6,a,i6,a,a)' )& 168 & ' mkmem = ',dtset%mkmem,' must be equal to nkpt ',dtset%nkpt,ch10,& 169 & ' keep all kpt in memory' 170 ABI_ERROR(message) 171 end if 172 173 if (dtset%usepaw /= 0) then 174 write(message, '(a,i0,a,a,a)' )& 175 & 'usepaw = ',dtset%usepaw,' must be equal to 0 ',ch10,& 176 & 'Not functional for PAW case yet.' 177 ABI_ERROR(message) 178 end if 179 180 cplex=2 181 fft_option = 0 ! just do direct fft 182 spin_symbol = (/'x','y','z'/) 183 spinor_sym = (/'u','d'/) 184 realimag = (/'Re','Im'/) 185 186 write(std_out,*) ' psps%mpsang,psps%mpssoang ', psps%mpsang,psps%mpssoang 187 188 !======================= main code ================================ 189 !----------------------------------------------------------------------------------------- 190 !----------------------------------------------------------------------------------------- 191 !first get normal contribution to current, as psi tau dpsidr + dpsidr tau psi 192 !where tau are 1/2 the pauli matrices 193 !----------------------------------------------------------------------------------------- 194 !----------------------------------------------------------------------------------------- 195 196 !init plane wave coeff counter 197 icg = 0 198 !init plane wave counter 199 ikg = 0 200 !init occupation/band counter 201 iocc = 1 202 203 !rspace point, cartesian direction, spin pol=x,y,z 204 ABI_MALLOC(spincurrent,(dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6),3,3)) 205 spincurrent = zero 206 207 ABI_MALLOC(dummy_denpot,(cplex*dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6))) 208 209 ABI_MALLOC(gbound,(2*dtset%mgfft+8,2)) 210 211 !allocate (density_matrix(2,dtset%ngfft(1)*dtset%ngfft(2)*dtset%ngfft(3),dtset%nspinor,& 212 !& dtset%ngfft(1)*dtset%ngfft(2)*dtset%ngfft(3),dtset%nspinor)) 213 !density_matrix= zero 214 ABI_MALLOC(density,(2,dtset%ngfft(1)*dtset%ngfft(2)*dtset%ngfft(3),dtset%nspinor,dtset%nspinor)) 215 density= zero 216 217 ABI_MALLOC(dpsidr,(2,dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6),dtset%nspinor,3)) 218 ABI_MALLOC(psi_r,(2,dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6),dtset%nspinor)) 219 220 !loop over kpoints 221 do ikpt=1,dtset%nkpt 222 223 ! number of plane waves for this kpt 224 npw = hdr%npwarr(ikpt) 225 226 ! allocate arrays dep on number of pw 227 ABI_MALLOC(kg_k,(3,npw)) 228 ABI_MALLOC(gpsi,(2,npw,dtset%nspinor,3)) 229 ABI_MALLOC(psi,(2,npw,dtset%nspinor)) 230 ABI_MALLOC(kgcart,(3,npw)) 231 232 ! get cartesian coordinates of k+G vectors around this kpoint 233 do ig=1,npw 234 kgcart(:,ig) = matmul(gprimd(:,:),dtset%kpt(:,ikpt)+kg(:,ikg+ig)) 235 kg_k (:,ig) = kg(:,ikg+ig) 236 end do 237 238 ! get gbound 239 call sphereboundary(gbound,dtset%istwfk(ikpt),kg_k,dtset%mgfft,npw) 240 241 ! loop over bands 242 do iband=1,dtset%nband(ikpt) 243 244 ! prefactor for sum over bands and kpoints 245 prefact_nk = hdr%occ(iocc) * dtset%wtk(ikpt) 246 247 ! initialize this wf 248 gpsi=zero 249 psi=zero 250 psi(:,1:npw,1) = cg(:,icg+1:icg+npw) 251 252 ! multiply psi by - i 2 pi G 253 do ig=1,npw 254 gpsi(1,ig,:,1) = two_pi * kgcart(1,ig)*psi(2,ig,:) 255 gpsi(2,ig,:,1) = -two_pi * kgcart(1,ig)*psi(1,ig,:) 256 gpsi(1,ig,:,2) = two_pi * kgcart(2,ig)*psi(2,ig,:) 257 gpsi(2,ig,:,2) = -two_pi * kgcart(2,ig)*psi(1,ig,:) 258 gpsi(1,ig,:,3) = two_pi * kgcart(3,ig)*psi(2,ig,:) 259 gpsi(2,ig,:,3) = -two_pi * kgcart(3,ig)*psi(1,ig,:) 260 end do 261 262 ! loop over spinorial components 263 do ispinor=1,dtset%nspinor 264 ! FT Gpsi_x to real space 265 call fourwf(cplex,dummy_denpot,gpsi(:,:,ispinor,1),dummy_fofgout,& 266 & dpsidr(:,:,:,:,ispinor,1),gbound,gbound,& 267 & hdr%istwfk(ikpt),kg_k,kg_k,dtset%mgfft,mpi_enreg,1,dtset%ngfft,npw,& 268 & npw,dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6),& 269 & fft_option,0,one,one,gpu_option=dtset%gpu_option) 270 271 ! FT Gpsi_y to real space 272 call fourwf(cplex,dummy_denpot,gpsi(:,:,ispinor,2),dummy_fofgout,& 273 & dpsidr(:,:,:,:,ispinor,2),gbound,gbound,& 274 & hdr%istwfk(ikpt),kg_k,kg_k,dtset%mgfft,mpi_enreg,1,dtset%ngfft,npw,& 275 & npw,dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6),& 276 & fft_option,0,one,one,gpu_option=dtset%gpu_option) 277 278 ! FT Gpsi_z to real space 279 call fourwf(cplex,dummy_denpot,gpsi(:,:,ispinor,3),dummy_fofgout,& 280 & dpsidr(:,:,:,:,ispinor,3),gbound,gbound,& 281 & hdr%istwfk(ikpt),kg_k,kg_k,dtset%mgfft,mpi_enreg,1,dtset%ngfft,npw,& 282 & npw,dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6),& 283 & fft_option,0,one,one,gpu_option=dtset%gpu_option) 284 285 ! FT psi to real space 286 call fourwf(cplex,dummy_denpot,psi(:,:,ispinor),dummy_fofgout,& 287 & psi_r(:,:,:,:,ispinor),gbound,gbound,& 288 & hdr%istwfk(ikpt),kg_k,kg_k,dtset%mgfft,mpi_enreg,1,dtset%ngfft,npw,& 289 & npw,dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6),& 290 & fft_option,0,one,one,gpu_option=dtset%gpu_option) 291 292 end do ! ispinor 293 294 ! dpsidr now contains the full derivative of psi wrt space (gradient) in cartesian coordinates 295 296 ! get 3 pauli matrix contributions to the current: x,y,z, cart dir, spin dir 297 do icartdir=1,3 298 299 ! x pauli spin matrix 300 ! sigma_x = | 0 1 | 301 ! | 1 0 | 302 spincurrent(:,:,:,icartdir,1) = spincurrent(:,:,:,icartdir,1) + prefact_nk * & 303 ! Re(psi_r(up)^* dpsidr(down)) 304 & real(psi_r(1,:,:,:,1)*dpsidr(1,:,:,:,2,icartdir) & 305 & + psi_r(2,:,:,:,1)*dpsidr(2,:,:,:,2,icartdir) & 306 ! Re(psi_r(down)^* dpsidr(up)) 307 & + psi_r(1,:,:,:,2)*dpsidr(1,:,:,:,1,icartdir) & 308 & + psi_r(2,:,:,:,2)*dpsidr(2,:,:,:,1,icartdir)) 309 310 ! y pauli spin matrix 311 ! sigma_y = | 0 -i | 312 ! | i 0 | 313 spincurrent(:,:,:,icartdir,2) = spincurrent(:,:,:,icartdir,2) + prefact_nk * & 314 ! Re(-i psi_r(up)^* dpsidr(down)) 315 & real(psi_r(1,:,:,:,1)*dpsidr(2,:,:,:,2,icartdir) & 316 & - psi_r(2,:,:,:,1)*dpsidr(1,:,:,:,2,icartdir) & 317 ! Re(i psi_r(down)^* dpsidr(up)) 318 & - psi_r(1,:,:,:,2)*dpsidr(2,:,:,:,1,icartdir) & 319 & + psi_r(2,:,:,:,2)*dpsidr(1,:,:,:,1,icartdir)) 320 321 ! z pauli spin matrix 322 ! sigma_z = | 1 0 | 323 ! | 0 -1 | 324 spincurrent(:,:,:,icartdir,3) = spincurrent(:,:,:,icartdir,3) + prefact_nk * & 325 ! Re(psi_r(up)^* dpsidr(up)) 326 & real(psi_r(1,:,:,:,1)*dpsidr(1,:,:,:,1,icartdir) & 327 & - psi_r(2,:,:,:,1)*dpsidr(2,:,:,:,1,icartdir) & 328 ! Re(-psi_r(down)^* dpsidr(down)) 329 & - psi_r(1,:,:,:,2)*dpsidr(1,:,:,:,2,icartdir) & 330 & + psi_r(2,:,:,:,2)*dpsidr(2,:,:,:,2,icartdir)) 331 end do ! end icartdir 332 333 ! 334 ! accumulate non local density matrix in real space 335 ! NOTE: if we are only using the local part of the current, this becomes the 336 ! density spinor matrix! (much lighter to calculate) rho(r, sigma, sigmaprime) 337 ! 338 do ispinor=1,dtset%nspinor 339 do i3=1,dtset%ngfft(3) 340 do i2=1,dtset%ngfft(2) 341 do i1=1,dtset%ngfft(1) 342 irealsp = i1 + (i2-1)*dtset%ngfft(1) + (i3-1)*dtset%ngfft(2)*dtset%ngfft(1) 343 344 do ispinorp=1,dtset%nspinor 345 density(1,irealsp,ispinor,ispinorp) = & 346 & density(1,irealsp,ispinor,ispinorp) + & 347 & prefact_nk * (psi_r(1,i1,i2,i3,ispinor)*psi_r(1,i1,i2,i3,ispinorp)& 348 & + psi_r(2,i1,i2,i3,ispinor)*psi_r(2,i1,i2,i3,ispinorp)) 349 density(2,irealsp,ispinor,ispinorp) = & 350 & density(2,irealsp,ispinor,ispinorp) + & 351 & prefact_nk * (psi_r(1,i1,i2,i3,ispinor)*psi_r(2,i1,i2,i3,ispinorp)& 352 & - psi_r(2,i1,i2,i3,ispinor)*psi_r(1,i1,i2,i3,ispinorp)) 353 354 ! do i3p=1,dtset%ngfft(3) 355 ! do i2p=1,dtset%ngfft(2) 356 ! do i1p=1,dtset%ngfft(1) 357 ! irealsp_p = i1p + (i2p-1)*dtset%ngfft(1) + (i3p-1)*dtset%ngfft(2)*dtset%ngfft(1) 358 ! 359 ! NOTE : sign changes in second terms below because rho = psi*(r) psi(rprime) 360 ! 361 ! density_matrix(1,irealsp,ispinor,irealsp_p,ispinorp) = & 362 ! & density_matrix(1,irealsp,ispinor,irealsp_p,ispinorp) + & 363 ! & prefact_nk * (psi_r(1,i1,i2,i3,ispinor)*psi_r(1,i1p,i2p,i3p,ispinorp)& 364 ! & + psi_r(2,i1,i2,i3,ispinor)*psi_r(2,i1p,i2p,i3p,ispinorp)) 365 ! density_matrix(2,irealsp,ispinor,irealsp_p,ispinorp) = & 366 ! & density_matrix(2,irealsp,ispinor,irealsp_p,ispinorp) + & 367 ! & prefact_nk * (psi_r(1,i1,i2,i3,ispinor)*psi_r(2,i1p,i2p,i3p,ispinorp)& 368 ! & - psi_r(2,i1,i2,i3,ispinor)*psi_r(1,i1p,i2p,i3p,ispinorp)) 369 ! end do 370 ! end do 371 ! end do ! end irealspprime 372 373 end do !end ispinorp do 374 375 end do 376 end do 377 end do ! end irealsp 378 end do !end ispinor do 379 380 ! update pw counter 381 icg=icg+npw 382 iocc=iocc+1 383 end do ! iband 384 385 ikg=ikg+npw 386 387 ! deallocate arrays dep on npw for this kpoint 388 ABI_FREE(kg_k) 389 ABI_FREE(gpsi) 390 ABI_FREE(psi) 391 ABI_FREE(kgcart) 392 393 end do ! ikpt 394 395 ABI_FREE(dpsidr) 396 ABI_FREE(psi_r) 397 ABI_FREE(dummy_denpot) 398 ABI_FREE(gbound) 399 400 !prefactor for contribution to spin current 401 !prefactor is 1/2 * 1/2 * 2 Re(.): 402 !1/2 from the formula for the current 403 !1/2 from the use of the normalized Pauli matrices 404 !2 from the complex conjugate part 405 !total = 1/2 406 spincurrent = half * spincurrent 407 408 !make array of positions for all points on grid 409 ABI_MALLOC(position_op,(3,dtset%ngfft(1),dtset%ngfft(2),dtset%ngfft(3))) 410 do i3=1,dtset%ngfft(3) 411 do i2=1,dtset%ngfft(2) 412 do i1=1,dtset%ngfft(1) 413 position_op(:,i1,i2,i3) = matmul(hdr%rprimd,(/i1-one,i2-one,i3-one/))& 414 & /(/dtset%ngfft(1),dtset%ngfft(2),dtset%ngfft(3)/) 415 end do 416 end do 417 end do 418 419 !----------------------------------------------------------------------------------------- 420 !----------------------------------------------------------------------------------------- 421 !add electric field term to current. Non local term in case of pseudopotential SO 422 !present theory is that it is equal to A(r,r') = (W_SO(r,r') + W_SO(r',r)) 423 !For the strictly local part of the current, this becomes 2 W_SO(r,r) 424 ! 425 !W_SO is the prefactor in the spinorbit part of the potential, such that it 426 !can be written V_SO = W_SO . p (momentum operator) 427 !decomposed from V_SO = v_SO(r,r') L.S = v_SO(r,r') (rxp).S = v_SO(r,r') (Sxr).p 428 !and ensuring symmetrization for the r operator wrt the two arguments of v_SO(r,r') 429 !Hence: 430 !W_SO(r,r) = v_SO(r,r) (Sxr) 431 !----------------------------------------------------------------------------------------- 432 !----------------------------------------------------------------------------------------- 433 434 !allocate (vso_realspace_nl(2,dtset%ngfft(1)*dtset%ngfft(2)*dtset%ngfft(3),dtset%nspinor,& 435 !& dtset%ngfft(1)*dtset%ngfft(2)*dtset%ngfft(3),dtset%nspinor)) 436 437 !call vso_realspace_nonlop(atindx,atindx1,dtfil,dtset,gmet,gprimd,hdr,kg,& 438 !& mpi_enreg,nattyp,ph1d,position_op,psps,rmet,ucvol,vso_realspace_nl,ylm,ylmgr) 439 !anticommutator of VSO with position operator 440 !--- not needed in local spin current case --- 441 442 ABI_MALLOC(vso_realspace,(2,dtset%ngfft(1)*dtset%ngfft(2)*dtset%ngfft(3),dtset%nspinor,dtset%nspinor,3)) 443 444 call vso_realspace_local(dtset,hdr,position_op,psps,vso_realspace) 445 446 447 !multiply by density (or density matrix for nonlocal case) 448 !and add to spin current 449 450 451 452 ABI_FREE(density) 453 454 realrecip = 0 ! real space for xsf output 455 ABI_MALLOC(xcart,(3,dtset%natom)) 456 call xred2xcart(dtset%natom,hdr%rprimd,xcart,hdr%xred) 457 458 !----------------------------------------------------------------------------------------- 459 !----------------------------------------------------------------------------------------- 460 !output 3 components of current for each real space point 461 !----------------------------------------------------------------------------------------- 462 !----------------------------------------------------------------------------------------- 463 do ispindir=1, 3 464 ! choose rescale_current such that the maximum current component printed out 465 ! is 1 percent of lattice distance 466 ! By default XCrysDen multiplies by 200 to get something comparable to a distance in real space. 467 rescale_current = maxval(abs(spincurrent(:, :, :, :, ispindir))) 468 if (abs(rescale_current) < tol8) then 469 rescale_current = one 470 else 471 rescale_current = 0.001_dp * sqrt(max(sum(hdr%rprimd(:,1)**2), & 472 & sum(hdr%rprimd(:,2)**2), sum(hdr%rprimd(:,3)**2))) 473 end if 474 475 filnam=trim(dtfil%fnameabo_spcur)//spin_symbol(ispindir)//".xsf" 476 if (open_file(filnam,message,newunit=spcur_unit,status='unknown') /= 0) then 477 ABI_ERROR(message) 478 end if 479 480 ! print header 481 write (spcur_unit,'(a)') '#' 482 write (spcur_unit,'(a)') '# Xcrysden format file' 483 write (spcur_unit,'(a)') '# spin current density, for all real space points' 484 write (spcur_unit,'(a,3(I5,1x))') '# fft grid is ', dtset%ngfft(1), dtset%ngfft(2), dtset%ngfft(3) 485 write (spcur_unit,'(a,a,a)') '# ', spin_symbol(ispindir), '-spin current, full vector ' 486 487 write (spcur_unit,'(a)') 'ATOMS' 488 do iatom = 1, dtset%natom 489 write (spcur_unit,'(I4, 2x, 3(E16.6, 1x))') int(dtset%znucl(dtset%typat(iatom))), xcart(:,iatom) 490 end do 491 492 do i3=1,dtset%ngfft(3) 493 do i2=1,dtset%ngfft(2) 494 do i1=1,dtset%ngfft(1) 495 write (spcur_unit,'(a, 3(E10.3),2x, 3(E20.10))') 'X ', & 496 & position_op(:, i1, i2, i3), spincurrent(i1, i2, i3, :, ispindir)*rescale_current 497 end do 498 end do 499 end do 500 close (spcur_unit) 501 502 end do ! end ispindir 503 504 !----------------------------------------------------------------------------------------- 505 !----------------------------------------------------------------------------------------- 506 !output 3 spin components of V_SO matrices, for each real space point 507 !----------------------------------------------------------------------------------------- 508 !----------------------------------------------------------------------------------------- 509 ABI_MALLOC(datagrid, (dtset%ngfft(1)*dtset%ngfft(2)*dtset%ngfft(3))) 510 511 do ispindir=1,3 512 do icplex=1,2 513 do ispinor=1,dtset%nspinor 514 do ispinorp=1,dtset%nspinor 515 516 ! for the moment only print out if non zero 517 if (abs(sum(vso_realspace(icplex, :, ispinor, ispinorp, ispindir))) < tol8) cycle 518 519 filnam=trim(dtfil%fnameabo_vso)//"_spin_"//spin_symbol(ispindir)//"_"//& 520 & spinor_sym(ispinor)//spinor_sym(ispinorp)//"_"//realimag(icplex)//".xsf" 521 522 if (open_file(filnam,message,newunit=spcur_unit,status='unknown') /= 0) then 523 ABI_ERROR(message) 524 end if 525 526 ! print header 527 write (spcur_unit,'(a)') '#' 528 write (spcur_unit,'(a)') '# Xcrysden format file' 529 write (spcur_unit,'(a)') '# spin-orbit potential (space diagonal), for all real space points' 530 write (spcur_unit,'(a)') '# Real part first, then imaginary part' 531 write (spcur_unit,'(a,3(I5,1x))') '# fft grid is ', dtset%ngfft(1), dtset%ngfft(2), dtset%ngfft(3) 532 write (spcur_unit,'(a,a,a)') '# ', spin_symbol(ispindir), '-spin contribution ' 533 write (spcur_unit,'(a,a,a)') '# ', spinor_sym(ispinor)//spinor_sym(ispinorp), & 534 '-spin element of the spinor 2x2 matrix ' 535 write (spcur_unit,'(a,a)') '# cart x * cart y * cart z ***',& 536 & ' up-up component * up-down * down-up * down-down ' 537 538 ! Build contiguous array 539 datagrid = vso_realspace(icplex,:,ispinor, ispinorp, ispindir) 540 541 call printxsf(dtset%ngfft(1),dtset%ngfft(2),dtset%ngfft(3),& 542 & datagrid,hdr%rprimd,(/zero,zero,zero/), dtset%natom, dtset%ntypat, & 543 & dtset%typat, xcart, dtset%znucl, spcur_unit,realrecip) 544 545 ! 546 ! NOTE: have chosen actual dims of grid (n123) instead of fft box, for which n45 547 ! may be different 548 ! 549 ! do i3_dum=1,dtset%ngfft(3)+1 550 ! i3 = mod(i3_dum-1,dtset%ngfft(3)) + 1 551 ! do i2_dum=1,dtset%ngfft(2)+1 552 ! i2 = mod(i2_dum-1,dtset%ngfft(2)) + 1 553 ! do i1_dum=1,dtset%ngfft(1)+1 554 ! i1 = mod(i1_dum-1,dtset%ngfft(1)) + 1 555 ! 556 ! irealsp = i1 + (i2-1)*dtset%ngfft(1) + (i3-1)*dtset%ngfft(2)*dtset%ngfft(1) 557 ! write (spcur_unit,'(E20.10,1x)')& 558 ! & vso_realspace(icplex,irealsp, ispinor, ispinorp, ispindir) 559 ! end do 560 ! end do 561 ! end do 562 563 close (spcur_unit) 564 565 end do ! ispinorp 566 end do ! ispinor 567 end do ! icplex 568 end do ! end ispindir 569 570 ABI_FREE(datagrid) 571 ABI_FREE(vso_realspace) 572 !deallocate (vso_realspace_nl) 573 ABI_FREE(position_op) 574 ABI_FREE(spincurrent) 575 ABI_FREE(xcart) 576 577 write(std_out,*) ' Exiting subroutine spin_current ' 578 579 end subroutine spin_current
m_spin_current/vso_realspace_local [ Functions ]
[ Top ] [ m_spin_current ] [ Functions ]
NAME
vso_realspace_local
FUNCTION
Calculate real space (local - (r,r)) values of the SO part of the pseudopotential. Reconstructed explicitly in the HGH/GTH case.
INPUTS
OUTPUT
SIDE EFFECTS
NOTES
SOURCE
600 subroutine vso_realspace_local(dtset,hdr,position_op,psps,vso_realspace) 601 602 !Arguments ------------------------------- 603 type(hdr_type),intent(inout) :: hdr 604 type(dataset_type),intent(in) :: dtset 605 type(pseudopotential_type),intent(in) :: psps 606 real(dp),intent(in) :: position_op(3,dtset%ngfft(1),dtset%ngfft(2),dtset%ngfft(3)) 607 real(dp),intent(out) :: vso_realspace(2,dtset%ngfft(1)*dtset%ngfft(2)*dtset%ngfft(3),& 608 & dtset%nspinor,dtset%nspinor,3) 609 610 !Local variables ------------------------- 611 !scalars 612 integer :: i,j,l, lmax,ipsp,iatom, ir1,ir2,ir3 613 integer :: rcexponent,irealsp 614 integer :: nradgrid,iradgrid 615 real(dp) :: gammai, gammaj, relative_position(3), radial_cutoff, norm_rel_pos 616 real(dp) :: expfact,lfact, vso_interpol, x,y,z 617 !arrays 618 real(dp) :: xcart(3,dtset%natom),splint_x(1),splint_y(1) 619 real(dp), allocatable :: radial_grid(:) 620 real(dp), allocatable :: prefact_ijl(:,:,:,:),tmpvso(:),tmpvso_pp(:) 621 real(dp), allocatable :: vso_radial(:,:),vso_radial_pp(:,:),tmp_spline(:) 622 real(dp), allocatable :: offdiag_l_fact(:,:,:),kpar_matrix(:,:) 623 624 ! ********************************************************************* 625 626 !recalculate xcart (option = 1) 627 call xred2xcart(dtset%natom,hdr%rprimd,xcart,hdr%xred) 628 629 lmax = psps%mpsang-1 630 631 !content of gth pseudo type: 632 !These are {rloc, C(1...4)} coefficients for psppar(0, :, :) indices, 633 !Followed by the h coefficients for psppar(1:2, 1:, :) indices. 634 !size (0:2, 0:4, npsp) 635 !potential radius r_l is in psppar(l+1,0,ipsp) 636 !real(dp), pointer :: psppar(:, :, :) 637 !The covalence radii for each pseudo (?) size (npsp) 638 !real(dp), pointer :: radii_cov(:) 639 !Cut-off radii for core part and long-range part. 640 !radii_cf(:, 1) is for the long-range cut-off and 641 !radii_cf(:, 2) is for the core cut-off. 642 !size (npsp, 2) 643 !real(dp), pointer :: radii_cf(:, :) 644 !Spin orbit coefficients in HGH/GTH formats: k11p 645 !etc... see psp3ini.F90 646 !dimension = num l channels, 3 coeffs, num psp = 647 !(1:lmax+1,1:3,npsp) 648 !real(dp), pointer :: psp_k_par(:, :, :) 649 650 !v_SO^l (r,r') = sum_i sum_j sum_m Y_{lm} (\hat{r}) p_i^l (r) k_{ij}^l p_j^l(r') Y^{*}_lm (\hat{r'}) 651 ! 652 !v_SO^l (r,r) = sum_ij p_i^l (r) k_{ij}^l p_j^l(r) sum_m Y_{lm} (\hat{r}) Y^{*}_lm (\hat{r}) 653 != (2l+1)/4\pi sum_ij p_i^l (r) k_{ij}^l p_j^l(r) (eq B.17 Patrick Rinke thesis) 654 !p are gaussian projectors (from HGH paper prb 58 3641) [[cite:Hartwigsen1998]] 655 !sum_l v_SO^l (r,r) is a purely radial quantity (function of |r|), so spline it 656 657 !maximum distance needed in unit cell 658 radial_cutoff = four * maxval(psps%gth_params%psppar(:, 0, :)) 659 660 !setup radial grid; Should we use a logarithmic grid? The spline functions can 661 !take it... 662 nradgrid = 201 ! this is heuristic 663 ABI_MALLOC(radial_grid,(nradgrid)) 664 do iradgrid=1,nradgrid 665 radial_grid(iradgrid) = (iradgrid-1)*radial_cutoff/(nradgrid-1) 666 end do 667 668 !calculate prefactors independent of r 669 ABI_MALLOC(prefact_ijl,(3,3,0:lmax,psps%npsp)) 670 ABI_MALLOC(offdiag_l_fact,(3,3,0:lmax)) 671 ABI_MALLOC(kpar_matrix,(3,3)) 672 673 !these factors complete the full 3x3 matrix of k (or h) parameters for the 674 !HGH pseudos 675 offdiag_l_fact = zero 676 !l=0 677 offdiag_l_fact(1,2,0) = -half*sqrt(three/five) 678 offdiag_l_fact(1,3,0) = half*sqrt(five/21._dp) 679 offdiag_l_fact(2,3,0) = -half*sqrt(100._dp/63._dp) 680 !l=1 681 offdiag_l_fact(1,2,1) = -half*sqrt(five/seven) 682 offdiag_l_fact(1,3,1) = sixth*sqrt(35._dp/11._dp) 683 offdiag_l_fact(2,3,1) = -sixth*14._dp/sqrt(11._dp) 684 !l=2 685 if (lmax >= 2) then 686 offdiag_l_fact(1,2,2) = -half*sqrt(seven/nine) 687 offdiag_l_fact(1,3,2) = half*sqrt(63._dp/143._dp) 688 offdiag_l_fact(2,3,2) = -half*18._dp /sqrt(143._dp) 689 end if 690 !l=3 691 if (lmax >= 3) then 692 offdiag_l_fact(1,2,3) = zero 693 offdiag_l_fact(1,3,3) = zero 694 offdiag_l_fact(2,3,3) = zero 695 end if 696 !get prefactors for evaluation of V_SO: terms that do not depend on r 697 prefact_ijl = zero 698 do l=0,lmax 699 ! first the diagonal i=j term 700 do i=1,3 701 call gamma_function(l+(4._dp*i-1._dp)*0.5_dp, gammai) 702 gammai = sqrt(gammai) 703 rcexponent = 2*l+2*i+2*i-1 704 do ipsp=1,psps%npsp 705 prefact_ijl(i,i,l,ipsp) = psps%gth_params%psp_k_par(l+1,i,ipsp) & 706 & / ( (psps%gth_params%psppar(l+1,0,ipsp))**(rcexponent) & 707 & * gammai * gammai) 708 end do 709 end do 710 ! now the off diagonal elements 711 do ipsp=1,psps%npsp 712 kpar_matrix(1,2) = offdiag_l_fact (1,2,l)* psps%gth_params%psp_k_par(l+1,2,ipsp) 713 kpar_matrix(2,1) = kpar_matrix(1,2) 714 kpar_matrix(1,3) = offdiag_l_fact (1,3,l)* psps%gth_params%psp_k_par(l+1,3,ipsp) 715 kpar_matrix(3,1) = kpar_matrix(1,3) 716 kpar_matrix(2,3) = offdiag_l_fact (2,3,l)* psps%gth_params%psp_k_par(l+1,3,ipsp) 717 kpar_matrix(3,2) = kpar_matrix(2,3) 718 end do 719 720 ! for the f case only the 1,1 matrix element is non 0 - it is done above and 721 ! all these terms are actually 0 722 if (l > 2) cycle 723 724 do i=1,3 725 call gamma_function(l+(4._dp*i-1._dp)*0.5_dp, gammai) 726 gammai = sqrt(gammai) 727 do j=1,3 728 if (j==i) cycle 729 rcexponent = 2*l+2*i+2*j-1 730 call gamma_function(l+(4._dp*j-1._dp)*0.5_dp,gammaj) 731 gammaj = sqrt(gammaj) 732 do ipsp=1,psps%npsp 733 prefact_ijl(i,j,l,ipsp) = kpar_matrix(i,j) & 734 & / ( (psps%gth_params%psppar(l+1,0,ipsp))**rcexponent & 735 & * gammai * gammaj ) 736 end do 737 end do 738 end do 739 end do 740 741 ABI_FREE(kpar_matrix) 742 ABI_FREE(offdiag_l_fact) 743 744 prefact_ijl = prefact_ijl * two 745 746 !calculate v_SO on radial grid 747 ! MGNAG Runtime Error: *** Arithmetic exception: Floating invalid operation - aborting 748 ABI_MALLOC(vso_radial,(nradgrid,psps%npsp)) 749 vso_radial = zero 750 do l=0,lmax 751 lfact=(2._dp*l+1._dp)/four/pi 752 do iradgrid=1,nradgrid 753 norm_rel_pos = radial_grid(iradgrid) 754 do ipsp=1,psps%npsp 755 expfact = exp(-norm_rel_pos**2 / (psps%gth_params%psppar(l+1,0,ipsp))**2) 756 757 do i=1,3 758 do j=1,3 759 rcexponent = 2*l +2*i+2*j-4 760 if(prefact_ijl(i,j,l,ipsp)/=0) then !vz_d 0**0 761 vso_radial(iradgrid,ipsp) = vso_radial(iradgrid,ipsp) + & 762 & prefact_ijl(i,j,l,ipsp)*(norm_rel_pos**rcexponent) * expfact 763 end if !vz_d 764 end do ! j 765 end do ! i 766 end do ! ipsp 767 end do ! iradgrid 768 end do ! lmax 769 770 !spline v_SO(radial coord): get second derivative coefficients 771 ABI_MALLOC(vso_radial_pp,(nradgrid,psps%npsp)) 772 773 ABI_MALLOC(tmp_spline,(nradgrid)) 774 ABI_MALLOC(tmpvso,(nradgrid)) 775 ABI_MALLOC(tmpvso_pp,(nradgrid)) 776 do ipsp=1,psps%npsp 777 tmpvso = vso_radial(:,ipsp) 778 call spline( radial_grid, tmpvso, nradgrid, zero, radial_grid(nradgrid), tmpvso_pp ) 779 vso_radial_pp(:,ipsp) = tmpvso_pp 780 end do 781 ABI_FREE(tmp_spline) 782 ABI_FREE(tmpvso) 783 ABI_FREE(tmpvso_pp) 784 785 !to optimize this I should precalculate the distances which are actually needed by 786 !symmetry, or only sum over irreducible points in space and use weights 787 788 !for each physical atom present in unit cell 789 vso_realspace = zero 790 do iatom=1,dtset%natom 791 ! atom type will be dtset%typat(iatom) 792 793 ! for each point on grid 794 do ir3=1,dtset%ngfft(3) 795 do ir2=1,dtset%ngfft(2) 796 do ir1=1,dtset%ngfft(1) 797 irealsp = ir1 + (ir2-1)*dtset%ngfft(1) + (ir3-1)*dtset%ngfft(2)*dtset%ngfft(1) 798 799 ! relative position from atom to point 800 relative_position = position_op(:,ir1,ir2,ir3) - xcart(:,iatom) 801 x=relative_position(1) 802 y=relative_position(2) 803 z=relative_position(3) 804 805 ! calculate norm^2 806 norm_rel_pos = relative_position(1)**2+relative_position(2)**2+relative_position(3)**2 807 808 ! if norm^2 is too large, skip this point 809 if (norm_rel_pos > radial_cutoff*radial_cutoff) cycle 810 811 ! calculate norm 812 splint_x(1) = sqrt(norm_rel_pos) 813 814 ! spline interpolate vso only depends on position (through pos - atomic position) 815 call splint (nradgrid,radial_grid,vso_radial(:,dtset%typat(iatom)),& 816 & vso_radial_pp(:,dtset%typat(iatom)),1,splint_x,splint_y) 817 vso_interpol=splint_y(1) 818 819 ! multiply by vectorial spin factor (S x r) 820 ! NOTE: this r is taken relative to atom center. It could be that the r operator should 821 ! applied in an absolute way wrt the origin... 822 ! 823 ! Is this correct: accumulated sum over atoms ? 824 vso_realspace(1,irealsp,:,:,1) = vso_realspace(1,irealsp,:,:,1) + & 825 & vso_interpol * reshape((/y, zero,zero,-y/),(/2,2/)) 826 vso_realspace(2,irealsp,:,:,1) = vso_realspace(2,irealsp,:,:,1) + & 827 & vso_interpol * reshape((/zero,z, -z, zero/),(/2,2/)) 828 829 vso_realspace(1,irealsp,:,:,2) = vso_realspace(1,irealsp,:,:,2) + & 830 & vso_interpol * reshape((/-x, z, z, x/),(/2,2/)) 831 vso_realspace(2,irealsp,:,:,2) = vso_realspace(2,irealsp,:,:,2) + & 832 & vso_interpol * reshape((/zero,zero,zero,zero/),(/2,2/)) 833 834 vso_realspace(1,irealsp,:,:,3) = vso_realspace(1,irealsp,:,:,3) + & 835 & vso_interpol * reshape((/zero,-y, -y, zero/),(/2,2/)) 836 vso_realspace(2,irealsp,:,:,3) = vso_realspace(2,irealsp,:,:,3) + & 837 & vso_interpol * reshape((/zero,-x, -x, zero/),(/2,2/)) 838 839 end do ! ir3 840 end do ! ir2 841 end do ! ir1 842 end do ! iatom 843 844 ABI_FREE(prefact_ijl) 845 ABI_FREE(vso_radial) 846 ABI_FREE(vso_radial_pp) 847 ABI_FREE(radial_grid) 848 849 end subroutine vso_realspace_local