Line data Source code
1 : !!****m* ABINIT/m_dfpt_rhotov
2 : !! NAME
3 : !! m_dfpt_rhotov
4 : !!
5 : !! FUNCTION
6 : !!
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1999-2026 ABINIT group (XG, DRH, MT, SPr)
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_dfpt_rhotov
23 :
24 : use defs_basis
25 : use m_abicore
26 : use m_errors
27 : use m_cgtools
28 :
29 : use defs_abitypes, only : MPI_type
30 : use m_time, only : timab
31 : use m_spacepar, only : hartrestr, hartre
32 : use m_dfpt_mkvxc, only : dfpt_mkvxc, dfpt_mkvxc_noncoll
33 : use m_dfpt_mkvxcstr, only : dfpt_mkvxcstr
34 : use m_dens, only : calcdenmagsph
35 : use m_numeric_tools, only : wrap2_zero_one
36 :
37 : implicit none
38 :
39 : private
40 : !!***
41 :
42 : public :: dfpt_rhotov
43 : !!***
44 :
45 : contains
46 : !!***
47 :
48 : !!****f* ABINIT/dfpt_rhotov
49 : !! NAME
50 : !! dfpt_rhotov
51 : !!
52 : !! FUNCTION
53 : !! This routine is called to compute, from a given 1st-order total density
54 : !! - the trial (local) 1st-order potential and/or the residual potential,
55 : !! - some contributions to the 2nd-order energy
56 : !!
57 : !! INPUTS
58 : !! cplex: if 1, real space 1-order WF on FFT grid are REAL; if 2, COMPLEX
59 : !! gsqcut=cutoff on (k+G)^2 (bohr^-2)
60 : !! icutcoul= type of Coulomb cutoff to apply
61 : !! idir=direction of atomic displacement (=1,2 or 3 : displacement of atom ipert along the 1st, 2nd or 3rd axis).
62 : !! ipert=type of the perturbation
63 : !! ixc= choice of exchange-correlation scheme
64 : !! kxc(nfft,nkxc)=exchange-correlation kernel
65 : !! mpi_enreg=information about MPI parallelization
66 : !! magpen=energy shift to apply on the spin degrees of freedom
67 : !! mpatpol(2)=initial and final atomic positions whose local magnetic moments will be penalized
68 : !! mpdir(3)=directions of the magnetic moments to be penalized
69 : !! natom=number of atoms in cell.
70 : !! nfft=(effective) number of FFT grid points (for this processor)
71 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
72 : !! nhat(nfft,nspden*nhatdim)= -PAW only- compensation density
73 : !! nhat1(cplex*nfft,2nspden*usepaw)= -PAW only- 1st-order compensation density
74 : !! nhat1gr(cplex*nfft,nspden,3*nhat1grdim)= -PAW only- gradients of 1st-order compensation density
75 : !! nhat1grdim= -PAW only- 1 if nhat1gr array is used ; 0 otherwise
76 : !! nkxc=second dimension of the array kxc, see rhotoxc.f for a description
77 : !! non_magnetic_xc= if true, handle density/potential as non-magnetic (even if it is)
78 : !! nspden=number of spin-density components
79 : !! ntypat=number of atom types
80 : !! n3xccc=dimension of xccc3d1 ; 0 if no XC core correction is used
81 : !! optene=0: the contributions to the 2nd order energy are not computed
82 : !! 1: the contributions to the 2nd order energy are computed
83 : !! optres=0: the trial potential residual is computed ; the input potential value is kept
84 : !! 1: the new value of the trial potential is computed in place of the input value
85 : !! qphon(3)=reduced coordinates for the phonon wavelength
86 : !! ratopt=1: spheres around atoms are build in real space
87 : !! 2: spheres around atoms are build in reciprocal space
88 : !! ratsm=smearing width for ratsph
89 : !! ratsph(ntypat)=radius of spheres around atoms
90 : !! rhog(2,nfft)=array for Fourier transform of GS electron density
91 : !! rhog1(2,nfft)=RF electron density in reciprocal space
92 : !! rhor(nfft,nspden)=array for GS electron density in electrons/bohr**3.
93 : !! rhor1(cplex*nfft,nspden)=RF electron density in real space (electrons/bohr**3).
94 : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
95 : !! typat(natom)=type of each atom
96 : !! ucvol=unit cell volume in ($\textrm{bohr}^{3}$)
97 : !! usepaw= 0 for non paw calculation; =1 for paw calculation
98 : !! usexcnhat= -PAW only- flag controling use of compensation density in Vxc
99 : !! vcutgeo(3)= array to describe the geometry of the Coulomb cutoff
100 : !! vpsp1(cplex*nfft)=first-order derivative of the ionic potential
101 : !! xccc3d1(cplex*n3xccc)=3D change in core charge density, see n3xccc
102 : !! xred(3,natom)=reduced dimensionless atomic coordinates
103 : !!
104 : !! OUTPUT
105 : !! vhartr1(cplex*nfft)=1-order Hartree potential (not output if size=0)
106 : !! vxc1(cplex*nfft,nspden)= 1st-order XC potential (not output if size=0)
107 : !! ==== if optene==1
108 : !! ehart01=inhomogeneous 1st-order Hartree part of 2nd-order total energy
109 : !! ehart1=1st-order Hartree part of 2nd-order total energy
110 : !! exc1=1st-order exchange-correlation part of 2nd-order total energy
111 : !! elpsp1=1st-order local pseudopot. part of 2nd-order total energy.
112 : !! emagpen1
113 : !! ==== if optres==0
114 : !! vresid1(cplex*nfft,nspden)=potential residual
115 : !! vres2=square of the norm of the residual
116 : !!
117 : !! SIDE EFFECTS
118 : !! ==== if optres==1
119 : !! vtrial1(cplex*nfft,nspden)= new value of 1st-order trial potential
120 : !!
121 : !! SOURCE
122 :
123 21014 : subroutine dfpt_rhotov(cplex,ehart01,ehart1,elmag1,elpsp1,emagpen1,exc1,gsqcut,icutcoul,idir,ipert,&
124 21014 : & ixc,kxc,magpen,mpatpol,mpdir,mpi_enreg,natom,nfft,ngfft,nhat,nhat1,nhat1gr,nhat1grdim,nkxc,nspden,ntypat,n3xccc,&
125 21014 : & non_magnetic_xc,optene,optres,qphon,ratsm,ratsph,rhog,rhog1,rhor,rhor1,rprimd,typat,ucvol,&
126 21014 : & usepaw,usexcnhat,vcutgeo,vhartr1,vpsp1,vresid1,vres2,vtrial1,vxc,vxc1,xccc3d1,ixcrot,xred,qgbt,use_gbt)
127 :
128 : !Arguments ------------------------------------
129 : !scalars
130 : integer,intent(in) :: cplex,icutcoul,idir,ipert,ixc,n3xccc,natom,nfft,nhat1grdim,nkxc,nspden
131 : integer,intent(in) :: ntypat,optene,optres,usepaw,usexcnhat,ixcrot
132 : logical,intent(in) :: non_magnetic_xc
133 : integer,intent(in) :: use_gbt
134 : real(dp),intent(in):: qgbt(3)
135 : real(dp),intent(in) :: gsqcut,magpen,ratsm,ucvol
136 : real(dp),intent(inout) :: ehart01,elpsp1,ehart1,exc1,elmag1,emagpen1
137 : real(dp),intent(out) :: vres2
138 : type(MPI_type),intent(in) :: mpi_enreg
139 : !arrays
140 : integer,intent(in) :: ngfft(18),typat(natom)
141 : integer,intent(in) :: mpatpol(2),mpdir(3)
142 : real(dp),intent(in) :: kxc(nfft,nkxc)
143 : real(dp),intent(in) :: vxc(nfft,nspden)
144 : real(dp),intent(in) :: nhat(nfft,nspden)
145 : real(dp),intent(in) :: nhat1(cplex*nfft,nspden) !vz_d
146 : real(dp),intent(in) :: nhat1gr(cplex*nfft,nspden,3*nhat1grdim)
147 : real(dp),intent(in) :: qphon(3),ratsph(ntypat),rhog(2,nfft)
148 : real(dp),intent(in) :: rhog1(2,nfft)
149 : real(dp),target,intent(in) :: rhor(nfft,nspden),rhor1(cplex*nfft,nspden)
150 : real(dp),intent(in) :: rprimd(3,3),vpsp1(cplex*nfft)
151 : real(dp),intent(in) :: xccc3d1(cplex*n3xccc)
152 : real(dp),intent(inout) :: vtrial1(cplex*nfft,nspden)
153 : real(dp),intent(out) :: vresid1(cplex*nfft,nspden)
154 : real(dp),target,intent(out) :: vhartr1(:),vxc1(:,:)
155 : real(dp),intent(in) :: vcutgeo(3)
156 : real(dp),intent(in) :: xred(3,natom)
157 :
158 : !Local variables-------------------------------
159 : !scalars
160 : integer :: idir_eff,ifft,ispden,nfftot,option
161 : integer :: optnc,nkxc_cur,prtopt
162 : logical :: vhartr1_allocated,vxc1_allocated
163 : real(dp) :: doti,elpsp10
164 : !arrays
165 : real(dp) :: tsec(20)
166 21014 : real(dp),allocatable :: rhor1_nohat(:,:),vhartr01(:),vxc1val(:,:)
167 21014 : real(dp),pointer :: rhor1_(:,:),vhartr1_(:),vxc1_(:,:),v1hspinfield(:,:)
168 21014 : real(dp),allocatable :: fatsph(:,:),intgden(:,:,:),rhomag(:,:),vmagpen1(:,:)
169 21014 : real(dp),allocatable :: taumr(:,:,:)
170 :
171 : ! *********************************************************************
172 :
173 21014 : call timab(157,1,tsec)
174 :
175 : !FR EB SPr
176 21014 : if (nspden==4) then
177 242 : if(usepaw==1) then
178 0 : ABI_ERROR('DFPT with nspden=4 works only for norm-conserving psp!')
179 : end if
180 : end if
181 :
182 : !Get size of FFT grid
183 21014 : nfftot=ngfft(1)*ngfft(2)*ngfft(3)
184 :
185 : !Eventually allocate temporary memory space
186 21014 : vhartr1_allocated=(size(vhartr1)>0)
187 21014 : if (vhartr1_allocated) then
188 21014 : vhartr1_ => vhartr1
189 : else
190 0 : ABI_MALLOC(vhartr1_,(cplex*nfft))
191 : end if
192 63042 : vxc1_allocated=(size(vxc1)>0)
193 21014 : if (vxc1_allocated) then
194 18718 : vxc1_ => vxc1
195 : else
196 9184 : ABI_MALLOC(vxc1_,(cplex*nfft,nspden))
197 : end if
198 :
199 : !If needed, store pseudo density without charge compensation
200 21014 : if (usepaw==1.and.usexcnhat==0) then
201 7172 : ABI_MALLOC(rhor1_,(cplex*nfft,nspden))
202 134952082 : rhor1_(:,:)=rhor1(:,:)-nhat1(:,:)
203 : else
204 19221 : rhor1_ => rhor1
205 : end if
206 :
207 : !Uniform Zeeman or scalar potential
208 21014 : if(ipert==natom+5.or.ipert==natom+6)then
209 82 : if (ipert==natom+5) idir_eff= idir
210 82 : if (ipert==natom+6) idir_eff= 4
211 328 : ABI_MALLOC(v1hspinfield,(cplex*nfft,nspden))
212 82 : call dfpt_v1hspinfield(nspden,nfft,cplex,idir_eff,v1hspinfield)
213 : end if
214 :
215 : !Preconditioned DFPT
216 21014 : if((ipert>natom+11.and.ipert<=2*natom+11).or.abs(magpen)>tol6) then
217 :
218 : !Compute the first-order magnetic moments.
219 62 : prtopt=1
220 310 : ABI_MALLOC(intgden,(cplex,nspden,natom))
221 186 : ABI_MALLOC(rhomag,(2,nspden))
222 248 : ABI_MALLOC(fatsph,(nfft,natom))
223 310 : ABI_MALLOC(taumr,(nfft,natom,3))
224 : call calcdenmagsph(mpi_enreg,natom,nfft,ngfft,nspden,&
225 : & ntypat,ratsm,ratsph,rhor1,rprimd,typat,xred,prtopt,cplex,&
226 62 : & qgbt,use_gbt,intgden=intgden,rhomag=rhomag,fatsph=fatsph,qphon=qphon,taumr=taumr)
227 : end if
228 :
229 21014 : if(ipert>natom+11.and.ipert<=2*natom+11)then
230 68 : ABI_MALLOC(v1hspinfield,(cplex*nfft,nspden))
231 : call dfpt_v1hspinfield_atsph(cplex,fatsph,idir,ipert,natom,nfft,nspden,&
232 17 : & qphon,taumr,v1hspinfield)
233 : end if
234 :
235 84056 : ABI_MALLOC(vmagpen1,(cplex*nfft,nspden))
236 296720089 : vmagpen1=zero
237 21014 : if (abs(magpen) > tol6) then
238 : call dfpt_v1magpen(cplex,emagpen1,fatsph,intgden,magpen,mpatpol,&
239 62 : & mpdir,natom,nfft,nspden,qphon,rhomag,taumr,vmagpen1)
240 : end if
241 :
242 : !------ Compute 1st-order Hartree potential (and energy) ----------------------
243 21014 : call hartre(cplex,gsqcut,icutcoul,0,mpi_enreg,nfft,ngfft,1,zero,rhog1,rprimd,vcutgeo,vhartr1_,qpt=qphon)
244 :
245 21014 : if (optene>0) then
246 18505 : call dotprod_vn(cplex,rhor1,ehart1,doti,nfft,nfftot,1,1,vhartr1_,ucvol)
247 : end if
248 :
249 21014 : if (optene>0) ehart01=zero
250 21014 : if(ipert==natom+3 .or. ipert==natom+4) then
251 5025 : ABI_MALLOC(vhartr01,(cplex*nfft))
252 1675 : call hartrestr(gsqcut,idir,ipert,mpi_enreg,natom,nfft,ngfft,rhog,rprimd,vhartr01)
253 1675 : if (optene>0) then
254 1507 : call dotprod_vn(cplex,rhor1,ehart01,doti,nfft,nfftot,1,1,vhartr01,ucvol)
255 1507 : ehart01=two*ehart01
256 1507 : ehart1=ehart1+ehart01
257 : end if
258 : ! Note that there is a factor 2.0_dp difference with the similar GS formula
259 11099319 : vhartr1_(:)=vhartr1_(:)+vhartr01(:)
260 :
261 1675 : ABI_FREE(vhartr01)
262 : end if
263 :
264 : !------ Compute 1st-order XC potential (and energy) ----------------------
265 : !(including the XC core correction)
266 :
267 : !Compute Vxc^(1) (with or without valence contribution according to options)
268 21014 : option=0;if (optene==0) option=1
269 21014 : if(ipert==natom+3.or.ipert==natom+4) then
270 : call dfpt_mkvxcstr(cplex,idir,ipert,kxc,mpi_enreg,natom,nfft,ngfft,nhat,&
271 : & nhat1,nkxc,non_magnetic_xc,nspden,n3xccc,option,qphon,rhor,rhor1,rprimd,&
272 1675 : & usepaw,usexcnhat,vxc1_,xccc3d1)
273 : else
274 : ! FR EB non-collinear magnetism
275 : ! the second nkxc should be nkxc_cur (see 67_common/nres2vres.F90)
276 19339 : if (nspden==4) then
277 242 : optnc=1
278 242 : nkxc_cur=nkxc ! TODO: remove nkxc_cur?
279 8435350 : vxc1_=zero
280 :
281 : call dfpt_mkvxc_noncoll(cplex,ixc,kxc,mpi_enreg,nfft,ngfft,nhat,usepaw,nhat1,usepaw,nhat1gr,nhat1grdim,nkxc,&
282 242 : & non_magnetic_xc,nspden,n3xccc,optnc,option,qphon,rhor,rhor1,rprimd,usexcnhat,vxc,vxc1_,xccc3d1,ixcrot=ixcrot)
283 :
284 : else
285 : call dfpt_mkvxc(cplex,ixc,kxc,mpi_enreg,nfft,ngfft,nhat1,usepaw,nhat1gr,nhat1grdim,nkxc,&
286 19097 : & non_magnetic_xc,nspden,n3xccc,option,qphon,rhor1,rprimd,usexcnhat,vxc1_,xccc3d1)
287 : end if !nspden==4
288 : end if
289 :
290 : !Compute local contribution to 2nd-order energy (includes Vxc and Vpsp and Vmag)
291 21014 : if (optene>0) then
292 18505 : if (usepaw==0) then
293 15059 : call dotprod_vn(cplex,rhor1,elpsp10,doti,nfft,nfftot,nspden,1,vxc1_,ucvol)
294 15059 : call dotprod_vn(cplex,rhor1,elpsp1 ,doti,nfft,nfftot,1 ,1,vpsp1,ucvol)
295 15059 : if (ipert==natom+5.or.ipert==natom+6.or.(ipert>natom+11.and.ipert<=2*natom+11)) then
296 87 : call dotprod_vn(cplex,rhor1,elmag1 ,doti,nfft,nfftot,nspden,1,v1hspinfield,ucvol)
297 : !A factor of four, present in <0|H^1|1> terms, compensates the missing
298 : !half factor
299 87 : elmag1=two*elmag1
300 : end if
301 : else
302 3446 : if (usexcnhat/=0) then
303 6237 : ABI_MALLOC(rhor1_nohat,(cplex*nfft,1))
304 28858687 : rhor1_nohat(:,1)=rhor1(:,1)-nhat1(:,1)
305 2079 : call dotprod_vn(cplex,rhor1 ,elpsp10,doti,nfft,nfftot,nspden,1,vxc1_,ucvol)
306 2079 : call dotprod_vn(cplex,rhor1_nohat,elpsp1 ,doti,nfft,nfftot,1 ,1,vpsp1,ucvol)
307 2079 : ABI_FREE(rhor1_nohat)
308 : else
309 1367 : call dotprod_vn(cplex,rhor1_,elpsp10,doti,nfft,nfftot,nspden,1,vxc1_,ucvol)
310 1367 : call dotprod_vn(cplex,rhor1_,elpsp1 ,doti,nfft,nfftot,1 ,1,vpsp1,ucvol)
311 : end if
312 : end if
313 :
314 : ! Note that there is a factor 2 difference with the similar GS formula
315 18505 : elpsp1=two*(elpsp1+elpsp10)
316 : end if
317 :
318 :
319 : !Compute XC valence contribution exc1 and complete eventually Vxc^(1)
320 21014 : if (optene>0) then
321 55515 : ABI_MALLOC(vxc1val,(cplex*nfft,nspden))
322 264248609 : vxc1val=zero
323 18505 : option=2
324 : !FR SPr EB non-collinear magnetism
325 18505 : if (nspden==4) then
326 215 : optnc=1
327 215 : nkxc_cur=nkxc
328 : call dfpt_mkvxc_noncoll(cplex,ixc,kxc,mpi_enreg,nfft,ngfft,nhat,usepaw,nhat1,usepaw,nhat1gr,nhat1grdim,nkxc,&
329 215 : & non_magnetic_xc,nspden,n3xccc,optnc,option,qphon,rhor,rhor1,rprimd,usexcnhat,vxc,vxc1val,xccc3d1,ixcrot=ixcrot)
330 : else
331 : call dfpt_mkvxc(cplex,ixc,kxc,mpi_enreg,nfft,ngfft,nhat1,usepaw,nhat1gr,nhat1grdim,nkxc,&
332 18290 : & non_magnetic_xc,nspden,n3xccc,option,qphon,rhor1,rprimd,usexcnhat,vxc1val,xccc3d1)
333 : end if !nspden==4
334 :
335 264248609 : vxc1_(:,:)=vxc1_(:,:)+vxc1val(:,:)
336 18505 : call dotprod_vn(cplex,rhor1_,exc1,doti,nfft,nfftot,nspden,1,vxc1val,ucvol)
337 18505 : ABI_FREE(vxc1val)
338 : end if
339 :
340 21014 : if (usepaw==1.and.usexcnhat==0) then
341 1793 : ABI_FREE(rhor1_)
342 : end if
343 :
344 : !DEBUG (do not take away)
345 : !Compute NSC energy ensc1 associated with rhor1 in vtrial1, for debugging purposes
346 : !call dotprod_vn(cplex,rhor1,ensc1,doti,nfft,nfftot,nspden,1,vtrial1,ucvol)
347 : !write(std_out,*)' ek0+eeig0+eloc0=',ek0+eeig0+eloc0
348 : !write(std_out,*)' ensc1=',ensc1
349 : !Compute NSC energy associated with vtrial1, for debugging purposes
350 : !call dotprod_vn(cplex,rhor1,ensc1,doti,mpi_enreg,nfft,nfftot,nspden,1,vtrial1,ucvol)
351 : !ensc1=ensc1+half*enl1
352 : !write(std_out,*)' dfpt_rhotov : check NSC energy, diff=',&
353 : !& ek0+edocc+eeig0+eloc0+enl0+ensc1
354 : !write(std_out,*)' evarNSC=',ek0+edocc+eeig0+eloc0+enl0
355 : !write(std_out,*)' ensc1,exc1=',ensc1,exc1
356 : !ENDDEBUG
357 :
358 : !Here, vhartr1 contains Hartree potential, vpsp1 contains local psp,
359 : !while vxc1 contain xc potential
360 :
361 : !------ Produce residual vector and square of norm of it -------------
362 : !(only if requested ; if optres==0)
363 21014 : if (optres==0) then
364 : !$OMP PARALLEL DO COLLAPSE(2)
365 38514 : do ispden=1,min(nspden,2)
366 258503231 : do ifft=1,cplex*nfft
367 258484850 : vresid1(ifft,ispden)=vhartr1_(ifft)+vxc1_(ifft,ispden)+vpsp1(ifft)-vtrial1(ifft,ispden)
368 : end do
369 : end do
370 18381 : if(nspden==4)then
371 : !$OMP PARALLEL DO COLLAPSE(2)
372 645 : do ispden=3,4
373 3822669 : do ifft=1,cplex*nfft
374 3822454 : vresid1(ifft,ispden)=vxc1_(ifft,ispden)-vtrial1(ifft,ispden)
375 : end do
376 : end do
377 : end if
378 :
379 18381 : if (ipert==natom+5.or.ipert==natom+6.or.(ipert>natom+11.and.ipert<=2*natom+11)) then
380 4746123 : vresid1 = vresid1 + v1hspinfield
381 : end if
382 :
383 18381 : if (abs(magpen) > tol6) then
384 891135 : vresid1 = vresid1 + vmagpen1
385 : end if
386 :
387 : ! Compute square norm vres2 of potential residual vresid
388 18381 : call sqnorm_v(cplex,nfft,vres2,nspden,optres,vresid1)
389 :
390 : else
391 :
392 : ! ------ Produce new value of trial potential-------------
393 : ! (only if requested ; if optres==1)
394 :
395 : !$OMP PARALLEL DO COLLAPSE(2)
396 5706 : do ispden=1,min(nspden,2)
397 33999304 : do ifft=1,cplex*nfft
398 33996671 : vtrial1(ifft,ispden)=vhartr1_(ifft)+vxc1_(ifft,ispden)+vpsp1(ifft)+vmagpen1(ifft,ispden)
399 : end do
400 : end do
401 2633 : if(nspden==4)then
402 : !$OMP PARALLEL DO COLLAPSE(2)
403 81 : do ispden=3,4
404 395127 : do ifft=1,cplex*nfft
405 395100 : vtrial1(ifft,ispden)=vxc1_(ifft,ispden)+vmagpen1(ifft,ispden)
406 : end do
407 : end do
408 : end if
409 :
410 2633 : if (ipert==natom+5.or.ipert==natom+6.or.(ipert>natom+11.and.ipert<=2*natom+11)) then
411 490808 : vtrial1 = vtrial1 + v1hspinfield
412 : end if
413 :
414 2633 : if (abs(magpen) > tol6) then
415 184375 : vtrial1 = vtrial1 + vmagpen1
416 : end if
417 :
418 : end if
419 :
420 : !Release temporary memory space
421 21014 : if (.not.vhartr1_allocated) then
422 0 : ABI_FREE(vhartr1_)
423 : end if
424 21014 : if (.not.vxc1_allocated) then
425 2296 : ABI_FREE(vxc1_)
426 : end if
427 :
428 21014 : if (ipert==natom+5.or.ipert==natom+6.or.(ipert>natom+11.and.ipert<=2*natom+11)) then
429 99 : ABI_FREE(v1hspinfield)
430 : end if
431 :
432 21014 : ABI_FREE(vmagpen1)
433 21014 : if((ipert>natom+11.and.ipert<=2*natom+11).or.abs(magpen)>tol6) then
434 62 : ABI_FREE(intgden)
435 62 : ABI_FREE(rhomag)
436 62 : ABI_FREE(fatsph)
437 62 : ABI_FREE(taumr)
438 : endif
439 :
440 21014 : call timab(157,2,tsec)
441 :
442 42028 : end subroutine dfpt_rhotov
443 : !!***
444 :
445 : !!****f* ABINIT/dfpt_v1hspinfield
446 : !! NAME
447 : !! dfpt_v1hspinfield
448 : !!
449 : !! FUNCTION
450 : !! Calculate 1st order spin magnetic (hspinfield) potential = -vec{\sigma}.\vec{b}, where
451 : !! sigma is the vector of Pauli matrices and \vec{b} is the unit
452 : !! vector indicating the perturbing field direction.
453 : !!
454 : !! INPUTS
455 : !! nspden = number of density matrix components
456 : !! nfft = numbder of fft grid points
457 : !! cplex = complex or real density matrix
458 : !! idir = direction of the perturbing field in Cartesian frame
459 : !! 1: along x
460 : !! 2: along y
461 : !! 3: along z
462 : !! 4: identity matrix at each fft point is returned (for density-density response)
463 : !!
464 : !! OUTPUT
465 : !! v1hspinfield(nfft*cplex,nspden)= 1st order spin magnetic (hspinfield) potential, or Identity matrix (electrostatic potential) for idir=4
466 : !!
467 : !! SIDE EFFECTS
468 : !!
469 : !! NOTES
470 : !! The definition of components of the potential matrix differ depending on cplex
471 : !! for nspden=4:
472 : !! For cplex=1, the potential is defined as (V_upup,V_dndn,Re[V_updn],Im[V_updn])
473 : !! For cplex=2, the definition is (V_upup,V_dndn,V_updn,i.V_updn)
474 : !!
475 : !! SOURCE
476 :
477 82 : subroutine dfpt_v1hspinfield(nspden,nfft,cplex,idir,v1hspinfield)
478 :
479 : !Arguments ------------------------------------
480 : integer , intent(in) :: idir,nfft,cplex,nspden
481 : real(dp), intent(inout) :: v1hspinfield(cplex*nfft,nspden)
482 :
483 : !Local variables-------------------------------
484 : integer :: ifft
485 : !character(len=500) :: msg
486 :
487 : ! *************************************************************************
488 :
489 : DBG_ENTER("COLL")
490 :
491 : ! if (option/=1 .and. option/=2 ) then
492 : ! write(msg,'(3a,i0)')&
493 : !& 'The argument option should be 1 or 2,',ch10,&
494 : !& 'however, option=',option
495 : ! ABI_BUG(msg)
496 : ! end if
497 : !
498 : ! if (sizein<1) then
499 : ! write(msg,'(3a,i0)')&
500 : !& 'The argument sizein should be a positive number,',ch10,&
501 : !& 'however, sizein=',sizein
502 : ! ABI_ERROR(msg)
503 : ! end if
504 :
505 : DBG_EXIT("COLL")
506 :
507 24 : select case(cplex)
508 : case(1)
509 24 : if (nspden==4) then
510 2 : if(idir==3)then ! Zeeman field along the 3rd axis (z)
511 0 : v1hspinfield(:,1)=-0.5d0
512 0 : v1hspinfield(:,2)=+0.5d0
513 0 : v1hspinfield(:,3)= 0.0d0
514 0 : v1hspinfield(:,4)= 0.0d0
515 2 : else if(idir==2)then ! Zeeman field along the 2nd axis (y)
516 0 : v1hspinfield(:,1)= 0.0d0
517 0 : v1hspinfield(:,2)= 0.0d0
518 0 : v1hspinfield(:,3)= 0.0d0
519 0 : v1hspinfield(:,4)=+0.5d0
520 2 : else if(idir==1)then ! Zeeman field along the 1st axis (x)
521 5122 : v1hspinfield(:,1)= 0.0d0
522 5122 : v1hspinfield(:,2)= 0.0d0
523 5122 : v1hspinfield(:,3)=-0.5d0
524 5122 : v1hspinfield(:,4)= 0.0d0
525 0 : else if(idir==4)then ! Scalar potential
526 0 : v1hspinfield(:,1)=-1.0d0
527 0 : v1hspinfield(:,2)=-1.0d0
528 0 : v1hspinfield(:,3)= 0.0d0
529 0 : v1hspinfield(:,4)= 0.0d0
530 : end if
531 22 : else if (nspden==2) then
532 22 : if (idir==4) then
533 0 : v1hspinfield(:,1)=-0.5d0
534 0 : v1hspinfield(:,2)=-0.5d0
535 : else
536 176022 : v1hspinfield(:,1)=-0.5d0
537 176022 : v1hspinfield(:,2)= 0.5d0
538 : end if
539 : else
540 0 : v1hspinfield(:,1)= 0.0d0
541 : end if
542 : case(2)
543 82 : if (nspden==2) then
544 0 : if (idir==4) then
545 0 : do ifft=1,nfft
546 0 : v1hspinfield(2*ifft-1,1) =-0.5e0
547 0 : v1hspinfield(2*ifft ,1) = 0.0e0
548 0 : v1hspinfield(2*ifft-1,2) =-0.5e0
549 0 : v1hspinfield(2*ifft ,2) = 0.0e0
550 : end do
551 : else
552 0 : do ifft=1,nfft
553 0 : v1hspinfield(2*ifft-1,1) =-0.5e0
554 0 : v1hspinfield(2*ifft ,1) = 0.0e0
555 0 : v1hspinfield(2*ifft-1,2) = 0.5e0
556 0 : v1hspinfield(2*ifft ,2) = 0.0e0
557 : end do
558 : end if
559 58 : else if (nspden==4) then
560 58 : select case(idir)
561 : case(1) !along x, v1=-sigma_x
562 177685 : do ifft=1,nfft
563 177664 : v1hspinfield(2*ifft-1,1)= 0.0e0 !Re[V^11]
564 177664 : v1hspinfield(2*ifft ,1)= 0.0e0 !Im[V^11]
565 177664 : v1hspinfield(2*ifft-1,2)= 0.0e0 !Re[V^22]
566 177664 : v1hspinfield(2*ifft ,2)= 0.0e0 !Im[V^22]
567 177664 : v1hspinfield(2*ifft-1,3)=-0.5e0 !Re[V^12]
568 177664 : v1hspinfield(2*ifft ,3)= 0.0e0 !Im[V^12]
569 177664 : v1hspinfield(2*ifft-1,4)= 0.0e0 !Re[i.V^21]=Im[V^12]
570 177685 : v1hspinfield(2*ifft ,4)=-0.5e0 !Im[i.V^21]=Re[V^12]
571 : end do
572 : case(2) !along y, v1 = -sigma_y
573 152075 : do ifft=1,nfft
574 152064 : v1hspinfield(2*ifft-1,1)= 0.0e0 !Re[V^11]
575 152064 : v1hspinfield(2*ifft ,1)= 0.0e0 !Im[V^11]
576 152064 : v1hspinfield(2*ifft-1,2)= 0.0e0 !Re[V^22]
577 152064 : v1hspinfield(2*ifft ,2)= 0.0e0 !Im[V^22]
578 152064 : v1hspinfield(2*ifft-1,3)= 0.0e0 !Re[V^12]
579 152064 : v1hspinfield(2*ifft ,3)=+0.5e0 !Im[V^12]
580 152064 : v1hspinfield(2*ifft-1,4)=+0.5e0 !Re[i.V^21]=Im[V^12]
581 152075 : v1hspinfield(2*ifft ,4)= 0.0e0 !Im[i.V^21]=Re[V^12]
582 : end do
583 : case(3)
584 207375 : do ifft=1,nfft
585 207360 : v1hspinfield(2*ifft-1,1)=-0.5e0 !Re[V^11]
586 207360 : v1hspinfield(2*ifft ,1)= 0.0e0 !Im[V^11]
587 207360 : v1hspinfield(2*ifft-1,2)= 0.5e0 !Re[V^22]
588 207360 : v1hspinfield(2*ifft ,2)= 0.0e0 !Im[V^22]
589 207360 : v1hspinfield(2*ifft-1,3)= 0.0e0 !Re[V^12]
590 207360 : v1hspinfield(2*ifft ,3)= 0.0e0 !Im[V^12]
591 207360 : v1hspinfield(2*ifft-1,4)= 0.0e0 !Re[i.V^21]
592 207375 : v1hspinfield(2*ifft ,4)= 0.0e0 !Im[i.V^21]
593 : end do
594 : case(4)
595 33861 : do ifft=1,nfft
596 33792 : v1hspinfield(2*ifft-1,1)=-0.5e0 !Re[V^11]
597 33792 : v1hspinfield(2*ifft ,1)= 0.0e0 !Im[V^11]
598 33792 : v1hspinfield(2*ifft-1,2)=-0.5e0 !Re[V^22]
599 33792 : v1hspinfield(2*ifft ,2)= 0.0e0 !Im[V^22]
600 33792 : v1hspinfield(2*ifft-1,3)= 0.0e0 !Re[V^12]
601 33792 : v1hspinfield(2*ifft ,3)= 0.0e0 !Im[V^12]
602 33792 : v1hspinfield(2*ifft-1,4)= 0.0e0 !Re[i.V^21]
603 33803 : v1hspinfield(2*ifft ,4)= 0.0e0 !Im[i.V^21]
604 : end do
605 : end select
606 : end if
607 : end select !cplex
608 :
609 82 : end subroutine dfpt_v1hspinfield
610 : !!***
611 :
612 : !!****f* ABINIT/dfpt_v1magpen
613 : !! NAME
614 : !! dfpt_v1magpen
615 : !!
616 : !! FUNCTION
617 : !! Calculate a 1st order penalty potential and energy in order to
618 : !! harden the spin degrees of freedom.
619 : !!
620 : !! INPUTS
621 : !! cplex = complex or real density matrix
622 : !! magpen=energy shift to apply on the spin degrees of freedom
623 : !! mpatpol(2)=initial and final atomic positions whose local magnetic moments will be penalized
624 : !! mpdir(3)=directions of the magnetic moments to be penalized
625 : !! mpi_enreg=information about MPI parallelization
626 : !! nfft = numbder of fft grid points
627 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
628 : !! nspden = number of density matrix components
629 : !! qphon(3)=reduced coordinates for the phonon wavelength
630 : !! taumr(nfft,natom,3)= array describing r-xred(iatom) at any point of the FFT grid
631 : !! xred(3,natom)=reduced dimensionless atomic coordinates
632 : !!
633 : !! OUTPUT
634 : !! vmagpen1(nfft*cplex,nspden)= 1st order magnetic penalty potential
635 : !! emagpen1= magnetic penalty energy contribution to second-order energy
636 : !!
637 : !! SIDE EFFECTS
638 : !!
639 : !! NOTES
640 : !! The definition of components of the potential matrix differ depending on cplex:
641 : !! For cplex=1, the potential is defined as (V_upup,V_dndn,Re[V_updn],Im[V_updn])
642 : !! For cplex=2, the definition is (V_upup,V_dndn,V_updn,i.V_updn)
643 : !!
644 : !! If magpen < 0 the penalty field is defined from the cell-integrated magnetic moments along the directions given by mpdir
645 : !! IF magpen > 0 the penalty field is defined from the atom spheres-integrated magnetic moments as given by mpatpol and mpdir
646 : !!
647 : !! SOURCE
648 :
649 62 : subroutine dfpt_v1magpen(cplex,emagpen1,fatsph,intgden,magpen,mpatpol,mpdir,&
650 62 : & natom,nfft,nspden,qphon,rhomag,taumr,vmagpen1)
651 :
652 : !Arguments
653 : !scalars:
654 : integer,intent(in) :: cplex,natom,nfft,nspden
655 : real(dp),intent(in) :: magpen
656 : real(dp),intent(out) :: emagpen1
657 : !arrays:
658 : integer,intent(in) :: mpatpol(2),mpdir(3)
659 : real(dp),intent(in) :: fatsph(nfft,natom)
660 : real(dp),intent(in) :: intgden(cplex,nspden,natom)
661 : real(dp), intent(in) :: qphon(3)
662 : real(dp),intent(in) :: rhomag(2,nspden)
663 : real(dp), intent(in) :: taumr(nfft,natom,3)
664 : real(dp),intent(out) :: vmagpen1(cplex*nfft,nspden)
665 :
666 : !Local variables-------------------------------
667 : !scalars:
668 : integer :: i,iatom,ifft,im,re
669 : real(dp) :: arg
670 : real(dp) :: phr1d_re,phr1d_im
671 : real(dp) :: Blocx_re,Blocy_re,Blocz_re
672 : real(dp) :: Blocx_im,Blocy_im,Blocz_im
673 : !arrays:
674 124 : real(dp) :: Bx(cplex),By(cplex),Bz(cplex)
675 124 : real(dp) :: Blocx(cplex*nfft),Blocy(cplex*nfft),Blocz(cplex*nfft)
676 62 : real(dp) :: rhomag_eff(2,nspden),intgden_eff(cplex,nspden,natom)
677 :
678 : ! *************************************************************************
679 :
680 176 : if (cplex==1.and.any(abs(qphon(:))>tol8)) then
681 0 : ABI_ERROR('Local Zeeman fields are cplex==2 at finite q vector')
682 : end if
683 :
684 : !Compute magnetic penalty from cell-integrated magnetic moments
685 62 : if (magpen < zero) then
686 :
687 0 : rhomag_eff(:,1)=rhomag(:,1)
688 0 : rhomag_eff(:,2:4)=half*rhomag(:,2:4) !Convert from mu_B to a.u.
689 0 : do i=1,3
690 0 : if (mpdir(i)==0) rhomag_eff(:,1+i) = zero
691 : end do
692 :
693 0 : if (cplex==1) then
694 0 : emagpen1=-one*magpen*(rhomag_eff(1,2)**2+rhomag_eff(1,3)**2+rhomag_eff(1,4)**2)
695 0 : else if (cplex==2) then
696 : emagpen1=-one*magpen*(rhomag_eff(1,2)**2+rhomag_eff(2,2)**2 &
697 : & + rhomag_eff(1,3)**2+rhomag_eff(2,3)**2 &
698 0 : & + rhomag_eff(1,4)**2+rhomag_eff(2,4)**2 )
699 : end if
700 :
701 0 : Bx(:)=-half*magpen*rhomag_eff(:,2)
702 0 : By(:)=-half*magpen*rhomag_eff(:,3)
703 0 : Bz(:)=-half*magpen*rhomag_eff(:,4)
704 0 : if (cplex==1) then
705 0 : do ifft=1,nfft
706 0 : vmagpen1(ifft,1)=Bz(1)
707 0 : vmagpen1(ifft,2)=-Bz(1)
708 0 : vmagpen1(ifft,3)=Bx(1)
709 0 : vmagpen1(ifft,4)=-By(1)
710 : end do
711 0 : else if (cplex==2) then
712 0 : do ifft=1,nfft
713 0 : vmagpen1(2*ifft-1,1)=Bz(1)
714 0 : vmagpen1(2*ifft ,1)=Bz(2)
715 0 : vmagpen1(2*ifft-1,2)=-Bz(1)
716 0 : vmagpen1(2*ifft ,2)=-Bz(2)
717 0 : vmagpen1(2*ifft-1,3)=Bx(1)+By(2)
718 0 : vmagpen1(2*ifft ,3)=Bx(2)-By(1)
719 0 : vmagpen1(2*ifft-1,4)=-Bx(2)-By(1)
720 0 : vmagpen1(2*ifft ,4)=Bx(1)-By(2)
721 : end do
722 : end if
723 :
724 : !Compute magnetic penalty from atom shperes-integrated magnetic moments
725 62 : else if (magpen > zero) then
726 62 : emagpen1=zero
727 268862 : Blocx=zero
728 268862 : Blocy=zero
729 268862 : Blocz=zero
730 :
731 1064 : intgden_eff(:,1,:)=intgden(:,1,:)
732 3440 : intgden_eff(:,2:4,:)=half*intgden(:,2:4,:) !Convert from mu_B to a.u.
733 186 : do iatom=mpatpol(1),mpatpol(2)
734 :
735 496 : do i=1,3
736 706 : if (mpdir(i)==0) intgden_eff(:,1+i,iatom) = zero
737 : end do
738 :
739 124 : if (cplex==1) then
740 : emagpen1=emagpen1+one*magpen*(intgden_eff(1,2,iatom)**2+ &
741 : & intgden_eff(1,3,iatom)**2+ &
742 38 : & intgden_eff(1,4,iatom)**2)
743 86 : else if (cplex==2) then
744 : emagpen1=emagpen1+one*magpen*(intgden_eff(1,2,iatom)**2+intgden_eff(2,2,iatom)**2 &
745 : & + intgden_eff(1,3,iatom)**2+intgden_eff(2,3,iatom)**2 &
746 86 : & + intgden_eff(1,4,iatom)**2+intgden_eff(2,4,iatom)**2 )
747 : end if
748 :
749 186 : if (cplex==1) then
750 97318 : do ifft=1,nfft
751 97280 : Blocx(ifft)=Blocx(ifft)+half*magpen*intgden_eff(1,2,iatom)*fatsph(ifft,iatom)
752 97280 : Blocy(ifft)=Blocy(ifft)+half*magpen*intgden_eff(1,3,iatom)*fatsph(ifft,iatom)
753 97318 : Blocz(ifft)=Blocz(ifft)+half*magpen*intgden_eff(1,4,iatom)*fatsph(ifft,iatom)
754 : end do
755 344 : else if (cplex==2.and.sum(qphon(:)**2) < tol8) then
756 97318 : do ifft=1,nfft
757 97280 : Blocx(2*ifft-1)=Blocx(2*ifft-1)+half*magpen*intgden_eff(1,2,iatom)*fatsph(ifft,iatom)
758 97280 : Blocy(2*ifft-1)=Blocy(2*ifft-1)+half*magpen*intgden_eff(1,3,iatom)*fatsph(ifft,iatom)
759 97280 : Blocz(2*ifft-1)=Blocz(2*ifft-1)+half*magpen*intgden_eff(1,4,iatom)*fatsph(ifft,iatom)
760 97280 : Blocx(2*ifft)=Blocx(2*ifft)+half*magpen*intgden_eff(2,2,iatom)*fatsph(ifft,iatom)
761 97280 : Blocy(2*ifft)=Blocy(2*ifft)+half*magpen*intgden_eff(2,3,iatom)*fatsph(ifft,iatom)
762 97318 : Blocz(2*ifft)=Blocz(2*ifft)+half*magpen*intgden_eff(2,4,iatom)*fatsph(ifft,iatom)
763 : end do
764 192 : else if (cplex==2.and.sum(qphon(:)**2) > tol8) then
765 122928 : do ifft=1,nfft
766 122880 : re=2*ifft-1
767 122880 : im=2*ifft
768 :
769 122880 : Blocx_re=+half*magpen*intgden_eff(1,2,iatom)*fatsph(ifft,iatom)
770 122880 : Blocy_re=+half*magpen*intgden_eff(1,3,iatom)*fatsph(ifft,iatom)
771 122880 : Blocz_re=+half*magpen*intgden_eff(1,4,iatom)*fatsph(ifft,iatom)
772 122880 : Blocx_im=+half*magpen*intgden_eff(2,2,iatom)*fatsph(ifft,iatom)
773 122880 : Blocy_im=+half*magpen*intgden_eff(2,3,iatom)*fatsph(ifft,iatom)
774 122880 : Blocz_im=+half*magpen*intgden_eff(2,4,iatom)*fatsph(ifft,iatom)
775 :
776 491520 : arg=two_pi*dot_product(qphon,-taumr(ifft,iatom,:))
777 122880 : phr1d_re=dcos(arg)
778 122880 : phr1d_im=dsin(arg)
779 :
780 122880 : Blocx(re)= Blocx(re)+phr1d_re*Blocx_re-phr1d_im*Blocx_im
781 122880 : Blocx(im)= Blocx(im)+phr1d_im*Blocx_re+phr1d_re*Blocx_im
782 122880 : Blocy(re)= Blocy(re)+phr1d_re*Blocy_re-phr1d_im*Blocy_im
783 122880 : Blocy(im)= Blocy(im)+phr1d_im*Blocy_re+phr1d_re*Blocy_im
784 122880 : Blocz(re)= Blocz(re)+phr1d_re*Blocz_re-phr1d_im*Blocz_im
785 122928 : Blocz(im)= Blocz(im)+phr1d_im*Blocz_re+phr1d_re*Blocz_im
786 : end do
787 : end if
788 :
789 : end do !iatom
790 :
791 62 : if (cplex==1) then
792 48659 : do ifft=1,nfft
793 48640 : vmagpen1(ifft,1)=Blocz(ifft)
794 48640 : vmagpen1(ifft,2)=-Blocz(ifft)
795 48640 : vmagpen1(ifft,3)=Blocx(ifft)
796 48659 : vmagpen1(ifft,4)=-Blocy(ifft)
797 : end do
798 43 : else if (cplex==2) then
799 110123 : do ifft=1,nfft
800 110080 : vmagpen1(2*ifft-1,1)=Blocz(2*ifft-1)
801 110080 : vmagpen1(2*ifft ,1)=Blocz(2*ifft)
802 110080 : vmagpen1(2*ifft-1,2)=-Blocz(2*ifft-1)
803 110080 : vmagpen1(2*ifft ,2)=-Blocz(2*ifft)
804 110080 : vmagpen1(2*ifft-1,3)=Blocx(2*ifft-1)+Blocy(2*ifft)
805 110080 : vmagpen1(2*ifft ,3)=Blocx(2*ifft)-Blocy(2*ifft-1)
806 110080 : vmagpen1(2*ifft-1,4)=-Blocx(2*ifft)-Blocy(2*ifft-1)
807 110123 : vmagpen1(2*ifft ,4)=Blocx(2*ifft-1)-Blocy(2*ifft)
808 : end do
809 : end if
810 :
811 : end if
812 :
813 62 : end subroutine dfpt_v1magpen
814 : !!***
815 :
816 : !!****f* ABINIT/dfpt_v1hspinfield_atsph
817 : !! NAME
818 : !! dfpt_v1hspinfield_atsph
819 : !!
820 : !! FUNCTION
821 : !! Calculate 1st order potential due to a local Zeeman field inside an
822 : !! atom centered sphere= -vec{\sigma}.\vec{b}*f_i(r), where
823 : !! sigma is the vector of Pauli matrices, \vec{b}(r) is the unit
824 : !! vector indicating the perturbing field direction and f_i(r) is the
825 : !! real-space function defining the sphere around atom i.
826 : !!
827 : !! INPUTS
828 : !! nspden = number of density matrix components
829 : !! nfft = numbder of fft grid points
830 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
831 : !! cplex = complex or real density matrix
832 : !! fatsph(nfft,natom)= functions defining the atomic spheres of integration in real space
833 : !! idir = direction of the perturbing field in Cartesian frame
834 : !! 1: along x
835 : !! 2: along y
836 : !! 3: along z
837 : !! 4: identity matrix at each fft point is returned (for density-density response)
838 : !! qphon(3)=reduced coordinates for the phonon wavelength
839 : !! taumr(nfft,natom,3)= array describing r-xred(iatom) at any point of the FFT grid
840 : !! xred(3,natom)=reduced dimensionless atomic coordinates
841 : !!
842 : !! OUTPUT
843 : !! v1hspinfield(nfft*cplex,nspden)= 1st order Zeeman potential, or Identity matrix (electrostatic potential) for idir=4
844 : !!
845 : !! SIDE EFFECTS
846 : !!
847 : !! NOTES
848 : !! The definition of components of the potential matrix differ depending on cplex
849 : !! for nspden=4:
850 : !! For cplex=1, the potential is defined as (V_upup,V_dndn,Re[V_updn],Im[V_updn])
851 : !! For cplex=2, the definition is (V_upup,V_dndn,V_updn,i.V_updn)
852 : !!
853 : !! SOURCE
854 :
855 17 : subroutine dfpt_v1hspinfield_atsph(cplex,fatsph,idir,ipert,natom,nfft,nspden,&
856 17 : & qphon,taumr,v1hspinfield)
857 :
858 : !Arguments ------------------------------------
859 : !scalars
860 : integer, intent(in) :: idir,ipert,nfft,cplex,natom,nspden
861 : !arrays
862 : real(dp), intent(in) :: fatsph(nfft,natom)
863 : real(dp), intent(in) :: qphon(3)
864 : real(dp), intent(inout):: v1hspinfield(cplex*nfft,nspden)
865 : real(dp), intent(in) :: taumr(nfft,natom,3)
866 :
867 : !Local variables-------------------------------
868 : !scalars
869 : integer :: ifft,iatom,im,re
870 : real(dp) :: arg
871 : real(dp) :: phr1d_re,phr1d_im
872 : real(dp) :: Bloc_re, Bloc_im
873 : character(len=500) :: msg
874 : !arrays
875 17 : real(dp) :: Bloc(cplex*nfft)
876 :
877 : ! *************************************************************************
878 :
879 47 : if (cplex==1.and.any(abs(qphon(:))>tol8)) then
880 0 : ABI_ERROR('Local Zeeman fields are cplex==2 at finite q vector')
881 : end if
882 :
883 17 : iatom=ipert-natom-11
884 :
885 : !Define the local magnetic field
886 17 : if (cplex==1) then
887 12805 : do ifft=1,nfft
888 12805 : Bloc(ifft)=-half*fatsph(ifft,iatom)
889 : end do
890 12 : else if (cplex==2) then
891 30732 : do ifft=1,nfft
892 30720 : re=2*ifft-1
893 30720 : im=2*ifft
894 122892 : if (sum(qphon(:)**2)<tol8) then
895 12800 : Bloc(re)=-half*fatsph(ifft,iatom)
896 12800 : Bloc(im)=zero
897 : else
898 17920 : Bloc_re=-half*fatsph(ifft,iatom)
899 17920 : Bloc_im=zero
900 71680 : arg=two_pi*dot_product(qphon,-taumr(ifft,iatom,:))
901 17920 : phr1d_re=dcos(arg)
902 17920 : phr1d_im=dsin(arg)
903 17920 : Bloc(re)=phr1d_re*Bloc_re-phr1d_im*Bloc_im
904 17920 : Bloc(im)=phr1d_im*Bloc_re+phr1d_re*Bloc_im
905 : end if
906 : end do
907 : end if
908 :
909 : !Build the first-order potential
910 5 : select case(cplex)
911 : case(1)
912 5 : if (nspden==4) then
913 5 : if(idir==3)then ! Zeeman field along the 3rd axis (z)
914 0 : do ifft=1,nfft
915 0 : v1hspinfield(ifft,1)=Bloc(ifft)
916 0 : v1hspinfield(ifft,2)=-Bloc(ifft)
917 0 : v1hspinfield(ifft,3)= 0.0d0
918 0 : v1hspinfield(ifft,4)= 0.0d0
919 : end do
920 5 : else if(idir==2)then ! Zeeman field along the 2nd axis (y)
921 0 : do ifft=1,nfft
922 0 : v1hspinfield(ifft,1)= 0.0d0
923 0 : v1hspinfield(ifft,2)= 0.0d0
924 0 : v1hspinfield(ifft,3)= 0.0d0
925 0 : v1hspinfield(ifft,4)=-Bloc(ifft)
926 : end do
927 : else ! Zeeman field along the 1st axis (x)
928 12805 : do ifft=1,nfft
929 12800 : v1hspinfield(ifft,1)= 0.0d0
930 12800 : v1hspinfield(ifft,2)= 0.0d0
931 12800 : v1hspinfield(ifft,3)=Bloc(ifft)
932 12805 : v1hspinfield(ifft,4)= 0.0d0
933 : end do
934 : end if
935 : else
936 0 : write(msg,*) 'Response to local Zeeman fields only implemented for nspden=4'
937 0 : ABI_BUG(msg)
938 : end if
939 : case(2)
940 17 : if (nspden==4) then
941 12 : select case(idir)
942 : case(1) !along x, v1=-sigma_x
943 30732 : do ifft=1,nfft
944 30720 : v1hspinfield(2*ifft-1,1)= 0.0e0 !Re[V^11]
945 30720 : v1hspinfield(2*ifft ,1)= 0.0e0 !Im[V^11]
946 30720 : v1hspinfield(2*ifft-1,2)= 0.0e0 !Re[V^22]
947 30720 : v1hspinfield(2*ifft ,2)= 0.0e0 !Im[V^22]
948 30720 : v1hspinfield(2*ifft-1,3)= Bloc(2*ifft-1) !Re[V^12]
949 30720 : v1hspinfield(2*ifft ,3)= Bloc(2*ifft) !Im[V^12]
950 30720 : v1hspinfield(2*ifft-1,4)=-Bloc(2*ifft) !Re[i.V^21]=Im[V^12]
951 30732 : v1hspinfield(2*ifft ,4)= Bloc(2*ifft-1) !Im[i.V^21]=Re[V^12]
952 : end do
953 : case(2) !along y, v1 = -sigma_y
954 0 : do ifft=1,nfft
955 0 : v1hspinfield(2*ifft-1,1)= 0.0e0 !Re[V^11]
956 0 : v1hspinfield(2*ifft ,1)= 0.0e0 !Im[V^11]
957 0 : v1hspinfield(2*ifft-1,2)= 0.0e0 !Re[V^22]
958 0 : v1hspinfield(2*ifft ,2)= 0.0e0 !Im[V^22]
959 0 : v1hspinfield(2*ifft-1,3)= Bloc(2*ifft) !Re[V^12]
960 0 : v1hspinfield(2*ifft ,3)=-Bloc(2*ifft-1) !Im[V^12]
961 0 : v1hspinfield(2*ifft-1,4)=-Bloc(2*ifft-1) !Re[i.V^21]=Im[V^12]
962 0 : v1hspinfield(2*ifft ,4)=-Bloc(2*ifft) !Im[i.V^21]=Re[V^12]
963 : end do
964 : case(3)
965 12 : do ifft=1,nfft
966 0 : v1hspinfield(2*ifft-1,1)= Bloc(2*ifft-1) !Re[V^11]
967 0 : v1hspinfield(2*ifft ,1)= Bloc(2*ifft) !Im[V^11]
968 0 : v1hspinfield(2*ifft-1,2)=-Bloc(2*ifft-1) !Re[V^22]
969 0 : v1hspinfield(2*ifft ,2)=-Bloc(2*ifft) !Im[V^22]
970 0 : v1hspinfield(2*ifft-1,3)= 0.0e0 !Re[V^12]
971 0 : v1hspinfield(2*ifft ,3)= 0.0e0 !Im[V^12]
972 0 : v1hspinfield(2*ifft-1,4)= 0.0e0 !Re[i.V^21]
973 0 : v1hspinfield(2*ifft ,4)= 0.0e0 !Im[i.V^21]
974 : end do
975 : end select
976 : else
977 0 : write(msg,*) 'Response to local spin fields only implemented for nspden=4'
978 0 : ABI_BUG(msg)
979 : end if
980 : end select !cplex
981 :
982 17 : end subroutine dfpt_v1hspinfield_atsph
983 : !!***
984 :
985 : end module m_dfpt_rhotov
986 : !!***
|
