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 : !!***
|
