TABLE OF CONTENTS


ABINIT/m_spin_current [ Modules ]

[ Top ] [ 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