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 : !!***
|
