Branch data Line data Source code
1 : : !!****m* ABINIT/m_spin_current
2 : : !! NAME
3 : : !! m_spin_current
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 2005-2026 ABINIT group (Mver)
10 : : !! This file is distributed under the terms of the
11 : : !! GNU General Public License, see ~abinit/COPYING
12 : : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : : !!
14 : : !! SOURCE
15 : :
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
43 : : !!***
44 : :
45 : : public :: spin_current
46 : : !!***
47 : :
48 : : contains
49 : : !!***
50 : :
51 : : !!****f* m_spin_current/spin_current
52 : : !! NAME
53 : : !! spin_current
54 : : !!
55 : : !! FUNCTION
56 : : !!
57 : : !! INPUTS
58 : : !! atindx(natom)=index table for atoms (see gstate.f)
59 : : !! atindx1(natom)=inverse of atindx
60 : : !! cg(2,mcg)=wavefunctions (may be read from disk instead of input)
61 : : !! dtfil <type(datafiles_type)>=variables related to files
62 : : !! dtset <type(dataset_type)>=all input variables in this dataset
63 : : !! eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree)
64 : : !! gmet = reciprocal space metric
65 : : !! gprimd = dimensionful reciprocal space vectors
66 : : !! hdr <type(hdr_type)>=the header of wf, den and pot files
67 : : !! kg(3,mpw*mkmem)=reduced (integer) coordinates of G vecs in basis sphere
68 : : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
69 : : !! mpi_enreg=information about MPI parallelization
70 : : !! nattyp(dtset%ntypat)=number of atoms of each type
71 : : !! nfftf = fft grid dimensions for fine grid
72 : : !! ph1d = phase factors in 1 radial dimension
73 : : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
74 : : !! | mpsang= 1+maximum angular momentum
75 : : !! rhog(2,nfftf)=Fourier transform of total electron density (including compensation density in PAW)
76 : : !! rhor(nfftf,nspden)=total electron density (including compensation density in PAW)
77 : : !! rmet = real space metric tensor
78 : : !! symrec(3,3,nsym)=symmetries in reciprocal space, reduced coordinates
79 : : !! ucvol = unit cell volume
80 : : !! wffnow=unit number for current wf disk file
81 : : !! ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point
82 : : !! ylmgr(mpw*mkmem,3,mpsang*mpsang*useylm)= gradients of real spherical harmonics
83 : : !!
84 : : !! OUTPUT
85 : : !! only output to file
86 : : !!
87 : : !! SOURCE
88 : :
89 : 1 : 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 : 1 : integer,allocatable :: gbound(:,:),kg_k(:,:)
131 : 1 : real(dp),allocatable :: dpsidr(:,:,:,:,:,:)
132 : 1 : real(dp),allocatable :: density(:,:,:,:)
133 : 1 : real(dp),allocatable :: dummy_denpot(:,:,:)
134 : 1 : real(dp),allocatable :: gpsi(:,:,:,:),kgcart(:,:)
135 : 1 : real(dp),allocatable :: position_op(:,:,:,:)
136 : 1 : real(dp),allocatable :: psi(:,:,:),psi_r(:,:,:,:,:)
137 : 1 : real(dp),allocatable :: spincurrent(:,:,:,:,:)
138 : 1 : real(dp),allocatable :: vso_realspace(:,:,:,:,:),datagrid(:)
139 : : real(dp) :: dummy_fofgout(0,0)
140 : 1 : 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 [ - + ]: 1 : if (dtset%nspinor /= 2) then
157 : 0 : write(message, '(a,i0)' )' nspinor must be 2, but it is ',dtset%nspinor
158 : 0 : ABI_ERROR(message)
159 : : end if
160 : :
161 [ - + ]: 1 : if (dtset%nsppol /= 1) then
162 : 0 : write(message, '(a,i0)' )' spin_current: nsppol must be 1 but it is ',dtset%nsppol
163 : 0 : ABI_ERROR(message)
164 : : end if
165 : :
166 [ - + ]: 1 : if (dtset%mkmem /= dtset%nkpt) then
167 : : write(message, '(a,i6,a,i6,a,a)' )&
168 : 0 : & ' mkmem = ',dtset%mkmem,' must be equal to nkpt ',dtset%nkpt,ch10,&
169 : 0 : & ' keep all kpt in memory'
170 : 0 : ABI_ERROR(message)
171 : : end if
172 : :
173 [ - + ]: 1 : if (dtset%usepaw /= 0) then
174 : : write(message, '(a,i0,a,a,a)' )&
175 : 0 : & 'usepaw = ',dtset%usepaw,' must be equal to 0 ',ch10,&
176 : 0 : & 'Not functional for PAW case yet.'
177 : 0 : ABI_ERROR(message)
178 : : end if
179 : :
180 : 1 : cplex=2
181 : 1 : fft_option = 0 ! just do direct fft
182 [ + + ]: 4 : spin_symbol = (/'x','y','z'/)
183 [ + + ]: 3 : spinor_sym = (/'u','d'/)
184 [ + + ]: 3 : realimag = (/'Re','Im'/)
185 : :
186 : 1 : 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 : 1 : icg = 0
198 : : !init plane wave counter
199 : 1 : ikg = 0
200 : : !init occupation/band counter
201 : 1 : iocc = 1
202 : :
203 : : !rspace point, cartesian direction, spin pol=x,y,z
204 [ - + + - : 7 : ABI_MALLOC(spincurrent,(dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6),3,3))
- + + - -
+ + - + -
+ - + - +
- - + -
+ ]
205 [ + + + + : 32548 : spincurrent = zero
+ + + + +
+ ]
206 : :
207 [ - + + - : 5 : ABI_MALLOC(dummy_denpot,(cplex*dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6)))
- + + - -
+ + - + -
+ - - + -
+ ]
208 : :
209 [ + - - + : 4 : 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 [ - + + - : 5 : ABI_MALLOC(density,(2,dtset%ngfft(1)*dtset%ngfft(2)*dtset%ngfft(3),dtset%nspinor,dtset%nspinor))
- + + - -
+ + - + -
+ - - + -
+ ]
215 [ + + + + : 40507 : density= zero
+ + + + ]
216 : :
217 [ - + + - : 7 : ABI_MALLOC(dpsidr,(2,dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6),dtset%nspinor,3))
- + + - -
+ + - - +
+ - + - +
- + - - +
- + ]
218 [ - + + - : 6 : ABI_MALLOC(psi_r,(2,dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6),dtset%nspinor))
- + + - -
+ + - - +
+ - + - +
- - + -
+ ]
219 : :
220 : : !loop over kpoints
221 [ + + ]: 2 : do ikpt=1,dtset%nkpt
222 : :
223 : : ! number of plane waves for this kpt
224 : 1 : npw = hdr%npwarr(ikpt)
225 : :
226 : : ! allocate arrays dep on number of pw
227 [ - + + - : 3 : ABI_MALLOC(kg_k,(3,npw))
+ - + - -
+ - + ]
228 [ - + + - : 5 : ABI_MALLOC(gpsi,(2,npw,dtset%nspinor,3))
- + + - +
- + - + -
- + - + ]
229 [ - + + - : 4 : ABI_MALLOC(psi,(2,npw,dtset%nspinor))
- + + - +
- + - - +
- + ]
230 [ - + + - : 3 : ABI_MALLOC(kgcart,(3,npw))
+ - + - -
+ - + ]
231 : :
232 : : ! get cartesian coordinates of k+G vectors around this kpoint
233 [ + + ]: 172 : do ig=1,npw
234 [ + + + + : 3249 : kgcart(:,ig) = matmul(gprimd(:,:),dtset%kpt(:,ikpt)+kg(:,ikg+ig))
+ + + + +
- ]
235 [ + + ]: 685 : kg_k (:,ig) = kg(:,ikg+ig)
236 : : end do
237 : :
238 : : ! get gbound
239 : 1 : call sphereboundary(gbound,dtset%istwfk(ikpt),kg_k,dtset%mgfft,npw)
240 : :
241 : : ! loop over bands
242 [ + + ]: 11 : do iband=1,dtset%nband(ikpt)
243 : :
244 : : ! prefactor for sum over bands and kpoints
245 : 10 : prefact_nk = hdr%occ(iocc) * dtset%wtk(ikpt)
246 : :
247 : : ! initialize this wf
248 [ + + + + : 30880 : gpsi=zero
+ + + + ]
249 [ + + + + : 10290 : psi=zero
+ + ]
250 [ + + + + ]: 5140 : psi(:,1:npw,1) = cg(:,icg+1:icg+npw)
251 : :
252 : : ! multiply psi by - i 2 pi G
253 [ + + ]: 1720 : do ig=1,npw
254 [ + + ]: 5130 : gpsi(1,ig,:,1) = two_pi * kgcart(1,ig)*psi(2,ig,:)
255 [ + + ]: 5130 : gpsi(2,ig,:,1) = -two_pi * kgcart(1,ig)*psi(1,ig,:)
256 [ + + ]: 5130 : gpsi(1,ig,:,2) = two_pi * kgcart(2,ig)*psi(2,ig,:)
257 [ + + ]: 5130 : gpsi(2,ig,:,2) = -two_pi * kgcart(2,ig)*psi(1,ig,:)
258 [ + + ]: 5130 : gpsi(1,ig,:,3) = two_pi * kgcart(3,ig)*psi(2,ig,:)
259 [ + + ]: 5140 : gpsi(2,ig,:,3) = -two_pi * kgcart(3,ig)*psi(1,ig,:)
260 : : end do
261 : :
262 : : ! loop over spinorial components
263 [ + + ]: 30 : 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 : 20 : & fft_option,0,one,one,gpu_option=dtset%gpu_option)
270 : :
271 : : ! FT Gpsi_y to real space
272 : : call fourwf(cplex,dummy_denpot,gpsi(:,:,ispinor,2),dummy_fofgout,&
273 : : & dpsidr(:,:,:,:,ispinor,2),gbound,gbound,&
274 : : & hdr%istwfk(ikpt),kg_k,kg_k,dtset%mgfft,mpi_enreg,1,dtset%ngfft,npw,&
275 : : & npw,dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6),&
276 : 20 : & fft_option,0,one,one,gpu_option=dtset%gpu_option)
277 : :
278 : : ! FT Gpsi_z to real space
279 : : call fourwf(cplex,dummy_denpot,gpsi(:,:,ispinor,3),dummy_fofgout,&
280 : : & dpsidr(:,:,:,:,ispinor,3),gbound,gbound,&
281 : : & hdr%istwfk(ikpt),kg_k,kg_k,dtset%mgfft,mpi_enreg,1,dtset%ngfft,npw,&
282 : : & npw,dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6),&
283 : 20 : & fft_option,0,one,one,gpu_option=dtset%gpu_option)
284 : :
285 : : ! FT psi to real space
286 : : call fourwf(cplex,dummy_denpot,psi(:,:,ispinor),dummy_fofgout,&
287 : : & psi_r(:,:,:,:,ispinor),gbound,gbound,&
288 : : & hdr%istwfk(ikpt),kg_k,kg_k,dtset%mgfft,mpi_enreg,1,dtset%ngfft,npw,&
289 : : & npw,dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6),&
290 : 30 : & fft_option,0,one,one,gpu_option=dtset%gpu_option)
291 : :
292 : : end do ! ispinor
293 : :
294 : : ! dpsidr now contains the full derivative of psi wrt space (gradient) in cartesian coordinates
295 : :
296 : : ! get 3 pauli matrix contributions to the current: x,y,z, cart dir, spin dir
297 [ + + ]: 40 : 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 [ + + + + : 108480 : & + 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 [ + + + + : 108480 : & + 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 [ + + + + : 108490 : & + 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 [ + + ]: 30 : do ispinor=1,dtset%nspinor
339 [ + + ]: 330 : do i3=1,dtset%ngfft(3)
340 [ + + ]: 4820 : do i2=1,dtset%ngfft(2)
341 [ + + ]: 72300 : do i1=1,dtset%ngfft(1)
342 : 67500 : irealsp = i1 + (i2-1)*dtset%ngfft(1) + (i3-1)*dtset%ngfft(2)*dtset%ngfft(1)
343 : :
344 [ + + ]: 207000 : 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 : 135000 : & + 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 : 202500 : & - 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 : 10 : icg=icg+npw
382 : 11 : iocc=iocc+1
383 : : end do ! iband
384 : :
385 : 1 : ikg=ikg+npw
386 : :
387 : : ! deallocate arrays dep on npw for this kpoint
388 [ - + ]: 1 : ABI_FREE(kg_k)
389 [ - + ]: 1 : ABI_FREE(gpsi)
390 [ - + ]: 1 : ABI_FREE(psi)
391 [ - + ]: 2 : ABI_FREE(kgcart)
392 : :
393 : : end do ! ikpt
394 : :
395 [ - + ]: 1 : ABI_FREE(dpsidr)
396 [ - + ]: 1 : ABI_FREE(psi_r)
397 [ - + ]: 1 : ABI_FREE(dummy_denpot)
398 [ - + ]: 1 : 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 [ + + + + : 32548 : spincurrent = half * spincurrent
+ + + + +
+ ]
407 : :
408 : : !make array of positions for all points on grid
409 [ - + + - : 5 : ABI_MALLOC(position_op,(3,dtset%ngfft(1),dtset%ngfft(2),dtset%ngfft(3)))
- + + - -
+ + - + -
+ - - + -
+ ]
410 [ + + ]: 16 : do i3=1,dtset%ngfft(3)
411 [ + + ]: 241 : do i2=1,dtset%ngfft(2)
412 [ + + ]: 3615 : do i1=1,dtset%ngfft(1)
413 : : position_op(:,i1,i2,i3) = matmul(hdr%rprimd,(/i1-one,i2-one,i3-one/))&
414 [ + + + + : 74475 : & /(/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 [ - + + - : 6 : ABI_MALLOC(vso_realspace,(2,dtset%ngfft(1)*dtset%ngfft(2)*dtset%ngfft(3),dtset%nspinor,dtset%nspinor,3))
- + + - -
+ + - + -
+ - + - -
+ - + ]
443 : :
444 : 1 : 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 [ - + ]: 1 : ABI_FREE(density)
453 : :
454 : 1 : realrecip = 0 ! real space for xsf output
455 [ - + + - : 3 : ABI_MALLOC(xcart,(3,dtset%natom))
+ - + - -
+ - + ]
456 : 1 : 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 [ + + ]: 4 : 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 [ + + + + : 32547 : rescale_current = maxval(abs(spincurrent(:, :, :, :, ispindir)))
+ + + + +
+ + - -
+ ]
468 [ + + ]: 3 : 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 [ + + + + : 10 : & sum(hdr%rprimd(:,2)**2), sum(hdr%rprimd(:,3)**2)))
+ + ]
473 : : end if
474 : :
475 [ + - ]: 3 : filnam=trim(dtfil%fnameabo_spcur)//spin_symbol(ispindir)//".xsf"
476 [ - + ]: 3 : if (open_file(filnam,message,newunit=spcur_unit,status='unknown') /= 0) then
477 : 0 : ABI_ERROR(message)
478 : : end if
479 : :
480 : : ! print header
481 : 3 : write (spcur_unit,'(a)') '#'
482 : 3 : write (spcur_unit,'(a)') '# Xcrysden format file'
483 : 3 : write (spcur_unit,'(a)') '# spin current density, for all real space points'
484 : 3 : write (spcur_unit,'(a,3(I5,1x))') '# fft grid is ', dtset%ngfft(1), dtset%ngfft(2), dtset%ngfft(3)
485 : 3 : write (spcur_unit,'(a,a,a)') '# ', spin_symbol(ispindir), '-spin current, full vector '
486 : :
487 : 3 : write (spcur_unit,'(a)') 'ATOMS'
488 [ + + ]: 6 : do iatom = 1, dtset%natom
489 : 6 : write (spcur_unit,'(I4, 2x, 3(E16.6, 1x))') int(dtset%znucl(dtset%typat(iatom))), xcart(:,iatom)
490 : : end do
491 : :
492 [ + + ]: 48 : do i3=1,dtset%ngfft(3)
493 [ + + ]: 723 : do i2=1,dtset%ngfft(2)
494 [ + + ]: 10845 : do i1=1,dtset%ngfft(1)
495 : 10125 : write (spcur_unit,'(a, 3(E10.3),2x, 3(E20.10))') 'X ', &
496 [ + + ]: 51300 : & position_op(:, i1, i2, i3), spincurrent(i1, i2, i3, :, ispindir)*rescale_current
497 : : end do
498 : : end do
499 : : end do
500 : 4 : 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 [ - + + - : 3 : ABI_MALLOC(datagrid, (dtset%ngfft(1)*dtset%ngfft(2)*dtset%ngfft(3)))
+ - + - -
+ - + ]
510 : :
511 [ + + ]: 4 : do ispindir=1,3
512 [ + + ]: 10 : do icplex=1,2
513 [ + + ]: 21 : do ispinor=1,dtset%nspinor
514 [ + + ]: 42 : do ispinorp=1,dtset%nspinor
515 : :
516 : : ! for the moment only print out if non zero
517 [ + + + + ]: 81024 : if (abs(sum(vso_realspace(icplex, :, ispinor, ispinorp, ispindir))) < tol8) cycle
518 : :
519 : : filnam=trim(dtfil%fnameabo_vso)//"_spin_"//spin_symbol(ispindir)//"_"//&
520 [ + - ]: 12 : & spinor_sym(ispinor)//spinor_sym(ispinorp)//"_"//realimag(icplex)//".xsf"
521 : :
522 [ - + ]: 12 : if (open_file(filnam,message,newunit=spcur_unit,status='unknown') /= 0) then
523 : 0 : ABI_ERROR(message)
524 : : end if
525 : :
526 : : ! print header
527 : 12 : write (spcur_unit,'(a)') '#'
528 : 12 : write (spcur_unit,'(a)') '# Xcrysden format file'
529 : 12 : write (spcur_unit,'(a)') '# spin-orbit potential (space diagonal), for all real space points'
530 : 12 : write (spcur_unit,'(a)') '# Real part first, then imaginary part'
531 : 12 : write (spcur_unit,'(a,3(I5,1x))') '# fft grid is ', dtset%ngfft(1), dtset%ngfft(2), dtset%ngfft(3)
532 : 12 : write (spcur_unit,'(a,a,a)') '# ', spin_symbol(ispindir), '-spin contribution '
533 : 12 : write (spcur_unit,'(a,a,a)') '# ', spinor_sym(ispinor)//spinor_sym(ispinorp), &
534 : 24 : '-spin element of the spinor 2x2 matrix '
535 : 12 : write (spcur_unit,'(a,a)') '# cart x * cart y * cart z ***',&
536 : 24 : & ' up-up component * up-down * down-up * down-down '
537 : :
538 : : ! Build contiguous array
539 [ - + + + ]: 40524 : 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 : 12 : & 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 : 36 : close (spcur_unit)
564 : :
565 : : end do ! ispinorp
566 : : end do ! ispinor
567 : : end do ! icplex
568 : : end do ! end ispindir
569 : :
570 [ - + ]: 1 : ABI_FREE(datagrid)
571 [ - + ]: 1 : ABI_FREE(vso_realspace)
572 : : !deallocate (vso_realspace_nl)
573 [ - + ]: 1 : ABI_FREE(position_op)
574 [ - + ]: 1 : ABI_FREE(spincurrent)
575 [ - + ]: 1 : ABI_FREE(xcart)
576 : :
577 : 1 : write(std_out,*) ' Exiting subroutine spin_current '
578 : :
579 : 1 : end subroutine spin_current
580 : : !!***
581 : :
582 : : !!****f* m_spin_current/vso_realspace_local
583 : : !! NAME
584 : : !! vso_realspace_local
585 : : !!
586 : : !! FUNCTION
587 : : !! Calculate real space (local - (r,r)) values of the SO part of the
588 : : !! pseudopotential. Reconstructed explicitly in the HGH/GTH case.
589 : : !!
590 : : !! INPUTS
591 : : !!
592 : : !! OUTPUT
593 : : !!
594 : : !! SIDE EFFECTS
595 : : !!
596 : : !! NOTES
597 : : !!
598 : : !! SOURCE
599 : :
600 : 1 : 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 : 2 : real(dp) :: xcart(3,dtset%natom),splint_x(1),splint_y(1)
619 : 1 : real(dp), allocatable :: radial_grid(:)
620 : 1 : real(dp), allocatable :: prefact_ijl(:,:,:,:),tmpvso(:),tmpvso_pp(:)
621 : 1 : real(dp), allocatable :: vso_radial(:,:),vso_radial_pp(:,:),tmp_spline(:)
622 : 1 : real(dp), allocatable :: offdiag_l_fact(:,:,:),kpar_matrix(:,:)
623 : :
624 : : ! *********************************************************************
625 : :
626 : : !recalculate xcart (option = 1)
627 : 1 : call xred2xcart(dtset%natom,hdr%rprimd,xcart,hdr%xred)
628 : :
629 : 1 : 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 [ + + + + : 7 : 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 : 1 : nradgrid = 201 ! this is heuristic
663 [ - + ]: 1 : ABI_MALLOC(radial_grid,(nradgrid))
664 [ + + ]: 202 : do iradgrid=1,nradgrid
665 : 202 : radial_grid(iradgrid) = (iradgrid-1)*radial_cutoff/(nradgrid-1)
666 : : end do
667 : :
668 : : !calculate prefactors independent of r
669 [ - + + - : 4 : ABI_MALLOC(prefact_ijl,(3,3,0:lmax,psps%npsp))
- + + - +
- + - - +
- + ]
670 [ - + + - : 3 : ABI_MALLOC(offdiag_l_fact,(3,3,0:lmax))
+ - + - -
+ - + ]
671 [ - + ]: 1 : 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 [ + + + + : 40 : offdiag_l_fact = zero
+ + ]
676 : : !l=0
677 : 1 : offdiag_l_fact(1,2,0) = -half*sqrt(three/five)
678 : 1 : offdiag_l_fact(1,3,0) = half*sqrt(five/21._dp)
679 : 1 : offdiag_l_fact(2,3,0) = -half*sqrt(100._dp/63._dp)
680 : : !l=1
681 : 1 : offdiag_l_fact(1,2,1) = -half*sqrt(five/seven)
682 : 1 : offdiag_l_fact(1,3,1) = sixth*sqrt(35._dp/11._dp)
683 : 1 : offdiag_l_fact(2,3,1) = -sixth*14._dp/sqrt(11._dp)
684 : : !l=2
685 [ + - ]: 1 : if (lmax >= 2) then
686 : 1 : offdiag_l_fact(1,2,2) = -half*sqrt(seven/nine)
687 : 1 : offdiag_l_fact(1,3,2) = half*sqrt(63._dp/143._dp)
688 : 1 : offdiag_l_fact(2,3,2) = -half*18._dp /sqrt(143._dp)
689 : : end if
690 : : !l=3
691 [ - + ]: 1 : if (lmax >= 3) then
692 : 0 : offdiag_l_fact(1,2,3) = zero
693 : 0 : offdiag_l_fact(1,3,3) = zero
694 : 0 : 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 [ + + + + : 41 : prefact_ijl = zero
+ + + + ]
698 [ + + ]: 4 : do l=0,lmax
699 : : ! first the diagonal i=j term
700 [ + + ]: 12 : do i=1,3
701 : 9 : call gamma_function(l+(4._dp*i-1._dp)*0.5_dp, gammai)
702 : 9 : gammai = sqrt(gammai)
703 : 9 : rcexponent = 2*l+2*i+2*i-1
704 [ + + ]: 21 : 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 : 18 : & * gammai * gammai)
708 : : end do
709 : : end do
710 : : ! now the off diagonal elements
711 [ + + ]: 6 : do ipsp=1,psps%npsp
712 : 3 : kpar_matrix(1,2) = offdiag_l_fact (1,2,l)* psps%gth_params%psp_k_par(l+1,2,ipsp)
713 : 3 : kpar_matrix(2,1) = kpar_matrix(1,2)
714 : 3 : kpar_matrix(1,3) = offdiag_l_fact (1,3,l)* psps%gth_params%psp_k_par(l+1,3,ipsp)
715 : 3 : kpar_matrix(3,1) = kpar_matrix(1,3)
716 : 3 : kpar_matrix(2,3) = offdiag_l_fact (2,3,l)* psps%gth_params%psp_k_par(l+1,3,ipsp)
717 : 6 : 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 [ + - ]: 3 : if (l > 2) cycle
723 : :
724 [ + + ]: 13 : do i=1,3
725 : 9 : call gamma_function(l+(4._dp*i-1._dp)*0.5_dp, gammai)
726 : 9 : gammai = sqrt(gammai)
727 [ + + ]: 39 : do j=1,3
728 [ + + ]: 27 : if (j==i) cycle
729 : 18 : rcexponent = 2*l+2*i+2*j-1
730 : 18 : call gamma_function(l+(4._dp*j-1._dp)*0.5_dp,gammaj)
731 : 18 : gammaj = sqrt(gammaj)
732 [ + + ]: 45 : 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 : 45 : & * gammai * gammaj )
736 : : end do
737 : : end do
738 : : end do
739 : : end do
740 : :
741 [ - + ]: 1 : ABI_FREE(kpar_matrix)
742 : 1 : ABI_FREE(offdiag_l_fact)
743 : :
744 [ + + + + : 41 : prefact_ijl = prefact_ijl * two
+ + + + ]
745 : :
746 : : !calculate v_SO on radial grid
747 : : ! MGNAG Runtime Error: *** Arithmetic exception: Floating invalid operation - aborting
748 [ - + + - : 3 : ABI_MALLOC(vso_radial,(nradgrid,psps%npsp))
+ - + - -
+ - + ]
749 [ + + + + ]: 203 : vso_radial = zero
750 [ + + ]: 4 : do l=0,lmax
751 : 606 : lfact=(2._dp*l+1._dp)/four/pi
752 [ + + ]: 607 : do iradgrid=1,nradgrid
753 : 603 : norm_rel_pos = radial_grid(iradgrid)
754 [ + + ]: 1209 : do ipsp=1,psps%npsp
755 : 603 : expfact = exp(-norm_rel_pos**2 / (psps%gth_params%psppar(l+1,0,ipsp))**2)
756 : :
757 [ + + ]: 3015 : do i=1,3
758 [ + + ]: 7839 : do j=1,3
759 : 5427 : rcexponent = 2*l +2*i+2*j-4
760 [ + + ]: 7236 : if(prefact_ijl(i,j,l,ipsp)/=0) then !vz_d 0**0
761 : : vso_radial(iradgrid,ipsp) = vso_radial(iradgrid,ipsp) + &
762 : 1005 : & 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 [ - + + - : 2 : ABI_MALLOC(vso_radial_pp,(nradgrid,psps%npsp))
- + - + ]
772 : :
773 [ - + ]: 1 : ABI_MALLOC(tmp_spline,(nradgrid))
774 [ - + ]: 1 : ABI_MALLOC(tmpvso,(nradgrid))
775 [ - + ]: 1 : ABI_MALLOC(tmpvso_pp,(nradgrid))
776 [ + + ]: 2 : do ipsp=1,psps%npsp
777 [ - + + + ]: 203 : tmpvso = vso_radial(:,ipsp)
778 : 1 : call spline( radial_grid, tmpvso, nradgrid, zero, radial_grid(nradgrid), tmpvso_pp )
779 [ + + ]: 203 : vso_radial_pp(:,ipsp) = tmpvso_pp
780 : : end do
781 [ - + ]: 1 : ABI_FREE(tmp_spline)
782 [ - + ]: 1 : ABI_FREE(tmpvso)
783 [ - + ]: 1 : 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 [ + + + + : 121522 : vso_realspace = zero
+ + + + +
+ ]
790 [ + + ]: 2 : do iatom=1,dtset%natom
791 : : ! atom type will be dtset%typat(iatom)
792 : :
793 : : ! for each point on grid
794 [ + + ]: 17 : do ir3=1,dtset%ngfft(3)
795 [ + + ]: 241 : do ir2=1,dtset%ngfft(2)
796 [ + + ]: 3615 : do ir1=1,dtset%ngfft(1)
797 : 3375 : irealsp = ir1 + (ir2-1)*dtset%ngfft(1) + (ir3-1)*dtset%ngfft(2)*dtset%ngfft(1)
798 : :
799 : : ! relative position from atom to point
800 [ + + ]: 13500 : relative_position = position_op(:,ir1,ir2,ir3) - xcart(:,iatom)
801 : 3375 : x=relative_position(1)
802 : 3375 : y=relative_position(2)
803 : 3375 : z=relative_position(3)
804 : :
805 : : ! calculate norm^2
806 : 3375 : 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 [ + + ]: 3375 : if (norm_rel_pos > radial_cutoff*radial_cutoff) cycle
810 : :
811 : : ! calculate norm
812 : 459 : 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 : 459 : & vso_radial_pp(:,dtset%typat(iatom)),1,splint_x,splint_y)
817 : 459 : 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 [ + + + + : 5049 : & vso_interpol * reshape((/y, zero,zero,-y/),(/2,2/))
+ + ]
826 : : vso_realspace(2,irealsp,:,:,1) = vso_realspace(2,irealsp,:,:,1) + &
827 [ + + + + : 5049 : & vso_interpol * reshape((/zero,z, -z, zero/),(/2,2/))
+ + ]
828 : :
829 : : vso_realspace(1,irealsp,:,:,2) = vso_realspace(1,irealsp,:,:,2) + &
830 [ + + + + : 5049 : & vso_interpol * reshape((/-x, z, z, x/),(/2,2/))
+ + ]
831 : : vso_realspace(2,irealsp,:,:,2) = vso_realspace(2,irealsp,:,:,2) + &
832 [ + + + + ]: 3213 : & vso_interpol * reshape((/zero,zero,zero,zero/),(/2,2/))
833 : :
834 : : vso_realspace(1,irealsp,:,:,3) = vso_realspace(1,irealsp,:,:,3) + &
835 [ + + + + : 5049 : & vso_interpol * reshape((/zero,-y, -y, zero/),(/2,2/))
+ + ]
836 : : vso_realspace(2,irealsp,:,:,3) = vso_realspace(2,irealsp,:,:,3) + &
837 [ + + + + : 8190 : & 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 [ - + ]: 1 : ABI_FREE(prefact_ijl)
845 [ - + ]: 1 : ABI_FREE(vso_radial)
846 [ - + ]: 1 : ABI_FREE(vso_radial_pp)
847 [ - + ]: 1 : ABI_FREE(radial_grid)
848 : :
849 : 1 : end subroutine vso_realspace_local
850 : : !!***
851 : :
852 : 171 : end module m_spin_current
853 : : !!***
|
