Branch data Line data Source code
1 : : !!****m* ABINIT/m_mklocl
2 : : !! NAME
3 : : !! m_mklocl
4 : : !!
5 : : !! FUNCTION
6 : : !! Routines related to the local part of the pseudopotentials.
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GMR, MM, DRH)
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_mklocl
23 : :
24 : : use defs_basis
25 : : use defs_wvltypes
26 : : use m_abicore
27 : : use m_errors
28 : : use m_xmpi
29 : : use m_dtset
30 : :
31 : : use defs_datatypes, only : pseudopotential_type
32 : : use defs_abitypes, only : MPI_type
33 : : use m_time, only : timab
34 : : use m_geometry, only : xred2xcart
35 : : use m_mpinfo, only : ptabs_fourdp
36 : : use m_pawtab, only : pawtab_type
37 : : use m_mklocl_realspace, only : mklocl_realspace, mklocl_wavelets
38 : : use m_fft, only : fourdp
39 : : use m_gtermcutoff,only : termcutoff
40 : :
41 : : use m_splines, only : splfit
42 : : use m_dfpt_mkvxc, only : dfpt_mkvxcgga_n0met
43 : :
44 : : #if defined HAVE_BIGDFT
45 : : use BigDFT_API, only : ELECTRONIC_DENSITY
46 : : use m_abi2big, only : wvl_rho_abi2big
47 : : #endif
48 : :
49 : : implicit none
50 : :
51 : : private
52 : : !!***
53 : :
54 : : public :: mklocl
55 : : public :: mklocl_recipspace
56 : : public :: dfpt_vlocal ! Local part of 1st-order potential due to atomic displacement.
57 : : public :: vlocalstr ! Compute strain derivatives of local ionic potential
58 : : public :: dfpt_vlocaldq ! Compute the first q-gradient of the 1st-order potential due to atomic displacement.
59 : : public :: dfpt_vlocaldqdq ! Compute the second q-gradient of the 1st-order potential due to atomic displacement.
60 : : public :: dfpt_vmetdqdq ! Compute the second q-gradient of the 1st-order potential due to a metric perturbation.
61 : : !!***
62 : :
63 : : contains
64 : : !!***
65 : :
66 : : !!****f* ABINIT/mklocl
67 : : !! NAME
68 : : !! mklocl
69 : : !!
70 : : !! FUNCTION
71 : : !! This method is a wrapper for mklocl_recipspace and mklocl_realspace.
72 : : !! It does some consistency checks before calling one of the two methods.
73 : : !!
74 : : !! Optionally compute :
75 : : !! option=1 : local ionic potential throughout unit cell
76 : : !! option=2 : contribution of local ionic potential to E gradient wrt xred
77 : : !! option=3 : contribution of local ionic potential to
78 : : !! stress tensor (only with reciprocal space computations)
79 : : !! option=4 : contribution of local ionic potential to
80 : : !! second derivative of E wrt xred (only with reciprocal space computations)
81 : : !!
82 : : !! INPUTS
83 : : !! if(option==3) eei=local pseudopotential part of total energy (hartree)
84 : : !! gmet(3,3)=reciprocal space metric ($\textrm{Bohr}^{-2}$).
85 : : !! gprimd(3,3)=reciprocal space dimensional primitive translations
86 : : !! gsqcut=cutoff on $|G|^2$: see setup1 for definition (doubled sphere).
87 : : !! mgfft=maximum size of 1D FFTs
88 : : !! mpi_enreg=information about MPI parallelization
89 : : !! natom=number of atoms in unit cell.
90 : : !! nattyp(ntypat)=number of atoms of each type in cell.
91 : : !! nfft=(effective) number of FFT grid points (for this processor)
92 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
93 : : !! nspden=number of spin-density components
94 : : !! ntypat=number of types of atoms.
95 : : !! option= (see above)
96 : : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
97 : : !! ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phase information.
98 : : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
99 : : !! qprtrb(3)= integer wavevector of possible perturbing potential
100 : : !! in basis of reciprocal lattice translations
101 : : !! rhog(2,nfft)=electron density rho(G) (electrons/$\textrm{Bohr}^3$)
102 : : !! (needed if option==2 or if option==4)
103 : : !! rhor(nfft,nspden)=electron density in electrons/bohr**3.
104 : : !! (needed if option==2 or if option==4)
105 : : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
106 : : !! ucvol=unit cell volume ($\textrm{Bohr}^3$).
107 : : !! vprtrb(2)=complex amplitude of possible perturbing potential; if nonzero,
108 : : !! perturbing potential is added of the form
109 : : !! $V(G)=(vprtrb(1)+I*vprtrb(2))/2$ at the values G=qprtrb and
110 : : !! $(vprtrb(1)-I*vprtrb(2))/2$ at $G=-qprtrb$ (integers)
111 : : !! xred(3,natom)=reduced dimensionless atomic coordinates
112 : : !!
113 : : !! OUTPUT
114 : : !! (if option==1) vpsp(nfft)=local crystal pseudopotential in real space.
115 : : !! (if option==2) grtn(3,natom)=grads of Etot wrt tn.
116 : : !! (if option==3) lpsstr(6)=components of local psp part of stress tensor
117 : : !! (Cartesian coordinates, symmetric tensor) in hartree/$\textrm{bohr}^3$
118 : : !! Store 6 unique components in order 11, 22, 33, 32, 31, 21
119 : : !! (if option==4) dyfrlo(3,3,natom)=d2 Eei/d tn(i)/d tn(j). (Hartrees)
120 : : !!
121 : : !! NOTES
122 : : !! Note that the present routine is tightly connected to the dfpt_vlocal.f routine,
123 : : !! that compute the derivative of the local ionic potential
124 : : !! with respect to one atomic displacement. The argument list
125 : : !! and the internal loops to be considered were sufficiently different
126 : : !! as to make the two routine different.
127 : : !!
128 : : !! SOURCE
129 : :
130 : 23956 : subroutine mklocl(dtset, dyfrlo,eei,gmet,gprimd,grtn,gsqcut,lpsstr,mgfft,&
131 : 23956 : & mpi_enreg,natom,nattyp,nfft,ngfft,nspden,ntypat,option,pawtab,ph1d,psps,qprtrb,&
132 : 23956 : & rhog,rhor,rprimd,ucvol,vprtrb,vpsp,wvl,wvl_den,xred)
133 : :
134 : : !Arguments ------------------------------------
135 : : !scalars
136 : : integer,intent(in) :: mgfft,natom,nfft,nspden,ntypat,option
137 : : real(dp),intent(in) :: eei,gsqcut,ucvol
138 : : type(MPI_type),intent(in) :: mpi_enreg
139 : : type(dataset_type),intent(in) :: dtset
140 : : type(pseudopotential_type),intent(in) :: psps
141 : : type(wvl_internal_type), intent(in) :: wvl
142 : : type(wvl_denspot_type), intent(inout) :: wvl_den
143 : : type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw)
144 : : !arrays
145 : : integer,intent(in) :: nattyp(ntypat),ngfft(18),qprtrb(3)
146 : : real(dp),intent(in) :: gmet(3,3),gprimd(3,3),ph1d(2,3*(2*mgfft+1)*natom)
147 : : real(dp),intent(in) :: rhog(2,nfft),rprimd(3,3)
148 : : real(dp),intent(in) :: vprtrb(2),xred(3,natom)
149 : : real(dp),intent(in),target :: rhor(nfft,nspden)
150 : : real(dp),intent(out) :: dyfrlo(3,3,natom),grtn(3,natom),lpsstr(6)
151 : : real(dp),intent(inout) :: vpsp(nfft)
152 : :
153 : : !Local variables-------------------------------
154 : : !scalars
155 : : character(len=500) :: message
156 : : !arrays
157 : 23956 : real(dp),allocatable :: xcart(:,:)
158 : : #if defined HAVE_BIGDFT
159 : : real(dp),pointer :: rhor_ptr(:,:)
160 : : #endif
161 : :
162 : : ! *************************************************************************
163 : :
164 [ - + ]: 23956 : if (option < 1 .or. option > 4) then
165 : : write(message,'(a,i0,a,a)')&
166 : 0 : & 'From the calling routine, option=',option,ch10,&
167 : 0 : & 'The only allowed values are between 1 and 4.'
168 : 0 : ABI_ERROR(message)
169 : : end if
170 [ - + - - ]: 23956 : if (option > 2 .and. .not.psps%vlspl_recipSpace) then
171 : : write(message,'(a,i0,a,a,a,a)')&
172 : 0 : & 'From the calling routine, option=',option,ch10,&
173 : 0 : & 'but the local part of the pseudo-potential is in real space.',ch10,&
174 : 0 : & 'Action: set icoulomb = 0 to turn-off real space computations.'
175 : 0 : ABI_ERROR(message)
176 : : end if
177 [ - + - - ]: 23956 : if (option > 2 .and. dtset%usewvl == 1) then
178 : : write(message,'(a,i0,a,a)')&
179 : 0 : & 'From the calling routine, option=',option,ch10,&
180 : 0 : & 'but this is not implemented yet from wavelets.'
181 : 0 : ABI_ERROR(message)
182 : : end if
183 : :
184 [ + - ]: 23956 : if (dtset%usewvl == 0) then
185 : : ! Plane wave case
186 [ + - ]: 23956 : if (psps%vlspl_recipSpace) then
187 : : call mklocl_recipspace(dyfrlo,eei,gmet,gprimd,grtn,gsqcut,&
188 : : & dtset%icutcoul,lpsstr,mgfft,mpi_enreg,psps%mqgrid_vl,natom,nattyp, &
189 : : & nfft,ngfft,dtset%nkpt,ntypat,option,ph1d,psps%qgrid_vl,qprtrb,dtset%rcut,&
190 : 23956 : & rhog,rprimd,ucvol,dtset%vcutgeo,psps%vlspl,vprtrb,vpsp)
191 : : else
192 : : call mklocl_realspace(grtn,dtset%icoulomb,mpi_enreg,natom,nattyp,nfft, &
193 : : & ngfft,dtset%nscforder,nspden,ntypat,option,pawtab,psps,rhog,rhor, &
194 : 0 : & rprimd,dtset%typat,ucvol,dtset%usewvl,vpsp,xred)
195 : : end if
196 : : else
197 : : ! Store xcart for each atom
198 [ # # # # : 0 : ABI_MALLOC(xcart,(3, dtset%natom))
# # # # #
# # # ]
199 : 0 : call xred2xcart(dtset%natom, rprimd, xcart, xred)
200 : : ! Eventually retrieve density
201 : : #if defined HAVE_BIGDFT
202 : : if (option>1.and.wvl_den%denspot%rhov_is/=ELECTRONIC_DENSITY) then
203 : : rhor_ptr => rhor ! Just to bypass intent(inout)
204 : : call wvl_rho_abi2big(1,rhor_ptr,wvl_den)
205 : : end if
206 : : #endif
207 : : ! Wavelets case
208 : : call mklocl_wavelets(dtset%efield, grtn, mpi_enreg, dtset%natom, &
209 : : & nfft, nspden, option, rprimd, vpsp, &
210 : 0 : & wvl_den, wvl, xcart)
211 [ # # ]: 0 : ABI_FREE(xcart)
212 : : end if
213 : :
214 : 23956 : end subroutine mklocl
215 : : !!***
216 : :
217 : : !!****f* ABINIT/mklocl_recipspace
218 : : !! NAME
219 : : !! mklocl_recipspace
220 : : !!
221 : : !! FUNCTION
222 : : !! Optionally compute :
223 : : !! option=1 : local ionic potential throughout unit cell
224 : : !! option=2 : contribution of local ionic potential to E gradient wrt xred
225 : : !! option=3 : contribution of local ionic potential to stress tensor
226 : : !! option=4 : contribution of local ionic potential to
227 : : !! second derivative of E wrt xred
228 : : !!
229 : : !! INPUTS
230 : : !! if(option==3) eei=local pseudopotential part of total energy (hartree)
231 : : !! gmet(3,3)=reciprocal space metric ($\textrm{Bohr}^{-2}$).
232 : : !! gprimd(3,3)=reciprocal space dimensional primitive translations
233 : : !! gsqcut=cutoff on $|G|^2$: see setup1 for definition (doubled sphere).
234 : : !! mgfft=maximum size of 1D FFTs
235 : : !! mpi_enreg=information about MPI parallelization
236 : : !! mqgrid=number of grid pts in q array for f(q) spline.
237 : : !! natom=number of atoms in unit cell.
238 : : !! nattyp(ntypat)=number of atoms of each type in cell.
239 : : !! nfft=(effective) number of FFT grid points (for this processor)
240 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
241 : : !! ntypat=number of types of atoms.
242 : : !! option= (see above)
243 : : !! ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phase information.
244 : : !! qgrid(mqgrid)=q grid for spline from 0 to qmax.
245 : : !! qprtrb(3)= integer wavevector of possible perturbing potential
246 : : !! in basis of reciprocal lattice translations
247 : : !! rhog(2,nfft)=electron density rho(G) (electrons/$\textrm{Bohr}^3$)
248 : : !! (needed if option==2 or if option==4)
249 : : !! ucvol=unit cell volume ($\textrm{Bohr}^3$).
250 : : !! vlspl(mqgrid,2,ntypat)=q^2 v(q) spline for each type of atom.
251 : : !! vprtrb(2)=complex amplitude of possible perturbing potential; if nonzero,
252 : : !! perturbing potential is added of the form
253 : : !! $V(G)=(vprtrb(1)+I*vprtrb(2))/2$ at the values G=qprtrb and
254 : : !! $(vprtrb(1)-I*vprtrb(2))/2$ at $G=-qprtrb$ (integers)
255 : : !!
256 : : !! OUTPUT
257 : : !! (if option==1) vpsp(nfft)=local crystal pseudopotential in real space.
258 : : !! (if option==2) grtn(3,natom)=grads of Etot wrt tn.
259 : : !! (if option==3) lpsstr(6)=components of local psp part of stress tensor
260 : : !! (Cartesian coordinates, symmetric tensor) in hartree/$\textrm{bohr}^3$
261 : : !! Store 6 unique components in order 11, 22, 33, 32, 31, 21
262 : : !! (if option==4) dyfrlo(3,3,natom)=d2 Eei/d tn(i)/d tn(j). (Hartrees)
263 : : !!
264 : : !! NOTES
265 : : !! Note that the present routine is tightly connected to the dfpt_vlocal.f routine,
266 : : !! that compute the derivative of the local ionic potential
267 : : !! with respect to one atomic displacement. The argument list
268 : : !! and the internal loops to be considered were sufficiently different
269 : : !! as to make the two routine different.
270 : : !!
271 : : !! SOURCE
272 : :
273 : 29276 : subroutine mklocl_recipspace(dyfrlo,eei,gmet,gprimd,grtn,gsqcut,icutcoul,lpsstr,mgfft,&
274 : 29276 : & mpi_enreg,mqgrid,natom,nattyp,nfft,ngfft,nkpt,ntypat,option,ph1d,qgrid,qprtrb,&
275 : 29276 : & rcut,rhog,rprimd,ucvol,vcutgeo,vlspl,vprtrb,vpsp)
276 : :
277 : : !Arguments ------------------------------------
278 : : !scalars
279 : : integer,intent(in) :: mgfft,mqgrid,natom,nfft,nkpt,ntypat,option,icutcoul
280 : : real(dp),intent(in) :: eei,gsqcut,rcut,rprimd(3,3),ucvol,vcutgeo(3)
281 : : type(MPI_type),intent(in) :: mpi_enreg
282 : : !arrays
283 : : integer,intent(in) :: nattyp(ntypat),ngfft(18),qprtrb(3)
284 : : real(dp),intent(in) :: gmet(3,3),gprimd(3,3),ph1d(2,3*(2*mgfft+1)*natom)
285 : : real(dp),intent(in) :: qgrid(mqgrid),rhog(2,nfft),vlspl(mqgrid,2,ntypat)
286 : : real(dp),intent(in) :: vprtrb(2)
287 : : real(dp),intent(out) :: dyfrlo(3,3,natom),grtn(3,natom),lpsstr(6) !vz_i
288 : : real(dp),intent(inout) :: vpsp(nfft) !vz_i
289 : :
290 : : !Local variables-------------------------------
291 : : !scalars
292 : : integer,parameter :: im=2,re=1
293 : : integer :: i1,i2,i3,ia,ia1,ia2,id1,id2,id3,ierr,ig1,ig2,ig3,ii,itypat
294 : : integer :: jj,me_fft,me_g0,n1,n2,n3,nproc_fft,shift1
295 : : integer :: shift2,shift3
296 : : #ifdef FC_NVHPC
297 : : !Silly trick to prevent NVHPC optimization issue
298 : : logical :: nothing=.false.
299 : : #endif
300 : : real(dp),parameter :: tolfix=1.0000001_dp
301 : : real(dp) :: aa,bb,cc,cutoff,dbl_ig1,dbl_ig2,dbl_ig3,dd,diff,dq,dq2div6,dqdiv6
302 : : real(dp) :: dqm1,ee,ff,gmag,gsquar!beta,gcart_para,gcart_perp
303 : : real(dp) :: ph12i,ph12r,ph1i,ph1r,ph2i,ph2r
304 : : real(dp) :: ph3i,ph3r,phimag_igia,phre_igia,rcut_loc,sfi,sfr
305 : : real(dp) :: svion,svioni,svionr,term,vion1,vion2,xnorm
306 : : character(len=500) :: message
307 : : !arrays
308 : 29276 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
309 : 29276 : integer, contiguous, pointer :: fftn3_distrib(:),ffti3_local(:)
310 : : real(dp) :: gcart(3),tsec(2)
311 : 29276 : real(dp),allocatable :: gcutoff(:)
312 : 29276 : real(dp),allocatable :: work1(:,:)
313 : : ! *************************************************************************
314 : :
315 : : !Define G^2 based on G space metric gmet.
316 : : ! gsq(i1,i2,i3)=dble(i1*i1)*gmet(1,1)+dble(i2*i2)*gmet(2,2)+&
317 : : !& dble(i3*i3)*gmet(3,3)+dble(2*i1*i2)*gmet(1,2)+&
318 : : !& dble(2*i2*i3)*gmet(2,3)+dble(2*i3*i1)*gmet(3,1)
319 : :
320 : : !Real and imaginary parts of phase--statment functions:
321 : : ! phr(x1,y1,x2,y2,x3,y3)=(x1*x2-y1*y2)*x3-(y1*x2+x1*y2)*y3
322 : : ! phi(x1,y1,x2,y2,x3,y3)=(x1*x2-y1*y2)*y3+(y1*x2+x1*y2)*x3
323 : : ! ph1(nri,i1,ia)=ph1d(nri,i1+1+n1+(ia-1)*(2*n1+1))
324 : : ! ph2(nri,i2,ia)=ph1d(nri,i2+1+n2+(ia-1)*(2*n2+1)+&
325 : : !& natom*(2*n1+1))
326 : : ! ph3(nri,i3,ia)=ph1d(nri,i3+1+n3+(ia-1)*(2*n3+1)+&
327 : : !& natom*(2*n1+1+2*n2+1))
328 : : ! phre(i1,i2,i3,ia)=phr(ph1(re,i1,ia),ph1(im,i1,ia),ph2(re,i2,ia),&
329 : : !& ph2(im,i2,ia),ph3(re,i3,ia),ph3(im,i3,ia))
330 : : ! phimag(i1,i2,i3,ia)=phi(ph1(re,i1,ia),ph1(im,i1,ia),ph2(re,i2,ia),&
331 : : !& ph2(im,i2,ia),ph3(re,i3,ia),ph3(im,i3,ia))
332 : :
333 : : !-----
334 : :
335 : : !Keep track of total time spent in mklocl
336 : 29276 : if(option==2)then
337 : 17081 : call timab(72,1,tsec)
338 : : end if
339 : 29276 : n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
340 : 29276 : me_fft=ngfft(11)
341 : 29276 : nproc_fft=ngfft(10)
342 : :
343 : : !Get the distrib associated with this fft_grid
344 : 29276 : call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
345 : :
346 : : !Zero out array to permit accumulation over atom types below:
347 [ + + ]: 29276 : if(option==1)then
348 [ - + + - : 20625 : ABI_MALLOC(work1,(2,nfft))
+ - + - -
+ - + ]
349 [ + + + + ]: 341606267 : work1(:,:)=zero
350 : : end if
351 : :
352 : : !
353 : 29276 : dq=(qgrid(mqgrid)-qgrid(1))/dble(mqgrid-1)
354 : 29276 : dqm1=1.0_dp/dq
355 : 29276 : dqdiv6=dq/6.0_dp
356 : 29276 : dq2div6=dq**2/6.0_dp
357 : 29276 : cutoff=gsqcut*tolfix
358 : 29276 : id1=n1/2+2
359 : 29276 : id2=n2/2+2
360 : 29276 : id3=n3/2+2
361 [ + + + + ]: 349392 : grtn(:,:)=zero
362 : 29276 : lpsstr(:)=zero
363 [ + + + + : 1069653 : dyfrlo(:,:,:)=zero
+ + ]
364 : 29276 : me_g0=0
365 : 29276 : ia1=1
366 : :
367 : : !Initialize Gcut-off array from m_gtermcutoff
368 : 29276 : call termcutoff(gcutoff,gsqcut,icutcoul,ngfft,nkpt,rcut,rprimd,vcutgeo)
369 : 29276 : rcut_loc = half*SQRT(DOT_PRODUCT(rprimd(:,3),rprimd(:,3)))
370 : :
371 [ + + ]: 72612 : do itypat=1,ntypat
372 : : ! ia1,ia2 sets range of loop over atoms:
373 : 43336 : ia2=ia1+nattyp(itypat)-1
374 : :
375 : 43336 : ii=0
376 [ + + ]: 831848 : do i3=1,n3
377 : 788512 : ig3=i3-(i3/id3)*n3-1
378 [ + + ]: 19070317 : do i2=1,n2
379 : 18238469 : ig2=i2-(i2/id2)*n2-1
380 [ + + ]: 19026981 : if(fftn2_distrib(i2) == me_fft ) then
381 [ + + ]: 636685195 : do i1=1,n1
382 : 618851054 : ig1=i1-(i1/id1)*n1-1
383 : :
384 : 618851054 : ii=ii+1
385 : : ! *** GET RID OF THIS THESE IF STATEMENTS (if they slow code)
386 : : ! Skip G=0:
387 : : ! if (ii==1) cycle
388 [ + + ]: 618851054 : if (ig1==0 .and. ig2==0 .and. ig3==0) me_g0=1
389 [ + + ]: 618851054 : if (ig1==0 .and. ig2==0 .and. ig3==0) cycle
390 : :
391 : 618808386 : gsquar=gsq_mk(ig1,ig2,ig3)
392 : : ! Skip G**2 outside cutoff:
393 [ + + ]: 636642527 : if (gsquar<=cutoff) then
394 : 258361020 : gmag=sqrt(gsquar)
395 : :
396 : : ! Compute vion(G) for given type of atom
397 : 258361020 : jj=1+int(gmag*dqm1)
398 : 258361020 : diff=gmag-qgrid(jj)
399 : :
400 : : ! Evaluate spline fit from q^2 V(q) to get V(q):
401 : : ! (p. 86 Numerical Recipes, Press et al;
402 : : ! NOTE error in book for sign
403 : : ! of "aa" term in derivative; also see splfit routine).
404 : :
405 : 258361020 : bb = diff*dqm1
406 : 258361020 : aa = 1.0_dp-bb
407 : 258361020 : cc = aa*(aa**2-1.0_dp)*dq2div6
408 : 258361020 : dd = bb*(bb**2-1.0_dp)*dq2div6
409 : :
410 : : vion1 = (aa*vlspl(jj,1,itypat)+bb*vlspl(jj+1,1,itypat) +&
411 : 258361020 : & cc*vlspl(jj,2,itypat)+dd*vlspl(jj+1,2,itypat) ) / gsquar * gcutoff(ii)
412 : :
413 [ + + ]: 258361020 : if(option==1)then
414 : :
415 : : ! Assemble structure factor over all atoms of given type:
416 : 58696177 : sfr=zero
417 : 58696177 : sfi=zero
418 [ + + ]: 193706198 : do ia=ia1,ia2
419 : 135010021 : sfr=sfr+phre_mk(ig1,ig2,ig3,ia)
420 : 193706198 : sfi=sfi-phimag_mk(ig1,ig2,ig3,ia)
421 : : end do
422 : : ! Multiply structure factor times vion:
423 : 58696177 : work1(re,ii)=work1(re,ii)+sfr*vion1
424 : 58696177 : work1(im,ii)=work1(im,ii)+sfi*vion1
425 : :
426 [ + + ]: 199664843 : else if(option==2 .or. option==4)then
427 : :
428 : : ! Compute Re and Im part of (2Pi)*Vion(G)*rho(G):
429 : 152111919 : svionr=(two_pi*vion1)*rhog(re,ii)
430 : 152111919 : svioni=(two_pi*vion1)*rhog(im,ii)
431 : :
432 : : ! Loop over all atoms of this type:
433 [ + + ]: 433890751 : do ia=ia1,ia2
434 : : #ifdef FC_NVHPC
435 : : !Silly trick to prevent NVHPC optimization issue
436 : : if(nothing) write(100,*) shift1,shift2,shift3
437 : : #endif
438 : :
439 : 281778832 : shift1=1+n1+(ia-1)*(2*n1+1)
440 : 281778832 : shift2=1+n2+(ia-1)*(2*n2+1)+natom*(2*n1+1)
441 : 281778832 : shift3=1+n3+(ia-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1)
442 : 281778832 : ph1r=ph1d(1,ig1+shift1)
443 : 281778832 : ph1i=ph1d(2,ig1+shift1)
444 : 281778832 : ph2r=ph1d(1,ig2+shift2)
445 : 281778832 : ph2i=ph1d(2,ig2+shift2)
446 : 281778832 : ph3r=ph1d(1,ig3+shift3)
447 : 281778832 : ph3i=ph1d(2,ig3+shift3)
448 : 281778832 : ph12r=ph1r*ph2r-ph1i*ph2i
449 : 281778832 : ph12i=ph1r*ph2i+ph1i*ph2r
450 : 281778832 : phre_igia=ph12r*ph3r-ph12i*ph3i
451 : 281778832 : phimag_igia=ph12r*ph3i+ph12i*ph3r
452 : :
453 [ + + ]: 433890751 : if(option==2)then
454 : :
455 : : ! Compute "Vion" part of gradient
456 : : ! svion=svioni*phre(ig1,ig2,ig3,ia)+svionr*phimag(ig1,ig2,ig3,ia)
457 : 277962554 : svion=svioni*phre_igia+svionr*phimag_igia
458 : :
459 : : ! Open loop over 3-index for speed:
460 : 277962554 : grtn(1,ia)=grtn(1,ia)-dble(ig1)*svion
461 : 277962554 : grtn(2,ia)=grtn(2,ia)-dble(ig2)*svion
462 : 277962554 : grtn(3,ia)=grtn(3,ia)-dble(ig3)*svion
463 : :
464 : : else
465 : :
466 : : ! Compute "Vion" part of the second derivative
467 : : ! svion=two_pi*
468 : : ! (svionr*phre(ig1,ig2,ig3,ia)-svioni*phimag(ig1,ig2,ig3,ia))
469 : 3816278 : svion=two_pi*(svionr*phre_igia-svioni*phimag_igia)
470 : :
471 : : ! Open loop over 3-index for speed
472 : 3816278 : dbl_ig1=dble(ig1) ; dbl_ig2=dble(ig2) ; dbl_ig3=dble(ig3)
473 : 3816278 : dyfrlo(1,1,ia)=dyfrlo(1,1,ia)-dbl_ig1*dbl_ig1*svion
474 : 3816278 : dyfrlo(1,2,ia)=dyfrlo(1,2,ia)-dbl_ig1*dbl_ig2*svion
475 : 3816278 : dyfrlo(1,3,ia)=dyfrlo(1,3,ia)-dbl_ig1*dbl_ig3*svion
476 : 3816278 : dyfrlo(2,2,ia)=dyfrlo(2,2,ia)-dbl_ig2*dbl_ig2*svion
477 : 3816278 : dyfrlo(2,3,ia)=dyfrlo(2,3,ia)-dbl_ig2*dbl_ig3*svion
478 : 3816278 : dyfrlo(3,3,ia)=dyfrlo(3,3,ia)-dbl_ig3*dbl_ig3*svion
479 : :
480 : : end if
481 : :
482 : : end do
483 : :
484 [ + - ]: 47552924 : else if(option==3)then
485 : : ! if(icutcoul .ne. 2) then
486 : : ! Also get (dV(q)/dq)/q:
487 : : ! (note correction of Numerical Recipes sign error
488 : : ! before (3._dp*aa**2-1._dp)
489 : : ! ee*dqm1 + ff*dqdiv6 is the best estimate of dV(q)/dq from splines
490 : 47552924 : ee= vlspl(jj+1,1,itypat)-vlspl(jj,1,itypat)
491 : : ff= (3._dp*bb**2-1._dp)*vlspl(jj+1,2,itypat) &
492 : 47552924 : & - (3._dp*aa**2-1._dp)*vlspl(jj,2,itypat)
493 : : vion2 = ( ( ee*dqm1 + ff*dqdiv6 )/gmag&
494 : 47552924 : & - 2.0_dp*vion1 ) / gsquar
495 : :
496 : : gcart(1)=gprimd(1,1)*dble(ig1)+gprimd(1,2)*dble(ig2)+&
497 : 47552924 : & gprimd(1,3)*dble(ig3)
498 : : gcart(2)=gprimd(2,1)*dble(ig1)+gprimd(2,2)*dble(ig2)+&
499 : 47552924 : & gprimd(2,3)*dble(ig3)
500 : : gcart(3)=gprimd(3,1)*dble(ig1)+gprimd(3,2)*dble(ig2)+&
501 : 47552924 : & gprimd(3,3)*dble(ig3)
502 : : ! Assemble structure over all atoms of given type
503 : 47552924 : sfr=zero
504 : 47552924 : sfi=zero
505 [ + + ]: 156538141 : do ia=ia1,ia2
506 : 108985217 : sfr=sfr+phre_mk(ig1,ig2,ig3,ia)
507 : 156538141 : sfi=sfi-phimag_mk(ig1,ig2,ig3,ia)
508 : : end do
509 : : ! Compute Re( rho^*(G)* sf ) * [(dV(G)/dG)/|G|]
510 : 47552924 : term=(rhog(re,ii)*sfr+rhog(im,ii)*sfi)*vion2
511 : : ! Compute contribution to stress tensor
512 : 47552924 : lpsstr(1)=lpsstr(1)-term*gcart(1)*gcart(1)
513 : 47552924 : lpsstr(2)=lpsstr(2)-term*gcart(2)*gcart(2)
514 : 47552924 : lpsstr(3)=lpsstr(3)-term*gcart(3)*gcart(3)
515 : 47552924 : lpsstr(4)=lpsstr(4)-term*gcart(3)*gcart(2)
516 : 47552924 : lpsstr(5)=lpsstr(5)-term*gcart(3)*gcart(1)
517 : 47552924 : lpsstr(6)=lpsstr(6)-term*gcart(2)*gcart(1)
518 : : ! else if (icutcoul .eq. 2) then
519 : : !! Also get (dV(q)/dq)/q:
520 : : !! (note correction of Numerical Recipes sign error
521 : : !! before (3._dp*aa**2-1._dp)
522 : : !! ee*dqm1 + ff*dqdiv6 is the best estimate of dV(q)/dq from splines
523 : : ! ee= vlspl(jj+1,1,itypat)-vlspl(jj,1,itypat)
524 : : ! ff= (3._dp*bb**2-1._dp)*vlspl(jj+1,2,itypat) &
525 : : !& - (3._dp*aa**2-1._dp)*vlspl(jj,2,itypat)
526 : : ! vion2 = ( ( ee*dqm1 + ff*dqdiv6 )/gmag&
527 : : !& - 2.0_dp*vion1 ) / gsquar
528 : : !
529 : : ! gcart(1)=gprimd(1,1)*dble(ig1)+gprimd(1,2)*dble(ig2)+&
530 : : !& gprimd(1,3)*dble(ig3)
531 : : ! gcart(2)=gprimd(2,1)*dble(ig1)+gprimd(2,2)*dble(ig2)+&
532 : : !& gprimd(2,3)*dble(ig3)
533 : : ! gcart(3)=gprimd(3,1)*dble(ig1)+gprimd(3,2)*dble(ig2)+&
534 : : !& gprimd(3,3)*dble(ig3)
535 : : !! Assemble structure over all atoms of given type
536 : : ! sfr=zero
537 : : ! sfi=zero
538 : : ! do ia=ia1,ia2
539 : : ! sfr=sfr+phre_mk(ig1,ig2,ig3,ia)
540 : : ! sfi=sfi-phimag_mk(ig1,ig2,ig3,ia)
541 : : ! end do
542 : : ! !Implement beta correction as in eq. 62 (PRB 96 075448 2017)
543 : : ! gcart_para = sqrt(gcart(1)**2+gcart(2)**2)
544 : : ! gcart_perp = gcart(3)
545 : : ! gsquar = gcart(1)**2+gcart(2)**2+gcart(3)**2
546 : : ! if(gcart_para .gt. tol12) then
547 : : ! beta = gsquar*rcut_loc/(two*gcart_para)* &
548 : : ! & exp(-gcart_para*rcut_loc)* &
549 : : ! &cos(gcart_perp*rcut_loc)/(one-exp(-gcart_para*rcut_loc)*cos(gcart_perp*rcut_loc))
550 : : ! else
551 : : ! beta = zero
552 : : ! end if
553 : : !! Compute Re( rho^*(G)* sf ) * [(dV(G)/dG)/|G|]
554 : : ! term=(rhog(re,ii)*sfr+rhog(im,ii)*sfi)*vion2
555 : : !! Compute contribution to stress tensor
556 : : ! lpsstr(1)=lpsstr(1)-term*(gcart(1)*gcart(1))*(1+beta)
557 : : ! lpsstr(2)=lpsstr(2)-term*(gcart(2)*gcart(2))*(1+beta)
558 : : ! lpsstr(3)=lpsstr(3)-term*(gcart(3)*gcart(3)-gsquar)
559 : : ! lpsstr(4)=lpsstr(4)-term*gcart(3)*gcart(2)
560 : : ! lpsstr(5)=lpsstr(5)-term*gcart(3)*gcart(1)
561 : : ! lpsstr(6)=lpsstr(6)-term*gcart(2)*gcart(1)
562 : : ! endif
563 : :
564 : : else
565 : 0 : write(message, '(a,i0,a)' )' mklocl: Option=',option,' not allowed.'
566 : 0 : ABI_BUG(message)
567 : : end if ! End option choice
568 : :
569 : : ! End skip G**2 outside cutoff:
570 : : end if
571 : :
572 : : ! End loop on n1, n2, n3. There is a "cycle" inside the loop
573 : : end do
574 : : end if ! this plane is for me_fft
575 : : end do
576 : : end do
577 : :
578 : : ! Symmetrize the dynamical matrix with respect to indices
579 [ + + ]: 123365 : do ia=ia1,ia2
580 : 80029 : dyfrlo(2,1,ia)=dyfrlo(1,2,ia)
581 : 80029 : dyfrlo(3,1,ia)=dyfrlo(1,3,ia)
582 : 123365 : dyfrlo(3,2,ia)=dyfrlo(2,3,ia)
583 : : end do
584 : :
585 : 72612 : ia1=ia2+1
586 : :
587 : : ! End loop on type of atoms
588 : : end do
589 : :
590 [ + + ]: 29276 : if(option==1)then
591 : : ! Dont't change work1 on g=0 if Poisson solver is used since work1
592 : : ! hold not the potential but the density generated by the pseudo.
593 [ + + ]: 6875 : if(me_g0 == 1) then
594 : : ! Set Vloc(G=0)=0:
595 : 6771 : work1(re,1)=zero
596 : 6771 : work1(im,1)=zero
597 : : end if
598 : : ! write(std_out,*) ' mklocl_recipspace : will add potential with strength vprtrb(:)=',vprtrb(:)
599 : :
600 : : ! Allow for the addition of a perturbing potential
601 [ + + ]: 6875 : if ((vprtrb(1)**2+vprtrb(2)**2) > 1.d-30) then
602 : : ! Find the linear indices which correspond with the input
603 : : ! wavevector qprtrb
604 : : ! The double modulus handles both i>=n and i<0, mapping into [0,n-1];
605 : : ! then add 1 to get range [1,n] for each
606 : 8 : i3=1+mod(n3+mod(qprtrb(3),n3),n3)
607 : 8 : i2=1+mod(n2+mod(qprtrb(2),n2),n2)
608 : 8 : i1=1+mod(n1+mod(qprtrb(1),n1),n1)
609 : : ! Compute the linear index in the 3 dimensional array
610 : 8 : ii=i1+n1*((ffti2_local(i2)-1)+(n2/nproc_fft)*(i3-1))
611 : : ! Add in the perturbation at G=qprtrb
612 : 8 : work1(re,ii)=work1(re,ii)+0.5_dp*vprtrb(1)
613 : 8 : work1(im,ii)=work1(im,ii)+0.5_dp*vprtrb(2)
614 : : ! Same thing for G=-qprtrb
615 : 8 : i3=1+mod(n3+mod(-qprtrb(3),n3),n3)
616 : 8 : i2=1+mod(n2+mod(-qprtrb(2),n2),n2)
617 : 8 : i1=1+mod(n1+mod(-qprtrb(1),n1),n1)
618 : : ! ii=i1+n1*((i2-1)+n2*(i3-1))
619 : 8 : work1(re,ii)=work1(re,ii)+0.5_dp*vprtrb(1)
620 : 8 : work1(im,ii)=work1(im,ii)-0.5_dp*vprtrb(2)
621 : : write(message, '(a,1p,2e12.4,a,0p,3i4,a)' )&
622 : 8 : & ' mklocl: perturbation of vprtrb=', vprtrb,&
623 : 16 : & ' and q=',qprtrb,' has been added'
624 : 8 : call wrtout(std_out,message,'COLL')
625 : : end if
626 : :
627 : : ! Transform back to real space
628 : 6875 : call fourdp(1,work1,vpsp,1,mpi_enreg,nfft,1,ngfft,0)
629 : :
630 : : ! Divide by unit cell volume
631 : 6875 : xnorm=1.0_dp/ucvol
632 [ + + ]: 113873339 : vpsp(:)=vpsp(:)*xnorm
633 : :
634 [ - + ]: 6875 : ABI_FREE(work1)
635 : :
636 : : end if
637 : :
638 [ - + ]: 29276 : ABI_FREE(gcutoff)
639 : :
640 [ + + ]: 29276 : if(option==2)then
641 : : ! Init mpi_comm
642 [ + + ]: 17081 : if(mpi_enreg%nproc_fft>1)then
643 : 904 : call timab(48,1,tsec)
644 : 904 : call xmpi_sum(grtn,mpi_enreg%comm_fft ,ierr)
645 : 904 : call timab(48,2,tsec)
646 : : end if
647 : 17081 : call timab(72,2,tsec)
648 : : end if
649 : :
650 [ + + ]: 29276 : if(option==3)then
651 : : ! Init mpi_comm
652 [ + + ]: 4915 : if(mpi_enreg%nproc_fft>1)then
653 : 200 : call timab(48,1,tsec)
654 : 200 : call xmpi_sum(lpsstr,mpi_enreg%comm_fft ,ierr)
655 : 200 : call timab(48,2,tsec)
656 : : end if
657 : :
658 : : ! Normalize and add term -eei/ucvol on diagonal
659 : : ! (see page 802 of notes)
660 : : ! if(icutcoul .ne. 2) then
661 : 4915 : lpsstr(1)=(lpsstr(1)-eei)/ucvol
662 : 4915 : lpsstr(2)=(lpsstr(2)-eei)/ucvol
663 : 4915 : lpsstr(3)=(lpsstr(3)-eei)/ucvol
664 : 4915 : lpsstr(4)=lpsstr(4)/ucvol
665 : 4915 : lpsstr(5)=lpsstr(5)/ucvol
666 : 4915 : lpsstr(6)=lpsstr(6)/ucvol
667 : : ! elseif (icutcoul .eq. 2) then
668 : : ! lpsstr(1)=(lpsstr(1)-eei)/ucvol
669 : : ! lpsstr(2)=(lpsstr(2)-eei)/ucvol
670 : : ! lpsstr(3)=(lpsstr(3)-eei)/ucvol
671 : : ! lpsstr(4)=lpsstr(4)/ucvol
672 : : ! lpsstr(5)=lpsstr(5)/ucvol
673 : : ! lpsstr(6)=lpsstr(6)/ucvol
674 : : !lpsstr=lpsstr/ucvol
675 : : ! endif
676 : :
677 : : end if
678 : :
679 [ + + + + : 58957 : if(option==4)then
- + ]
680 : : ! Init mpi_comm
681 [ - + ]: 405 : if(mpi_enreg%nproc_fft>1)then
682 : 0 : call timab(48,1,tsec)
683 : 0 : call xmpi_sum(dyfrlo,mpi_enreg%comm_fft ,ierr)
684 : 0 : call timab(48,2,tsec)
685 : : end if
686 : : end if
687 : :
688 : : contains
689 : :
690 : : !Real and imaginary parts of phase--statment functions:
691 : 243995238 : function phr_mk(x1,y1,x2,y2,x3,y3)
692 : :
693 : : real(dp) :: phr_mk,x1,x2,x3,y1,y2,y3
694 : 243995238 : phr_mk=(x1*x2-y1*y2)*x3-(y1*x2+x1*y2)*y3
695 : : end function phr_mk
696 : :
697 : 243995238 : function phi_mk(x1,y1,x2,y2,x3,y3)
698 : :
699 : : real(dp):: phi_mk,x1,x2,x3,y1,y2,y3
700 : 243995238 : phi_mk=(x1*x2-y1*y2)*y3+(y1*x2+x1*y2)*x3
701 : : end function phi_mk
702 : :
703 : 975980952 : function ph1_mk(nri,ig1,ia)
704 : :
705 : : real(dp):: ph1_mk
706 : : integer :: nri,ig1,ia
707 : 975980952 : ph1_mk=ph1d(nri,ig1+1+n1+(ia-1)*(2*n1+1))
708 : 975980952 : end function ph1_mk
709 : :
710 : 975980952 : function ph2_mk(nri,ig2,ia)
711 : :
712 : : real(dp):: ph2_mk
713 : : integer :: nri,ig2,ia
714 : 975980952 : ph2_mk=ph1d(nri,ig2+1+n2+(ia-1)*(2*n2+1)+natom*(2*n1+1))
715 : 975980952 : end function ph2_mk
716 : :
717 : 975980952 : function ph3_mk(nri,ig3,ia)
718 : :
719 : : real(dp):: ph3_mk
720 : : integer :: nri,ig3,ia
721 : 975980952 : ph3_mk=ph1d(nri,ig3+1+n3+(ia-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1))
722 : 975980952 : end function ph3_mk
723 : :
724 : 243995238 : function phre_mk(ig1,ig2,ig3,ia)
725 : :
726 : : real(dp):: phre_mk
727 : : integer :: ig1,ig2,ig3,ia
728 : : phre_mk=phr_mk(ph1_mk(re,ig1,ia),ph1_mk(im,ig1,ia),&
729 : 243995238 : & ph2_mk(re,ig2,ia),ph2_mk(im,ig2,ia),ph3_mk(re,ig3,ia),ph3_mk(im,ig3,ia))
730 : 243995238 : end function phre_mk
731 : :
732 : 243995238 : function phimag_mk(ig1,ig2,ig3,ia)
733 : :
734 : : real(dp) :: phimag_mk
735 : : integer :: ig1,ig2,ig3,ia
736 : : phimag_mk=phi_mk(ph1_mk(re,ig1,ia),ph1_mk(im,ig1,ia),&
737 : 243995238 : & ph2_mk(re,ig2,ia),ph2_mk(im,ig2,ia),ph3_mk(re,ig3,ia),ph3_mk(im,ig3,ia))
738 : 243995238 : end function phimag_mk
739 : :
740 : 618808386 : function gsq_mk(i1,i2,i3)
741 : :
742 : : real(dp) :: gsq_mk
743 : : integer :: i1,i2,i3
744 : : gsq_mk=dble(i1*i1)*gmet(1,1)+dble(i2*i2)*gmet(2,2)+&
745 : : & dble(i3*i3)*gmet(3,3)+dble(2*i1*i2)*gmet(1,2)+&
746 : 618808386 : & dble(2*i2*i3)*gmet(2,3)+dble(2*i3*i1)*gmet(3,1)
747 : 618808386 : end function gsq_mk
748 : :
749 : : end subroutine mklocl_recipspace
750 : : !!***
751 : :
752 : : !!****f* ABINIT/dfpt_vlocal
753 : : !! NAME
754 : : !! dfpt_vlocal
755 : : !!
756 : : !! FUNCTION
757 : : !! Compute local part of 1st-order potential from the appropriate
758 : : !! atomic pseudopotential with structure and derivative factor.
759 : : !! In case of derivative with respect to k or
760 : : !! electric (magnetic Zeeman) field perturbation, the 1st-order local potential vanishes.
761 : : !!
762 : : !! INPUTS
763 : : !! atindx(natom)=index table for atoms (see gstate.f)
764 : : !! cplex: if 1, real space 1-order functions on FFT grid
765 : : !! are REAL, if 2, COMPLEX
766 : : !! gmet(3,3)=reciprocal space metric (Bohr**-2)
767 : : !! gsqcut=cutoff G**2 for included G s in fft box.
768 : : !! icutcoul= type of Coulomb cutoff to apply
769 : : !! idir=direction of atomic displacement (=1,2 or 3 : displacement of
770 : : !! atom ipert along the 1st, 2nd or 3rd axis).
771 : : !! ipert=number of the atom being displaced in the frozen-phonon
772 : : !! mpi_enreg=information about MPI parallelization
773 : : !! mqgrid=dimension of q grid for pseudopotentials
774 : : !! natom=number of atoms in cell.
775 : : !! nattyp(ntypat)=number of atoms of each type in cell.
776 : : !! nfft=(effective) number of FFT grid points (for this processor)
777 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
778 : : !! ntypat=number of types of atoms in cell.
779 : : !! n1,n2,n3=fft grid.
780 : : !! ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phase information.
781 : : !! qgrid(mqgrid)=grid of q points from 0 to qmax.
782 : : !! qphon(3)=wavevector of the phonon
783 : : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
784 : : !! ucvol=unit cell volume (Bohr**3).
785 : : !! vcutgeo(3)= array to describe the geometry of the Coulomb cutoff
786 : : !! vlspl(mqgrid,2,ntypat)=spline fit of q^2 V(q) for each type of atom.
787 : : !! xred(3,natom)=reduced atomic coordinates
788 : : !!
789 : : !! OUTPUT
790 : : !! vpsp1(cplex*nfft)=first-order local crystal pseudopotential in real space
791 : : !! (including the minus sign, forgotten in the paper non-linear..
792 : : !!
793 : : !! SOURCE
794 : :
795 : 11870 : subroutine dfpt_vlocal(atindx,cplex,gmet,gsqcut,icutcoul,idir,ipert,&
796 : 11870 : & mpi_enreg,mqgrid,natom,nattyp,nfft,ngfft,nkpt,&
797 : 11870 : & ntypat,n1,n2,n3,ph1d,qgrid,qphon,rcut,rprimd,ucvol,vcutgeo,vlspl,vpsp1,xred,&
798 : : & zion) !Optional
799 : :
800 : : !Arguments -------------------------------
801 : : !scalars
802 : : integer,intent(in) :: cplex,icutcoul,idir,ipert,mqgrid,n1,n2,n3,natom,nfft,nkpt,ntypat
803 : : real(dp),intent(in) :: gsqcut,rcut,ucvol
804 : : type(MPI_type),intent(in) :: mpi_enreg
805 : : !arrays
806 : : integer,intent(in) :: atindx(natom),nattyp(ntypat),ngfft(18)
807 : : real(dp),intent(in) :: gmet(3,3),ph1d(2,(2*n1+1+2*n2+1+2*n3+1)*natom)
808 : : real(dp),intent(in) :: qgrid(mqgrid),qphon(3),rprimd(3,3),vcutgeo(3),vlspl(mqgrid,2,ntypat)
809 : : real(dp),intent(in) :: xred(3,natom)
810 : : real(dp),intent(in),optional :: zion(ntypat)
811 : : real(dp),intent(out) :: vpsp1(cplex*nfft)
812 : :
813 : : !Local variables -------------------------
814 : : !scalars
815 : : integer :: i1,i2,i3,ia1,iatom,id1,id2,id3,ig1,ig2,ig3,ii,ii1,im=2
816 : : integer :: itypat,jj,re=1
817 : : real(dp),parameter :: tolfix=1.000000001_dp
818 : : real(dp) :: aa,bb,cc,cutoff,dd,diff,dq,dq2div6,dqdiv6,dqm1,facg0,gmag,gq1
819 : : real(dp) :: gq2,gq3,gsquar,phqim,phqre
820 : : real(dp) :: qxred2pi,sfi,sfr,vion1,xnorm
821 : : logical :: qeq0
822 : : !arrays
823 : 11870 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
824 : 11870 : integer, contiguous, pointer :: fftn3_distrib(:),ffti3_local(:)
825 : : real(dp) :: gq(3)
826 : 11870 : real(dp),allocatable :: gcutoff(:)
827 : 11870 : real(dp),allocatable :: work1(:,:)
828 : : ! *********************************************************************
829 : :
830 : 11870 : iatom=ipert
831 : :
832 : : if(iatom==natom+1 .or. iatom==natom+2 .or. iatom==natom+10 .or. iatom==natom+11 &
833 [ + + + + : 23740 : & .or. iatom==natom+5 .or. iatom==natom+6 .or. iatom==natom+6 .or. (iatom>natom+11.and.iatom<=2*natom+11))then
+ + + + +
+ + + + +
+ - - + ]
834 : :
835 : : ! (In case of d/dk or an electric field, or magnetic (Zeeman) field->[natom+5] SPr deb )
836 [ + + ]: 4537224 : vpsp1(1:cplex*nfft)=zero
837 : :
838 : : else
839 : :
840 : : ! (In case of a phonon perturbation)
841 [ - + + - : 33681 : ABI_MALLOC(work1,(2,nfft))
+ - + - -
+ - + ]
842 [ + + + + ]: 350520388 : work1(1:2,1:nfft)=0.0_dp
843 : :
844 : 11227 : dq=(qgrid(mqgrid)-qgrid(1))/dble(mqgrid-1)
845 : 11227 : dqm1=1.0_dp/dq
846 : 11227 : dqdiv6=dq/6.0_dp
847 : 11227 : dq2div6=dq**2/6.0_dp
848 : 11227 : cutoff=gsqcut*tolfix
849 : 11227 : id1=n1/2+2
850 : 11227 : id2=n2/2+2
851 : 11227 : id3=n3/2+2
852 : :
853 : : ! Get the distrib associated with this fft_grid
854 : 11227 : call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
855 : :
856 : : ! This is to allow q=0
857 : 11227 : qeq0=.false.
858 [ + + ]: 11227 : if(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15)qeq0=.true.
859 : :
860 : : ! Determination of the atom type
861 : 11227 : ia1=0
862 : 11227 : itypat=0
863 [ + + ]: 27447 : do ii=1,ntypat
864 : 16220 : ia1=ia1+nattyp(ii)
865 [ + + + + ]: 32440 : if(atindx(iatom)<=ia1.and.itypat==0)itypat=ii
866 : : end do
867 : :
868 : : ! Determination of phase qxred*
869 : : qxred2pi=2.0_dp*pi*(qphon(1)*xred(1,iatom)+ &
870 : : & qphon(2)*xred(2,iatom)+ &
871 : 11227 : & qphon(3)*xred(3,iatom) )
872 : 11227 : phqre=cos(qxred2pi)
873 : 11227 : phqim=sin(qxred2pi)
874 : 11227 : ii=0
875 : :
876 : : ! Initialize Gcut-off array from m_gtermcutoff
877 : 11227 : call termcutoff(gcutoff,gsqcut,icutcoul,ngfft,nkpt,rcut,rprimd,vcutgeo,qpt=qphon)
878 : :
879 [ + + ]: 228538 : do i3=1,n3
880 : 217311 : ig3=i3-(i3/id3)*n3-1
881 : 217311 : gq3=dble(ig3)+qphon(3)
882 : 217311 : gq(3)=gq3
883 [ + + ]: 4639861 : do i2=1,n2
884 [ + - ]: 4628634 : if (fftn2_distrib(i2)==mpi_enreg%me_fft) then
885 : 4411323 : ig2=i2-(i2/id2)*n2-1
886 : 4411323 : gq2=dble(ig2)+qphon(2)
887 : 4411323 : gq(2)=gq2
888 : :
889 : : ! Note the lower limit of the next loop
890 : 4411323 : ii1=1
891 : 4411323 : facg0=zero
892 [ + + + + : 4411323 : if(i3==1 .and. i2==1 .and. qeq0 .and. ig2==0 .and. ig3==0)then
+ - ]
893 : 6696 : ii1=2
894 : 6696 : ii=ii+1
895 : : end if
896 [ + + + + : 4411323 : if(i3==1 .and. i2==1 .and. (.not.qeq0) .and. icutcoul==55 .and. ig2==0 .and. ig3==0)then
+ + + - ]
897 : 36 : facg0 = four_pi * zion(itypat) / (two_pi)**2
898 : : end if
899 [ + + ]: 121241014 : do i1=ii1,n1
900 : 116829691 : ig1=i1-(i1/id1)*n1-1
901 : 116829691 : gq1=dble(ig1)+qphon(1)
902 : 116829691 : gq(1)=gq1
903 : 116829691 : ii=ii+1
904 : 116829691 : gsquar=gsq_vl3(gq1,gq2,gq3)
905 : : ! Skip G**2 outside cutoff:
906 [ + + ]: 121241014 : if (gsquar<=cutoff) then
907 : 41986670 : gmag=sqrt(gsquar)
908 : :
909 : : ! Compute vion(G) for given type of atom
910 : 41986670 : jj=1+int(gmag*dqm1)
911 : 41986670 : diff=gmag-qgrid(jj)
912 : :
913 : : ! Evaluate spline fit from q^2 V(q) to get V(q):
914 : : ! (p. 86 Numerical Recipes, Press et al; NOTE error in book for sign
915 : : ! of "aa" term in derivative; also see splfit routine.
916 : : ! This bug fixed here 27 Jan 1992.)
917 : :
918 : 41986670 : bb = diff*dqm1
919 : 41986670 : aa = 1.0_dp-bb
920 : 41986670 : cc = aa*(aa**2-1.0_dp)*dq2div6
921 : 41986670 : dd = bb*(bb**2-1.0_dp)*dq2div6
922 : : vion1 = (aa*vlspl(jj,1,itypat)+bb*vlspl(jj+1,1,itypat) + &
923 : : & cc*vlspl(jj,2,itypat)+dd*vlspl(jj+1,2,itypat) + facg0 ) &
924 : 41986670 : & / gsquar*gcutoff(ii)
925 : :
926 : : ! Phase G*xred (complex conjugate) * -i *2pi*(g+q)*vion
927 : 41986670 : sfr=-phimag_vl3(ig1,ig2,ig3,iatom)*2.0_dp*pi*gq(idir)*vion1
928 : 41986670 : sfi=-phre_vl3(ig1,ig2,ig3,iatom)*2.0_dp*pi*gq(idir)*vion1
929 : :
930 : : ! Phase q*xred (complex conjugate)
931 : 41986670 : work1(re,ii)=sfr*phqre+sfi*phqim
932 : 41986670 : work1(im,ii)=-sfr*phqim+sfi*phqre
933 : :
934 : 41986670 : facg0=zero
935 : : end if
936 : :
937 : : end do
938 : : end if
939 : : end do
940 : : end do
941 : :
942 : : ! Transform back to real space
943 : 11227 : call fourdp(cplex,work1,vpsp1,1,mpi_enreg,nfft,1,ngfft,0)
944 : :
945 : 11227 : xnorm=1.0_dp/ucvol
946 [ + + ]: 150557733 : vpsp1(1:cplex*nfft)=vpsp1(1:cplex*nfft)*xnorm
947 : :
948 [ - + ]: 11227 : ABI_FREE(work1)
949 [ - + ]: 11227 : ABI_FREE(gcutoff)
950 : :
951 : : ! End the condition of non-electric-field
952 : : end if
953 : :
954 : : contains
955 : :
956 : : !Real and imaginary parts of phase.
957 : 41986670 : function phr_vl3(x1,y1,x2,y2,x3,y3)
958 : :
959 : : real(dp) :: phr_vl3
960 : : real(dp),intent(in) :: x1,x2,x3,y1,y2,y3
961 : 41986670 : phr_vl3=(x1*x2-y1*y2)*x3-(y1*x2+x1*y2)*y3
962 : : end function phr_vl3
963 : :
964 : 41986670 : function phi_vl3(x1,y1,x2,y2,x3,y3)
965 : :
966 : : real(dp) :: phi_vl3
967 : : real(dp),intent(in) :: x1,x2,x3,y1,y2,y3
968 : 41986670 : phi_vl3=(x1*x2-y1*y2)*y3+(y1*x2+x1*y2)*x3
969 : : end function phi_vl3
970 : :
971 : : ! Warning : this function differ from similar ones for ground-state calculations : note the atindx !!
972 : 167946680 : function ph1_vl3(nri,ig1,ia)
973 : :
974 : : real(dp) :: ph1_vl3
975 : : integer,intent(in) :: nri,ig1,ia
976 : 167946680 : ph1_vl3=ph1d(nri,ig1+1+n1+(atindx(ia)-1)*(2*n1+1))
977 : 167946680 : end function ph1_vl3
978 : :
979 : : ! Warning : this function differ from similar ones for ground-state calculations : note the atindx !!
980 : 167946680 : function ph2_vl3(nri,ig2,ia)
981 : :
982 : : real(dp) :: ph2_vl3
983 : : integer,intent(in) :: nri,ig2,ia
984 : 167946680 : ph2_vl3=ph1d(nri,ig2+1+n2+(atindx(ia)-1)*(2*n2+1)+natom*(2*n1+1))
985 : 167946680 : end function ph2_vl3
986 : :
987 : : ! Warning : this function differ from similar ones for ground-state calculations : note the atindx !!
988 : 167946680 : function ph3_vl3(nri,ig3,ia)
989 : :
990 : : real(dp) :: ph3_vl3
991 : : integer,intent(in) :: nri,ig3,ia
992 : 167946680 : ph3_vl3=ph1d(nri,ig3+1+n3+(atindx(ia)-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1))
993 : 167946680 : end function ph3_vl3
994 : :
995 : 41986670 : function phre_vl3(ig1,ig2,ig3,ia)
996 : :
997 : : real(dp) :: phre_vl3
998 : : integer,intent(in) :: ig1,ig2,ig3,ia
999 : : phre_vl3=phr_vl3(ph1_vl3(re,ig1,ia),ph1_vl3(im,ig1,ia),&
1000 : 41986670 : & ph2_vl3(re,ig2,ia),ph2_vl3(im,ig2,ia),ph3_vl3(re,ig3,ia),ph3_vl3(im,ig3,ia))
1001 : 41986670 : end function phre_vl3
1002 : :
1003 : 41986670 : function phimag_vl3(ig1,ig2,ig3,ia)
1004 : :
1005 : : real(dp) :: phimag_vl3
1006 : : integer,intent(in) :: ig1,ig2,ig3,ia
1007 : : phimag_vl3=phi_vl3(ph1_vl3(re,ig1,ia),ph1_vl3(im,ig1,ia),&
1008 : 41986670 : & ph2_vl3(re,ig2,ia),ph2_vl3(im,ig2,ia),ph3_vl3(re,ig3,ia),ph3_vl3(im,ig3,ia))
1009 : 41986670 : end function phimag_vl3
1010 : :
1011 : 116829691 : function gsq_vl3(g1,g2,g3)
1012 : :
1013 : : real(dp) :: gsq_vl3
1014 : : real(dp),intent(in) :: g1,g2,g3 ! Note that they are real, unlike in other similar function definitions
1015 : : !Define G^2 based on G space metric gmet.
1016 : : gsq_vl3=g1*g1*gmet(1,1)+g2*g2*gmet(2,2)+&
1017 : : & g3*g3*gmet(3,3)+2.0_dp*g1*g2*gmet(1,2)+&
1018 : 116829691 : & 2.0_dp*g2*g3*gmet(2,3)+2.0_dp*g3*g1*gmet(3,1)
1019 : 116829691 : end function gsq_vl3
1020 : :
1021 : : end subroutine dfpt_vlocal
1022 : : !!***
1023 : :
1024 : : !!****f* ABINIT/vlocalstr
1025 : : !! NAME
1026 : : !! vlocalstr
1027 : : !!
1028 : : !! FUNCTION
1029 : : !! Compute strain derivatives of local ionic potential
1030 : : !! second derivative of E wrt xred
1031 : : !!
1032 : : !! INPUTS
1033 : : !! gmet(3,3)=reciprocal space metric ($\textrm{Bohr}^{-2}$).
1034 : : !! gprimd(3,3)=reciprocal space dimensional primitive translations
1035 : : !! gsqcut=cutoff on $|G|^2$: see setup1 for definition (doubled sphere).
1036 : : !! istr=1,...6 specifies cartesian strain component 11,22,33,32,31,21
1037 : : !! mgfft=maximum size of 1D FFTs
1038 : : !! mpi_enreg=information about MPI parallelization
1039 : : !! mqgrid=number of grid pts in q array for f(q) spline.
1040 : : !! natom=number of atoms in unit cell.
1041 : : !! nattyp(ntypat)=number of atoms of each type in cell.
1042 : : !! nfft=(effective) number of FFT grid points (for this processor)
1043 : : !! ngfft(18)=contain all needed information about 3D FFT,
1044 : : !! see ~abinit/doc/variables/vargs.htm#ngfft
1045 : : !! ntypat=number of types of atoms.
1046 : : !! ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phase information.
1047 : : !! qgrid(mqgrid)=q grid for spline from 0 to qmax.
1048 : : !! ucvol=unit cell volume ($\textrm{Bohr}^3$).
1049 : : !! vlspl(mqgrid,2,ntypat)=q^2 v(q) spline for each type of atom.
1050 : : !! [g0term]= optional, if present an alternative treatment of the G=0 term,
1051 : : !! adopoted for the flexoelectric tensor calculation, is performed.
1052 : : !!
1053 : : !! OUTPUT
1054 : : !! vpsp1(nfft)=first-order local crystal pseudopotential in real space.
1055 : : !!
1056 : : !! NOTES
1057 : : !! * Note that the present routine is tightly connected to the dfpt_vlocal.f routine,
1058 : : !! that compute the derivative of the local ionic potential
1059 : : !! with respect to one atomic displacement. The argument list
1060 : : !! and the internal loops to be considered were sufficiently different
1061 : : !! as to make the two routines different.
1062 : : !! * The routine was adapted from mklocl.F90
1063 : : !!
1064 : : !! SOURCE
1065 : :
1066 : 672 : subroutine vlocalstr(gmet,gprimd,gsqcut,istr,mgfft,mpi_enreg,&
1067 : 672 : & mqgrid,natom,nattyp,nfft,ngfft,ntypat,ph1d,qgrid,&
1068 : 672 : & ucvol,vlspl,vpsp1,g0term)
1069 : :
1070 : : !Arguments ------------------------------------
1071 : : !scalars
1072 : : integer,intent(in) :: istr,mgfft,mqgrid,natom,nfft,ntypat
1073 : : integer,optional,intent(in) :: g0term
1074 : : real(dp),intent(in) :: gsqcut,ucvol
1075 : : type(MPI_type),intent(in) :: mpi_enreg
1076 : : !arrays
1077 : : integer,intent(in) :: nattyp(ntypat),ngfft(18)
1078 : : real(dp),intent(in) :: gmet(3,3),gprimd(3,3),ph1d(2,3*(2*mgfft+1)*natom)
1079 : : real(dp),intent(in) :: qgrid(mqgrid),vlspl(mqgrid,2,ntypat)
1080 : : real(dp),intent(out) :: vpsp1(nfft)
1081 : :
1082 : : !Local variables-------------------------------
1083 : : !scalars
1084 : : integer,parameter :: im=2,re=1
1085 : : integer :: g0term_
1086 : : integer :: i1,i2,i3,ia,ia1,ia2,id1,id2,id3,ig1,ig2,ig3,ii,itypat,jj
1087 : : integer :: ka,kb,n1,n2,n3
1088 : : real(dp),parameter :: tolfix=1.0000001_dp
1089 : : real(dp) :: aa,bb,cc,cutoff,dd,dgsquards,diff
1090 : : real(dp) :: dq,dq2div6,dqdiv6,dqm1,ee,ff,gmag,gsquar
1091 : : real(dp) :: sfi,sfr,term,vion1,vion2,vlocg0
1092 : : real(dp) :: xnorm
1093 : : character(len=500) :: message
1094 : : !arrays
1095 : : integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/)
1096 : 672 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
1097 : 672 : integer, contiguous, pointer :: fftn3_distrib(:),ffti3_local(:)
1098 : : real(dp) :: dgmetds(3,3)
1099 : 672 : real(dp),allocatable :: work1(:,:)
1100 : : ! *************************************************************************
1101 : :
1102 : : !Define G^2 based on G space metric gmet.
1103 : : ! gsq_vl(i1,i2,i3)=dble(i1*i1)*gmet(1,1)+dble(i2*i2)*gmet(2,2)+&
1104 : : !& dble(i3*i3)*gmet(3,3)+dble(2*i1*i2)*gmet(1,2)+&
1105 : : !& dble(2*i2*i3)*gmet(2,3)+dble(2*i3*i1)*gmet(3,1)
1106 : :
1107 : : !Define dG^2/ds based on G space metric derivative dgmetds.
1108 : : ! dgsqds_vl(i1,i2,i3)=dble(i1*i1)*dgmetds(1,1)+dble(i2*i2)*dgmetds(2,2)+&
1109 : : !& dble(i3*i3)*dgmetds(3,3)+&
1110 : : !& dble(i1*i2)*(dgmetds(1,2)+dgmetds(2,1))+&
1111 : : !& dble(i1*i3)*(dgmetds(1,3)+dgmetds(3,1))+&
1112 : : !& dble(i2*i3)*(dgmetds(2,3)+dgmetds(3,2))
1113 : :
1114 : : !Real and imaginary parts of phase--statment functions:
1115 : : ! phr_vl(x1,y1,x2,y2,x3,y3)=(x1*x2-y1*y2)*x3-(y1*x2+x1*y2)*y3
1116 : : ! phi_vl(x1,y1,x2,y2,x3,y3)=(x1*x2-y1*y2)*y3+(y1*x2+x1*y2)*x3
1117 : : ! ph1_vl(nri,i1,ia)=ph1d(nri,i1+1+n1+(ia-1)*(2*n1+1))
1118 : : ! ph2_vl(nri,i2,ia)=ph1d(nri,i2+1+n2+(ia-1)*(2*n2+1)+&
1119 : : !& natom*(2*n1+1))
1120 : : ! ph3_vl(nri,i3,ia)=ph1d(nri,i3+1+n3+(ia-1)*(2*n3+1)+&
1121 : : !& natom*(2*n1+1+2*n2+1))
1122 : : ! phre_vl(i1,i2,i3,ia)=phr_vl(ph1_vl(re,i1,ia),ph1_vl(im,i1,ia),ph2_vl(re,i2,ia),&
1123 : : !& ph2_vl(im,i2,ia),ph3_vl(re,i3,ia),ph3_vl(im,i3,ia))
1124 : : ! phimag_vl(i1,i2,i3,ia)=phi_vl(ph1_vl(re,i1,ia),ph1_vl(im,i1,ia),ph2_vl(re,i2,ia),&
1125 : : !& ph2_vl(im,i2,ia),ph3_vl(re,i3,ia),ph3_vl(im,i3,ia))
1126 : :
1127 : : !-----
1128 : : !Compute derivative of metric tensor wrt strain component istr
1129 : 672 : if(istr<1 .or. istr>6)then
1130 : : write(message, '(a,i10,a,a,a)' )&
1131 : 0 : & ' Input istr=',istr,' not allowed.',ch10,&
1132 : 0 : & ' Possible values are 1,2,3,4,5,6 only.'
1133 : 0 : ABI_BUG(message)
1134 : : end if
1135 : :
1136 : 672 : ka=idx(2*istr-1);kb=idx(2*istr)
1137 [ + + ]: 2688 : do ii = 1,3
1138 [ + + ]: 8736 : dgmetds(:,ii)=-(gprimd(ka,:)*gprimd(kb,ii)+gprimd(kb,:)*gprimd(ka,ii))
1139 : : end do
1140 : : !For historical reasons:
1141 [ + + + + ]: 8736 : dgmetds(:,:)=0.5_dp*dgmetds(:,:)
1142 : :
1143 : 672 : n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
1144 : :
1145 : : !Get the distrib associated with this fft_grid
1146 : 672 : call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
1147 : :
1148 : : !Zero out array to permit accumulation over atom types below:
1149 [ - + + - : 2016 : ABI_MALLOC(work1,(2,nfft))
+ - + - -
+ - + ]
1150 [ + + + + ]: 8074584 : work1(:,:)=0.0_dp
1151 : : !
1152 : 672 : dq=(qgrid(mqgrid)-qgrid(1))/dble(mqgrid-1)
1153 : 672 : dqm1=1.0_dp/dq
1154 : 672 : dqdiv6=dq/6.0_dp
1155 : 672 : dq2div6=dq**2/6.0_dp
1156 : 672 : cutoff=gsqcut*tolfix
1157 : 672 : id1=n1/2+2
1158 : 672 : id2=n2/2+2
1159 : 672 : id3=n3/2+2
1160 : :
1161 : 672 : ia1=1
1162 [ + + ]: 1757 : do itypat=1,ntypat
1163 : : ! ia1,ia2 sets range of loop over atoms:
1164 : 1085 : ia2=ia1+nattyp(itypat)-1
1165 : :
1166 : 1085 : ii=0
1167 [ + + ]: 19859 : do i3=1,n3
1168 : 18774 : ig3=i3-(i3/id3)*n3-1
1169 [ + + ]: 303415 : do i2=1,n2
1170 [ + - ]: 302330 : if (fftn2_distrib(i2)==mpi_enreg%me_fft) then
1171 : 283556 : ig2=i2-(i2/id2)*n2-1
1172 [ + + ]: 4723124 : do i1=1,n1
1173 : 4439568 : ig1=i1-(i1/id1)*n1-1
1174 : 4439568 : ii=ii+1
1175 : : ! *** GET RID OF THIS THESE IF STATEMENTS (if they slow code)
1176 : : ! Skip G=0:
1177 : : ! if (ii==1) cycle
1178 [ + + ]: 4439568 : if (ig1==0 .and. ig2==0 .and. ig3==0) cycle
1179 : 4438483 : gsquar=gsq_vl(ig1,ig2,ig3)
1180 : :
1181 : : ! Skip G**2 outside cutoff:
1182 [ + + ]: 4722039 : if (gsquar<=cutoff) then
1183 : 1420174 : gmag=sqrt(gsquar)
1184 : 1420174 : dgsquards=dgsqds_vl(ig1,ig2,ig3)
1185 : : ! Compute vion(G) for given type of atom
1186 : 1420174 : jj=1+int(gmag*dqm1)
1187 : 1420174 : diff=gmag-qgrid(jj)
1188 : :
1189 : : ! Evaluate spline fit from q^2 V(q) to get V(q):
1190 : : ! (p. 86 Numerical Recipes, Press et al;
1191 : : ! NOTE error in book for sign
1192 : : ! of "aa" term in derivative; also see splfit routine).
1193 : :
1194 : 1420174 : bb = diff*dqm1
1195 : 1420174 : aa = 1.0_dp-bb
1196 : 1420174 : cc = aa*(aa**2-1.0_dp)*dq2div6
1197 : 1420174 : dd = bb*(bb**2-1.0_dp)*dq2div6
1198 : :
1199 : : vion1 = (aa*vlspl(jj,1,itypat)+bb*vlspl(jj+1,1,itypat) +&
1200 : : & cc*vlspl(jj,2,itypat)+dd*vlspl(jj+1,2,itypat) ) &
1201 : 1420174 : & / gsquar
1202 : :
1203 : : ! Also get (dV(q)/dq)/q:
1204 : : ! (note correction of Numerical Recipes sign error
1205 : : ! before (3._dp*aa**2-1._dp)
1206 : 1420174 : ee= vlspl(jj+1,1,itypat)-vlspl(jj,1,itypat)
1207 : : ff= (3._dp*bb**2-1._dp)*vlspl(jj+1,2,itypat) &
1208 : 1420174 : & - (3._dp*aa**2-1._dp)*vlspl(jj,2,itypat)
1209 : : vion2 = ( ( ee*dqm1 + ff*dqdiv6 )/gmag&
1210 : 1420174 : & - 2.0_dp*vion1 ) / gsquar
1211 : :
1212 : :
1213 : : ! Assemble structure factor over all atoms of given type:
1214 : 1420174 : sfr=0.0_dp
1215 : 1420174 : sfi=0.0_dp
1216 [ + + ]: 3609368 : do ia=ia1,ia2
1217 : 2189194 : sfr=sfr+phre_vl(ig1,ig2,ig3,ia)
1218 : 3609368 : sfi=sfi-phimag_vl(ig1,ig2,ig3,ia)
1219 : : end do
1220 : :
1221 : 1420174 : term=dgsquards*vion2
1222 : : ! Add potential for diagonal strain components
1223 [ + + ]: 1420174 : if(istr <=3) then
1224 : 718770 : term=term-vion1
1225 : : end if
1226 : :
1227 : : ! Multiply structure factor times vion derivatives:
1228 : 1420174 : work1(re,ii)=work1(re,ii)+sfr*term
1229 : 1420174 : work1(im,ii)=work1(im,ii)+sfi*term
1230 : :
1231 : : ! End skip G**2 outside cutoff:
1232 : : end if
1233 : : ! End loop on n1, n2, n3. There is a "cycle" inside the loop
1234 : : end do
1235 : : end if
1236 : : end do
1237 : : end do
1238 : :
1239 : 1757 : ia1=ia2+1
1240 : :
1241 : : ! End loop on type of atoms
1242 : : end do
1243 : :
1244 : :
1245 : : !Set Vloc(G=0)=0:
1246 : 672 : work1(re,1)=0.0_dp
1247 : 672 : work1(im,1)=0.0_dp
1248 : :
1249 : : !Alternative treatment of Vloc(G=0) for the flexoelectric tensor calculation
1250 [ + - ]: 672 : g0term_=0; if (present(g0term)) g0term_=g0term
1251 [ + + ]: 672 : if (g0term_==1) then
1252 : 126 : vlocg0=zero
1253 [ + + ]: 126 : if (istr<=3) then
1254 : 126 : ia1=1
1255 [ + + ]: 126 : do itypat=1,ntypat
1256 : : ! ia1,ia2 sets range of loop over atoms:
1257 : :
1258 : 63 : ia2=ia1+nattyp(itypat)-1
1259 [ + + ]: 252 : do ia=ia1,ia2
1260 : 189 : vlocg0=vlocg0+vlspl(1,2,itypat)
1261 : : end do
1262 : : end do
1263 : 63 : work1(re,1)=-half*vlocg0
1264 : : end if
1265 : : end if
1266 : :
1267 : : !Transform back to real space
1268 : 672 : call fourdp(1,work1,vpsp1,1,mpi_enreg,nfft,1,ngfft,0)
1269 : :
1270 : : !Divide by unit cell volume
1271 : 672 : xnorm=1.0_dp/ucvol
1272 [ + + ]: 2691976 : vpsp1(:)=vpsp1(:)*xnorm
1273 : :
1274 [ - + ]: 672 : ABI_FREE(work1)
1275 : :
1276 : : contains
1277 : :
1278 : : !Real and imaginary parts of phase.
1279 : 2189194 : function phr_vl(x1,y1,x2,y2,x3,y3)
1280 : :
1281 : : real(dp) :: phr_vl,x1,x2,x3,y1,y2,y3
1282 : 2189194 : phr_vl=(x1*x2-y1*y2)*x3-(y1*x2+x1*y2)*y3
1283 : : end function phr_vl
1284 : :
1285 : 2189194 : function phi_vl(x1,y1,x2,y2,x3,y3)
1286 : :
1287 : : real(dp):: phi_vl,x1,x2,x3,y1,y2,y3
1288 : 2189194 : phi_vl=(x1*x2-y1*y2)*y3+(y1*x2+x1*y2)*x3
1289 : : end function phi_vl
1290 : :
1291 : 8756776 : function ph1_vl(nri,ig1,ia)
1292 : :
1293 : : real(dp):: ph1_vl
1294 : : integer :: nri,ig1,ia
1295 : 8756776 : ph1_vl=ph1d(nri,ig1+1+n1+(ia-1)*(2*n1+1))
1296 : 8756776 : end function ph1_vl
1297 : :
1298 : 8756776 : function ph2_vl(nri,ig2,ia)
1299 : :
1300 : : real(dp):: ph2_vl
1301 : : integer :: nri,ig2,ia
1302 : 8756776 : ph2_vl=ph1d(nri,ig2+1+n2+(ia-1)*(2*n2+1)+natom*(2*n1+1))
1303 : 8756776 : end function ph2_vl
1304 : :
1305 : 8756776 : function ph3_vl(nri,ig3,ia)
1306 : :
1307 : : real(dp):: ph3_vl
1308 : : integer :: nri,ig3,ia
1309 : 8756776 : ph3_vl=ph1d(nri,ig3+1+n3+(ia-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1))
1310 : 8756776 : end function ph3_vl
1311 : :
1312 : 2189194 : function phre_vl(ig1,ig2,ig3,ia)
1313 : :
1314 : : real(dp):: phre_vl
1315 : : integer :: ig1,ig2,ig3,ia
1316 : : phre_vl=phr_vl(ph1_vl(re,ig1,ia),ph1_vl(im,ig1,ia),&
1317 : 2189194 : & ph2_vl(re,ig2,ia),ph2_vl(im,ig2,ia),ph3_vl(re,ig3,ia),ph3_vl(im,ig3,ia))
1318 : 2189194 : end function phre_vl
1319 : :
1320 : 2189194 : function phimag_vl(ig1,ig2,ig3,ia)
1321 : :
1322 : : real(dp) :: phimag_vl
1323 : : integer :: ig1,ig2,ig3,ia
1324 : : phimag_vl=phi_vl(ph1_vl(re,ig1,ia),ph1_vl(im,ig1,ia),&
1325 : 2189194 : & ph2_vl(re,ig2,ia),ph2_vl(im,ig2,ia),ph3_vl(re,ig3,ia),ph3_vl(im,ig3,ia))
1326 : 2189194 : end function phimag_vl
1327 : :
1328 : 4438483 : function gsq_vl(i1,i2,i3)
1329 : :
1330 : : real(dp) :: gsq_vl
1331 : : integer :: i1,i2,i3
1332 : : !Define G^2 based on G space metric gmet.
1333 : : gsq_vl=dble(i1*i1)*gmet(1,1)+dble(i2*i2)*gmet(2,2)+&
1334 : : & dble(i3*i3)*gmet(3,3)+dble(2*i1*i2)*gmet(1,2)+&
1335 : 4438483 : & dble(2*i2*i3)*gmet(2,3)+dble(2*i3*i1)*gmet(3,1)
1336 : 4438483 : end function gsq_vl
1337 : :
1338 : 1420174 : function dgsqds_vl(i1,i2,i3)
1339 : :
1340 : : real(dp) :: dgsqds_vl
1341 : : integer :: i1,i2,i3
1342 : : !Define dG^2/ds based on G space metric derivative dgmetds.
1343 : : dgsqds_vl=dble(i1*i1)*dgmetds(1,1)+dble(i2*i2)*dgmetds(2,2)+&
1344 : : & dble(i3*i3)*dgmetds(3,3)+&
1345 : : & dble(i1*i2)*(dgmetds(1,2)+dgmetds(2,1))+&
1346 : : & dble(i1*i3)*(dgmetds(1,3)+dgmetds(3,1))+&
1347 : 1420174 : & dble(i2*i3)*(dgmetds(2,3)+dgmetds(3,2))
1348 : 1420174 : end function dgsqds_vl
1349 : :
1350 : : end subroutine vlocalstr
1351 : : !!***
1352 : :
1353 : : !!****f* ABINIT/dfpt_vlocaldq
1354 : : !! NAME
1355 : : !! dfpt_vlocaldq
1356 : : !!
1357 : : !! FUNCTION
1358 : : !! Compute q-gradient (at q=0) of the local part of 1st-order
1359 : : !! atomic displacement local potential or pseudocore density from the appropriate
1360 : : !! atomic pseudopotential data with structure and derivative factor.
1361 : : !!
1362 : : !! INPUTS
1363 : : !! atindx(natom)=index table for atoms (see gstate.f)
1364 : : !! cplex: if 1, real space 1-order functions on FFT grid
1365 : : !! are REAL, if 2, COMPLEX
1366 : : !! gmet(3,3)=reciprocal space metric (Bohr**-2)
1367 : : !! gsqcut=cutoff G**2 for included G s in fft box.
1368 : : !! idir=direction of atomic displacement (=1,2 or 3 : displacement of
1369 : : !! atom ipert along the 1st, 2nd or 3rd axis).
1370 : : !! ipert=number of the atom being displaced in the frozen-phonon
1371 : : !! mpi_enreg=information about MPI parallelization
1372 : : !! mqgrid=dimension of q grid for pseudopotentials
1373 : : !! natom=number of atoms in cell.
1374 : : !! nattyp(ntypat)=number of atoms of each type in cell.
1375 : : !! nfft=(effective) number of FFT grid points (for this processor)
1376 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/input_variables/vargs.htm#ngfft
1377 : : !! ntypat=number of types of atoms in cell.
1378 : : !! n1,n2,n3=fft grid.
1379 : : !! ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phase information.
1380 : : !! qdir=direction of the q-gradient
1381 : : !! qgrid(mqgrid)=grid of q points from 0 to qmax.
1382 : : !! qphon(3)=wavevector of the phonon
1383 : : !! ucvol=unit cell volume (Bohr**3).
1384 : : !! vlspl(mqgrid,2,ntypat)=spline fit of q^2 V(q) for each type of atom.
1385 : : !! if optnc=1-> spline fit of n_c(q) for each type of atom.
1386 : : !! optnc= (OPTIONAL) if present and equal to 1, does the calculation assuming a
1387 : : !! pseudocore density
1388 : : !!
1389 : : !! OUTPUT
1390 : : !! vpsp1dq(cplex*nfft)=q-gradient (at q=0) of the first-order local
1391 : : !! crystal pseudopotential in real space
1392 : : !! (including the minus sign, forgotten in the paper non-linear..
1393 : : !!
1394 : : !! NOTES
1395 : : !! * IMPORTANT: the formalism followed in this routine
1396 : : !! assumes a phase factor for the perturbation that
1397 : : !! is different to the one used elsewhere in the code (See M.Stengel paper):
1398 : : !!
1399 : : !! here: e^{i q (R_l + \tau_{\kappa})}
1400 : : !! rest of ABINIT: e^{i q R_l}
1401 : : !!
1402 : : !! **A -i factor has been factorized out in all the contributions of the first
1403 : : !! q-gradient of the atomic displacement Hamiltonian. This is lately included
1404 : : !! in the matrix element calculation.
1405 : : !!
1406 : : !! SOURCE
1407 : :
1408 : 799 : subroutine dfpt_vlocaldq(atindx,cplex,gmet,gsqcut,idir,ipert,&
1409 : 799 : & mpi_enreg,mqgrid,natom,nattyp,nfft,ngfft,&
1410 : 799 : & ntypat,n1,n2,n3,ph1d,qdir,qgrid,qphon,ucvol,vlspl,vpsp1dq,&
1411 : : & optnc) !optional
1412 : :
1413 : : !Arguments -------------------------------
1414 : : !scalars
1415 : : integer,intent(in) :: cplex,idir,ipert,mqgrid,n1,n2,n3,natom,nfft,ntypat
1416 : : integer,intent(in) :: qdir
1417 : : integer,optional,intent(in) :: optnc
1418 : : real(dp),intent(in) :: gsqcut,ucvol
1419 : : type(MPI_type),intent(in) :: mpi_enreg
1420 : : !arrays
1421 : : integer,intent(in) :: atindx(natom),nattyp(ntypat),ngfft(18)
1422 : : real(dp),intent(in) :: gmet(3,3),ph1d(2,(2*n1+1+2*n2+1+2*n3+1)*natom)
1423 : : real(dp),intent(in) :: qgrid(mqgrid),qphon(3),vlspl(mqgrid,2,ntypat)
1424 : : real(dp),intent(out) :: vpsp1dq(cplex*nfft)
1425 : :
1426 : : !Local variables -------------------------
1427 : : !scalars
1428 : : integer :: i1,i2,i3,ia1,iatom,id1,id2,id3,ig1,ig2,ig3,ii,ii1
1429 : : integer :: itypat, optnc_
1430 : : integer, parameter :: im=2, re=1
1431 : : real(dp),parameter :: tolfix=1.000000001_dp
1432 : : real(dp) :: cutoff,gfact,gmag,gq1
1433 : : real(dp) :: gq2,gq3,gsquar
1434 : : real(dp) :: sfi,sfr,xnorm
1435 : : logical :: qeq0
1436 : : character(len=500) :: msg
1437 : : !arrays
1438 : 799 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
1439 : 799 : integer, contiguous, pointer :: fftn3_distrib(:),ffti3_local(:)
1440 : : real(dp) :: gq(3),gvec(3),vion1(1),vion1dq(1)
1441 : 799 : real(dp),allocatable :: work1(:,:)
1442 : : ! *********************************************************************
1443 : :
1444 : 799 : iatom=ipert
1445 : :
1446 : 27 : optnc_=0; if (present(optnc)) optnc_=optnc
1447 : : if(iatom==natom+1 .or. iatom==natom+2 .or. iatom==natom+10 .or. iatom==natom+11 &
1448 [ + + - + : 826 : & .or. iatom==natom+5 .or. iatom==natom+6 .or. (iatom>natom+11.and.iatom<=2*natom+11))then
- + - + -
+ - + - +
- + - - ]
1449 : :
1450 : : ! (In case of d/dk or an electric field, or magnetic (Zeeman) field->[natom+5] SPr deb )
1451 [ # # ]: 0 : vpsp1dq(1:cplex*nfft)=zero
1452 : :
1453 : : else
1454 : :
1455 : : ! (In case of a phonon perturbation)
1456 [ - + + - : 2397 : ABI_MALLOC(work1,(2,nfft))
+ - + - -
+ - + ]
1457 [ + + + + ]: 10394170 : work1(1:2,1:nfft)=0.0_dp
1458 : :
1459 : 799 : cutoff=gsqcut*tolfix
1460 : 799 : id1=n1/2+2
1461 : 799 : id2=n2/2+2
1462 : 799 : id3=n3/2+2
1463 : :
1464 : : ! Get the distrib associated with this fft_grid
1465 : 799 : call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
1466 : :
1467 : : ! This is to allow q=0
1468 : 799 : qeq0=.false.
1469 [ - + ]: 799 : if(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15) then
1470 : : qeq0=.true.
1471 : : else
1472 : 0 : msg='This routine cannot be used for q/=0'
1473 : 0 : ABI_BUG(msg)
1474 : : end if
1475 : :
1476 : : ! Determination of the atom type
1477 : 799 : ia1=0
1478 : 799 : itypat=0
1479 [ + + ]: 1728 : do ii=1,ntypat
1480 : 929 : ia1=ia1+nattyp(ii)
1481 [ + + + + ]: 1858 : if(atindx(iatom)<=ia1.and.itypat==0)itypat=ii
1482 : : end do
1483 : :
1484 : : ii=0
1485 : :
1486 [ + + ]: 13216 : do i3=1,n3
1487 : 12417 : ig3=i3-(i3/id3)*n3-1
1488 : 12417 : gq3=dble(ig3)+qphon(3)
1489 : 12417 : gq(3)=gq3
1490 [ + + ]: 213223 : do i2=1,n2
1491 [ + - ]: 212424 : if (fftn2_distrib(i2)==mpi_enreg%me_fft) then
1492 : 200007 : ig2=i2-(i2/id2)*n2-1
1493 : 200007 : gq2=dble(ig2)+qphon(2)
1494 : 200007 : gq(2)=gq2
1495 : :
1496 : : ! Note the lower limit of the next loop
1497 : 200007 : ii1=1
1498 [ + + + - : 200007 : if(i3==1 .and. i2==1 .and. qeq0 .and. ig2==0 .and. ig3==0)then
+ - ]
1499 : 799 : ii1=2
1500 : 799 : ii=ii+1
1501 : : end if
1502 [ + + ]: 3663665 : do i1=ii1,n1
1503 : 3463658 : ig1=i1-(i1/id1)*n1-1
1504 : 3463658 : gq1=dble(ig1)+qphon(1)
1505 : 3463658 : gq(1)=gq1
1506 : 3463658 : ii=ii+1
1507 : 3463658 : gsquar=gsq_vl3(gq1,gq2,gq3)
1508 : : ! Skip G**2 outside cutoff:
1509 [ + + ]: 3663665 : if (gsquar<=cutoff) then
1510 : 958320 : gmag=sqrt(gsquar)
1511 : :
1512 : : ! Evaluate spline fit to get V(q) and V(q)':
1513 [ + + ]: 1916640 : call splfit(qgrid,vion1dq,vlspl(:,:,itypat),1,(/gmag/),vion1,mqgrid,1)
1514 : :
1515 [ + + ]: 958320 : if (optnc_==0) then
1516 [ + + ]: 1814616 : vion1=vion1/gsquar
1517 [ + + ]: 1814616 : vion1dq=(vion1dq-2.0_dp*gmag*vion1)/gsquar
1518 : : end if
1519 : :
1520 [ + + ]: 3833280 : gvec=(/ig1,ig2,ig3/)
1521 [ + + ]: 3833280 : gfact=dot_product(gmet(qdir,:),gvec(:))/gmag
1522 : :
1523 : : ! Phase G*xred (complex conjugate) *2*pi*(g_idir)*vion1dq*gfact
1524 : 958320 : sfr=phre_vl3(ig1,ig2,ig3,iatom)*two_pi*gq(idir)*vion1dq(1)*gfact
1525 : 958320 : sfi=-phimag_vl3(ig1,ig2,ig3,iatom)*two_pi*gq(idir)*vion1dq(1)*gfact
1526 : : ! Phase G*xred (complex conjugate) *2*pi*(\delta_{idir,qdir})*vion1
1527 [ + + ]: 958320 : if (idir==qdir) then
1528 : 390504 : sfr=sfr + phre_vl3(ig1,ig2,ig3,iatom)*two_pi*vion1(1)
1529 : 390504 : sfi=sfi - phimag_vl3(ig1,ig2,ig3,iatom)*two_pi*vion1(1)
1530 : : end if
1531 : :
1532 : 958320 : work1(re,ii)=sfr
1533 : 958320 : work1(im,ii)=sfi
1534 : : end if
1535 : :
1536 : : end do
1537 : : end if
1538 : : end do
1539 : : end do
1540 : :
1541 : : ! Consider the macroscopic term in the case of the pseudocore density
1542 [ + + + + ]: 799 : if (optnc_==1.and.idir==qdir) then
1543 : 11 : work1(re,1)= two_pi*vlspl(1,1,itypat)
1544 : : end if
1545 : :
1546 : : ! Transform back to real space
1547 : 799 : call fourdp(cplex,work1,vpsp1dq,1,mpi_enreg,nfft,1,ngfft,0)
1548 : :
1549 : 799 : xnorm=1.0_dp/ucvol
1550 [ + + ]: 6929713 : vpsp1dq(1:cplex*nfft)=vpsp1dq(1:cplex*nfft)*xnorm
1551 : :
1552 [ - + ]: 799 : ABI_FREE(work1)
1553 : :
1554 : : ! End the condition of non-electric-field
1555 : : end if
1556 : :
1557 : : contains
1558 : :
1559 : : !Real and imaginary parts of phase.
1560 : 958320 : function phr_vl3(x1,y1,x2,y2,x3,y3)
1561 : : real(dp) :: phr_vl3
1562 : : real(dp),intent(in) :: x1,x2,x3,y1,y2,y3
1563 : 958320 : phr_vl3=(x1*x2-y1*y2)*x3-(y1*x2+x1*y2)*y3
1564 : : end function phr_vl3
1565 : :
1566 : 958320 : function phi_vl3(x1,y1,x2,y2,x3,y3)
1567 : : real(dp) :: phi_vl3
1568 : : real(dp),intent(in) :: x1,x2,x3,y1,y2,y3
1569 : 958320 : phi_vl3=(x1*x2-y1*y2)*y3+(y1*x2+x1*y2)*x3
1570 : : end function phi_vl3
1571 : :
1572 : : ! Warning : this function differ from similar ones for ground-state calculations : note the atindx !!
1573 : 3833280 : function ph1_vl3(nri,ig1,ia)
1574 : : real(dp) :: ph1_vl3
1575 : : integer,intent(in) :: nri,ig1,ia
1576 : 3833280 : ph1_vl3=ph1d(nri,ig1+1+n1+(atindx(ia)-1)*(2*n1+1))
1577 : 3833280 : end function ph1_vl3
1578 : :
1579 : : ! Warning : this function differ from similar ones for ground-state calculations : note the atindx !!
1580 : 3833280 : function ph2_vl3(nri,ig2,ia)
1581 : : real(dp) :: ph2_vl3
1582 : : integer,intent(in) :: nri,ig2,ia
1583 : 3833280 : ph2_vl3=ph1d(nri,ig2+1+n2+(atindx(ia)-1)*(2*n2+1)+natom*(2*n1+1))
1584 : 3833280 : end function ph2_vl3
1585 : :
1586 : : ! Warning : this function differ from similar ones for ground-state calculations : note the atindx !!
1587 : 3833280 : function ph3_vl3(nri,ig3,ia)
1588 : : real(dp) :: ph3_vl3
1589 : : integer,intent(in) :: nri,ig3,ia
1590 : 3833280 : ph3_vl3=ph1d(nri,ig3+1+n3+(atindx(ia)-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1))
1591 : 3833280 : end function ph3_vl3
1592 : :
1593 : 958320 : function phre_vl3(ig1,ig2,ig3,ia)
1594 : : real(dp) :: phre_vl3
1595 : : integer,intent(in) :: ig1,ig2,ig3,ia
1596 : : phre_vl3=phr_vl3(ph1_vl3(re,ig1,ia),ph1_vl3(im,ig1,ia),&
1597 : 958320 : & ph2_vl3(re,ig2,ia),ph2_vl3(im,ig2,ia),ph3_vl3(re,ig3,ia),ph3_vl3(im,ig3,ia))
1598 : 958320 : end function phre_vl3
1599 : :
1600 : 958320 : function phimag_vl3(ig1,ig2,ig3,ia)
1601 : : real(dp) :: phimag_vl3
1602 : : integer,intent(in) :: ig1,ig2,ig3,ia
1603 : : phimag_vl3=phi_vl3(ph1_vl3(re,ig1,ia),ph1_vl3(im,ig1,ia),&
1604 : 958320 : & ph2_vl3(re,ig2,ia),ph2_vl3(im,ig2,ia),ph3_vl3(re,ig3,ia),ph3_vl3(im,ig3,ia))
1605 : 958320 : end function phimag_vl3
1606 : :
1607 : 3463658 : function gsq_vl3(g1,g2,g3)
1608 : : real(dp) :: gsq_vl3
1609 : : real(dp),intent(in) :: g1,g2,g3 ! Note that they are real, unlike in other similar function definitions
1610 : : !Define G^2 based on G space metric gmet.
1611 : : gsq_vl3=g1*g1*gmet(1,1)+g2*g2*gmet(2,2)+&
1612 : : & g3*g3*gmet(3,3)+2.0_dp*g1*g2*gmet(1,2)+&
1613 : 3463658 : & 2.0_dp*g2*g3*gmet(2,3)+2.0_dp*g3*g1*gmet(3,1)
1614 : 3463658 : end function gsq_vl3
1615 : :
1616 : : end subroutine dfpt_vlocaldq
1617 : : !!***
1618 : :
1619 : : !!****f* ABINIT/dfpt_vlocaldqdq
1620 : : !! NAME
1621 : : !! dfpt_vlocaldqdq
1622 : : !!
1623 : : !! FUNCTION
1624 : : !! Compute 2nd q-gradient (at q=0) of the local part of 1st-order
1625 : : !! atomic displacement potential from the appropriate
1626 : : !! atomic pseudopotential with structure and derivative factor.
1627 : : !!
1628 : : !! INPUTS
1629 : : !! atindx(natom)=index table for atoms (see gstate.f)
1630 : : !! cplex: if 1, real space 1-order functions on FFT grid
1631 : : !! are REAL, if 2, COMPLEX
1632 : : !! gmet(3,3)=reciprocal space metric (Bohr**-2)
1633 : : !! gsqcut=cutoff G**2 for included G s in fft box.
1634 : : !! idir=direction of atomic displacement (=1,2 or 3 : displacement of
1635 : : !! atom ipert along the 1st, 2nd or 3rd axis).
1636 : : !! ipert=number of the atom being displaced in the frozen-phonon
1637 : : !! mpi_enreg=information about MPI parallelization
1638 : : !! mqgrid=dimension of q grid for pseudopotentials
1639 : : !! natom=number of atoms in cell.
1640 : : !! nattyp(ntypat)=number of atoms of each type in cell.
1641 : : !! nfft=(effective) number of FFT grid points (for this processor)
1642 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/input_variables/vargs.htm#ngfft
1643 : : !! ntypat=number of types of atoms in cell.
1644 : : !! n1,n2,n3=fft grid.
1645 : : !! ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phase information.
1646 : : !! qdir1=direction of the first q-gradient
1647 : : !! qdir2=direction of the second q-gradient
1648 : : !! qgrid(mqgrid)=grid of q points from 0 to qmax.
1649 : : !! qphon(3)=wavevector of the phonon
1650 : : !! ucvol=unit cell volume (Bohr**3).
1651 : : !! vlspl(mqgrid,2,ntypat)=spline fit of q^2 V(q) for each type of atom.
1652 : : !!
1653 : : !! OUTPUT
1654 : : !! vpsp1dqdq(cplex*nfft)=2nd q-gradient (at q=0) of the first-order local
1655 : : !! crystal pseudopotential in real space
1656 : : !!
1657 : : !! NOTES
1658 : : !! * IMPORTANT: the formalism followed in this routine
1659 : : !! assumes a phase factor for the perturbation that
1660 : : !! is different to the one used elsewhere in the code (See M.Stengel paper):
1661 : : !!
1662 : : !! here: e^{i q (R_l + \tau_{\kappa})}
1663 : : !! rest of ABINIT: e^{i q R_l}
1664 : : !!
1665 : : !! **A -i factor has been factorized out in all the contributions of the second
1666 : : !! q-gradient of the atomic displacement Hamiltonian. This is lately included
1667 : : !! in the whole frozen contribution to the q-gradient of the
1668 : : !! 2nd order energy wrt an atomic displacement and a strain:
1669 : : !! \Delta E^{\tau_{\kappa\alpha}^* (\beta)}_{m\kvec,\gamma\delta}
1670 : : !!
1671 : : !!
1672 : : !! SOURCE
1673 : :
1674 : 486 : subroutine dfpt_vlocaldqdq(atindx,cplex,gmet,gsqcut,idir,ipert,&
1675 : 486 : & mpi_enreg,mqgrid,natom,nattyp,nfft,ngfft,&
1676 : 486 : & ntypat,n1,n2,n3,ph1d,qdir1,qdir2,qgrid,qphon,ucvol,vlspl,vpsp1dqdq)
1677 : :
1678 : : !Arguments -------------------------------
1679 : : !scalars
1680 : : integer,intent(in) :: cplex,idir,ipert,mqgrid,n1,n2,n3,natom,nfft,ntypat
1681 : : integer,intent(in) :: qdir1,qdir2
1682 : : real(dp),intent(in) :: gsqcut,ucvol
1683 : : type(MPI_type),intent(in) :: mpi_enreg
1684 : : !arrays
1685 : : integer,intent(in) :: atindx(natom),nattyp(ntypat),ngfft(18)
1686 : : real(dp),intent(in) :: gmet(3,3),ph1d(2,(2*n1+1+2*n2+1+2*n3+1)*natom)
1687 : : real(dp),intent(in) :: qgrid(mqgrid),qphon(3),vlspl(mqgrid,2,ntypat)
1688 : : real(dp),intent(out) :: vpsp1dqdq(cplex*nfft)
1689 : :
1690 : : !Local variables -------------------------
1691 : : !scalars
1692 : : integer :: alpha, delta, gamma
1693 : : integer :: i1,i2,i3,ia1,iatom,id1,id2,id3,ig1,ig2,ig3,ii,ii1,im=2
1694 : : integer :: itypat,re=1
1695 : : real(dp),parameter :: tolfix=1.000000001_dp
1696 : : real(dp) :: cutoff,delad,delag,gfact,gfact1,gfact2,gmag,gq1
1697 : : real(dp) :: gq2,gq3,gsquar
1698 : : real(dp) :: sfi,sfr,term1,term2,xnorm
1699 : : logical :: qeq0
1700 : : character(len=500) :: msg
1701 : : !arrays
1702 : 486 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
1703 : 486 : integer, contiguous, pointer :: fftn3_distrib(:),ffti3_local(:)
1704 : : real(dp) :: gq(3),gvec(3),vion1(1),vion1dq(1),vion1dqdq(1)
1705 : 486 : real(dp),allocatable :: work1(:,:)
1706 : : ! *********************************************************************
1707 : :
1708 : 486 : iatom=ipert
1709 : :
1710 : : if(iatom==natom+1 .or. iatom==natom+2 .or. iatom==natom+10 .or. iatom==natom+11 &
1711 [ - + - + : 486 : & .or. iatom==natom+5 .or. iatom==natom+6 .or. (iatom>natom+11.and.iatom<=2*natom+11))then
- + - + -
+ - + - +
- - ]
1712 : :
1713 : : ! (In case of d/dk or an electric field, or magnetic (Zeeman) field->[natom+5] SPr deb )
1714 [ # # ]: 0 : vpsp1dqdq(1:cplex*nfft)=zero
1715 : :
1716 : : else
1717 : :
1718 : 486 : alpha=idir; delta=qdir2; gamma=qdir1
1719 : :
1720 : : !Kronecker deltas
1721 : 486 : delad=0.0_dp; delag=0.0_dp
1722 [ + + ]: 486 : if (alpha==delta) delad=1.0_dp
1723 [ + + ]: 486 : if (alpha==gamma) delag=1.0_dp
1724 : :
1725 : : ! (In case of a phonon perturbation)
1726 [ - + + - : 1458 : ABI_MALLOC(work1,(2,nfft))
+ - + - -
+ - + ]
1727 [ + + + + ]: 4921236 : work1(1:2,1:nfft)=0.0_dp
1728 : :
1729 : 486 : cutoff=gsqcut*tolfix
1730 : 486 : id1=n1/2+2
1731 : 486 : id2=n2/2+2
1732 : 486 : id3=n3/2+2
1733 : :
1734 : : ! Get the distrib associated with this fft_grid
1735 : 486 : call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
1736 : :
1737 : : ! This is to allow q=0
1738 : 486 : qeq0=.false.
1739 [ - + ]: 486 : if(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15) then
1740 : : qeq0=.true.
1741 : : else
1742 : 0 : msg='This routine cannot be used for q/=0'
1743 : 0 : ABI_BUG(msg)
1744 : : end if
1745 : :
1746 : : ! Determination of the atom type
1747 : 486 : ia1=0
1748 : 486 : itypat=0
1749 [ + + ]: 972 : do ii=1,ntypat
1750 : 486 : ia1=ia1+nattyp(ii)
1751 [ - + - + ]: 972 : if(atindx(iatom)<=ia1.and.itypat==0)itypat=ii
1752 : : end do
1753 : :
1754 : : ii=0
1755 : :
1756 [ + + ]: 7776 : do i3=1,n3
1757 : 7290 : ig3=i3-(i3/id3)*n3-1
1758 : 7290 : gq3=dble(ig3)+qphon(3)
1759 : 7290 : gq(3)=gq3
1760 [ + + ]: 117126 : do i2=1,n2
1761 [ + - ]: 116640 : if (fftn2_distrib(i2)==mpi_enreg%me_fft) then
1762 : 109350 : ig2=i2-(i2/id2)*n2-1
1763 : 109350 : gq2=dble(ig2)+qphon(2)
1764 : 109350 : gq(2)=gq2
1765 : :
1766 : : ! Note the lower limit of the next loop
1767 : 109350 : ii1=1
1768 [ + + + - : 109350 : if(i3==1 .and. i2==1 .and. qeq0 .and. ig2==0 .and. ig3==0)then
+ - ]
1769 : 486 : ii1=2
1770 : 486 : ii=ii+1
1771 : : end if
1772 [ + + ]: 1749114 : do i1=ii1,n1
1773 : 1639764 : ig1=i1-(i1/id1)*n1-1
1774 : 1639764 : gq1=dble(ig1)+qphon(1)
1775 : 1639764 : gq(1)=gq1
1776 : 1639764 : ii=ii+1
1777 : 1639764 : gsquar=gsq_vl3(gq1,gq2,gq3)
1778 : : ! Skip G**2 outside cutoff:
1779 [ + + ]: 1749114 : if (gsquar<=cutoff) then
1780 : 363528 : gmag=sqrt(gsquar)
1781 : :
1782 : : ! Evaluate spline fit to get first V(q) and V(q)' and later V(q)'':
1783 [ + + ]: 727056 : call splfit(qgrid,vion1dq,vlspl(:,:,itypat),1,(/gmag/),vion1,mqgrid,1)
1784 [ + + ]: 727056 : vion1=vion1/gsquar
1785 [ + + ]: 727056 : vion1dq=(vion1dq-2.0_dp*gmag*vion1)/gsquar
1786 : :
1787 [ + + ]: 727056 : call splfit(qgrid,vion1dqdq,vlspl(:,:,itypat),2,(/gmag/),vion1,mqgrid,1)
1788 [ + + ]: 727056 : vion1dqdq=(vion1dqdq-4.0_dp*gmag*vion1dq-2.0_dp*vion1)/gsquar
1789 : :
1790 [ + + ]: 1454112 : gvec=(/ig1,ig2,ig3/)
1791 [ + + ]: 1454112 : gfact1=dot_product(gmet(gamma,:),gvec(:))
1792 [ + + ]: 1454112 : gfact2=dot_product(gmet(delta,:),gvec(:))
1793 : 363528 : gfact=gvec(alpha)*gfact1*gfact2/gsquar
1794 : :
1795 : 363528 : term1=delag*gfact2+delad*gfact1+gvec(alpha)*gmet(gamma,delta)
1796 : 363528 : term1=term1-gfact
1797 : 363528 : term1=term1*vion1dq(1)/gmag
1798 : :
1799 : 363528 : term2=vion1dqdq(1)*gfact
1800 : :
1801 : : ! structure factors
1802 : 363528 : sfr=phre_vl3(ig1,ig2,ig3,iatom)
1803 : 363528 : sfi=-phimag_vl3(ig1,ig2,ig3,iatom)
1804 : :
1805 : : ! Multiply structure factor times vion derivatives:
1806 : 363528 : work1(re,ii)=sfr*(term1+term2)*two_pi
1807 : 363528 : work1(im,ii)=sfi*(term1+term2)*two_pi
1808 : :
1809 : : end if
1810 : :
1811 : : end do
1812 : : end if
1813 : : end do
1814 : : end do
1815 : :
1816 : : ! Transform back to real space
1817 : 486 : call fourdp(cplex,work1,vpsp1dqdq,1,mpi_enreg,nfft,1,ngfft,0)
1818 : :
1819 : 486 : xnorm=1.0_dp/ucvol
1820 [ + + ]: 3280986 : vpsp1dqdq(1:cplex*nfft)=vpsp1dqdq(1:cplex*nfft)*xnorm
1821 : :
1822 [ - + ]: 486 : ABI_FREE(work1)
1823 : :
1824 : : ! End the condition of non-electric-field
1825 : : end if
1826 : :
1827 : : contains
1828 : :
1829 : : !Real and imaginary parts of phase.
1830 : 363528 : function phr_vl3(x1,y1,x2,y2,x3,y3)
1831 : : real(dp) :: phr_vl3
1832 : : real(dp),intent(in) :: x1,x2,x3,y1,y2,y3
1833 : 363528 : phr_vl3=(x1*x2-y1*y2)*x3-(y1*x2+x1*y2)*y3
1834 : : end function phr_vl3
1835 : :
1836 : 363528 : function phi_vl3(x1,y1,x2,y2,x3,y3)
1837 : : real(dp) :: phi_vl3
1838 : : real(dp),intent(in) :: x1,x2,x3,y1,y2,y3
1839 : 363528 : phi_vl3=(x1*x2-y1*y2)*y3+(y1*x2+x1*y2)*x3
1840 : : end function phi_vl3
1841 : :
1842 : : ! Warning : this function differ from similar ones for ground-state calculations : note the atindx !!
1843 : 1454112 : function ph1_vl3(nri,ig1,ia)
1844 : : real(dp) :: ph1_vl3
1845 : : integer,intent(in) :: nri,ig1,ia
1846 : 1454112 : ph1_vl3=ph1d(nri,ig1+1+n1+(atindx(ia)-1)*(2*n1+1))
1847 : 1454112 : end function ph1_vl3
1848 : :
1849 : : ! Warning : this function differ from similar ones for ground-state calculations : note the atindx !!
1850 : 1454112 : function ph2_vl3(nri,ig2,ia)
1851 : : real(dp) :: ph2_vl3
1852 : : integer,intent(in) :: nri,ig2,ia
1853 : 1454112 : ph2_vl3=ph1d(nri,ig2+1+n2+(atindx(ia)-1)*(2*n2+1)+natom*(2*n1+1))
1854 : 1454112 : end function ph2_vl3
1855 : :
1856 : : ! Warning : this function differ from similar ones for ground-state calculations : note the atindx !!
1857 : 1454112 : function ph3_vl3(nri,ig3,ia)
1858 : : real(dp) :: ph3_vl3
1859 : : integer,intent(in) :: nri,ig3,ia
1860 : 1454112 : ph3_vl3=ph1d(nri,ig3+1+n3+(atindx(ia)-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1))
1861 : 1454112 : end function ph3_vl3
1862 : :
1863 : 363528 : function phre_vl3(ig1,ig2,ig3,ia)
1864 : : real(dp) :: phre_vl3
1865 : : integer,intent(in) :: ig1,ig2,ig3,ia
1866 : : phre_vl3=phr_vl3(ph1_vl3(re,ig1,ia),ph1_vl3(im,ig1,ia),&
1867 : 363528 : & ph2_vl3(re,ig2,ia),ph2_vl3(im,ig2,ia),ph3_vl3(re,ig3,ia),ph3_vl3(im,ig3,ia))
1868 : 363528 : end function phre_vl3
1869 : :
1870 : 363528 : function phimag_vl3(ig1,ig2,ig3,ia)
1871 : : real(dp) :: phimag_vl3
1872 : : integer,intent(in) :: ig1,ig2,ig3,ia
1873 : : phimag_vl3=phi_vl3(ph1_vl3(re,ig1,ia),ph1_vl3(im,ig1,ia),&
1874 : 363528 : & ph2_vl3(re,ig2,ia),ph2_vl3(im,ig2,ia),ph3_vl3(re,ig3,ia),ph3_vl3(im,ig3,ia))
1875 : 363528 : end function phimag_vl3
1876 : :
1877 : 1639764 : function gsq_vl3(g1,g2,g3)
1878 : : real(dp) :: gsq_vl3
1879 : : real(dp),intent(in) :: g1,g2,g3 ! Note that they are real, unlike in other similar function definitions
1880 : : !Define G^2 based on G space metric gmet.
1881 : : gsq_vl3=g1*g1*gmet(1,1)+g2*g2*gmet(2,2)+&
1882 : : & g3*g3*gmet(3,3)+2.0_dp*g1*g2*gmet(1,2)+&
1883 : 1639764 : & 2.0_dp*g2*g3*gmet(2,3)+2.0_dp*g3*g1*gmet(3,1)
1884 : 1639764 : end function gsq_vl3
1885 : :
1886 : : end subroutine dfpt_vlocaldqdq
1887 : : !!***
1888 : :
1889 : :
1890 : : !!****f* ABINIT/dfpt_vmetdqdq
1891 : : !! NAME
1892 : : !! dfpt_vmetdqdq
1893 : : !!
1894 : : !! FUNCTION
1895 : : !! Compute second q-gradient (at q=0) of the local part of 1st-order
1896 : : !! metric potential from the appropriate atomic pseudopotential
1897 : : !! with structure and derivative factor. Additionaly, compute the
1898 : : !! second q-gradient (at q=0) of the Hartree and XC (if GGA) potentials of the metric
1899 : : !! perturbation.
1900 : : !! Cartesian coordinates are employed to define the direction of the
1901 : : !! metric perturbation and the two q-gradients.
1902 : : !!
1903 : : !! INPUTS
1904 : : !! cplex: if 1, real space 1-order functions on FFT grid
1905 : : !! are REAL, if 2, COMPLEX
1906 : : !! gmet(3,3)=reciprocal space metric (Bohr**-2)
1907 : : !! gsqcut=cutoff G**2 for included G s in fft box.
1908 : : !! gprimd(3,3)=dimensional reciprocal space primitive translations
1909 : : !! idir= strain perturbation direction
1910 : : !! ipert=number of the atom being displaced in the frozen-phonon
1911 : : !! kxc(nfft,nkxc)=exchange and correlation kernel
1912 : : !! mpi_enreg=information about MPI parallelization
1913 : : !! mqgrid=dimension of q grid for pseudopotentials
1914 : : !! natom=number of atoms in cell.
1915 : : !! nattyp(ntypat)=number of atoms of each type in cell.
1916 : : !! nfft=(effective) number of FFT grid points (for this processor)
1917 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/input_variables/vargs.htm#ngfft
1918 : : !! nkxc=second dimension of the kxc array. If /=0, the XC kernel must be computed.
1919 : : !! nspden=number of spin-density components
1920 : : !! ntypat=number of types of atoms in cell.
1921 : : !! n1,n2,n3=fft grid.
1922 : : !! opthartdqdq= if 1 activates the calculation 2nd q-gradient of the Hartree potential
1923 : : !! ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phase information.
1924 : : !! qdir=direction of the q-gradient
1925 : : !! qgrid(mqgrid)=grid of q points from 0 to qmax.
1926 : : !! qphon(3)=wavevector of the phonon
1927 : : !! rhog(2,nfft)=array for Fourier transform of GS electron density
1928 : : !! rhor(nfftf,nspden)=array for GS electron density in electrons/bohr**3.
1929 : : !! ucvol=unit cell volume (Bohr**3).
1930 : : !! vlspl(mqgrid,2,ntypat)=spline fit of q^2 V(q) for each type of atom.
1931 : : !!
1932 : : !! OUTPUT
1933 : : !! vhart1dqdq(cplex*nfft)=2nd q-gradient (at q=0) of the GS density Hartree potential from the metric perturbation
1934 : : !! vpsp1dqdq(cplex*nfft)=2nd q-gradient (at q=0) of the first-order metric local
1935 : : !! crystal pseudopotential in real space
1936 : : !! vxc1dqdq(cplex*nfft)=2nd q-gradient (at q=0) of the GS density XC potential from the metric perturbation (only finite if GGA)
1937 : : !!
1938 : : !! NOTES
1939 : : !! ** IMPORTANT: the formalism followed in this routine
1940 : : !! assumes a phase factor for the perturbation that
1941 : : !! is different to the one used elsewhere in the code (See M.Stengel paper):
1942 : : !!
1943 : : !! here: e^{i q (R_l + \tau_{\kappa})}
1944 : : !! rest of ABINIT: e^{i q R_l}
1945 : : !!
1946 : : !! **Since the 2nd derivative w.r.t q-vector is calculated along cartesian
1947 : : !! directions, the 1/twopi**2 factor (that in the rest of the code is applied
1948 : : !! in the reduced to cartesian derivative conversion process) is here
1949 : : !! explicictly included in the formulas.
1950 : : !!
1951 : : !! **Notice that idir=1-9, in contrast to the strain perturbation (idir=1-6),
1952 : : !! because this term is not symmetric w.r.t permutations of the two strain
1953 : : !! indices.
1954 : : !!
1955 : : !! **A -i factor has been factorized out in all the contributions of the second
1956 : : !! q-gradient of the metric Hamiltonian. This is lately included in the contribution
1957 : : !! of the corresponing term (T4) to the flexoelectric tensor in dfpt_flexoout.F90
1958 : : !!
1959 : : !! SOURCE
1960 : :
1961 : 729 : subroutine dfpt_vmetdqdq(cplex,gmet,gprimd,gsqcut,idir,ipert,&
1962 : 729 : & kxc,mpi_enreg,mqgrid,natom,nattyp,nfft,ngfft,&
1963 : 729 : & ntypat,n1,n2,n3,nkxc,nspden,opthartdqdq,ph1d,qdir,qgrid,qphon,rhog,rhor,&
1964 : 729 : & ucvol,vlspl,vhart1dqdq,vpsp1dqdq,vxc1dqdq)
1965 : :
1966 : : !Arguments -------------------------------
1967 : : !scalars
1968 : : integer,intent(in) :: cplex,idir,ipert,mqgrid,n1,n2,n3,natom,nfft,nkxc,ntypat
1969 : : integer,intent(in) :: nspden,opthartdqdq,qdir
1970 : : real(dp),intent(in) :: gsqcut,ucvol
1971 : : type(MPI_type),intent(in) :: mpi_enreg
1972 : : !arrays
1973 : : integer,intent(in) :: nattyp(ntypat),ngfft(18)
1974 : : real(dp),intent(in) :: gmet(3,3),gprimd(3,3), kxc(nfft,nkxc)
1975 : : real(dp),intent(in) :: ph1d(2,(2*n1+1+2*n2+1+2*n3+1)*natom)
1976 : : real(dp),intent(in) :: qgrid(mqgrid),qphon(3),rhog(2,nfft),rhor(cplex*nfft,nspden)
1977 : : real(dp),intent(in) :: vlspl(mqgrid,2,ntypat)
1978 : : real(dp),intent(out) :: vhart1dqdq(cplex*nfft),vpsp1dqdq(cplex*nfft)
1979 : : real(dp),intent(out) :: vxc1dqdq(cplex*nfft)
1980 : :
1981 : : !Local variables -------------------------
1982 : : !scalars
1983 : : integer :: beta, delta, gamma
1984 : : integer :: ia,i1,i2,i3,ia1,ia2,id1,id2,id3,ig1,ig2,ig3,ii,ii1
1985 : : integer :: itypat,jj
1986 : : integer, parameter :: im=2, re=1
1987 : : real(dp),parameter :: tolfix=1.000000001_dp
1988 : : real(dp) :: cutoff,delbd,delbg,deldg,gfact,gmag,gq1
1989 : : real(dp) :: gq2,gq3,gsquar,pisqrinv
1990 : : real(dp) :: sfi,sfr,term1,term2,tmpre,tmpim,uogsquar,work1re,xnorm
1991 : : logical :: qeq0
1992 : : character(len=500) :: msg
1993 : : !arrays
1994 : : integer,save :: idx(18)=(/1,1,2,2,3,3,3,2,3,1,2,1,2,3,1,3,1,2/)
1995 : 729 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
1996 : 729 : integer, contiguous, pointer :: fftn3_distrib(:),ffti3_local(:)
1997 : : real(dp) :: gq(3),gqc(3),vion1(1),vion1dq(1),vion1dqdq(1)
1998 : 729 : real(dp),allocatable :: work1(:,:)
1999 : : ! *********************************************************************
2000 : :
2001 [ + + - + ]: 729 : if(ipert/=natom+3 .and. ipert/=natom+4)then
2002 : :
2003 [ # # ]: 0 : vpsp1dqdq(1:cplex*nfft)=zero
2004 : :
2005 : : else
2006 : :
2007 : 729 : beta=idx(2*idir-1); delta=idx(2*idir); gamma=qdir
2008 : :
2009 : : !Kronecker deltas
2010 : 729 : delbd=0.0_dp; delbg=0.0_dp; deldg=0.0_dp
2011 [ + + ]: 729 : if (beta==delta) delbd=1.0_dp
2012 [ + + ]: 729 : if (beta==gamma) delbg=1.0_dp
2013 [ + + ]: 729 : if (delta==gamma) deldg=1.0_dp
2014 : :
2015 [ - + + - : 2187 : ABI_MALLOC(work1,(2,nfft))
+ - + - -
+ - + ]
2016 [ + + + + ]: 7381854 : work1(1:2,1:nfft)=0.0_dp
2017 : :
2018 : 729 : cutoff=gsqcut*tolfix
2019 : 729 : id1=n1/2+2
2020 : 729 : id2=n2/2+2
2021 : 729 : id3=n3/2+2
2022 : :
2023 : : !Get the distrib associated with this fft_grid
2024 : 729 : call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
2025 : :
2026 : : ! This is to allow q=0
2027 : 729 : qeq0=.false.
2028 [ - + ]: 729 : if(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15) then
2029 : : qeq0=.true.
2030 : : else
2031 : 0 : msg='This routine cannot be used for q/=0'
2032 : 0 : ABI_BUG(msg)
2033 : : end if
2034 : :
2035 : 729 : ia1=1
2036 [ + + ]: 1458 : do itypat=1,ntypat
2037 : : ! ia1,ia2 sets range of loop over atoms:
2038 : 729 : ia2=ia1+nattyp(itypat)-1
2039 : :
2040 : 729 : ii=0
2041 : :
2042 [ + + ]: 11664 : do i3=1,n3
2043 : 10935 : ig3=i3-(i3/id3)*n3-1
2044 : 10935 : gq3=dble(ig3)+qphon(3)
2045 : 10935 : gq(3)=gq3
2046 [ + + ]: 175689 : do i2=1,n2
2047 [ + - ]: 174960 : if (fftn2_distrib(i2)==mpi_enreg%me_fft) then
2048 : 164025 : ig2=i2-(i2/id2)*n2-1
2049 : 164025 : gq2=dble(ig2)+qphon(2)
2050 : 164025 : gq(2)=gq2
2051 : :
2052 : : ! Note the lower limit of the next loop
2053 : 164025 : ii1=1
2054 [ + + + - : 164025 : if(i3==1 .and. i2==1 .and. qeq0 .and. ig2==0 .and. ig3==0)then
+ - ]
2055 : 729 : ii1=2
2056 : 729 : ii=ii+1
2057 : : end if
2058 [ + + ]: 2623671 : do i1=ii1,n1
2059 : 2459646 : ig1=i1-(i1/id1)*n1-1
2060 : 2459646 : gq1=dble(ig1)+qphon(1)
2061 : 2459646 : gq(1)=gq1
2062 : 2459646 : ii=ii+1
2063 : 2459646 : gsquar=gsq_vl(ig1,ig2,ig3)
2064 : : ! Skip G**2 outside cutoff:
2065 [ + + ]: 2623671 : if (gsquar<=cutoff) then
2066 : 545292 : gmag=sqrt(gsquar)
2067 : :
2068 : : ! Obtain G in cartesian coordinates
2069 : 545292 : gqc(1)=gprimd(1,1)*gq(1)+gprimd(1,2)*gq(2)+gprimd(1,3)*gq(3)
2070 : 545292 : gqc(2)=gprimd(2,1)*gq(1)+gprimd(2,2)*gq(2)+gprimd(2,3)*gq(3)
2071 : 545292 : gqc(3)=gprimd(3,1)*gq(1)+gprimd(3,2)*gq(2)+gprimd(3,3)*gq(3)
2072 : :
2073 : : ! Evaluate spline fit to get first V(q) and V(q)' and later V(q)'':
2074 [ + + ]: 1090584 : call splfit(qgrid,vion1dq,vlspl(:,:,itypat),1,(/gmag/),vion1,mqgrid,1)
2075 [ + + ]: 1090584 : vion1=vion1/gsquar
2076 [ + + ]: 1090584 : vion1dq=(vion1dq-2.0_dp*gmag*vion1)/gsquar
2077 : :
2078 [ + + ]: 1090584 : call splfit(qgrid,vion1dqdq,vlspl(:,:,itypat),2,(/gmag/),vion1,mqgrid,1)
2079 [ + + ]: 1090584 : vion1dqdq=(vion1dqdq-4.0_dp*gmag*vion1dq-2.0_dp*vion1)/gsquar
2080 : :
2081 : : ! Assemble structure factor over all atoms of given type:
2082 : 545292 : sfr=0.0_dp
2083 : 545292 : sfi=0.0_dp
2084 [ + + ]: 1635876 : do ia=ia1,ia2
2085 : 1090584 : sfr=sfr+phre_vl(ig1,ig2,ig3,ia)
2086 : 1635876 : sfi=sfi-phimag_vl(ig1,ig2,ig3,ia)
2087 : : end do
2088 : :
2089 : 545292 : gfact=gqc(beta)*gqc(delta)*gqc(gamma)/gsquar
2090 : :
2091 : 545292 : term1=delbd*gqc(gamma)+delbg*gqc(delta)+deldg*gqc(beta)
2092 : 545292 : term1=term1-gfact
2093 : 545292 : term1=term1*vion1dq(1)/gmag
2094 : :
2095 : 545292 : term2=vion1dqdq(1)*gfact
2096 : :
2097 : : ! Multiply structure factor times vion derivatives:
2098 : 545292 : work1(re,ii)=work1(re,ii)+sfr*(term1+term2)
2099 : 545292 : work1(im,ii)=work1(im,ii)+sfi*(term1+term2)
2100 : :
2101 : : ! End skip G**2 outside cutoff:
2102 : : end if
2103 : :
2104 : : end do
2105 : : end if
2106 : : end do
2107 : : end do
2108 : :
2109 : 1458 : ia1=ia2+1
2110 : :
2111 : : ! End loop on type of atoms
2112 : : end do
2113 : :
2114 : : ! Set Vloc(G=0)=0:
2115 : 729 : work1(re,1)=0.0_dp
2116 : 729 : work1(im,1)=0.0_dp
2117 : :
2118 : : ! Transform back to real space
2119 : 729 : call fourdp(cplex,work1,vpsp1dqdq,1,mpi_enreg,nfft,1,ngfft,0)
2120 : :
2121 : 729 : xnorm=1.0_dp/ucvol/two_pi
2122 [ + + ]: 4921479 : vpsp1dqdq(1:cplex*nfft)=vpsp1dqdq(1:cplex*nfft)*xnorm
2123 : :
2124 [ + + + + ]: 7381854 : work1=0.0_dp
2125 : :
2126 : : ! Calculate the GS density Hartree contribution
2127 [ + - ]: 729 : if (opthartdqdq==1) then
2128 : :
2129 : : pisqrinv=1.0_dp/pi**2
2130 : :
2131 : : ii=0
2132 [ + + ]: 11664 : do i3=1,n3
2133 : 10935 : ig3=i3-(i3/id3)*n3-1
2134 : 10935 : gq3=dble(ig3)+qphon(3)
2135 : 10935 : gq(3)=gq3
2136 [ + + ]: 175689 : do i2=1,n2
2137 [ + - ]: 174960 : if (fftn2_distrib(i2)==mpi_enreg%me_fft) then
2138 : 164025 : ig2=i2-(i2/id2)*n2-1
2139 : 164025 : gq2=dble(ig2)+qphon(2)
2140 : 164025 : gq(2)=gq2
2141 : :
2142 : : ! Note the lower limit of the next loop
2143 : 164025 : ii1=1
2144 [ + + + - : 164025 : if(i3==1 .and. i2==1 .and. qeq0 .and. ig2==0 .and. ig3==0)then
+ - + - ]
2145 : 729 : ii1=2
2146 : 729 : ii=ii+1
2147 : : end if
2148 [ + + ]: 2623671 : do i1=ii1,n1
2149 : 2459646 : ig1=i1-(i1/id1)*n1-1
2150 : 2459646 : gq1=dble(ig1)+qphon(1)
2151 : 2459646 : gq(1)=gq1
2152 : 2459646 : ii=ii+1
2153 : 2459646 : gsquar=gsq_vl(ig1,ig2,ig3)
2154 : : ! Skip G**2 outside cutoff:
2155 [ + + ]: 2623671 : if (gsquar<=cutoff) then
2156 : :
2157 : : ! Precalculate quotient of G powers
2158 : 545292 : uogsquar= 1.0_dp/gsquar
2159 : :
2160 : : ! Obtain G in cartesian coordinates
2161 : 545292 : gqc(1)=gprimd(1,1)*gq(1)+gprimd(1,2)*gq(2)+gprimd(1,3)*gq(3)
2162 : 545292 : gqc(2)=gprimd(2,1)*gq(1)+gprimd(2,2)*gq(2)+gprimd(2,3)*gq(3)
2163 : 545292 : gqc(3)=gprimd(3,1)*gq(1)+gprimd(3,2)*gq(2)+gprimd(3,3)*gq(3)
2164 : :
2165 : 545292 : term1=4.0_dp*gqc(beta)*gqc(gamma)*gqc(delta)*uogsquar*uogsquar
2166 : 545292 : term2=delbd*gqc(gamma)+delbg*gqc(delta)+deldg*gqc(beta)
2167 : 545292 : term2=-term2*uogsquar
2168 : :
2169 : 545292 : work1re=pisqrinv*uogsquar*(term1+term2)
2170 : 545292 : work1(re,ii)=rhog(re,ii)*work1re
2171 : 545292 : work1(im,ii)=rhog(im,ii)*work1re
2172 : :
2173 : : ! End skip G**2 outside cutoff:
2174 : : end if
2175 : :
2176 : : end do
2177 : : end if
2178 : : end do
2179 : : end do
2180 : :
2181 : : ! Set V(G=0)=0:
2182 : 729 : work1(re,1)=0.0_dp
2183 : 729 : work1(im,1)=0.0_dp
2184 : :
2185 : : ! Transform back to real space
2186 : 729 : call fourdp(cplex,work1,vhart1dqdq,1,mpi_enreg,nfft,1,ngfft,0)
2187 : :
2188 : : ! End the calculation of the Hartree contribution
2189 : : end if
2190 : :
2191 : 729 : ABI_FREE(work1)
2192 : :
2193 : : ! Calculate the GS density XC contribution (if GGA)
2194 [ + + ]: 4921479 : vxc1dqdq(:)=zero
2195 [ + + ]: 729 : if (nkxc == 7) then
2196 : : call dfpt_mkvxcgga_n0met(beta,1,delta,gamma,gprimd,kxc,mpi_enreg, &
2197 : 243 : & nfft,ngfft,nkxc,nspden,rhor,vxc1dqdq)
2198 : :
2199 : : !Fictitious i factor temporarily applied.
2200 : : !It is later canceled by the (-i) factor of the total matrix element
2201 [ + + ]: 820368 : do ii=1,nfft
2202 : 820125 : jj=ii*2
2203 : 820125 : tmpre=vxc1dqdq(jj-1); tmpim=vxc1dqdq(jj)
2204 : 820368 : vxc1dqdq(jj-1)=-tmpim; vxc1dqdq(jj)=tmpre
2205 : : end do
2206 : : end if
2207 : :
2208 : : !End the condition of non-electric-field
2209 : : end if
2210 : :
2211 : : contains
2212 : :
2213 : : !Real and imaginary parts of phase.
2214 : 1090584 : function phr_vl(x1,y1,x2,y2,x3,y3)
2215 : : real(dp) :: phr_vl,x1,x2,x3,y1,y2,y3
2216 : 1090584 : phr_vl=(x1*x2-y1*y2)*x3-(y1*x2+x1*y2)*y3
2217 : : end function phr_vl
2218 : :
2219 : 1090584 : function phi_vl(x1,y1,x2,y2,x3,y3)
2220 : : real(dp):: phi_vl,x1,x2,x3,y1,y2,y3
2221 : 1090584 : phi_vl=(x1*x2-y1*y2)*y3+(y1*x2+x1*y2)*x3
2222 : : end function phi_vl
2223 : :
2224 : 4362336 : function ph1_vl(nri,ig1,ia)
2225 : : real(dp):: ph1_vl
2226 : : integer :: nri,ig1,ia
2227 : 4362336 : ph1_vl=ph1d(nri,ig1+1+n1+(ia-1)*(2*n1+1))
2228 : 4362336 : end function ph1_vl
2229 : :
2230 : 4362336 : function ph2_vl(nri,ig2,ia)
2231 : : real(dp):: ph2_vl
2232 : : integer :: nri,ig2,ia
2233 : 4362336 : ph2_vl=ph1d(nri,ig2+1+n2+(ia-1)*(2*n2+1)+natom*(2*n1+1))
2234 : 4362336 : end function ph2_vl
2235 : :
2236 : 4362336 : function ph3_vl(nri,ig3,ia)
2237 : : real(dp):: ph3_vl
2238 : : integer :: nri,ig3,ia
2239 : 4362336 : ph3_vl=ph1d(nri,ig3+1+n3+(ia-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1))
2240 : 4362336 : end function ph3_vl
2241 : :
2242 : 1090584 : function phre_vl(ig1,ig2,ig3,ia)
2243 : : real(dp):: phre_vl
2244 : : integer :: ig1,ig2,ig3,ia
2245 : : phre_vl=phr_vl(ph1_vl(re,ig1,ia),ph1_vl(im,ig1,ia),&
2246 : 1090584 : & ph2_vl(re,ig2,ia),ph2_vl(im,ig2,ia),ph3_vl(re,ig3,ia),ph3_vl(im,ig3,ia))
2247 : 1090584 : end function phre_vl
2248 : :
2249 : 1090584 : function phimag_vl(ig1,ig2,ig3,ia)
2250 : : real(dp) :: phimag_vl
2251 : : integer :: ig1,ig2,ig3,ia
2252 : : phimag_vl=phi_vl(ph1_vl(re,ig1,ia),ph1_vl(im,ig1,ia),&
2253 : 1090584 : & ph2_vl(re,ig2,ia),ph2_vl(im,ig2,ia),ph3_vl(re,ig3,ia),ph3_vl(im,ig3,ia))
2254 : 1090584 : end function phimag_vl
2255 : :
2256 : 4919292 : function gsq_vl(i1,i2,i3)
2257 : : real(dp) :: gsq_vl
2258 : : integer :: i1,i2,i3
2259 : : !Define G^2 based on G space metric gmet.
2260 : : gsq_vl=dble(i1*i1)*gmet(1,1)+dble(i2*i2)*gmet(2,2)+&
2261 : : & dble(i3*i3)*gmet(3,3)+dble(2*i1*i2)*gmet(1,2)+&
2262 : 4919292 : & dble(2*i2*i3)*gmet(2,3)+dble(2*i3*i1)*gmet(3,1)
2263 : 4919292 : end function gsq_vl
2264 : :
2265 : : end subroutine dfpt_vmetdqdq
2266 : :
2267 : : end module m_mklocl
2268 : : !!***
|
