Branch data Line data Source code
1 : : !!****m* ABINIT/m_atm2fft
2 : : !! NAME
3 : : !! m_atm2fft
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !! COPYRIGHT
8 : : !! Copyright (C) 1998-2026 ABINIT group (FJ, MT)
9 : : !! This file is distributed under the terms of the
10 : : !! GNU General Public License, see ~abinit/COPYING
11 : : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : : !!
13 : : !! SOURCE
14 : :
15 : : #if defined HAVE_CONFIG_H
16 : : #include "config.h"
17 : : #endif
18 : :
19 : : #include "abi_common.h"
20 : :
21 : : ! nvtx related macro definition
22 : : #include "nvtx_macros.h"
23 : :
24 : : module m_atm2fft
25 : :
26 : : use defs_basis
27 : : use m_abicore
28 : : use m_errors
29 : : use m_xmpi
30 : : use m_dtset
31 : :
32 : : use defs_abitypes, only : mpi_type
33 : : use m_time, only : timab
34 : : use defs_datatypes,only : pseudopotential_type
35 : : use m_distribfft, only : distribfft_type
36 : : use m_gtermcutoff, only : termcutoff
37 : : use m_pawtab, only : pawtab_type
38 : : use m_fft, only : zerosym, fourdp
39 : : use m_mpinfo, only : set_mpi_enreg_fft, unset_mpi_enreg_fft, initmpi_seq
40 : :
41 : : #if defined(HAVE_GPU_MARKERS)
42 : : use m_nvtx_data
43 : : #endif
44 : :
45 : : implicit none
46 : :
47 : : private
48 : : !!***
49 : :
50 : : public :: atm2fft
51 : : public :: dfpt_atm2fft
52 : : !!***
53 : :
54 : : contains
55 : : !!***
56 : :
57 : : !!****f* ABINIT/atm2fft
58 : : !! NAME
59 : : !! atm2fft
60 : : !!
61 : : !! FUNCTION
62 : : !! This routine sums atomic functions (density, kinetic density or potential) defined
63 : : !! (in rec. space) on a radial grid to get global quantities on the
64 : : !! fine FFT grid. It can also compute contribution to energy derivatives
65 : : !! of these atomic functions.
66 : : !!
67 : : !! Possible options:
68 : : !! optn=1: compute a sum of local atomic [kinetic] densities or contrib. to energy derivatives
69 : : !! optv=1: compute a sum of local atomic potentials or contrib. to energy derivatives
70 : : !!
71 : : !! optatm =1: computes sum of atomic potentials/densities
72 : : !! optgr =1: computes contribution of atomic pot./dens. to forces
73 : : !! optstr =1: computes contribution of atomic pot./dens. to stress tensor
74 : : !! optdyfr=1: computes contribution of atomic pot./dens. to frozen part of dyn. matrix
75 : : !!
76 : : !! INPUTS
77 : : !! atindx1(natom)=index table for atoms, inverse of atindx
78 : : !! distribfft<type(distribfft_type)>=--optional-- contains infos related to FFT parallelism
79 : : !! gauss(2,ntypat)= params for gaussian atm density (optn2=3) for each atom type
80 : : !! gmet(3,3)=reciprocal space metric
81 : : !! gprimd(3,3)=reciprocal space dimensional primitive translations
82 : : !! gsqcut=cutoff on |G|^2: see setup1 for definition (doubled sphere)
83 : : !! kxc(2,nfft)=exchange and correlation kernel
84 : : !! me_g0=--optional-- 1 if the current process treat the g=0 plane-wave (only needed when comm_fft is present)
85 : : !! mgfft=maximum size of 1D FFTs
86 : : !! comm_fft=--optional-- MPI communicator over FFT components
87 : : !! mqgrid=number of grid pts in q array for f(q) spline.
88 : : !! natom=number of atoms in unit cell.
89 : : !! nattyp(ntypat)=number of atoms of each type in cell.
90 : : !! nfft=(effective) number of FFT grid points (for this processor)
91 : : !! ngfft(18)=contain all needed information about 3D FFT
92 : : !! ntypat=number of types of atoms.
93 : : !! optatm,optdyfr,optgr,optn,optn2,optstr,optv= (see NOTES below)
94 : : !! paral_kgb=--optional-- 1 if "band-FFT" parallelism is activated (only needed when comm_fft is present)
95 : : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
96 : : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
97 : : !! ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phase information
98 : : !! qgrid(mqgrid)=q grid for spline from 0 to qmax.
99 : : !! qprtrb(3)= integer wavevector of possible perturbing potential
100 : : !! in basis of reciprocal lattice translations
101 : : !! rhog(2,nfft)=electron density rho(G) in reciprocal space
102 : : !! (used only if optv=1 and (optgr=1 or optstr=1 or optdyfr=1))
103 : : !! ucvol=unit cell volume
104 : : !! usepaw= 0 for non paw calculation; =1 for paw calculation
105 : : !! vspl(mqgrid,2,ntypat)=q^2 v(q) spline of an atomic potential
106 : : !! (used only if optv=1)
107 : : !! vprtrb(2)=complex amplitude of possible perturbing potential; if nonzero,
108 : : !! perturbing potential is added of the form V(G)=(vprtrb(1)+I*vprtrb(2))/2
109 : : !! at the values G=qprtrb and (vprtrb(1)-I*vprtrb(2))/2 at G=-qprtrb
110 : : !! vg(2,nfft)= potential V(G) in reciprocal space
111 : : !! (used only if optn=1 and (optgr=1 or optstr=1 or optdyfr=1))
112 : : !! vg1(2,nfft)= 1st-order potential V(G) in reciprocal space
113 : : !! (used only if opteltfr==1)
114 : : !! vg1_core(2,nfft)= 1-order potential V(G) in reciprocal space with only core contribution
115 : : !! (used only if opteltfr==1)
116 : : !! OUTPUT
117 : : !! ======= if optv==1 =======
118 : : !! ============================
119 : : !! --- if optatm==1
120 : : !! atmvloc(nfft)=sum of local atomic potentials in real space
121 : : !! --- if optgr==1
122 : : !! grv(3,natom)=contribution of atomic potentials to forces
123 : : !! --- if optstr==1
124 : : !! strv(6)=contribution of atomic potentials to stress tensor
125 : : !! cart. coordinates, symmetric tensor, 6 comp. in order 11,22,33,32,31,21
126 : : !! --- if optdyfr==1
127 : : !! dyfrv(3,3,natom)=contribution of atomic potentials to frozen part of dyn. matrix
128 : : !!
129 : : !! ======= if optn==1 =======
130 : : !! ============================
131 : : !! --- if optatm==1
132 : : !! atmrho(nfft)=sum of atomic densities in real space
133 : : !! --- if optgr==1
134 : : !! grn(3,natom)=contribution of atomic densities to forces
135 : : !! --- if optstr==1
136 : : !! strn(6)=contribution of atomic densities to stress tensor
137 : : !! cart. coordinates, symmetric tensor, 6 comp. in order 11,22,33,32,31,21
138 : : !! --- if optdyfr==1
139 : : !! dyfrn(3,3,natom)=contribution of atomic densities to frozen part of dyn. matrix
140 : : !! --- if opteltfr==1
141 : : !! eltfrn(6+3*natom,6)=contribution of atomic density to frozen part of stress tensor
142 : : !!
143 : : !! NOTES
144 : : !! Details on possible options:
145 : : !! ============================
146 : : !! optv: controls the computation of a local potential as sum of atomic potentials
147 : : !! Vloc(r)=Sum_R[V^AT(r-R)]
148 : : !! or its contributions to energy derivatives, i.e. derivatives of Int[Vloc(r).rho(r).dr]
149 : : !! rho(r) is stored in reciprocal space in array rhog()
150 : : !! V^AT is stored in reciprocal space in array vspl (in practice vspl(q)=q^2.V^AT(q))
151 : : !!
152 : : !! optn: controls the computation of a density as sum of atomic densities
153 : : !! n(r)=Sum_R[n^AT(r-R)]
154 : : !! or its contributions to energy derivatives, i.e. derivatives of Int[n(r).V(r).dr]
155 : : !! V(r) is stored in reciprocal space in array vg()
156 : : !! n^AT is stored in reciprocal space:
157 : : !! if optn2=1: n^AT is the atomic PAW PS core density stored in array pawtab%tcorespl()
158 : : !! 2: n^AT is the atomic PAW PS valence density stored in array pawtab%tvalespl()
159 : : !! 3: n^AT is a gaussian density: n(g)=gauss(1,ityp)*exp[-(gauss(2,ityp)*G)^2]
160 : : !! 4: n^AT is the atomic PAW PS core kinetic density stored in array pawtab%ttaucorespl()
161 : : !! Note: optv and optn can be activated together
162 : : !!
163 : : !! Options controlling which contrib. to Etot derivatives are computed:
164 : : !! optatm =1: computes Vloc(r) or n(r) as sum of atomic potentials/densities
165 : : !! optgr =1: computes contribution of atomic Vloc(r) or n(r) to forces
166 : : !! optstr =1: computes contribution of atomic Vloc(r) or n(r) to stress tensor
167 : : !! optdyfr =1: computes contribution of atomic Vloc(r) or n(r) to fr part of dyn. matrix
168 : : !! opteltfr=1: computes contribution of atomic Vloc(r) or n(r) to elastic tensor
169 : : !! Note: optatm, optgr, optstr, optelfr and optdyfr can be activated together
170 : : !!
171 : : !! Typical uses:
172 : : !! =============
173 : : !! Computation of:
174 : : !! - local potential: optv=1, optatm=1
175 : : !! - contrib. of local potential to Etot derivatives: optv=1, rhog=total valence density
176 : : !! optgr=1 or optstr=1 or optdyfr=1
177 : : !! - PS core density: optn=1, optn2=1, optatm=1
178 : : !! - contrib. of NLCC to Etot derivatives: optn=1, optn2=1, vg=XC potential
179 : : !! optgr=1 or optstr=1 or optdyfr=1
180 : : !! - sum of atomic valence densities: optn=1, optn2=2 or 3, optatm=1
181 : : !! - correction of forces due to potential residual: optn=1, optn2=2 or 3, optgr=1
182 : : !! vg=potential residual
183 : : !! etc...
184 : : !!
185 : : !! SOURCE
186 : :
187 : 16294 : subroutine atm2fft(atindx1,atmrho,atmvloc,dyfrn,dyfrv,eltfrn,gauss,gmet,gprimd,&
188 : 16294 : & grn,grv,gsqcut,mgfft,mqgrid,natom,nattyp,nfft,ngfft,ntypat,&
189 : : & optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,&
190 : 16294 : & psps,pawtab,ph1d,qgrid,qprtrb,rcut,rhog,rprimd,strn,strv,ucvol,usepaw,vg,vg1,vg1_core,vprtrb,vspl,&
191 : : & is2_in,comm_fft,me_g0,paral_kgb,distribfft) ! optional arguments
192 : :
193 : : !Arguments ------------------------------------
194 : : !scalars
195 : : integer,intent(in) :: mgfft,mqgrid,natom,nfft,ntypat,optatm,optdyfr,opteltfr
196 : : integer,intent(in) :: optgr,optn,optn2,optstr,optv,usepaw
197 : : integer,optional,intent(in) :: is2_in,me_g0,comm_fft,paral_kgb
198 : : real(dp),intent(in) :: rcut,gsqcut,ucvol
199 : : type(pseudopotential_type),target,intent(in) :: psps
200 : : type(distribfft_type),optional,intent(in),target :: distribfft
201 : : !arrays
202 : : integer,intent(in) :: atindx1(natom),nattyp(ntypat),ngfft(18),qprtrb(3)
203 : : real(dp),intent(in) :: gauss(2,ntypat*(optn2/3)),gmet(3,3),gprimd(3,3)
204 : : real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*natom),qgrid(mqgrid)
205 : : real(dp),intent(in) :: rhog(2,nfft*optv*max(optgr,optstr,optdyfr,opteltfr))
206 : : real(dp),intent(inout) :: rprimd(3,3)
207 : : real(dp),intent(in) :: vg(2,nfft*optn*max(optgr,optstr,optdyfr,opteltfr))
208 : : real(dp),intent(in) :: vg1(2,nfft*optn*opteltfr),vg1_core(2,nfft*optn*opteltfr)
209 : : real(dp),intent(in) :: vprtrb(2),vspl(mqgrid,2,ntypat*optv)
210 : : real(dp),intent(out) :: atmrho(nfft*optn)
211 : : real(dp),intent(inout) :: atmvloc(nfft*optv)
212 : : real(dp),intent(out) :: dyfrn(3,3,natom*optn*optdyfr),dyfrv(3,3,natom*optv*optdyfr)
213 : : real(dp),intent(out) :: eltfrn(6+3*natom,6)
214 : : real(dp),intent(inout) :: grn(3,natom*optn*optgr)
215 : : real(dp),intent(out) :: grv(3,natom*optv*optgr),strn(6*optn*optstr)
216 : : real(dp),intent(out) :: strv(6*optv*optstr)
217 : : type(pawtab_type),target,intent(in) :: pawtab(ntypat*usepaw)
218 : :
219 : : !Local variables ------------------------------
220 : : !scalars
221 : : integer,parameter :: im=2,re=1,icutcoul=3
222 : : integer :: i1,i2,i3,ia,ia1,ia2,id1,id2,id3,ierr,ig1,ig1_,ig2,ig2_,ig3,ig3_,ii,is1,is2
223 : : integer :: itypat,jj,js,ka,kb,kd,kg,me_fft,my_comm_fft,ndir,n1,n2,n3,nproc_fft,paral_kgb_fft
224 : : integer :: shift1,shift2,shift3
225 : : logical :: have_g0
226 : : #ifdef FC_NVHPC
227 : : !Silly trick to prevent NVHPC optimization issue
228 : : logical :: nothing=.false.
229 : : #endif
230 : : real(dp),parameter :: tolfix=1.0000001_dp, vcutgeo(3)=zero
231 : : real(dp) :: aa,alf2pi2,bb,cc,cutoff,dbl_ig1,dbl_ig2,dbl_ig3,dd,dg1,dg2,d2g,diff
232 : : real(dp) :: dn_at,d2n_at,d2n_at2,dq,dq2div6,dqdiv6,dqm1,dv_at,ee,ff,gauss1,gauss2,gg,gmag,gsquar,n_at
233 : : real(dp) :: ph12i,ph12r,ph1i,ph1r,ph2i,ph2r,ph3i,ph3r,sfi,sfr,term,term1,term2,tmpni,tmpnr
234 : : real(dp) :: tmpvi,tmpvr,v_at,xnorm
235 : : character(len=500) :: message
236 : : type(distribfft_type),pointer :: my_distribfft
237 : 16294 : type(distribfft_type),target :: my_distribfft_
238 : 16294 : type(mpi_type) :: mpi_enreg_fft
239 : : !arrays
240 : 16294 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
241 : 16294 : real(dp), contiguous, pointer :: tvalespl(:,:),tcorespl(:,:),ttaucorespl(:,:)
242 : : real(dp), pointer :: dncdq0, dtaucdq0, dnvdq0
243 : : integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/)
244 : : integer :: delta(6)=(/1,1,1,0,0,0/)
245 : : real(dp) :: dgm(3,3,6),d2gm(3,3,6,6),gcart(3),tsec(2)
246 : 16294 : real(dp),allocatable :: dyfrn_indx(:,:,:),dyfrv_indx(:,:,:),grn_indx(:,:)
247 : 16294 : real(dp),allocatable :: grv_indx(:,:),phim_igia(:),phre_igia(:),workn(:,:)
248 : 16294 : real(dp),allocatable :: gcutoff(:), workv(:,:)
249 : : ! *************************************************************************
250 : :
251 : : DBG_ENTER("COLL")
252 : : !Check optional arguments
253 : 16294 : if (present(comm_fft)) then
254 [ - + ]: 16128 : if ((.not.present(paral_kgb)).or.(.not.present(me_g0))) then
255 : 0 : ABI_BUG(' Need paral_kgb and me_g0 with comm_fft !')
256 : : end if
257 : : end if
258 : :
259 : 16294 : n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
260 : 16294 : me_fft=ngfft(11)
261 : 16294 : nproc_fft=ngfft(10)
262 : :
263 : : !Get the distrib associated with this fft_grid
264 [ + + ]: 16294 : if (present(distribfft)) then
265 : : my_distribfft => distribfft
266 : : else
267 : 166 : my_distribfft => my_distribfft_
268 : 166 : call my_distribfft%init_seq('f',n2,n3,'fourdp')
269 : : end if
270 [ + + ]: 16294 : if (n2==my_distribfft%n2_coarse) then
271 : 1936 : fftn2_distrib => my_distribfft%tab_fftdp2_distrib
272 : 1936 : ffti2_local => my_distribfft%tab_fftdp2_local
273 [ + - ]: 14358 : else if (n2 == my_distribfft%n2_fine) then
274 : 14358 : fftn2_distrib => my_distribfft%tab_fftdp2dg_distrib
275 : 14358 : ffti2_local => my_distribfft%tab_fftdp2dg_local
276 : : else
277 : 0 : ABI_BUG("Unable to find an allocated distrib for this fft grid")
278 : : end if
279 : :
280 [ + + ]: 16294 : if (present(is2_in)) then
281 [ - + ]: 72 : if(is2_in<1.or.is2_in>6) then
282 : 0 : ABI_BUG("is2_in must be between 1 and 6")
283 : : else
284 : : ndir = 1
285 : : end if
286 : : ndir = -1
287 : : end if
288 : :
289 : : !Zero out arrays to permit accumulation over atom types
290 [ + + + + ]: 16294 : if (optv==1.and.optatm==1) then
291 [ - + + - : 9807 : ABI_MALLOC(workv,(2,nfft))
+ - + - -
+ - + ]
292 [ + + + + ]: 304390292 : workv(:,:)=zero
293 : : end if
294 [ + + + + ]: 16294 : if (optn==1.and.optatm==1) then
295 [ - + + - : 9462 : ABI_MALLOC(workn,(2,nfft))
+ - + - -
+ - + ]
296 [ + + + + ]: 264198229 : workn(:,:)=zero
297 : : end if
298 [ + + + + ]: 16294 : if (optv==1.and.optgr==1) then
299 [ - + + - : 20082 : ABI_MALLOC(grv_indx,(3,natom))
+ - + - -
+ - + ]
300 [ + + + + ]: 92474 : grv_indx(:,:)=zero
301 : : end if
302 [ + + + + ]: 16294 : if (optn==1.and.optgr==1) then
303 [ - + + - : 25755 : ABI_MALLOC(grn_indx,(3,natom))
+ - + - -
+ - + ]
304 [ + + + + ]: 116993 : grn_indx(:,:)=zero
305 : : end if
306 [ + + + + ]: 16294 : if (optv==1.and.optdyfr==1) then
307 [ - + + - : 135 : ABI_MALLOC(dyfrv_indx,(3,3,natom))
+ - + - -
+ - + ]
308 [ + + + + : 1787 : dyfrv_indx(:,:,:)=zero
+ + ]
309 : : end if
310 [ + + + + ]: 16294 : if (optn==1.and.optdyfr==1) then
311 [ - + + - : 123 : ABI_MALLOC(dyfrn_indx,(3,3,natom))
+ - + - -
+ - + ]
312 [ + + + + : 1705 : dyfrn_indx(:,:,:)=zero
+ + ]
313 : : end if
314 [ + + + + : 24664 : if (optv==1.and.optstr==1) strv(:)=zero
+ + ]
315 [ + + + + : 22114 : if (optn==1.and.optstr==1) strn(:)=zero
+ + ]
316 [ + + + + : 23854 : if (opteltfr==1) eltfrn(:,:) = zero
+ + ]
317 : :
318 : : !Compute 1st and 2nd derivatives of metric tensor wrt all strain components
319 : : !and store for use in inner loop below for elastic tensor.
320 [ + + ]: 16294 : if (opteltfr==1) then
321 : 72 : dgm(:,:,:)=zero
322 : 72 : d2gm(:,:,:,:)=zero
323 : : ! Loop over 2nd strain index
324 [ + + ]: 504 : do is2=1,6
325 : 432 : kg=idx(2*is2-1);kd=idx(2*is2)
326 [ + + ]: 1728 : do jj = 1,3
327 [ + + ]: 5616 : dgm(:,jj,is2)=-(gprimd(kg,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(kg,jj))
328 : : end do
329 : :
330 : : ! Loop over 1st strain index, upper triangle only
331 [ + + ]: 2016 : do is1=1,is2
332 : 1512 : ka=idx(2*is1-1);kb=idx(2*is1)
333 [ + + ]: 6480 : do jj = 1,3
334 [ + + ]: 4536 : if(ka==kg) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
335 [ + + ]: 8640 : & +gprimd(kb,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(kb,jj)
336 [ + + ]: 4536 : if(ka==kd) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
337 [ + + ]: 5184 : & +gprimd(kb,:)*gprimd(kg,jj)+gprimd(kg,:)*gprimd(kb,jj)
338 [ + + ]: 4536 : if(kb==kg) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
339 [ + + ]: 6048 : & +gprimd(ka,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(ka,jj)
340 [ + + ]: 4536 : if(kb==kd) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
341 [ + + ]: 10152 : & +gprimd(ka,:)*gprimd(kg,jj)+gprimd(kg,:)*gprimd(ka,jj)
342 : : end do
343 : : end do !is1
344 : : end do !is2
345 : : end if
346 : :
347 : :
348 : 16294 : dq=(qgrid(mqgrid)-qgrid(1))/dble(mqgrid-1)
349 : 16294 : dqm1=1.0_dp/dq
350 : 16294 : dqdiv6=dq/6.0_dp
351 : 16294 : dq2div6=dq**2/6.0_dp
352 : 16294 : cutoff=gsqcut*tolfix
353 : 16294 : id1=n1/2+2
354 : 16294 : id2=n2/2+2
355 : 16294 : id3=n3/2+2
356 : :
357 [ - + + - : 48882 : ABI_MALLOC(phre_igia,(natom))
+ - + - -
+ - + ]
358 [ - + + - : 32588 : ABI_MALLOC(phim_igia,(natom))
- + - + ]
359 : :
360 : : !Initialize Gcut-off array from m_termcutoff
361 : : !ABI_MALLOC(gcutoff,(nfft))
362 : 16294 : call termcutoff(gcutoff,gsqcut,icutcoul,ngfft,1,rcut,rprimd,vcutgeo)
363 : :
364 : 16294 : ia1=1
365 [ + + ]: 40479 : do itypat=1,ntypat
366 : : ! ia1,ia2 sets range of loop over atoms:
367 : 24185 : ia2=ia1+nattyp(itypat)-1
368 : 24185 : ii=0
369 : :
370 [ + + ]: 24185 : if (optn2==3)then
371 : 2065 : gauss1=gauss(1,itypat)
372 : 2065 : gauss2=gauss(2,itypat)
373 : 2065 : alf2pi2=(two_pi*gauss2)**2
374 : : end if
375 : :
376 [ + + ]: 24185 : if (usepaw == 1) then
377 : 22392 : tcorespl => pawtab(itypat)%tcorespl
378 : 22392 : tvalespl => pawtab(itypat)%tvalespl
379 : 22392 : ttaucorespl => pawtab(itypat)%tcoretauspl
380 : 22392 : dncdq0 => pawtab(itypat)%dncdq0
381 : 22392 : dnvdq0 => pawtab(itypat)%dnvdq0
382 : 22392 : dtaucdq0 => pawtab(itypat)%dtaucdq0
383 : : else
384 : 1793 : tcorespl => psps%nctab(itypat)%tcorespl
385 : 1793 : tvalespl => psps%nctab(itypat)%tvalespl
386 : 1793 : ttaucorespl => psps%nctab(itypat)%ttaucorespl
387 : 1793 : dncdq0 => psps%nctab(itypat)%dncdq0
388 : 1793 : dnvdq0 => psps%nctab(itypat)%dnvdq0
389 : 1793 : dtaucdq0 => psps%nctab(itypat)%dtaucdq0
390 : : end if
391 : :
392 [ + + ]: 762515 : do i3=1,n3
393 : 738330 : ig3=i3-(i3/id3)*n3-1
394 [ - + ]: 738330 : ig3_=ig3;if (ig3_==(n3/2+1)) ig3_=0
395 [ + + ]: 26456967 : do i2=1,n2
396 : 25694452 : ig2=i2-(i2/id2)*n2-1
397 [ - + ]: 25694452 : ig2_=ig2;if (ig2_==(n2/2+1)) ig2_=0
398 [ + + ]: 26432782 : if(fftn2_distrib(i2)==me_fft) then
399 [ + + ]: 1613118068 : do i1=1,n1
400 : 1589338440 : ig1=i1-(i1/id1)*n1-1
401 [ - + ]: 1589338440 : ig1_=ig1;if (ig1_==(n1/2+1)) ig1_=0
402 : 1589338440 : ii=ii+1
403 : 1589338440 : gsquar=gsq_atm(ig1,ig2,ig3)
404 : :
405 : : ! Skip G**2 outside cutoff:
406 [ + + ]: 1613118068 : if (gsquar<=cutoff) then
407 : :
408 : 665896771 : gmag=sqrt(gsquar)
409 [ + + + + : 665896771 : have_g0=(ig1==0.and.ig2==0.and.ig3==0)
+ + ]
410 : :
411 : 665896771 : jj=1+int(gmag*dqm1)
412 : 665896771 : diff=gmag-qgrid(jj)
413 : :
414 : : ! Compute structure factor for all atoms of given type:
415 [ + + ]: 2635595242 : do ia=ia1,ia2
416 : : #ifdef FC_NVHPC
417 : : !Silly trick to prevent NVHPC optimization issue
418 : : if(nothing) write(100,*) shift1,shift2,shift3
419 : : #endif
420 : 1969698471 : shift1=1+n1+(ia-1)*(2*n1+1)
421 : 1969698471 : shift2=1+n2+(ia-1)*(2*n2+1)+natom*(2*n1+1)
422 : 1969698471 : shift3=1+n3+(ia-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1)
423 : 1969698471 : ph1r=ph1d(1,ig1+shift1);ph1i=ph1d(2,ig1+shift1)
424 : 1969698471 : ph2r=ph1d(1,ig2+shift2);ph2i=ph1d(2,ig2+shift2)
425 : 1969698471 : ph3r=ph1d(1,ig3+shift3);ph3i=ph1d(2,ig3+shift3)
426 : 1969698471 : ph12r=ph1r*ph2r-ph1i*ph2i
427 : 1969698471 : ph12i=ph1r*ph2i+ph1i*ph2r
428 : 1969698471 : phre_igia(ia)=ph12r*ph3r-ph12i*ph3i
429 : 2635595242 : phim_igia(ia)=ph12r*ph3i+ph12i*ph3r
430 : : end do
431 : :
432 : : ! Assemble structure factors for this type of atom= sum[exp(-i.piG.R)]
433 [ + + + + : 665896771 : if (optatm==1.or.optstr==1.or.opteltfr==1) then
+ + ]
434 : : sfr=zero;sfi=zero
435 [ + + ]: 433579940 : do ia=ia1,ia2
436 : 330770412 : sfr=sfr+phre_igia(ia)
437 : 996667183 : sfi=sfi-phim_igia(ia)
438 : : end do
439 : : end if
440 : :
441 : : ! Compute V^AT(G) and/or n^AT(G) for given type of atom
442 : : ! Evaluate spline fit: p. 86 Numerical Recipes, Press et al;
443 : : ! NOTE: error in book for sign of "aa" term in derivative;
444 : : ! ! also see splfit routine.
445 [ + + + + ]: 665896771 : if (optv==1.or.optn2/=3) then
446 : 657459598 : bb = diff*dqm1
447 : 657459598 : aa = 1.0_dp-bb
448 : 657459598 : cc = aa*(aa**2-1.0_dp)*dq2div6
449 : 657459598 : dd = bb*(bb**2-1.0_dp)*dq2div6
450 : : end if
451 [ + + ]: 665896771 : if (optv==1) then
452 [ + + ]: 391803386 : if (have_g0) then
453 : : v_at=zero
454 : : else
455 : : v_at=(aa*vspl(jj,1,itypat)+bb*vspl(jj+1,1,itypat)+&
456 : : & cc*vspl(jj,2,itypat)+dd*vspl(jj+1,2,itypat)) &
457 : 391787692 : & /gsquar * gcutoff(ii)
458 : : end if
459 : : end if
460 [ + + ]: 665896771 : if (optn==1) then
461 [ + + ]: 371214334 : if (optn2==1) then
462 : 97584115 : n_at=(aa*tcorespl(jj,1)+bb*tcorespl(jj+1,1)+cc*tcorespl(jj,2)+dd*tcorespl(jj+1,2))
463 [ + + ]: 273630219 : else if (optn2==2) then
464 : 263865722 : n_at=(aa*tvalespl(jj,1)+bb*tvalespl(jj+1,1)+cc*tvalespl(jj,2)+dd*tvalespl(jj+1,2))
465 [ + + ]: 9764497 : else if (optn2==3) then
466 : 8437173 : n_at=gauss1*exp(-gsquar*alf2pi2)
467 [ + - ]: 1327324 : else if (optn2==4) then
468 : 1327324 : n_at=(aa*ttaucorespl(jj,1)+bb*ttaucorespl(jj+1,1)+cc*ttaucorespl(jj,2)+dd*ttaucorespl(jj+1,2))
469 : : else
470 : : n_at=zero
471 : : end if
472 : : !DEBUG MJV for mGGA NC potentials
473 : : !write (1003, *) jj, n_at, '#optn==1 jj n_at xccc/xcctau in recip space'
474 : : end if
475 : :
476 : : ! Compute sum of local atomic potentials or densities
477 : : ! ---------------------------------------------------
478 [ + + ]: 665896771 : if(optatm==1) then
479 : : ! Accumulate V^AT(G)*SF(G) or n^AT(G)*SF(G)
480 [ + + ]: 67676929 : if (optv==1) then
481 : 53554149 : workv(re,ii)=workv(re,ii)+sfr*v_at
482 : 53554149 : workv(im,ii)=workv(im,ii)+sfi*v_at
483 : : end if
484 [ + + ]: 67676929 : if (optn==1) then
485 : 41865241 : workn(re,ii)=workn(re,ii)+sfr*n_at
486 : 41865241 : workn(im,ii)=workn(im,ii)+sfi*n_at
487 : : end if
488 : :
489 : : ! Compute contrib. to forces and/or frozen part of dyn. matrix
490 : : ! -------------------------------------------------------------
491 [ + + + + ]: 598219842 : else if (optgr==1.or.optdyfr==1) then
492 : 563087243 : dbl_ig1=dble(ig1_);dbl_ig2=dble(ig2_);dbl_ig3=dble(ig3_)
493 : : ! Compute (2Pi)*V^AT(G)*rho(G) or (2Pi)*n^AT(G)*V(G)
494 [ + + ]: 563087243 : if (optv==1) then
495 : 303984459 : tmpvr=(two_pi*v_at)*rhog(re,ii)
496 : 303984459 : tmpvi=(two_pi*v_at)*rhog(im,ii)
497 : : end if
498 [ + + ]: 563087243 : if (optn==1) then
499 : 315594398 : tmpnr=(two_pi*n_at)*vg(re,ii)
500 : 315594398 : tmpni=(two_pi*n_at)*vg(im,ii)
501 : : end if
502 : : ! === contrib. to forces
503 [ + + ]: 563087243 : if (optgr==1) then
504 : : ! Accumulate -(2Pi.G)*V^AT(G)*rho(G)*SF(G)
505 : : ! or -(2Pi)*n^AT(G)*V(G)*SF(G) into forces
506 [ + + ]: 562513849 : if (optv==1) then
507 [ + + ]: 1246859326 : do ia=ia1,ia2
508 : 943448261 : term=tmpvi*phre_igia(ia)+tmpvr*phim_igia(ia)
509 : 943448261 : grv_indx(1,ia)=grv_indx(1,ia)-dbl_ig1*term
510 : 943448261 : grv_indx(2,ia)=grv_indx(2,ia)-dbl_ig2*term
511 : 1246859326 : grv_indx(3,ia)=grv_indx(3,ia)-dbl_ig3*term
512 : : end do
513 : : end if
514 [ + + ]: 562513849 : if (optn==1) then
515 [ + + ]: 1274441444 : do ia=ia1,ia2
516 : 959369124 : term=tmpni*phre_igia(ia)+tmpnr*phim_igia(ia)
517 : 959369124 : grn_indx(1,ia)=grn_indx(1,ia)-dbl_ig1*term
518 : 959369124 : grn_indx(2,ia)=grn_indx(2,ia)-dbl_ig2*term
519 : 1274441444 : grn_indx(3,ia)=grn_indx(3,ia)-dbl_ig3*term
520 : : end do
521 : : end if
522 : : end if
523 : : ! === contrib. to frozen part of dyn. matrix
524 [ + + ]: 563087243 : if (optdyfr==1) then
525 : : ! Accumulate -(2Pi^2.Gi.Gj)*V^AT(G)*rho(G)*SF(G)
526 : : ! or -(2Pi^2.Gi.Gj)*n^AT(G)*V(G)*SF(G) into dyn. matrix
527 [ + - ]: 573394 : if (optv==1) then
528 [ + + ]: 2226100 : do ia=ia1,ia2
529 : 1652706 : term=two_pi*(tmpvr*phre_igia(ia)-tmpvi*phim_igia(ia))
530 : 1652706 : dyfrv_indx(1,1,ia)=dyfrv_indx(1,1,ia)-dbl_ig1*dbl_ig1*term
531 : 1652706 : dyfrv_indx(1,2,ia)=dyfrv_indx(1,2,ia)-dbl_ig1*dbl_ig2*term
532 : 1652706 : dyfrv_indx(1,3,ia)=dyfrv_indx(1,3,ia)-dbl_ig1*dbl_ig3*term
533 : 1652706 : dyfrv_indx(2,2,ia)=dyfrv_indx(2,2,ia)-dbl_ig2*dbl_ig2*term
534 : 1652706 : dyfrv_indx(2,3,ia)=dyfrv_indx(2,3,ia)-dbl_ig2*dbl_ig3*term
535 : 2226100 : dyfrv_indx(3,3,ia)=dyfrv_indx(3,3,ia)-dbl_ig3*dbl_ig3*term
536 : : end do
537 : : end if
538 [ + + ]: 573394 : if (optn==1) then
539 [ + + ]: 2075882 : do ia=ia1,ia2
540 : 1553804 : term=two_pi*(tmpnr*phre_igia(ia)-tmpni*phim_igia(ia))
541 : 1553804 : dyfrn_indx(1,1,ia)=dyfrn_indx(1,1,ia)-dbl_ig1*dbl_ig1*term
542 : 1553804 : dyfrn_indx(1,2,ia)=dyfrn_indx(1,2,ia)-dbl_ig1*dbl_ig2*term
543 : 1553804 : dyfrn_indx(1,3,ia)=dyfrn_indx(1,3,ia)-dbl_ig1*dbl_ig3*term
544 : 1553804 : dyfrn_indx(2,2,ia)=dyfrn_indx(2,2,ia)-dbl_ig2*dbl_ig2*term
545 : 1553804 : dyfrn_indx(2,3,ia)=dyfrn_indx(2,3,ia)-dbl_ig2*dbl_ig3*term
546 : 2075882 : dyfrn_indx(3,3,ia)=dyfrn_indx(3,3,ia)-dbl_ig3*dbl_ig3*term
547 : : end do
548 : : end if
549 : : end if
550 : : end if
551 : :
552 : : ! Compute (dV^AT(q)/dq)/q and/or (dn^AT(q)/dq)/q
553 : : ! For stress tensor and/or elastic tensor
554 : : ! ---------------------------------
555 [ + + + + ]: 665896771 : if (optstr==1.or.opteltfr==1) then
556 : : ! Note: correction of Numerical Recipes sign error before (3._dp*aa**2-1._dp)
557 : : ! ee*dqm1 + ff*dqdiv6 is the best estimate of dV(q)/dq from splines
558 [ + + ]: 35132599 : if (optv==1) then
559 [ + + ]: 34264778 : if (have_g0) then
560 : : dv_at=zero
561 : : else
562 : 34262965 : ee=vspl(jj+1,1,itypat)-vspl(jj,1,itypat)
563 : : ff=(3._dp*bb**2-1._dp)*vspl(jj+1,2,itypat)&
564 : 34262965 : & -(3._dp*aa**2-1._dp)*vspl(jj ,2,itypat)
565 : 34262965 : dv_at=((ee*dqm1+ff*dqdiv6)/gmag-2.0_dp*v_at)/gsquar
566 : : end if
567 : : end if
568 [ + + ]: 35132599 : if (optn==1) then
569 [ + + ]: 13754695 : if (have_g0) then
570 [ + + ]: 1396 : if (optn2==1) then
571 : 1367 : dn_at=dncdq0
572 [ - + ]: 29 : else if (optn2==2) then
573 : 0 : dn_at=dnvdq0
574 [ - + ]: 29 : else if (optn2==3) then
575 : 0 : dn_at=-two*gauss1*alf2pi2
576 [ + - ]: 29 : else if (optn2==4) then
577 : 29 : dn_at=dtaucdq0
578 : : end if
579 [ + + ]: 1396 : if (opteltfr==1) then
580 : 132 : d2n_at = 0
581 : : end if
582 : : else
583 [ + + ]: 13753299 : if (optn2==1) then
584 : 13322939 : ee=tcorespl(jj+1,1)-tcorespl(jj,1)
585 : : ff=(3._dp*bb**2-1._dp)*tcorespl(jj+1,2) &
586 : 13322939 : & -(3._dp*aa**2-1._dp)*tcorespl(jj,2)
587 : : ! Also get nc''(q)
588 [ + + ]: 13322939 : if (opteltfr==1) gg=aa*tcorespl(jj,2)+bb*tcorespl(jj+1,2)
589 [ - + ]: 430360 : else if (optn2==2) then
590 : 0 : ee=tvalespl(jj+1,1)-tvalespl(jj,1)
591 : : ff=(3._dp*bb**2-1._dp)*tvalespl(jj+1,2) &
592 : 0 : & -(3._dp*aa**2-1._dp)*tvalespl(jj,2)
593 : : ! Also get nc''(q)
594 [ # # ]: 0 : if (opteltfr==1) then
595 : 0 : gg=aa*tvalespl(jj,2)+bb*tvalespl(jj+1,2)
596 : : end if
597 [ + - ]: 430360 : else if (optn2==3) then
598 : 13753299 : dn_at=-two*gauss1*alf2pi2*exp(-gsquar*alf2pi2)
599 [ + - + - ]: 430360 : else if (optn2==4.and.usepaw==1) then
600 : 430360 : ee=ttaucorespl(jj+1,1)-ttaucorespl(jj,1)
601 : : ff=(3._dp*bb**2-1._dp)*ttaucorespl(jj+1,2) &
602 : 430360 : & -(3._dp*aa**2-1._dp)*ttaucorespl(jj,2)
603 : : ! Also get nc''(q)
604 [ - + ]: 430360 : if (opteltfr==1) gg=aa*ttaucorespl(jj,2)+bb*ttaucorespl(jj+1,2)
605 : : end if
606 : 13753299 : dn_at = (ee*dqm1+ff*dqdiv6)/gmag
607 [ + + ]: 13753299 : if (opteltfr==1) then
608 : 426120 : d2n_at = (gg-dn_at)/gsquar
609 : : d2n_at2 = gg/gmag**3
610 : : end if
611 : : end if
612 : : end if
613 : :
614 : : ! Compute G in cartesian coordinates
615 : : gcart(1)=gprimd(1,1)*dble(ig1)+gprimd(1,2)*dble(ig2)+&
616 : 35132599 : & gprimd(1,3)*dble(ig3)
617 : : gcart(2)=gprimd(2,1)*dble(ig1)+gprimd(2,2)*dble(ig2)+&
618 : 35132599 : & gprimd(2,3)*dble(ig3)
619 : : gcart(3)=gprimd(3,1)*dble(ig1)+gprimd(3,2)*dble(ig2)+&
620 : 35132599 : & gprimd(3,3)*dble(ig3)
621 : : end if
622 : :
623 : : ! Compute contrib. to stress tensor
624 : : ! ---------------------------------
625 [ + + ]: 665896771 : if (optstr==1)then
626 : : ! Accumulate -dV^AT/dG*rho(G)*SF(G)*Gi.Gj/G
627 : : ! or -dn^AT/dG*V(G)*SF(G)*Gi.Gj/G
628 : : ! into stress tensor
629 [ + + ]: 34706347 : if (optv==1) then
630 : 34264778 : term=(rhog(re,ii)*sfr+rhog(im,ii)*sfi)
631 : 34264778 : strv(1)=strv(1)-term*(dv_at*gcart(1)*gcart(1)+v_at)
632 : 34264778 : strv(2)=strv(2)-term*(dv_at*gcart(2)*gcart(2)+v_at)
633 : 34264778 : strv(3)=strv(3)-term*(dv_at*gcart(3)*gcart(3)+v_at)
634 : 34264778 : strv(4)=strv(4)-term*dv_at*gcart(3)*gcart(2)
635 : 34264778 : strv(5)=strv(5)-term*dv_at*gcart(3)*gcart(1)
636 : 34264778 : strv(6)=strv(6)-term*dv_at*gcart(2)*gcart(1)
637 : : end if
638 [ + + ]: 34706347 : if (optn==1) then
639 : 13328443 : term=(vg(re,ii)*sfr+vg(im,ii)*sfi)*dn_at
640 : 13328443 : strn(1)=strn(1)-term*gcart(1)*gcart(1)
641 : 13328443 : strn(2)=strn(2)-term*gcart(2)*gcart(2)
642 : 13328443 : strn(3)=strn(3)-term*gcart(3)*gcart(3)
643 : 13328443 : strn(4)=strn(4)-term*gcart(3)*gcart(2)
644 : 13328443 : strn(5)=strn(5)-term*gcart(3)*gcart(1)
645 : 13328443 : strn(6)=strn(6)-term*gcart(2)*gcart(1)
646 : : end if
647 : : end if
648 : :
649 : : ! Compute contrib. to elastic tensor
650 : : ! ---------------------------------
651 [ + + ]: 665896771 : if (opteltfr==1) then
652 : 426252 : dbl_ig1=dble(ig1_);dbl_ig2=dble(ig2_);dbl_ig3=dble(ig3_)
653 : : ! if (ig1==0 .and. ig2==0 .and. ig3==0) cycle
654 : : ! Compute G*dG/Deps_{\gamme\delta}
655 : 426252 : dg2=0.5_dp*dgsqds_atm(ig1,ig2,ig3,is2_in)
656 : :
657 : 426252 : term = vg (re,ii)*sfr + vg (im,ii)*sfi
658 : 426252 : term1 = vg1(re,ii)*sfr + vg1(im,ii)*sfi
659 : 426252 : term2 = vg1_core(re,ii)*sfr + vg1_core(im,ii)*sfi
660 : :
661 [ + + ]: 2983764 : do is1=1,6
662 : :
663 : : ! Compute G*dG/Deps_{\alpha\beta}
664 : 2557512 : dg1=0.5_dp*dgsqds_atm(ig1,ig2,ig3,is1)
665 : : ! Compute G^2*d2G/Deps_{alphabetagammadelta}
666 : 2557512 : d2g=(0.25_dp*d2gsqds_atm(ig1,ig2,ig3,is1,is2_in))
667 : :
668 : 2557512 : eltfrn(is1,is2_in) = eltfrn(is1,is2_in) + (term*(d2n_at*dg1*dg2 + d2g*dn_at))
669 : 2557512 : eltfrn(is1,is2_in) = eltfrn(is1,is2_in) + 0.5*(term1*dn_at*dg1)
670 : 2557512 : eltfrn(is2_in,is1) = eltfrn(is2_in,is1) + 0.5*(term1*dn_at*dg1)
671 : :
672 [ + + ]: 2557512 : if(is2_in<=3)then
673 : 1278756 : eltfrn(is1,is2_in) = eltfrn(is1,is2_in) - term*dn_at*dg1
674 : : end if
675 [ + + ]: 2557512 : if(is1<=3)then
676 : 1278756 : eltfrn(is1,is2_in) = eltfrn(is1,is2_in) - term*dn_at*dg2
677 : : end if
678 [ + + + + ]: 2983764 : if(is2_in<=3.and.is1<=3)then
679 : 639378 : eltfrn(is1,is2_in) = eltfrn(is1,is2_in) - (term1-term)*n_at
680 : : end if
681 : : end do
682 : :
683 : : ! internal strain
684 [ + + ]: 1261104 : do ia=ia1,ia2
685 : 834852 : js=7+3*(ia-1)
686 : : ! Compute -2pi*G*i*vxcis_core(G)*nat(G)*(exp(-iGr))
687 : 834852 : term=(vg1_core(im,ii)*phre_igia(ia)+vg1_core(re,ii)*phim_igia(ia))*n_at*two_pi
688 : 834852 : eltfrn(js ,is2_in) = eltfrn(js ,is2_in) - dbl_ig1*term
689 : 834852 : eltfrn(js+1,is2_in) = eltfrn(js+1,is2_in) - dbl_ig2*term
690 : 834852 : eltfrn(js+2,is2_in) = eltfrn(js+2,is2_in) - dbl_ig3*term
691 : :
692 : : ! Compute -2pi*G*i*vxc(G)*(dnat*dG/deps-delta*nat(G))*(exp(-iGr))
693 : : term=(vg(im,ii)*phre_igia(ia)+vg(re,ii)*phim_igia(ia))*&
694 : 834852 : & (dn_at*dg2-delta(is2_in)*n_at)*two_pi
695 : 834852 : eltfrn(js ,is2_in) = eltfrn(js ,is2_in) - dbl_ig1*term
696 : 834852 : eltfrn(js+1,is2_in) = eltfrn(js+1,is2_in) - dbl_ig2*term
697 : 1261104 : eltfrn(js+2,is2_in) = eltfrn(js+2,is2_in) - dbl_ig3*term
698 : : end do
699 : : end if
700 : : ! End skip G**2 outside cutoff:
701 : : end if
702 : : ! End loop on n1, n2, n3
703 : : end do
704 : : end if ! this plane is for me_fft
705 : : end do
706 : : end do
707 : :
708 : : ! Symmetrize the dynamical matrix with respect to indices
709 [ + + ]: 24185 : if (optdyfr==1) then
710 [ + - ]: 82 : if (optv==1) then
711 [ + + ]: 216 : do ia=ia1,ia2
712 : 134 : dyfrv_indx(2,1,ia)=dyfrv_indx(1,2,ia)
713 : 134 : dyfrv_indx(3,1,ia)=dyfrv_indx(1,3,ia)
714 : 216 : dyfrv_indx(3,2,ia)=dyfrv_indx(2,3,ia)
715 : : end do
716 : : end if
717 [ + + ]: 82 : if (optn==1) then
718 [ + + ]: 206 : do ia=ia1,ia2
719 : 128 : dyfrn_indx(2,1,ia)=dyfrn_indx(1,2,ia)
720 : 128 : dyfrn_indx(3,1,ia)=dyfrn_indx(1,3,ia)
721 : 206 : dyfrn_indx(3,2,ia)=dyfrn_indx(2,3,ia)
722 : : end do
723 : : end if
724 : : end if
725 : :
726 : 40479 : ia1=ia2+1
727 : :
728 : : ! End loop on type of atoms
729 : : end do
730 : :
731 [ - + ]: 16294 : ABI_FREE(phre_igia)
732 [ - + ]: 16294 : ABI_FREE(phim_igia)
733 [ - + ]: 16294 : ABI_FREE(gcutoff)
734 : :
735 : : !Get local potential or density back to real space
736 [ + + ]: 16294 : if(optatm==1)then
737 : : ! Allow for the addition of a perturbing potential
738 [ + + ]: 3651 : if (optv==1) then
739 [ - + ]: 3269 : if ((vprtrb(1)**2+vprtrb(2)**2) > 1.d-30) then
740 : : ! Find the linear indices which correspond with the input wavevector qprtrb
741 : : ! The double modulus handles both i>=n and i<0, mapping into [0,n-1];
742 : : ! then add 1 to get range [1,n] for each
743 : 0 : i3=1+mod(n3+mod(qprtrb(3),n3),n3)
744 : 0 : i2=1+mod(n2+mod(qprtrb(2),n2),n2)
745 : 0 : i1=1+mod(n1+mod(qprtrb(1),n1),n1)
746 : : ! Compute the linear index in the 3 dimensional array
747 : 0 : ii=i1+n1*((ffti2_local(i2)-1)+(n2/nproc_fft)*(i3-1))
748 : : ! Add in the perturbation at G=qprtrb
749 : 0 : workv(re,ii)=workv(re,ii)+0.5_dp*vprtrb(1)
750 : 0 : workv(im,ii)=workv(im,ii)+0.5_dp*vprtrb(2)
751 : : ! Same thing for G=-qprtrb
752 : 0 : i3=1+mod(n3+mod(-qprtrb(3),n3),n3)
753 : 0 : i2=1+mod(n2+mod(-qprtrb(2),n2),n2)
754 : 0 : i1=1+mod(n1+mod(-qprtrb(1),n1),n1)
755 : : ! ii=i1+n1*((i2-1)+n2*(i3-1))
756 : 0 : workv(re,ii)=workv(re,ii)+0.5_dp*vprtrb(1)
757 : 0 : workv(im,ii)=workv(im,ii)-0.5_dp*vprtrb(2)
758 : : write(message, '(a,1p,2e12.4,a,0p,3i4,a)' )&
759 : 0 : & ' atm2fft: perturbation of vprtrb=', vprtrb,&
760 : 0 : & ' and q=',qprtrb,' has been added'
761 : 0 : call wrtout(std_out,message,'COLL')
762 : : end if
763 : : end if
764 : :
765 [ + + - + ]: 3651 : if (optv==1.or.optn==1) then
766 : : ! Create fake mpi_enreg to wrap fourdp
767 : 3651 : call initmpi_seq(mpi_enreg_fft)
768 [ - + - + : 3651 : ABI_FREE(mpi_enreg_fft%distribfft)
- + - + -
+ - + - +
- + - + -
+ - + - +
- + ]
769 [ + + ]: 3651 : if (present(comm_fft)) then
770 : 3530 : call set_mpi_enreg_fft(mpi_enreg_fft,comm_fft,my_distribfft,me_g0,paral_kgb)
771 : 3530 : my_comm_fft=comm_fft;paral_kgb_fft=paral_kgb
772 : : else
773 : 121 : my_comm_fft=xmpi_comm_self;paral_kgb_fft=0;
774 : 121 : mpi_enreg_fft%distribfft => my_distribfft
775 : : end if
776 : : ! Non-symetrized non-zero elements have to be nullified
777 : : ! Transform back to real space; divide by unit cell volume
778 : 3651 : xnorm=one/ucvol
779 [ + + ]: 3651 : if (optv==1) then
780 : 3269 : call zerosym(workv,2,n1,n2,n3,comm_fft=my_comm_fft,distribfft=my_distribfft)
781 : 3269 : call fourdp(1,workv,atmvloc,1,mpi_enreg_fft,nfft,1,ngfft,0)
782 [ + + ]: 101465610 : atmvloc(:)=atmvloc(:)*xnorm
783 [ - + ]: 3269 : ABI_FREE(workv)
784 : : end if
785 [ + + ]: 3651 : if (optn==1) then
786 : 3154 : call zerosym(workn,2,n1,n2,n3,comm_fft=my_comm_fft,distribfft=my_distribfft)
787 : 3154 : call fourdp(1,workn,atmrho,1,mpi_enreg_fft,nfft,1,ngfft,0)
788 [ + + ]: 88068179 : atmrho(:)=atmrho(:)*xnorm
789 [ - + ]: 3154 : ABI_FREE(workn)
790 : : end if
791 : : ! Destroy fake mpi_enreg
792 : 3651 : call unset_mpi_enreg_fft(mpi_enreg_fft)
793 : : end if
794 : :
795 : : end if
796 : :
797 : : !Additional treatment in case of parallelization
798 [ + + ]: 16294 : if (present(comm_fft)) then
799 [ + + + + : 16128 : if ((xmpi_comm_size(comm_fft)>1).and.&
+ + + - -
+ ]
800 : : & (optgr==1.or.optstr==1.or.optdyfr==1.or.opteltfr==1)) then
801 : 2764 : call timab(48,1,tsec)
802 [ + + ]: 2764 : if (optv==1) then
803 [ + + ]: 1048 : if (optgr==1)then
804 : 812 : call xmpi_sum(grv_indx,comm_fft,ierr)
805 : : end if
806 [ + + ]: 1048 : if (optstr==1)then
807 : 236 : call xmpi_sum(strv,comm_fft,ierr)
808 : : end if
809 [ - + ]: 1048 : if (optdyfr==1)then
810 : 0 : call xmpi_sum(dyfrv_indx,comm_fft,ierr)
811 : : end if
812 : : end if
813 [ + + ]: 2764 : if (optn==1) then
814 [ + + ]: 2748 : if (optgr==1)then
815 : 2520 : call xmpi_sum(grn_indx,comm_fft,ierr)
816 : : end if
817 [ + + ]: 2748 : if (optstr==1)then
818 : 228 : call xmpi_sum(strn,comm_fft,ierr)
819 : : end if
820 [ - + ]: 2748 : if (optdyfr==1)then
821 : 0 : call xmpi_sum(dyfrn_indx,comm_fft,ierr)
822 : : end if
823 : : end if
824 : 2764 : call timab(48,2,tsec)
825 : : end if
826 : : end if
827 : :
828 : : !Forces: re-order atoms
829 [ + + ]: 16294 : if (optgr==1) then
830 [ + + ]: 11099 : if (optv==1) then
831 [ + + ]: 28139 : do ia=1,natom
832 [ + + ]: 92474 : grv(1:3,atindx1(ia))=grv_indx(1:3,ia)
833 : : end do
834 [ - + ]: 6694 : ABI_FREE(grv_indx)
835 : : end if
836 [ + + ]: 11099 : if (optn==1) then
837 [ + + ]: 35687 : do ia=1,natom
838 [ + + ]: 116993 : grn(1:3,atindx1(ia))=grn_indx(1:3,ia)
839 : : end do
840 [ - + ]: 8585 : ABI_FREE(grn_indx)
841 : : end if
842 : : end if
843 : :
844 : : !Elastic tensor: Fill in lower triangle
845 : : ! if (opteltfr==1) then
846 : : ! do is2=2,6
847 : : ! do is1=1,is2-1
848 : : ! eltfrn(is2,is1)=eltfrn(is1,is2)
849 : : ! end do
850 : : ! end do
851 : : ! end if
852 : :
853 : : !Normalize stress tensor:
854 [ + + ]: 16294 : if (optstr==1) then
855 [ + + ]: 1427 : if (optv==1) then
856 [ + + ]: 9765 : strv(:)=strv(:)/ucvol
857 : : end if
858 [ + + + + ]: 7247 : if (optn==1) strn(:)=strn(:)/ucvol
859 : : end if
860 : :
861 : : !Dynamical matrix: re-order atoms
862 [ + + ]: 16294 : if (optdyfr==1) then
863 [ + - ]: 45 : if (optv==1) then
864 [ + + ]: 179 : do ia=1,natom
865 [ + + + + ]: 1787 : dyfrv(1:3,1:3,atindx1(ia))=dyfrv_indx(1:3,1:3,ia)
866 : : end do
867 [ - + ]: 45 : ABI_FREE(dyfrv_indx)
868 : : end if
869 [ + + ]: 45 : if (optn==1) then
870 [ + + ]: 169 : do ia=1,natom
871 [ + + + + ]: 1705 : dyfrn(1:3,1:3,atindx1(ia))=dyfrn_indx(1:3,1:3,ia)
872 : : end do
873 [ - + ]: 41 : ABI_FREE(dyfrn_indx)
874 : : end if
875 : : end if
876 : :
877 [ + + + + : 32754 : if (.not.present(distribfft)) then
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + ]
878 : 166 : call my_distribfft%free()
879 : : end if
880 : :
881 : : DBG_EXIT("COLL")
882 : :
883 : : contains
884 : :
885 : 1589338440 : function gsq_atm(i1,i2,i3)
886 : :
887 : : real(dp) :: gsq_atm
888 : : integer,intent(in) :: i1,i2,i3
889 : : gsq_atm=dble(i1*i1)*gmet(1,1)+dble(i2*i2)*gmet(2,2)+dble(i3*i3)*gmet(3,3) &
890 : 1589338440 : & +two*(dble(i1*i2)*gmet(1,2)+dble(i2*i3)*gmet(2,3)+dble(i3*i1)*gmet(3,1))
891 : 1589338440 : end function gsq_atm
892 : :
893 : 2983764 : function dgsqds_atm(i1,i2,i3,is)
894 : : !Define dG^2/ds based on G space metric derivative
895 : : real(dp) :: dgsqds_atm
896 : : integer,intent(in) :: i1,i2,i3,is
897 : : dgsqds_atm=dble(i1*i1)*dgm(1,1,is)+dble(i2*i2)*dgm(2,2,is)+&
898 : : & dble(i3*i3)*dgm(3,3,is)+&
899 : : & dble(i1*i2)*(dgm(1,2,is)+dgm(2,1,is))+&
900 : : & dble(i1*i3)*(dgm(1,3,is)+dgm(3,1,is))+&
901 : 2983764 : & dble(i2*i3)*(dgm(2,3,is)+dgm(3,2,is))
902 : 2983764 : end function dgsqds_atm
903 : :
904 : 2557512 : function d2gsqds_atm(i1,i2,i3,is1,is2)
905 : : ! Define 2dG^2/ds1ds2 based on G space metric derivative
906 : : real(dp) :: d2gsqds_atm
907 : : integer,intent(in) :: i1,i2,i3,is1,is2
908 : : d2gsqds_atm=dble(i1*i1)*d2gm(1,1,is1,is2)+&
909 : : & dble(i2*i2)*d2gm(2,2,is1,is2)+dble(i3*i3)*d2gm(3,3,is1,is2)+&
910 : : & dble(i1*i2)*(d2gm(1,2,is1,is2)+d2gm(2,1,is1,is2))+&
911 : : & dble(i1*i3)*(d2gm(1,3,is1,is2)+d2gm(3,1,is1,is2))+&
912 : 2557512 : & dble(i2*i3)*(d2gm(2,3,is1,is2)+d2gm(3,2,is1,is2))
913 : 2557512 : end function d2gsqds_atm
914 : :
915 : : end subroutine atm2fft
916 : : !!***
917 : :
918 : : !!****f* ABINIT/dfpt_atm2fft
919 : : !! NAME
920 : : !! dfpt_atm2fft
921 : : !!
922 : : !! FUNCTION
923 : : !! This routine sums 1st-order atomic functions (density or potential)
924 : : !! defined (in rec. space) on a radial grid to get global 1st-order
925 : : !! quantities on the fine FFT grid.
926 : : !!
927 : : !! Possible options:
928 : : !! optn=1: compute a sum of local 1st-order atomic densities
929 : : !! optv=1: compute a sum of local 1st-order atomic potentials
930 : : !!
931 : : !! INPUTS
932 : : !! atindx(natom)=index table for atoms ordered by type
933 : : !! cplex: if 1, real space 1-order functions on FFT grid
934 : : !! distribfft<type(distribfft_type)>=--optional-- contains infos related to FFT parallelism
935 : : !! eei=local pseudopotential part of total energy
936 : : !! gauss(2,ntypat)= params for gaussian atm density (optn2=3) for each atom type
937 : : !! gmet(3,3)=reciprocal space metric
938 : : !! gprimd(3,3)=reciprocal space dimensional primitive translations
939 : : !! gsqcut=cutoff on |G|^2: see setup1 for definition (doubled sphere)
940 : : !! idir=direction of atomic displacement (in case of phonons perturb.)
941 : : !! used only if ndir=1 (see below)
942 : : !! ipert=nindex of perturbation
943 : : !! me_g0=--optional-- 1 if the current process treat the g=0 plane-wave (only needed when comm_fft is present)
944 : : !! mgfft=maximum size of 1D FFTs
945 : : !! comm_fft=--optional-- MPI communicator over FFT components
946 : : !! mqgrid=number of grid pts in q array for f(q) spline.
947 : : !! natom=number of atoms in unit cell.
948 : : !! ndir=number of directions of atomic displacement (in case of phonon):
949 : : !! can be 1 (idir direction in then used) or 3 (all directions)
950 : : !! 6 cartesian strain component 11,22,33,32,31,21 (in case strain perturbation)
951 : : !! nfft=(effective) number of FFT grid points (for this processor)
952 : : !! ngfft(18)=contain all needed information about 3D FFT
953 : : !! nattyp(ntypat)=array describing how many atoms of each type in cell
954 : : !! ntypat=number of types of atoms.
955 : : !! optn,optn2,optv= (see NOTES below)
956 : : !! paral_kgb=--optional-- 1 if "band-FFT" parallelism is activated (only needed when comm_fft is present)
957 : : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
958 : : !! ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phase information
959 : : !! qgrid(mqgrid)=q grid for spline from 0 to qmax.
960 : : !! qphon(3)=wavevector of the phonon
961 : : !! typat(natom)=type of each atom
962 : : !! ucvol=unit cell volume
963 : : !! usepaw= 0 for non paw calculation; =1 for paw calculation
964 : : !! vspl(mqgrid,2,ntypat)=q^2 v(q) spline of an atomic potential
965 : : !! (used only if optv=1)
966 : : !! xred(3,natom)=reduced atomic coordinates
967 : : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
968 : : !!
969 : : !! OUTPUT
970 : : !! ======= if optv==1 =======
971 : : !! atmvlocr1(cplex*nfft,ndir)=sum of local 1st-order atomic potentials in real space
972 : : !! atmvlocg1(2,nfft,ndir)=sum of local 1st-order atomic potentials in G-space
973 : : !! ======= if optn==1 =======
974 : : !! --- if optatm==1
975 : : !! atmrhor1(cplex*nfft,ndir)=sum of 1st-order atomic densities in real space
976 : : !! atmrhog1(2,nfft,ndir)=sum of 1st-order atomic densities in G-space
977 : : !!
978 : : !! NOTES
979 : : !! Details on possible options:
980 : : !! ============================
981 : : !! optv: controls the computation of a local 1st-order potential as sum of atomic potentials
982 : : !! Vloc(r)=Sum_R[V1^AT(r-R)]
983 : : !! optn: controls the computation of a 1st-order density as sum of atomic densities
984 : : !! n(r)=Sum_R[n1^AT(r-R)]
985 : : !! n^AT is stored in reciprocal space:
986 : : !! if optn2=1: n^AT is the atomic PAW PS core density stored in array pawtab%tcorespl()
987 : : !! 2: n^AT is the atomic PAW PS valence density stored in array pawtab%tvalespl()
988 : : !! 3: n^AT is a gaussian density: n(g)=gauss(1,ityp)*exp[-(gauss(2,ityp)*G)^2]
989 : : !! Note: optv and optn can be activated together
990 : : !!
991 : : !! Typical uses:
992 : : !! =============
993 : : !! Computation of:
994 : : !! - 1st-order local potential: optv=1
995 : : !! - 1st-order PS core density: optn=1, optn2=1
996 : : !!
997 : : !! SOURCE
998 : :
999 : 10532 : subroutine dfpt_atm2fft(atindx,cplex,gmet,gprimd,gsqcut,idir,ipert,&
1000 : : & mgfft,mqgrid,natom,ndir,nfft,ngfft,ntypat,&
1001 : 5266 : & ph1d,qgrid,qphon,typat,ucvol,usepaw,xred,psps,pawtab,&
1002 : 8933 : & atmrhor1,atmrhog1,atmvlocr1,atmvlocg1,distribfft,gauss,comm_fft,me_g0,optn_in,&
1003 : 3850 : & optn2_in,optv_in,paral_kgb,vspl) ! optional arguments
1004 : :
1005 : : !Arguments ------------------------------------
1006 : : !scalars
1007 : : integer,intent(in) :: cplex,idir,ipert,mgfft,mqgrid,natom,ndir,nfft,ntypat,usepaw
1008 : : integer,optional,intent(in) :: optn_in,optn2_in,optv_in
1009 : : integer,optional,intent(in) :: me_g0,comm_fft,paral_kgb
1010 : : real(dp),intent(in) :: gsqcut,ucvol
1011 : : type(pseudopotential_type),target,intent(in) :: psps
1012 : : type(distribfft_type),optional,intent(in),target :: distribfft
1013 : : !arrays
1014 : : integer,intent(in) :: atindx(natom),ngfft(18),typat(natom)
1015 : : real(dp),intent(in) :: gmet(3,3),gprimd(3,3)
1016 : : real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*natom),qgrid(mqgrid),qphon(3)
1017 : : real(dp),intent(in) :: xred(3,natom)
1018 : : real(dp),optional,intent(in) :: gauss(2,ntypat),vspl(mqgrid,2,ntypat)
1019 : : real(dp),optional,intent(out) :: atmrhor1(cplex*nfft,ndir)
1020 : : real(dp),optional,intent(out) :: atmrhog1(2,nfft,ndir)
1021 : : real(dp),optional,intent(out) :: atmvlocr1(cplex*nfft,ndir)
1022 : : real(dp),optional,intent(out) :: atmvlocg1(2,nfft,ndir)
1023 : : type(pawtab_type),target,intent(in) :: pawtab(ntypat*usepaw)
1024 : :
1025 : : !Local variables ------------------------------
1026 : : !scalars
1027 : : integer,parameter :: im=2,re=1
1028 : : integer :: i1,i2,i3,ia,ia1,ia2,iatm,iatom,id,id1,id2,id3
1029 : : integer :: ig1,ig1max,ig1min,ig2,ig2max,ig2min,ig3,ig3max,ig3min
1030 : : integer :: ii,itypat,jj,me_fft,my_comm_fft,n1,n2,n3,nattyp,nproc_fft,ntype,paral_kgb_fft
1031 : : integer :: optn,optv,optn2,shift1,shift2,shift3,type1,type2
1032 : : logical :: have_g0,qeq0,qeq05
1033 : : #ifdef FC_NVHPC
1034 : : !Silly trick to prevent NVHPC optimization issue
1035 : : logical :: nothing=.false.
1036 : : #endif
1037 : : real(dp),parameter :: tolfix=1.0000001_dp
1038 : : real(dp) :: aa,alf2pi2,bb,cc,cutoff,dd,diff,dq,dq2div6,dqdiv6,dqm1,ee,ff
1039 : : real(dp) :: gauss1,gauss2,gmag,gq1,gq2,gq3,gsquar,n_at,dn_at,ph12i,ph12r,ph1i
1040 : : real(dp) :: ph1r,ph2i,ph2r,ph3i,ph3r,phqim,phqre,qxred2pi
1041 : : real(dp) :: sfi,sfqi,sfqr,sfr,term_n,term_v,v_at,dv_at,xnorm
1042 : : type(distribfft_type),pointer :: my_distribfft
1043 : 5266 : type(distribfft_type),target :: my_distribfft_
1044 : 5266 : type(mpi_type) :: mpi_enreg_fft
1045 : : !arrays
1046 : 10532 : integer :: eps1(6)=(/1,2,3,2,3,1/),eps2(6)=(/1,2,3,3,1,2/),jdir(ndir)
1047 : 5266 : integer, contiguous, pointer :: fftn2_distrib(:)
1048 : 5266 : real(dp), contiguous, pointer :: tvalespl(:,:),tcorespl(:,:)
1049 : 5266 : real(dp), contiguous, pointer :: ttaucorespl(:,:)
1050 : : real(dp), pointer :: dncdq0, dtaucdq0, dnvdq0
1051 : : real(dp) :: gq(6),gcart(3)
1052 : 5266 : real(dp),allocatable :: phim_igia(:),phre_igia(:),workn(:,:,:),workv(:,:,:)
1053 : : ! *************************************************************************
1054 : :
1055 : : DBG_ENTER("COLL")
1056 : : ABI_NVTX_START_RANGE(NVTX_DFPT_ATM2FFT)
1057 : :
1058 : : ! Check optional arguments
1059 : 5266 : if (present(comm_fft)) then
1060 [ - + ]: 169 : if ((.not.present(paral_kgb)).or.(.not.present(me_g0))) then
1061 : 0 : ABI_BUG('Need paral_kgb and me_g0 with comm_fft !')
1062 : : end if
1063 : : end if
1064 : :
1065 [ - + ]: 5266 : if (present(gauss))then
1066 [ # # ]: 0 : if (.not.present(optn2_in)) then
1067 : : optn2 = 3
1068 : : else
1069 [ # # ]: 0 : if(optn2_in/=3)then
1070 : 0 : ABI_BUG('optn2 must be set to 3!')
1071 : : else
1072 : : optn2 = optn2_in
1073 : : end if
1074 : : end if
1075 : : end if
1076 : :
1077 [ + + ]: 5266 : if (present(atmrhor1))then
1078 [ + + ]: 5169 : if(.not.present(optn_in))then
1079 : : optn = 1
1080 : : else
1081 : 3963 : optn = optn_in
1082 : : end if
1083 [ - + ]: 5169 : if (.not.present(optn2_in)) then
1084 : 0 : ABI_BUG('rho1 calculation need optn2 !')
1085 : : else
1086 : 5169 : optn2 = optn2_in
1087 : : end if
1088 : : else
1089 : : optn = 0
1090 : : optn2 = 0
1091 : : end if
1092 : :
1093 [ + + ]: 5266 : if (present(atmvlocr1))then
1094 [ + + ]: 3764 : if(.not.present(optv_in))then
1095 : : optv = 1
1096 : : else
1097 [ - + ]: 3010 : if(optv_in/=1)then
1098 : 0 : ABI_BUG('optv_in must be set to 1!')
1099 : : else
1100 : : optv = optv_in
1101 : : end if
1102 : : end if
1103 [ - + ]: 3764 : if(.not.present(vspl))then
1104 : 0 : ABI_BUG('vloc1 calculation need vspl!')
1105 : : end if
1106 : : else
1107 : : optv = 0
1108 : : end if
1109 : :
1110 [ + + + + : 5266 : if(ipert==natom+1.or.ipert==natom+2.or.ipert==natom+6.or.ipert==natom+10.or.ipert==natom+11) then
+ + + + +
+ ]
1111 : :
1112 : : ! (In case of d/dk or an electric/magnetic field)
1113 [ + + ]: 830 : if (optn==1) then
1114 [ + + + + ]: 7722862 : atmrhor1(1:cplex*nfft,1:ndir)=zero
1115 [ - + - - : 803 : if (present(atmrhog1)) atmrhog1 = zero
- - - - ]
1116 : : end if
1117 [ + + ]: 830 : if (optv==1) then
1118 [ + + + + ]: 23008230 : atmvlocr1(1:cplex*nfft,1:ndir)=zero
1119 [ - + - - : 759 : if (present(atmvlocg1)) atmvlocg1 = zero
- - - - ]
1120 : : end if
1121 : :
1122 : : else
1123 : :
1124 : : ! Useful quantities
1125 [ + + + + ]: 4436 : if (ipert/=natom+3.and.ipert/=natom+4) then
1126 : 3714 : iatom=ipert;iatm=atindx(iatom)
1127 : 3714 : itypat=typat(iatom)
1128 : : else
1129 : : !sum of all (strain perturbation)
1130 : : iatom = 1
1131 : : iatm = 1
1132 : : itypat = 1
1133 : : end if
1134 : :
1135 : 4436 : n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
1136 : 4436 : me_fft=ngfft(11)
1137 : 4436 : nproc_fft=ngfft(10)
1138 [ + + ]: 4436 : if (ndir==1) then
1139 : 4339 : jdir(1)=idir
1140 : : else
1141 [ + + ]: 466 : do id=1,ndir
1142 : 466 : jdir(id)=id
1143 : : end do
1144 : : end if
1145 : :
1146 : : ! Get the distrib associated with this fft_grid
1147 [ + + ]: 4436 : if (present(distribfft)) then
1148 : : my_distribfft => distribfft
1149 : : else
1150 : 4267 : my_distribfft => my_distribfft_
1151 : 4267 : call my_distribfft_%init_seq('f',n2,n3,'fourdp')
1152 : : end if
1153 [ + + ]: 4436 : if (n2==my_distribfft%n2_coarse) then
1154 : 60 : fftn2_distrib => my_distribfft%tab_fftdp2_distrib
1155 [ + - ]: 4376 : else if (n2 == my_distribfft%n2_fine) then
1156 : 4376 : fftn2_distrib => my_distribfft%tab_fftdp2dg_distrib
1157 : : else
1158 : 0 : ABI_BUG("Unable to find an allocated distrib for this fft grid")
1159 : : end if
1160 : :
1161 : 4436 : qeq0=(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15)
1162 : : qeq05=(abs(abs(qphon(1))-half)<tol12.or. &
1163 : : & abs(abs(qphon(2))-half)<tol12.or. &
1164 [ + + + + : 4436 : & abs(abs(qphon(3))-half)<tol12)
+ + ]
1165 : :
1166 [ - + ]: 4436 : if (nproc_fft>1.and.qeq05) then
1167 : 0 : ABI_ERROR('not compatible with FFT parallelism')
1168 : : end if
1169 : :
1170 [ - + ]: 4436 : if (optn2==3)then
1171 : 0 : gauss1=gauss(1,itypat)
1172 : 0 : gauss2=gauss(2,itypat)
1173 : 0 : alf2pi2=(two_pi*gauss2)**2
1174 : : end if
1175 : :
1176 : 4436 : dq=(qgrid(mqgrid)-qgrid(1))/dble(mqgrid-1)
1177 : 4436 : dqm1=one/dq
1178 : 4436 : dqdiv6=dq/six
1179 : 4436 : dq2div6=dq**2/six
1180 : 4436 : cutoff=gsqcut*tolfix
1181 : 4436 : id1=n1/2+2;id2=n2/2+2;id3=n3/2+2
1182 : 4436 : ig1max=-1;ig2max=-1;ig3max=-1
1183 : 4436 : ig1min=n1;ig2min=n2;ig3min=n3
1184 : :
1185 : : ! Determination of phase qxred*
1186 : : qxred2pi=two_pi*(qphon(1)*xred(1,iatom)+ &
1187 : : & qphon(2)*xred(2,iatom)+ &
1188 : 4436 : & qphon(3)*xred(3,iatom) )
1189 : 4436 : phqre=cos(qxred2pi)
1190 : 4436 : phqim=sin(qxred2pi)
1191 : :
1192 : : ! Zero out temporary arrays
1193 [ + + ]: 4436 : if (optv==1) then
1194 [ - + + - : 12020 : ABI_MALLOC(workv,(2,nfft,ndir))
- + + - +
- + - - +
- + ]
1195 [ + + + + : 195570090 : workv(:,:,:)=zero
+ + ]
1196 : : end if
1197 [ + + ]: 4436 : if (optn==1) then
1198 [ - + + - : 17260 : ABI_MALLOC(workn,(2,nfft,ndir))
- + + - +
- + - - +
- + ]
1199 [ + + + + : 270846053 : workn(:,:,:)=zero
+ + ]
1200 : : end if
1201 : :
1202 [ + + + + ]: 4436 : if (ipert==natom+3.or.ipert==natom+4) then
1203 : 722 : ntype = ntypat
1204 : 722 : ia1=1
1205 : 722 : type1 = 1
1206 : 722 : type2 = ntype
1207 [ - + + - : 2166 : ABI_MALLOC(phre_igia,(natom))
+ - + - -
+ - + ]
1208 [ - + + - : 1444 : ABI_MALLOC(phim_igia,(natom))
- + - + ]
1209 : : else
1210 : 3714 : type1 = itypat
1211 : 3714 : type2 = itypat
1212 : 3714 : ntype = 1
1213 [ - + ]: 3714 : ABI_MALLOC(phre_igia,(iatm:iatm))
1214 [ - + ]: 3714 : ABI_MALLOC(phim_igia,(iatm:iatm))
1215 : : end if
1216 : :
1217 : :
1218 [ + + ]: 9508 : do itypat=type1,type2
1219 : : ! ia1,ia2 sets range of loop over atoms:
1220 [ + + + + ]: 5072 : if (ipert==natom+3.or.ipert==natom+4) then
1221 [ + + + + ]: 6148 : nattyp = count(typat(:)==itypat)
1222 : 1358 : ia2=ia1+nattyp-1
1223 : : else
1224 : : ia1 = iatm
1225 : : ia2 = iatm
1226 : : end if
1227 : :
1228 [ + + ]: 5072 : if (usepaw == 1) then
1229 : 4834 : tcorespl => pawtab(itypat)%tcorespl
1230 : 4834 : tvalespl => pawtab(itypat)%tvalespl
1231 : 4834 : ttaucorespl => pawtab(itypat)%tcoretauspl
1232 : 4834 : dncdq0 => pawtab(itypat)%dncdq0
1233 : 4834 : dnvdq0 => pawtab(itypat)%dnvdq0
1234 : 4834 : dtaucdq0 => pawtab(itypat)%dtaucdq0
1235 : : else
1236 : 238 : tcorespl => psps%nctab(itypat)%tcorespl
1237 : 238 : tvalespl => psps%nctab(itypat)%tvalespl
1238 : 238 : ttaucorespl => psps%nctab(itypat)%ttaucorespl
1239 : 238 : dncdq0 => psps%nctab(itypat)%dncdq0
1240 : 238 : dnvdq0 => psps%nctab(itypat)%dnvdq0
1241 : 238 : dtaucdq0 => psps%nctab(itypat)%dtaucdq0
1242 : : end if
1243 : :
1244 : 5072 : ii=0
1245 [ + + ]: 156089 : do i3=1,n3
1246 : 151017 : ig3=i3-(i3/id3)*n3-1
1247 : 151017 : gq3=dble(ig3)+qphon(3)
1248 : 151017 : gq(3)=gq3
1249 : :
1250 [ + + ]: 3711232 : do i2=1,n2
1251 [ + - ]: 3706160 : if (fftn2_distrib(i2)==me_fft) then
1252 : 3555143 : ig2=i2-(i2/id2)*n2-1
1253 : 3555143 : gq2=dble(ig2)+qphon(2)
1254 : 3555143 : gq(2)=gq2
1255 : :
1256 [ + + ]: 103149548 : do i1=1,n1
1257 : 99594405 : ig1=i1-(i1/id1)*n1-1
1258 : 99594405 : gq1=dble(ig1)+qphon(1)
1259 : 99594405 : gq(1)=gq1
1260 : :
1261 : 99594405 : ii=ii+1
1262 : : ! gsquar=gsq(gq1,gq2,gq3)
1263 : : gsquar=gq1*gq1*gmet(1,1)+gq2*gq2*gmet(2,2)+gq3*gq3*gmet(3,3) &
1264 : 99594405 : & +two*(gq1*gq2*gmet(1,2)+gq2*gq3*gmet(2,3)+gq3*gq1*gmet(3,1))
1265 : :
1266 : :
1267 : : ! Skip G**2 outside cutoff:
1268 [ + + ]: 103149548 : if (gsquar<=cutoff) then
1269 : :
1270 : : ! Identify min/max indexes (to cancel unbalanced contributions later)
1271 [ + + ]: 33867197 : if (qeq05) then
1272 : 1050227 : ig1max=max(ig1max,ig1);ig1min=min(ig1min,ig1)
1273 : 1050227 : ig2max=max(ig2max,ig2);ig2min=min(ig2min,ig2)
1274 : 1050227 : ig3max=max(ig3max,ig3);ig3min=min(ig3min,ig3)
1275 : : end if
1276 : :
1277 : 33867197 : gmag=sqrt(gsquar)
1278 [ + + + + ]: 33867197 : have_g0=(ig1==0.and.ig2==0.and.ig3==0.and.qeq0)
1279 : :
1280 : 33867197 : jj=1+int(gmag*dqm1)
1281 : 33867197 : diff=gmag-qgrid(jj)
1282 : :
1283 : : ! Compute structure factor
1284 [ + + ]: 79490842 : phre_igia(:) = zero
1285 [ + + ]: 79490842 : phim_igia(:) = zero
1286 : :
1287 [ + + ]: 71411794 : do ia=ia1,ia2
1288 : : #ifdef FC_NVHPC
1289 : : !Silly trick to prevent NVHPC optimization issue
1290 : : if(nothing) write(100,*) shift1,shift2,shift3
1291 : : #endif
1292 : 37544597 : shift1=1+n1+(ia-1)*(2*n1+1)
1293 : 37544597 : shift2=1+n2+(ia-1)*(2*n2+1)+natom*(2*n1+1)
1294 : 37544597 : shift3=1+n3+(ia-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1)
1295 : 37544597 : ph1r=ph1d(re,ig1+shift1);ph1i=ph1d(im,ig1+shift1)
1296 : 37544597 : ph2r=ph1d(re,ig2+shift2);ph2i=ph1d(im,ig2+shift2)
1297 : 37544597 : ph3r=ph1d(re,ig3+shift3);ph3i=ph1d(im,ig3+shift3)
1298 : 37544597 : ph12r=ph1r*ph2r-ph1i*ph2i
1299 : 37544597 : ph12i=ph1r*ph2i+ph1i*ph2r
1300 : 37544597 : phre_igia(ia)=ph12r*ph3r-ph12i*ph3i
1301 : 71411794 : phim_igia(ia)=ph12r*ph3i+ph12i*ph3r
1302 : : end do
1303 : : ! Compute V^AT(g+q) and/or n^AT(g+q) for given type of atom
1304 : : ! Evaluate spline fit: p. 86 Numerical Recipes, Press et al;
1305 : : ! Note the error in book for sign of "aa" term in derivative.
1306 [ + - ]: 33867197 : if (optv==1.or.optn2/=3) then
1307 : 33867197 : bb = diff*dqm1
1308 : 33867197 : aa = one-bb
1309 : 33867197 : cc = aa*(aa**2-one)*dq2div6
1310 : 33867197 : dd = bb*(bb**2-one)*dq2div6
1311 : : end if
1312 : :
1313 [ + + ]: 33867197 : if (optv==1) then
1314 [ + + ]: 22633499 : if (have_g0) then
1315 : : v_at=zero
1316 : : dv_at=zero
1317 : : else
1318 : : v_at=(aa*vspl(jj,1,itypat)+bb*vspl(jj+1,1,itypat)+&
1319 : 22630136 : & cc*vspl(jj,2,itypat)+dd*vspl(jj+1,2,itypat))/gsquar
1320 : :
1321 : : ! Also get (dV(q)/dq)/q:
1322 : : ! (note correction of Numerical Recipes sign error
1323 : : ! before (3._dp*aa**2-1._dp)
1324 : 22630136 : ee= vspl(jj+1,1,itypat)-vspl(jj,1,itypat)
1325 : : ff= (3._dp*bb**2-1._dp)*vspl(jj+1,2,itypat) &
1326 : 22630136 : & - (3._dp*aa**2-1._dp)*vspl(jj,2,itypat)
1327 : : dv_at = ( ( ee*dqm1 + ff*dqdiv6 )/gmag&
1328 : 22630136 : & - 2.0_dp*v_at) / gsquar
1329 : : end if
1330 : : end if ! end optv
1331 : :
1332 [ + + ]: 33867197 : if (optn==1) then
1333 [ + + ]: 32780978 : if (optn2==1) then
1334 : 31459217 : n_at=(aa*tcorespl(jj,1)+bb*tcorespl(jj+1,1)+cc*tcorespl(jj,2)+dd*tcorespl(jj+1,2))
1335 [ + - ]: 1321761 : else if (optn2==2) then
1336 : 1321761 : n_at=(aa*tvalespl(jj,1)+bb*tvalespl(jj+1,1)+cc*tvalespl(jj,2)+dd*tvalespl(jj+1,2))
1337 [ # # ]: 0 : else if (optn2==3) then
1338 : 0 : n_at=gauss1*exp(-gsquar*alf2pi2)
1339 : : else
1340 : : n_at=zero
1341 : : end if
1342 : :
1343 : : ! Also get (dn^AT(q)/dq)/q:
1344 [ + + ]: 32780978 : if (have_g0) then
1345 [ + + ]: 4442 : if (optn2==1) then
1346 : 4364 : dn_at=dncdq0
1347 [ + - ]: 78 : else if (optn2==2) then
1348 : 78 : dn_at=dnvdq0
1349 [ # # ]: 0 : else if (optn2==3) then
1350 : 0 : dn_at=-two*gauss1*alf2pi2
1351 : : end if
1352 : : else
1353 [ + + ]: 32776536 : if (optn2==1) then
1354 : 31454853 : ee=tcorespl(jj+1,1)-tcorespl(jj,1)
1355 : : ff=(3._dp*bb**2-1._dp)*tcorespl(jj+1,2) &
1356 : 31454853 : & -(3._dp*aa**2-1._dp)*tcorespl(jj,2)
1357 [ + - ]: 1321683 : else if (optn2==2) then
1358 : 1321683 : ee=tvalespl(jj+1,1)-tvalespl(jj,1)
1359 : : ff=(3._dp*bb**2-1._dp)*tvalespl(jj+1,2) &
1360 : 1321683 : & -(3._dp*aa**2-1._dp)*tvalespl(jj,2)
1361 : : else if (optn2==3) then
1362 : : dn_at=-two*gauss1*alf2pi2*exp(-gsquar*alf2pi2)
1363 : : else
1364 : : end if
1365 : 32776536 : dn_at=(ee*dqm1+ff*dqdiv6)/gmag
1366 : : end if
1367 : : end if ! end optn
1368 : :
1369 [ + + ]: 69790258 : do id=1,ndir
1370 : : sfr=zero;sfi=zero
1371 [ + + ]: 76340722 : do ia=ia1,ia2
1372 [ + + + + ]: 76340722 : if (ipert==natom+3.or.ipert==natom+4) then
1373 : : ! sum[Exp(-i.2pi.g.xred)]
1374 : 9854664 : sfr=sfr+phre_igia(ia)
1375 : 9854664 : sfi=sfi-phim_igia(ia)
1376 : : else
1377 : : ! Exp(-i.2pi.g.xred) * -i.2pi.(g+q)
1378 : 30562997 : sfr=-(two_pi*gq(jdir(id))*phim_igia(ia))
1379 : 30562997 : sfi=-(two_pi*gq(jdir(id))*phre_igia(ia))
1380 : : end if
1381 : : end do
1382 : :
1383 [ + + + + ]: 69790258 : if (ipert/=natom+3.and.ipert/=natom+4) then
1384 : : ! Phonons case
1385 : :
1386 : : ! Exp(-i.2pi.q.xred) => -i.2pi.(g+q).Exp(-i.2pi.(g+q).xred)
1387 : 30562997 : sfqr= sfr*phqre+sfi*phqim
1388 : 30562997 : sfqi=-sfr*phqim+sfi*phqre
1389 : :
1390 [ + + ]: 30562997 : if (optv == 1) then
1391 : 19755551 : workv(re,ii,id) = workv(re,ii,id) + sfqr*v_at
1392 : 19755551 : workv(im,ii,id) = workv(im,ii,id) + sfqi*v_at
1393 : : end if
1394 [ + + ]: 30562997 : if (optn == 1) then
1395 : 28419206 : workn(re,ii,id) = workn(re,ii,id) + sfqr*n_at
1396 : 28419206 : workn(im,ii,id) = workn(im,ii,id) + sfqi*n_at
1397 : : end if
1398 : :
1399 : : else
1400 : : ! Strain case
1401 : :
1402 : : ! Compute G in cartesian coordinates
1403 : : gcart(1)=gprimd(1,1)*dble(ig1)+gprimd(1,2)*dble(ig2)+&
1404 : 5360064 : & gprimd(1,3)*dble(ig3)
1405 : : gcart(2)=gprimd(2,1)*dble(ig1)+gprimd(2,2)*dble(ig2)+&
1406 : 5360064 : & gprimd(2,3)*dble(ig3)
1407 : : gcart(3)=gprimd(3,1)*dble(ig1)+gprimd(3,2)*dble(ig2)+&
1408 : 5360064 : & gprimd(3,3)*dble(ig3)
1409 : :
1410 : : ! Accumulate -dV^AT/dG*rho(G)*SF(G)*Gi.Gj/G
1411 : : ! or -dn^AT/dG*V(G)*SF(G)*Gi.Gj/G
1412 [ + + ]: 5360064 : if (optv==1) then
1413 [ + + ]: 4933812 : if(jdir(id)<=3) then
1414 : 2466906 : term_v = dv_at*gcart(eps1(jdir(id)))*gcart(eps2(jdir(id))) + v_at
1415 : : else
1416 : 2466906 : term_v = dv_at*gcart(eps1(jdir(id)))*gcart(eps2(jdir(id)))
1417 : : end if
1418 : 4933812 : workv(re,ii,id) = workv(re,ii,id) - (sfr*term_v)
1419 : 4933812 : workv(im,ii,id) = workv(im,ii,id) - (sfi*term_v)
1420 : : end if
1421 : :
1422 [ + + ]: 5360064 : if (optn==1) then
1423 [ + + ]: 4361772 : if(jdir(id)<=3) then
1424 : 2180886 : term_n = dn_at*gcart(eps1(jdir(id)))*gcart(eps2(jdir(id))) + n_at
1425 : : else
1426 : 2180886 : term_n = dn_at*gcart(eps1(jdir(id)))*gcart(eps2(jdir(id)))
1427 : : end if
1428 : :
1429 : 4361772 : workn(re,ii,id) = workn(re,ii,id) - (sfr*term_n)
1430 : 4361772 : workn(im,ii,id) = workn(im,ii,id) - (sfi*term_n)
1431 : : end if
1432 : :
1433 : : end if
1434 : : ! End loop on ndir
1435 : :
1436 : : end do
1437 : : ! End skip G**2 outside cutoff
1438 : : end if
1439 : : ! End loop on n1, n2, n3
1440 : : end do
1441 : : end if ! this plane is selected
1442 : : end do
1443 : : end do
1444 : 9508 : ia1=ia2+1
1445 : : end do ! end loop itype
1446 : :
1447 [ - + ]: 4436 : ABI_FREE(phre_igia)
1448 [ - + ]: 4436 : ABI_FREE(phim_igia)
1449 : :
1450 [ + + + + ]: 4436 : if(ipert==natom+3.or.ipert==natom+4) then
1451 : : ! Set Vloc(G=0)=0:
1452 [ + + ]: 722 : if (optv==1) then
1453 [ + + ]: 1430 : workv(re,1,:)=zero
1454 [ + + ]: 1430 : workv(im,1,:)=zero
1455 : : end if
1456 : : end if
1457 : :
1458 : : ! Identify unbalanced g-vectors
1459 [ + + ]: 4436 : if (qeq05) then !This doesn't work in parallel
1460 [ + + ]: 281 : ig1=-1;if (mod(n1,2)==0) ig1=1+n1/2
1461 [ + + ]: 281 : ig2=-1;if (mod(n2,2)==0) ig2=1+n2/2
1462 [ + + ]: 281 : ig3=-1;if (mod(n3,2)==0) ig3=1+n3/2
1463 [ + + ]: 281 : if (abs(abs(qphon(1))-half)<tol12) then
1464 [ - + ]: 58 : if (abs(ig1min)<abs(ig1max)) ig1=abs(ig1max)
1465 [ + + ]: 58 : if (abs(ig1min)>abs(ig1max)) ig1=n1-abs(ig1min)
1466 : : end if
1467 [ + + ]: 281 : if (abs(abs(qphon(2))-half)<tol12) then
1468 [ - + ]: 227 : if (abs(ig2min)<abs(ig2max)) ig2=abs(ig2max)
1469 [ + + ]: 227 : if (abs(ig2min)>abs(ig2max)) ig2=n2-abs(ig2min)
1470 : : end if
1471 [ + + ]: 281 : if (abs(abs(qphon(3))-half)<tol12) then
1472 [ - + ]: 38 : if (abs(ig3min)<abs(ig3max)) ig3=abs(ig3max)
1473 [ + - ]: 38 : if (abs(ig3min)>abs(ig3max)) ig3=n3-abs(ig3min)
1474 : : end if
1475 : : end if
1476 : :
1477 : :
1478 : : ! Get 1st-order potential/density back to real space
1479 : : ! Non-symetrized non-zero elements have to be nullified
1480 : : ! Divide by unit cell volume
1481 : :
1482 [ + - ]: 4436 : if (optv==1.or.optn==1) then
1483 : :
1484 : 4436 : xnorm=one/ucvol
1485 : : ! Create fake mpi_enreg to wrap fourdp
1486 : 4436 : call initmpi_seq(mpi_enreg_fft)
1487 [ - + - + : 4436 : ABI_FREE(mpi_enreg_fft%distribfft)
- + - + -
+ - + - +
- + - + -
+ - + - +
- + ]
1488 [ + + ]: 4436 : if (present(comm_fft)) then
1489 : 169 : call set_mpi_enreg_fft(mpi_enreg_fft,comm_fft,my_distribfft,me_g0,paral_kgb)
1490 : 169 : my_comm_fft=comm_fft;paral_kgb_fft=paral_kgb
1491 : : else
1492 : 4267 : my_comm_fft=xmpi_comm_self;paral_kgb_fft=0;
1493 : 4267 : mpi_enreg_fft%distribfft => my_distribfft
1494 : : end if
1495 : :
1496 [ + + ]: 4436 : if (optv==1) then
1497 [ + + ]: 6282 : do id=1,ndir
1498 : : ! Eliminate unbalanced g-vectors
1499 [ + + ]: 3277 : if (qeq0) then !q=0
1500 : 3031 : call zerosym(workv(:,:,id),2,n1,n2,n3,comm_fft=my_comm_fft,distribfft=my_distribfft)
1501 [ + + ]: 246 : else if (qeq05) then !q=1/2; this doesn't work in parallel
1502 : 132 : call zerosym(workv(:,:,id),2,n1,n2,n3,ig1=ig1,ig2=ig2,ig3=ig3)
1503 : : end if
1504 [ + - ]: 3277 : call fourdp(cplex,workv(:,:,id),atmvlocr1(:,id),1,mpi_enreg_fft,nfft,1,ngfft,0)
1505 [ + + ]: 68451066 : atmvlocr1(:,id)=atmvlocr1(:,id)*xnorm
1506 : : end do
1507 : :
1508 : : !if (present(atmvlocg1)) atmvlocg1 = workv
1509 [ - + ]: 3005 : ABI_FREE(workv)
1510 : : end if
1511 : :
1512 [ + + ]: 4436 : if (optn==1) then
1513 [ + + ]: 8630 : do id=1,ndir
1514 : : ! Eliminate unbalanced g-vectors
1515 [ + + ]: 4315 : if (qeq0) then !q=0
1516 : 3830 : call zerosym(workn(:,:,id),2,n1,n2,n3,comm_fft=my_comm_fft,distribfft=my_distribfft)
1517 [ + + ]: 485 : else if (qeq05) then !q=1/2; this doesn't work in parallel
1518 : 261 : call zerosym(workn(:,:,id),2,n1,n2,n3,ig1=ig1,ig2=ig2,ig3=ig3)
1519 : : end if
1520 [ + - ]: 4315 : call fourdp(cplex,workn(:,:,id),atmrhor1(:,id),1,mpi_enreg_fft,nfft,1,ngfft,0)
1521 [ + + ]: 97217575 : atmrhor1(:,id)=atmrhor1(:,id)*xnorm
1522 : : end do
1523 : : !if (present(atmrhog1)) atmrhog1 = workn
1524 [ - + ]: 4315 : ABI_FREE(workn)
1525 : : end if
1526 : :
1527 : : ! Destroy fake mpi_enreg
1528 : 4436 : call unset_mpi_enreg_fft(mpi_enreg_fft)
1529 : : end if
1530 : :
1531 [ + + ]: 4436 : if (.not.present(distribfft)) then
1532 : 4267 : call my_distribfft%free()
1533 : : end if
1534 : :
1535 : : ! End the condition of non-electric-field
1536 : : end if
1537 : :
1538 : : ABI_NVTX_END_RANGE()
1539 : : DBG_EXIT("COLL")
1540 : :
1541 [ + + + + : 23315 : end subroutine dfpt_atm2fft
+ + - + -
+ + + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + ]
1542 : : !!***
1543 : :
1544 : : end module m_atm2fft
1545 : : !!***
|
