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 23934 : subroutine mklocl(dtset, dyfrlo,eei,gmet,gprimd,grtn,gsqcut,lpsstr,mgfft,&
131 23934 : & mpi_enreg,natom,nattyp,nfft,ngfft,nspden,ntypat,option,pawtab,ph1d,psps,qprtrb,&
132 23934 : & 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 23934 : real(dp),allocatable :: xcart(:,:)
158 : #if defined HAVE_BIGDFT
159 : real(dp),pointer :: rhor_ptr(:,:)
160 : #endif
161 :
162 : ! *************************************************************************
163 :
164 23934 : 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 23934 : 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 23934 : 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 23934 : if (dtset%usewvl == 0) then
185 : ! Plane wave case
186 23934 : 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 23934 : & 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 23934 : 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 29246 : subroutine mklocl_recipspace(dyfrlo,eei,gmet,gprimd,grtn,gsqcut,icutcoul,lpsstr,mgfft,&
274 29246 : & mpi_enreg,mqgrid,natom,nattyp,nfft,ngfft,nkpt,ntypat,option,ph1d,qgrid,qprtrb,&
275 29246 : & 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 29246 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
309 29246 : integer, contiguous, pointer :: fftn3_distrib(:),ffti3_local(:)
310 : real(dp) :: gcart(3),tsec(2)
311 29246 : real(dp),allocatable :: gcutoff(:)
312 29246 : 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 29246 : if(option==2)then
337 17073 : call timab(72,1,tsec)
338 : end if
339 29246 : n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
340 29246 : me_fft=ngfft(11)
341 29246 : nproc_fft=ngfft(10)
342 :
343 : !Get the distrib associated with this fft_grid
344 29246 : 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 29246 : if(option==1)then
348 20583 : ABI_MALLOC(work1,(2,nfft))
349 341143209 : work1(:,:)=zero
350 : end if
351 :
352 : !
353 29246 : dq=(qgrid(mqgrid)-qgrid(1))/dble(mqgrid-1)
354 29246 : dqm1=1.0_dp/dq
355 29246 : dqdiv6=dq/6.0_dp
356 29246 : dq2div6=dq**2/6.0_dp
357 29246 : cutoff=gsqcut*tolfix
358 29246 : id1=n1/2+2
359 29246 : id2=n2/2+2
360 29246 : id3=n3/2+2
361 349194 : grtn(:,:)=zero
362 29246 : lpsstr(:)=zero
363 1069077 : dyfrlo(:,:,:)=zero
364 29246 : me_g0=0
365 29246 : ia1=1
366 :
367 : !Initialize Gcut-off array from m_gtermcutoff
368 29246 : call termcutoff(gcutoff,gsqcut,icutcoul,ngfft,nkpt,rcut,rprimd,vcutgeo)
369 29246 : rcut_loc = half*SQRT(DOT_PRODUCT(rprimd(:,3),rprimd(:,3)))
370 :
371 72544 : do itypat=1,ntypat
372 : ! ia1,ia2 sets range of loop over atoms:
373 43298 : ia2=ia1+nattyp(itypat)-1
374 :
375 43298 : ii=0
376 831004 : do i3=1,n3
377 787706 : ig3=i3-(i3/id3)*n3-1
378 19051151 : do i2=1,n2
379 18220147 : ig2=i2-(i2/id2)*n2-1
380 19007853 : if(fftn2_distrib(i2) == me_fft ) then
381 636226579 : do i1=1,n1
382 618410760 : ig1=i1-(i1/id1)*n1-1
383 :
384 618410760 : 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 618410760 : if (ig1==0 .and. ig2==0 .and. ig3==0) me_g0=1
389 618410760 : if (ig1==0 .and. ig2==0 .and. ig3==0) cycle
390 :
391 618368130 : gsquar=gsq_mk(ig1,ig2,ig3)
392 : ! Skip G**2 outside cutoff:
393 636183949 : if (gsquar<=cutoff) then
394 258249564 : gmag=sqrt(gsquar)
395 :
396 : ! Compute vion(G) for given type of atom
397 258249564 : jj=1+int(gmag*dqm1)
398 258249564 : 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 258249564 : bb = diff*dqm1
406 258249564 : aa = 1.0_dp-bb
407 258249564 : cc = aa*(aa**2-1.0_dp)*dq2div6
408 258249564 : dd = bb*(bb**2-1.0_dp)*dq2div6
409 :
410 : vion1 = (aa*vlspl(jj,1,itypat)+bb*vlspl(jj+1,1,itypat) +&
411 258249564 : & cc*vlspl(jj,2,itypat)+dd*vlspl(jj+1,2,itypat) ) / gsquar * gcutoff(ii)
412 :
413 258249564 : if(option==1)then
414 :
415 : ! Assemble structure factor over all atoms of given type:
416 58652841 : sfr=zero
417 58652841 : sfi=zero
418 193614974 : do ia=ia1,ia2
419 134962133 : sfr=sfr+phre_mk(ig1,ig2,ig3,ia)
420 193614974 : sfi=sfi-phimag_mk(ig1,ig2,ig3,ia)
421 : end do
422 : ! Multiply structure factor times vion:
423 58652841 : work1(re,ii)=work1(re,ii)+sfr*vion1
424 58652841 : work1(im,ii)=work1(im,ii)+sfi*vion1
425 :
426 199596723 : else if(option==2 .or. option==4)then
427 :
428 : ! Compute Re and Im part of (2Pi)*Vion(G)*rho(G):
429 152077859 : svionr=(two_pi*vion1)*rhog(re,ii)
430 152077859 : svioni=(two_pi*vion1)*rhog(im,ii)
431 :
432 : ! Loop over all atoms of this type:
433 433820355 : 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 281742496 : shift1=1+n1+(ia-1)*(2*n1+1)
440 281742496 : shift2=1+n2+(ia-1)*(2*n2+1)+natom*(2*n1+1)
441 281742496 : shift3=1+n3+(ia-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1)
442 281742496 : ph1r=ph1d(1,ig1+shift1)
443 281742496 : ph1i=ph1d(2,ig1+shift1)
444 281742496 : ph2r=ph1d(1,ig2+shift2)
445 281742496 : ph2i=ph1d(2,ig2+shift2)
446 281742496 : ph3r=ph1d(1,ig3+shift3)
447 281742496 : ph3i=ph1d(2,ig3+shift3)
448 281742496 : ph12r=ph1r*ph2r-ph1i*ph2i
449 281742496 : ph12i=ph1r*ph2i+ph1i*ph2r
450 281742496 : phre_igia=ph12r*ph3r-ph12i*ph3i
451 281742496 : phimag_igia=ph12r*ph3i+ph12i*ph3r
452 :
453 433820355 : 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 277926218 : svion=svioni*phre_igia+svionr*phimag_igia
458 :
459 : ! Open loop over 3-index for speed:
460 277926218 : grtn(1,ia)=grtn(1,ia)-dble(ig1)*svion
461 277926218 : grtn(2,ia)=grtn(2,ia)-dble(ig2)*svion
462 277926218 : 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 47518864 : 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 47518864 : ee= vlspl(jj+1,1,itypat)-vlspl(jj,1,itypat)
491 : ff= (3._dp*bb**2-1._dp)*vlspl(jj+1,2,itypat) &
492 47518864 : & - (3._dp*aa**2-1._dp)*vlspl(jj,2,itypat)
493 : vion2 = ( ( ee*dqm1 + ff*dqdiv6 )/gmag&
494 47518864 : & - 2.0_dp*vion1 ) / gsquar
495 :
496 : gcart(1)=gprimd(1,1)*dble(ig1)+gprimd(1,2)*dble(ig2)+&
497 47518864 : & gprimd(1,3)*dble(ig3)
498 : gcart(2)=gprimd(2,1)*dble(ig1)+gprimd(2,2)*dble(ig2)+&
499 47518864 : & gprimd(2,3)*dble(ig3)
500 : gcart(3)=gprimd(3,1)*dble(ig1)+gprimd(3,2)*dble(ig2)+&
501 47518864 : & gprimd(3,3)*dble(ig3)
502 : ! Assemble structure over all atoms of given type
503 47518864 : sfr=zero
504 47518864 : sfi=zero
505 156467745 : do ia=ia1,ia2
506 108948881 : sfr=sfr+phre_mk(ig1,ig2,ig3,ia)
507 156467745 : sfi=sfi-phimag_mk(ig1,ig2,ig3,ia)
508 : end do
509 : ! Compute Re( rho^*(G)* sf ) * [(dV(G)/dG)/|G|]
510 47518864 : term=(rhog(re,ii)*sfr+rhog(im,ii)*sfi)*vion2
511 : ! Compute contribution to stress tensor
512 47518864 : lpsstr(1)=lpsstr(1)-term*gcart(1)*gcart(1)
513 47518864 : lpsstr(2)=lpsstr(2)-term*gcart(2)*gcart(2)
514 47518864 : lpsstr(3)=lpsstr(3)-term*gcart(3)*gcart(3)
515 47518864 : lpsstr(4)=lpsstr(4)-term*gcart(3)*gcart(2)
516 47518864 : lpsstr(5)=lpsstr(5)-term*gcart(3)*gcart(1)
517 47518864 : 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 123285 : do ia=ia1,ia2
580 79987 : dyfrlo(2,1,ia)=dyfrlo(1,2,ia)
581 79987 : dyfrlo(3,1,ia)=dyfrlo(1,3,ia)
582 123285 : dyfrlo(3,2,ia)=dyfrlo(2,3,ia)
583 : end do
584 :
585 72544 : ia1=ia2+1
586 :
587 : ! End loop on type of atoms
588 : end do
589 :
590 29246 : 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 6861 : if(me_g0 == 1) then
594 : ! Set Vloc(G=0)=0:
595 6757 : work1(re,1)=zero
596 6757 : 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 6861 : 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 6861 : call fourdp(1,work1,vpsp,1,mpi_enreg,nfft,1,ngfft,0)
629 :
630 : ! Divide by unit cell volume
631 6861 : xnorm=1.0_dp/ucvol
632 113718977 : vpsp(:)=vpsp(:)*xnorm
633 :
634 6861 : ABI_FREE(work1)
635 :
636 : end if
637 :
638 29246 : ABI_FREE(gcutoff)
639 :
640 29246 : if(option==2)then
641 : ! Init mpi_comm
642 17073 : 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 17073 : call timab(72,2,tsec)
648 : end if
649 :
650 29246 : if(option==3)then
651 : ! Init mpi_comm
652 4907 : 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 4907 : lpsstr(1)=(lpsstr(1)-eei)/ucvol
662 4907 : lpsstr(2)=(lpsstr(2)-eei)/ucvol
663 4907 : lpsstr(3)=(lpsstr(3)-eei)/ucvol
664 4907 : lpsstr(4)=lpsstr(4)/ucvol
665 4907 : lpsstr(5)=lpsstr(5)/ucvol
666 4907 : 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 58897 : 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 243911014 : function phr_mk(x1,y1,x2,y2,x3,y3)
692 :
693 : real(dp) :: phr_mk,x1,x2,x3,y1,y2,y3
694 243911014 : phr_mk=(x1*x2-y1*y2)*x3-(y1*x2+x1*y2)*y3
695 : end function phr_mk
696 :
697 243911014 : function phi_mk(x1,y1,x2,y2,x3,y3)
698 :
699 : real(dp):: phi_mk,x1,x2,x3,y1,y2,y3
700 243911014 : phi_mk=(x1*x2-y1*y2)*y3+(y1*x2+x1*y2)*x3
701 : end function phi_mk
702 :
703 975644056 : function ph1_mk(nri,ig1,ia)
704 :
705 : real(dp):: ph1_mk
706 : integer :: nri,ig1,ia
707 975644056 : ph1_mk=ph1d(nri,ig1+1+n1+(ia-1)*(2*n1+1))
708 975644056 : end function ph1_mk
709 :
710 975644056 : function ph2_mk(nri,ig2,ia)
711 :
712 : real(dp):: ph2_mk
713 : integer :: nri,ig2,ia
714 975644056 : ph2_mk=ph1d(nri,ig2+1+n2+(ia-1)*(2*n2+1)+natom*(2*n1+1))
715 975644056 : end function ph2_mk
716 :
717 975644056 : function ph3_mk(nri,ig3,ia)
718 :
719 : real(dp):: ph3_mk
720 : integer :: nri,ig3,ia
721 975644056 : ph3_mk=ph1d(nri,ig3+1+n3+(ia-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1))
722 975644056 : end function ph3_mk
723 :
724 243911014 : 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 243911014 : & ph2_mk(re,ig2,ia),ph2_mk(im,ig2,ia),ph3_mk(re,ig3,ia),ph3_mk(im,ig3,ia))
730 243911014 : end function phre_mk
731 :
732 243911014 : 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 243911014 : & ph2_mk(re,ig2,ia),ph2_mk(im,ig2,ia),ph3_mk(re,ig3,ia),ph3_mk(im,ig3,ia))
738 243911014 : end function phimag_mk
739 :
740 618368130 : 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 618368130 : & dble(2*i2*i3)*gmet(2,3)+dble(2*i3*i1)*gmet(3,1)
747 618368130 : 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 : !!***
|
