TABLE OF CONTENTS
ABINIT/spin_current [ Functions ]
NAME
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 .
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
SIDE EFFECTS
NOTES
PARENTS
afterscfloop
CHILDREN
fourwf,printxsf,sphereboundary,vso_realspace_local,xred2xcart
SOURCE
56 #if defined HAVE_CONFIG_H 57 #include "config.h" 58 #endif 59 60 #include "abi_common.h" 61 62 subroutine spin_current(cg,dtfil,dtset,gprimd,hdr,kg,mcg,mpi_enreg,psps) 63 64 use defs_basis 65 use defs_datatypes 66 use defs_abitypes 67 use m_errors 68 use m_profiling_abi 69 70 use m_io_tools, only : open_file 71 use m_pptools, only : printxsf 72 use m_geometry, only : xred2xcart 73 74 !This section has been created automatically by the script Abilint (TD). 75 !Do not modify the following lines by hand. 76 #undef ABI_FUNC 77 #define ABI_FUNC 'spin_current' 78 use interfaces_52_fft_mpi_noabirule 79 use interfaces_53_ffts 80 use interfaces_67_common, except_this_one => spin_current 81 !End of the abilint section 82 83 implicit none 84 85 !Arguments ------------------------------------ 86 !scalars 87 !integer,intent(in) :: nfftf 88 !real(dp),intent(in) :: ucvol 89 integer,intent(in) :: mcg 90 type(MPI_type),intent(in) :: mpi_enreg 91 type(datafiles_type),intent(in) :: dtfil 92 type(dataset_type),intent(in) :: dtset 93 type(hdr_type),intent(inout) :: hdr 94 type(pseudopotential_type),intent(in) :: psps 95 !type(wffile_type),intent(in) :: wffnow 96 !arrays 97 !integer,intent(in) :: atindx(dtset%natom),atindx1(dtset%natom) 98 integer,intent(in) :: kg(3,dtset%mpw*dtset%mkmem) 99 !integer,intent(in) :: nattyp(dtset%ntypat) 100 !integer,intent(in) :: symrec(3,3,dtset%nsym) 101 real(dp),intent(in) :: cg(2,mcg) 102 !real(dp),intent(in) :: eigen(dtset%mband*dtset%nkpt*dtset%nsppol),gmet(3,3) 103 real(dp),intent(in) :: gprimd(3,3) 104 !real(dp),intent(in) :: rhog(2,nfftf),rhor(nfftf,dtset%nspden) 105 !real(dp),intent(in) :: rmet(3,3) 106 !real(dp),intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm) 107 !real(dp),intent(in) :: ylmgr(dtset%mpw*dtset%mkmem,3,psps%mpsang*psps%mpsang*psps%useylm) 108 !real(dp),intent(inout) :: ph1d(2,3*(2*dtset%mgfft+1)*dtset%natom) 109 110 !Local variables------------------------------- 111 !scalars 112 integer :: cplex,fft_option,i1 113 integer :: i2,i3,iband,icartdir,icg,ig 114 integer :: ikg,ikpt,iocc,irealsp,ispindir,ispinor,ispinorp 115 integer :: npw 116 integer :: icplex 117 integer :: realrecip 118 integer :: iatom,spcur_unit 119 real(dp) :: prefact_nk 120 real(dp) :: rescale_current 121 character(len=500) :: message 122 character(len=fnlen) :: filnam 123 !arrays 124 integer,allocatable :: gbound(:,:),kg_k(:,:) 125 real(dp),allocatable :: dpsidr(:,:,:,:,:,:) 126 real(dp),allocatable :: density(:,:,:,:) 127 real(dp),allocatable :: dummy_denpot(:,:,:) 128 real(dp),allocatable :: gpsi(:,:,:,:),kgcart(:,:) 129 real(dp),allocatable :: position_op(:,:,:,:) 130 real(dp),allocatable :: psi(:,:,:),psi_r(:,:,:,:,:) 131 real(dp),allocatable :: spincurrent(:,:,:,:,:) 132 real(dp),allocatable :: vso_realspace(:,:,:,:,:),datagrid(:) 133 real(dp) :: dummy_fofgout(0,0) 134 real(dp),allocatable :: xcart(:,:) 135 character :: spin_symbol(3) 136 character :: spinor_sym(2) 137 character(len=2) :: realimag(2) 138 !no_abirules 139 !real(dp),allocatable :: density_matrix(:,:,:,:,:) 140 !real(dp),allocatable :: vso_realspace_nl(:,:,:,:,:) 141 142 ! ************************************************************************* 143 !source 144 145 !write(std_out,*) ' Entering subroutine spin_current ' 146 !write(std_out,*) ' dtset%ngfft = ', dtset%ngfft 147 !write(std_out,*) ' hdr%istwfk = ', hdr%istwfk 148 149 !===================== init and checks ============================ 150 !check if nspinor is 2 151 if (dtset%nspinor /= 2) then 152 write(message, '(a,i0)' )' nspinor must be 2, but it is ',dtset%nspinor 153 MSG_ERROR(message) 154 end if 155 156 if (dtset%nsppol /= 1) then 157 write(message, '(a,i0)' )' spin_current: nsppol must be 1 but it is ',dtset%nsppol 158 MSG_ERROR(message) 159 end if 160 161 if (dtset%mkmem /= dtset%nkpt) then 162 write(message, '(a,i6,a,i6,a,a)' )& 163 & ' mkmem = ',dtset%mkmem,' must be equal to nkpt ',dtset%nkpt,ch10,& 164 & ' keep all kpt in memory' 165 MSG_ERROR(message) 166 end if 167 168 if (dtset%usepaw /= 0) then 169 write(message, '(a,i0,a,a,a)' )& 170 & 'usepaw = ',dtset%usepaw,' must be equal to 0 ',ch10,& 171 & 'Not functional for PAW case yet.' 172 MSG_ERROR(message) 173 end if 174 175 cplex=2 176 fft_option = 0 ! just do direct fft 177 spin_symbol = (/'x','y','z'/) 178 spinor_sym = (/'u','d'/) 179 realimag = (/'Re','Im'/) 180 181 182 write(std_out,*) ' psps%mpsang,psps%mpssoang ', psps%mpsang,psps%mpssoang 183 184 !======================= main code ================================ 185 !----------------------------------------------------------------------------------------- 186 !----------------------------------------------------------------------------------------- 187 !first get normal contribution to current, as psi tau dpsidr + dpsidr tau psi 188 !where tau are 1/2 the pauli matrices 189 !----------------------------------------------------------------------------------------- 190 !----------------------------------------------------------------------------------------- 191 192 !init plane wave coeff counter 193 icg = 0 194 !init plane wave counter 195 ikg = 0 196 !init occupation/band counter 197 iocc = 1 198 199 !rspace point, cartesian direction, spin pol=x,y,z 200 ABI_ALLOCATE(spincurrent,(dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6),3,3)) 201 spincurrent = zero 202 203 ABI_ALLOCATE(dummy_denpot,(cplex*dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6))) 204 205 ABI_ALLOCATE(gbound,(2*dtset%mgfft+8,2)) 206 207 !allocate (density_matrix(2,dtset%ngfft(1)*dtset%ngfft(2)*dtset%ngfft(3),dtset%nspinor,& 208 !& dtset%ngfft(1)*dtset%ngfft(2)*dtset%ngfft(3),dtset%nspinor)) 209 !density_matrix= zero 210 ABI_ALLOCATE(density,(2,dtset%ngfft(1)*dtset%ngfft(2)*dtset%ngfft(3),dtset%nspinor,dtset%nspinor)) 211 density= zero 212 213 ABI_ALLOCATE(dpsidr,(2,dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6),dtset%nspinor,3)) 214 ABI_ALLOCATE(psi_r,(2,dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6),dtset%nspinor)) 215 216 !loop over kpoints 217 do ikpt=1,dtset%nkpt 218 219 220 ! number of plane waves for this kpt 221 npw = hdr%npwarr(ikpt) 222 223 ! allocate arrays dep on number of pw 224 ABI_ALLOCATE(kg_k,(3,npw)) 225 ABI_ALLOCATE(gpsi,(2,npw,dtset%nspinor,3)) 226 ABI_ALLOCATE(psi,(2,npw,dtset%nspinor)) 227 ABI_ALLOCATE(kgcart,(3,npw)) 228 229 ! get cartesian coordinates of k+G vectors around this kpoint 230 do ig=1,npw 231 kgcart(:,ig) = matmul(gprimd(:,:),dtset%kpt(:,ikpt)+kg(:,ikg+ig)) 232 kg_k (:,ig) = kg(:,ikg+ig) 233 end do 234 235 ! get gbound 236 call sphereboundary(gbound,dtset%istwfk(ikpt),kg_k,dtset%mgfft,npw) 237 238 ! loop over bands 239 do iband=1,dtset%nband(ikpt) 240 241 ! prefactor for sum over bands and kpoints 242 prefact_nk = hdr%occ(iocc) * dtset%wtk(ikpt) 243 244 ! initialize this wf 245 gpsi=zero 246 psi=zero 247 psi(:,1:npw,1) = cg(:,icg+1:icg+npw) 248 249 ! multiply psi by - i 2 pi G 250 do ig=1,npw 251 gpsi(1,ig,:,1) = two_pi * kgcart(1,ig)*psi(2,ig,:) 252 gpsi(2,ig,:,1) = -two_pi * kgcart(1,ig)*psi(1,ig,:) 253 gpsi(1,ig,:,2) = two_pi * kgcart(2,ig)*psi(2,ig,:) 254 gpsi(2,ig,:,2) = -two_pi * kgcart(2,ig)*psi(1,ig,:) 255 gpsi(1,ig,:,3) = two_pi * kgcart(3,ig)*psi(2,ig,:) 256 gpsi(2,ig,:,3) = -two_pi * kgcart(3,ig)*psi(1,ig,:) 257 end do 258 259 ! loop over spinorial components 260 do ispinor=1,dtset%nspinor 261 ! FT Gpsi_x to real space 262 call fourwf(cplex,dummy_denpot,gpsi(:,:,ispinor,1),dummy_fofgout,& 263 & dpsidr(:,:,:,:,ispinor,1),gbound,gbound,& 264 & hdr%istwfk(ikpt),kg_k,kg_k,dtset%mgfft,mpi_enreg,1,dtset%ngfft,npw,& 265 & npw,dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6),& 266 & fft_option,dtset%paral_kgb,0,one,one,use_gpu_cuda=dtset%use_gpu_cuda) 267 268 ! FT Gpsi_y to real space 269 call fourwf(cplex,dummy_denpot,gpsi(:,:,ispinor,2),dummy_fofgout,& 270 & dpsidr(:,:,:,:,ispinor,2),gbound,gbound,& 271 & hdr%istwfk(ikpt),kg_k,kg_k,dtset%mgfft,mpi_enreg,1,dtset%ngfft,npw,& 272 & npw,dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6),& 273 & fft_option,dtset%paral_kgb,0,one,one,use_gpu_cuda=dtset%use_gpu_cuda) 274 275 ! FT Gpsi_z to real space 276 call fourwf(cplex,dummy_denpot,gpsi(:,:,ispinor,3),dummy_fofgout,& 277 & dpsidr(:,:,:,:,ispinor,3),gbound,gbound,& 278 & hdr%istwfk(ikpt),kg_k,kg_k,dtset%mgfft,mpi_enreg,1,dtset%ngfft,npw,& 279 & npw,dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6),& 280 & fft_option,dtset%paral_kgb,0,one,one,use_gpu_cuda=dtset%use_gpu_cuda) 281 282 ! FT psi to real space 283 call fourwf(cplex,dummy_denpot,psi(:,:,ispinor),dummy_fofgout,& 284 & psi_r(:,:,:,:,ispinor),gbound,gbound,& 285 & hdr%istwfk(ikpt),kg_k,kg_k,dtset%mgfft,mpi_enreg,1,dtset%ngfft,npw,& 286 & npw,dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6),& 287 & fft_option,dtset%paral_kgb,0,one,one,use_gpu_cuda=dtset%use_gpu_cuda) 288 289 end do ! ispinor 290 291 ! dpsidr now contains the full derivative of psi wrt space (gradient) in cartesian coordinates 292 293 ! get 3 pauli matrix contributions to the current: x,y,z, cart dir, spin dir 294 do icartdir=1,3 295 296 ! x pauli spin matrix 297 ! sigma_x = | 0 1 | 298 ! | 1 0 | 299 spincurrent(:,:,:,icartdir,1) = spincurrent(:,:,:,icartdir,1) + prefact_nk * & 300 ! Re(psi_r(up)^* dpsidr(down)) 301 & real(psi_r(1,:,:,:,1)*dpsidr(1,:,:,:,2,icartdir) & 302 & + psi_r(2,:,:,:,1)*dpsidr(2,:,:,:,2,icartdir) & 303 ! Re(psi_r(down)^* dpsidr(up)) 304 & + psi_r(1,:,:,:,2)*dpsidr(1,:,:,:,1,icartdir) & 305 & + psi_r(2,:,:,:,2)*dpsidr(2,:,:,:,1,icartdir)) 306 307 ! y pauli spin matrix 308 ! sigma_y = | 0 -i | 309 ! | i 0 | 310 spincurrent(:,:,:,icartdir,2) = spincurrent(:,:,:,icartdir,2) + prefact_nk * & 311 ! Re(-i psi_r(up)^* dpsidr(down)) 312 & real(psi_r(1,:,:,:,1)*dpsidr(2,:,:,:,2,icartdir) & 313 & - psi_r(2,:,:,:,1)*dpsidr(1,:,:,:,2,icartdir) & 314 ! Re(i psi_r(down)^* dpsidr(up)) 315 & - psi_r(1,:,:,:,2)*dpsidr(2,:,:,:,1,icartdir) & 316 & + psi_r(2,:,:,:,2)*dpsidr(1,:,:,:,1,icartdir)) 317 318 ! z pauli spin matrix 319 ! sigma_z = | 1 0 | 320 ! | 0 -1 | 321 spincurrent(:,:,:,icartdir,3) = spincurrent(:,:,:,icartdir,3) + prefact_nk * & 322 ! Re(psi_r(up)^* dpsidr(up)) 323 & real(psi_r(1,:,:,:,1)*dpsidr(1,:,:,:,1,icartdir) & 324 & - psi_r(2,:,:,:,1)*dpsidr(2,:,:,:,1,icartdir) & 325 ! Re(-psi_r(down)^* dpsidr(down)) 326 & - psi_r(1,:,:,:,2)*dpsidr(1,:,:,:,2,icartdir) & 327 & + psi_r(2,:,:,:,2)*dpsidr(2,:,:,:,2,icartdir)) 328 end do ! end icartdir 329 330 ! 331 ! accumulate non local density matrix in real space 332 ! NOTE: if we are only using the local part of the current, this becomes the 333 ! density spinor matrix! (much lighter to calculate) rho(r, sigma, sigmaprime) 334 ! 335 do ispinor=1,dtset%nspinor 336 do i3=1,dtset%ngfft(3) 337 do i2=1,dtset%ngfft(2) 338 do i1=1,dtset%ngfft(1) 339 irealsp = i1 + (i2-1)*dtset%ngfft(1) + (i3-1)*dtset%ngfft(2)*dtset%ngfft(1) 340 341 do ispinorp=1,dtset%nspinor 342 density(1,irealsp,ispinor,ispinorp) = & 343 & density(1,irealsp,ispinor,ispinorp) + & 344 & prefact_nk * (psi_r(1,i1,i2,i3,ispinor)*psi_r(1,i1,i2,i3,ispinorp)& 345 & + psi_r(2,i1,i2,i3,ispinor)*psi_r(2,i1,i2,i3,ispinorp)) 346 density(2,irealsp,ispinor,ispinorp) = & 347 & density(2,irealsp,ispinor,ispinorp) + & 348 & prefact_nk * (psi_r(1,i1,i2,i3,ispinor)*psi_r(2,i1,i2,i3,ispinorp)& 349 & - psi_r(2,i1,i2,i3,ispinor)*psi_r(1,i1,i2,i3,ispinorp)) 350 351 ! do i3p=1,dtset%ngfft(3) 352 ! do i2p=1,dtset%ngfft(2) 353 ! do i1p=1,dtset%ngfft(1) 354 ! irealsp_p = i1p + (i2p-1)*dtset%ngfft(1) + (i3p-1)*dtset%ngfft(2)*dtset%ngfft(1) 355 ! 356 ! NOTE : sign changes in second terms below because rho = psi*(r) psi(rprime) 357 ! 358 ! density_matrix(1,irealsp,ispinor,irealsp_p,ispinorp) = & 359 ! & density_matrix(1,irealsp,ispinor,irealsp_p,ispinorp) + & 360 ! & prefact_nk * (psi_r(1,i1,i2,i3,ispinor)*psi_r(1,i1p,i2p,i3p,ispinorp)& 361 ! & + psi_r(2,i1,i2,i3,ispinor)*psi_r(2,i1p,i2p,i3p,ispinorp)) 362 ! density_matrix(2,irealsp,ispinor,irealsp_p,ispinorp) = & 363 ! & density_matrix(2,irealsp,ispinor,irealsp_p,ispinorp) + & 364 ! & prefact_nk * (psi_r(1,i1,i2,i3,ispinor)*psi_r(2,i1p,i2p,i3p,ispinorp)& 365 ! & - psi_r(2,i1,i2,i3,ispinor)*psi_r(1,i1p,i2p,i3p,ispinorp)) 366 ! end do 367 ! end do 368 ! end do ! end irealspprime 369 370 end do !end ispinorp do 371 372 end do 373 end do 374 end do ! end irealsp 375 end do !end ispinor do 376 377 ! update pw counter 378 icg=icg+npw 379 iocc=iocc+1 380 end do ! iband 381 382 ikg=ikg+npw 383 384 ! deallocate arrays dep on npw for this kpoint 385 ABI_DEALLOCATE(kg_k) 386 ABI_DEALLOCATE(gpsi) 387 ABI_DEALLOCATE(psi) 388 ABI_DEALLOCATE(kgcart) 389 390 end do ! ikpt 391 392 ABI_DEALLOCATE(dpsidr) 393 ABI_DEALLOCATE(psi_r) 394 ABI_DEALLOCATE(dummy_denpot) 395 ABI_DEALLOCATE(gbound) 396 397 !prefactor for contribution to spin current 398 !prefactor is 1/2 * 1/2 * 2 Re(.): 399 !1/2 from the formula for the current 400 !1/2 from the use of the normalized Pauli matrices 401 !2 from the complex conjugate part 402 !total = 1/2 403 spincurrent = half * spincurrent 404 405 !make array of positions for all points on grid 406 ABI_ALLOCATE(position_op,(3,dtset%ngfft(1),dtset%ngfft(2),dtset%ngfft(3))) 407 do i3=1,dtset%ngfft(3) 408 do i2=1,dtset%ngfft(2) 409 do i1=1,dtset%ngfft(1) 410 position_op(:,i1,i2,i3) = matmul(hdr%rprimd,(/i1-one,i2-one,i3-one/))& 411 & /(/dtset%ngfft(1),dtset%ngfft(2),dtset%ngfft(3)/) 412 end do 413 end do 414 end do 415 416 !----------------------------------------------------------------------------------------- 417 !----------------------------------------------------------------------------------------- 418 !add electric field term to current. Non local term in case of pseudopotential SO 419 !present theory is that it is equal to A(r,r') = (W_SO(r,r') + W_SO(r',r)) 420 !For the strictly local part of the current, this becomes 2 W_SO(r,r) 421 ! 422 !W_SO is the prefactor in the spinorbit part of the potential, such that it 423 !can be written V_SO = W_SO . p (momentum operator) 424 !decomposed from V_SO = v_SO(r,r') L.S = v_SO(r,r') (rxp).S = v_SO(r,r') (Sxr).p 425 !and ensuring symmetrization for the r operator wrt the two arguments of v_SO(r,r') 426 !Hence: 427 !W_SO(r,r) = v_SO(r,r) (Sxr) 428 !----------------------------------------------------------------------------------------- 429 !----------------------------------------------------------------------------------------- 430 431 !allocate (vso_realspace_nl(2,dtset%ngfft(1)*dtset%ngfft(2)*dtset%ngfft(3),dtset%nspinor,& 432 !& dtset%ngfft(1)*dtset%ngfft(2)*dtset%ngfft(3),dtset%nspinor)) 433 434 !call vso_realspace_nonlop(atindx,atindx1,dtfil,dtset,gmet,gprimd,hdr,kg,& 435 !& mpi_enreg,nattyp,ph1d,position_op,psps,rmet,ucvol,vso_realspace_nl,ylm,ylmgr) 436 !anticommutator of VSO with position operator 437 !--- not needed in local spin current case --- 438 439 ABI_ALLOCATE(vso_realspace,(2,dtset%ngfft(1)*dtset%ngfft(2)*dtset%ngfft(3),dtset%nspinor,dtset%nspinor,3)) 440 441 call vso_realspace_local(dtset,hdr,position_op,psps,vso_realspace) 442 443 444 !multiply by density (or density matrix for nonlocal case) 445 !and add to spin current 446 447 448 449 ABI_DEALLOCATE(density) 450 451 realrecip = 0 ! real space for xsf output 452 ABI_ALLOCATE(xcart,(3,dtset%natom)) 453 call xred2xcart(dtset%natom,hdr%rprimd,xcart,hdr%xred) 454 455 !----------------------------------------------------------------------------------------- 456 !----------------------------------------------------------------------------------------- 457 !output 3 components of current for each real space point 458 !----------------------------------------------------------------------------------------- 459 !----------------------------------------------------------------------------------------- 460 do ispindir=1, 3 461 ! choose rescale_current such that the maximum current component printed out 462 ! is 1 percent of lattice distance 463 ! By default XCrysDen multiplies by 200 to get something comparable to a distance in real space. 464 rescale_current = maxval(abs(spincurrent(:, :, :, :, ispindir))) 465 if (abs(rescale_current) < tol8) then 466 rescale_current = one 467 else 468 rescale_current = 0.001_dp * sqrt(max(sum(hdr%rprimd(:,1)**2), & 469 & sum(hdr%rprimd(:,2)**2), sum(hdr%rprimd(:,3)**2))) 470 end if 471 472 filnam=trim(dtfil%fnameabo_spcur)//spin_symbol(ispindir)//".xsf" 473 if (open_file(filnam,message,newunit=spcur_unit,status='unknown') /= 0) then 474 MSG_ERROR(message) 475 end if 476 477 ! print header 478 write (spcur_unit,'(a)') '#' 479 write (spcur_unit,'(a)') '# Xcrysden format file' 480 write (spcur_unit,'(a)') '# spin current density, for all real space points' 481 write (spcur_unit,'(a,3(I5,1x))') '# fft grid is ', dtset%ngfft(1), dtset%ngfft(2), dtset%ngfft(3) 482 write (spcur_unit,'(a,a,a)') '# ', spin_symbol(ispindir), '-spin current, full vector ' 483 484 write (spcur_unit,'(a)') 'ATOMS' 485 do iatom = 1, dtset%natom 486 write (spcur_unit,'(I4, 2x, 3(E16.6, 1x))') int(dtset%znucl(dtset%typat(iatom))), xcart(:,iatom) 487 end do 488 489 do i3=1,dtset%ngfft(3) 490 do i2=1,dtset%ngfft(2) 491 do i1=1,dtset%ngfft(1) 492 write (spcur_unit,'(a, 3(E10.3),2x, 3(E20.10))') 'X ', & 493 & position_op(:, i1, i2, i3), spincurrent(i1, i2, i3, :, ispindir)*rescale_current 494 end do 495 end do 496 end do 497 close (spcur_unit) 498 499 end do ! end ispindir 500 501 !----------------------------------------------------------------------------------------- 502 !----------------------------------------------------------------------------------------- 503 !output 3 spin components of V_SO matrices, for each real space point 504 !----------------------------------------------------------------------------------------- 505 !----------------------------------------------------------------------------------------- 506 ABI_MALLOC(datagrid, (dtset%ngfft(1)*dtset%ngfft(2)*dtset%ngfft(3))) 507 508 do ispindir=1,3 509 do icplex=1,2 510 do ispinor=1,dtset%nspinor 511 do ispinorp=1,dtset%nspinor 512 513 ! for the moment only print out if non zero 514 if (abs(sum(vso_realspace(icplex, :, ispinor, ispinorp, ispindir))) < tol8) cycle 515 516 filnam=trim(dtfil%fnameabo_vso)//"_spin_"//spin_symbol(ispindir)//"_"//& 517 & spinor_sym(ispinor)//spinor_sym(ispinorp)//"_"//realimag(icplex)//".xsf" 518 519 if (open_file(filnam,message,newunit=spcur_unit,status='unknown') /= 0) then 520 MSG_ERROR(message) 521 end if 522 523 ! print header 524 write (spcur_unit,'(a)') '#' 525 write (spcur_unit,'(a)') '# Xcrysden format file' 526 write (spcur_unit,'(a)') '# spin-orbit potential (space diagonal), for all real space points' 527 write (spcur_unit,'(a)') '# Real part first, then imaginary part' 528 write (spcur_unit,'(a,3(I5,1x))') '# fft grid is ', dtset%ngfft(1), dtset%ngfft(2), dtset%ngfft(3) 529 write (spcur_unit,'(a,a,a)') '# ', spin_symbol(ispindir), '-spin contribution ' 530 write (spcur_unit,'(a,a,a)') '# ', spinor_sym(ispinor)//spinor_sym(ispinorp), '-spin element of the spinor 2x2 matrix ' 531 532 write (spcur_unit,'(a,a)') '# cart x * cart y * cart z ***',& 533 & ' up-up component * up-down * down-up * down-down ' 534 535 ! Build contiguous array 536 datagrid = vso_realspace(icplex,:,ispinor, ispinorp, ispindir) 537 538 call printxsf(dtset%ngfft(1),dtset%ngfft(2),dtset%ngfft(3),& 539 & datagrid,hdr%rprimd,(/zero,zero,zero/), dtset%natom, dtset%ntypat, & 540 & dtset%typat, xcart, dtset%znucl, spcur_unit,realrecip) 541 542 ! 543 ! NOTE: have chosen actual dims of grid (n123) instead of fft box, for which n45 544 ! may be different 545 ! 546 ! do i3_dum=1,dtset%ngfft(3)+1 547 ! i3 = mod(i3_dum-1,dtset%ngfft(3)) + 1 548 ! do i2_dum=1,dtset%ngfft(2)+1 549 ! i2 = mod(i2_dum-1,dtset%ngfft(2)) + 1 550 ! do i1_dum=1,dtset%ngfft(1)+1 551 ! i1 = mod(i1_dum-1,dtset%ngfft(1)) + 1 552 ! 553 ! irealsp = i1 + (i2-1)*dtset%ngfft(1) + (i3-1)*dtset%ngfft(2)*dtset%ngfft(1) 554 ! write (spcur_unit,'(E20.10,1x)')& 555 ! & vso_realspace(icplex,irealsp, ispinor, ispinorp, ispindir) 556 ! end do 557 ! end do 558 ! end do 559 560 close (spcur_unit) 561 562 end do ! ispinorp 563 end do ! ispinor 564 end do ! icplex 565 end do ! end ispindir 566 567 ABI_FREE(datagrid) 568 ABI_DEALLOCATE(vso_realspace) 569 !deallocate (vso_realspace_nl) 570 ABI_DEALLOCATE(position_op) 571 ABI_DEALLOCATE(spincurrent) 572 ABI_DEALLOCATE(xcart) 573 574 write(std_out,*) ' Exiting subroutine spin_current ' 575 576 end subroutine spin_current