TABLE OF CONTENTS


ABINIT/m_spin_current [ Modules ]

[ Top ] [ Modules ]

NAME

  m_spin_current

FUNCTION

COPYRIGHT

 Copyright (C) 2005-2022 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,use_gpu_cuda=dtset%use_gpu_cuda)
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,use_gpu_cuda=dtset%use_gpu_cuda)
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,use_gpu_cuda=dtset%use_gpu_cuda)
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,use_gpu_cuda=dtset%use_gpu_cuda)
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