Branch data Line data Source code
1 : : !!****m* ABINIT/m_dfpt_elt
2 : : !! NAME
3 : : !! m_dfpt_elt
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 1998-2026 ABINIT group (DRH, DCA, XG, GM, AR, MB)
10 : : !! This file is distributed under the terms of the
11 : : !! GNU General Public License, see ~abinit/COPYING
12 : : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : : !!
14 : : !! SOURCE
15 : :
16 : : #if defined HAVE_CONFIG_H
17 : : #include "config.h"
18 : : #endif
19 : :
20 : : #include "abi_common.h"
21 : :
22 : : module m_dfpt_elt
23 : :
24 : : use defs_basis
25 : : use m_abicore
26 : : use m_errors
27 : : use m_xmpi
28 : : use m_dtset
29 : :
30 : : use defs_datatypes, only : pseudopotential_type
31 : : use defs_abitypes, only : MPI_type
32 : : use m_time, only : timab
33 : : use m_special_funcs, only : abi_derfc
34 : : use m_geometry, only : metric
35 : : use m_cgtools, only : dotprod_vn
36 : : use m_pawtab, only : pawtab_type,pawtab_free,pawtab_nullify
37 : : use m_pawrad, only : pawrad_type,pawrad_init,pawrad_free
38 : : use m_pawpsp, only : pawpsp_cg
39 : : use m_paw_numeric, only : paw_spline
40 : : use m_spacepar, only : redgr
41 : : use m_atm2fft, only : atm2fft, dfpt_atm2fft
42 : : use m_mkcore, only : dfpt_mkcore
43 : : use m_dfpt_mkvxcstr, only : dfpt_mkvxcstr
44 : : use m_paral_atom, only : get_my_atmtab, free_my_atmtab
45 : : use m_mpinfo, only : ptabs_fourdp, proc_distrb_cycle, proc_distrb_nband
46 : : use m_fftcore, only : sphereboundary
47 : : use m_fft, only : fourdp
48 : : use m_gtermcutoff, only : termcutoff
49 : :
50 : : implicit none
51 : :
52 : : private
53 : : !!***
54 : :
55 : : public :: dfpt_eltfrxc
56 : : public :: dfpt_eltfrloc
57 : : public :: dfpt_eltfrkin
58 : : public :: dfpt_eltfrhar
59 : : public :: elt_ewald
60 : : public :: dfpt_ewald
61 : : public :: dfpt_ewalddq
62 : : public :: dfpt_ewalddqdq
63 : : !!***
64 : :
65 : : contains
66 : : !!***
67 : :
68 : : !!****f* ABINIT/dfpt_eltfrxc
69 : : !! NAME
70 : : !! dfpt_eltfrxc
71 : : !!
72 : : !! FUNCTION
73 : : !! Compute the 2nd derivatives of exchange-correlation energy
74 : : !! with respect to all pairs of strain and strain-atomic displacement
75 : : !! for the frozen wavefunction contribution to the elastic
76 : : !! and internal strain tensors
77 : : !!
78 : : !! INPUTS
79 : : !! atindx(natom)=index table for atoms ordered by type
80 : : !! dtset <type(dataset_type)>=all input variables for this dataset
81 : : !! | natom=number of atoms in unit cell
82 : : !! | nfft=(effective) number of FFT grid points (for this processor)
83 : : !! | nspden=number of spin-density components
84 : : !! | ntypat=number of types of atoms in cell.
85 : : !! | typat(natom)=integer type for each atom in cell
86 : : !! enxc=exchange and correlation energy (hartree)
87 : : !! gsqcut=Fourier cutoff on G^2 for "large sphere" of radius double that of the basis sphere
88 : : !! kxc(nfft,nkxc)=exchange and correlation kernel
89 : : !! mgfft=maximum size of 1D FFTs
90 : : !! mpi_enreg=information about MPI parallelization
91 : : !! ngfft(18)=contain all needed information about 3D FFT,
92 : : !! see ~abinit/doc/variables/vargs.htm#ngfft
93 : : !! ngfftf(18)= -PAW ONLY- contain all needed information about 3D FFT for the fine grid
94 : : !! (ngfftf=ngfft for norm-conserving potential runs)
95 : : !! nkxc=2nd dimension of kxc
96 : : !! n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
97 : : !! n3xccc=dimension of xccc3d (0 if no core charge, nfft otherwise)
98 : : !! nhat(nfft,nspden*nhatdim)= -PAW only- compensation density
99 : : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
100 : : !! ph1d(2,3*(2*mgfft+1)*natom)=1-dim phase (structure factor) information
101 : : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
102 : : !! rhor(nfft,nspden)=electron density in r space
103 : : !! (if spin polarized, array contains total density in first half and
104 : : !! spin-up density in second half)
105 : : !! (for non-collinear magnetism, first element: total density,
106 : : !! 3 next ones: mx,my,mz)
107 : : !! rprimd(3,3)=dimensional primitive translation vectors (bohr)
108 : : !! usexcnhat= -PAW only- 1 if nhat density has to be taken into account in Vxc
109 : : !! vxc(nfft,nspden)=xc potential (spin up in first half and spin down in
110 : : !! second half if nspden=2)
111 : : !! xcccrc(ntypat)=XC core correction cutoff radius (bohr) for each atom type
112 : : !! xccc1d(n1xccc,6,ntypat)=1D core charge function and five derivatives,
113 : : !! for each type of atom, from psp
114 : : !! xccc3d(n3xccc)=3D core electron density for XC core correction, bohr^-3
115 : : !! xred(3,natom)=reduced coordinates for atoms in unit cell
116 : : !!
117 : : !! OUTPUT
118 : : !! eltfrxc(6+3*natom,6) = xc frozen wavefunction contribution to the
119 : : !! elastic tensor
120 : : !!
121 : : !! SIDE EFFECTS
122 : : !!
123 : : !! NOTES
124 : : !! Much of the code in versions of this routine prior to 4.4.5
125 : : !! has been transfered to its child eltxccore.
126 : : !!
127 : : !! SOURCE
128 : :
129 : 37 : subroutine dfpt_eltfrxc(atindx,dtset,eltfrxc,enxc,gsqcut,kxc,mpi_enreg,mgfft,&
130 : 37 : & nattyp,nfft,ngfft,ngfftf,nhat,nkxc,n3xccc,pawtab,ph1d,psps,rhor,rprimd,&
131 : 37 : & usexcnhat,vxc,xccc3d,xred)
132 : :
133 : : !Arguments ------------------------------------
134 : : !type
135 : : !scalars
136 : : integer,intent(in) :: mgfft,n3xccc,nfft,nkxc,usexcnhat
137 : : real(dp),intent(in) :: enxc,gsqcut
138 : : type(MPI_type),intent(in) :: mpi_enreg
139 : : type(dataset_type),intent(in) :: dtset
140 : : type(pseudopotential_type),intent(inout) :: psps
141 : : !arrays
142 : : integer,intent(in) :: atindx(dtset%natom),nattyp(dtset%ntypat),ngfft(18)
143 : : integer,intent(in) :: ngfftf(18)
144 : : real(dp),intent(in) :: nhat(nfft,dtset%nspden*psps%usepaw)
145 : : real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*dtset%natom)
146 : : real(dp),intent(in) :: vxc(nfft,dtset%nspden),xccc3d(n3xccc)
147 : : real(dp),intent(in) :: xred(3,dtset%natom)
148 : : real(dp),intent(in),target :: rhor(nfft,dtset%nspden)
149 : : real(dp),intent(inout) :: kxc(nfft,nkxc)
150 : : real(dp),intent(out) :: eltfrxc(6+3*dtset%natom,6),rprimd(3,3)
151 : : type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*dtset%usepaw)
152 : :
153 : : !Local variables-------------------------------
154 : : !scalars
155 : : integer,parameter :: mshift=401
156 : : integer :: cplex,fgga,ia,idir,ielt,ieltx,ierr,ifft,ii,ipert,is1,is2,ispden,ispden_c,jj,ka,kb
157 : : integer :: kd,kg,n1,n1xccc,n2,n3,n3xccc_loc,optatm,optdyfr,opteltfr,optgr
158 : : integer :: option,optn,optn2,optstr,optv
159 : : logical :: nmxc
160 : : real(dp) :: d2eacc,d2ecdgs2,d2exdgs2,d2gsds1ds2,d2gstds1ds2,decdgs,dexdgs
161 : : real(dp) :: dgsds10,dgsds20,dgstds10,dgstds20,rstep,spnorm,tmp0,tmp0t
162 : : real(dp) :: ucvol,valuei,yp1,ypn
163 : 37 : type(pawrad_type) :: core_mesh
164 : : !arrays
165 : : integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/)
166 : : real(dp) :: corstr(6),dummy6(0),dummy_in(0,0)
167 : : real(dp) :: dummy_out1(0),dummy_out2(0),dummy_out3(0),dummy_out4(0),dummy_out5(0),dummy_out6(0)
168 : 74 : real(dp) :: eltfrxc_test1(6+3*dtset%natom,6),eltfrxc_test2(6+3*dtset%natom,6)
169 : : real(dp) :: gmet(3,3),gprimd(3,3),qphon(3),rmet(3,3),tsec(2)
170 : : real(dp) :: strn_dummy6(0), strv_dummy6(0)
171 : 37 : real(dp),allocatable :: d2gm(:,:,:,:),dgm(:,:,:),eltfrxc_tmp(:,:)
172 : 37 : real(dp),allocatable :: eltfrxc_tmp2(:,:),elt_work(:,:),rho0_redgr(:,:,:)
173 : 37 : real(dp),allocatable :: vxc10(:,:),vxc10_core(:),vxc10_coreg(:,:)
174 : 37 : real(dp),allocatable :: vxc1is_core(:),vxc1is_coreg(:,:),vxc_core(:)
175 : 37 : real(dp),allocatable :: vxc_coreg(:,:),work(:),workgr(:,:),xccc1d(:,:,:)
176 : 37 : real(dp),allocatable :: xccc3d1(:),xccc3d1_temp(:,:),xcccrc(:)
177 : 37 : real(dp),pointer :: rhor_(:,:)
178 : 37 : type(pawtab_type),allocatable :: pawtab_test(:)
179 : :
180 : : ! *************************************************************************
181 : :
182 : : !Initialize variables
183 : 37 : cplex=1
184 : 37 : qphon(:)=zero
185 : 37 : n1=ngfft(1)
186 : 37 : n2=ngfft(2)
187 : 37 : n3=ngfft(3)
188 : :
189 : 37 : n1xccc = psps%n1xccc
190 : 37 : if(psps%usepaw==0)then
191 [ - + + - : 75 : ABI_MALLOC(xcccrc,(dtset%ntypat))
+ - + - -
+ - + ]
192 [ + + + - : 100 : ABI_MALLOC(xccc1d,(n1xccc,6,dtset%ntypat))
- + + - +
- + + - +
- + ]
193 [ + - - + : 89 : xcccrc = psps%xcccrc
- - - - +
+ ]
194 [ + - + - : 375473 : xccc1d = psps%xccc1d
+ - - + -
- - - - -
- - + + +
+ + + ]
195 : : end if
196 : :
197 [ + + + + ]: 37 : if (usexcnhat==0.and.dtset%usepaw==1) then
198 [ - + + - : 8 : ABI_MALLOC(rhor_,(nfft,dtset%nspden))
- + + - +
- + - - +
- + ]
199 [ + - + + : 31112 : rhor_(:,:) = rhor(:,:)-nhat(:,:)
+ + + + +
+ ]
200 : : else
201 : 35 : rhor_ => rhor
202 : : end if
203 : :
204 : : !HACK - should be fixed globally
205 [ + + ]: 37 : if(n1xccc==0) then
206 : 9 : n3xccc_loc=0
207 : : else
208 : 28 : n3xccc_loc=n3xccc
209 : : end if
210 : :
211 [ + + ]: 37 : fgga=0 ; if(nkxc==7.or.nkxc==19) fgga=1
212 [ + + + - ]: 37 : nmxc=(dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4)
213 : :
214 [ + - - + : 185 : ABI_MALLOC(eltfrxc_tmp,(6+3*dtset%natom,6))
+ - + - +
- + - - +
- + ]
215 [ - + + - : 74 : ABI_MALLOC(eltfrxc_tmp2,(6+3*dtset%natom,6))
- + - + ]
216 [ - + + - : 148 : ABI_MALLOC(vxc10,(nfft,dtset%nspden))
- + + - +
- + - - +
- + ]
217 [ - + + - : 111 : ABI_MALLOC(xccc3d1,(cplex*nfft))
+ - + - -
+ - + ]
218 : :
219 [ + + ]: 37 : if(n1xccc/=0) then
220 [ - + + - : 56 : ABI_MALLOC(vxc_core,(nfft))
- + - + ]
221 [ - + + - : 56 : ABI_MALLOC(vxc10_core,(nfft))
- + - + ]
222 [ - + + - : 56 : ABI_MALLOC(vxc1is_core,(nfft))
- + - + ]
223 : :
224 [ + + ]: 28 : if(dtset%nspden==1) then
225 [ + + ]: 175399 : vxc_core(:)=vxc(:,1)
226 : : else
227 [ + + ]: 11666 : vxc_core(:)=0.5_dp*(vxc(:,1)+vxc(:,2))
228 : : end if
229 : : end if
230 : :
231 : : !Compute gmet, gprimd and ucvol from rprimd
232 : 37 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
233 : :
234 : : !For GGA case, prepare quantities needed to evaluate contributions
235 : : !arising from the strain dependence of the gradient operator itself
236 : :
237 [ + + ]: 37 : if(fgga==1) then
238 [ - + + - : 16 : ABI_MALLOC(rho0_redgr,(3,nfft,dtset%nspden))
- + + - +
- + - - +
- + ]
239 [ - + + - : 8 : ABI_MALLOC(work,(nfft))
- + - + ]
240 [ - + + - : 12 : ABI_MALLOC(workgr,(nfft,3))
+ - + - -
+ - + ]
241 : :
242 : : ! Set up metric tensor derivatives
243 [ - + ]: 4 : ABI_MALLOC(dgm,(3,3,6))
244 [ - + ]: 4 : ABI_MALLOC(d2gm,(3,3,6,6))
245 : : ! Loop over 2nd strain index
246 [ + + ]: 28 : do is2=1,6
247 : 24 : kg=idx(2*is2-1);kd=idx(2*is2)
248 [ + + ]: 96 : do jj = 1,3
249 [ + + ]: 312 : dgm(:,jj,is2)=-(gprimd(kg,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(kg,jj))
250 : : end do
251 : :
252 : : ! Loop over 1st strain index
253 [ + + ]: 172 : do is1=1,6
254 : 144 : ka=idx(2*is1-1);kb=idx(2*is1)
255 [ + + + + ]: 1872 : d2gm(:,:,is1,is2)=0._dp
256 [ + + ]: 576 : do jj = 1,3
257 [ + + ]: 432 : if(ka==kg) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
258 [ + + ]: 672 : & +gprimd(kb,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(kb,jj)
259 [ + + ]: 432 : if(ka==kd) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
260 [ + + ]: 480 : & +gprimd(kb,:)*gprimd(kg,jj)+gprimd(kg,:)*gprimd(kb,jj)
261 [ + + ]: 432 : if(kb==kg) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
262 [ + + ]: 480 : & +gprimd(ka,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(ka,jj)
263 [ + + ]: 432 : if(kb==kd) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
264 [ + + ]: 816 : & +gprimd(ka,:)*gprimd(kg,jj)+gprimd(kg,:)*gprimd(ka,jj)
265 : : end do
266 [ + + + + ]: 1896 : d2gm(:,:,is1,is2)=0.5_dp*d2gm(:,:,is1,is2)
267 : : end do
268 : : end do
269 : :
270 : : ! Compute the reduced gradients of the zero-order charge density.
271 : : ! Note that in the spin-polarized case, we are computing the reduced
272 : : ! gradients of 2 X the spin-up or spin-down charge. This simplifies
273 : : ! subsequent code for the non-spin-polarized case.
274 [ + + ]: 4 : if(dtset%nspden==1) then
275 [ + + ]: 12878 : work(:)=rhor_(:,1)
276 : : else
277 [ + + ]: 5833 : work(:)=2.0_dp*rhor_(:,2)
278 : : end if
279 [ + + ]: 4 : if(n1xccc/=0) then
280 [ + + ]: 13834 : work(:)=work(:)+xccc3d(:)
281 : : end if
282 : 4 : call redgr (work,workgr,mpi_enreg,nfft,ngfft)
283 [ + + ]: 18711 : do ifft=1,nfft
284 [ + + ]: 74832 : rho0_redgr(:,ifft,1)=workgr(ifft,:)
285 : : end do
286 [ + + ]: 4 : if(dtset%nspden==2) then
287 [ + + ]: 5833 : work(:)=2.0_dp*(rhor_(:,1)-rhor_(:,2))
288 [ + - ]: 1 : if(n1xccc/=0) then
289 [ + + ]: 5833 : work(:)=work(:)+xccc3d(:)
290 : : end if
291 : 1 : call redgr(work,workgr,mpi_enreg,nfft,ngfft)
292 [ + + ]: 5833 : do ifft=1,nfft
293 [ + + ]: 23329 : rho0_redgr(:,ifft,2)=workgr(ifft,:)
294 : : end do
295 : : end if
296 [ - + ]: 4 : ABI_FREE(work)
297 [ - + ]: 4 : ABI_FREE(workgr)
298 : : end if !GGA
299 : :
300 : : !Null the elastic tensor accumulator
301 [ + + + + : 9849 : eltfrxc(:,:)=zero;eltfrxc_tmp(:,:)=zero;eltfrxc_tmp2(:,:) = zero
+ + + + +
+ + + ]
302 : :
303 : : !Normalization factor
304 [ + + ]: 37 : if(dtset%nspden==1) then
305 : : spnorm=one
306 : : else
307 : 2 : spnorm=half
308 : : end if
309 : :
310 : : !Big loop over 2nd strain index
311 [ + + ]: 259 : do is2=1,6
312 : :
313 : : ! Translate strain index as needed by dfpt_mkcore below.
314 [ + + ]: 222 : if(is2<=3) then
315 : 111 : ipert=dtset%natom+3
316 : 111 : idir=is2
317 : : else
318 : 111 : ipert=dtset%natom+4
319 : 111 : idir=is2-3
320 : : end if
321 : :
322 : : ! Generate first-order core charge for is2 strain if core charges are present.
323 [ + + ]: 222 : if(n1xccc/=0)then
324 : :
325 [ + + - + ]: 168 : if (psps%usepaw==1 .or. psps%nc_xccc_gspace==1) then
326 : : ! Calculation in Reciprocal space for paw or NC with nc_xccc_gspace
327 [ - + + - : 144 : ABI_MALLOC(xccc3d1_temp,(cplex*nfft,1))
- + - + ]
328 [ + + + + ]: 653328 : xccc3d1_temp = zero
329 : :
330 : : call dfpt_atm2fft(atindx,cplex,gmet,gprimd,gsqcut,is2,ipert,&
331 : : & mgfft,psps%mqgrid_vl,dtset%natom,1,nfft,ngfftf,dtset%ntypat,&
332 : : & ph1d,psps%qgrid_vl,qphon,dtset%typat,ucvol,psps%usepaw,xred,psps,pawtab,&
333 : : & atmrhor1=xccc3d1_temp,optn2_in=1,&
334 : : & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,&
335 : 72 : & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft)
336 [ + + ]: 653256 : xccc3d1(:) = xccc3d1_temp(:,1)
337 [ - + ]: 72 : ABI_FREE(xccc3d1_temp)
338 : :
339 : : else
340 : : ! Calculation in direct space for norm conserving:
341 : : call dfpt_mkcore(cplex,idir,ipert,dtset%natom,dtset%ntypat,n1,n1xccc,&
342 : : & n2,n3,qphon,rprimd,dtset%typat,ucvol,&
343 : 96 : & xcccrc,xccc1d,xccc3d1,xred)
344 : : end if
345 : : else
346 [ + + ]: 229116 : xccc3d1(:)=zero
347 : : end if
348 : :
349 : : ! Compute the first-order potentials.
350 : : ! Standard first-order potential for LDA and GGA with core charge
351 [ + + + + ]: 222 : if(fgga==0 .or. (fgga==1 .and. n1xccc/=0)) then
352 : 210 : option=0
353 : : call dfpt_mkvxcstr(cplex,idir,ipert,kxc,mpi_enreg,dtset%natom,nfft,ngfft,nhat,&
354 : : & dummy_in,nkxc,nmxc,dtset%nspden,n3xccc_loc,option,qphon,rhor,rhor,&
355 : 210 : & rprimd,dtset%usepaw,usexcnhat,vxc10,xccc3d1)
356 [ + + ]: 210 : if(n1xccc/=0)then
357 [ + + ]: 168 : if(dtset%nspden==1) then
358 [ + + ]: 1052394 : vxc10_core(:)=vxc10(:,1)
359 [ + + ]: 1052394 : vxc1is_core(:)=vxc10(:,1)
360 : : else
361 [ + + ]: 69996 : vxc10_core(:)=0.5_dp*(vxc10(:,1)+vxc10(:,2))
362 [ + + ]: 69996 : vxc1is_core(:)=0.5_dp*(vxc10(:,1)+vxc10(:,2))
363 : : end if
364 : : end if
365 : : end if
366 : :
367 : : ! For GGA, first-order potential with doubled gradient operator strain
368 : : ! derivative terms needed for elastic tensor but not internal strain.
369 [ + + ]: 222 : if(fgga==1) then
370 : 24 : option=2
371 : : call dfpt_mkvxcstr(cplex,idir,ipert,kxc,mpi_enreg,dtset%natom,nfft,ngfft,nhat,&
372 : : & dummy_in,nkxc,nmxc,dtset%nspden,n3xccc_loc,option,qphon,rhor,rhor,&
373 : 24 : & rprimd,dtset%usepaw,usexcnhat,vxc10,xccc3d1)
374 [ + + ]: 24 : if(n1xccc/=0)then
375 [ + + ]: 12 : if(dtset%nspden==1) then
376 [ + + ]: 48006 : vxc10_core(:)=vxc10(:,1)
377 : : else
378 [ + + ]: 34998 : vxc10_core(:)=0.5_dp*(vxc10(:,1)+vxc10(:,2))
379 : : end if
380 : : end if
381 : : end if
382 : :
383 : :
384 : : ! Additional term for diagonal strains.
385 [ + + ]: 222 : if(is2<=3) then
386 [ + + + + ]: 710862 : vxc10(:,:)=vxc10(:,:)+vxc(:,:)
387 [ + + ]: 111 : if(n1xccc/=0) then
388 [ + + ]: 561195 : vxc10_core(:)=vxc10_core(:)+2.0_dp*vxc_core(:)
389 [ + + ]: 561195 : vxc1is_core(:)=vxc1is_core(:)+vxc_core(:)
390 : : end if
391 : : end if
392 : :
393 : : ! For GGA, compute the contributions from the strain derivatives acting
394 : : ! on the gradient operators.
395 [ + + ]: 222 : if(fgga==1) then
396 : :
397 [ + + ]: 24 : if (dtset%nspden==1) then
398 [ + + ]: 77268 : do ifft=1,nfft
399 : : ! Collect the needed derivatives of Exc. The factors introduced
400 : : ! deal with the difference between density as used here and
401 : : ! spin density as used with these kxc terms in other contexts.
402 : 77250 : dexdgs =half *kxc(ifft,2)
403 : 77250 : d2exdgs2=quarter*kxc(ifft,4)
404 : : ! Loop over 1st strain index
405 [ + + ]: 540768 : do is1=1,6
406 : : ! The notation here is .gs... for the derivatives of the squared-
407 : : ! gradient of (2X) each spin density, and .gst... for the total density.
408 : : dgsds10=zero;dgsds20=zero;d2gsds1ds2=zero
409 [ + + ]: 1854000 : do jj=1,3
410 [ + + ]: 6025500 : do ii=1,3
411 : 4171500 : tmp0=rho0_redgr(ii,ifft,1)*rho0_redgr(jj,ifft,1)
412 : 4171500 : dgsds10=dgsds10+dgm(ii,jj,is1)*tmp0
413 : 4171500 : dgsds20=dgsds20+dgm(ii,jj,is2)*tmp0
414 : 5562000 : d2gsds1ds2=d2gsds1ds2+d2gm(ii,jj,is1,is2)*tmp0
415 : : end do
416 : : end do
417 : : ! Volume derivative terms added
418 [ + + ]: 463500 : if(is1<=3) d2gsds1ds2=d2gsds1ds2+dgsds20
419 [ + + ]: 463500 : if(is2<=3) d2gsds1ds2=d2gsds1ds2+dgsds10
420 : : ! Add the gradient derivative terms to eltfrxc.
421 : 540750 : eltfrxc(is1,is2)=eltfrxc(is1,is2)+d2exdgs2*dgsds10*dgsds20+dexdgs*d2gsds1ds2
422 : : end do !is1
423 : : end do !ifft
424 : :
425 : : else ! nspden==2
426 : :
427 [ + + ]: 18 : do ispden=1,dtset%nspden
428 : 12 : ispden_c=dtset%nspden-ispden+1
429 : :
430 [ + + ]: 70002 : do ifft=1,nfft
431 : :
432 : : ! Collect the needed derivatives of Exc. The factors introduced
433 : : ! deal with the difference between density as used here and
434 : : ! spin density as used with these kxc terms in other contexts.
435 : 69984 : dexdgs =quarter *kxc(ifft,3+ispden)
436 : 69984 : d2exdgs2=quarter*eighth*kxc(ifft,7+ispden)
437 : 69984 : decdgs =eighth *kxc(ifft,10)
438 : 69984 : d2ecdgs2=eighth*eighth *kxc(ifft,13)
439 : :
440 : : ! Loop over 1st strain index
441 [ + + ]: 489900 : do is1=1,6
442 : :
443 : : ! The notation here is .gs... for the derivatives of the squared-
444 : : ! gradient of (2X) each spin density, and .gst... for the total
445 : : ! density. Note the hack that the the total density is given
446 : : ! by the same expression for either the non-polarized or spin-
447 : : ! polarized case, implemented with the "complementary" index ispden_c
448 : : ! in the expression for tmp0t below.
449 : : dgsds10=zero;dgsds20=zero;d2gsds1ds2=zero
450 : : dgstds10=zero;dgstds20=zero;d2gstds1ds2=zero
451 [ + + ]: 1679616 : do jj=1,3
452 [ + + ]: 5458752 : do ii=1,3
453 : 3779136 : tmp0=rho0_redgr(ii,ifft,ispden)*rho0_redgr(jj,ifft,ispden)
454 : : tmp0t=(rho0_redgr(ii,ifft,ispden)+rho0_redgr(ii,ifft,ispden_c))&
455 : 3779136 : & *(rho0_redgr(jj,ifft,ispden)+rho0_redgr(jj,ifft,ispden_c))
456 : 3779136 : dgsds10=dgsds10+dgm(ii,jj,is1)*tmp0
457 : 3779136 : dgsds20=dgsds20+dgm(ii,jj,is2)*tmp0
458 : 3779136 : dgstds10=dgstds10+dgm(ii,jj,is1)*tmp0t
459 : 3779136 : dgstds20=dgstds20+dgm(ii,jj,is2)*tmp0t
460 : 3779136 : d2gsds1ds2=d2gsds1ds2+d2gm(ii,jj,is1,is2)*tmp0
461 : 5038848 : d2gstds1ds2=d2gstds1ds2+d2gm(ii,jj,is1,is2)*tmp0t
462 : : end do
463 : : end do
464 : : ! Volume derivative terms added
465 [ + + ]: 419904 : if(is1<=3) then
466 : 209952 : d2gsds1ds2=d2gsds1ds2+dgsds20
467 : 209952 : d2gstds1ds2=d2gstds1ds2+dgstds20
468 : : end if
469 [ + + ]: 419904 : if(is2<=3) then
470 : 209952 : d2gsds1ds2=d2gsds1ds2+dgsds10
471 : 209952 : d2gstds1ds2=d2gstds1ds2+dgstds10
472 : : end if
473 : :
474 : : ! Add the gradient derivative terms to eltfrxc.
475 : : eltfrxc(is1,is2)=eltfrxc(is1,is2)+spnorm*&
476 : : & (d2exdgs2*(dgsds10 *dgsds20) + dexdgs*d2gsds1ds2&
477 : 489888 : & +d2ecdgs2*(dgstds10*dgstds20)+ decdgs*d2gstds1ds2)
478 : :
479 : : end do !is1
480 : : end do !ifft
481 : : end do !ispden
482 : :
483 : : end if ! nspden
484 : :
485 : : end if !GGA
486 : :
487 : : ! Compute valence electron 1st-order charge contributions. Recall that
488 : : ! the diagonal strain derivatives of the valence charge are minus the
489 : : ! zero-order density. The explicit symmetrization avoids the need
490 : : ! to store vxc10 for strain indices other than is2.
491 : :
492 : : call dotprod_vn(1,rhor_,d2eacc,valuei,nfft,nfft,dtset%nspden,1,&
493 [ - + ]: 222 : & vxc10,ucvol)
494 [ + + ]: 888 : do is1=1,3
495 : 666 : eltfrxc_tmp(is1,is2)=eltfrxc_tmp(is1,is2)-0.5_dp*d2eacc
496 : 888 : eltfrxc_tmp(is2,is1)=eltfrxc_tmp(is2,is1)-0.5_dp*d2eacc
497 : : end do
498 : :
499 : : ! Compute additional core contributions from is1 perturbation
500 : : ! Internal strain terms calculated here.
501 [ + + ]: 222 : if(n1xccc/=0) then
502 : :
503 [ + + - + ]: 168 : if (psps%usepaw==1 .or. psps%nc_xccc_gspace==1) then
504 : : ! Calculation in Reciprocal space for paw or NC with nc_xccc_gspace
505 : 72 : optatm=0;optdyfr=0;optgr=0;optstr=0;optv=0;optn=n3xccc/nfft;optn2=1;opteltfr=1
506 [ - + + - : 216 : ABI_MALLOC(vxc10_coreg,(2,nfft))
+ - + - -
+ - + ]
507 [ - + + - : 144 : ABI_MALLOC(vxc_coreg,(2,nfft))
- + - + ]
508 [ - + + - : 216 : ABI_MALLOC(vxc1is_coreg,(2,nfft))
+ - + - -
+ - + ]
509 : :
510 [ + + + + : 5878872 : vxc10_coreg(:,:)=zero;vxc10_coreg(:,:)=zero;vxc1is_coreg(:,:)=zero;
+ + + + +
+ + + ]
511 : :
512 : : ! Fourier transform of Vxc_core/vxc10_core to use in atm2fft (reciprocal space calculation)
513 : 72 : call fourdp(1,vxc10_coreg,vxc10_core,-1,mpi_enreg,nfft,1, ngfft,0)
514 : 72 : call fourdp(1,vxc_coreg,vxc_core,-1,mpi_enreg,nfft,1, ngfft, 0)
515 : 72 : call fourdp(1,vxc1is_coreg,vxc1is_core,-1,mpi_enreg,nfft,1, ngfft, 0)
516 : :
517 : : call atm2fft(atindx,dummy_out1,dummy_out2,dummy_out3,dummy_out4,eltfrxc_tmp2,dummy_in,gmet,gprimd,&
518 : : & dummy_out5,dummy_out6,gsqcut,mgfft,psps%mqgrid_vl,dtset%natom,nattyp,nfft,ngfft,dtset%ntypat,&
519 : : & optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab,ph1d,psps%qgrid_vl,dtset%qprtrb,&
520 : : & dtset%rcut,dummy_in,rprimd,strn_dummy6,strv_dummy6,ucvol,psps%usepaw,&
521 : : & vxc_coreg,vxc10_coreg,vxc1is_coreg,dtset%vprtrb,psps%vlspl,is2_in=is2,&
522 : : & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,&
523 : 72 : & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft)
524 : :
525 : : ! The indexing array atindx is used to reestablish the correct order of atoms
526 [ + - - + : 288 : ABI_MALLOC(elt_work,(6+3*dtset%natom,6))
+ - + - +
- + - - +
- + ]
527 [ + + + + ]: 3096 : elt_work(1:6,1:6)=eltfrxc_tmp2(1:6,1:6)
528 [ + + ]: 324 : do ia=1,dtset%natom
529 : 252 : ielt=7+3*(ia-1)
530 : 252 : ieltx=7+3*(atindx(ia)-1)
531 [ + + + + ]: 6372 : elt_work(ielt:ielt+2,1:6)=eltfrxc_tmp2(ieltx:ieltx+2,1:6)
532 : : end do
533 [ + + + + ]: 7632 : eltfrxc_tmp2(:,:)=elt_work(:,:)
534 [ - + ]: 72 : ABI_FREE(elt_work)
535 : :
536 : :
537 [ - + ]: 72 : ABI_FREE(vxc10_coreg)
538 [ - + ]: 72 : ABI_FREE(vxc_coreg)
539 [ - + ]: 72 : ABI_FREE(vxc1is_coreg)
540 [ + + + + ]: 7632 : eltfrxc(:,:)= eltfrxc(:,:) + eltfrxc_tmp2(:,:)
541 : :
542 : : else
543 : :
544 : : call eltxccore(eltfrxc,is2,mpi_enreg%my_natom,dtset%natom,nfft,dtset%ntypat,&
545 : : & n1,n1xccc,n2,n3,rprimd,dtset%typat,ucvol,vxc_core,vxc10_core,vxc1is_core,&
546 : 96 : & xcccrc,xccc1d,xred,mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
547 : :
548 : : !DEBUG
549 : : ! TEST ZONE (DO NOT REMOVE) USE TO RECIPROCAL SPACE IN NC CASE
550 [ - + ]: 96 : if (dtset%userid==567) then
551 [ # # # # ]: 0 : eltfrxc_test1(:,is2)=zero;eltfrxc_test2(:,is2)=zero
552 : : call eltxccore(eltfrxc_test1,is2,mpi_enreg%my_natom,dtset%natom,nfft,dtset%ntypat,&
553 : : & n1,n1xccc,n2,n3,rprimd,dtset%typat,ucvol,vxc_core,vxc10_core,vxc1is_core,&
554 : 0 : & xcccrc,xccc1d,xred,mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
555 : : ! if (is2==1) print*,"elt-frxc from eltxccore",is2,eltfrxc_test1(1,1)*ucvol/dble(nfft)
556 [ # # # # : 0 : ABI_MALLOC(pawtab_test,(dtset%ntypat))
# # # # #
# # # # #
# # ]
557 : 0 : call pawtab_nullify(pawtab_test)
558 [ # # ]: 0 : do jj=1,dtset%ntypat
559 : 0 : pawtab_test(jj)%mqgrid=psps%mqgrid_vl
560 [ # # # # : 0 : ABI_MALLOC(pawtab_test(jj)%tcorespl,(pawtab_test(jj)%mqgrid,2))
# # # # #
# # # #
# ]
561 : 0 : rstep=xcccrc(jj)/dble(n1xccc-1)
562 : 0 : call pawrad_init(mesh=core_mesh,mesh_size=n1xccc,mesh_type=1,rstep=rstep)
563 : : call pawpsp_cg(pawtab_test(jj)%dncdq0,pawtab_test(jj)%d2ncdq0,psps%mqgrid_vl,psps%qgrid_vl,&
564 : 0 : & pawtab_test(jj)%tcorespl(:,1),core_mesh,xccc1d(:,1,jj),yp1,ypn)
565 : 0 : call paw_spline(psps%qgrid_vl,pawtab_test(jj)%tcorespl(:,1),psps%mqgrid_vl,yp1,ypn,pawtab_test(jj)%tcorespl(:,2))
566 : : ! if (is2==1) then
567 : : ! do ii=1,n1xccc;write(100+jj,*) (ii-1)*rstep,xccc1d(ii,1,jj);enddo
568 : : ! do ii=1,psps%mqgrid_vl;write(200+jj,*) psps%qgrid_vl(ii),pawtab_test(jj)%tcorespl(ii,1);enddo
569 : : ! end if
570 : : end do
571 [ # # # # : 0 : ABI_MALLOC(vxc10_coreg,(2,nfft))
# # # # #
# # # ]
572 [ # # # # : 0 : ABI_MALLOC(vxc_coreg,(2,nfft))
# # # # ]
573 [ # # # # : 0 : ABI_MALLOC(vxc1is_coreg,(2,nfft))
# # # # ]
574 [ # # # # : 0 : vxc10_coreg(:,:)=zero;vxc10_coreg(:,:)=zero;vxc1is_coreg(:,:)=zero;
# # # # #
# # # ]
575 : 0 : call fourdp(1,vxc10_coreg,vxc10_core,-1,mpi_enreg,nfft,1, ngfft, 0)
576 : 0 : call fourdp(1,vxc_coreg,vxc_core,-1,mpi_enreg,nfft,1, ngfft,0)
577 : 0 : call fourdp(1,vxc1is_coreg,vxc1is_core,-1,mpi_enreg,nfft,1, ngfft, 0)
578 : 0 : optatm=0;optdyfr=0;optgr=0;optstr=0;optv=0;optn=1;optn2=1;opteltfr=1;corstr=zero
579 : : call atm2fft(atindx,dummy_out1,dummy_out2,dummy_out3,dummy_out4,eltfrxc_test2,dummy_in,gmet,gprimd,&
580 : : & dummy_out5,dummy_out6,gsqcut,mgfft,psps%mqgrid_vl,dtset%natom,nattyp,nfft,ngfft,dtset%ntypat,&
581 : : & optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab_test,ph1d,psps%qgrid_vl,dtset%qprtrb,&
582 : : & dtset%rcut,dummy_in,rprimd,corstr,dummy6,ucvol,psps%usepaw,&
583 : 0 : & vxc_coreg,vxc10_coreg,vxc1is_coreg,dtset%vprtrb,psps%vlspl,is2_in=is2)
584 [ # # ]: 0 : ABI_FREE(vxc10_coreg)
585 [ # # ]: 0 : ABI_FREE(vxc_coreg)
586 [ # # ]: 0 : ABI_FREE(vxc1is_coreg)
587 : 0 : call pawrad_free(core_mesh)
588 : 0 : call pawtab_free(pawtab_test)
589 [ # # # # : 0 : ABI_FREE(pawtab_test)
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # ]
590 [ # # # # ]: 0 : eltfrxc(:,:)= eltfrxc(:,:)+eltfrxc_test2(:,:)
591 : : ! if (is2==1) print*,"cor-str from atm2fft",is2,corstr*ucvol
592 : : ! if (is2==1) print*,"elt-frxc from atm2fft ",is2,eltfrxc_test2(1,1)
593 : : end if
594 : : !DEBUG
595 : :
596 : : end if
597 : : end if
598 : :
599 : : ! Additional term for diagonal strains
600 [ + + ]: 259 : if(is2<=3) then
601 [ + + ]: 444 : do is1=1,3
602 : 444 : eltfrxc_tmp(is1,is2)=eltfrxc_tmp(is1,is2)+enxc
603 : : end do
604 : : end if
605 : : end do !is2 outermost strain loop
606 : :
607 : : !Accumulate eltfrxc accross processors
608 : 37 : call timab(48,1,tsec)
609 : 37 : call xmpi_sum(eltfrxc,mpi_enreg%comm_fft,ierr)
610 : 37 : call timab(48,2,tsec)
611 : :
612 : : !Normalize accumulated 2nd derivatives in NC case
613 [ + + ]: 37 : if(psps%usepaw==1)then
614 [ + + + + ]: 1272 : eltfrxc(:,:)=eltfrxc_tmp(:,:)+eltfrxc
615 : : else
616 [ + + + + ]: 2011 : eltfrxc(:,:)=eltfrxc_tmp(:,:)+eltfrxc*ucvol/dble(nfft)
617 : : end if
618 : :
619 [ - + ]: 37 : ABI_FREE(eltfrxc_tmp)
620 [ - + ]: 37 : ABI_FREE(eltfrxc_tmp2)
621 [ - + ]: 37 : ABI_FREE(vxc10)
622 [ - + ]: 37 : ABI_FREE(xccc3d1)
623 [ + + ]: 37 : if(psps%usepaw==0)then
624 [ - + ]: 25 : ABI_FREE(xccc1d)
625 [ - + ]: 25 : ABI_FREE(xcccrc)
626 : : end if
627 [ + + + + ]: 37 : if (usexcnhat==0.and.dtset%usepaw==1) then
628 [ - + ]: 2 : ABI_FREE(rhor_)
629 : : end if
630 : :
631 [ + + ]: 37 : if(n1xccc/=0) then
632 [ - + ]: 28 : ABI_FREE(vxc_core)
633 [ - + ]: 28 : ABI_FREE(vxc10_core)
634 [ - + ]: 28 : ABI_FREE(vxc1is_core)
635 : : end if
636 : :
637 [ + + ]: 37 : if(fgga==1) then
638 [ - + ]: 4 : ABI_FREE(rho0_redgr)
639 [ - + ]: 4 : ABI_FREE(dgm)
640 [ - + ]: 4 : ABI_FREE(d2gm)
641 : : end if
642 : :
643 [ + + - + : 74 : end subroutine dfpt_eltfrxc
- + - + -
+ - + - +
- + - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- + - + -
+ - + -
+ ]
644 : : !!***
645 : :
646 : : !!****f* ABINIT/eltxccore
647 : : !! NAME
648 : : !! eltxccore
649 : : !!
650 : : !! FUNCTION
651 : : !! Compute the core charge contributions to the 2nd derivatives
652 : : !! of the exchange-correlation energy with respect to all pairs of
653 : : !! strain or strain and atomic displacement for the frozen wavefunction
654 : : !! contribution to the elastic tensor. 1st-order potentials representing
655 : : !! the perturbation by one strain are supplied, and the routine loops
656 : : !! over the second strain and over all atomic displacements.
657 : : !!
658 : : !! COPYRIGHT
659 : : !! Copyright (C) 1998-2026 ABINIT group (DRH, DCA, XG, GMR)
660 : : !! This file is distributed under the terms of the
661 : : !! GNU General Public License, see ~abinit/COPYING
662 : : !! or http://www.gnu.org/copyleft/gpl.txt .
663 : : !!
664 : : !! INPUTS
665 : : !! mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
666 : : !! comm_atom=--optional-- MPI communicator over atoms
667 : : !! my_natom=number of atoms treated by current processor
668 : : !! natom=number of atoms in cell.
669 : : !! nfft=number of fft grid points
670 : : !! ntypat=number of types of atoms in cell.
671 : : !! n1,n2,n3=fft grid dimensions.
672 : : !! n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
673 : : !! rprimd(3,3)=dimensional primitive translation vectors (bohr)
674 : : !! typat(natom)=integer type for each atom in cell
675 : : !! ucvol=unit cell volume (bohr**3).
676 : : !! vxc_core(nfft)=spin-averaged xc potential
677 : : !! vxc10_core(nfft)=spin-averaged 1st-order xc potential for elastic tensor
678 : : !! vxc1is_core(nfft)=spin-averaged 1st-order xc potential for internal strain
679 : : !! xcccrc(ntypat)=XC core correction cutoff radius (bohr) for each atom type
680 : : !! xccc1d(n1xccc,6,ntypat)=1D core charge function and five derivatives,
681 : : !! for each type of atom, from psp
682 : : !! xred(3,natom)=reduced coordinates for atoms in unit cell
683 : : !!
684 : : !! OUTPUT
685 : : !! eltfrxc(6+3*natom,6) = xc frozen wavefunction contribution to the
686 : : !! elastic tensor
687 : : !!
688 : : !! SIDE EFFECTS
689 : : !! eltfrxc(6+3*natom,6) = xc frozen wavefunction contribution to the
690 : : !! elastic and internal-strain tensor. One column is incremented
691 : : !! by the core contribution.
692 : : !!
693 : : !! NOTES
694 : : !! Note that this routine is related to the mkcore.f routine
695 : : !!
696 : : !! SOURCE
697 : :
698 : 192 : subroutine eltxccore(eltfrxc,is2_in,my_natom,natom,nfft,ntypat,&
699 : 96 : & n1,n1xccc,n2,n3,rprimd,typat,ucvol,vxc_core,vxc10_core,vxc1is_core,&
700 : 96 : & xcccrc,xccc1d,xred, &
701 [ + - - + : 96 : & mpi_atmtab,comm_atom) ! optional arguments (parallelism)
- - ]
702 : :
703 : : !Arguments ------------------------------------
704 : : !scalars
705 : : integer,intent(in) :: is2_in,n1,n1xccc,n2,n3,my_natom,natom,nfft,ntypat
706 : : integer,optional,intent(in) :: comm_atom
707 : : real(dp),intent(in) :: ucvol
708 : : !arrays
709 : : integer,intent(in) :: typat(natom)
710 : : integer,optional,target,intent(in) :: mpi_atmtab(:)
711 : : real(dp),intent(in) :: vxc10_core(nfft),vxc1is_core(nfft)
712 : : real(dp),intent(in) :: vxc_core(nfft),xccc1d(n1xccc,6,ntypat)
713 : : real(dp),intent(in) :: xcccrc(ntypat),xred(3,natom)
714 : : real(dp),intent(inout) :: eltfrxc(6+3*natom,6),rprimd(3,3)
715 : :
716 : : !Local variables-------------------------------
717 : : !scalars
718 : : integer,parameter :: mshift=401
719 : : integer :: i1,i2,i3,iat,iatom,ierr,ifft,is1,is2,ishift,ishift1,ishift2
720 : : integer :: ishift3,ixp,jj,js,ka,kb,kd,kg,mu,my_comm_atom,nu
721 : : logical :: my_atmtab_allocated,paral_atom
722 : : real(dp) :: aa,bb,cc,d2rss,dd,delta,delta2div6,deltam1,diff
723 : : real(dp) :: difmag,difmag2,difmag2_fact,difmag2_part,drss1,drss2,func1
724 : : real(dp) :: func2,range,range2,rangem1,rdiff1,rdiff2,rdiff3
725 : : real(dp) :: term1,term2,yy
726 : : character(len=500) :: message
727 : : !arrays
728 : : integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/)
729 : : integer :: igrid(3),ii(mshift,3),irange(3),ngfft(3)
730 : 96 : integer,pointer :: my_atmtab(:)
731 : 192 : real(dp) :: drm(3,3,6),eltfrxc_core(6+3*natom,6),lencp(3),rmet(3,3),rrdiff(mshift,3)
732 : : real(dp) :: scale(3),tau(3),ts2(3),tsec(2),tt(3)
733 : 96 : real(dp),allocatable :: d2rm(:,:,:,:)
734 : :
735 : : ! *************************************************************************
736 : :
737 : : !Compute lengths of cross products for pairs of primitive
738 : : !translation vectors (used in setting index search range below)
739 : : lencp(1)=cross_elt(rprimd(1,2),rprimd(2,2),rprimd(3,2),&
740 : 96 : & rprimd(1,3),rprimd(2,3),rprimd(3,3))
741 : : lencp(2)=cross_elt(rprimd(1,3),rprimd(2,3),rprimd(3,3),&
742 : 96 : & rprimd(1,1),rprimd(2,1),rprimd(3,1))
743 : : lencp(3)=cross_elt(rprimd(1,1),rprimd(2,1),rprimd(3,1),&
744 : 96 : & rprimd(1,2),rprimd(2,2),rprimd(3,2))
745 : :
746 : : !Set up parallelism over atoms
747 [ + - ]: 96 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
748 [ - + ]: 96 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
749 [ + - ]: 96 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
750 : 96 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
751 : :
752 : : !Compute factor R1.(R2xR3)/|R2xR3| etc for 1, 2, 3
753 : : !(recall ucvol=R1.(R2xR3))
754 [ + + ]: 384 : scale(:)=ucvol/lencp(:)
755 : :
756 : : !Compute metric tensor in real space rmet
757 [ + + ]: 384 : do nu=1,3
758 : : rmet(:,nu)=rprimd(1,:)*rprimd(1,nu)+rprimd(2,:)*rprimd(2,nu)+&
759 [ + + ]: 1248 : & rprimd(3,:)*rprimd(3,nu)
760 : : end do
761 : :
762 : : !Compute 1st and 2nd derivatives of metric tensor wrt all strain components
763 : : !and store for use in inner loop below.
764 : :
765 [ - + ]: 96 : ABI_MALLOC(d2rm,(3,3,6,6))
766 : :
767 : : !Loop over 2nd strain index
768 [ + + ]: 672 : do is2=1,6
769 : 576 : kg=idx(2*is2-1);kd=idx(2*is2)
770 [ + + ]: 2304 : do jj = 1,3
771 [ + + ]: 7488 : drm(:,jj,is2)=rprimd(kg,:)*rprimd(kd,jj)+rprimd(kd,:)*rprimd(kg,jj)
772 : : end do
773 : :
774 : : ! Loop over 1st strain index
775 [ + + ]: 4128 : do is1=1,6
776 : :
777 : 3456 : ka=idx(2*is1-1);kb=idx(2*is1)
778 [ + + + + ]: 44928 : d2rm(:,:,is1,is2)=0._dp
779 [ + + ]: 14400 : do jj = 1,3
780 [ + + ]: 10368 : if(ka==kg) d2rm(:,jj,is1,is2)=d2rm(:,jj,is1,is2)&
781 [ + + ]: 16128 : & +rprimd(kb,:)*rprimd(kd,jj)+rprimd(kd,:)*rprimd(kb,jj)
782 [ + + ]: 10368 : if(ka==kd) d2rm(:,jj,is1,is2)=d2rm(:,jj,is1,is2)&
783 [ + + ]: 11520 : & +rprimd(kb,:)*rprimd(kg,jj)+rprimd(kg,:)*rprimd(kb,jj)
784 [ + + ]: 10368 : if(kb==kg) d2rm(:,jj,is1,is2)=d2rm(:,jj,is1,is2)&
785 [ + + ]: 11520 : & +rprimd(ka,:)*rprimd(kd,jj)+rprimd(kd,:)*rprimd(ka,jj)
786 [ + + ]: 10368 : if(kb==kd) d2rm(:,jj,is1,is2)=d2rm(:,jj,is1,is2)&
787 [ + + ]: 19584 : & +rprimd(ka,:)*rprimd(kg,jj)+rprimd(kg,:)*rprimd(ka,jj)
788 : : end do
789 : : end do !is1
790 : : end do !is2
791 : :
792 : 96 : ngfft(1)=n1
793 : 96 : ngfft(2)=n2
794 : 96 : ngfft(3)=n3
795 : 96 : delta=1.0_dp/(n1xccc-1)
796 : 96 : deltam1=n1xccc-1
797 : 96 : delta2div6=delta**2/6.0_dp
798 : :
799 : : !Loop over atoms in unit cell
800 [ + + + + ]: 7476 : eltfrxc_core(:,:)=zero
801 : :
802 [ + + ]: 282 : do iat=1,my_natom
803 [ - + ]: 186 : iatom=iat;if (paral_atom) iatom=my_atmtab(iat)
804 : 186 : js=7+3*(iatom-1)
805 : : ! Set search range (density cuts off perfectly beyond range)
806 : : ! Cycle if no range.
807 : 186 : range=0.0_dp
808 : 186 : range=xcccrc(typat(iatom))
809 [ + + ]: 186 : if(range<1.d-16) cycle
810 : :
811 : 180 : range2=range**2
812 : 180 : rangem1=1.0_dp/range
813 : :
814 : : ! Consider each component in turn
815 [ + + ]: 720 : do mu=1,3
816 : 540 : tau(mu)=mod(xred(mu,iatom)+1._dp-aint(xred(mu,iatom)),1._dp)
817 : :
818 : : ! Use tau to find nearest grid point along R(mu)
819 : : ! (igrid=0 is the origin; shift by 1 to agree with usual index)
820 : 540 : igrid(mu)=nint(tau(mu)*dble(ngfft(mu)))
821 : :
822 : : ! Use range to compute an index range along R(mu)
823 : : ! (add 1 to make sure it covers full range)
824 : 540 : irange(mu)=1+nint((range/scale(mu))*dble(ngfft(mu)))
825 : :
826 : : ! Check that the largest range is smallest than the maximum
827 : : ! allowed one
828 [ - + ]: 540 : if(2*irange(mu)+1 > mshift)then
829 : 0 : write(message, '(a,i0,a)' )' The range around atom',iatom,' is too large.'
830 : 0 : ABI_BUG(message)
831 : : end if
832 : :
833 : : ! Set up a counter that explore the relevant range
834 : : ! of points around the atom
835 : 540 : ishift=0
836 [ + + ]: 19320 : do ixp=igrid(mu)-irange(mu),igrid(mu)+irange(mu)
837 : 18600 : ishift=ishift+1
838 : 18600 : ii(ishift,mu)=1+mod(ngfft(mu)+mod(ixp,ngfft(mu)),ngfft(mu))
839 : 19140 : rrdiff(ishift,mu)=dble(ixp)/dble(ngfft(mu))-tau(mu)
840 : : end do
841 : :
842 : : ! End loop on mu
843 : : end do
844 : :
845 : : ! Conduct triple loop over restricted range of grid points for iatom
846 : :
847 [ + + ]: 6468 : do ishift3=1,1+2*irange(3)
848 : : ! map back to [1,ngfft(3)] for usual fortran index in unit cell
849 : 6192 : i3=ii(ishift3,3)
850 : : ! find vector from atom location to grid point (reduced)
851 : 6192 : rdiff3=rrdiff(ishift3,3)
852 : :
853 [ + + ]: 226506 : do ishift2=1,1+2*irange(2)
854 : 220128 : i2=ii(ishift2,2)
855 : 220128 : rdiff2=rrdiff(ishift2,2)
856 : : ! Prepare the computation of difmag2
857 : : difmag2_part=rmet(3,3)*rdiff3**2+rmet(2,2)*rdiff2**2&
858 : 220128 : & +2.0_dp*rmet(3,2)*rdiff3*rdiff2
859 : 220128 : difmag2_fact=2.0_dp*(rmet(3,1)*rdiff3+rmet(2,1)*rdiff2)
860 : :
861 [ + + ]: 8302464 : do ishift1=1,1+2*irange(1)
862 : 8076144 : rdiff1=rrdiff(ishift1,1)
863 : :
864 : : ! Compute (rgrid-tau-Rprim)**2
865 : 8076144 : difmag2= difmag2_part+rdiff1*(difmag2_fact+rmet(1,1)*rdiff1)
866 : :
867 : : ! Only accept contribution inside defined range
868 [ + + ]: 8296272 : if (difmag2<range2) then
869 : :
870 : : ! Prepare computation of core charge function and derivatives,
871 : : ! using splines
872 : 2641878 : difmag=sqrt(difmag2)
873 [ + + ]: 2641878 : if (difmag>=1.0d-10) then
874 : 2641764 : i1=ii(ishift1,1)
875 : 2641764 : yy=difmag*rangem1
876 : :
877 : : ! Compute index of yy over 1 to n1xccc scale
878 : 2641764 : jj=1+int(yy*(n1xccc-1))
879 : 2641764 : diff=yy-(jj-1)*delta
880 : :
881 : : ! Will evaluate spline fit (p. 86 Numerical Recipes, Press et al;
882 : : ! NOTE error in book for sign of "aa" term in derivative;
883 : : ! also see splfit routine).
884 : 2641764 : bb = diff*deltam1
885 : 2641764 : aa = 1.0_dp-bb
886 : 2641764 : cc = aa*(aa**2-1.0_dp)*delta2div6
887 : 2641764 : dd = bb*(bb**2-1.0_dp)*delta2div6
888 : :
889 : : ! Evaluate spline fit of 1st der of core charge density
890 : : ! from xccc1d(:,2,:) and (:,4,:)
891 : : func1=aa*xccc1d(jj,2,typat(iatom))+bb*xccc1d(jj+1,2,typat(iatom)) +&
892 : 2641764 : & cc*xccc1d(jj,4,typat(iatom))+dd*xccc1d(jj+1,4,typat(iatom))
893 : 2641764 : term1=func1*rangem1
894 : : ! Evaluate spline fit of 2nd der of core charge density
895 : : ! from xccc1d(:,3,:) and (:,5,:)
896 : : func2=aa*xccc1d(jj,3,typat(iatom))+bb*xccc1d(jj+1,3,typat(iatom)) +&
897 : 2641764 : & cc*xccc1d(jj,5,typat(iatom))+dd*xccc1d(jj+1,5,typat(iatom))
898 : 2641764 : term2=func2*rangem1**2
899 : :
900 : 2641764 : ifft=i1+n1*(i2-1+n2*(i3-1))
901 [ + + ]: 10567056 : tt(:)=rmet(:,1)*rdiff1+rmet(:,2)*rdiff2+rmet(:,3)*rdiff3
902 : :
903 : : ! Add contributions to 2nd derivative tensor
904 : : drss2=&
905 : : & (rdiff1*(drm(1,1,is2_in)*rdiff1+drm(1,2,is2_in)*rdiff2&
906 : : & +drm(1,3,is2_in)*rdiff3)&
907 : : & +rdiff2*(drm(2,1,is2_in)*rdiff1+drm(2,2,is2_in)*rdiff2&
908 : : & +drm(2,3,is2_in)*rdiff3)&
909 : : & +rdiff3*(drm(3,1,is2_in)*rdiff1+drm(3,2,is2_in)*rdiff2&
910 : 2641764 : & +drm(3,3,is2_in)*rdiff3))
911 : :
912 : : ! Loop over 1st strain index
913 [ + + ]: 18492348 : do is1=1,6
914 : :
915 : : drss1=&
916 : : & (rdiff1*(drm(1,1,is1)*rdiff1+drm(1,2,is1)*rdiff2&
917 : : & +drm(1,3,is1)*rdiff3)&
918 : : & +rdiff2*(drm(2,1,is1)*rdiff1+drm(2,2,is1)*rdiff2&
919 : : & +drm(2,3,is1)*rdiff3)&
920 : : & +rdiff3*(drm(3,1,is1)*rdiff1+drm(3,2,is1)*rdiff2&
921 : 15850584 : & +drm(3,3,is1)*rdiff3))
922 : :
923 : : d2rss=&
924 : : & (rdiff1*(d2rm(1,1,is1,is2_in)*rdiff1+d2rm(1,2,is1,is2_in)*rdiff2&
925 : : & +d2rm(1,3,is1,is2_in)*rdiff3)&
926 : : & +rdiff2*(d2rm(2,1,is1,is2_in)*rdiff1+d2rm(2,2,is1,is2_in)*rdiff2&
927 : : & +d2rm(2,3,is1,is2_in)*rdiff3)&
928 : : & +rdiff3*(d2rm(3,1,is1,is2_in)*rdiff1+d2rm(3,2,is1,is2_in)*rdiff2&
929 : 15850584 : & +d2rm(3,3,is1,is2_in)*rdiff3))
930 : :
931 : : ! Vall(0) X Rhocore(2) term
932 : : eltfrxc_core(is1,is2_in)=eltfrxc_core(is1,is2_in)+0.25_dp*&
933 : : & (vxc_core(ifft)*(term1*(d2rss/difmag&
934 : : & -drss1*drss2/difmag**3)&
935 : 15850584 : & +term2*drss1*drss2/difmag**2))
936 : :
937 : : ! Vall(1) X Rhocore(1) term
938 : : eltfrxc_core(is1,is2_in)=eltfrxc_core(is1,is2_in)+0.25_dp*&
939 : 15850584 : & vxc10_core(ifft)*drss1*term1/difmag
940 : : eltfrxc_core(is2_in,is1)=eltfrxc_core(is2_in,is1)+0.25_dp*&
941 : 18492348 : & vxc10_core(ifft)*drss1*term1/difmag
942 : :
943 : : ! End loop in is1
944 : : end do
945 : : ! Internal strain contributions
946 : : ts2(:)=drm(:,1,is2_in)*rdiff1+drm(:,2,is2_in)*rdiff2&
947 [ + + ]: 10567056 : & +drm(:,3,is2_in)*rdiff3
948 : :
949 : : eltfrxc_core(js:js+2,is2_in)=eltfrxc_core(js:js+2,is2_in)&
950 : : & -(vxc1is_core(ifft)*term1/difmag&
951 : : & +0.5_dp*vxc_core(ifft)*(term2-term1/difmag)*drss2/difmag**2)*tt(:)&
952 [ + + ]: 10567056 : & -(vxc_core(ifft)*term1/difmag)*ts2(:)
953 : :
954 : : ! End of the condition for the distance not to vanish
955 : : end if
956 : :
957 : : ! End of condition to be inside the range
958 : : end if
959 : :
960 : : ! End loop on ishift1
961 : : end do
962 : :
963 : : ! End loop on ishift2
964 : : end do
965 : :
966 : : ! End loop on ishift3
967 : : end do
968 : :
969 : : ! End loop on atoms
970 : : end do
971 : :
972 : : !In case of parallelism over atoms: communicate
973 [ - + ]: 96 : if (paral_atom) then
974 : 0 : call timab(48,1,tsec)
975 : 0 : call xmpi_sum(eltfrxc_core,my_comm_atom,ierr)
976 : 0 : call timab(48,2,tsec)
977 : : end if
978 : :
979 : : !Add core contribution to XC elastic tensor
980 [ + + + + ]: 7476 : eltfrxc(:,:)=eltfrxc(:,:)+eltfrxc_core(:,:)
981 : :
982 : : !Destroy atom table used for parallelism
983 : 96 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
984 : :
985 [ + - - + ]: 192 : ABI_FREE(d2rm)
986 : :
987 : : contains
988 : :
989 : 288 : function cross_elt(xx,yy,zz,aa,bb,cc)
990 : : !Define magnitude of cross product of two vectors
991 : : real(dp) :: cross_elt
992 : : real(dp),intent(in) :: xx,yy,zz,aa,bb,cc
993 : 288 : cross_elt=sqrt((yy*cc-zz*bb)**2+(zz*aa-xx*cc)**2+(xx*bb-yy*aa)**2)
994 : 288 : end function cross_elt
995 : :
996 : : end subroutine eltxccore
997 : : !!***
998 : :
999 : : !!****f* ABINIT/dfpt_eltfrloc
1000 : : !! NAME
1001 : : !! dfpt_eltfrloc
1002 : : !!
1003 : : !! FUNCTION
1004 : : !! Compute the frozen-wavefunction local pseudopotential contribution
1005 : : !! to the elastic tensor and the internal strain (derivative wrt one
1006 : : !! cartesian strain component and one reduced-coordinate atomic displacement).
1007 : : !!
1008 : : !! INPUTS
1009 : : !! atindx(natom)=index table for atoms (see gstate.f)
1010 : : !! gmet(3,3)=reciprocal space metric (bohr^-2)
1011 : : !! gprimd(3,3)=dimensional primitive translations for reciprocal space (bohr**-1)
1012 : : !! gsqcut=cutoff on G^2 based on ecut
1013 : : !! mgfft=maximum size of 1D FFTs
1014 : : !! mpi_enreg=information about MPI parallelization
1015 : : !! mqgrid=dimensioned number of q grid points for local psp spline
1016 : : !! natom=number of atoms in unit cell
1017 : : !! nattyp(ntypat)=number of atoms of each type
1018 : : !! nfft=(effective) number of FFT grid points (for this processor)
1019 : : !! ngfft(18)=contain all needed information about 3D FFT,
1020 : : !! see ~abinit/doc/variables/vargs.htm#ngfft
1021 : : !! ntypat=number of types of atoms
1022 : : !! ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information
1023 : : !! qgrid(mqgrid)=q point array for local psp spline fits
1024 : : !! rhog(2,nfft)=electron density in G space
1025 : : !! vlspl(mqgrid,2,ntypat)=q^2 v(q) spline for each type of atom.
1026 : : !!
1027 : : !! OUTPUT
1028 : : !! eltfrloc(6+3*natom,6)=non-symmetrized local pseudopotenial contribution
1029 : : !! to the elastic tensor and internal strain.
1030 : : !!
1031 : : !! SOURCE
1032 : :
1033 : 37 : subroutine dfpt_eltfrloc(atindx,eltfrloc,gmet,gprimd,gsqcut,mgfft,&
1034 : 37 : & mpi_enreg,mqgrid,natom,nattyp,nfft,ngfft,ntypat,ph1d,qgrid,rhog,vlspl)
1035 : :
1036 : : !Arguments ------------------------------------
1037 : : !scalars
1038 : : integer,intent(in) :: mgfft,mqgrid,natom,nfft,ntypat
1039 : : real(dp),intent(in) :: gsqcut
1040 : : type(MPI_type),intent(in) :: mpi_enreg
1041 : : !arrays
1042 : : integer,intent(in) :: atindx(natom),nattyp(ntypat),ngfft(18)
1043 : : real(dp),intent(in) :: gmet(3,3),gprimd(3,3),ph1d(2,3*(2*mgfft+1)*natom)
1044 : : real(dp),intent(in) :: qgrid(mqgrid),rhog(2,nfft),vlspl(mqgrid,2,ntypat)
1045 : : real(dp),intent(out) :: eltfrloc(6+3*natom,6)
1046 : :
1047 : : !Local variables-------------------------------
1048 : : !scalars
1049 : : integer,parameter :: im=2,re=1
1050 : : integer :: i1,i2,i3,ia,ia1,ia2,id1,id2,id3,ielt,ieltx,ierr,ig1,ig2,ig3,ii
1051 : : integer :: is1,is2,itypat,jj,ka,kb,kd,kg,me_fft,n1,n2,n3,nproc_fft
1052 : : real(dp),parameter :: tolfix=1.0000001_dp
1053 : : real(dp) :: aa,bb,cc,cutoff,d2g
1054 : : real(dp) :: dd,dg1,dg2,diff,dq
1055 : : real(dp) :: dq2div6,dqdiv6,dqm1,ee,ff,gmag,gsquar
1056 : : real(dp) :: sfi,sfr,term,term1
1057 : : !real(dp) :: ph1_elt,ph2_elt,ph3_elt,phi_elt,phr_elt
1058 : : real(dp) :: term2,term3,term4,term5,vion1,vion2,vion3
1059 : : !arrays
1060 : : integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/)
1061 : 37 : integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
1062 : 37 : integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
1063 : : real(dp) :: dgm(3,3,6),tsec(2)
1064 : 37 : real(dp),allocatable :: d2gm(:,:,:,:),elt_work(:,:)
1065 : :
1066 : : ! *************************************************************************
1067 : :
1068 : 37 : n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
1069 : 37 : me_fft=ngfft(11) ; nproc_fft=ngfft(10)
1070 : :
1071 : : !Get the distrib associated with this fft_grid
1072 : 37 : call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
1073 : :
1074 : : !-----
1075 : : !Compute 1st and 2nd derivatives of metric tensor wrt all strain components
1076 : : !and store for use in inner loop below.
1077 [ - + ]: 37 : ABI_MALLOC(d2gm,(3,3,6,6))
1078 : :
1079 : : !Loop over 2nd strain index
1080 [ + + ]: 259 : do is2=1,6
1081 : 222 : kg=idx(2*is2-1);kd=idx(2*is2)
1082 [ + + ]: 888 : do jj = 1,3
1083 [ + + ]: 2886 : dgm(:,jj,is2)=-(gprimd(kg,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(kg,jj))
1084 : : end do
1085 : : ! Loop over 1st strain index, upper triangle only
1086 [ + + ]: 1036 : do is1=1,is2
1087 : 777 : ka=idx(2*is1-1);kb=idx(2*is1)
1088 [ + + + + ]: 10101 : d2gm(:,:,is1,is2)=0._dp
1089 [ + + ]: 3330 : do jj = 1,3
1090 [ + + ]: 2331 : if(ka==kg) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
1091 [ + + ]: 4440 : & +gprimd(kb,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(kb,jj)
1092 [ + + ]: 2331 : if(ka==kd) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
1093 [ + + ]: 2664 : & +gprimd(kb,:)*gprimd(kg,jj)+gprimd(kg,:)*gprimd(kb,jj)
1094 [ + + ]: 2331 : if(kb==kg) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
1095 [ + + ]: 3108 : & +gprimd(ka,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(ka,jj)
1096 [ + + ]: 2331 : if(kb==kd) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
1097 [ + + ]: 5217 : & +gprimd(ka,:)*gprimd(kg,jj)+gprimd(kg,:)*gprimd(ka,jj)
1098 : : end do
1099 : : end do !is1
1100 : : end do !is2
1101 : :
1102 : : !Zero out array to permit accumulation over atom types below:
1103 [ + + + + ]: 3283 : eltfrloc(:,:)=0.0_dp
1104 : :
1105 : 37 : dq=(qgrid(mqgrid)-qgrid(1))/dble(mqgrid-1)
1106 : 37 : dqm1=1.0_dp/dq
1107 : 37 : dqdiv6=dq/6.0_dp
1108 : 37 : dq2div6=dq**2/6.0_dp
1109 : 37 : cutoff=gsqcut*tolfix
1110 : 37 : id1=n1/2+2
1111 : 37 : id2=n2/2+2
1112 : 37 : id3=n3/2+2
1113 : :
1114 : 37 : ia1=1
1115 [ + + ]: 98 : do itypat=1,ntypat
1116 : : ! ia1,ia2 sets range of loop over atoms:
1117 : 61 : ia2=ia1+nattyp(itypat)-1
1118 : 61 : ii=0
1119 [ + + ]: 1411 : do i3=1,n3
1120 : 1350 : ig3=i3-(i3/id3)*n3-1
1121 [ + + ]: 24317 : do i2=1,n2
1122 [ + - ]: 24256 : if (fftn2_distrib(i2)==me_fft) then
1123 : 22906 : ig2=i2-(i2/id2)*n2-1
1124 [ + + ]: 419836 : do i1=1,n1
1125 : 396930 : ig1=i1-(i1/id1)*n1-1
1126 : :
1127 : 396930 : ii=ii+1
1128 : : ! Skip G=0:
1129 [ + + ]: 396930 : if (ig1==0 .and. ig2==0 .and. ig3==0) cycle
1130 : :
1131 : : ! Skip G**2 outside cutoff:
1132 : 396869 : gsquar=gsq_elt(ig1,ig2,ig3)
1133 [ + + ]: 419775 : if (gsquar<=cutoff) then
1134 : 132058 : gmag=sqrt(gsquar)
1135 : :
1136 : : ! Compute vion(G) for given type of atom
1137 : 132058 : jj=1+int(gmag*dqm1)
1138 : 132058 : diff=gmag-qgrid(jj)
1139 : :
1140 : : ! Evaluate spline fit from q^2 V(q) to get V(q):
1141 : : ! (p. 86 Numerical Recipes, Press et al; NOTE error in book for sign
1142 : : ! of "aa" term in derivative; also see splfit routine).
1143 : 132058 : bb = diff*dqm1
1144 : 132058 : aa = 1.0_dp-bb
1145 : 132058 : cc = aa*(aa**2-1.0_dp)*dq2div6
1146 : 132058 : dd = bb*(bb**2-1.0_dp)*dq2div6
1147 : : term1 = (aa*vlspl(jj,1,itypat)+bb*vlspl(jj+1,1,itypat) +&
1148 : 132058 : & cc*vlspl(jj,2,itypat)+dd*vlspl(jj+1,2,itypat))
1149 : 132058 : vion1=term1 / gsquar
1150 : :
1151 : : ! Also get dV(q)/dq:
1152 : : ! (note correction of Numerical Recipes sign error
1153 : : ! before (3._dp*aa**2-1._dp)
1154 : 132058 : ee= vlspl(jj+1,1,itypat)-vlspl(jj,1,itypat)
1155 : : ff= (3._dp*bb**2-1._dp)*vlspl(jj+1,2,itypat) &
1156 : 132058 : & - (3._dp*aa**2-1._dp)*vlspl(jj,2,itypat)
1157 : 132058 : term2 = ee*dqm1 + ff*dqdiv6
1158 : 132058 : vion2 = term2/gsquar - 2._dp*term1/(gsquar*gmag)
1159 : :
1160 : : ! Also get V''(q)
1161 : 132058 : term3=aa*vlspl(jj,2,itypat)+bb*vlspl(jj+1,2,itypat)
1162 : 132058 : vion3 = (term3 - 4.0_dp*term2/gmag + 6._dp*term1/gsquar)/gsquar
1163 : :
1164 : : ! Assemble structure factor over all atoms of given type:
1165 : 132058 : sfr=zero;sfi=zero
1166 [ + + ]: 361046 : do ia=ia1,ia2
1167 : 228988 : sfr=sfr+phre_elt(ig1,ig2,ig3,ia)
1168 : 361046 : sfi=sfi-phimag_elt(ig1,ig2,ig3,ia)
1169 : : end do
1170 : 132058 : term=(rhog(re,ii)*sfr+rhog(im,ii)*sfi)
1171 : :
1172 : : ! Loop over 2nd strain index
1173 [ + + ]: 924406 : do is2=1,6
1174 : 792348 : dg2=0.5_dp*dgsqds_elt(ig1,ig2,ig3,is2)/gmag
1175 : : ! Loop over 1st strain index, upper triangle only
1176 [ + + ]: 3565566 : do is1=1,is2
1177 : 2773218 : dg1=0.5_dp*dgsqds_elt(ig1,ig2,ig3,is1)/gmag
1178 : 2773218 : d2g=(0.25_dp*d2gsqds_elt(ig1,ig2,ig3,is1,is2)-dg1*dg2)/gmag
1179 : : eltfrloc(is1,is2)=eltfrloc(is1,is2)+&
1180 : 2773218 : & term*(vion3*dg1*dg2+vion2*d2g)
1181 [ + + ]: 2773218 : if(is2<=3)&
1182 : 792348 : & eltfrloc(is1,is2)=eltfrloc(is1,is2)-term*vion2*dg1
1183 [ + + ]: 2773218 : if(is1<=3)&
1184 : 1980870 : & eltfrloc(is1,is2)=eltfrloc(is1,is2)-term*vion2*dg2
1185 [ + + + + ]: 2773218 : if(is1<=3 .and. is2<=3)&
1186 : 1584696 : & eltfrloc(is1,is2)=eltfrloc(is1,is2)+term*vion1
1187 : : end do !is1
1188 : :
1189 : : ! Internal strain section - loop over current atoms
1190 [ + + ]: 2298334 : do ia=ia1,ia2
1191 [ + + ]: 1373928 : if(is2 <=3) then
1192 : 686964 : term4=vion2*dg2-vion1
1193 : : else
1194 : 686964 : term4=vion2*dg2
1195 : : end if
1196 : : term5=-two_pi*(rhog(re,ii)*phimag_elt(ig1,ig2,ig3,ia)&
1197 : 1373928 : & +rhog(im,ii)*phre_elt(ig1,ig2,ig3,ia))*term4
1198 : 1373928 : eltfrloc(7+3*(ia-1),is2)=eltfrloc(7+3*(ia-1),is2)+term5*dble(ig1)
1199 : 1373928 : eltfrloc(8+3*(ia-1),is2)=eltfrloc(8+3*(ia-1),is2)+term5*dble(ig2)
1200 : 2166276 : eltfrloc(9+3*(ia-1),is2)=eltfrloc(9+3*(ia-1),is2)+term5*dble(ig3)
1201 : : end do
1202 : :
1203 : : end do !is2
1204 : :
1205 : : ! End skip G**2 outside cutoff:
1206 : : end if
1207 : :
1208 : : ! End loop on n1, n2, n3. There is a "cycle" inside the loop
1209 : : end do
1210 : : end if
1211 : : end do
1212 : : end do
1213 : :
1214 : : ! End loop on type of atoms
1215 : 98 : ia1=ia2+1
1216 : : end do
1217 : : !Init mpi_comm
1218 : 37 : call timab(48,1,tsec)
1219 : 37 : call xmpi_sum(eltfrloc,mpi_enreg%comm_fft,ierr)
1220 : 37 : call timab(48,2,tsec)
1221 : :
1222 : : !Fill in lower triangle
1223 [ + + ]: 222 : do is2=2,6
1224 [ + + ]: 777 : do is1=1,is2-1
1225 : 740 : eltfrloc(is2,is1)=eltfrloc(is1,is2)
1226 : : end do
1227 : : end do
1228 : :
1229 : : !The indexing array atindx is used to reestablish the correct
1230 : : !order of atoms
1231 [ + - - + : 185 : ABI_MALLOC(elt_work,(6+3*natom,6))
+ - + - +
- + - - +
- + ]
1232 [ + + + + ]: 1591 : elt_work(1:6,1:6)=eltfrloc(1:6,1:6)
1233 [ + + ]: 131 : do ia=1,natom
1234 : 94 : ielt=7+3*(ia-1)
1235 : 94 : ieltx=7+3*(atindx(ia)-1)
1236 [ + + + + ]: 2387 : elt_work(ielt:ielt+2,1:6)=eltfrloc(ieltx:ieltx+2,1:6)
1237 : : end do
1238 [ + + + + ]: 3283 : eltfrloc(:,:)=elt_work(:,:)
1239 : :
1240 [ - + ]: 37 : ABI_FREE(d2gm)
1241 [ - + ]: 74 : ABI_FREE(elt_work)
1242 : :
1243 : : contains
1244 : :
1245 : : !Real and imaginary parts of phase.
1246 : 1602916 : function phr_elt(x1,y1,x2,y2,x3,y3)
1247 : :
1248 : : real(dp) :: phr_elt,x1,x2,x3,y1,y2,y3
1249 : 1602916 : phr_elt=(x1*x2-y1*y2)*x3-(y1*x2+x1*y2)*y3
1250 : : end function phr_elt
1251 : :
1252 : 1602916 : function phi_elt(x1,y1,x2,y2,x3,y3)
1253 : :
1254 : : real(dp):: phi_elt,x1,x2,x3,y1,y2,y3
1255 : 1602916 : phi_elt=(x1*x2-y1*y2)*y3+(y1*x2+x1*y2)*x3
1256 : : end function phi_elt
1257 : :
1258 : 6411664 : function ph1_elt(nri,ig1,ia)
1259 : :
1260 : : real(dp):: ph1_elt
1261 : : integer :: nri,ig1,ia
1262 : 6411664 : ph1_elt=ph1d(nri,ig1+1+n1+(ia-1)*(2*n1+1))
1263 : 6411664 : end function ph1_elt
1264 : :
1265 : 6411664 : function ph2_elt(nri,ig2,ia)
1266 : :
1267 : : real(dp):: ph2_elt
1268 : : integer :: nri,ig2,ia
1269 : 6411664 : ph2_elt=ph1d(nri,ig2+1+n2+(ia-1)*(2*n2+1)+natom*(2*n1+1))
1270 : 6411664 : end function ph2_elt
1271 : :
1272 : 6411664 : function ph3_elt(nri,ig3,ia)
1273 : :
1274 : : real(dp):: ph3_elt
1275 : : integer :: nri,ig3,ia
1276 : 6411664 : ph3_elt=ph1d(nri,ig3+1+n3+(ia-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1))
1277 : 6411664 : end function ph3_elt
1278 : :
1279 : 1602916 : function phre_elt(ig1,ig2,ig3,ia)
1280 : :
1281 : : real(dp):: phre_elt
1282 : : integer :: ig1,ig2,ig3,ia
1283 : : phre_elt=phr_elt(ph1_elt(re,ig1,ia),ph1_elt(im,ig1,ia),&
1284 : 1602916 : & ph2_elt(re,ig2,ia),ph2_elt(im,ig2,ia),ph3_elt(re,ig3,ia),ph3_elt(im,ig3,ia))
1285 : 1602916 : end function phre_elt
1286 : :
1287 : 1602916 : function phimag_elt(ig1,ig2,ig3,ia)
1288 : :
1289 : : real(dp) :: phimag_elt
1290 : : integer :: ig1,ig2,ig3,ia
1291 : : phimag_elt=phi_elt(ph1_elt(re,ig1,ia),ph1_elt(im,ig1,ia),&
1292 : 1602916 : & ph2_elt(re,ig2,ia),ph2_elt(im,ig2,ia),ph3_elt(re,ig3,ia),ph3_elt(im,ig3,ia))
1293 : 1602916 : end function phimag_elt
1294 : :
1295 : 396869 : function gsq_elt(i1,i2,i3)
1296 : :
1297 : : real(dp) :: gsq_elt
1298 : : integer :: i1,i2,i3
1299 : : !Define G^2 based on G space metric gmet.
1300 : : gsq_elt=dble(i1*i1)*gmet(1,1)+dble(i2*i2)*gmet(2,2)+&
1301 : : & dble(i3*i3)*gmet(3,3)+dble(2*i1*i2)*gmet(1,2)+&
1302 : 396869 : & dble(2*i2*i3)*gmet(2,3)+dble(2*i3*i1)*gmet(3,1)
1303 : 396869 : end function gsq_elt
1304 : :
1305 : 3565566 : function dgsqds_elt(i1,i2,i3,is)
1306 : :
1307 : : real(dp) :: dgsqds_elt
1308 : : integer :: i1,i2,i3,is
1309 : : !Define dG^2/ds based on G space metric derivative
1310 : : dgsqds_elt=dble(i1*i1)*dgm(1,1,is)+dble(i2*i2)*dgm(2,2,is)+&
1311 : : & dble(i3*i3)*dgm(3,3,is)+&
1312 : : & dble(i1*i2)*(dgm(1,2,is)+dgm(2,1,is))+&
1313 : : & dble(i1*i3)*(dgm(1,3,is)+dgm(3,1,is))+&
1314 : 3565566 : & dble(i2*i3)*(dgm(2,3,is)+dgm(3,2,is))
1315 : 3565566 : end function dgsqds_elt
1316 : :
1317 : 2773218 : function d2gsqds_elt(i1,i2,i3,is1,is2)
1318 : :
1319 : : real(dp) :: d2gsqds_elt
1320 : : integer :: i1,i2,i3,is1,is2
1321 : : !Define 2dG^2/ds1ds2 based on G space metric derivative
1322 : : d2gsqds_elt=dble(i1*i1)*d2gm(1,1,is1,is2)+&
1323 : : & dble(i2*i2)*d2gm(2,2,is1,is2)+dble(i3*i3)*d2gm(3,3,is1,is2)+&
1324 : : & dble(i1*i2)*(d2gm(1,2,is1,is2)+d2gm(2,1,is1,is2))+&
1325 : : & dble(i1*i3)*(d2gm(1,3,is1,is2)+d2gm(3,1,is1,is2))+&
1326 : 2773218 : & dble(i2*i3)*(d2gm(2,3,is1,is2)+d2gm(3,2,is1,is2))
1327 : 2773218 : end function d2gsqds_elt
1328 : :
1329 : : end subroutine dfpt_eltfrloc
1330 : : !!***
1331 : :
1332 : : !!****f* ABINIT/dfpt_eltfrkin
1333 : : !! NAME
1334 : : !! dfpt_eltfrkin
1335 : : !!
1336 : : !! FUNCTION
1337 : : !! Compute the frozen-wavefunction kinetic enegy contribution to the
1338 : : !! elastic tensor
1339 : : !!
1340 : : !! INPUTS
1341 : : !! cg(2,mpw*nspinor*mband_mem*mkmem*nsppol)=<G|Cnk>=Fourier coefficients of wavefunction
1342 : : !! ecut=cut-off energy for plane wave basis sphere (Ha)
1343 : : !! ecutsm=smearing energy for plane wave kinetic energy (Ha) (NOT NEEDED !)
1344 : : !! effmass_free=effective mass for electrons (1. in common case)
1345 : : !! istwfk(nkpt)=input option parameter that describes the storage of wfs
1346 : : !! kg(3,mpw*mkmem)=work array for coordinates of G vectors in basis
1347 : : !! kptns(3,nkpt)=coordinates of k points in terms of reciprocal space
1348 : : !! primitive translations
1349 : : !! mband=maximum number of bands
1350 : : !! mband_mem=maximum number of bands in memory
1351 : : !! mgfft=maximum size of 1D FFTs
1352 : : !! mkmem=number of k points treated by this node.
1353 : : !! mpi_enreg=information about MPI parallelization
1354 : : !! mpw=maximum dimension for number of planewaves
1355 : : !! nband(nkpt*nsppol)=number of bands being considered per k point
1356 : : !! nkpt=number of k points
1357 : : !! ngfft(18)=contain all needed information about 3D FFT, i
1358 : : !! see ~abinit/doc/variables/vargs.htm#ngfft
1359 : : !! npwarr(nkpt)=number of planewaves at each k point, and boundary
1360 : : !! nspinor=number of spinorial components of the wavefunctions
1361 : : !! nsppol=1 for unpolarized, 2 for polarized
1362 : : !! occ(mband*nkpt*nsppol)=occupation numbers of bands (usually 2)
1363 : : !! at each k point
1364 : : !! rprimd(3,3)=dimensional real space primitive translations (bohr)
1365 : : !! wtk(nkpt)=k point weights
1366 : : !!
1367 : : !! OUTPUT
1368 : : !! eltfrkin(6,6)=non-symmetrized kinetic energy contribution to the
1369 : : !! elastic tensor
1370 : : !!
1371 : : !! SOURCE
1372 : :
1373 : 37 : subroutine dfpt_eltfrkin(cg,eltfrkin,ecut,ecutsm,effmass_free,&
1374 : 37 : & istwfk,kg,kptns,mband,mband_mem,mgfft,mkmem,mpi_enreg,&
1375 : 37 : & mpw,nband,nkpt,ngfft,npwarr,nspinor,nsppol,occ,rprimd,wtk)
1376 : :
1377 : : !Arguments ------------------------------------
1378 : : !scalars
1379 : : integer,intent(in) :: mband,mband_mem,mgfft,mkmem,mpw,nkpt,nspinor,nsppol
1380 : : real(dp),intent(in) :: ecut,ecutsm,effmass_free
1381 : : type(MPI_type),intent(in) :: mpi_enreg
1382 : : !arrays
1383 : : integer,intent(in) :: istwfk(nkpt),kg(3,mpw*mkmem),nband(nkpt*nsppol)
1384 : : integer,intent(in) :: ngfft(18),npwarr(nkpt)
1385 : : real(dp),intent(in) :: cg(2,mpw*nspinor*mband_mem*mkmem*nsppol),kptns(3,nkpt)
1386 : : real(dp),intent(in) :: occ(mband*nkpt*nsppol),rprimd(3,3),wtk(nkpt)
1387 : : real(dp),intent(out) :: eltfrkin(6,6)
1388 : :
1389 : : !Local variables-------------------------------
1390 : : !scalars
1391 : : integer :: bdtot_index,iband,icg,ierr,ii,ikg
1392 : : integer :: ikpt,index,ipw,isppol,istwf_k,jj,master,me,n1,n2
1393 : : integer :: n3,nband_k,nkinout,npw_k,spaceComm
1394 : : integer :: nband_me, iband_me
1395 : : real(dp) :: ucvol
1396 : : !arrays
1397 : 37 : integer,allocatable :: gbound(:,:),kg_k(:,:)
1398 : : real(dp) :: gmet(3,3),gprimd(3,3),kpoint(3),rmet(3,3),tsec(2)
1399 : 37 : real(dp),allocatable :: cwavef(:,:),ekinout(:)
1400 : 37 : real(dp),allocatable :: eltfrkink(:,:)
1401 : :
1402 : : ! *************************************************************************
1403 : :
1404 : : DBG_ENTER("COLL")
1405 : :
1406 : : !Default for sequential use
1407 : 37 : master=0
1408 : : !Init mpi_comm
1409 : 37 : spaceComm=mpi_enreg%comm_cell
1410 : 37 : me=mpi_enreg%me_kpt
1411 : :
1412 : : !Compute gmet, gprimd and ucvol from rprimd
1413 : 37 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
1414 : :
1415 : 37 : eltfrkin(:,:)=0.0_dp
1416 : 37 : bdtot_index=0
1417 : 37 : icg=0
1418 : :
1419 : 37 : n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
1420 [ - + + - : 111 : ABI_MALLOC(kg_k,(3,mpw))
+ - + - -
+ - + ]
1421 [ - + + - : 111 : ABI_MALLOC(cwavef,(2,mpw*nspinor))
+ - + - -
+ - + ]
1422 [ - + ]: 37 : ABI_MALLOC(eltfrkink,(6,6))
1423 : :
1424 : : !Define k-points distribution
1425 : :
1426 : : !LOOP OVER SPINS
1427 [ + + ]: 76 : do isppol=1,nsppol
1428 : 39 : ikg=0
1429 : :
1430 : : ! Loop over k points
1431 [ + + ]: 5756 : do ikpt=1,nkpt
1432 : :
1433 : 5680 : nband_k=nband(ikpt+(isppol-1)*nkpt)
1434 : 5680 : istwf_k=istwfk(ikpt)
1435 : 5680 : npw_k=npwarr(ikpt)
1436 : :
1437 : : ! Skip this k-point if not the proper processor
1438 [ + + ]: 5680 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then
1439 : 32 : bdtot_index=bdtot_index+nband_k
1440 : 32 : cycle
1441 : : end if
1442 : :
1443 : : ! find number of bands I will actually treat
1444 : 5648 : nband_me = proc_distrb_nband(mpi_enreg%proc_distrb,ikpt,nband_k,isppol,me)
1445 : :
1446 [ + - - + : 22592 : ABI_MALLOC(gbound,(2*mgfft+8,2))
+ - + - +
- + - - +
- + ]
1447 [ + + ]: 22592 : kpoint(:)=kptns(:,ikpt)
1448 : :
1449 [ + + + + ]: 2776608 : kg_k(:,:) = 0
1450 : :
1451 : :
1452 : : !$OMP PARALLEL DO PRIVATE(ipw) SHARED(ikg,kg,kg_k,npw_k)
1453 [ + + ]: 662007 : do ipw=1,npw_k
1454 : 656359 : kg_k(1,ipw)=kg(1,ipw+ikg)
1455 : 656359 : kg_k(2,ipw)=kg(2,ipw+ikg)
1456 : 662007 : kg_k(3,ipw)=kg(3,ipw+ikg)
1457 : : end do
1458 : :
1459 : 5648 : call sphereboundary(gbound,istwf_k,kg_k,mgfft,npw_k)
1460 : :
1461 : 5648 : index=1+icg
1462 : :
1463 [ + + + + ]: 242864 : eltfrkink(:,:)=0.0_dp
1464 : :
1465 : 5648 : nkinout=6*6
1466 [ - + ]: 5648 : ABI_MALLOC(ekinout,(nkinout))
1467 [ + + ]: 208976 : ekinout(:)=zero
1468 : :
1469 : : iband_me = 0
1470 [ + + ]: 29080 : do iband=1,nband_k
1471 : :
1472 [ - + ]: 23432 : if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= me) cycle
1473 : 23432 : iband_me = iband_me + 1
1474 : :
1475 [ + + + + ]: 8715776 : cwavef(:,1:npw_k*nspinor)=cg(:,1+(iband_me-1)*npw_k*nspinor+icg:iband_me*npw_k*nspinor+icg)
1476 : :
1477 : : call d2kindstr2(cwavef,ecut,ecutsm,effmass_free,ekinout,gmet,gprimd,&
1478 : 23432 : & istwf_k,kg_k,kpoint,npw_k,nspinor)
1479 : :
1480 [ + + + + ]: 1013224 : eltfrkink(:,:)=eltfrkink(:,:)+ occ(iband+bdtot_index)* reshape(ekinout(:), (/6,6/) )
1481 : :
1482 : : end do !iband
1483 : :
1484 [ - + ]: 5648 : ABI_FREE(ekinout)
1485 : :
1486 [ + + + + ]: 242864 : eltfrkin(:,:)=eltfrkin(:,:)+wtk(ikpt)*eltfrkink(:,:)
1487 : :
1488 [ - + ]: 5648 : ABI_FREE(gbound)
1489 : :
1490 : 5648 : bdtot_index=bdtot_index+nband_k
1491 : :
1492 [ + - ]: 5687 : if (mkmem/=0) then
1493 : : ! Handle case in which kg, cg, are kept in core
1494 : 5648 : icg=icg+npw_k*nspinor*nband_me
1495 : 5648 : ikg=ikg+npw_k
1496 : : end if
1497 : :
1498 : : end do
1499 : : end do ! End loops on isppol and ikpt
1500 : :
1501 : : !Fill in lower triangle
1502 [ + + ]: 222 : do jj=2,6
1503 [ + + ]: 777 : do ii=1,jj-1
1504 : 740 : eltfrkin(jj,ii)=eltfrkin(ii,jj)
1505 : : end do
1506 : : end do
1507 : :
1508 : : !Accumulate eltfrkin on all proc.
1509 : 37 : call timab(48,1,tsec)
1510 : 37 : call xmpi_sum(eltfrkin,spaceComm,ierr)
1511 : 37 : call timab(48,2,tsec)
1512 : :
1513 [ - + ]: 37 : ABI_FREE(cwavef)
1514 [ - + ]: 37 : ABI_FREE(eltfrkink)
1515 [ - + ]: 37 : ABI_FREE(kg_k)
1516 : :
1517 : : DBG_EXIT("COLL")
1518 : :
1519 : : contains
1520 : : !!***
1521 : :
1522 : : !!****f* ABINIT/d2kindstr2
1523 : : !! NAME
1524 : : !! d2kindstr2
1525 : : !!
1526 : : !! FUNCTION
1527 : : !! compute expectation value of the second derivatives of the kinetic energy
1528 : : !! wrt strain for one band and kpoint
1529 : : !!
1530 : : !! INPUTS
1531 : : !! cwavef(2,npw*nspinor)=wavefunction for current band
1532 : : !! ecut=cut-off energy for plane wave basis sphere (Ha)
1533 : : !! ecutsm=smearing energy for plane wave kinetic energy (Ha)
1534 : : !! effmass_free=effective mass for electrons (1. in common case)
1535 : : !! gmet(3,3)=reciprocal lattice metric tensor ($\textrm{Bohr}^{-2}$)
1536 : : !! gprimd(3,3)=primitive vectors in reciprocal space
1537 : : !! istwfk=information about wavefunction storage
1538 : : !! kg_k(3,npw)=integer coordinates of planewaves in basis sphere.
1539 : : !! kpt(3)=reduced coordinates of k point
1540 : : !! npw=number of plane waves at kpt.
1541 : : !! nspinor=number of spinorial components of the wavefunction
1542 : : !!
1543 : : !! OUTPUT
1544 : : !! ekinout(36)=expectation values of the second strain derivatives
1545 : : !! of the (modified) kinetic energy
1546 : : !!
1547 : : !! NOTES
1548 : : !! Usually, the kinetic energy expression is $(1/2) (2 \pi)^2 (k+G)^2 $
1549 : : !! However, the present implementation allows for a modification
1550 : : !! of this kinetic energy, in order to obtain smooth total energy
1551 : : !! curves with respect to the cut-off energy or the cell size and shape.
1552 : : !! Thus the usual expression is kept if it is lower then ecut-ecutsm,
1553 : : !! zero is returned beyond ecut, and in between, the kinetic
1554 : : !! energy is DIVIDED by a smearing factor (to make it infinite at the
1555 : : !! cut-off energy). The smearing factor is $x^2 (3-2x)$, where
1556 : : !! x = (ecut- unmodified energy)/ecutsm.
1557 : : !!
1558 : : !! SOURCE
1559 : :
1560 : 23432 : subroutine d2kindstr2(cwavef,ecut,ecutsm,effmass_free,ekinout,gmet,gprimd,&
1561 : 23432 : & istwfk,kg_k,kpt,npw,nspinor)
1562 : :
1563 : : !Arguments ------------------------------------
1564 : : !scalars
1565 : : integer,intent(in) :: istwfk,npw,nspinor
1566 : : real(dp),intent(in) :: ecut,ecutsm,effmass_free
1567 : : !arrays
1568 : : integer,intent(in) :: kg_k(3,npw)
1569 : : real(dp),intent(in) :: cwavef(2,npw*nspinor),gmet(3,3),gprimd(3,3),kpt(3)
1570 : : real(dp),intent(inout) :: ekinout(36) !vz_i
1571 : :
1572 : : !Local variables-------------------------------
1573 : : !scalars
1574 : : integer,parameter :: im=2,re=1
1575 : : integer :: ig,igs,ii,ispinor,istr1,istr2,ka,kb,kd,kg
1576 : : real(dp) :: d2fkin,d2fsm,d2kinacc,d2kpg2,dfkin,dfsm,dkpg21,dkpg22,ecutsm_inv
1577 : : real(dp) :: fsm,gpk1,gpk2,gpk3,htpisq,kpg2,term,xx
1578 : : !arrays
1579 : : integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/)
1580 : : real(dp) :: d2gm(3,3),dgm01(3,3),dgm10(3,3)
1581 : :
1582 : : ! *************************************************************************
1583 : : !
1584 : : !htpisq is (1/2) (2 Pi) **2:
1585 : 23432 : htpisq=0.5_dp*(two_pi)**2
1586 : :
1587 : 23432 : ecutsm_inv=0.0_dp
1588 [ + + ]: 23432 : if(ecutsm>1.0d-20)ecutsm_inv=1/ecutsm
1589 : :
1590 : : !Loop over 2nd strain index
1591 [ + + ]: 164024 : do istr2=1,6
1592 : : ! Loop over 1st strain index, upper triangle only
1593 [ + + ]: 656096 : do istr1=1,istr2
1594 : :
1595 : 492072 : ka=idx(2*istr1-1);kb=idx(2*istr1);kg=idx(2*istr2-1);kd=idx(2*istr2)
1596 : :
1597 [ + + ]: 1968288 : do ii = 1,3
1598 [ + + ]: 5904864 : dgm01(:,ii)=-(gprimd(ka,:)*gprimd(kb,ii)+gprimd(kb,:)*gprimd(ka,ii))
1599 [ + + ]: 6396936 : dgm10(:,ii)=-(gprimd(kg,:)*gprimd(kd,ii)+gprimd(kd,:)*gprimd(kg,ii))
1600 : : end do
1601 : :
1602 : 492072 : d2gm(:,:)=0._dp
1603 [ + + ]: 1968288 : do ii = 1,3
1604 [ + + ]: 1476216 : if(ka==kg) d2gm(:,ii)=d2gm(:,ii)&
1605 [ + + ]: 2811840 : & +gprimd(kb,:)*gprimd(kd,ii)+gprimd(kd,:)*gprimd(kb,ii)
1606 [ + + ]: 1476216 : if(ka==kd) d2gm(:,ii)=d2gm(:,ii)&
1607 [ + + ]: 1687104 : & +gprimd(kb,:)*gprimd(kg,ii)+gprimd(kg,:)*gprimd(kb,ii)
1608 [ + + ]: 1476216 : if(kb==kg) d2gm(:,ii)=d2gm(:,ii)&
1609 [ + + ]: 1968288 : & +gprimd(ka,:)*gprimd(kd,ii)+gprimd(kd,:)*gprimd(ka,ii)
1610 [ + + ]: 1476216 : if(kb==kd) d2gm(:,ii)=d2gm(:,ii)&
1611 [ + + ]: 3303912 : & +gprimd(ka,:)*gprimd(kg,ii)+gprimd(kg,:)*gprimd(ka,ii)
1612 : : end do
1613 [ + + + + ]: 6396936 : d2gm(:,:)=0.5_dp*d2gm(:,:)
1614 : :
1615 : : d2kinacc=0._dp
1616 : :
1617 : : ! loop on spinor index
1618 [ + + ]: 984144 : do ispinor=1,nspinor
1619 : 492072 : igs=(ispinor-1)*npw
1620 : : ! loop on plane waves
1621 [ + + ]: 61830552 : do ig=1,npw
1622 : 60846408 : gpk1=dble(kg_k(1,ig))+kpt(1)
1623 : 60846408 : gpk2=dble(kg_k(2,ig))+kpt(2)
1624 : 60846408 : gpk3=dble(kg_k(3,ig))+kpt(3)
1625 : : kpg2=htpisq*&
1626 : : & ( gmet(1,1)*gpk1**2+ &
1627 : : & gmet(2,2)*gpk2**2+ &
1628 : : & gmet(3,3)*gpk3**2 &
1629 : : & +2.0_dp*(gpk1*gmet(1,2)*gpk2+ &
1630 : : & gpk1*gmet(1,3)*gpk3+ &
1631 : 60846408 : & gpk2*gmet(2,3)*gpk3 ) )
1632 : : dkpg21=htpisq*&
1633 : : & ( dgm01(1,1)*gpk1**2+ &
1634 : : & dgm01(2,2)*gpk2**2+ &
1635 : : & dgm01(3,3)*gpk3**2 &
1636 : : & +2.0_dp*(gpk1*dgm01(1,2)*gpk2+ &
1637 : : & gpk1*dgm01(1,3)*gpk3+ &
1638 : 60846408 : & gpk2*dgm01(2,3)*gpk3 ) )
1639 : : dkpg22=htpisq*&
1640 : : & ( dgm10(1,1)*gpk1**2+ &
1641 : : & dgm10(2,2)*gpk2**2+ &
1642 : : & dgm10(3,3)*gpk3**2 &
1643 : : & +2.0_dp*(gpk1*dgm10(1,2)*gpk2+ &
1644 : : & gpk1*dgm10(1,3)*gpk3+ &
1645 : 60846408 : & gpk2*dgm10(2,3)*gpk3 ) )
1646 : : d2kpg2=htpisq*&
1647 : : & ( d2gm(1,1)*gpk1**2+ &
1648 : : & d2gm(2,2)*gpk2**2+ &
1649 : : & d2gm(3,3)*gpk3**2 &
1650 : : & +2.0_dp*(gpk1*d2gm(1,2)*gpk2+ &
1651 : : & gpk1*d2gm(1,3)*gpk3+ &
1652 : 60846408 : & gpk2*d2gm(2,3)*gpk3 ) )
1653 : :
1654 [ + - ]: 60846408 : if(kpg2>ecut-tol12)then
1655 : : dfkin=0._dp
1656 : : d2fkin=0._dp
1657 [ + + ]: 60846408 : elseif(kpg2>ecut-ecutsm)then
1658 : : ! This kinetic cutoff smoothing function and its xx derivatives
1659 : : ! were produced with Mathematica and the fortran code has been
1660 : : ! numerically checked against Mathematica.
1661 : 6141576 : xx=(ecut-kpg2)*ecutsm_inv
1662 : 6141576 : fsm=1.0_dp/(xx**2*(3+xx*(1+xx*(-6+3*xx))))
1663 : 6141576 : dfsm=-3.0_dp*(-1+xx)**2*xx*(2+5*xx)*fsm**2
1664 : : d2fsm=6.0_dp*xx**2*(9+xx*(8+xx*(-52+xx*(-3+xx*(137+xx*&
1665 : 6141576 : & (-144+45*xx))))))*fsm**3
1666 : 6141576 : dfkin=fsm-ecutsm_inv*kpg2*dfsm
1667 : 6141576 : d2fkin=ecutsm_inv*(-2.0_dp*dfsm+ecutsm_inv*kpg2*d2fsm)
1668 : : else
1669 : : dfkin=1._dp
1670 : : d2fkin=0._dp
1671 : : end if
1672 : :
1673 : : ! accumulate kinetic energy 2nd derivative with wavefunction components
1674 : 60846408 : term=d2fkin*dkpg21*dkpg22 + dfkin*d2kpg2
1675 [ - + - - ]: 60846408 : if(istwfk==2 .and. ig/=1)term=2.0_dp*term
1676 [ - + ]: 60846408 : if(istwfk>2)term=2.0_dp*term
1677 : 61338480 : d2kinacc=d2kinacc + term*(cwavef(re,ig+igs)**2 + cwavef(im,ig+igs)**2)
1678 : :
1679 : : end do !ig
1680 : : end do !ispinor
1681 : :
1682 : 632664 : ekinout(istr1+6*(istr2-1))=d2kinacc/effmass_free
1683 : :
1684 : : end do !istr1
1685 : : end do !istr2
1686 : :
1687 : 23432 : end subroutine d2kindstr2
1688 : : !!***
1689 : :
1690 : : end subroutine dfpt_eltfrkin
1691 : : !!***
1692 : :
1693 : : !!****f* ABINIT/dfpt_eltfrhar
1694 : : !! NAME
1695 : : !! dfpt_eltfrhar
1696 : : !!
1697 : : !! FUNCTION
1698 : : !! Compute the frozen-wavefunction hartree enegy contribution to the elastic tensor
1699 : : !!
1700 : : !! INPUTS
1701 : : !! rprimd(3,3)=dimensional primitive translation vectors (bohr)
1702 : : !! gsqcut =Fourier cutoff on G^2 for "large sphere" of radius double
1703 : : !! that of the basis sphere--appropriate for charge density rho(G),
1704 : : !! Hartree potential, and pseudopotentials
1705 : : !! mpi_enreg=information about MPI parallelization
1706 : : !! nfft =(effective) number of FFT grid points (for this processor)
1707 : : !! ngfft(18)=contain all needed information about 3D FFT,
1708 : : !! see ~abinit/doc/variables/vargs.htm#ngfft
1709 : : !! rhog(2,nfft)=total electron density in G space
1710 : : !!
1711 : : !! OUTPUT
1712 : : !! eltfrhar(6,6)=non-symmetrized kinetic energy contribution to the
1713 : : !! elastic tensor
1714 : : !! NOTES
1715 : : !! *based largely on hartre.f
1716 : : !!
1717 : : !! SOURCE
1718 : :
1719 : 37 : subroutine dfpt_eltfrhar(eltfrhar,rprimd,gsqcut,mpi_enreg,nfft,ngfft,rhog)
1720 : :
1721 : : !Arguments ------------------------------------
1722 : : !scalars
1723 : : integer,intent(in) :: nfft
1724 : : real(dp),intent(in) :: gsqcut
1725 : : type(MPI_type),intent(in) :: mpi_enreg
1726 : : !arrays
1727 : : integer,intent(in) :: ngfft(18)
1728 : : real(dp),intent(in) :: rhog(2,nfft),rprimd(3,3)
1729 : : real(dp),intent(out) :: eltfrhar(6,6)
1730 : :
1731 : : !Local variables-------------------------------
1732 : : !scalars
1733 : : integer,parameter :: im=2,re=1
1734 : : integer :: i1,i2,i23,i3,id2,id3,ierr,ig,ig2,ig3,ii,ii1,ing,istr1,istr2,jj
1735 : : integer :: ka,kb,kd,kg,me_fft,n1,n2,n3,nproc_fft
1736 : : real(dp),parameter :: tolfix=1.000000001_dp
1737 : : real(dp) :: cutoff,d2eacc,d2etot,d2gs,deacc01,deacc10,dgs01,dgs10,eacc,fact,gs
1738 : : real(dp) :: term,ucvol
1739 : : !arrays
1740 : : integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/)
1741 : : integer :: id(3)
1742 : 37 : integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
1743 : 37 : integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
1744 : : real(dp) :: d2gm(3,3),dgm01(3,3),dgm10(3,3),gmet(3,3),gprimd(3,3),gqr(3)
1745 : : real(dp) :: rmet(3,3),tsec(2)
1746 : 37 : real(dp),allocatable :: gq(:,:)
1747 : :
1748 : : ! *************************************************************************
1749 : :
1750 : : !Compute gmet, gprimd and ucvol from rprimd
1751 : 37 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
1752 : :
1753 : 37 : eltfrhar(:,:)=0.0_dp
1754 : :
1755 : 37 : n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
1756 : 37 : me_fft=ngfft(11)
1757 : 37 : nproc_fft=ngfft(10)
1758 : :
1759 : : !Get the distrib associated with this fft_grid
1760 : 37 : call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
1761 : :
1762 : : !Initialize a few quantities
1763 : 37 : fact=0.5_dp*ucvol/pi
1764 : 37 : cutoff=gsqcut*tolfix
1765 : :
1766 : : !In order to speed the routine, precompute the components of g+q
1767 : : !Also check if the booked space was large enough...
1768 [ - + + - : 111 : ABI_MALLOC(gq,(3,max(n1,n2,n3)))
+ - + - -
+ + - ]
1769 [ + + ]: 148 : do ii=1,3
1770 : 111 : id(ii)=ngfft(ii)/2+2
1771 [ + + ]: 2130 : do ing=1,ngfft(ii)
1772 : 1982 : ig=ing-(ing/id(ii))*ngfft(ii)-1
1773 : 2093 : gq(ii,ing)=ig
1774 : : end do
1775 : : end do
1776 : :
1777 : : !Loop over 2nd strain index
1778 [ + + ]: 259 : do istr2=1,6
1779 : : ! Loop over 1st strain index, upper triangle only
1780 [ + + ]: 1036 : do istr1=1,istr2
1781 : :
1782 : 777 : ka=idx(2*istr1-1);kb=idx(2*istr1);kg=idx(2*istr2-1);kd=idx(2*istr2)
1783 : :
1784 [ + + ]: 3108 : do ii = 1,3
1785 [ + + ]: 9324 : dgm01(:,ii)=-(gprimd(ka,:)*gprimd(kb,ii)+gprimd(kb,:)*gprimd(ka,ii))
1786 [ + + ]: 10101 : dgm10(:,ii)=-(gprimd(kg,:)*gprimd(kd,ii)+gprimd(kd,:)*gprimd(kg,ii))
1787 : : end do
1788 : :
1789 : 777 : d2gm(:,:)=0._dp
1790 [ + + ]: 3108 : do ii = 1,3
1791 [ + + ]: 2331 : if(ka==kg) d2gm(:,ii)=d2gm(:,ii)&
1792 [ + + ]: 4440 : & +gprimd(kb,:)*gprimd(kd,ii)+gprimd(kd,:)*gprimd(kb,ii)
1793 [ + + ]: 2331 : if(ka==kd) d2gm(:,ii)=d2gm(:,ii)&
1794 [ + + ]: 2664 : & +gprimd(kb,:)*gprimd(kg,ii)+gprimd(kg,:)*gprimd(kb,ii)
1795 [ + + ]: 2331 : if(kb==kg) d2gm(:,ii)=d2gm(:,ii)&
1796 [ + + ]: 3108 : & +gprimd(ka,:)*gprimd(kd,ii)+gprimd(kd,:)*gprimd(ka,ii)
1797 [ + + ]: 2331 : if(kb==kd) d2gm(:,ii)=d2gm(:,ii)&
1798 [ + + ]: 5217 : & +gprimd(ka,:)*gprimd(kg,ii)+gprimd(kg,:)*gprimd(ka,ii)
1799 : : end do
1800 [ + + + + ]: 10101 : d2gm(:,:)=0.5_dp*d2gm(:,:)
1801 : :
1802 : : ! initialize energy accumulator
1803 : 777 : eacc=0._dp
1804 : 777 : deacc01=0._dp
1805 : 777 : deacc10=0._dp
1806 : 777 : d2eacc=0._dp
1807 : :
1808 : 777 : id2=n2/2+2
1809 : 777 : id3=n3/2+2
1810 : : ! Triple loop on each dimension
1811 [ + + ]: 17115 : do i3=1,n3
1812 : 16338 : ig3=i3-(i3/id3)*n3-1
1813 : 16338 : gqr(3)=gq(3,i3)
1814 [ + + ]: 292299 : do i2=1,n2
1815 [ + - ]: 291522 : if (fftn2_distrib(i2)==me_fft) then
1816 : 275184 : gqr(2)=gq(2,i2)
1817 : 275184 : ig2=i2-(i2/id2)*n2-1
1818 : 275184 : i23=n1*(ffti2_local(i2)-1 +(n2/nproc_fft)*(i3-1))
1819 : : ! Do the test that eliminates the Gamma point outside
1820 : : ! of the inner loop
1821 : 275184 : ii1=1
1822 [ + + ]: 275184 : if(i23==0 .and. ig2==0 .and. ig3==0)then
1823 : 777 : ii1=2
1824 : : end if
1825 : :
1826 : : ! Final inner loop on the first dimension
1827 : : ! (note the lower limit)
1828 [ + + ]: 5003901 : do i1=ii1,n1
1829 : 4728717 : gqr(1)=gq(1,i1)
1830 : : gs=(gmet(1,1)*gqr(1)*gqr(1)+gmet(2,2)*gqr(2)*gqr(2)+&
1831 : : & gmet(3,3)*gqr(3)*gqr(3)+2._dp*&
1832 : : & (gmet(1,2)*gqr(1)*gqr(2) + gmet(1,3)*gqr(1)*gqr(3)+&
1833 : 4728717 : & gmet(2,3)*gqr(2)*gqr(3)) )
1834 : 4728717 : ii=i1+i23
1835 [ + + ]: 5003901 : if(gs<=cutoff)then
1836 : : dgs01=(dgm01(1,1)*gqr(1)*gqr(1)+dgm01(2,2)*gqr(2)*gqr(2)+&
1837 : : & dgm01(3,3)*gqr(3)*gqr(3)+2._dp*&
1838 : : & (dgm01(1,2)*gqr(1)*gqr(2) + dgm01(1,3)*gqr(1)*gqr(3)+&
1839 : 1552110 : & dgm01(2,3)*gqr(2)*gqr(3)) )
1840 : : dgs10=(dgm10(1,1)*gqr(1)*gqr(1)+dgm10(2,2)*gqr(2)*gqr(2)+&
1841 : : & dgm10(3,3)*gqr(3)*gqr(3)+2._dp*&
1842 : : & (dgm10(1,2)*gqr(1)*gqr(2) + dgm10(1,3)*gqr(1)*gqr(3)+&
1843 : 1552110 : & dgm10(2,3)*gqr(2)*gqr(3)) )
1844 : : d2gs =(d2gm(1,1)*gqr(1)*gqr(1)+d2gm(2,2)*gqr(2)*gqr(2)+&
1845 : : & d2gm(3,3)*gqr(3)*gqr(3)+2._dp*&
1846 : : & (d2gm(1,2)*gqr(1)*gqr(2) + d2gm(1,3)*gqr(1)*gqr(3)+&
1847 : 1552110 : & d2gm(2,3)*gqr(2)*gqr(3)) )
1848 : :
1849 : 1552110 : term=(rhog(re,ii)**2+rhog(im,ii)**2)/gs
1850 : 1552110 : eacc=eacc+term
1851 : 1552110 : deacc01=deacc01+dgs01*term/gs
1852 : 1552110 : deacc10=deacc10+dgs10*term/gs
1853 : 1552110 : d2eacc=d2eacc+(-d2gs+2._dp*dgs01*dgs10/gs)*term/gs
1854 : : end if
1855 : :
1856 : : ! End loop on i1
1857 : : end do
1858 : : end if
1859 : : ! End loop on i2
1860 : : end do
1861 : : ! End loop on i3
1862 : : end do
1863 : :
1864 : : ! Add contributions taking account diagonal strain terms (from ucvol
1865 : : ! derivatives)
1866 : 777 : d2etot=d2eacc
1867 [ + + ]: 777 : if(istr1<=3) d2etot=d2etot+deacc10
1868 [ + + ]: 777 : if(istr2<=3) d2etot=d2etot+deacc01
1869 [ + + ]: 777 : if(istr1<=3 .and. istr2<=3) d2etot=d2etot+eacc
1870 : :
1871 : 999 : eltfrhar(istr1,istr2)=fact*d2etot
1872 : :
1873 : : ! End loop on istr1
1874 : : end do
1875 : : ! End loop in istr2
1876 : : end do
1877 : :
1878 [ - + ]: 37 : ABI_FREE(gq)
1879 : :
1880 : : !Init mpi_comm
1881 : 37 : call timab(48,1,tsec)
1882 : 37 : call xmpi_sum(eltfrhar,mpi_enreg%comm_fft,ierr)
1883 : 37 : call timab(48,2,tsec)
1884 : :
1885 : : !Fill in lower triangle
1886 [ + + ]: 222 : do jj=2,6
1887 [ + + ]: 777 : do ii=1,jj-1
1888 : 740 : eltfrhar(jj,ii)=eltfrhar(ii,jj)
1889 : : end do
1890 : : end do
1891 : 37 : end subroutine dfpt_eltfrhar
1892 : : !!***
1893 : :
1894 : : !!****f* ABINIT/elt_ewald
1895 : : !!
1896 : : !! NAME
1897 : : !! elt_ewald
1898 : : !!
1899 : : !! FUNCTION
1900 : : !! Compute 2nd derivatives of Ewald energy wrt strain for frozen wavefunction
1901 : : !! contributions to elastic tensor
1902 : : !!
1903 : : !! INPUTS
1904 : : !! gmet(3,3)=metric tensor in reciprocal space (bohr^-2)
1905 : : !! gprimd(3,3)=dimensional primitive translations for reciprocal space (bohr^-1)
1906 : : !! mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
1907 : : !! comm_atom=--optional-- MPI communicator over atoms
1908 : : !! my_natom=number of atoms treated by current processor
1909 : : !! natom=number of atoms in unit cell
1910 : : !! ntypat=numbe of type of atoms
1911 : : !! rmet(3,3)=metric tensor in real space (bohr^2)
1912 : : !! rprimd(3,3)=dimensional primitive translation vectors (bohr)
1913 : : !! typat(natom)=integer label of each type of atom (1,2,...)
1914 : : !! ucvol=unit cell volume (bohr^3)
1915 : : !! xred(3,natom)=relative coords of atoms in unit cell (dimensionless)
1916 : : !! zion(ntypat)=charge on each type of atom (real number)
1917 : : !!
1918 : : !! OUTPUT
1919 : : !! elteew(6+3*natom,6)=2nd derivatives of Ewald energy wrt strain
1920 : : !!
1921 : : !! SOURCE
1922 : :
1923 : 74 : subroutine elt_ewald(elteew,gmet,gprimd,my_natom,natom,ntypat,rmet,rprimd,&
1924 : 37 : & typat,ucvol,xred,zion, &
1925 [ + - + + : 37 : & mpi_atmtab,comm_atom) ! optional arguments (parallelism)
- + ]
1926 : :
1927 : : !Arguments ------------------------------------
1928 : : !scalars
1929 : : integer,intent(in) :: my_natom,natom,ntypat
1930 : : real(dp),intent(in) :: ucvol
1931 : : !arrays
1932 : : integer,intent(in) :: typat(natom)
1933 : : integer,optional,intent(in) :: comm_atom
1934 : : integer,optional,target,intent(in) :: mpi_atmtab(:)
1935 : : real(dp),intent(in) :: gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3)
1936 : : real(dp),intent(in) :: xred(3,natom),zion(ntypat)
1937 : : real(dp),intent(out) :: elteew(6+3*natom,6)
1938 : :
1939 : : !Local variables-------------------------------
1940 : : !scalars
1941 : : integer :: ia,ia0,ib,ierr,ig1,ig2,ig3,ir1,ir2,ir3,is1,is2,jj,js,ka,kb,kd,kg,my_comm_atom,newg,newr,ng,nr
1942 : : logical :: my_atmtab_allocated,paral_atom
1943 : : real(dp) :: arg,ch,chsq,cos_term,d2derfc,d2gss,d2r,d2rs,dderfc,derfc_arg
1944 : : real(dp) :: dgss1,dgss2,direct,dr1,dr2,drs1,drs2,eew,eta,fac,fraca1,fraca2
1945 : : real(dp) :: fraca3,fracb1,fracb2,fracb3,gsq,gsum,r1,r2,r3,recip,reta
1946 : : real(dp) :: rmagn,rsq,sin_term,sumg,summi,summr,sumr,t1,term
1947 : : character(len=500) :: message
1948 : : !arrays
1949 : : integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/)
1950 : 37 : integer,pointer :: my_atmtab(:)
1951 : : real(dp) :: d2gm(3,3,6,6),d2ris(3),d2rm(3,3,6,6),dgm(3,3,6),dris(3),drm(3,3,6)
1952 : : real(dp) :: t2(3),ts2(3),tsec(2),tt(3)
1953 : 37 : real(dp),allocatable :: d2sumg(:,:),d2sumr(:,:),drhoisi(:,:),drhoisr(:,:)
1954 : 37 : real(dp),allocatable :: mpibuf(:)
1955 : :
1956 : : ! *************************************************************************
1957 : :
1958 : : !DEBUG
1959 : : !write(std_out,*)' elt_ewald : enter '
1960 : : !stop
1961 : : !ENDDEBUG
1962 : :
1963 : : !Compute 1st and 2nd derivatives of metric tensor wrt all strain components
1964 : : !and store for use in inner loop below.
1965 : :
1966 : : !Set up parallelism over atoms
1967 [ + + ]: 37 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
1968 [ + + ]: 37 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
1969 [ + - ]: 37 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
1970 : 37 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
1971 : :
1972 : : !Loop over 2nd strain index
1973 [ + + ]: 259 : do is2=1,6
1974 : 222 : kg=idx(2*is2-1);kd=idx(2*is2)
1975 [ + + ]: 888 : do jj = 1,3
1976 [ + + ]: 2664 : drm(:,jj,is2)=rprimd(kg,:)*rprimd(kd,jj)+rprimd(kd,:)*rprimd(kg,jj)
1977 [ + + ]: 2886 : dgm(:,jj,is2)=-(gprimd(kg,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(kg,jj))
1978 : : end do
1979 : :
1980 : : ! Loop over 1st strain index, upper triangle only
1981 [ + + ]: 1036 : do is1=1,is2
1982 : :
1983 : 777 : ka=idx(2*is1-1);kb=idx(2*is1)
1984 [ + + + + ]: 10101 : d2rm(:,:,is1,is2)=zero
1985 [ + + + + ]: 10101 : d2gm(:,:,is1,is2)=zero
1986 [ + + ]: 3330 : do jj = 1,3
1987 [ + + ]: 2331 : if(ka==kg) d2rm(:,jj,is1,is2)=d2rm(:,jj,is1,is2)&
1988 [ + + ]: 4440 : & +rprimd(kb,:)*rprimd(kd,jj)+rprimd(kd,:)*rprimd(kb,jj)
1989 [ + + ]: 2331 : if(ka==kd) d2rm(:,jj,is1,is2)=d2rm(:,jj,is1,is2)&
1990 [ + + ]: 2664 : & +rprimd(kb,:)*rprimd(kg,jj)+rprimd(kg,:)*rprimd(kb,jj)
1991 [ + + ]: 2331 : if(kb==kg) d2rm(:,jj,is1,is2)=d2rm(:,jj,is1,is2)&
1992 [ + + ]: 3108 : & +rprimd(ka,:)*rprimd(kd,jj)+rprimd(kd,:)*rprimd(ka,jj)
1993 [ + + ]: 2331 : if(kb==kd) d2rm(:,jj,is1,is2)=d2rm(:,jj,is1,is2)&
1994 [ + + ]: 4440 : & +rprimd(ka,:)*rprimd(kg,jj)+rprimd(kg,:)*rprimd(ka,jj)
1995 : :
1996 [ + + ]: 2331 : if(ka==kg) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
1997 [ + + ]: 4440 : & +gprimd(kb,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(kb,jj)
1998 [ + + ]: 2331 : if(ka==kd) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
1999 [ + + ]: 2664 : & +gprimd(kb,:)*gprimd(kg,jj)+gprimd(kg,:)*gprimd(kb,jj)
2000 [ + + ]: 2331 : if(kb==kg) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
2001 [ + + ]: 3108 : & +gprimd(ka,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(ka,jj)
2002 [ + + ]: 2331 : if(kb==kd) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
2003 [ + + ]: 5217 : & +gprimd(ka,:)*gprimd(kg,jj)+gprimd(kg,:)*gprimd(ka,jj)
2004 : : end do
2005 : : end do !is1
2006 : : end do !is2
2007 : :
2008 : : !Add up total charge and sum of $charge^2$ in cell
2009 : : chsq=zero
2010 : : ch=zero
2011 [ + + ]: 131 : do ia=1,natom
2012 : 94 : ch=ch+zion(typat(ia))
2013 : 131 : chsq=chsq+zion(typat(ia))**2
2014 : : end do
2015 : :
2016 : : !Compute eta, the Ewald summation convergence parameter,
2017 : : !for approximately optimized summations:
2018 : : direct=rmet(1,1)+rmet(1,2)+rmet(1,3)+rmet(2,1)+&
2019 : 37 : & rmet(2,2)+rmet(2,3)+rmet(3,1)+rmet(3,2)+rmet(3,3)
2020 : : recip=gmet(1,1)+gmet(1,2)+gmet(1,3)+gmet(2,1)+&
2021 : 37 : & gmet(2,2)+gmet(2,3)+gmet(3,1)+gmet(3,2)+gmet(3,3)
2022 : : !Here, a bias is introduced, because G-space summation scales
2023 : : !better than r space summation ! Note : debugging is the most
2024 : : !easier at fixed eta.
2025 : : ! eta=pi*200._dp/33.0_dp*sqrt(1.69_dp*recip/direct)
2026 : 37 : eta=1.0_dp
2027 : :
2028 : : !Conduct reciprocal space summations
2029 : 37 : fac=pi**2/eta ; gsum=zero
2030 [ + - - + : 185 : ABI_MALLOC(d2sumg,(6+3*natom,6))
+ - + - +
- + - - +
- + ]
2031 [ - + + - : 111 : ABI_MALLOC(drhoisr,(3,natom))
+ - + - -
+ - + ]
2032 [ - + + - : 74 : ABI_MALLOC(drhoisi,(3,natom))
- + - + ]
2033 [ + + + + ]: 3283 : d2sumg(:,:)=zero
2034 : :
2035 : : !Sum over G space, done shell after shell until all
2036 : : !contributions are too small.
2037 : 37 : ng=0
2038 : : do
2039 : 930 : ng=ng+1
2040 : 930 : newg=0
2041 : :
2042 [ + + ]: 28672 : do ig3=-ng,ng
2043 [ + + ]: 1174074 : do ig2=-ng,ng
2044 [ + + ]: 57209902 : do ig1=-ng,ng
2045 : :
2046 : : ! Exclude shells previously summed over
2047 : : if(abs(ig1)==ng .or. abs(ig2)==ng .or. abs(ig3)==ng&
2048 [ + + + + : 57182160 : & .or. ng==1 ) then
+ + + + ]
2049 : :
2050 : : ! gsq is G dot G = |G|^2
2051 : : gsq=gmet(1,1)*dble(ig1*ig1)+gmet(2,2)*dble(ig2*ig2)+&
2052 : : & gmet(3,3)*dble(ig3*ig3)+2._dp*(gmet(2,1)*dble(ig1*ig2)+&
2053 : 6546985 : & gmet(3,1)*dble(ig1*ig3)+gmet(3,2)*dble(ig3*ig2))
2054 : :
2055 : : ! Skip g=0:
2056 [ + + ]: 6546985 : if (gsq>1.0d-20) then
2057 : 6546948 : arg=fac*gsq
2058 : :
2059 : : ! Larger arg gives 0 contribution because of exp(-arg)
2060 [ + + ]: 6546948 : if (arg <= 80._dp) then
2061 : : ! When any term contributes then include next shell
2062 : 1275492 : newg=1
2063 : 1275492 : term=exp(-arg)/gsq
2064 : 1275492 : summr = zero
2065 : 1275492 : summi = zero
2066 : : ! Note that if reduced atomic coordinates xred drift outside
2067 : : ! of unit cell (outside [0,1)) it is irrelevant in the following
2068 : : ! term, which only computes a phase.
2069 [ + + ]: 5371384 : do ia=1,natom
2070 : 4095892 : arg=two_pi*(ig1*xred(1,ia)+ig2*xred(2,ia)+ig3*xred(3,ia))
2071 : : ! Sum real and imaginary parts (avoid complex variables)
2072 : 4095892 : cos_term=cos(arg)
2073 : 4095892 : sin_term=sin(arg)
2074 : 4095892 : summr=summr+zion(typat(ia))*cos_term
2075 : 4095892 : summi=summi+zion(typat(ia))*sin_term
2076 : 4095892 : drhoisr(1,ia)=-two_pi*zion(typat(ia))*sin_term*dble(ig1)
2077 : 4095892 : drhoisi(1,ia)= two_pi*zion(typat(ia))*cos_term*dble(ig1)
2078 : 4095892 : drhoisr(2,ia)=-two_pi*zion(typat(ia))*sin_term*dble(ig2)
2079 : 4095892 : drhoisi(2,ia)= two_pi*zion(typat(ia))*cos_term*dble(ig2)
2080 : 4095892 : drhoisr(3,ia)=-two_pi*zion(typat(ia))*sin_term*dble(ig3)
2081 : 5371384 : drhoisi(3,ia)= two_pi*zion(typat(ia))*cos_term*dble(ig3)
2082 : : end do
2083 : :
2084 : : ! The following two checks avoid an annoying
2085 : : ! underflow error message
2086 [ + + ]: 1275492 : if (abs(summr)<1.d-16) summr=zero
2087 [ + + ]: 1275492 : if (abs(summi)<1.d-16) summi=zero
2088 : :
2089 : : ! The product of term and summr**2 or summi**2 below
2090 : : ! can underflow if not for checks above
2091 : 1275492 : t1=term*(summr*summr+summi*summi)
2092 : 1275492 : gsum=gsum+t1
2093 : : ! Loop over 2nd strain index
2094 [ + + ]: 8928444 : do is2=1,6
2095 : : dgss2=dgm(1,1,is2)*dble(ig1*ig1)+dgm(2,2,is2)*dble(ig2*ig2)+&
2096 : : & dgm(3,3,is2)*dble(ig3*ig3)+2._dp*(dgm(2,1,is2)*dble(ig1*ig2)+&
2097 : 7652952 : & dgm(3,1,is2)*dble(ig1*ig3)+dgm(3,2,is2)*dble(ig3*ig2))
2098 : : ! Loop over 1st strain index, upper triangle only
2099 [ + + ]: 34438284 : do is1=1,is2
2100 : : dgss1=dgm(1,1,is1)*dble(ig1*ig1)+dgm(2,2,is1)*dble(ig2*ig2)+&
2101 : : & dgm(3,3,is1)*dble(ig3*ig3)+2._dp*(dgm(2,1,is1)*dble(ig1*ig2)+&
2102 : 26785332 : & dgm(3,1,is1)*dble(ig1*ig3)+dgm(3,2,is1)*dble(ig3*ig2))
2103 : :
2104 : : d2gss=d2gm(1,1,is1,is2)*dble(ig1*ig1)+&
2105 : : & d2gm(2,2,is1,is2)*dble(ig2*ig2)+&
2106 : : & d2gm(3,3,is1,is2)*dble(ig3*ig3)+2._dp*&
2107 : : & (d2gm(2,1,is1,is2)*dble(ig1*ig2)+&
2108 : : & d2gm(3,1,is1,is2)*dble(ig1*ig3)+&
2109 : 26785332 : & d2gm(3,2,is1,is2)*dble(ig3*ig2))
2110 : :
2111 : : d2sumg(is1,is2)=d2sumg(is1,is2)+&
2112 : : & t1*((fac**2 + 2.0_dp*fac/gsq + 2.0_dp/(gsq**2))*dgss1*dgss2 -&
2113 : 26785332 : & 0.5_dp*(fac + 1.0_dp/gsq)*d2gss)
2114 [ + + ]: 26785332 : if(is1<=3) d2sumg(is1,is2)=d2sumg(is1,is2)+&
2115 : 19132380 : & t1*(fac + 1.0_dp/gsq)*dgss2
2116 [ + + ]: 26785332 : if(is2<=3) d2sumg(is1,is2)=d2sumg(is1,is2)+&
2117 : 7652952 : & t1*(fac + 1.0_dp/gsq)*dgss1
2118 [ + + ]: 34438284 : if(is1<=3 .and. is2<=3) d2sumg(is1,is2)=d2sumg(is1,is2)+t1
2119 : :
2120 : : end do !is1
2121 : :
2122 : : ! Internal strain contributions
2123 [ + + ]: 33503796 : do ia=1,natom
2124 : 24575352 : js=7+3*(ia-1)
2125 [ + + ]: 98301408 : t2(:)=2.0_dp*term*(summr*drhoisr(:,ia)+summi*drhoisi(:,ia))
2126 : : d2sumg(js:js+2,is2)=d2sumg(js:js+2,is2)-&
2127 [ + + ]: 98301408 : & (fac + 1.0_dp/gsq)*dgss2*t2(:)
2128 [ + + + + ]: 69091332 : if(is2<=3) d2sumg(js:js+2,is2)=d2sumg(js:js+2,is2)-t2(:)
2129 : : end do
2130 : : end do !is2
2131 : :
2132 : : end if ! End condition of not larger than 80.0
2133 : : end if ! End skip g=0
2134 : : end if ! End triple loop over G s and associated new shell condition
2135 : : end do
2136 : : end do
2137 : : end do
2138 : :
2139 : : ! Check if new shell must be calculated
2140 [ + + ]: 930 : if (newg==0) exit
2141 : : end do ! End the loop on ng (new shells). Note that there is one exit from this loop.
2142 : :
2143 : 37 : sumg=gsum/(two_pi*ucvol)
2144 [ + + + + ]: 3283 : d2sumg(:,:)=d2sumg(:,:)/(two_pi*ucvol)
2145 : :
2146 [ - + ]: 37 : ABI_FREE(drhoisr)
2147 [ - + ]: 37 : ABI_FREE(drhoisi)
2148 : : !Stress tensor is now computed elsewhere (ewald2) hence do not need
2149 : : !length scale gradients (used to compute them here).
2150 : :
2151 : : !Conduct real space summations
2152 : 37 : reta=sqrt(eta)
2153 : 37 : fac=2._dp*sqrt(eta/pi)
2154 [ - + + - : 74 : ABI_MALLOC(d2sumr,(6+3*natom,6))
- + - + ]
2155 [ + + + + ]: 3283 : sumr=zero;d2sumr(:,:)=zero
2156 : :
2157 : : !In the following a summation is being conducted over all
2158 : : !unit cells (ir1, ir2, ir3) so it is appropriate to map all
2159 : : !reduced coordinates xred back into [0,1).
2160 : : !
2161 : : !Loop on shells in r-space as was done in g-space
2162 : 37 : nr=0
2163 : : do
2164 : 79 : nr=nr+1
2165 : 79 : newr=0
2166 : : !
2167 [ + + ]: 410 : do ir3=-nr,nr
2168 [ + + ]: 1913 : do ir2=-nr,nr
2169 [ + + ]: 9173 : do ir1=-nr,nr
2170 [ + + + + : 8842 : if( abs(ir3)==nr .or. abs(ir2)==nr .or. abs(ir1)==nr .or. nr==1 )then
+ + + + ]
2171 : :
2172 [ + + ]: 16773 : do ia0=1,my_natom
2173 [ + + ]: 11058 : ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
2174 : 11058 : js=7+3*(ia-1)
2175 : : ! Map reduced coordinate xred(mu,ia) into [0,1)
2176 : 11058 : fraca1=xred(1,ia)-aint(xred(1,ia))+0.5_dp-sign(0.5_dp,xred(1,ia))
2177 : 11058 : fraca2=xred(2,ia)-aint(xred(2,ia))+0.5_dp-sign(0.5_dp,xred(2,ia))
2178 : 11058 : fraca3=xred(3,ia)-aint(xred(3,ia))+0.5_dp-sign(0.5_dp,xred(3,ia))
2179 [ + + ]: 48641 : do ib=1,natom
2180 : 30244 : fracb1=xred(1,ib)-aint(xred(1,ib))+0.5_dp-sign(0.5_dp,xred(1,ib))
2181 : 30244 : fracb2=xred(2,ib)-aint(xred(2,ib))+0.5_dp-sign(0.5_dp,xred(2,ib))
2182 : 30244 : fracb3=xred(3,ib)-aint(xred(3,ib))+0.5_dp-sign(0.5_dp,xred(3,ib))
2183 : 30244 : r1=dble(ir1)+fracb1-fraca1
2184 : 30244 : r2=dble(ir2)+fracb2-fraca2
2185 : 30244 : r3=dble(ir3)+fracb3-fraca3
2186 : : rsq=rmet(1,1)*r1*r1+rmet(2,2)*r2*r2+rmet(3,3)*r3*r3+&
2187 : 30244 : & 2.0_dp*(rmet(2,1)*r2*r1+rmet(3,2)*r3*r2+rmet(3,1)*r1*r3)
2188 : :
2189 : : ! Avoid zero denominators in 'term':
2190 [ + + ]: 41302 : if (rsq>=1.0d-24) then
2191 : :
2192 : : ! Note: erfc(8) is about 1.1e-29,
2193 : : ! so do not bother with larger arg.
2194 : : ! Also: exp(-64) is about 1.6e-28,
2195 : : ! so do not bother with larger arg**2 in exp.
2196 : 30166 : term=zero
2197 [ + + ]: 30166 : if (eta*rsq<64.0_dp) then
2198 : 1380 : newr=1
2199 : 1380 : rmagn=sqrt(rsq)
2200 : 1380 : arg=reta*rmagn
2201 : : ! derfc computes the complementary error function
2202 : : ! dderfc is the derivative of the complementary error function
2203 : : ! d2derfc is the 2nd derivative of the complementary error function
2204 : 1380 : dderfc=-fac*exp(-eta*rsq)
2205 : 1380 : d2derfc=-2._dp*eta*rmagn*dderfc
2206 : 1380 : derfc_arg = abi_derfc(arg)
2207 : 1380 : term=derfc_arg/rmagn
2208 : 1380 : sumr=sumr+zion(typat(ia))*zion(typat(ib))*term
2209 [ + + ]: 5520 : tt(:)=rmet(:,1)*r1+rmet(:,2)*r2+rmet(:,3)*r3
2210 [ + + ]: 5520 : dris(:)=tt(:)/rmagn
2211 : : ! Loop over 2nd strain index
2212 [ + + ]: 9660 : do is2=1,6
2213 : : drs2=drm(1,1,is2)*r1*r1+drm(2,2,is2)*r2*r2+&
2214 : : & drm(3,3,is2)*r3*r3+&
2215 : : & 2.0_dp*(drm(2,1,is2)*r2*r1+drm(3,2,is2)*r3*r2+&
2216 : 8280 : & drm(3,1,is2)*r1*r3)
2217 : 8280 : dr2=0.5_dp*drs2/rmagn
2218 : : ! Loop over 1st strain index, upper triangle only
2219 [ + + ]: 37260 : do is1=1,is2
2220 : : drs1=drm(1,1,is1)*r1*r1+drm(2,2,is1)*r2*r2+&
2221 : : & drm(3,3,is1)*r3*r3+&
2222 : : & 2.0_dp*(drm(2,1,is1)*r2*r1+drm(3,2,is1)*r3*r2+&
2223 : 28980 : & drm(3,1,is1)*r1*r3)
2224 : 28980 : dr1=0.5_dp*drs1/rmagn
2225 : : d2rs=d2rm(1,1,is1,is2)*r1*r1+d2rm(2,2,is1,is2)*r2*r2+&
2226 : : & d2rm(3,3,is1,is2)*r3*r3+&
2227 : : & 2.0_dp*(d2rm(2,1,is1,is2)*r2*r1+d2rm(3,2,is1,is2)*r3*r2+&
2228 : 28980 : & d2rm(3,1,is1,is2)*r1*r3)
2229 : 28980 : d2r=(0.25_dp*d2rs-dr1*dr2)/rmagn
2230 : : d2sumr(is1,is2)=d2sumr(is1,is2)+&
2231 : : & zion(typat(ia))*zion(typat(ib))*&
2232 : : & ((d2derfc-2.0_dp*dderfc/rmagn+2.0_dp*derfc_arg/rsq)*dr1*dr2+&
2233 : 37260 : & (dderfc-derfc_arg/rmagn)*d2r)/rmagn
2234 : : end do !is1
2235 : : ! Internal strain contribution
2236 [ + + ]: 33120 : ts2(:)=drm(:,1,is2)*r1+drm(:,2,is2)*r2+drm(:,3,is2)*r3
2237 [ + + ]: 33120 : d2ris(:)=ts2(:)/rmagn-0.5_dp*drs2*tt(:)/(rsq*rmagn)
2238 : :
2239 : : d2sumr(js:js+2,is2)=d2sumr(js:js+2,is2)-&
2240 : : & 2.0_dp*zion(typat(ia))*zion(typat(ib))*&
2241 : : & ((d2derfc-2.0_dp*dderfc/rmagn+2.0_dp*derfc_arg/rsq)*dr2*dris(:)+&
2242 [ + + ]: 34500 : & (dderfc-derfc_arg/rmagn)*d2ris(:))/rmagn
2243 : : end do !is2
2244 : : end if
2245 : : end if ! End avoid zero denominators in'term'
2246 : :
2247 : : end do ! end loop over ib:
2248 : : end do ! end loop over ia:
2249 : : end if ! end triple loop over real space points and associated condition of new shell
2250 : : end do
2251 : : end do
2252 : : end do
2253 : :
2254 : : ! Check if new shell must be calculated
2255 [ + + ]: 79 : if(newr==0) exit
2256 : : end do ! End loop on nr (new shells). Note that there is an exit within the loop
2257 : :
2258 : : !In case of parallelism over atoms: communicate
2259 [ + + ]: 37 : if (paral_atom) then
2260 : 6 : call timab(48,1,tsec)
2261 [ + - - + : 18 : ABI_MALLOC(mpibuf,((6+3*natom)*6+1))
+ - + - +
- - + -
+ ]
2262 [ + + ]: 12 : mpibuf(1:(6+3*natom)*6)=reshape(d2sumr(:,:),shape=(/((6+3*natom)*6)/))
2263 : 6 : mpibuf((6+3*natom)*6+1)=sumr
2264 : 6 : call xmpi_sum(mpibuf,my_comm_atom,ierr)
2265 : 6 : sumr=mpibuf((6+3*natom)*6+1)
2266 [ + + ]: 18 : d2sumr(:,:)=reshape(mpibuf(1:(6+3*natom)*6),shape=(/(6+3*natom),6/))
2267 : 6 : ABI_FREE(mpibuf)
2268 : 12 : call timab(48,2,tsec)
2269 : : end if
2270 : :
2271 : 37 : sumr=0.5_dp*sumr
2272 [ + + + + ]: 3283 : d2sumr(:,:)=0.5_dp*d2sumr(:,:)
2273 : 37 : fac=pi*ch**2/(2.0_dp*eta*ucvol)
2274 : :
2275 : : !Finally assemble Ewald energy, eew
2276 : 37 : eew=sumg+sumr-chsq*reta/sqrt(pi)-fac
2277 : :
2278 [ + + + + ]: 3283 : elteew(:,:)=d2sumg(:,:)+d2sumr(:,:)
2279 : :
2280 : : !Additional term for all strains diagonal (from "fac" term in eew)
2281 [ + + + + ]: 481 : elteew(1:3,1:3)=elteew(1:3,1:3)-fac
2282 : :
2283 : : !Fill in lower triangle
2284 [ + + ]: 222 : do is2=2,6
2285 [ + + ]: 777 : do is1=1,is2-1
2286 : 740 : elteew(is2,is1)=elteew(is1,is2)
2287 : : end do
2288 : : end do
2289 : :
2290 [ - + ]: 37 : ABI_FREE(d2sumg)
2291 [ - + ]: 37 : ABI_FREE(d2sumr)
2292 : :
2293 : : !Output the final values of ng and nr
2294 : 37 : write(message, '(a,i4,a,i4)' )' elt_ewald : nr and ng are ',nr,' and ',ng
2295 : 37 : call wrtout(std_out,message,'COLL')
2296 : :
2297 : : !Destroy atom table used for parallelism
2298 : 37 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
2299 : :
2300 [ + - - + ]: 74 : end subroutine elt_ewald
2301 : : !!***
2302 : :
2303 : : !!****f* ABINIT/dfpt_ewald
2304 : : !!
2305 : : !! NAME
2306 : : !! dfpt_ewald
2307 : : !!
2308 : : !! FUNCTION
2309 : : !! Compute ewald contribution to the dynamical matrix, at a given q wavevector.
2310 : : !! Note: the q=0 part should be subtracted, by another call to
2311 : : !! the present routine, with q=0. The present routine correspond
2312 : : !! to the quantity C_bar defined in Eq.(24) or (27) in Phys. Rev. B 55, 10355 (1997) [[cite:Gonze1997a]].
2313 : : !! The two calls correspond to Eq.(23) of the same paper.
2314 : : !! If q=0 is asked, sumg0 should be put to 0. Otherwise, it should be put to 1.
2315 : : !!
2316 : : !! INPUTS
2317 : : !! gmet(3,3)=metric tensor in reciprocal space (length units **-2)
2318 : : !! mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
2319 : : !! comm_atom=--optional-- MPI communicator over atoms
2320 : : !! icutcoul= type of Coulomb cutoff to apply
2321 : : !! my_natom=number of atoms treated by current processor
2322 : : !! natom=number of atoms in unit cell
2323 : : !! qphon(3)=phonon wavevector (same system of coordinates as the
2324 : : !! reciprocal lattice vectors)
2325 : : !! rmet(3,3)=metric tensor in real space (length units squared)
2326 : : !! sumg0: if=1, the sum in reciprocal space must include g=0,
2327 : : !! if=0, this contribution must be skipped (q=0 singularity)
2328 : : !! typat(natom)=integer label of each type of atom (1,2,...)
2329 : : !! ucvol=unit cell volume in (whatever length scale units)**3
2330 : : !! vcutgeo(3)= array to describe the geometry of the Coulomb cutoff
2331 : : !! xred(3,natom)=relative coords of atoms in unit cell (dimensionless)
2332 : : !! zion(ntypat)=charge on each type of atom (real number)
2333 : : !!
2334 : : !! OUTPUT
2335 : : !! dyew(2,3,natom,3,natom)= Ewald part of the dynamical matrix,
2336 : : !! second energy derivative wrt xred(3,natom), Hartrees.
2337 : : !!
2338 : : !! SOURCE
2339 : :
2340 : 1800 : subroutine dfpt_ewald(dyew,gmet,gsqcut,icutcoul,my_natom,natom,ngfft,nkpt,qphon,rcut, &
2341 : 900 : & rmet,rprimd,sumg0,typat,ucvol,vcutgeo,xred,zion, &
2342 [ + - + + : 900 : & mpi_atmtab,comm_atom ) ! optional arguments (parallelism))
- + ]
2343 : :
2344 : : !Arguments -------------------------------
2345 : : !scalars
2346 : : integer,intent(in) :: icutcoul,my_natom,natom,nkpt,sumg0
2347 : : real(dp),intent(in) :: gsqcut,rcut,ucvol
2348 : : !arrays
2349 : : integer,intent(in) :: ngfft(18),typat(natom)
2350 : : integer,optional,intent(in) :: comm_atom
2351 : : integer,optional,target,intent(in) :: mpi_atmtab(:)
2352 : : real(dp),intent(in) :: gmet(3,3),qphon(3),rmet(3,3),rprimd(3,3),vcutgeo(3),xred(3,natom),zion(*)
2353 : : real(dp),intent(out) :: dyew(2,3,natom,3,natom)
2354 : :
2355 : : !Local variables -------------------------
2356 : : !nr, ng affect convergence of sums (nr=3,ng=5 is not good enough):
2357 : : !scalars
2358 : : integer,parameter :: im=2,ng=10,nr=6,re=1
2359 : : integer :: ia,ia0,ib,ierr,ig,ig1,ig2,ig23,ig3,ii,ing
2360 : : integer :: ir1,ir2,ir3,mu,my_comm_atom,nh,nu
2361 : : logical :: my_atmtab_allocated,paral_atom,computeit
2362 : : real(dp) :: arg,arga,argb,c1i,c1r,da1,da2,da3,derfc_arg
2363 : : real(dp) :: direct,dot1,dot2,dot3,dotr1,dotr2,dotr3
2364 : : real(dp) :: eta,fac,gdot12,gdot13,gdot23,gsq,gsum,norm1
2365 : : real(dp) :: r1,r2,r3,rdot12,rdot13,rdot23,recip,reta
2366 : : real(dp) :: reta3m,rmagn,rsq,term,term1,term2
2367 : : real(dp) :: term3,facg0
2368 : : character(len=500) :: message
2369 : : !arrays
2370 : : integer :: id(3)
2371 : 900 : integer, allocatable :: inv_ig(:,:)
2372 : : real(dp) :: tsec(2)
2373 : 900 : integer,pointer :: my_atmtab(:)
2374 : : real(dp) :: gpq(3),rq(3)
2375 : 900 : real(dp),allocatable :: gcutoff(:)
2376 : :
2377 : : ! *************************************************************************
2378 : :
2379 : : !Set up parallelism over atoms
2380 [ + + ]: 900 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
2381 [ + + ]: 900 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
2382 [ + - ]: 900 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
2383 : 900 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
2384 : :
2385 : : !Compute eta for approximately optimized summations:
2386 : : direct=rmet(1,1)+rmet(1,2)+rmet(1,3)+rmet(2,1)+&
2387 : 900 : & rmet(2,2)+rmet(2,3)+rmet(3,1)+rmet(3,2)+rmet(3,3)
2388 : : recip=gmet(1,1)+gmet(1,2)+gmet(1,3)+gmet(2,1)+&
2389 : 900 : & gmet(2,2)+gmet(2,3)+gmet(3,1)+gmet(3,2)+gmet(3,3)
2390 : 900 : eta=pi*(dble(ng)/dble(nr))*sqrt(1.69_dp*recip/direct)
2391 : :
2392 : : !Test Ewald s summation
2393 : : !eta=1.2_dp*eta
2394 : :
2395 : : !Initialize Gcut-off array from m_gtermcutoff
2396 : : call termcutoff(gcutoff,gsqcut,icutcoul,ngfft,nkpt,rcut,rprimd,vcutgeo,&
2397 : 900 : & qpt=qphon)
2398 : :
2399 : : !Need a way to reverse mapping the indexes inside termcutoff and Ewald rotines
2400 : 900 : nh=MAX(ngfft(1),ngfft(2),ngfft(3))/2
2401 [ - + + - : 2700 : ABI_MALLOC(inv_ig,(3,-nh:nh))
+ - + - -
+ - + ]
2402 [ + + ]: 3600 : do ii=1,3
2403 : 2700 : id(ii)=ngfft(ii)/2+2
2404 [ + + ]: 50628 : do ing=1,ngfft(ii)
2405 : 47028 : ig=ing-(ing/id(ii))*ngfft(ii)-1
2406 : : ! Create reverse mapping
2407 : 49728 : inv_ig(ii, ig) = ing
2408 : : end do
2409 : : end do
2410 : :
2411 : : !Sum terms over g space:
2412 : 900 : fac=pi**2/eta
2413 : 900 : gsum=zero
2414 : 900 : da1=zero
2415 : 900 : da2=zero
2416 : 900 : da3=zero
2417 [ + + + + : 180724 : dyew(:,:,:,:,:)=zero
+ + + + +
+ ]
2418 : : ii=0
2419 [ + + ]: 19800 : do ig3=-ng,ng
2420 [ + + ]: 416700 : do ig2=-ng,ng
2421 [ + + ]: 8750700 : do ig1=-ng,ng
2422 : 8334900 : ii=ii+1
2423 : 8334900 : gpq(1)=dble(ig1)+qphon(1)
2424 : 8334900 : gpq(2)=dble(ig2)+qphon(2)
2425 : 8334900 : gpq(3)=dble(ig3)+qphon(3)
2426 : 8334900 : gdot12=gmet(2,1)*gpq(1)*gpq(2)
2427 : 8334900 : gdot13=gmet(3,1)*gpq(1)*gpq(3)
2428 : 8334900 : gdot23=gmet(3,2)*gpq(2)*gpq(3)
2429 : 8334900 : dot1=gmet(1,1)*gpq(1)**2+gdot12+gdot13
2430 : 8334900 : dot2=gmet(2,2)*gpq(2)**2+gdot12+gdot23
2431 : 8334900 : dot3=gmet(3,3)*gpq(3)**2+gdot13+gdot23
2432 : 8334900 : gsq=dot1+dot2+dot3
2433 : 8334900 : facg0= zero
2434 : 8334900 : computeit= .true.
2435 : : ! Skip q=0:
2436 [ + + ]: 8334900 : if (gsq<1.0d-20) then
2437 : 690 : computeit= .false.
2438 [ - + ]: 690 : if (sumg0==1) then
2439 : : write(message,'(5a)')&
2440 : 0 : & 'The phonon wavelength should not be zero : ',ch10,&
2441 : 0 : & 'there are non-analytical terms that the code cannot handle.',ch10,&
2442 : 0 : & 'Action : subtract this wavelength from the input.'
2443 : 0 : ABI_ERROR(message)
2444 : : end if
2445 : : else
2446 [ + - - - : 8334210 : if (icutcoul==55.and.(ig1==0 .and. ig2==0 .and. ig3==0).and. sumg0==0) then
- - ]
2447 : : facg0= one
2448 : : end if
2449 : : ! Endif g/=0 :
2450 : : end if
2451 : :
2452 : 396900 : if (computeit) then
2453 : 8334210 : arg=fac*gsq
2454 : : ! Larger arg gives 0 contribution:
2455 [ + + ]: 8334210 : if (arg <= 80._dp) then
2456 : :
2457 : : ! Apply cutoff
2458 : : if ((abs(ig1).lt.ngfft(1)/2).and.&
2459 [ + + + + : 962086 : & (abs(ig2).lt.ngfft(2)/2).and.&
+ + ]
2460 : : & (abs(ig3).lt.ngfft(3)/2)) then
2461 : : ! Use inv_ig to map back to ing indices
2462 : 702357 : ig23=ngfft(1)*(inv_ig(2,ig2)-1 + ngfft(2)*(inv_ig(3,ig3)-1))
2463 : 702357 : ii=inv_ig(1,ig1)+ig23
2464 : 702357 : term=exp(-arg-facg0)/gsq * gcutoff(ii)
2465 [ + - ]: 259729 : else if (icutcoul.ne.3) then
2466 : : term=zero
2467 : : else
2468 : 259729 : term=exp(-arg)/gsq
2469 : : end if
2470 : :
2471 [ + + ]: 3179405 : do ia0=1,my_natom
2472 [ + + ]: 2217319 : ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
2473 : 2217319 : arga=two_pi*(gpq(1)*xred(1,ia)+gpq(2)*xred(2,ia)+gpq(3)*xred(3,ia))
2474 [ + + ]: 9399142 : do ib=1,ia
2475 : 6219737 : argb=two_pi*(gpq(1)*xred(1,ib)+gpq(2)*xred(2,ib)+gpq(3)*xred(3,ib))
2476 : 6219737 : arg=arga-argb
2477 : 6219737 : c1r=cos(arg)*term
2478 : 6219737 : c1i=sin(arg)*term
2479 : :
2480 [ + + ]: 27096267 : do mu=1,3
2481 [ + + ]: 62197370 : do nu=1,mu
2482 : 37318422 : dyew(re,mu,ia,nu,ib)=dyew(re,mu,ia,nu,ib)+gpq(mu)*gpq(nu)*c1r
2483 : 55977633 : dyew(im,mu,ia,nu,ib)=dyew(im,mu,ia,nu,ib)+gpq(mu)*gpq(nu)*c1i
2484 : : end do
2485 : : end do
2486 : :
2487 : : end do
2488 : : end do
2489 : : end if
2490 : : ! Endif computeit:
2491 : : end if
2492 : : ! End triple loop over G s:
2493 : : end do
2494 : : end do
2495 : : end do
2496 : :
2497 [ - + ]: 900 : ABI_FREE(gcutoff)
2498 [ - + ]: 900 : ABI_FREE(inv_ig)
2499 : :
2500 : : !End G summation by accounting for some common factors.
2501 : : !(for the charges:see end of routine)
2502 : 900 : norm1=4.0_dp*pi/ucvol
2503 [ + + ]: 2654 : do ia0=1,my_natom
2504 [ + + ]: 1754 : ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
2505 [ + + ]: 6060 : do ib=1,ia
2506 [ + + ]: 15378 : do mu=1,3
2507 [ + + ]: 34060 : do nu=1,mu
2508 [ + + ]: 71526 : dyew(:,mu,ia,nu,ib)=dyew(:,mu,ia,nu,ib)*norm1
2509 : : end do
2510 : : end do
2511 : : end do
2512 : : end do
2513 : :
2514 : : !Do sums over real space:
2515 : 900 : reta=sqrt(eta)
2516 : 900 : reta3m=-eta*reta
2517 : 900 : fac=4._dp/3.0_dp/sqrt(pi)
2518 [ + + ]: 12600 : do ir3=-nr,nr
2519 [ + + ]: 164700 : do ir2=-nr,nr
2520 [ + + ]: 2141100 : do ir1=-nr,nr
2521 : 1977300 : arg=-two_pi*(qphon(1)*ir1+qphon(2)*ir2+qphon(3)*ir3)
2522 : 1977300 : c1r=cos(arg)*reta3m
2523 : 1977300 : c1i=sin(arg)*reta3m
2524 [ + + ]: 5982938 : do ia0=1,my_natom
2525 [ + + ]: 3853538 : ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
2526 [ + + ]: 13313820 : do ib=1,ia
2527 : 7482982 : r1=dble(ir1)+xred(1,ia)-xred(1,ib)
2528 : 7482982 : r2=dble(ir2)+xred(2,ia)-xred(2,ib)
2529 : 7482982 : r3=dble(ir3)+xred(3,ia)-xred(3,ib)
2530 : 7482982 : rdot12=rmet(2,1)*r1*r2
2531 : 7482982 : rdot13=rmet(3,1)*r1*r3
2532 : 7482982 : rdot23=rmet(3,2)*r2*r3
2533 : 7482982 : dotr1=rmet(1,1)*r1**2+rdot12+rdot13
2534 : 7482982 : dotr2=rmet(2,2)*r2**2+rdot12+rdot23
2535 : 7482982 : dotr3=rmet(3,3)*r3**2+rdot13+rdot23
2536 : 7482982 : rsq=dotr1+dotr2+dotr3
2537 : 7482982 : rmagn=sqrt(rsq)
2538 : : ! Avoid zero denominators in term :
2539 [ + + ]: 11336520 : if (rmagn>=1.0d-12) then
2540 : 7481228 : arg=reta*rmagn
2541 : 7481228 : term=zero
2542 [ + + ]: 7481228 : if (arg<8.0_dp) then
2543 : : ! Note: erfc(8) is about 1.1e-29,
2544 : : ! so don t bother with larger arg.
2545 : : ! Also: exp(-64) is about 1.6e-28,
2546 : : ! so don t bother with larger arg**2 in exp.
2547 : 883974 : derfc_arg = abi_derfc(arg)
2548 : 883974 : term=derfc_arg/arg**3
2549 : 883974 : term1=2.0_dp/sqrt(pi)*exp(-arg**2)/arg**2
2550 : 883974 : term2=-(term+term1)
2551 : 883974 : term3=(3*term+term1*(3.0_dp+2.0_dp*arg**2))/rsq
2552 : 883974 : rq(1)=rmet(1,1)*r1+rmet(1,2)*r2+rmet(1,3)*r3
2553 : 883974 : rq(2)=rmet(2,1)*r1+rmet(2,2)*r2+rmet(2,3)*r3
2554 : 883974 : rq(3)=rmet(3,1)*r1+rmet(3,2)*r2+rmet(3,3)*r3
2555 [ + + ]: 3535896 : do mu=1,3
2556 [ + + ]: 8839740 : do nu=1,mu
2557 : : dyew(re,mu,ia,nu,ib)=dyew(re,mu,ia,nu,ib)+&
2558 : 5303844 : & c1r*(rq(mu)*rq(nu)*term3+rmet(mu,nu)*term2)
2559 : : dyew(im,mu,ia,nu,ib)=dyew(im,mu,ia,nu,ib)+&
2560 : 7955766 : & c1i*(rq(mu)*rq(nu)*term3+rmet(mu,nu)*term2)
2561 : : end do
2562 : : end do
2563 : : end if
2564 : : else
2565 [ - + ]: 1754 : if (ia/=ib)then
2566 : : write(message,'(a,a,a,a,a,i5,a,i5,a)')&
2567 : 0 : & 'The distance between two atoms vanishes.',ch10,&
2568 : 0 : & 'This is not allowed.',ch10,&
2569 : 0 : & 'Action: check the input for the atoms number',ia,' and',ib,'.'
2570 : 0 : ABI_ERROR(message)
2571 : : else
2572 [ + + ]: 7016 : do mu=1,3
2573 [ + + ]: 17540 : do nu=1,mu
2574 : : dyew(re,mu,ia,nu,ib)=dyew(re,mu,ia,nu,ib)+&
2575 : 15786 : & fac*reta3m*rmet(mu,nu)
2576 : : end do
2577 : : end do
2578 : : end if
2579 : : end if
2580 : :
2581 : : end do ! End loop over ib:
2582 : : end do ! End loop over ia:
2583 : : end do ! End triple loop over real space points:
2584 : : end do
2585 : : end do
2586 : :
2587 : : !Take account of the charges
2588 : : !write(std_out,*)' '
2589 [ + + ]: 2654 : do ia0=1,my_natom
2590 [ + + ]: 1754 : ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
2591 [ + + ]: 6060 : do ib=1,ia
2592 [ + + ]: 15378 : do mu=1,3
2593 [ + + ]: 34060 : do nu=1,mu
2594 [ + + ]: 71526 : do ii=1,2
2595 : : ! write(std_out,*)dyew(ii,mu,ia,nu,ib)
2596 : : dyew(ii,mu,ia,nu,ib)=dyew(ii,mu,ia,nu,ib)*&
2597 : 61308 : & zion(typat(ia))*zion(typat(ib))
2598 : : end do
2599 : : end do
2600 : : end do
2601 : : end do
2602 : : end do
2603 : :
2604 : : !Symmetrize with respect to the directions
2605 [ + + ]: 2654 : do ia0=1,my_natom
2606 [ + + ]: 1754 : ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
2607 [ + + ]: 6060 : do ib=1,ia
2608 [ + + ]: 15378 : do mu=1,3
2609 [ + + ]: 34060 : do nu=1,mu
2610 : 20436 : dyew(re,nu,ia,mu,ib)=dyew(re,mu,ia,nu,ib)
2611 : 30654 : dyew(im,nu,ia,mu,ib)=dyew(im,mu,ia,nu,ib)
2612 : : end do
2613 : : end do
2614 : : end do
2615 : : end do
2616 : :
2617 : : !In case of parallelism over atoms: communicate
2618 [ + + ]: 900 : if (paral_atom) then
2619 : 44 : call timab(48,1,tsec)
2620 : 44 : call xmpi_sum(dyew,my_comm_atom,ierr)
2621 : 44 : call timab(48,2,tsec)
2622 : : end if
2623 : :
2624 : : !Fill the upper part of the matrix, with the hermitian conjugate
2625 [ + + ]: 2776 : do ia=1,natom
2626 [ + + ]: 6586 : do ib=1,ia
2627 [ + + ]: 17116 : do nu=1,3
2628 [ + + ]: 49530 : do mu=1,3
2629 : 34290 : dyew(re,mu,ib,nu,ia)=dyew(re,mu,ia,nu,ib)
2630 : 45720 : dyew(im,mu,ib,nu,ia)=-dyew(im,mu,ia,nu,ib)
2631 : : end do
2632 : : end do
2633 : : end do
2634 : : end do
2635 : :
2636 : : !Destroy atom table used for parallelism
2637 : 900 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
2638 : :
2639 [ + - ]: 900 : end subroutine dfpt_ewald
2640 : : !!***
2641 : :
2642 : : !!****f* ABINIT/dfpt_ewalddq
2643 : : !!
2644 : : !! NAME
2645 : : !! dfpt_ewalddq
2646 : : !!
2647 : : !! FUNCTION
2648 : : !! Compute the first q-gradient of Ewald contribution to the dynamical matrix, at a given q wavevector.
2649 : : !! If q=0 is asked, sumg0 should be put to 0. Otherwise, it should be put to 1.
2650 : : !!
2651 : : !! COPYRIGHT
2652 : : !! Copyright (C) 1998-2026 ABINIT group (MR, MS)
2653 : : !! This file is distributed under the terms of the
2654 : : !! GNU General Public License, see ~abinit/COPYING
2655 : : !! or http://www.gnu.org/copyleft/gpl.txt .
2656 : : !!
2657 : : !! INPUTS
2658 : : !! gmet(3,3)=metric tensor in reciprocal space (length units **-2)
2659 : : !! mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
2660 : : !! comm_atom=--optional-- MPI communicator over atoms
2661 : : !! my_natom=number of atoms treated by current processor
2662 : : !! natom=number of atoms in unit cell
2663 : : !! qphon(3)=phonon wavevector (same system of coordinates as the
2664 : : !! reciprocal lattice vectors)
2665 : : !! rmet(3,3)=metric tensor in real space (length units squared)
2666 : : !! sumg0: if=1, the sum in reciprocal space must include g=0,
2667 : : !! if=0, this contribution must be skipped (q=0 singularity)
2668 : : !! typat(natom)=integer label of each type of atom (1,2,...)
2669 : : !! ucvol=unit cell volume in (whatever length scale units)**3
2670 : : !! xred(3,natom)=relative coords of atoms in unit cell (dimensionless)
2671 : : !! zion(ntypat)=charge on each type of atom (real number)
2672 : : !!
2673 : : !! OUTPUT
2674 : : !! dyewdq(2,3,natom,3,natom,3)= First q-gradient of Ewald part of the dynamical matrix,
2675 : : !! second energy derivative wrt xred(3,natom), Hartrees.
2676 : : !!
2677 : : !! SOURCE
2678 : :
2679 : 6 : subroutine dfpt_ewalddq(dyewdq,gmet,my_natom,natom,qphon,rmet,sumg0,typat,ucvol,xred,zion, &
2680 [ - + - - : 3 : & mpi_atmtab,comm_atom ) ! optional arguments (parallelism))
- - ]
2681 : :
2682 : : !Arguments -------------------------------
2683 : : !scalars
2684 : : integer,intent(in) :: my_natom,natom,sumg0
2685 : : real(dp),intent(in) :: ucvol
2686 : : !arrays
2687 : : integer,intent(in) :: typat(natom)
2688 : : integer,optional,intent(in) :: comm_atom
2689 : : integer,optional,target,intent(in) :: mpi_atmtab(:)
2690 : : real(dp),intent(in) :: gmet(3,3),qphon(3),rmet(3,3),xred(3,natom),zion(*)
2691 : : real(dp),intent(out) :: dyewdq(2,3,natom,3,natom,3)
2692 : :
2693 : : !Local variables -------------------------
2694 : : !nr, ng affect convergence of sums (nr=3,ng=5 is not good enough):
2695 : : !scalars
2696 : : integer,parameter :: im=2,nng=10,nnr=6,re=1
2697 : : integer ::ia,ia0,ib,ierr,ig1,ig2,ig3,ii,iq,ir1,ir2,ir3,mu,my_comm_atom,newg,newr,ng,nr,nu
2698 : : logical :: my_atmtab_allocated,paral_atom
2699 : : real(dp) :: arg,arga,argb,c1i,c1r,delag,delbg,derfc_arg
2700 : : real(dp) :: direct,dot1,dot2,dot3,dotr1,dotr2,dotr3
2701 : : real(dp) :: eta,fac,fac2,gdot12,gdot13,gdot23,gsq,gpqdq,gterms,norm1
2702 : : real(dp) :: r1,r2,r3,rdot12,rdot13,rdot23,recip,reta
2703 : : real(dp) :: reta3m,rmagn,rsq,term,term1,term2,term3
2704 : : character(len=500) :: message
2705 : : !arrays
2706 : : real(dp) :: tsec(2)
2707 : 3 : integer,pointer :: my_atmtab(:)
2708 : : real(dp) :: dakk(3),gpq(3),rq(3)
2709 : :
2710 : : ! *************************************************************************
2711 : :
2712 : : !Set up parallelism over atoms
2713 [ + - ]: 3 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
2714 [ - + ]: 3 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
2715 [ + - ]: 3 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
2716 : 3 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
2717 : :
2718 : : !Compute eta for approximately optimized summations:
2719 : : direct=rmet(1,1)+rmet(1,2)+rmet(1,3)+rmet(2,1)+&
2720 : 3 : & rmet(2,2)+rmet(2,3)+rmet(3,1)+rmet(3,2)+rmet(3,3)
2721 : : recip=gmet(1,1)+gmet(1,2)+gmet(1,3)+gmet(2,1)+&
2722 : 3 : & gmet(2,2)+gmet(2,3)+gmet(3,1)+gmet(3,2)+gmet(3,3)
2723 : 3 : eta=pi*(dble(nng)/dble(nnr))*sqrt(1.69_dp*recip/direct)
2724 : :
2725 : : !Test Ewald s summation
2726 : : !eta=1.2_dp*eta
2727 : :
2728 : : !Sum over G space, done shell after shell until all
2729 : : !contributions are too small.
2730 : 3 : fac=pi**2.d0/eta
2731 : 3 : fac2=2.d0*fac
2732 [ + + + + : 1164 : dyewdq(:,:,:,:,:,:)=zero
+ + + + +
+ + + ]
2733 : : ng=0
2734 : : do
2735 : 18 : ng=ng+1
2736 : 18 : newg=0
2737 : :
2738 [ + + ]: 162 : do ig3=-ng,ng
2739 [ + + ]: 1524 : do ig2=-ng,ng
2740 [ + + ]: 15762 : do ig1=-ng,ng
2741 : :
2742 : : ! Exclude shells previously summed over
2743 : : if(abs(ig1)==ng .or. abs(ig2)==ng .or. abs(ig3)==ng&
2744 [ + + + + : 15618 : & .or. ng==1 ) then
+ + + + ]
2745 : :
2746 : 6591 : gpq(1)=dble(ig1)+qphon(1)
2747 : 6591 : gpq(2)=dble(ig2)+qphon(2)
2748 : 6591 : gpq(3)=dble(ig3)+qphon(3)
2749 : 6591 : gdot12=gmet(2,1)*gpq(1)*gpq(2)
2750 : 6591 : gdot13=gmet(3,1)*gpq(1)*gpq(3)
2751 : 6591 : gdot23=gmet(3,2)*gpq(2)*gpq(3)
2752 : 6591 : dot1=gmet(1,1)*gpq(1)**2+gdot12+gdot13
2753 : 6591 : dot2=gmet(2,2)*gpq(2)**2+gdot12+gdot23
2754 : 6591 : dot3=gmet(3,3)*gpq(3)**2+gdot13+gdot23
2755 : 6591 : gsq=dot1+dot2+dot3
2756 : : ! Skip q=0:
2757 [ + + ]: 6591 : if (gsq<1.0d-20) then
2758 [ - + ]: 3 : if (sumg0==1) then
2759 : : write(message,'(3a)')&
2760 : 0 : & 'The G=0 term has no contributions at first order in q: ',ch10,&
2761 : 0 : & 'Action : sumg0=0 '
2762 : 0 : ABI_ERROR(message)
2763 : : end if
2764 : : else
2765 : 6588 : arg=fac*gsq
2766 : : ! Larger arg gives 0 contribution:
2767 [ + + ]: 6588 : if (arg <= 80._dp) then
2768 : : ! When any term contributes then include next shell
2769 : 1230 : newg=1
2770 : 1230 : term=exp(-arg)/gsq
2771 [ + + ]: 3690 : do ia0=1,my_natom
2772 [ - + ]: 2460 : ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
2773 : 2460 : arga=two_pi*(gpq(1)*xred(1,ia)+gpq(2)*xred(2,ia)+gpq(3)*xred(3,ia))
2774 [ + + ]: 7380 : do ib=1,ia
2775 : 3690 : argb=two_pi*(gpq(1)*xred(1,ib)+gpq(2)*xred(2,ib)+gpq(3)*xred(3,ib))
2776 : 3690 : arg=arga-argb
2777 : 3690 : c1r=cos(arg)*term
2778 : 3690 : c1i=sin(arg)*term
2779 : :
2780 [ + + ]: 17220 : do iq=1,3
2781 : 11070 : gpqdq=gmet(iq,1)*gpq(1)+gmet(iq,2)*gpq(2)+gmet(iq,3)*gpq(3)
2782 [ + + ]: 47970 : do mu=1,3
2783 [ + + ]: 33210 : delag=zero; if(iq==mu) delag=one
2784 [ + + ]: 110700 : do nu=1,mu
2785 [ + + ]: 66420 : delbg=zero; if(iq==nu) delbg=one
2786 : 66420 : term1=delag*gpq(nu)+delbg*gpq(mu)
2787 : 66420 : term2=gpq(mu)*gpq(nu)*gpqdq
2788 : 66420 : term3=fac2*term2
2789 : 66420 : term2=two*term2/gsq
2790 : 66420 : gterms=term1-term2-term3
2791 : 66420 : dyewdq(re,mu,ia,nu,ib,iq)=dyewdq(re,mu,ia,nu,ib,iq)+gterms*c1r
2792 : 99630 : dyewdq(im,mu,ia,nu,ib,iq)=dyewdq(im,mu,ia,nu,ib,iq)+gterms*c1i
2793 : : end do
2794 : : end do
2795 : : end do
2796 : :
2797 : : end do
2798 : : end do
2799 : : end if
2800 : : ! Endif g/=0 :
2801 : : end if
2802 : : end if
2803 : : ! End triple loop over G s:
2804 : : end do
2805 : : end do
2806 : : end do
2807 : :
2808 : : ! Check if new shell must be calculated
2809 [ + + ]: 18 : if (newg==0) exit
2810 : : end do ! End the loop on ng (new shells). Note that there is one exit from this loop.
2811 : :
2812 : : !End G summation by accounting for some common factors.
2813 : : !(for the charges:see end of routine)
2814 : 3 : norm1=4.0_dp*pi/ucvol
2815 [ + + ]: 9 : do ia0=1,my_natom
2816 [ - + ]: 6 : ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
2817 [ + + ]: 18 : do ib=1,ia
2818 [ + + ]: 42 : do iq=1,3
2819 [ + + ]: 117 : do mu=1,3
2820 [ + + ]: 270 : do nu=1,mu
2821 [ + + ]: 567 : dyewdq(:,mu,ia,nu,ib,iq)=dyewdq(:,mu,ia,nu,ib,iq)*norm1
2822 : : end do
2823 : : end do
2824 : : end do
2825 : : end do
2826 : : end do
2827 : :
2828 : : !Do sums over real space:
2829 : 3 : reta=sqrt(eta)
2830 : 3 : reta3m=-eta*reta
2831 : 3 : fac=4._dp/3.0_dp/sqrt(pi)
2832 : 3 : nr=0
2833 : : do
2834 : 18 : nr=nr+1
2835 : 18 : newr=0
2836 : :
2837 [ + + ]: 162 : do ir3=-nr,nr
2838 [ + + ]: 1524 : do ir2=-nr,nr
2839 [ + + ]: 15762 : do ir1=-nr,nr
2840 [ + + + + : 15618 : if( abs(ir3)==nr .or. abs(ir2)==nr .or. abs(ir1)==nr .or. nr==1 )then
+ + + + ]
2841 : :
2842 : 6591 : arg=two_pi*(qphon(1)*ir1+qphon(2)*ir2+qphon(3)*ir3)
2843 : 6591 : c1r=cos(arg)*reta3m
2844 : 6591 : c1i=sin(arg)*reta3m
2845 [ + + ]: 19773 : do ia0=1,my_natom
2846 [ - + ]: 13182 : ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
2847 [ + + ]: 47211 : do ib=1,ia
2848 : 19773 : r1=dble(ir1)+xred(1,ib)-xred(1,ia)
2849 : 19773 : r2=dble(ir2)+xred(2,ib)-xred(2,ia)
2850 : 19773 : r3=dble(ir3)+xred(3,ib)-xred(3,ia)
2851 [ + + ]: 79092 : dakk(:)=two_pi*(/r1,r2,r3/)
2852 : 19773 : rdot12=rmet(2,1)*r1*r2
2853 : 19773 : rdot13=rmet(3,1)*r1*r3
2854 : 19773 : rdot23=rmet(3,2)*r2*r3
2855 : 19773 : dotr1=rmet(1,1)*r1**2+rdot12+rdot13
2856 : 19773 : dotr2=rmet(2,2)*r2**2+rdot12+rdot23
2857 : 19773 : dotr3=rmet(3,3)*r3**2+rdot13+rdot23
2858 : 19773 : rsq=dotr1+dotr2+dotr3
2859 : 19773 : rmagn=sqrt(rsq)
2860 : : ! Avoid zero denominators in term :
2861 [ + + ]: 32955 : if (rmagn>=1.0d-12) then
2862 : 19767 : arg=reta*rmagn
2863 : 19767 : term=zero
2864 [ + + ]: 19767 : if (arg<8.0_dp) then
2865 : : ! Note: erfc(8) is about 1.1e-29,
2866 : : ! so don t bother with larger arg.
2867 : : ! Also: exp(-64) is about 1.6e-28,
2868 : : ! so don t bother with larger arg**2 in exp.
2869 : 4200 : newr=1
2870 : 4200 : derfc_arg = abi_derfc(arg)
2871 : 4200 : term=derfc_arg/arg**3
2872 : 4200 : term1=2.0_dp/sqrt(pi)*exp(-arg**2)/arg**2
2873 : 4200 : term2=-(term+term1)
2874 : 4200 : term3=(3*term+term1*(3.0_dp+2.0_dp*arg**2))/rsq
2875 : 4200 : rq(1)=rmet(1,1)*r1+rmet(1,2)*r2+rmet(1,3)*r3
2876 : 4200 : rq(2)=rmet(2,1)*r1+rmet(2,2)*r2+rmet(2,3)*r3
2877 : 4200 : rq(3)=rmet(3,1)*r1+rmet(3,2)*r2+rmet(3,3)*r3
2878 [ + + ]: 16800 : do iq=1,3
2879 [ + + ]: 54600 : do mu=1,3
2880 : : ! do nu=1,3
2881 [ + + ]: 126000 : do nu=1,mu
2882 : : dyewdq(re,mu,ia,nu,ib,iq)=dyewdq(re,mu,ia,nu,ib,iq)-&
2883 : 75600 : & c1i*dakk(iq)*(rq(mu)*rq(nu)*term3+rmet(mu,nu)*term2)
2884 : : dyewdq(im,mu,ia,nu,ib,iq)=dyewdq(im,mu,ia,nu,ib,iq)+&
2885 : 113400 : & c1r*dakk(iq)*(rq(mu)*rq(nu)*term3+rmet(mu,nu)*term2)
2886 : : end do
2887 : : end do
2888 : : end do
2889 : : end if
2890 : : else
2891 [ + - ]: 6 : if (ia/=ib)then
2892 : : write(message,'(a,a,a,a,a,i5,a,i5,a)')&
2893 : 0 : & 'The distance between two atoms vanishes.',ch10,&
2894 : 0 : & 'This is not allowed.',ch10,&
2895 : 0 : & 'Action: check the input for the atoms number',ia,' and',ib,'.'
2896 : 0 : ABI_ERROR(message)
2897 : : end if
2898 : : end if
2899 : :
2900 : : end do ! End loop over ib:
2901 : : end do ! End loop over ia:
2902 : : end if
2903 : : end do ! End triple loop over real space points:
2904 : : end do
2905 : : end do
2906 : :
2907 : : ! Check if new shell must be calculated
2908 [ + + ]: 18 : if(newr==0) exit
2909 : : end do ! End loop on nr (new shells). Note that there is an exit within the loop
2910 : :
2911 : : !Take account of the charges
2912 : : !write(std_out,*)' '
2913 [ + + ]: 9 : do ia0=1,my_natom
2914 [ - + ]: 6 : ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
2915 [ + + ]: 18 : do ib=1,ia
2916 [ + + ]: 42 : do iq=1,3
2917 [ + + ]: 117 : do mu=1,3
2918 [ + + ]: 270 : do nu=1,mu
2919 [ + + ]: 567 : do ii=1,2
2920 : : dyewdq(ii,mu,ia,nu,ib,iq)=dyewdq(ii,mu,ia,nu,ib,iq)*&
2921 : 486 : & zion(typat(ia))*zion(typat(ib))
2922 : : end do
2923 : : end do
2924 : : end do
2925 : : end do
2926 : : end do
2927 : : end do
2928 : :
2929 : : !Symmetrize with respect to the directions
2930 [ + + ]: 9 : do ia0=1,my_natom
2931 [ - + ]: 6 : ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
2932 [ + + ]: 18 : do ib=1,ia
2933 [ + + ]: 42 : do iq=1,3
2934 [ + + ]: 117 : do mu=1,3
2935 [ + + ]: 270 : do nu=1,mu
2936 : 162 : dyewdq(re,nu,ia,mu,ib,iq)=dyewdq(re,mu,ia,nu,ib,iq)
2937 : 243 : dyewdq(im,nu,ia,mu,ib,iq)=dyewdq(im,mu,ia,nu,ib,iq)
2938 : : end do
2939 : : end do
2940 : : end do
2941 : : end do
2942 : : end do
2943 : :
2944 : : !In case of parallelism over atoms: communicate
2945 [ - + ]: 3 : if (paral_atom) then
2946 : 0 : call timab(48,1,tsec)
2947 : 0 : call xmpi_sum(dyewdq,my_comm_atom,ierr)
2948 : 0 : call timab(48,2,tsec)
2949 : : end if
2950 : :
2951 : : !Fill the upper part of the matrix, with the hermitian conjugate
2952 [ + + ]: 9 : do ia=1,natom
2953 [ + + ]: 18 : do ib=1,ia
2954 [ + + ]: 42 : do iq=1,3
2955 [ + + ]: 117 : do nu=1,3
2956 [ + + ]: 351 : do mu=1,3
2957 : 243 : dyewdq(re,mu,ib,nu,ia,iq)=dyewdq(re,mu,ia,nu,ib,iq)
2958 : 324 : dyewdq(im,mu,ib,nu,ia,iq)=-dyewdq(im,mu,ia,nu,ib,iq)
2959 : : end do
2960 : : end do
2961 : : end do
2962 : : end do
2963 : : end do
2964 : :
2965 : : !Destroy atom table used for parallelism
2966 : 3 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
2967 : :
2968 [ + - ]: 3 : end subroutine dfpt_ewalddq
2969 : : !!***
2970 : :
2971 : : !!****f* ABINIT/dfpt_ewalddqdq
2972 : : !!
2973 : : !! NAME
2974 : : !! dfpt_ewalddqdq
2975 : : !!
2976 : : !! FUNCTION
2977 : : !! Compute the second q-gradient of Ewald contribution to the dynamical matrix, at a given q wavevector,
2978 : : !! sumed over the second atomic sublattice.
2979 : : !! If q=0 is asked, sumg0 should be put to 0. Otherwise, it should be put to 1.
2980 : : !!
2981 : : !! COPYRIGHT
2982 : : !! Copyright (C) 1998-2026 ABINIT group (MR, MS)
2983 : : !! This file is distributed under the terms of the
2984 : : !! GNU General Public License, see ~abinit/COPYING
2985 : : !! or http://www.gnu.org/copyleft/gpl.txt .
2986 : : !!
2987 : : !! INPUTS
2988 : : !! gmet(3,3)=metric tensor in reciprocal space (length units **-2)
2989 : : !! mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
2990 : : !! comm_atom=--optional-- MPI communicator over atoms
2991 : : !! my_natom=number of atoms treated by current processor
2992 : : !! natom=number of atoms in unit cell
2993 : : !! qphon(3)=phonon wavevector (same system of coordinates as the
2994 : : !! reciprocal lattice vectors)
2995 : : !! rmet(3,3)=metric tensor in real space (length units squared)
2996 : : !! sumg0: if=1, the sum in reciprocal space must include g=0,
2997 : : !! if=0, this contribution must be skipped (q=0 singularity)
2998 : : !! typat(natom)=integer label of each type of atom (1,2,...)
2999 : : !! ucvol=unit cell volume in (whatever length scale units)**3
3000 : : !! xred(3,natom)=relative coords of atoms in unit cell (dimensionless)
3001 : : !! zion(ntypat)=charge on each type of atom (real number)
3002 : : !!
3003 : : !! OUTPUT
3004 : : !! dyewdqdq(2,3,natom,3,3,3)= First q-gradient of Ewald part of the dynamical matrix,
3005 : : !! sumed over second atomic sublattice.
3006 : : !!
3007 : : !! SOURCE
3008 : :
3009 : 6 : subroutine dfpt_ewalddqdq(dyewdqdq,gmet,my_natom,natom,qphon,rmet,sumg0,typat,ucvol,xred,zion, &
3010 [ - + - - : 3 : & mpi_atmtab,comm_atom ) ! optional arguments (parallelism))
- - ]
3011 : :
3012 : : !Arguments -------------------------------
3013 : : !scalars
3014 : : integer,intent(in) :: my_natom,natom,sumg0
3015 : : real(dp),intent(in) :: ucvol
3016 : : !arrays
3017 : : integer,intent(in) :: typat(natom)
3018 : : integer,optional,intent(in) :: comm_atom
3019 : : integer,optional,target,intent(in) :: mpi_atmtab(:)
3020 : : real(dp),intent(in) :: gmet(3,3),qphon(3),rmet(3,3),xred(3,natom),zion(*)
3021 : : real(dp),intent(out) :: dyewdqdq(2,3,natom,3,3,3)
3022 : :
3023 : : !Local variables -------------------------
3024 : : !nr, ng affect convergence of sums (nr=3,ng=5 is not good enough):
3025 : : !scalars
3026 : : integer,parameter :: im=2,nng=10,nnr=6,re=1
3027 : : integer :: ia,ia0,ib,ierr,ig1,ig2,ig3,ii,iq1,iq2,ir1,ir2,ir3,mu,my_comm_atom,newg,newr,ng,nr,nu
3028 : : logical :: my_atmtab_allocated,paral_atom
3029 : : real(dp) :: arg,arga,argb,c1i,c1r,delad,delag,delbd,delbg,derfc_arg
3030 : : real(dp) :: direct,dot1,dot2,dot3,dotr1,dotr2,dotr3
3031 : : real(dp) :: eta,fac,fac2,fac8,fac2sqr,gdot12,gdot13,gdot23,gsq,gsqsq,gpqdq1,gpqdq2,gterms,g0term,norm1
3032 : : real(dp) :: r1,r2,r3,rdot12,rdot13,rdot23,recip,reta
3033 : : real(dp) :: reta3m,rmagn,rsq,term,term1,term2,term3
3034 : : character(len=500) :: message
3035 : : !arrays
3036 : 3 : integer,pointer :: my_atmtab(:)
3037 : : real(dp) :: dakk(3),gpq(3),rq(3)
3038 : : real(dp) :: tsec(2)
3039 : 3 : real(dp),allocatable :: work(:,:,:,:,:,:,:)
3040 : :
3041 : : ! *************************************************************************
3042 : :
3043 : : !Set up parallelism over atoms
3044 [ + - ]: 3 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
3045 [ - + ]: 3 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
3046 [ + - ]: 3 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
3047 : 3 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
3048 : :
3049 : : !Compute eta for approximately optimized summations:
3050 : : direct=rmet(1,1)+rmet(1,2)+rmet(1,3)+rmet(2,1)+&
3051 : 3 : & rmet(2,2)+rmet(2,3)+rmet(3,1)+rmet(3,2)+rmet(3,3)
3052 : : recip=gmet(1,1)+gmet(1,2)+gmet(1,3)+gmet(2,1)+&
3053 : 3 : & gmet(2,2)+gmet(2,3)+gmet(3,1)+gmet(3,2)+gmet(3,3)
3054 : 3 : eta=pi*(dble(nng)/dble(nnr))*sqrt(1.69_dp*recip/direct)
3055 : : ! eta=1.0_dp
3056 : :
3057 : : !Test Ewald s summation
3058 : : !eta=1.2_dp*eta
3059 : :
3060 : : !Sum terms over g space:
3061 : 3 : fac=pi**2.0_dp/eta
3062 : 3 : fac2=2.0_dp*fac
3063 : 3 : fac8=4.0_dp*fac2
3064 : 3 : fac2sqr=fac2*fac2
3065 [ - + + - : 21 : ABI_MALLOC(work,(2,3,natom,3,natom,3,3))
+ - - + +
- + - + -
+ - + - -
+ - + ]
3066 [ + + + + : 3495 : work(:,:,:,:,:,:,:)=zero
+ + + + +
+ + + +
+ ]
3067 : : ng=0
3068 : : do
3069 : 18 : ng=ng+1
3070 : 18 : newg=0
3071 : :
3072 [ + + ]: 162 : do ig3=-ng,ng
3073 [ + + ]: 1524 : do ig2=-ng,ng
3074 [ + + ]: 15762 : do ig1=-ng,ng
3075 : :
3076 : : ! Exclude shells previously summed over
3077 : : if(abs(ig1)==ng .or. abs(ig2)==ng .or. abs(ig3)==ng&
3078 [ + + + + : 15618 : & .or. ng==1 ) then
+ + + + ]
3079 : :
3080 : 6591 : gpq(1)=dble(ig1)+qphon(1)
3081 : 6591 : gpq(2)=dble(ig2)+qphon(2)
3082 : 6591 : gpq(3)=dble(ig3)+qphon(3)
3083 : 6591 : gdot12=gmet(2,1)*gpq(1)*gpq(2)
3084 : 6591 : gdot13=gmet(3,1)*gpq(1)*gpq(3)
3085 : 6591 : gdot23=gmet(3,2)*gpq(2)*gpq(3)
3086 : 6591 : dot1=gmet(1,1)*gpq(1)**2+gdot12+gdot13
3087 : 6591 : dot2=gmet(2,2)*gpq(2)**2+gdot12+gdot23
3088 : 6591 : dot3=gmet(3,3)*gpq(3)**2+gdot13+gdot23
3089 : 6591 : gsq=dot1+dot2+dot3
3090 : 6591 : gsqsq=gsq*gsq
3091 : : ! Skip q=0:
3092 [ + + ]: 6591 : if (gsq<1.0d-20) then
3093 : :
3094 : : ! At second order in q there is a nonvanishing G=0 contribution in the longwave limit
3095 [ + - ]: 3 : if (sumg0==1) then
3096 [ + + ]: 9 : do ia0=1,my_natom
3097 [ - + ]: 6 : ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
3098 [ + + ]: 18 : do ib=1,ia
3099 [ + + ]: 42 : do iq2=1,3
3100 [ + + ]: 117 : do iq1=1,3
3101 [ + + ]: 351 : do mu=1,3
3102 [ + + ]: 243 : delag=zero; if(iq1==mu) delag=one
3103 [ + + ]: 243 : delad=zero; if(iq2==mu) delad=one
3104 [ + + ]: 810 : do nu=1,mu
3105 [ + + ]: 486 : delbg=zero; if(iq1==nu) delbg=one
3106 [ + + ]: 486 : delbd=zero; if(iq2==nu) delbd=one
3107 : 486 : g0term=-fac*(delad*delbg+delbd*delag)
3108 : 729 : work(re,mu,ia,nu,ib,iq1,iq2)=work(re,mu,ia,nu,ib,iq1,iq2)+g0term
3109 : : end do
3110 : : end do
3111 : : end do
3112 : : end do
3113 : : end do
3114 : : end do
3115 : : end if
3116 : :
3117 : : else
3118 : 6588 : arg=fac*gsq
3119 : : ! Larger arg gives 0 contribution:
3120 [ + + ]: 6588 : if (arg <= 80._dp) then
3121 : 1230 : newg=1
3122 : 1230 : term=exp(-arg)/gsq
3123 [ + + ]: 3690 : do ia0=1,my_natom
3124 [ - + ]: 2460 : ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
3125 : 2460 : arga=two_pi*(gpq(1)*xred(1,ia)+gpq(2)*xred(2,ia)+gpq(3)*xred(3,ia))
3126 [ + + ]: 7380 : do ib=1,ia
3127 : 3690 : argb=two_pi*(gpq(1)*xred(1,ib)+gpq(2)*xred(2,ib)+gpq(3)*xred(3,ib))
3128 : 3690 : arg=arga-argb
3129 : 3690 : c1r=cos(arg)*term
3130 : 3690 : c1i=sin(arg)*term
3131 : :
3132 [ + + ]: 17220 : do iq2=1,3
3133 : 11070 : gpqdq2=gmet(iq2,1)*gpq(1)+gmet(iq2,2)*gpq(2)+gmet(iq2,3)*gpq(3)
3134 [ + + ]: 47970 : do iq1=1,3
3135 : 33210 : gpqdq1=gmet(iq1,1)*gpq(1)+gmet(iq1,2)*gpq(2)+gmet(iq1,3)*gpq(3)
3136 [ + + ]: 143910 : do mu=1,3
3137 [ + + ]: 99630 : delag=zero; if(iq1==mu) delag=one
3138 [ + + ]: 99630 : delad=zero; if(iq2==mu) delad=one
3139 [ + + ]: 332100 : do nu=1,mu
3140 [ + + ]: 199260 : delbg=zero; if(iq1==nu) delbg=one
3141 [ + + ]: 199260 : delbd=zero; if(iq2==nu) delbd=one
3142 : :
3143 : 199260 : term1=gpqdq2*(delag*gpq(nu)+delbg*gpq(mu))
3144 : 199260 : term1=term1+gpqdq1*(delad*gpq(nu)+delbd*gpq(mu))
3145 : 199260 : term1=term1+gpq(mu)*gpq(nu)*gmet(iq1,iq2)
3146 : 199260 : term1=-term1*(fac2+2.0_dp/gsq)
3147 : :
3148 : 199260 : term2=delag*delbd + delbg*delad
3149 : :
3150 : 199260 : term3=gpqdq1*gpqdq2*gpq(mu)*gpq(nu)
3151 : 199260 : term3=term3*(fac8/gsq + fac2sqr + 8.0_dp/gsqsq)
3152 : :
3153 : 199260 : gterms=term1+term2+term3
3154 : 199260 : work(re,mu,ia,nu,ib,iq1,iq2)=work(re,mu,ia,nu,ib,iq1,iq2)+gterms*c1r
3155 : 298890 : work(im,mu,ia,nu,ib,iq1,iq2)=work(im,mu,ia,nu,ib,iq1,iq2)+gterms*c1i
3156 : : end do
3157 : : end do
3158 : : end do
3159 : : end do
3160 : : end do
3161 : : end do
3162 : : end if
3163 : : ! Endif g/=0 :
3164 : : end if
3165 : : end if
3166 : : ! End triple loop over G s:
3167 : : end do
3168 : : end do
3169 : : end do
3170 : :
3171 : : ! Check if new shell must be calculated
3172 [ + + ]: 18 : if (newg==0) exit
3173 : : end do ! End the loop on ng (new shells). Note that there is one exit from this loop.
3174 : :
3175 : : !End G summation by accounting for some common factors.
3176 : : !(for the charges:see end of routine)
3177 : 3 : norm1=4.0_dp*pi/ucvol
3178 [ + + ]: 9 : do ia0=1,my_natom
3179 [ - + ]: 6 : ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
3180 [ + + ]: 18 : do ib=1,ia
3181 [ + + ]: 42 : do iq2=1,3
3182 [ + + ]: 117 : do iq1=1,3
3183 [ + + ]: 351 : do mu=1,3
3184 [ + + ]: 810 : do nu=1,mu
3185 [ + + ]: 1701 : work(:,mu,ia,nu,ib,iq1,iq2)=work(:,mu,ia,nu,ib,iq1,iq2)*norm1
3186 : : end do
3187 : : end do
3188 : : end do
3189 : : end do
3190 : : end do
3191 : : end do
3192 : :
3193 : : !Do sums over real space:
3194 : 3 : reta=sqrt(eta)
3195 : 3 : reta3m=eta*reta
3196 : 3 : fac=4._dp/3.0_dp/sqrt(pi)
3197 : 3 : nr=0
3198 : : do
3199 : 18 : nr=nr+1
3200 : 18 : newr=0
3201 : :
3202 [ + + ]: 162 : do ir3=-nr,nr
3203 [ + + ]: 1524 : do ir2=-nr,nr
3204 [ + + ]: 15762 : do ir1=-nr,nr
3205 [ + + + + : 15618 : if( abs(ir3)==nr .or. abs(ir2)==nr .or. abs(ir1)==nr .or. nr==1 )then
+ + + + ]
3206 : :
3207 : 6591 : arg=two_pi*(qphon(1)*ir1+qphon(2)*ir2+qphon(3)*ir3)
3208 : 6591 : c1r=cos(arg)*reta3m
3209 : 6591 : c1i=sin(arg)*reta3m
3210 [ + + ]: 19773 : do ia0=1,my_natom
3211 [ - + ]: 13182 : ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
3212 [ + + ]: 47211 : do ib=1,ia
3213 : 19773 : r1=dble(ir1)+xred(1,ib)-xred(1,ia)
3214 : 19773 : r2=dble(ir2)+xred(2,ib)-xred(2,ia)
3215 : 19773 : r3=dble(ir3)+xred(3,ib)-xred(3,ia)
3216 [ + + ]: 79092 : dakk(:)=two_pi*(/r1,r2,r3/)
3217 : 19773 : rdot12=rmet(2,1)*r1*r2
3218 : 19773 : rdot13=rmet(3,1)*r1*r3
3219 : 19773 : rdot23=rmet(3,2)*r2*r3
3220 : 19773 : dotr1=rmet(1,1)*r1**2+rdot12+rdot13
3221 : 19773 : dotr2=rmet(2,2)*r2**2+rdot12+rdot23
3222 : 19773 : dotr3=rmet(3,3)*r3**2+rdot13+rdot23
3223 : 19773 : rsq=dotr1+dotr2+dotr3
3224 : 19773 : rmagn=sqrt(rsq)
3225 : : ! Avoid zero denominators in term :
3226 [ + + ]: 32955 : if (rmagn>=1.0d-12) then
3227 : 19767 : arg=reta*rmagn
3228 : 19767 : term=zero
3229 [ + + ]: 19767 : if (arg<8.0_dp) then
3230 : : ! Note: erfc(8) is about 1.1e-29,
3231 : : ! so don t bother with larger arg.
3232 : : ! Also: exp(-64) is about 1.6e-28,
3233 : : ! so don t bother with larger arg**2 in exp.
3234 : 4200 : newr=1
3235 : 4200 : derfc_arg = abi_derfc(arg)
3236 : 4200 : term=derfc_arg/arg**3
3237 : 4200 : term1=2.0_dp/sqrt(pi)*exp(-arg**2)/arg**2
3238 : 4200 : term2=-(term+term1)
3239 : 4200 : term3=(3*term+term1*(3.0_dp+2.0_dp*arg**2))/rsq
3240 : 4200 : rq(1)=rmet(1,1)*r1+rmet(1,2)*r2+rmet(1,3)*r3
3241 : 4200 : rq(2)=rmet(2,1)*r1+rmet(2,2)*r2+rmet(2,3)*r3
3242 : 4200 : rq(3)=rmet(3,1)*r1+rmet(3,2)*r2+rmet(3,3)*r3
3243 [ + + ]: 16800 : do iq2=1,3
3244 [ + + ]: 54600 : do iq1=1,3
3245 [ + + ]: 163800 : do mu=1,3
3246 : : ! do nu=1,3
3247 [ + + ]: 378000 : do nu=1,mu
3248 : : work(re,mu,ia,nu,ib,iq1,iq2)=work(re,mu,ia,nu,ib,iq1,iq2)+&
3249 : 226800 : & c1r*dakk(iq1)*dakk(iq2)*(rq(mu)*rq(nu)*term3+rmet(mu,nu)*term2)
3250 : : work(im,mu,ia,nu,ib,iq1,iq2)=work(im,mu,ia,nu,ib,iq1,iq2)+&
3251 : 340200 : & c1i*dakk(iq1)*dakk(iq2)*(rq(mu)*rq(nu)*term3+rmet(mu,nu)*term2)
3252 : : end do
3253 : : end do
3254 : : end do
3255 : : end do
3256 : : end if
3257 : : else
3258 [ + - ]: 6 : if (ia/=ib)then
3259 : : write(message,'(a,a,a,a,a,i5,a,i5,a)')&
3260 : 0 : & 'The distance between two atoms vanishes.',ch10,&
3261 : 0 : & 'This is not allowed.',ch10,&
3262 : 0 : & 'Action: check the input for the atoms number',ia,' and',ib,'.'
3263 : 0 : ABI_ERROR(message)
3264 : : end if
3265 : : end if
3266 : :
3267 : : end do ! End loop over ib:
3268 : : end do ! End loop over ia:
3269 : : end if
3270 : : end do ! End triple loop over real space points:
3271 : : end do
3272 : : end do
3273 : :
3274 : : ! Check if new shell must be calculated
3275 [ + + ]: 18 : if(newr==0) exit
3276 : : end do ! End loop on nr (new shells). Note that there is an exit within the loop
3277 : :
3278 : : !Take account of the charges
3279 : : !write(std_out,*)' '
3280 [ + + ]: 9 : do ia0=1,my_natom
3281 [ - + ]: 6 : ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
3282 [ + + ]: 18 : do ib=1,ia
3283 [ + + ]: 42 : do iq2=1,3
3284 [ + + ]: 117 : do iq1=1,3
3285 [ + + ]: 351 : do mu=1,3
3286 [ + + ]: 810 : do nu=1,mu
3287 [ + + ]: 1701 : do ii=1,2
3288 : : work(ii,mu,ia,nu,ib,iq1,iq2)=work(ii,mu,ia,nu,ib,iq1,iq2)*&
3289 : 1458 : & zion(typat(ia))*zion(typat(ib))
3290 : : end do
3291 : : end do
3292 : : end do
3293 : : end do
3294 : : end do
3295 : : end do
3296 : : end do
3297 : :
3298 : : !Symmetrize with respect to the directions
3299 [ + + ]: 9 : do ia0=1,my_natom
3300 [ - + ]: 6 : ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
3301 [ + + ]: 18 : do ib=1,ia
3302 [ + + ]: 42 : do iq2=1,3
3303 [ + + ]: 117 : do iq1=1,3
3304 [ + + ]: 351 : do mu=1,3
3305 [ + + ]: 810 : do nu=1,mu
3306 : 486 : work(re,nu,ia,mu,ib,iq1,iq2)=work(re,mu,ia,nu,ib,iq1,iq2)
3307 : 729 : work(im,nu,ia,mu,ib,iq1,iq2)=work(im,mu,ia,nu,ib,iq1,iq2)
3308 : : end do
3309 : : end do
3310 : : end do
3311 : : end do
3312 : : end do
3313 : : end do
3314 : :
3315 : : !In case of parallelism over atoms: communicate
3316 [ - + ]: 3 : if (paral_atom) then
3317 : 0 : call timab(48,1,tsec)
3318 : 0 : call xmpi_sum(work,my_comm_atom,ierr)
3319 : 0 : call timab(48,2,tsec)
3320 : : end if
3321 : :
3322 : : !Fill the upper part of the matrix, with the hermitian conjugate
3323 [ + + ]: 9 : do ia=1,natom
3324 [ + + ]: 18 : do ib=1,ia
3325 [ + + ]: 42 : do iq2=1,3
3326 [ + + ]: 117 : do iq1=1,3
3327 [ + + ]: 351 : do nu=1,3
3328 [ + + ]: 1053 : do mu=1,3
3329 : 729 : work(re,mu,ib,nu,ia,iq1,iq2)=work(re,mu,ia,nu,ib,iq1,iq2)
3330 : 972 : work(im,mu,ib,nu,ia,iq1,iq2)=-work(im,mu,ia,nu,ib,iq1,iq2)
3331 : : end do
3332 : : end do
3333 : : end do
3334 : : end do
3335 : : end do
3336 : : end do
3337 : :
3338 : : !Perform the summation over the second atomic sublattice
3339 [ + + + + : 1740 : dyewdqdq(:,:,:,:,:,:)=zero
+ + + + +
+ + + ]
3340 [ + + ]: 9 : do ia=1,natom
3341 [ + + ]: 27 : do iq2=1,3
3342 [ + + ]: 78 : do iq1=1,3
3343 [ + + ]: 234 : do nu=1,3
3344 [ + + ]: 702 : do mu=1,3
3345 [ + + ]: 1620 : do ib=1,natom
3346 : : dyewdqdq(re,mu,ia,nu,iq1,iq2)=dyewdqdq(re,mu,ia,nu,iq1,iq2) + &
3347 : 972 : & work(re,mu,ia,nu,ib,iq1,iq2)
3348 : : dyewdqdq(im,mu,ia,nu,iq1,iq2)=dyewdqdq(im,mu,ia,nu,iq1,iq2) + &
3349 : 1458 : & work(im,mu,ia,nu,ib,iq1,iq2)
3350 : : end do
3351 : : end do
3352 : : end do
3353 : : end do
3354 : : end do
3355 : : end do
3356 [ - + ]: 3 : ABI_FREE(work)
3357 : :
3358 : : !Destroy atom table used for parallelism
3359 : 3 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
3360 : :
3361 [ + - ]: 3 : end subroutine dfpt_ewalddqdq
3362 : : !!***
3363 : :
3364 : : end module m_dfpt_elt
3365 : : !!***
|
