Branch data 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 : 21140 : subroutine dfpt_rhotov(cplex,ehart01,ehart1,elmag1,elpsp1,emagpen1,exc1,gsqcut,icutcoul,idir,ipert,&
124 : 21140 : & ixc,kxc,magpen,mpatpol,mpdir,mpi_enreg,natom,nfft,ngfft,nhat,nhat1,nhat1gr,nhat1grdim,nkxc,nspden,ntypat,n3xccc,&
125 : 21140 : & non_magnetic_xc,optene,optres,qphon,ratsm,ratsph,rhog,rhog1,rhor,rhor1,rprimd,typat,ucvol,&
126 [ - + - + ]: 21140 : & 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 : 21140 : real(dp),allocatable :: rhor1_nohat(:,:),vhartr01(:),vxc1val(:,:)
167 : 21140 : real(dp),pointer :: rhor1_(:,:),vhartr1_(:),vxc1_(:,:),v1hspinfield(:,:)
168 : 21140 : real(dp),allocatable :: fatsph(:,:),intgden(:,:,:),rhomag(:,:),vmagpen1(:,:)
169 : 21140 : real(dp),allocatable :: taumr(:,:,:)
170 : :
171 : : ! *********************************************************************
172 : :
173 : 21140 : call timab(157,1,tsec)
174 : :
175 : : !FR EB SPr
176 [ + + ]: 21140 : 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 : 21140 : nfftot=ngfft(1)*ngfft(2)*ngfft(3)
184 : :
185 : : !Eventually allocate temporary memory space
186 : 21140 : vhartr1_allocated=(size(vhartr1)>0)
187 [ + - ]: 21140 : if (vhartr1_allocated) then
188 : 21140 : vhartr1_ => vhartr1
189 : : else
190 [ # # # # : 0 : ABI_MALLOC(vhartr1_,(cplex*nfft))
# # # # #
# # # ]
191 : : end if
192 [ + + ]: 63420 : vxc1_allocated=(size(vxc1)>0)
193 [ + + ]: 21140 : if (vxc1_allocated) then
194 : 18802 : vxc1_ => vxc1
195 : : else
196 [ - + + - : 9352 : ABI_MALLOC(vxc1_,(cplex*nfft,nspden))
- + + - +
- + - - +
- + ]
197 : : end if
198 : :
199 : : !If needed, store pseudo density without charge compensation
200 [ + + + + ]: 21140 : if (usepaw==1.and.usexcnhat==0) then
201 [ - + + - : 7508 : ABI_MALLOC(rhor1_,(cplex*nfft,nspden))
- + + - +
- + - - +
- + ]
202 [ + - + + : 302952418 : rhor1_(:,:)=rhor1(:,:)-nhat1(:,:)
+ + + + +
+ ]
203 : : else
204 : 19263 : rhor1_ => rhor1
205 : : end if
206 : :
207 : : !Uniform Zeeman or scalar potential
208 [ + + + + ]: 21140 : 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 [ + + - + : 21140 : 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 [ + + + - ]: 21140 : 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 [ - + + - : 84560 : ABI_MALLOC(vmagpen1,(cplex*nfft,nspden))
- + + - +
- + - - +
+ - ]
236 [ + + + + ]: 390010111 : vmagpen1=zero
237 [ + + ]: 21140 : 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 [ - + ]: 21140 : call hartre(cplex,gsqcut,icutcoul,0,mpi_enreg,nfft,ngfft,1,zero,rhog1,rprimd,vcutgeo,vhartr1_,qpt=qphon)
244 : :
245 [ + + ]: 21140 : if (optene>0) then
246 [ - + ]: 18617 : call dotprod_vn(cplex,rhor1,ehart1,doti,nfft,nfftot,1,1,vhartr1_,ucvol)
247 : : end if
248 : :
249 [ + + ]: 21140 : if (optene>0) ehart01=zero
250 [ + + + + ]: 21140 : 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 [ + + ]: 21140 : option=0;if (optene==0) option=1
269 [ + + + + ]: 21140 : 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 [ + + ]: 19465 : 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 [ - + ]: 19223 : & 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 [ + + ]: 21140 : if (optene>0) then
292 [ + + ]: 18617 : 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 [ + + ]: 3558 : if (usexcnhat/=0) then
303 [ - + + - : 6357 : ABI_MALLOC(rhor1_nohat,(cplex*nfft,1))
+ - + - -
+ + - ]
304 [ + + ]: 33282407 : rhor1_nohat(:,1)=rhor1(:,1)-nhat1(:,1)
305 [ - + ]: 2119 : call dotprod_vn(cplex,rhor1 ,elpsp10,doti,nfft,nfftot,nspden,1,vxc1_,ucvol)
306 : 2119 : call dotprod_vn(cplex,rhor1_nohat,elpsp1 ,doti,nfft,nfftot,1 ,1,vpsp1,ucvol)
307 [ - + ]: 2119 : ABI_FREE(rhor1_nohat)
308 : : else
309 [ - + - + ]: 1439 : call dotprod_vn(cplex,rhor1_,elpsp10,doti,nfft,nfftot,nspden,1,vxc1_,ucvol)
310 [ - + ]: 1439 : 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 : 18617 : elpsp1=two*(elpsp1+elpsp10)
316 : : end if
317 : :
318 : :
319 : : !Compute XC valence contribution exc1 and complete eventually Vxc^(1)
320 [ + + ]: 21140 : if (optene>0) then
321 [ - + + - : 55851 : ABI_MALLOC(vxc1val,(cplex*nfft,nspden))
- + + - -
+ + - ]
322 [ + + + + ]: 345096233 : vxc1val=zero
323 : 18617 : option=2
324 : : !FR SPr EB non-collinear magnetism
325 [ + + ]: 18617 : 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 : 18402 : & non_magnetic_xc,nspden,n3xccc,option,qphon,rhor1,rprimd,usexcnhat,vxc1val,xccc3d1)
333 : : end if !nspden==4
334 : :
335 [ + + + + ]: 345096233 : vxc1_(:,:)=vxc1_(:,:)+vxc1val(:,:)
336 [ - + ]: 18617 : call dotprod_vn(cplex,rhor1_,exc1,doti,nfft,nfftot,nspden,1,vxc1val,ucvol)
337 [ - + ]: 18617 : ABI_FREE(vxc1val)
338 : : end if
339 : :
340 [ + + + + ]: 21140 : if (usepaw==1.and.usexcnhat==0) then
341 [ - + ]: 1877 : 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 [ + + ]: 21140 : if (optres==0) then
364 : : !$OMP PARALLEL DO COLLAPSE(2)
365 [ + + ]: 38778 : do ispden=1,min(nspden,2)
366 [ + + ]: 339350855 : do ifft=1,cplex*nfft
367 : 339332362 : vresid1(ifft,ispden)=vhartr1_(ifft)+vxc1_(ifft,ispden)+vpsp1(ifft)-vtrial1(ifft,ispden)
368 : : end do
369 : : end do
370 [ + + ]: 18493 : 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 [ + + + + : 18493 : 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 [ + + ]: 18493 : if (abs(magpen) > tol6) then
384 [ + + + + ]: 891135 : vresid1 = vresid1 + vmagpen1
385 : : end if
386 : :
387 : : ! Compute square norm vres2 of potential residual vresid
388 : 18493 : 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 [ + + ]: 5736 : do ispden=1,min(nspden,2)
397 [ + + ]: 46441702 : do ifft=1,cplex*nfft
398 : 46439055 : vtrial1(ifft,ispden)=vhartr1_(ifft)+vxc1_(ifft,ispden)+vpsp1(ifft)
399 : : end do
400 : : end do
401 [ + + ]: 2647 : 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)
406 : : end do
407 : : end do
408 : : end if
409 : :
410 [ + + + + : 2647 : 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 [ + + ]: 2647 : if (abs(magpen) > tol6) then
415 [ + + + + ]: 184375 : vtrial1 = vtrial1 + vmagpen1
416 : : end if
417 : :
418 : : end if
419 : :
420 : : !Release temporary memory space
421 [ - + ]: 21140 : if (.not.vhartr1_allocated) then
422 [ # # ]: 0 : ABI_FREE(vhartr1_)
423 : : end if
424 [ + + ]: 21140 : if (.not.vxc1_allocated) then
425 [ - + ]: 2338 : ABI_FREE(vxc1_)
426 : : end if
427 : :
428 [ + + + + : 21140 : 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 [ - + ]: 21140 : ABI_FREE(vmagpen1)
433 [ + + - + : 21140 : 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 : 21140 : call timab(157,2,tsec)
441 : :
442 [ - + - + : 42280 : 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 : : !!***
|
