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