Branch data Line data Source code
1 : : !!****m* ABINIT/m_ewald
2 : : !! NAME
3 : : !! m_ewald
4 : : !!
5 : : !! FUNCTION
6 : : !! This module gathers routines to compute the Ewald energy and its derivatives
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 2014-2026 ABINIT group (DCA, XG, JJC, GMR)
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_ewald
23 : :
24 : : use defs_basis
25 : : use m_abicore
26 : : use m_errors
27 : : use m_splines
28 : : use m_time
29 : : use m_xmpi
30 : :
31 : : use m_gtermcutoff, only : termcutoff
32 : : use m_special_funcs, only : abi_derfc
33 : : use m_matrix, only : matr3inv
34 : :
35 : : implicit none
36 : :
37 : : private
38 : :
39 : : public :: ewald ! Compute Ewald energy and derivatives with respect to xred
40 : : public :: ewald2 ! Derivative of the Ewald energy with respect to strain.
41 : : public :: ewald9 ! Compute ewald contribution to the dynamical matrix, at a given
42 : : ! q wavevector, including anisotropic dielectric tensor and effective charges
43 : : public :: ewald9_2D! Compute ewald contribution to the dynamical matrix, at a given
44 : : ! q wavevector, in the case of a 2D material with an external
45 : : ! dielectric environment
46 : :
47 : : contains
48 : : !!***
49 : :
50 : : !!****f* m_ewald/ewald
51 : : !!
52 : : !! NAME
53 : : !! ewald
54 : : !!
55 : : !! FUNCTION
56 : : !! Compute Ewald energy and derivatives with respect to dimensionless
57 : : !! reduced atom coordinates xred.
58 : : !!
59 : : !! INPUTS
60 : : !! gmet(3,3)=metric tensor in reciprocal space (bohr^-2)
61 : : !! natom=number of atoms in unit cell
62 : : !! ntypat=numbe of type of atoms
63 : : !! rmet(3,3)=metric tensor in real space (bohr^2)
64 : : !! typat(natom)=integer label of each type of atom (1,2,...)
65 : : !! ucvol=unit cell volume (bohr^3)
66 : : !! xred(3,natom)=relative coords of atoms in unit cell (dimensionless)
67 : : !! zion(ntypat)=charge on each type of atom (real number)
68 : : !!
69 : : !! OUTPUT
70 : : !! eew=final ewald energy in hartrees
71 : : !! grewtn(3,natom)=grads of eew wrt xred(3,natom), hartrees.
72 : : !!
73 : : !! SOURCE
74 : :
75 : 7686 : subroutine ewald(eew,gmet,grewtn,gsqcut,icutcoul,natom,ngfft,nkpt,ntypat,rcut,&
76 : 7686 : rmet,rprimd,typat,ucvol,vcutgeo,xred,zion)
77 : :
78 : : !Arguments ------------------------------------
79 : : !scalars
80 : : integer,intent(in) :: icutcoul,natom,nkpt,ntypat
81 : : real(dp),intent(in) :: gsqcut,rcut,ucvol
82 : : real(dp),intent(out) :: eew
83 : : !arrays
84 : : integer,intent(in) :: ngfft(18),typat(natom)
85 : : real(dp),intent(in) :: gmet(3,3),rmet(3,3),rprimd(3,3),xred(3,natom),vcutgeo(3),zion(ntypat)
86 : : real(dp),intent(out) :: grewtn(3,natom)
87 : :
88 : : !Local variables-------------------------------
89 : : !scalars
90 : : integer :: ia,ib,ig1,ig2,ig3,ig23,ii,ir1,ir2,ir3,newg,newr,ng,nr
91 : : real(dp) :: arg,c1i,ch,chsq,derfc_arg,direct,drdta1,drdta2,drdta3,eta,fac
92 : : real(dp) :: fraca1,fraca2,fraca3,fracb1,fracb2,fracb3,gsq,gsum,phi,phr,r1
93 : : real(dp) :: minexparg
94 : : real(dp) :: r1a1d,r2,r2a2d,r3,r3a3d,recip,reta,rmagn,rsq,sumg,summi,summr,sumr
95 : : real(dp) :: t1,term ,zcut !, gcart_para, gcart_perp
96 : : !character(len=500) :: msg
97 : : !arrays
98 : 7686 : real(dp),allocatable :: gcutoff(:)
99 : : ! *************************************************************************
100 : :
101 : : !This is the minimum argument of an exponential, with some safety
102 : 7686 : minexparg=log(tiny(0._dp))+five
103 : :
104 : : !Add up total charge and sum of $charge^2$ in cell
105 : :
106 : 7686 : chsq=0._dp
107 : 7686 : ch=0._dp
108 [ + + ]: 30112 : do ia=1,natom
109 : 22426 : ch=ch+zion(typat(ia))
110 : 30112 : chsq=chsq+zion(typat(ia))**2
111 : : end do
112 : :
113 : : !Compute eta, the Ewald summation convergence parameter,
114 : : ! for approximately optimized summations:
115 : : direct=rmet(1,1)+rmet(1,2)+rmet(1,3)+rmet(2,1)+&
116 : 7686 : rmet(2,2)+rmet(2,3)+rmet(3,1)+rmet(3,2)+rmet(3,3)
117 : : recip=gmet(1,1)+gmet(1,2)+gmet(1,3)+gmet(2,1)+&
118 : 7686 : gmet(2,2)+gmet(2,3)+gmet(3,1)+gmet(3,2)+gmet(3,3)
119 : :
120 : : !A bias is introduced, because G-space summation scales
121 : : !better than r space summation ! Note: debugging is the most easier at fixed eta.
122 : : zcut=SQRT(DOT_PRODUCT(rprimd(:,3),rprimd(:,3)))/2.0_dp
123 [ - + ]: 7686 : if(icutcoul.eq.1) then
124 [ # # ]: 0 : eta=SQRT(16.0_dp/SQRT(DOT_PRODUCT(rprimd(:,1),rprimd(:,1))))
125 : : ! else if (icutcoul.eq.2) then
126 : : ! zcut=SQRT(DOT_PRODUCT(rprimd(:,3),rprimd(:,3)))/2.0_dp
127 : : ! eta=217.6_dp/zcut**2.0_dp
128 : : ! eta=1.0_dp/zcut**2.0_dp
129 : : ! eta=SQRT(16.0_dp/SQRT(DOT_PRODUCT(rprimd(:,1),rprimd(:,1))))
130 : : ! else if (icutcoul.eq.2) then
131 : : ! zcut=SQRT(DOT_PRODUCT(rprimd(:,3),rprimd(:,3)))/2.0_dp
132 : : ! eta=SQRT(8.0_dp/zcut)
133 : : ! eta=SQRT(16.0_dp/SQRT(DOT_PRODUCT(rprimd(:,1),rprimd(:,1))))
134 : : ! else if (icutcoul.eq.2) then
135 : : ! zcut=SQRT(DOT_PRODUCT(rprimd(:,3),rprimd(:,3)))/2.0_dp
136 : : ! eta=SQRT(8.0_dp/zcut)
137 : : else
138 : 7686 : eta=pi*200.0_dp/33.0_dp*sqrt(1.69_dp*recip/direct)
139 : : end if
140 : :
141 : : !Conduct reciprocal space summations
142 : 7686 : fac=pi**2/eta
143 : 7686 : gsum=0._dp
144 [ + + + + ]: 97390 : grewtn(:,:)=0.0_dp
145 : :
146 : : !Initialize Gcut-off array from m_gtermcutoff
147 : : !ABI_MALLOC(gcutoff,(ngfft(1)*ngfft(2)*ngfft(3)))
148 : 7686 : call termcutoff(gcutoff,gsqcut,icutcoul,ngfft,nkpt,rcut,rprimd,vcutgeo)
149 : :
150 : : !if (icutcoul.eq.3) then
151 : : !Sum over G space, done shell after shell until all
152 : : !contributions are too small.
153 : 7686 : ng=0
154 : : do
155 : 319904 : ng=ng+1
156 : 319904 : newg=0
157 : : ! Instead of this warning that most normal users do not understand (because they are doing GS calculations, and not RF calculations),
158 : : ! one should optimize this routine. But usually this is a very small fraction of any ABINIT run.
159 : : ! if (ng > 20 .and. mod(ng,10)==0) then
160 : : ! write (msg,'(3a,I10)') "Very large box of G neighbors in ewald: you probably do not want to do this.", ch10,&
161 : : !& " If you have a metal consider setting dipdip 0. ng = ", ng
162 : : ! ABI_WARNING(msg)
163 : : ! end if
164 : 319904 : ii=1
165 [ + + ]: 15947432 : do ig3=-ng,ng
166 [ + + ]: 1243477224 : do ig2=-ng,ng
167 [ + + ]: >18139*10^7 : do ig1=-ng,ng
168 : : ! Exclude shells previously summed over
169 [ + + + + : >18137*10^7 : if(abs(ig1)==ng .or. abs(ig2)==ng .or. abs(ig3)==ng .or. ng==1 ) then
+ + + + ]
170 : :
171 : : ! gsq is G dot G = |G|^2
172 : : gsq=gmet(1,1)*dble(ig1*ig1)+gmet(2,2)*dble(ig2*ig2)+&
173 : : & gmet(3,3)*dble(ig3*ig3)+2._dp*(gmet(2,1)*dble(ig1*ig2)+&
174 : 7180215334 : & gmet(3,1)*dble(ig1*ig3)+gmet(3,2)*dble(ig3*ig2))
175 : :
176 : : ! Skip g=0:
177 [ + + ]: 7180215334 : if (gsq>1.0d-20) then
178 : 7180207648 : arg=fac*gsq
179 : :
180 : : ! Larger arg gives 0 contribution because of exp(-arg)
181 [ + + ]: 7180207648 : if (arg <= -minexparg ) then
182 : : ! When any term contributes then include next shell
183 : 1847389240 : newg=1
184 : :
185 : : if((abs(ig1).lt.ngfft(1)).and.&
186 [ + + + + : 1847389240 : &(abs(ig2).lt.ngfft(2)).and.&
+ + ]
187 : : &(abs(ig3).lt.ngfft(3))) then
188 : 569930776 : ig23=ngfft(1)*(abs(ig2)+ngfft(2)*(abs(ig3)))
189 : 569930776 : ii=abs(ig1)+ig23+1
190 : : !term= ( exp(-arg) + gcutoff(ii) - 1.0_dp )/gsq
191 : : !term=exp(-arg)/gsq*gcutoff(ii)
192 : : !term= ( exp(-arg) + gcutoff(ii) - 1.0_dp)/gsq
193 : 569930776 : term=exp(-arg)/gsq*gcutoff(ii)
194 [ + + ]: 1277458464 : else if (icutcoul.ne.3) then
195 : : term=zero !exp(-arg)/gsq
196 : : else
197 : 1277236418 : term=exp(-arg)/gsq
198 : : endif
199 : :
200 : 1847389240 : summr = 0.0_dp
201 : 1847389240 : summi = 0.0_dp
202 : :
203 : :
204 : : ! XG 20180531 : the two do-loops on ia should be merged, in order to spare
205 : : ! the waste of computing twice the sin and cos.
206 : :
207 : : ! Note that if reduced atomic coordinates xred drift outside
208 : : ! of unit cell (outside [0,1)) it is irrelevant in the following
209 : : ! term, which only computes a phase.
210 [ + + ]: 7604674612 : do ia=1,natom
211 : 5757285372 : arg=two_pi*(ig1*xred(1,ia)+ig2*xred(2,ia)+ig3*xred(3,ia))
212 : : ! Sum real and imaginary parts (avoid complex variables)
213 : 5757285372 : summr=summr+zion(typat(ia))*cos(arg)
214 : 7604674612 : summi=summi+zion(typat(ia))*sin(arg)
215 : : end do
216 : :
217 : : ! The following two checks avoid an annoying underflow error msg
218 [ + + ]: 1847389240 : if (abs(summr)<1.d-16) summr=0.0_dp
219 [ + + ]: 1847389240 : if (abs(summi)<1.d-16) summi=0.0_dp
220 : :
221 : : ! The product of term and summr**2 or summi**2 below
222 : : ! can underflow if not for checks above
223 : 1847389240 : t1=term*(summr*summr+summi*summi)
224 : 1847389240 : gsum=gsum+t1
225 : :
226 [ + + ]: 7604674612 : do ia=1,natom
227 : : ! Again only phase is computed so xred may fall outside [0,1).
228 : 5757285372 : arg=two_pi*(ig1*xred(1,ia)+ig2*xred(2,ia)+ig3*xred(3,ia))
229 : 5757285372 : phr= cos(arg)
230 : 5757285372 : phi=-sin(arg)
231 : : ! (note: do not need real part, commented out)
232 : : ! c1r=(phr*summr-phi*summi)*(term*zion(typat(ia)))
233 : 5757285372 : c1i=(phi*summr+phr*summi)*(term*zion(typat(ia)))
234 : : ! compute coordinate gradients
235 : 5757285372 : grewtn(1,ia)=grewtn(1,ia)-c1i*ig1
236 : 5757285372 : grewtn(2,ia)=grewtn(2,ia)-c1i*ig2
237 : 7604674612 : grewtn(3,ia)=grewtn(3,ia)-c1i*ig3
238 : : end do
239 : :
240 : : end if ! End condition of not larger than -minexparg
241 : : end if ! End skip g=0
242 : : end if ! End triple loop over G s and associated new shell condition
243 : :
244 : : end do
245 : : end do
246 : : end do
247 : :
248 : : ! Check if new shell must be calculated
249 [ + + ]: 319904 : if (newg==0) exit
250 : :
251 : : end do ! End the loop on ng (new shells). Note that there is one exit from this loop.
252 : : !endif
253 : :
254 : 7686 : sumg=gsum/(two_pi*ucvol)
255 : :
256 : : !Stress tensor is now computed elsewhere (ewald2) hence do not need
257 : : !length scale gradients (used to compute them here).
258 : :
259 : : !normalize coordinate gradients by unit cell volume ucvol
260 : 7686 : term=-2._dp/ucvol
261 [ + + + + ]: 97390 : grewtn(:,:)=grewtn(:,:)*term
262 : : !call DSCAL(3*natom,term,grewtn,1)
263 : :
264 : : !Conduct real space summations
265 : 7686 : reta=sqrt(eta)
266 : 7686 : fac=2._dp*sqrt(eta/pi)
267 : 7686 : sumr=0.0_dp
268 : :
269 : : !In the following a summation is being conducted over all
270 : : !unit cells (ir1, ir2, ir3) so it is appropriate to map all
271 : : !reduced coordinates xred back into [0,1).
272 : : !
273 : : !Loop on shells in r-space as was done in g-space
274 : 7686 : nr=0
275 : : do
276 : 22367 : nr=nr+1
277 : 22367 : newr=0
278 : : ! Instead of this warning that most normal users do not understand (because they are doing GS calculations, and not RF calculations),
279 : : ! one should optimize this routine. But usually this is a very small fraction of any ABINIT run.
280 : : ! if (nr > 20 .and. mod(nr,10)==0) then
281 : : ! write (msg,'(3a,I10)') "Very large box of R neighbors in ewald: you probably do not want to do this.", ch10,&
282 : : !& " If you have a metal consider setting dipdip 0. nr = ", nr
283 : : ! ABI_WARNING(msg)
284 : : ! end if
285 : : !
286 [ + + ]: 135718 : do ir3=-nr,nr
287 [ + + ]: 787349 : do ir2=-nr,nr
288 [ + + ]: 4944269 : do ir1=-nr,nr
289 [ + + + + : 4830918 : if( abs(ir3)==nr .or. abs(ir2)==nr .or. abs(ir1)==nr .or. nr==1 )then
+ + + + ]
290 : :
291 [ + + ]: 11533170 : do ia=1,natom
292 : : ! Map reduced coordinate xred(mu,ia) into [0,1)
293 : 8796974 : fraca1=xred(1,ia)-aint(xred(1,ia))+0.5_dp-sign(0.5_dp,xred(1,ia))
294 : 8796974 : fraca2=xred(2,ia)-aint(xred(2,ia))+0.5_dp-sign(0.5_dp,xred(2,ia))
295 : 8796974 : fraca3=xred(3,ia)-aint(xred(3,ia))+0.5_dp-sign(0.5_dp,xred(3,ia))
296 : 8796974 : drdta1=0.0_dp
297 : 8796974 : drdta2=0.0_dp
298 : 8796974 : drdta3=0.0_dp
299 : :
300 [ + + ]: 103325178 : do ib=1,natom
301 : : ! fraca and fracb should be precomputedi and become arrays with natom dimension.
302 : : ! Also the combination with dble(ir1), dble(ir2), dble(ir3) or fraca should be done outside of the ib loop.
303 : 94528204 : fracb1=xred(1,ib)-aint(xred(1,ib))+0.5_dp-sign(0.5_dp,xred(1,ib))
304 : 94528204 : fracb2=xred(2,ib)-aint(xred(2,ib))+0.5_dp-sign(0.5_dp,xred(2,ib))
305 : 94528204 : fracb3=xred(3,ib)-aint(xred(3,ib))+0.5_dp-sign(0.5_dp,xred(3,ib))
306 : 94528204 : r1=dble(ir1)+fracb1-fraca1
307 : 94528204 : r2=dble(ir2)+fracb2-fraca2
308 : 94528204 : r3=dble(ir3)+fracb3-fraca3
309 : : rsq=rmet(1,1)*r1*r1+rmet(2,2)*r2*r2+rmet(3,3)*r3*r3+&
310 : 94528204 : & 2.0_dp*(rmet(2,1)*r2*r1+rmet(3,2)*r3*r2+rmet(3,1)*r1*r3)
311 : :
312 : : ! Avoid zero denominators in 'term':
313 [ + + ]: 103325178 : if (rsq>=1.0d-24) then
314 : :
315 : : ! Note: erfc(8) is about 1.1e-29, so do not bother with larger arg.
316 : : ! Also: exp(-64) is about 1.6e-28, so do not bother with larger arg**2 in exp.
317 : 94505778 : term=0._dp
318 [ + + ]: 94505778 : if (eta*rsq<64.0_dp) then
319 : 5563276 : newr=1
320 : 5563276 : rmagn=sqrt(rsq)
321 : 5563276 : arg=reta*rmagn
322 : : ! derfc is the real(dp) complementary error function
323 : 5563276 : derfc_arg = abi_derfc(arg)
324 : 5563276 : term=derfc_arg/rmagn
325 : 5563276 : sumr=sumr+zion(typat(ia))*zion(typat(ib))*term
326 : : term=zion(typat(ia))*zion(typat(ib))*&
327 : 5563276 : & (term+fac*exp(-eta*rsq))/rsq
328 : : ! Length scale grads now handled with stress tensor in ewald2
329 : 5563276 : r1a1d=rmet(1,1)*r1+rmet(1,2)*r2+rmet(1,3)*r3
330 : 5563276 : r2a2d=rmet(2,1)*r1+rmet(2,2)*r2+rmet(2,3)*r3
331 : 5563276 : r3a3d=rmet(3,1)*r1+rmet(3,2)*r2+rmet(3,3)*r3
332 : : ! Compute terms related to coordinate gradients
333 : 5563276 : drdta1=drdta1+term*r1a1d
334 : 5563276 : drdta2=drdta2+term*r2a2d
335 : 5563276 : drdta3=drdta3+term*r3a3d
336 : : end if
337 : : end if ! End avoid zero denominators in'term'
338 : : end do ! end loop over ib:
339 : :
340 : 8796974 : grewtn(1,ia)=grewtn(1,ia)+drdta1
341 : 8796974 : grewtn(2,ia)=grewtn(2,ia)+drdta2
342 : 12976261 : grewtn(3,ia)=grewtn(3,ia)+drdta3
343 : : end do ! end loop over ia:
344 : : end if
345 : : end do ! end triple loop over real space points and associated condition of new shell
346 : : end do
347 : : end do
348 : :
349 : : ! Check if new shell must be calculated
350 [ + + ]: 22367 : if(newr==0) exit
351 : : end do ! End loop on nr (new shells). Note that there is an exit within the loop
352 : : !
353 : 7686 : sumr=0.5_dp*sumr
354 : 7686 : fac=pi*ch**2.0_dp/(2.0_dp*eta*ucvol)
355 : :
356 : : !Finally assemble Ewald energy, eew
357 [ + + ]: 7686 : if(icutcoul.ne.3) then
358 : : !eew=sumg+sumr-chsq*reta/sqrt(pi)-fac
359 : 10 : eew=sumg+sumr-chsq*reta/sqrt(pi)
360 : : else
361 : 7676 : eew=sumg+sumr-chsq*reta/sqrt(pi)-fac
362 : : end if
363 : :
364 [ - + ]: 7686 : ABI_FREE(gcutoff)
365 : :
366 : : !DEBUG
367 : : !write(std_out,*)'eew=sumg+sumr-chsq*reta/sqrt(pi)-fac'
368 : : !write(std_out,*)eew,sumg,sumr,chsq*reta/sqrt(pi),fac
369 : : !ENDDEBUG
370 : :
371 : : !Length scale grads handled with stress tensor, ewald2
372 : :
373 : : !Output the final values of ng and nr
374 : : ! write(msg, '(a,a,i4,a,i4)' )ch10,' ewald : nr and ng are ',nr,' and ',ng
375 : : ! call wrtout(std_out,msg,'COLL')
376 : :
377 : 7686 : end subroutine ewald
378 : : !!***
379 : :
380 : : !----------------------------------------------------------------------
381 : :
382 : : !!****f* m_ewald/ewald2
383 : : !!
384 : : !! NAME
385 : : !! ewald2
386 : : !!
387 : : !! FUNCTION
388 : : !! Compute the part of the stress tensor coming from the Ewald energy
389 : : !! which is calculated by derivating the Ewald energy with respect to strain.
390 : : !! See Nielsen and Martin, Phys. Rev. B 32, 3792 (1985) [[cite:Nielsen1985a]].
391 : : !! Definition of stress tensor is $(1/ucvol)*d(Etot)/d(strain(a,b))$.
392 : : !!
393 : : !! INPUTS
394 : : !! gmet(3,3)=metric tensor in reciprocal space (bohr^-2)
395 : : !! natom=number of atoms in umit cell
396 : : !! ntypat=number of type of atoms
397 : : !! rmet(3,3)=metric tensor in real space (bohr^2) (inverse transpose of gmet)
398 : : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
399 : : !! typat(natom)=integer label of each type of atom (1,2,...)
400 : : !! ucvol=unit cell volume (bohr^3)
401 : : !! xred(3,natom)=relative coords of atoms in unit cell (dimensionless)
402 : : !! zion(ntypat)=charge on each type of atom (real number)
403 : : !!
404 : : !! OUTPUT
405 : : !! $stress(6)=(1/ucvol)*gradient$ of Ewald energy with respect to strain,
406 : : !! in hartrees/bohr^3
407 : : !! Cartesian components of stress are provided for this symmetric
408 : : !! tensor in the order 11 22 33 32 31 21.
409 : : !!
410 : : !! SOURCE
411 : :
412 : 6331 : subroutine ewald2(gmet,natom,ntypat,rmet,rprimd,stress,typat,ucvol,xred,zion)
413 : :
414 : : !Arguments ------------------------------------
415 : : !scalars
416 : : integer,intent(in) :: natom,ntypat
417 : : real(dp),intent(in) :: ucvol
418 : : !arrays
419 : : integer,intent(in) :: typat(natom)
420 : : real(dp),intent(in) :: gmet(3,3),rmet(3,3),rprimd(3,3),xred(3,natom)
421 : : real(dp),intent(in) :: zion(ntypat)
422 : : real(dp),intent(out) :: stress(6)
423 : :
424 : : !Local variables-------------------------------
425 : : !scalars
426 : : integer :: ia,ib,ig1,ig2,ig3,ir1,ir2,ir3,newg,newr,ng,nr
427 : : real(dp) :: arg1,arg2,arg3,ch,dderfc,derfc_arg,direct,eta,fac,fraca1
428 : : real(dp) :: fraca2,fraca3,fracb1,fracb2,fracb3,g1,g2,g3,gsq,r1,r1c,r2,r2c
429 : : real(dp) :: minexparg
430 : : real(dp) :: r3,r3c,recip,reta,rmagn,rsq,summi,summr,t1,t2,t3,t4,t5,t6,term1
431 : : real(dp) :: term2,term3,term4
432 : : !arrays
433 : : real(dp) :: gprimd(3,3),strg(6),strr(6)
434 : : ! *************************************************************************
435 : :
436 : : !Define dimensional reciprocal space primitive translations gprimd
437 : : !(inverse transpose of rprimd)
438 : 6331 : call matr3inv(rprimd,gprimd)
439 : :
440 : : !This is the minimum argument of an exponential, with some safety
441 : 6331 : minexparg=log(tiny(0._dp))+five
442 : :
443 : : !Add up total charge and sum of charge^2 in cell
444 : 6331 : ch=0._dp
445 [ + + ]: 25039 : do ia=1,natom
446 : 25039 : ch=ch+zion(typat(ia))
447 : : end do
448 : :
449 : : !Compute eta, the Ewald summation convergence parameter,
450 : : !for approximately optimized summations:
451 : : direct=rmet(1,1)+rmet(1,2)+rmet(1,3)+rmet(2,1)+&
452 : 6331 : & rmet(2,2)+rmet(2,3)+rmet(3,1)+rmet(3,2)+rmet(3,3)
453 : : recip=gmet(1,1)+gmet(1,2)+gmet(1,3)+gmet(2,1)+&
454 : 6331 : & gmet(2,2)+gmet(2,3)+gmet(3,1)+gmet(3,2)+gmet(3,3)
455 : : !Here, a bias is introduced, because G-space summation scales
456 : : !better than r space summation !
457 : 6331 : eta=pi*200.0_dp/33.0_dp*sqrt(1.69_dp*recip/direct)
458 : :
459 : 6331 : fac=pi**2/eta
460 : :
461 : : !Conduct reciprocal space summations
462 : 6331 : strg(1:6)=0.0_dp
463 : :
464 : : !Sum over G space, done shell after shell until all
465 : : !contributions are too small
466 : 6331 : ng=0
467 : : do
468 : 270962 : ng=ng+1
469 : 270962 : newg=0
470 : :
471 [ + + ]: 13836894 : do ig3=-ng,ng
472 [ + + ]: 1111516032 : do ig2=-ng,ng
473 [ + + ]: >17024*10^7 : do ig1=-ng,ng
474 : :
475 : : ! Exclude shells previously summed over
476 [ + + + + : >17023*10^7 : if(abs(ig1)==ng .or. abs(ig2)==ng .or. abs(ig3)==ng .or. ng==1 ) then
+ + + + ]
477 : :
478 : : ! Compute Cartesian components of each G
479 : : ! TODO : make this a blas call, and batch things up
480 : 6425457671 : g1=gprimd(1,1)*ig1+gprimd(1,2)*ig2+gprimd(1,3)*ig3
481 : 6425457671 : g2=gprimd(2,1)*ig1+gprimd(2,2)*ig2+gprimd(2,3)*ig3
482 : 6425457671 : g3=gprimd(3,1)*ig1+gprimd(3,2)*ig2+gprimd(3,3)*ig3
483 : : ! Compute |G|^2 (no pi factors)
484 : 6425457671 : gsq=(g1**2+g2**2+g3**2)
485 : :
486 : : ! skip g=0:
487 [ + + ]: 6425457671 : if (gsq>1.0d-20) then
488 : 6425451340 : arg1=fac*gsq
489 : :
490 : : ! larger arg1 gives 0 contribution because of exp(-arg1)
491 [ + + ]: 6425451340 : if (arg1<= -minexparg) then
492 : : ! When any term contributes then include next shell
493 : 1646238332 : newg=1
494 : 1646238332 : term1=exp(-arg1)/arg1
495 : 1646238332 : summr = 0.0_dp
496 : 1646238332 : summi = 0.0_dp
497 [ + + ]: 6756439430 : do ia=1,natom
498 : 5110201098 : arg2=two_pi*(ig1*xred(1,ia)+ig2*xred(2,ia)+ig3*xred(3,ia))
499 : : ! Sum real and imaginary parts (avoid complex variables)
500 : 5110201098 : summr=summr+zion(typat(ia))*cos(arg2)
501 : 6756439430 : summi=summi+zion(typat(ia))*sin(arg2)
502 : : end do
503 : :
504 : : ! Avoid underflow error messages
505 [ + + ]: 1646238332 : if (abs(summr)<1.d-16) summr=0.0_dp
506 [ + + ]: 1646238332 : if (abs(summi)<1.d-16) summi=0.0_dp
507 : :
508 : 1646238332 : term2=(2._dp/gsq)*(1._dp+arg1)
509 : 1646238332 : t1=term2*g1*g1-1._dp
510 : 1646238332 : t2=term2*g2*g2-1._dp
511 : 1646238332 : t3=term2*g3*g3-1._dp
512 : 1646238332 : t4=term2*g2*g3
513 : 1646238332 : t5=term2*g1*g3
514 : 1646238332 : t6=term2*g1*g2
515 : 1646238332 : term3=term1*(summr*summr+summi*summi)
516 : 1646238332 : strg(1)=strg(1)+t1*term3
517 : 1646238332 : strg(2)=strg(2)+t2*term3
518 : 1646238332 : strg(3)=strg(3)+t3*term3
519 : 1646238332 : strg(4)=strg(4)+t4*term3
520 : 1646238332 : strg(5)=strg(5)+t5*term3
521 : 1646238332 : strg(6)=strg(6)+t6*term3
522 : :
523 : : end if ! End condition not being larger than -minexparg
524 : : end if ! End skip g=0
525 : :
526 : : end if ! End triple loop and condition of new shell
527 : : end do
528 : : end do
529 : : end do
530 : :
531 : : ! Check if new shell must be calculated
532 [ + + ]: 270962 : if (newg==0) exit
533 : : end do ! End loop on new shell. Note that there is an "exit" instruction within the loop
534 : :
535 : :
536 : : !Conduct real space summations
537 : 6331 : reta=sqrt(eta)
538 : 6331 : strr(1:6)=0.0_dp
539 : :
540 : : !Loop on shells in r-space as was done in g-space
541 : 6331 : nr=0
542 : : do
543 : 18221 : nr=nr+1
544 : 18221 : newr=0
545 : :
546 [ + + ]: 110270 : do ir3=-nr,nr
547 [ + + ]: 638827 : do ir2=-nr,nr
548 [ + + ]: 4019471 : do ir1=-nr,nr
549 [ + + + + : 3927422 : if( abs(ir3)==nr .or. abs(ir2)==nr .or. abs(ir1)==nr .or. nr==1 )then
+ + + + ]
550 : :
551 [ + + ]: 9448833 : do ia=1,natom
552 : : ! Convert reduced atomic coordinates to [0,1)
553 : 7229980 : fraca1=xred(1,ia)-aint(xred(1,ia))+0.5_dp-sign(0.5_dp,xred(1,ia))
554 : 7229980 : fraca2=xred(2,ia)-aint(xred(2,ia))+0.5_dp-sign(0.5_dp,xred(2,ia))
555 : 7229980 : fraca3=xred(3,ia)-aint(xred(3,ia))+0.5_dp-sign(0.5_dp,xred(3,ia))
556 [ + + ]: 92869881 : do ib=1,natom
557 : 82241036 : fracb1=xred(1,ib)-aint(xred(1,ib))+0.5_dp-sign(0.5_dp,xred(1,ib))
558 : 82241036 : fracb2=xred(2,ib)-aint(xred(2,ib))+0.5_dp-sign(0.5_dp,xred(2,ib))
559 : 82241036 : fracb3=xred(3,ib)-aint(xred(3,ib))+0.5_dp-sign(0.5_dp,xred(3,ib))
560 : 82241036 : r1=ir1+fracb1-fraca1
561 : 82241036 : r2=ir2+fracb2-fraca2
562 : 82241036 : r3=ir3+fracb3-fraca3
563 : : ! Convert from reduced to cartesian coordinates
564 : 82241036 : r1c=rprimd(1,1)*r1+rprimd(1,2)*r2+rprimd(1,3)*r3
565 : 82241036 : r2c=rprimd(2,1)*r1+rprimd(2,2)*r2+rprimd(2,3)*r3
566 : 82241036 : r3c=rprimd(3,1)*r1+rprimd(3,2)*r2+rprimd(3,3)*r3
567 : : ! Compute |r|^2
568 : 82241036 : rsq=r1c**2+r2c**2+r3c**2
569 : 82241036 : rmagn=sqrt(rsq)
570 : :
571 : : ! Avoid zero denominators in 'term':
572 [ + + ]: 89471016 : if (rmagn>=1.0d-12) then
573 : :
574 : : ! Note: erfc(8) is about 1.1e-29, so do not bother with larger arg.
575 : : ! Also: exp(-64) is about 1.6e-28, so do not bother with larger arg**2 in exp.
576 : 82222328 : arg3=reta*rmagn
577 [ + + ]: 82222328 : if (arg3<8.0_dp) then
578 : 4512404 : newr=1
579 : : ! derfc computes the complementary error function
580 : : ! dderfc is the derivative of the complementary error function
581 : 4512404 : dderfc=(-2/sqrt(pi))*exp(-eta*rsq)
582 : 4512404 : derfc_arg = abi_derfc(arg3)
583 : 4512404 : term3=dderfc-derfc_arg/arg3
584 : 4512404 : term4=zion(typat(ia))*zion(typat(ib))*term3
585 : 4512404 : strr(1)=strr(1)+term4*r1c*r1c/rsq
586 : 4512404 : strr(2)=strr(2)+term4*r2c*r2c/rsq
587 : 4512404 : strr(3)=strr(3)+term4*r3c*r3c/rsq
588 : 4512404 : strr(4)=strr(4)+term4*r2c*r3c/rsq
589 : 4512404 : strr(5)=strr(5)+term4*r1c*r3c/rsq
590 : 4512404 : strr(6)=strr(6)+term4*r1c*r2c/rsq
591 : : end if ! End the condition of not being to large
592 : : end if ! End avoid zero denominator
593 : :
594 : : end do ! End loop over ib:
595 : : end do ! End loop over ia:
596 : :
597 : : end if ! End triple loop overs real space points, and associated new shell condition
598 : : end do
599 : : end do
600 : : end do
601 : :
602 : : ! Check if new shell must be calculated
603 [ + + ]: 18221 : if(newr==0) exit
604 : : end do ! End loop on new shells
605 : :
606 : : !Finally assemble stress tensor coming from Ewald energy, stress
607 : : !(note division by unit cell volume in accordance with definition
608 : : !found in Nielsen and Martin, Phys. Rev. B 32, 3792 (1985) [[cite:Nielsen1985a]]
609 : :
610 : 6331 : fac = pi/(2._dp*ucvol*eta)
611 : 6331 : stress(1)=(0.5_dp*reta*strr(1)+fac*(strg(1)+(ch**2)))/ucvol
612 : 6331 : stress(2)=(0.5_dp*reta*strr(2)+fac*(strg(2)+(ch**2)))/ucvol
613 : 6331 : stress(3)=(0.5_dp*reta*strr(3)+fac*(strg(3)+(ch**2)))/ucvol
614 : 6331 : stress(4)=(0.5_dp*reta*strr(4)+fac*strg(4))/ucvol
615 : 6331 : stress(5)=(0.5_dp*reta*strr(5)+fac*strg(5))/ucvol
616 : 6331 : stress(6)=(0.5_dp*reta*strr(6)+fac*strg(6))/ucvol
617 : :
618 : 6331 : end subroutine ewald2
619 : : !!***
620 : :
621 : : !!****f* m_ewald/ewald9
622 : : !! NAME
623 : : !! ewald9
624 : : !!
625 : : !! FUNCTION
626 : : !! Compute Ewald contribution to the dynamical matrix, at a given
627 : : !! q wavevector, including anisotropic dielectric tensor and effective charges
628 : : !! See Phys. Rev. B 55, 10355 (1997) [[cite:Gonze1997a]], equations (72) to (75).
629 : : !! This has been generalized to quadrupoles.
630 : : !! Delivers the left hand side of Eq.(72), possibly generalized.
631 : : !!
632 : : !! INPUTS
633 : : !! acell = lengths by which lattice vectors are multiplied
634 : : !! dielt(3,3)=dielectric tensor
635 : : !! gmet(3,3) = metric in reciprocal space.
636 : : !! gprim(3,3)=dimensionless primitive translations in reciprocal space
637 : : !! natom=number of atoms in unit cell
638 : : !! qphon(3)=phonon wavevector (same system of coordinates as the reciprocal lattice vectors)
639 : : !! rmet = metric in real space
640 : : !! rprim(3,3)=dimensionless primitive translations in real space
641 : : !! sumg0: if=1, the sum in reciprocal space must include g=0,
642 : : !! if=0, this contribution must be skipped (q=0 singularity)
643 : : !! ucvol=unit cell volume in (whatever length scale units)**3
644 : : !! xred(3,natom)=relative coords of atoms in unit cell (dimensionless)
645 : : !! zeff(3,3,natom)=effective charge on each atom, versus electric
646 : : !! field and atomic displacement
647 : : !! qdrp_cart(3,3,3,natom)=Quadrupole tensor on each atom in cartesian cordinates
648 : : !! option= 0: use old implementation;
649 : : !! 1: reduce the smalest argument of the exponentials to be evaluated,
650 : : !! set eta to 1 and skip real space sum, leads to a significant speedup
651 : : !! [dipquad] = if 1, atmfrc has been build without dipole-quadrupole part
652 : : !! [quadquad] = if 1, atmfrc has been build without quadrupole-quadrupole part
653 : : !!
654 : : !! OUTPUT
655 : : !! dyew(2,3,natom,3,natom)= Ewald part of the dynamical matrix,
656 : : !! second energy derivative wrt xred(3,natom) in Hartrees
657 : : !! Set to zero if all(zeff == zero)
658 : : !! eta: parameter used to split R and G-space summation
659 : : !!
660 : : !! NOTES
661 : : !! 1. The q=0 part should be subtracted, by another call to
662 : : !! the present routine, with q=0. The present routine correspond
663 : : !! to the quantity written A-bar in the explanatory notes.
664 : : !! If q=0 is asked, sumg0 should be put to 0. Otherwise, it should be put to 1.
665 : : !! 2. Because this routine can be used many times in the
666 : : !! evaluation of phonons in ppddb9, it has been
667 : : !! optimized carefully. There is still possibility
668 : : !! for improvement, by using bloking on G and R!
669 : : !! 3. There can be small numerical variations due to the
670 : : !! fact that the input dielectric tensor is usually
671 : : !! not perfectly symmetric.
672 : : !!
673 : : !! SOURCE
674 : :
675 : 175670 : subroutine ewald9(acell,dielt,dyew,gmet,gprim,natom,qphon,rmet,rprim,sumg0,ucvol,xred,zeff, qdrp_cart, eta, &
676 : : option, dipquad, quadquad) ! optional
677 : :
678 : : !Arguments -------------------------------
679 : : !scalars
680 : : integer,intent(in) :: natom,sumg0
681 : : integer,optional,intent(in) :: option, dipquad, quadquad
682 : : real(dp),intent(in) :: ucvol
683 : : !arrays
684 : : real(dp),intent(in) :: acell(3),dielt(3,3),gmet(3,3),gprim(3,3),qphon(3)
685 : : real(dp),intent(in) :: rmet(3,3),rprim(3,3),xred(3,natom),zeff(3,3,natom)
686 : : real(dp),intent(in) :: qdrp_cart(3,3,3,natom)
687 : : real(dp),intent(out) :: dyew(2,3,natom,3,natom), eta
688 : :
689 : : !Local variables -------------------------
690 : : !scalars
691 : : integer,parameter :: mr=10000
692 : : integer :: ia,ib,ig1,ig2,ig3,ii,ll,kk,ir,ir1,ir2,ir3,jj
693 : : integer :: info,lwork,mu,newg,newr,ng,nr,nu,ng_expxq
694 : : integer :: ewald_option
695 : : integer :: dipquad_,quadquad_
696 : : logical :: do_quadrupole
697 : : logical, save :: firstcall = .TRUE.
698 : : real(dp),parameter :: fac=4.0_dp/3.0_dp/sqrt(pi)
699 : : real(dp),parameter :: fact2=2.0_dp/sqrt(pi)
700 : : real(dp),parameter :: y2max=64.0_dp, y2min=1.0d-24
701 : : real(dp) :: cddi,cddr,cqdi,cqdr,cqqi,cqqr,g3,g4
702 : : real(dp) :: arg1,arg2,arg3,arga,c123r,c123i,c23i,c23r,detdlt,inv_detdlt
703 : : real(dp) :: direct,fact1,fact3,gsq,recip,reta,reta3,inv4eta
704 : : real(dp) :: minexparg,sigma_max
705 : : real(dp) :: term1,term2,term3,term4,term5,y2,yy,invy,invy2,derfc_yy
706 : : character(len=700) :: msg
707 : : !arrays
708 : : real(dp) :: c1i(2*mr+1),c1r(2*mr+1),c2i(2*mr+1),c2r(2*mr+1),c3i(2*mr+1)
709 : 351340 : real(dp) :: c3r(2*mr+1),cosqxred(natom),wdielt(3,3),eig_dielt(3),gpq(3),gpqfac(3,3),gpqgpq(3,3)
710 : 351340 : real(dp) :: invdlt(3,3),ircar(3),ircax(3),rr(3),sinqxred(natom)
711 : 351340 : real(dp) :: xredcar(3,natom),xredcax(3,natom),xredicar(3),xredicax(3),xx(3)
712 : : real(dp) :: gprimbyacell(3,3) !,tsec(2)
713 : 175670 : real(dp),allocatable :: dyddt(:,:,:,:,:), dydqt(:,:,:,:,:,:), dyqqt(:,:,:,:,:,:,:)
714 : 175670 : real(dp),allocatable :: work(:)
715 [ + + ]: 175670 : complex(dp) :: exp2piqx(natom)
716 : 175670 : complex(dp),allocatable :: expx1(:,:), expx2(:,:), expx3(:,:)
717 : : ! *********************************************************************
718 : :
719 : : ! This routine is expensive so skip the calculation and return zeros if zeff == zero.
720 : : ! Typically this happens when the DDB file does not contains zeff but dipdip = 1 is used (default).
721 [ + + + + : 16412068 : if (all(zeff == zero).and.all(qdrp_cart == zero)) then
+ + + + +
+ + + + +
+ + + + +
+ ]
722 [ + + + + : 3128282 : dyew = zero; return
+ + + + +
+ ]
723 : : end if
724 [ + + + + : 13098440 : do_quadrupole = any(qdrp_cart /= zero)
+ + + + +
+ ]
725 : :
726 : : ! Keep track of total time spent.
727 : : !call timab(1749, 1, tsec)
728 : :
729 : : ! Initialize dipquad and quadquad options
730 [ + + ]: 135760 : dipquad_=0; if(present(dipquad)) dipquad_=dipquad
731 [ + + ]: 135760 : quadquad_=0; if(present(quadquad)) quadquad_=quadquad
732 : :
733 : : ! Deactivate real space sums for quadrupolar fields or for dipdip = -1
734 [ + + ]: 135760 : ewald_option = 0; if (present(option)) ewald_option = option
735 [ + + + + ]: 135760 : if (do_quadrupole.and.(dipquad_==1.or.quadquad_==1)) ewald_option = 1
736 : : !ewald_option = 0
737 : :
738 : : !This is the minimum argument of an exponential, with some safety
739 : 135616 : minexparg=log(tiny(0._dp))+five
740 [ + + ]: 135760 : if (ewald_option == 1) minexparg=-20.0_dp
741 : :
742 : : ! initialize complex phase factors
743 [ + + ]: 473720 : do ia = 1, natom
744 : : arga = two_pi*( (qphon(1))*xred(1,ia)&
745 : : +(qphon(2))*xred(2,ia)&
746 : 337960 : +(qphon(3))*xred(3,ia) )
747 : 473720 : exp2piqx(ia) = exp(arga*j_dpc)
748 : : end do
749 : 135760 : ng_expxq = 1000
750 [ - + + - : 407280 : ABI_MALLOC(expx1, (-ng_expxq:ng_expxq, natom))
+ - + - -
+ - + ]
751 [ - + + - : 271520 : ABI_MALLOC(expx2, (-ng_expxq:ng_expxq, natom))
- + - + ]
752 [ - + + - : 271520 : ABI_MALLOC(expx3, (-ng_expxq:ng_expxq, natom))
- + - + ]
753 [ + + ]: 473720 : do ia = 1, natom
754 [ + + ]: 676731680 : do ig1 = -ng_expxq, ng_expxq
755 : 676257960 : expx1(ig1, ia) = exp(ig1*two_pi*xred(1,ia)*j_dpc)
756 : 676257960 : expx2(ig1, ia) = exp(ig1*two_pi*xred(2,ia)*j_dpc)
757 : 676595920 : expx3(ig1, ia) = exp(ig1*two_pi*xred(3,ia)*j_dpc)
758 : : end do
759 : : end do
760 : :
761 : 135760 : gprimbyacell = gprim
762 [ + + ]: 543040 : gprimbyacell(:,1) = gprimbyacell(:,1) / acell(1)
763 [ + + ]: 543040 : gprimbyacell(:,2) = gprimbyacell(:,2) / acell(2)
764 [ + + ]: 543040 : gprimbyacell(:,3) = gprimbyacell(:,3) / acell(3)
765 : :
766 : : ! compute eta for approximately optimized summations:
767 : : direct=rmet(1,1)+rmet(1,2)+rmet(1,3)+rmet(2,1)+&
768 : 135760 : rmet(2,2)+rmet(2,3)+rmet(3,1)+rmet(3,2)+rmet(3,3)
769 : : recip=gmet(1,1)+gmet(1,2)+gmet(1,3)+gmet(2,1)+&
770 : 135760 : gmet(2,2)+gmet(2,3)+gmet(3,1)+gmet(3,2)+gmet(3,3)
771 : :
772 : 135760 : eta=pi*100.0_dp/33.0_dp*sqrt(1.69_dp*recip/direct)
773 : :
774 : : ! Compute a material-dependent width for the Gaussians that hopefully
775 : : ! will make the Ewald real-space summation unnecessary.
776 [ + + ]: 135760 : if (ewald_option == 1) then
777 : :
778 : 145 : wdielt(:,:)=dielt(:,:)
779 : :
780 : : ! Diagonalize dielectric matrix
781 : 145 : lwork=-1
782 [ - + ]: 145 : ABI_MALLOC(work,(10))
783 : 145 : call dsyev('N','U',3, wdielt, 3, eig_dielt, work, lwork,info)
784 : 145 : lwork=nint(work(1))
785 : 145 : ABI_FREE(work)
786 : :
787 [ - + + - : 435 : ABI_MALLOC(work,(lwork))
+ - + - -
+ - + ]
788 : 145 : call dsyev('V','U',3, wdielt, 3, eig_dielt, work, lwork,info)
789 [ - + ]: 145 : ABI_FREE(work)
790 : :
791 : : ! This is a tentative maximum value for the gaussian width in real space
792 : 145 : sigma_max=three
793 : :
794 : : ! Set eta taking into account that the eps_inf is used as a metric in reciprocal space
795 [ - + + - : 725 : eta=sqrt(maxval(eig_dielt))/sigma_max
+ + ]
796 : :
797 [ + + ]: 145 : if (firstcall) then
798 : 3 : firstcall = .FALSE.
799 : 3 : write(msg, '(4a,f9.4,9a)' ) ch10,&
800 : 3 : ' Warning : due to the use of quadrupolar fields, the width of the reciprocal space gaussians', ch10, &
801 : 3 : ' in ewald9 has been set to eta= ', eta, ' 1/bohr and the real-space sums have been neglected.', ch10, &
802 : 3 : ' One should check whether this choice leads to correct results for the specific system under study', &
803 : 3 : ' and q-point grid.',ch10, &
804 : 3 : ' It is recommended to check that calculations with dipdip=1 and -1 (both with dipquad=0 and quadquad=0)', ch10, &
805 : 6 : ' lead to identical results. Otherwise increase the resolution of the q-point grid and repeat this test.', ch10
806 [ + + ]: 9 : call wrtout([ab_out,std_out], msg)
807 : : end if
808 : :
809 : : !Internally eta is the square of the gaussians width
810 : 145 : eta=eta*eta
811 : : end if
812 : :
813 : 135760 : inv4eta = one / four / eta
814 : :
815 [ - + + - : 678800 : ABI_MALLOC(dyddt,(2,3,natom,3,natom))
+ - - + +
- + - + -
- + - + ]
816 [ - + + - : 678800 : ABI_MALLOC(dydqt,(2,3,natom,3,natom,3))
- + + - +
- + - + -
- + - + ]
817 [ - + + - : 678800 : ABI_MALLOC(dyqqt,(2,3,natom,3,natom,3,3))
- + + - +
- + - + -
- + - + ]
818 : :
819 [ + + + + : 36459200 : dyddt = zero
+ + + + +
+ ]
820 [ + + + + : 109513360 : dydqt = zero
+ + + + +
+ + + ]
821 [ + + + + : 328675840 : dyqqt = zero
+ + + + +
+ + + +
+ ]
822 : :
823 : : !Sum terms over g space:
824 : : ng=0
825 : : do
826 : 1827227 : ng=ng+1
827 : :
828 : : ! if needed, update the complex phases for larger G vectors
829 [ - + ]: 1827227 : if (ng > ng_expxq) then
830 : : !write(std_out,*)"have to realloc"
831 [ # # ]: 0 : ABI_FREE(expx1)
832 [ # # ]: 0 : ABI_FREE(expx2)
833 [ # # ]: 0 : ABI_FREE(expx3)
834 : :
835 : 0 : ng_expxq = ng_expxq*2
836 : : ! TODO: half of this space is not needed, as it contains the complex conjugate of the other half.
837 : : ! present duplication avoids if statements inside the loop, however
838 [ # # # # : 0 : ABI_MALLOC(expx1, (-ng_expxq:ng_expxq, natom))
# # # # #
# # # # #
# # ]
839 [ # # # # : 0 : ABI_MALLOC(expx2, (-ng_expxq:ng_expxq, natom))
# # # # #
# # # ]
840 [ # # # # : 0 : ABI_MALLOC(expx3, (-ng_expxq:ng_expxq, natom))
# # # # #
# # # ]
841 [ # # ]: 0 : do ia = 1, natom
842 [ # # ]: 0 : do ig1 = -ng_expxq, ng_expxq
843 : 0 : expx1(ig1, ia) = exp(ig1*two_pi*xred(1,ia)*j_dpc)
844 : 0 : expx2(ig1, ia) = exp(ig1*two_pi*xred(2,ia)*j_dpc)
845 : 0 : expx3(ig1, ia) = exp(ig1*two_pi*xred(3,ia)*j_dpc)
846 : : end do
847 : : end do
848 : : end if
849 : :
850 : 1827227 : newg=0
851 [ + + ]: 32786322 : do ig3=-ng,ng
852 [ + + ]: 777120741 : do ig2=-ng,ng
853 [ + + ]: 26476385801 : do ig1=-ng,ng
854 [ + + + + : 26445426706 : if(abs(ig1)==ng .or. abs(ig2)==ng .or. abs(ig3)==ng .or. ng==1 )then
+ + + + ]
855 : :
856 : : gpq(1)=(ig1+qphon(1))*gprimbyacell(1,1)+(ig2+qphon(2))*&
857 : 4109250950 : gprimbyacell(1,2)+(ig3+qphon(3))*gprimbyacell(1,3)
858 : : gpq(2)=(ig1+qphon(1))*gprimbyacell(2,1)+(ig2+qphon(2))*&
859 : 4109250950 : gprimbyacell(2,2)+(ig3+qphon(3))*gprimbyacell(2,3)
860 : : gpq(3)=(ig1+qphon(1))*gprimbyacell(3,1)+(ig2+qphon(2))*&
861 : 4109250950 : gprimbyacell(3,2)+(ig3+qphon(3))*gprimbyacell(3,3)
862 : 4109250950 : gsq=zero
863 [ + + ]: 16437003800 : do jj=1,3
864 [ + + ]: 53420262350 : do ii=1,3
865 : 36983258550 : gpqgpq(ii,jj)=gpq(ii)*gpq(jj)
866 : 49311011400 : gsq=gsq+gpqgpq(ii,jj)*dielt(ii,jj)
867 : : end do
868 : : end do
869 : :
870 : : ! Skip q=0:
871 [ + + ]: 4109250950 : if (gsq<1.0d-20) then
872 [ - + ]: 19491 : if (sumg0==1) then
873 : : write(msg,'(5a)' )&
874 : 0 : 'The phonon wavelength should not be zero :',ch10,&
875 : 0 : 'there are non-analytical terms that cannot be treated.',ch10,&
876 : 0 : 'Action: subtract this wavelength from the input file.'
877 : 0 : ABI_ERROR(msg)
878 : : end if
879 : :
880 : : else
881 : :
882 : 4109231459 : arg1=(two_pi**2)*gsq* inv4eta
883 : :
884 : : ! Larger arg gives 0 contribution:
885 [ + + ]: 4109231459 : if (arg1<= -minexparg ) then
886 : 1073049269 : newg=1
887 : :
888 : : ! Here calculate the term
889 : 1073049269 : term1=exp(-arg1)/gsq
890 [ + + ]: 4292197076 : do jj=1,3
891 [ + + ]: 13949640497 : do ii=1,3
892 : 12876591228 : gpqfac(ii,jj)=gpqgpq(ii,jj)*term1
893 : : end do
894 : : end do
895 : :
896 : : ! MJV: replaced old calls to cos and sin.
897 : : ! Checked for 10 tests in v2 that max error is about 6.e-15, usually < 2.e-15
898 [ + + ]: 4452116262 : do ia=1,natom
899 : 3379066993 : cosqxred(ia)= real(exp2piqx(ia)*expx1(ig1, ia)*expx2(ig2, ia)*expx3(ig3, ia))
900 : 4452116262 : sinqxred(ia)=aimag(exp2piqx(ia)*expx1(ig1, ia)*expx2(ig2, ia)*expx3(ig3, ia))
901 : : end do
902 : :
903 : : ! First, the diagonal terms
904 [ + + ]: 4292197076 : do nu=1,3
905 [ + + ]: 14429398055 : do ia=1,natom
906 [ + + ]: 33630750744 : do mu=nu,3
907 : 30411602937 : dyddt(1,mu,ia,nu,ia)=dyddt(1,mu,ia,nu,ia)+gpqfac(mu,nu)
908 : : end do
909 : : end do
910 : : end do
911 : :
912 : : ! Then, the non-diagonal ones
913 [ + + ]: 3379066993 : do ib=2,natom
914 [ + + ]: 9784825564 : do ia=1,ib-1
915 : : ! phase factor dipole-dipole
916 : 6405758571 : cddr=cosqxred(ia)*cosqxred(ib)+sinqxred(ia)*sinqxred(ib)
917 : 6405758571 : cddi=sinqxred(ia)*cosqxred(ib)-cosqxred(ia)*sinqxred(ib)
918 : :
919 : : ! Dipole-dipole contribution
920 [ + + ]: 27929052008 : do nu=1,3
921 [ + + ]: 64057585710 : do mu=nu,3
922 : 38434551426 : dyddt(1,mu,ia,nu,ib)=dyddt(1,mu,ia,nu,ib)+gpqfac(mu,nu)*cddr
923 : 57651827139 : dyddt(2,mu,ia,nu,ib)=dyddt(2,mu,ia,nu,ib)+gpqfac(mu,nu)*cddi
924 : : end do
925 : : end do
926 : : end do
927 : : end do
928 : :
929 [ + + ]: 1073049269 : if (do_quadrupole) then
930 [ + + ]: 273367160 : do ib=1,natom
931 [ + + ]: 647448417 : do ia=1,natom
932 : :
933 : : ! phase factor for dipole-quadrupole
934 : 374081257 : cqdr=cosqxred(ia)*sinqxred(ib)-sinqxred(ia)*cosqxred(ib)
935 : 374081257 : cqdi=cosqxred(ia)*cosqxred(ib)+sinqxred(ia)*sinqxred(ib)
936 : :
937 : : ! phase factor quadrupole-quadrupole
938 : 374081257 : cqqr=cosqxred(ia)*cosqxred(ib)+sinqxred(ia)*sinqxred(ib)
939 : 374081257 : cqqi=sinqxred(ia)*cosqxred(ib)-cosqxred(ia)*sinqxred(ib)
940 : :
941 : : ! Dipole-quadrupole contribution
942 [ + + ]: 1496325028 : do ii=1,3
943 [ + + ]: 4863056341 : do jj=1,3
944 [ + + ]: 14589169023 : do kk=1,3
945 : 10100193939 : g3=gpq(ii)*gpq(jj)*gpq(kk)
946 : 10100193939 : dydqt(1,ii,ia,jj,ib,kk)=dydqt(1,ii,ia,jj,ib,kk)+g3*term1*cqdr
947 : 13466925252 : dydqt(2,ii,ia,jj,ib,kk)=dydqt(2,ii,ia,jj,ib,kk)+g3*term1*cqdi
948 : : end do ! kk
949 : : end do ! jj
950 : : end do ! ii
951 : :
952 : : ! Quadrupole-quadrupole contribution
953 [ + + ]: 1679441775 : do ii=1,3
954 [ + + ]: 4863056341 : do jj=1,3
955 [ + + ]: 14589169023 : do kk=1,3
956 [ + + ]: 43767507069 : do ll=1,3
957 : 30300581817 : g4 = gpq(ii)*gpq(jj)*gpq(kk)*gpq(ll)
958 : 30300581817 : dyqqt(1,ii,ia,jj,ib,kk,ll)=dyqqt(1,ii,ia,jj,ib,kk,ll)+g4*term1*cqqr
959 : 40400775756 : dyqqt(2,ii,ia,jj,ib,kk,ll)=dyqqt(2,ii,ia,jj,ib,kk,ll)+g4*term1*cqqi
960 : : end do
961 : : end do ! kk
962 : : end do ! jj
963 : : end do ! ii
964 : : end do ! ia
965 : : end do ! ib
966 : : end if
967 : :
968 : : end if ! endif exp() argument is smaller than -minexparg
969 : : end if ! Endif g/=0 :
970 : : end if ! End triple summation over Gs:
971 : : end do
972 : : end do
973 : : end do
974 : :
975 : : ! Check if new shell must be calculated
976 [ + + ]: 1827227 : if(newg==0)exit
977 : : end do
978 : :
979 : : ! Multiplies by common factor
980 : 135760 : fact1=4.0_dp*pi/ucvol
981 [ + + ]: 473720 : do ib=1,natom
982 [ + + ]: 1225560 : do ia=1,ib
983 [ + + ]: 3345320 : do nu=1,3
984 [ + + ]: 7518400 : do mu=nu,3
985 : 4511040 : dyddt(1,mu,ia,nu,ib)=dyddt(1,mu,ia,nu,ib)*fact1
986 : 6766560 : dyddt(2,mu,ia,nu,ib)=dyddt(2,mu,ia,nu,ib)*fact1
987 : : end do
988 : : end do
989 : : end do
990 : : end do
991 [ + + ]: 135760 : if (do_quadrupole) then
992 [ + + + + : 2998350 : dydqt=dydqt*fact1/two * two_pi
+ + + + +
+ + + ]
993 [ + + + + : 9001431 : dyqqt=dyqqt*fact1/four * two_pi ** 2
+ + + + +
+ + + +
+ ]
994 : : end if
995 : :
996 : 135760 : reta=sqrt(eta)
997 : 135760 : reta3=-eta*reta
998 : :
999 : : !Calculating the inverse (transpose) of the dielectric tensor
1000 : 135760 : call matr3inv(dielt,invdlt)
1001 : :
1002 : : !Calculating the determinant of the dielectric tensor
1003 : : detdlt=dielt(1,1)*dielt(2,2)*dielt(3,3)+dielt(1,3)*dielt(2,1)*&
1004 : : dielt(3,2)+dielt(1,2)*dielt(2,3)*dielt(3,1)-dielt(1,3)*&
1005 : : dielt(2,2)*dielt(3,1)-dielt(1,1)*dielt(2,3)*dielt(3,2)-&
1006 : 135760 : dielt(1,2)*dielt(2,1)*dielt(3,3)
1007 : :
1008 [ - + ]: 135760 : if(detdlt<tol6)then
1009 : : write(msg, '(a,es16.6,11a)' )&
1010 : 0 : 'The determinant of the dielectrix matrix, detdlt=',detdlt,' is smaller than 1.0d-6.',ch10,&
1011 : 0 : 'The use of the dipole-dipole model for interatomic force constants is not possible.',ch10,&
1012 : 0 : 'It is likely that you have not treated the electric field perturbations,',ch10,&
1013 : 0 : 'because you not are dealing with an insulator, so that',ch10,&
1014 : 0 : 'your dielectric matrix was simply set to zero in the Derivative DataBase.',ch10,&
1015 : 0 : 'Action: set the input variable dipdip to 0 .'
1016 : 0 : ABI_ERROR(msg)
1017 : : end if
1018 : :
1019 : 135760 : inv_detdlt = one / sqrt(detdlt)
1020 : 135760 : fact3=reta3 * inv_detdlt
1021 : :
1022 [ + + ]: 135760 : if (ewald_option /= 1) then
1023 : : ! Preparing the loop on real space
1024 [ + + ]: 473285 : do ia=1,natom
1025 [ + + ]: 1486295 : do ii=1,3
1026 : : xredcar(ii,ia)=(xred(1,ia)*acell(1)*rprim(ii,1)+&
1027 : : xred(2,ia)*acell(2)*rprim(ii,2)+&
1028 : 1350680 : xred(3,ia)*acell(3)*rprim(ii,3) )*reta
1029 : : end do
1030 : : end do
1031 [ + + ]: 473285 : do ia=1,natom
1032 [ + + ]: 1486295 : do ii=1,3
1033 : : xredcax(ii,ia)= invdlt(1,ii)*xredcar(ii,ia)+&
1034 : : invdlt(2,ii)*xredcar(ii,ia)+&
1035 : 1350680 : invdlt(3,ii)*xredcar(ii,ia)
1036 : : end do
1037 : : end do
1038 : :
1039 : : ! Prepare the evaluation of exp(iq*R)
1040 [ + + ]: 2712571230 : do ir=-mr,mr
1041 : 2712435615 : arg1=-two_pi*qphon(1)*ir
1042 : 2712435615 : arg2=-two_pi*qphon(2)*ir
1043 : 2712435615 : arg3=-two_pi*qphon(3)*ir
1044 : 2712435615 : c1r(ir+mr+1)=cos(arg1)
1045 : 2712435615 : c1i(ir+mr+1)=sin(arg1)
1046 : 2712435615 : c2r(ir+mr+1)=cos(arg2)
1047 : 2712435615 : c2i(ir+mr+1)=sin(arg2)
1048 : 2712435615 : c3r(ir+mr+1)=cos(arg3)
1049 : 2712571230 : c3i(ir+mr+1)=sin(arg3)
1050 : : end do
1051 : :
1052 [ + - ]: 1173289 : do nr=1,mr
1053 : 1173289 : newr=0
1054 : :
1055 : : ! Begin big loop on real space vectors
1056 [ + + ]: 15000954 : do ir3=-nr,nr
1057 [ + + ]: 235663251 : do ir2=-nr,nr
1058 : :
1059 : : ! Here, construct the cosine and sine of q*R for components 2 and 3
1060 : 220662297 : c23r = c2r(ir2+mr+1) * c3r(ir3+mr+1) - c2i(ir2+mr+1) * c3i(ir3+mr+1)
1061 : 220662297 : c23i = c2i(ir2+mr+1) * c3r(ir3+mr+1) + c2r(ir2+mr+1) * c3i(ir3+mr+1)
1062 : :
1063 : : ! Also multiplies by fact3, because it is a rather economical place to do so
1064 : 220662297 : c23r=c23r * fact3
1065 : 220662297 : c23i=c23i * fact3
1066 : :
1067 [ + + ]: 4575695659 : do ir1=-nr,nr
1068 [ + + + + : 4561867994 : if( abs(ir3)==nr .or. abs(ir2)==nr .or. abs(ir1)==nr .or. nr==1 )then
+ + + + ]
1069 : :
1070 : : ! This is the real part and imaginary part of the phase factor exp(iq*R)
1071 : 1167563729 : c123r = c1r(ir1+mr+1) * c23r - c1i(ir1+mr+1) * c23i
1072 : 1167563729 : c123i = c1i(ir1+mr+1) * c23r + c1r(ir1+mr+1) * c23i
1073 : :
1074 [ + + ]: 4670254916 : do ii=1,3
1075 : : ircar(ii)= ( ir1*acell(1)*rprim(ii,1)+&
1076 : : ir2*acell(2)*rprim(ii,2)+&
1077 : 4670254916 : ir3*acell(3)*rprim(ii,3) ) * reta
1078 : : end do
1079 [ + + ]: 4670254916 : do ii=1,3
1080 : : ircax(ii)= invdlt(1,ii)*ircar(ii)+&
1081 : : invdlt(2,ii)*ircar(ii)+&
1082 : 4670254916 : invdlt(3,ii)*ircar(ii)
1083 : : end do
1084 : :
1085 : : ! Here loops on atoms
1086 [ + + ]: 4002401223 : do ib=1,natom
1087 [ + + ]: 11339349976 : do ii=1,3
1088 : 8504512482 : xredicar(ii)=ircar(ii)-xredcar(ii,ib)
1089 : 11339349976 : xredicax(ii)=ircax(ii)-xredcax(ii,ib)
1090 : : end do
1091 [ + + ]: 12704604322 : do ia=1,ib
1092 [ + + ]: 22114244524 : do ii=1,3
1093 : 16585683393 : rr(ii)=xredicar(ii)+xredcar(ii,ia)
1094 : 22114244524 : xx(ii)=xredicax(ii)+xredcax(ii,ia)
1095 : : end do
1096 : :
1097 : 5528561131 : y2=rr(1)*xx(1)+rr(2)*xx(2)+rr(3)*xx(3)
1098 : :
1099 : : ! The atoms should not be too far of each other
1100 [ + + ]: 8363398625 : if (y2 < y2max) then
1101 : : ! Note: erfc(8) is about 1.1e-29, so dont bother with larger y.
1102 : : ! Also: exp(-64) is about 1.6e-28, do dont bother with larger y**2 in exp.
1103 : :
1104 : : ! Avoid zero denominators in term:
1105 [ + + ]: 1271968614 : if (y2 >= y2min) then
1106 : 1271630944 : newr=1
1107 : 1271630944 : yy=sqrt(y2)
1108 : 1271630944 : invy=1.0_dp/yy
1109 : 1271630944 : invy2=invy**2
1110 : 1271630944 : derfc_yy = abi_derfc(yy)
1111 : 1271630944 : term2=derfc_yy*invy*invy2
1112 : 1271630944 : term3=fact2*exp(-y2)*invy2
1113 : 1271630944 : term4=-(term2+term3)
1114 : 1271630944 : term5=(3.0_dp*term2+term3*(3.0_dp+2.0_dp*y2))*invy2
1115 [ + + ]: 5086523776 : do nu=1,3
1116 [ + + ]: 12716309440 : do mu=nu,3
1117 : 7629785664 : dyddt(1,mu,ia,nu,ib)=dyddt(1,mu,ia,nu,ib)+c123r*(xx(nu)*xx(mu)*term5+term4*invdlt(nu,mu))
1118 : 11444678496 : dyddt(2,mu,ia,nu,ib)=dyddt(2,mu,ia,nu,ib)+c123i*(xx(nu)*xx(mu)*term5+term4*invdlt(nu,mu))
1119 : : end do
1120 : : end do
1121 : : else
1122 : : ! If zero denominator, the atoms should be identical
1123 [ - + ]: 337670 : if (ia/=ib)then
1124 : : write(msg, '(5a,i0,a,i0,a)' )&
1125 : 0 : 'The distance between two atoms seem to vanish.',ch10,&
1126 : 0 : 'This is not allowed.',ch10,&
1127 : 0 : 'Action: check the input for the atoms number',ia,' and',ib,'.'
1128 : 0 : ABI_ERROR(msg)
1129 : : else
1130 : : ! This is the correction when the atoms are identical
1131 [ + + ]: 1350680 : do nu=1,3
1132 [ + + ]: 4389710 : do mu=1,3
1133 : : dyddt(1,mu,ia,nu,ib)=dyddt(1,mu,ia,nu,ib)+&
1134 : 4052040 : fac*reta3*invdlt(nu,mu) * inv_detdlt
1135 : : end do
1136 : : end do
1137 : : end if
1138 : : end if ! End the condition for avoiding zero denominators
1139 : : end if ! End the condition of too large distance between atoms
1140 : : end do
1141 : : end do ! End loop over ia and ib :
1142 : : end if ! End triple loop over real space points:
1143 : : end do ! ir1
1144 : : end do ! ir2
1145 : : end do ! ir3
1146 : :
1147 : : ! Check if new shell must be calculated
1148 [ + + ]: 1173289 : if(newr==0)exit
1149 [ - + ]: 1173289 : if(newr==1 .and. nr==mr) ABI_BUG('mr is too small')
1150 : : end do
1151 : : end if ! check if should compute real part
1152 : :
1153 : : ! Now, symmetrizes
1154 [ + + ]: 337960 : do ib=1,natom-1
1155 [ + + ]: 944560 : do nu=1,3
1156 [ + + ]: 2050440 : do ia=ib+1,natom
1157 [ + + ]: 4331520 : do mu=nu,3
1158 : 2483280 : dyddt(1,mu,ia,nu,ib)= dyddt(1,mu,ib,nu,ia)
1159 : 3724920 : dyddt(2,mu,ia,nu,ib)=-dyddt(2,mu,ib,nu,ia)
1160 : : end do
1161 : : end do
1162 : : end do
1163 : : end do
1164 : :
1165 [ + + ]: 473720 : do ib=1,natom
1166 [ + + ]: 1149640 : do nu=2,3
1167 [ + + ]: 3345320 : do ia=1,natom
1168 [ + + ]: 6504520 : do mu=1,nu-1
1169 : 3497160 : dyddt(1,mu,ia,nu,ib)=dyddt(1,nu,ia,mu,ib)
1170 : 5828600 : dyddt(2,mu,ia,nu,ib)=dyddt(2,nu,ia,mu,ib)
1171 : : end do
1172 : : end do
1173 : : end do
1174 : : end do
1175 : :
1176 : : !Tests
1177 : : !write(std_out,*)' ewald9 : take into account the effective charges '
1178 [ + + + + : 36459200 : dyew = zero
+ + + + +
+ ]
1179 [ + + ]: 473720 : do ib=1,natom
1180 [ + + ]: 1487600 : do nu=1,3
1181 [ + + ]: 4849000 : do ia=1,natom
1182 [ + + ]: 15002520 : do mu=1,3
1183 [ + + ]: 45463080 : do ii=1,3
1184 [ + + ]: 136389240 : do jj=1,3
1185 : : ! dipole-dipole correction
1186 : : dyew(1,mu,ia,nu,ib)=dyew(1,mu,ia,nu,ib) + &
1187 : 94423320 : zeff(ii,mu,ia)*zeff(jj,nu,ib)*dyddt(1,ii,ia,jj,ib)
1188 : : dyew(2,mu,ia,nu,ib)=dyew(2,mu,ia,nu,ib) + &
1189 : 94423320 : zeff(ii,mu,ia)*zeff(jj,nu,ib)*dyddt(2,ii,ia,jj,ib)
1190 [ + + ]: 125897760 : if (do_quadrupole) then
1191 [ + + ]: 10101996 : do kk=1,3
1192 [ + + ]: 7576497 : if (dipquad_==1) then
1193 : : ! dipole-quadrupole correction
1194 : : dyew(1,mu,ia,nu,ib)=dyew(1,mu,ia,nu,ib) + &
1195 : : (zeff(ii,nu,ib)*qdrp_cart(kk,jj,mu,ia) - &
1196 : 139968 : zeff(ii,mu,ia)*qdrp_cart(kk,jj,nu,ib)) * dydqt(1,ii,ia,jj,ib,kk)
1197 : : dyew(2,mu,ia,nu,ib)=dyew(2,mu,ia,nu,ib) + &
1198 : : (zeff(ii,nu,ib)*qdrp_cart(kk,jj,mu,ia) - &
1199 : 139968 : zeff(ii,mu,ia)*qdrp_cart(kk,jj,nu,ib)) * dydqt(2,ii,ia,jj,ib,kk)
1200 : : end if
1201 : :
1202 : : ! quadrupole-quadrupole correction
1203 [ + + ]: 10101996 : if (quadquad_==1) then
1204 [ + + ]: 559872 : do ll=1,3
1205 : : dyew(1,mu,ia,nu,ib)=dyew(1,mu,ia,nu,ib) + &
1206 : 419904 : (qdrp_cart(ll,ii,mu,ia)*qdrp_cart(kk,jj,nu,ib)) * dyqqt(1,ii,ia,jj,ib,kk,ll)
1207 : : dyew(2,mu,ia,nu,ib)=dyew(2,mu,ia,nu,ib) + &
1208 : 559872 : (qdrp_cart(ll,ii,mu,ia)*qdrp_cart(kk,jj,nu,ib)) * dyqqt(2,ii,ia,jj,ib,kk,ll)
1209 : : end do
1210 : : end if
1211 : : end do
1212 : : end if
1213 : :
1214 : : end do
1215 : : end do
1216 : : end do
1217 : : end do
1218 : : end do
1219 : : end do
1220 : :
1221 [ - + ]: 135760 : ABI_FREE(expx1)
1222 [ - + ]: 135760 : ABI_FREE(expx2)
1223 [ - + ]: 135760 : ABI_FREE(expx3)
1224 [ - + ]: 135760 : ABI_FREE(dyddt)
1225 [ - + ]: 135760 : ABI_FREE(dydqt)
1226 [ - + ]: 135760 : ABI_FREE(dyqqt)
1227 : : !call timab(1749, 2, tsec)
1228 : :
1229 : : end subroutine ewald9
1230 : : !!***
1231 : :
1232 : : !!****f* m_ewald/ewald9_2D
1233 : : !!
1234 : : !! NAME
1235 : : !! ewald9_2D
1236 : : !!
1237 : : !! FUNCTION
1238 : : !! Compute the long-range electrostatics contribution to interatomic force constants
1239 : : !! in the bi-dimensional (2D) case, considering the 2D is embedded in a dielectric environment
1240 : : !! and has a given dielectric thickness. The singularity of the Coulomb potential is treated
1241 : : !! using the Ewald summation approach. It is possible to input a more complicate model
1242 : : !! with two consecutive dielectric slabs.
1243 : : !!
1244 : : !! INPUTS
1245 : : !! natom=number of atoms in unit cell
1246 : : !! acell(3)=length of unit cell vectors
1247 : : !! xred(3,natom)=reduced coordinates of the atoms
1248 : : !! rprim(3,3)=unit cell vectors (unscaled)
1249 : : !! dielt(3,3)=dielectric tensor of the 2D
1250 : : !! dyew(2,3,natom,3,natom)=long-range electrostatics IFCs following Ewald
1251 : : !! qphon(3)=phonon wavevector in reduced coordinates
1252 : : !! zeff(3,3,natom)=Born effective charge tensor
1253 : : !! qdrp_cart(3,3,3,natom)=Dynamical quadrupoles
1254 : : !! dielt_env=dielectric constant of the embedding environment (1 in vacuum)
1255 : : !! thick(2)=dielectric thicknesses of the slab, first value correspond to the outer
1256 : : !! dielectric, second to the inner dielectric slab (if any)
1257 : : !! sys_dim= system dimensionality (indicates axis without periodicity)
1258 : : !!
1259 : : !! OUTPUT
1260 : : !! dyew(2,3,natom,3,natom)=long-range electrostatics IFCs following Ewald
1261 : : !!
1262 : : !! SOURCE
1263 : :
1264 : 561 : subroutine ewald9_2D(natom,acell,xred,rprim,dielt,dyew,qphon,zeff,qdrp_cart,dielt_env,thick,sys_dim)
1265 : :
1266 : : !Arguments -------------------------------
1267 : : !scalars
1268 : : real(dp), intent(in) :: dielt_env
1269 : : integer :: natom, sys_dim
1270 : : !arrays
1271 : : real(dp),intent(in) :: acell(3),thick(2),xred(3,natom),dielt(3,3),qphon(3)
1272 : : real(dp),intent(in) :: rprim(3,3),zeff(3,3,natom),qdrp_cart(3,3,3,natom)
1273 : : real(dp),intent(out) :: dyew(2,3,natom,3,natom)
1274 : : character(len=700) :: msg
1275 : :
1276 : : !Local variables -------------------------
1277 : : !scalars
1278 : : integer :: gmax,idir1,idir2,ibz1,ibz2,ipert1,ipert2,inner_thick,ndir,mdir
1279 : : real(dp) :: detdlt, delta_perp, lambda, dielt_perp,dielt_perp1,dielt_perp2, eta,eta1,xi, dielt_eff
1280 : : real(dp) :: dielt_eff1,dielt_eff2, norm_kvec, norm_kvec0, phi
1281 : : real(dp) :: ewald_fun, ewald_fun1, ewald_fun2, rflct_coeff, out_thick
1282 : : real(dp) :: rflct_coeff1, rflct_coeff2, rprimd_perp
1283 : : real(dp) :: fac_erfc, fac_ewald1, fac_ewald2, fac_ewald2b,fac_exp
1284 : : real(dp) :: fac_exp1, trans_fun, fac_gauss, fac_mirror, fac_mirror1, fac_real
1285 : : real(dp) :: mean2_perp, mirror_diff, mirror_parapara, mirror_paraperp, mirror_perpperp
1286 : : real(dp) :: rvec_norm, sqrt_norm, ucsurf, xmean
1287 : : logical, save :: firstcall = .TRUE.
1288 : : !arrays
1289 : : integer :: periodic_dir(3)
1290 : 1122 : real(dp) :: dyew_real(2,3,natom,3,natom),dyew_rec(2,3,natom,3,natom),kvec(2),invdlt_para(2,2)
1291 : : real(dp) :: rprimd_para(2,2), gprimd_para(2,2), gvec(2)
1292 : : real(dp) :: norm_dielt(2), dielt_para(2,2), invdlt(3,3), qvec(3), qvec_para(2)
1293 : 1122 : real(dp) :: diff_xcart(3), xcart_para(2,natom),xcart_perp(natom)
1294 : 1122 : real(dp) :: kvec_para(2),zeff_para(2,3,natom), zeff_perp(3,natom)
1295 : 1122 : real(dp) :: xcart(3,natom), rprimd(3,3), gprimd(3,3), rvec_dielt(3)
1296 : 1122 : real(dp) :: qdrp_parapara(2,2,3,natom),qdrp_perpperp(3,natom), qdrp_paraperp(2,3,natom)
1297 : : real(dp) :: rho_gerade1(2),rho_gerade2(2), rho_ungerade1(2), rho_ungerade2(2)
1298 : : ! *************************************************************************
1299 : :
1300 : 561 : periodic_dir(:) = 0
1301 [ - + ]: 561 : if (sys_dim ==2) then ! 2D along x
1302 : 0 : periodic_dir(2) =1 ; periodic_dir(3) = 1
1303 [ - + ]: 561 : elseif (sys_dim==3) then ! 2D along y
1304 : 0 : periodic_dir(1) = 1 ; periodic_dir(3) = 1
1305 [ + - ]: 561 : elseif (sys_dim==4) then ! 2D along z
1306 : 561 : periodic_dir(1) = 1 ; periodic_dir(2) = 1
1307 : : end if
1308 : :
1309 : 561 : rprimd=zero
1310 [ + + ]: 2244 : do idir1=1,3
1311 [ + + ]: 7293 : do idir2=1,3
1312 : 6732 : rprimd(idir1,idir2) = rprim(idir1,idir2)*acell(idir2)
1313 : : end do
1314 : : end do
1315 : :
1316 : :
1317 [ + + + + ]: 8017 : xcart=zero
1318 [ + + ]: 2425 : do ipert1=1,natom
1319 [ + + + + : 30385 : xcart(:,ipert1)=matmul(rprimd, xred(:,ipert1))
+ + ]
1320 : : end do
1321 : :
1322 : 561 : call matr3inv(rprimd,gprimd)
1323 [ + + + + : 26657 : zeff_para=zero ; zeff_perp=zero
+ + + + +
+ ]
1324 [ + + ]: 2244 : qvec(:) = qphon(1)*gprimd(:,1)+qphon(2)*gprimd(:,2)
1325 : : ndir=0
1326 [ + + ]: 2244 : do idir1=1,3
1327 [ + + ]: 2244 : if (periodic_dir(idir1)==1) then
1328 : 1122 : ndir=ndir+1
1329 : 1122 : qvec_para(ndir)=qvec(idir1)
1330 [ + + + + ]: 16034 : zeff_para(ndir,:,:) = zeff(idir1,:,:)
1331 [ + + ]: 4850 : xcart_para(ndir,:) = xcart(idir1,:)
1332 : : else
1333 [ + + + + ]: 8017 : zeff_perp(:,:) = zeff(idir1,:,:)
1334 [ + + ]: 2425 : xcart_perp(:) = xcart(idir1,:)
1335 [ + - ]: 561 : if (qvec(idir1)>tol6) then !Check if phonon mode is not out-of-plane
1336 : : write(msg, '(a,es16.6,5a)')&
1337 : 0 : 'The phonon wavevector along the confined direction is',qvec(idir1),' 1/Bohr >1.0d-6',ch10,&
1338 : 0 : 'The phonon wavevector should be purely along the periodic direction', ch10, &
1339 : 0 : 'when using Ewald summation in 2D. Please check your input file and structure'
1340 : 0 : ABI_ERROR(msg)
1341 : : end if
1342 : : end if
1343 : : end do
1344 [ + + ]: 2425 : xmean = sum(xcart_perp)/natom
1345 [ + + ]: 2425 : xcart_perp(:)=xcart_perp(:)-xmean
1346 : :
1347 : : !Calculating the inverse (transpose) of the dielectric tensor
1348 : 561 : call matr3inv(dielt,invdlt)
1349 : : !Calculating the determinant of the dielectric tensor
1350 : : detdlt=dielt(1,1)*dielt(2,2)*dielt(3,3)+dielt(1,3)*dielt(2,1)*&
1351 : : & dielt(3,2)+dielt(1,2)*dielt(2,3)*dielt(3,1)-dielt(1,3)*&
1352 : : & dielt(2,2)*dielt(3,1)-dielt(1,1)*dielt(2,3)*dielt(3,2)-&
1353 : 561 : & dielt(1,2)*dielt(2,1)*dielt(3,3)
1354 : :
1355 [ - + ]: 561 : if(detdlt<tol6)then
1356 : : write(msg, '(a,es16.6,11a)' )&
1357 : 0 : 'The determinant of the dielectrix matrix, detdlt=',detdlt,' is smaller than 1.0d-6.',ch10,&
1358 : 0 : 'The use of the dipole-dipole model for interatomic force constants is not possible.',ch10,&
1359 : 0 : 'It is likely that you have not treated the electric field perturbations,',ch10,&
1360 : 0 : 'because you not are dealing with an insulator, so that',ch10,&
1361 : 0 : 'your dielectric matrix was simply set to zero in the Derivative DataBase.',ch10,&
1362 : 0 : 'Action: set the input variable dipdip to 0 .'
1363 : 0 : ABI_ERROR(msg)
1364 : : end if
1365 : :
1366 : : ! The dielectric tensor must be diagonal in the confined direction
1367 : 561 : rprimd_para = zero ; rprimd_perp = zero
1368 : 561 : dielt_para = zero ; dielt_perp = zero
1369 : 561 : invdlt_para = zero ; gprimd_para = zero
1370 : 561 : ndir=0
1371 [ + + ]: 2244 : do idir1=1,3
1372 [ + + ]: 1683 : if (periodic_dir(idir1)==1) then
1373 : 1122 : ndir=ndir+1
1374 : : end if
1375 : 1683 : mdir=0
1376 [ + + ]: 7293 : do idir2=1,3
1377 [ + + ]: 5049 : if (periodic_dir(idir2)==1) then
1378 : 3366 : mdir=mdir+1
1379 : : end if
1380 [ + + + + : 5049 : if ((periodic_dir(idir1)==1 .and. periodic_dir(idir2)==0) .or. &
+ + + + ]
1381 : 1683 : (periodic_dir(idir1)==0 .and. periodic_dir(idir2)==1)) then
1382 [ + - - + ]: 2244 : if (abs(dielt(idir1,idir2))>tol6 .or. abs(rprimd(idir1,idir2))>tol6) then !No cross in-plane out-of-plane are allowed
1383 : : write(msg, '(7a)' )&
1384 : 0 : 'The dielectric matrix shows off-diagonal components in the confined direction larger than 1d-6',ch10,&
1385 : 0 : 'This is forbidden when considering the Ewald summation for 2D systems. Please check if your', ch10, &
1386 : 0 : 'confined direction is correctly specified in the anaddb input or if the vacuum size if sufficiently', ch10, &
1387 : 0 : 'large in the confined direction to avoid spurious interactions between unit cells'
1388 : 0 : ABI_ERROR(msg)
1389 : : end if
1390 [ + + + + ]: 32068 : qdrp_paraperp(ndir,:,:) = qdrp_cart(idir1,idir2,:,:)
1391 [ + + + - ]: 2805 : elseif (periodic_dir(idir1)==0 .and. periodic_dir(idir2)==0) then
1392 : 561 : dielt_perp = dielt(idir1,idir2)
1393 : 561 : rprimd_perp = rprimd(idir1,idir2)
1394 [ + + + + ]: 8017 : qdrp_perpperp(:,:) = qdrp_cart(idir1,idir2,:,:)
1395 : : else
1396 : 2244 : dielt_para(ndir,mdir) = dielt(idir1,idir2)
1397 : 2244 : rprimd_para(ndir,mdir)=rprimd(idir1,idir2)
1398 : 2244 : gprimd_para(ndir,mdir) = gprimd(idir1,idir2)
1399 : 2244 : invdlt_para(ndir,mdir)=invdlt(idir1,idir2)
1400 [ + + + + ]: 32068 : qdrp_parapara(ndir,mdir,:,:) = qdrp_cart(idir1,idir2,:,:)
1401 : : end if
1402 : : end do
1403 : : end do
1404 : :
1405 : :
1406 : : ! Consistency check for dielectric thicknesses
1407 : 561 : inner_thick = thick(2)
1408 : 561 : out_thick = thick(1)
1409 : :
1410 [ + + ]: 2425 : do ipert1=1,natom
1411 [ + + ]: 9121 : do ipert2=1,natom
1412 [ + - - + ]: 8560 : if (xcart_perp(ipert1)-xcart_perp(ipert2)>out_thick .or. xcart_perp(ipert1)+xcart_perp(ipert2)>out_thick) then
1413 : : write(msg, '(5a)' )&
1414 : 0 : 'Some atoms seem to be located away from the dielectric slab.',ch10,&
1415 : 0 : 'The present model only allows for this specific scenario. Please increases slightly',ch10,&
1416 : 0 : 'the dielectric slab thickness'
1417 : 0 : ABI_ERROR(msg)
1418 : : end if
1419 : : end do
1420 : : end do
1421 [ + + ]: 561 : if (inner_thick>zero) then
1422 [ - + ]: 190 : if (inner_thick>out_thick) then
1423 : : write(msg, '(3a)' )&
1424 : 0 : 'When considering consecutive dielectric slab models, inner thickness',ch10,&
1425 : 0 : 'should be smaller than the outer thickness. Please check your input file'
1426 : 0 : ABI_ERROR(msg)
1427 : : else
1428 [ + + ]: 570 : do ipert1=1,natom
1429 [ + + ]: 1330 : do ipert2=1,natom
1430 [ + - - + ]: 1140 : if (xcart_perp(ipert1)-xcart_perp(ipert2)>inner_thick .or. xcart_perp(ipert1)+xcart_perp(ipert2)>inner_thick) then
1431 : : write(msg, '(5a)' )&
1432 : 0 : 'Some atoms seem to be located away from the inner dielectric slab.',ch10,&
1433 : 0 : 'The present model only allows for this specific scenario. Please increases slightly',ch10,&
1434 : 0 : 'the inner dielectric slab thickness'
1435 : 0 : ABI_ERROR(msg)
1436 : : end if
1437 : : end do
1438 : : end do
1439 : : end if
1440 : : end if
1441 : :
1442 : : ! First needs to determine the Gaussian broadening intrinsic to the Ewald summation. In 2D, both the real and
1443 : : ! reciprocal summation are related to the complementary error function. We want to restrict the real-part to
1444 : : ! the first Wigner cell. We use the fact that sqrt(1-e^{-x^2}) < erf(x) < sqrt(1-e^{-4x^2/pi})
1445 : : ! and invert those relationships to estimate the broadening required to restrict the real-part summation of
1446 : : !the Ewald summation; here fixes the threshold to 1e-9 for contribution from later unit cells
1447 : : rvec_dielt = zero
1448 : 561 : ndir=0
1449 [ + + ]: 1683 : do idir1=1,2
1450 [ + + + + : 10659 : norm_dielt(idir1) = dot_product(rprimd_para(idir1,:),matmul(invdlt_para(:,:),rprimd_para(idir1,:)))
+ + ]
1451 : : end do
1452 : :
1453 [ - + + - : 2244 : lambda = dsqrt(maxval(norm_dielt))/dsqrt(-two*dlog(one-(one-tol9)**2))
+ + ]
1454 : :
1455 : : ! Now that we have computed the value of lambda, we can compute the Ewald summation, starting from the
1456 : : ! reciprocal sum. We need first to estimate the max. number of reciprocal vectors we need to consider
1457 : : ! for the Ewald summation. We use here a stricter tolerance than for the determination of lambda
1458 : : ! Note that the worst case scenario is always when considering Rka=Rk'b in this case
1459 : :
1460 [ - + + - : 2244 : gmax = int(dsqrt(-two*dlog(one-(one-tol12)**2))/lambda/(two_pi/dsqrt(minval(norm_dielt))))
+ + ]
1461 [ + + ]: 561 : if (firstcall) then
1462 : 2 : firstcall = .FALSE.
1463 : 2 : write(msg, '(6a,f9.4,3a,i3,1a)' ) ch10,&
1464 : 2 : ' Ewald treatment of 2D long-range electrostatics interatomic force constants', ch10, &
1465 : 2 : ' To restrict the real-part summation of Ewald to the first unit cell, the Gaussian broadening', ch10, &
1466 : 2 : ' has been set to ', lambda, ' 1/Bohr. For the reciprocal sum, this corresponds to max.', ch10, &
1467 : 4 : ' ',2*gmax-1, ' Brillouin zone repetitions in either in-plane directions'
1468 [ + + ]: 6 : call wrtout([ab_out,std_out], msg)
1469 : : end if
1470 : :
1471 [ + + + + : 208897 : dyew_rec = zero
+ + + + +
+ ]
1472 : : ! If one dielectric slab model, same dielectric for both regions
1473 : : ! Otherwise, inner dielectric ~1 and the other has been computed
1474 : : ! accordingly in the anaddb driver
1475 [ + + ]: 561 : if (inner_thick> zero) then
1476 : : dielt_perp1 = one
1477 : : else
1478 : 371 : dielt_perp1=dielt_perp
1479 : : end if
1480 : 561 : dielt_perp2 = dielt_perp
1481 [ + + ]: 8216 : do ibz1 = -gmax,gmax
1482 [ + + ]: 114681 : do ibz2 = -gmax,gmax
1483 [ + + ]: 319395 : gvec(:) = ibz1*gprimd_para(:,1)+ibz2*gprimd_para(:,2)
1484 [ + + ]: 319395 : kvec(:) = gvec(:) + qvec_para(:)
1485 [ + + ]: 319395 : kvec(:) = kvec(:)*two_pi
1486 [ + + + + ]: 745255 : kvec_para(:) = matmul(dielt_para,kvec)
1487 [ + + ]: 319395 : norm_kvec0 = dot_product(kvec,kvec)
1488 [ + + ]: 319395 : norm_kvec = dot_product(kvec,kvec_para)
1489 [ + + ]: 114120 : if (abs(norm_kvec)>tol6) then !Remove G=q=0 case
1490 : 106456 : eta = dsqrt(norm_kvec/dielt_perp)
1491 : 106456 : eta1 = dsqrt(norm_kvec/dielt_perp1)
1492 : 106456 : xi = dsqrt(norm_kvec/dielt_perp2)
1493 : : ! Effective dielectric constants (depends on direction)
1494 : 106456 : dielt_eff = dsqrt(norm_kvec*dielt_perp/norm_kvec0)
1495 : 106456 : dielt_eff1 = dsqrt(norm_kvec*dielt_perp1/norm_kvec0)
1496 : 106456 : dielt_eff2 = dsqrt(norm_kvec*dielt_perp2/norm_kvec0)
1497 : : ! Reflection coefficient at the dielectric interfaces
1498 : 106456 : rflct_coeff = (dielt_eff-dielt_env)/(dielt_eff+dielt_env)
1499 : 106456 : rflct_coeff2 = (dielt_eff2-dielt_env)/(dielt_eff2+dielt_env)
1500 : 106456 : trans_fun = (one+rflct_coeff2*dexp(-xi*(out_thick-inner_thick)))
1501 : 106456 : trans_fun = trans_fun/(one-rflct_coeff2*dexp(-xi*(out_thick-inner_thick)))
1502 : 106456 : rflct_coeff1 = (dielt_eff1*trans_fun-dielt_eff2)/(dielt_eff1*trans_fun+dielt_eff2)
1503 : : ! Dipole-dipole charges prefactors
1504 : 106456 : fac_exp = rflct_coeff*dexp(-eta*out_thick)
1505 : 106456 : fac_exp1 = rflct_coeff1*dexp(-eta1*inner_thick)
1506 : 106456 : fac_mirror = two*fac_exp/(one-fac_exp**2)
1507 : 106456 : fac_mirror1 = two*fac_exp1/(one-fac_exp1**2)
1508 : : ! Ewald factor in error function
1509 : 106456 : fac_ewald1= eta*dsqrt(dielt_perp)*lambda/dsqrt(two)
1510 [ + + ]: 486310 : do ipert1=1,natom
1511 [ + + ]: 1913786 : do ipert2=1,natom
1512 : 1427476 : delta_perp = (xcart_perp(ipert2)-xcart_perp(ipert1))
1513 : 1427476 : mean2_perp = (xcart_perp(ipert2)+xcart_perp(ipert1))
1514 : 1427476 : fac_ewald2= delta_perp/lambda/sqrt(two*dielt_perp)
1515 : 1427476 : fac_ewald2b= delta_perp/lambda/sqrt(two*dielt_perp1)
1516 : : ! Ewald function and derivatives (eta factorized)
1517 : : ewald_fun = half*(dexp(-eta*delta_perp)*(one-erf(fac_ewald1-fac_ewald2)))+ &
1518 : 1427476 : half*(dexp(eta*delta_perp)*(one-erf(fac_ewald1+fac_ewald2)))
1519 : : ewald_fun1 = half*(-dexp(-eta*delta_perp)*(one-erf(fac_ewald1-fac_ewald2)))+ &
1520 : 1427476 : half*(dexp(eta*delta_perp)*(one-erf(fac_ewald1+fac_ewald2)))
1521 : : ! For second derivative, there is in principle a Gaussian term as well
1522 : : ! However, by an appropriate choice of the electrostatic gauge (mean average
1523 : : ! potential), we can neglect it. This approximation has been validated
1524 : : ! with respect to real-space dipoles and exact calculated points
1525 : : ewald_fun2 = half*(dexp(-eta1*delta_perp)*(one-erf(fac_ewald1-fac_ewald2b)))+ &
1526 : 1427476 : half*(dexp(eta1*delta_perp)*(one-erf(fac_ewald1+fac_ewald2b)))
1527 : 1427476 : ewald_fun2=-ewald_fun2*eta1**2
1528 : :
1529 : : ! Phase factor and mirror terms
1530 [ + + ]: 4282428 : phi = dot_product(kvec,xcart_para(:,ipert1)-xcart_para(:,ipert2))
1531 : 1427476 : mirror_parapara = fac_mirror*(dcosh(eta*mean2_perp)+fac_exp*dcosh(eta*delta_perp))
1532 [ + + ]: 1427476 : if (inner_thick >tol6) then
1533 : 91940 : mirror_perpperp = fac_mirror1*(dcosh(eta1*mean2_perp)-fac_exp1*dcosh(eta1*delta_perp))
1534 : : else
1535 : 1335536 : mirror_perpperp = fac_mirror*(dcosh(eta*mean2_perp)-fac_exp1*dcosh(eta*delta_perp))
1536 : : end if
1537 : 1427476 : mirror_paraperp = fac_mirror*fac_exp*dsinh(eta*delta_perp)
1538 : 1427476 : mirror_diff = fac_mirror*dsinh(eta*mean2_perp)
1539 : : ! Then compute the charge prefactor
1540 [ + + ]: 6089758 : do idir1=1,3
1541 [ + + ]: 18557188 : do idir2=1,3
1542 [ + + + + : 115625556 : rho_gerade1(1) = -half*dot_product(kvec,matmul(qdrp_parapara(:,:,idir1,ipert1),kvec))
+ + ]
1543 [ + + + + : 115625556 : rho_gerade2(1) = -half*dot_product(kvec,matmul(qdrp_parapara(:,:,idir2,ipert2),kvec))
+ + ]
1544 : 12847284 : rho_gerade1(1) = rho_gerade1(1)+half*eta1**2*qdrp_perpperp(idir1,ipert1)
1545 : 12847284 : rho_gerade2(1) = rho_gerade2(1)+half*eta1**2*qdrp_perpperp(idir2,ipert2)
1546 [ + + ]: 38541852 : rho_gerade1(2) = -dot_product(kvec,zeff_para(:,idir1,ipert1))
1547 [ + + ]: 38541852 : rho_gerade2(2) = -dot_product(kvec,zeff_para(:,idir2,ipert2))
1548 : 12847284 : rho_ungerade1(1) = -zeff_perp(idir1,ipert1)
1549 [ + + ]: 38541852 : rho_ungerade1(2) = eta1*dot_product(kvec,qdrp_paraperp(:,idir1,ipert1))
1550 : 12847284 : rho_ungerade2(1) = -zeff_perp(idir2,ipert2)
1551 [ + + ]: 38541852 : rho_ungerade2(2) = eta1*dot_product(kvec,qdrp_paraperp(:,idir2,ipert2))
1552 : : ! First, add the source charges (gerade gerade)
1553 : : dyew_rec(1,idir1,ipert1,idir2,ipert2)= dyew_rec(1,idir1,ipert1,idir2,ipert2)+&
1554 : : (rho_gerade1(1)*rho_gerade2(1)+rho_gerade1(2)*rho_gerade2(2))&
1555 : 12847284 : *ewald_fun/eta/dielt_perp*cos(phi)
1556 : : dyew_rec(1,idir1,ipert1,idir2,ipert2)= dyew_rec(1,idir1,ipert1,idir2,ipert2)-&
1557 : : (rho_gerade1(2)*rho_gerade2(1)-rho_gerade1(1)*rho_gerade2(2))&
1558 : 12847284 : *ewald_fun/eta/dielt_perp*sin(phi)
1559 : : dyew_rec(2,idir1,ipert1,idir2,ipert2)= dyew_rec(2,idir1,ipert1,idir2,ipert2)+&
1560 : : (rho_gerade1(1)*rho_gerade2(1)+rho_gerade1(2)*rho_gerade2(2))&
1561 : 12847284 : *ewald_fun/eta/dielt_perp*sin(phi)
1562 : : dyew_rec(2,idir1,ipert1,idir2,ipert2)= dyew_rec(2,idir1,ipert1,idir2,ipert2)+&
1563 : : (rho_gerade1(2)*rho_gerade2(1)-rho_gerade1(1)*rho_gerade2(2))&
1564 : 12847284 : *ewald_fun/eta/dielt_perp*cos(phi)
1565 : : ! Second, add the source charges (ungerade ungerade)
1566 : : dyew_rec(1,idir1,ipert1,idir2,ipert2)= dyew_rec(1,idir1,ipert1,idir2,ipert2)+&
1567 : : (rho_ungerade1(1)*rho_ungerade2(1)+rho_ungerade1(2)*rho_ungerade2(2))&
1568 : 12847284 : *ewald_fun2/eta1/dielt_perp1*cos(phi)
1569 : : dyew_rec(1,idir1,ipert1,idir2,ipert2)= dyew_rec(1,idir1,ipert1,idir2,ipert2)-&
1570 : : (rho_ungerade1(2)*rho_ungerade2(1)-rho_ungerade1(1)*rho_ungerade2(2))&
1571 : 12847284 : *ewald_fun2/eta1/dielt_perp1*sin(phi)
1572 : : dyew_rec(2,idir1,ipert1,idir2,ipert2)= dyew_rec(2,idir1,ipert1,idir2,ipert2)+&
1573 : : (rho_ungerade1(1)*rho_ungerade2(1)+rho_ungerade1(2)*rho_ungerade2(2))&
1574 : 12847284 : *ewald_fun2/eta1/dielt_perp1*sin(phi)
1575 : : dyew_rec(2,idir1,ipert1,idir2,ipert2)= dyew_rec(2,idir1,ipert1,idir2,ipert2)+&
1576 : : (rho_ungerade1(2)*rho_ungerade2(1)-rho_ungerade1(1)*rho_ungerade2(2))&
1577 : 12847284 : *ewald_fun2/eta1/dielt_perp1*cos(phi)
1578 : : ! Third, add the source charge (gerade ungerade)
1579 : : ! For sake of consistenty, only used when there is only one dielectric thickness
1580 [ + + ]: 12847284 : if (inner_thick <tol6) then
1581 : : dyew_rec(1,idir1,ipert1,idir2,ipert2)= dyew_rec(1,idir1,ipert1,idir2,ipert2)-&
1582 : : (rho_gerade1(1)*rho_ungerade2(1)+rho_ungerade1(1)*rho_gerade2(1)&
1583 : : +rho_gerade1(2)*rho_ungerade2(2)+rho_ungerade1(2)*rho_gerade2(2)) &
1584 : 12019824 : *ewald_fun1/eta/dielt_perp*cos(phi)
1585 : : dyew_rec(1,idir1,ipert1,idir2,ipert2)= dyew_rec(1,idir1,ipert1,idir2,ipert2)+&
1586 : : (rho_gerade1(2)*rho_ungerade2(1)+rho_ungerade1(2)*rho_gerade2(1)&
1587 : : -rho_gerade1(1)*rho_ungerade2(2)-rho_ungerade1(1)*rho_gerade2(2)) &
1588 : 12019824 : *ewald_fun1/eta/dielt_perp*sin(phi)
1589 : : dyew_rec(2,idir1,ipert1,idir2,ipert2)= dyew_rec(2,idir1,ipert1,idir2,ipert2)-&
1590 : : (rho_gerade1(1)*rho_ungerade2(1)+rho_ungerade1(1)*rho_gerade2(1)&
1591 : : +rho_gerade1(2)*rho_ungerade2(2)+rho_ungerade1(2)*rho_gerade2(2)) &
1592 : 12019824 : *ewald_fun1/eta/dielt_perp*sin(phi)
1593 : : dyew_rec(2,idir1,ipert1,idir2,ipert2)= dyew_rec(2,idir1,ipert1,idir2,ipert2)-&
1594 : : (rho_gerade1(2)*rho_ungerade2(1)+rho_ungerade1(2)*rho_gerade2(1)&
1595 : : -rho_gerade1(1)*rho_ungerade2(2)-rho_ungerade1(1)*rho_gerade2(2)) &
1596 : 12019824 : *ewald_fun1/eta/dielt_perp*cos(phi)
1597 : : end if
1598 : : ! Now add the interactions with the mirror charges... para para
1599 : : dyew_rec(1,idir1,ipert1,idir2,ipert2)= dyew_rec(1,idir1,ipert1,idir2,ipert2)+&
1600 : : (rho_gerade1(1)*rho_gerade2(1)+rho_gerade1(2)*rho_gerade2(2))&
1601 : 12847284 : *mirror_parapara/eta/dielt_perp*cos(phi)
1602 : : dyew_rec(1,idir1,ipert1,idir2,ipert2)= dyew_rec(1,idir1,ipert1,idir2,ipert2)-&
1603 : : (rho_gerade1(2)*rho_gerade2(1)-rho_gerade1(1)*rho_gerade2(2))&
1604 : 12847284 : *mirror_parapara/eta/dielt_perp*sin(phi)
1605 : : dyew_rec(2,idir1,ipert1,idir2,ipert2)= dyew_rec(2,idir1,ipert1,idir2,ipert2)+&
1606 : : (rho_gerade1(1)*rho_gerade2(1)+rho_gerade1(2)*rho_gerade2(2))&
1607 : 12847284 : *mirror_parapara/eta/dielt_perp*sin(phi)
1608 : : dyew_rec(2,idir1,ipert1,idir2,ipert2)= dyew_rec(2,idir1,ipert1,idir2,ipert2)+&
1609 : : (rho_gerade1(2)*rho_gerade2(1)-rho_gerade1(1)*rho_gerade2(2))&
1610 : 12847284 : *mirror_parapara/eta/dielt_perp*cos(phi)
1611 : : ! Now with perp perp
1612 : : dyew_rec(1,idir1,ipert1,idir2,ipert2)= dyew_rec(1,idir1,ipert1,idir2,ipert2)+&
1613 : : (rho_ungerade1(1)*rho_ungerade2(1)+rho_ungerade1(2)*rho_ungerade2(2))&
1614 : 12847284 : *mirror_perpperp*eta1/dielt_perp1*cos(phi)
1615 : : dyew_rec(1,idir1,ipert1,idir2,ipert2)= dyew_rec(1,idir1,ipert1,idir2,ipert2)-&
1616 : : (rho_ungerade1(2)*rho_ungerade2(1)-rho_ungerade1(1)*rho_ungerade2(2))&
1617 : 12847284 : *mirror_perpperp*eta1/dielt_perp1*sin(phi)
1618 : : dyew_rec(2,idir1,ipert1,idir2,ipert2)= dyew_rec(2,idir1,ipert1,idir2,ipert2)+&
1619 : : (rho_ungerade1(1)*rho_ungerade2(1)+rho_ungerade1(2)*rho_ungerade2(2))&
1620 : 12847284 : *mirror_perpperp*eta1/dielt_perp1*sin(phi)
1621 : : dyew_rec(2,idir1,ipert1,idir2,ipert2)= dyew_rec(2,idir1,ipert1,idir2,ipert2)+&
1622 : : (rho_ungerade1(2)*rho_ungerade2(1)-rho_ungerade1(1)*rho_ungerade2(2))&
1623 : 12847284 : *mirror_perpperp*eta1/dielt_perp1*cos(phi)
1624 : : ! Third, add the source charge (gerade ungerade)
1625 [ + + ]: 17129712 : if (inner_thick <tol6) then
1626 : : dyew_rec(1,idir1,ipert1,idir2,ipert2)= dyew_rec(1,idir1,ipert1,idir2,ipert2)+&
1627 : : (rho_gerade1(1)*rho_ungerade2(1)+rho_ungerade1(1)*rho_gerade2(1)&
1628 : : +rho_gerade1(2)*rho_ungerade2(2)+rho_ungerade1(2)*rho_gerade2(2)) &
1629 : 12019824 : *mirror_paraperp/eta/dielt_perp*cos(phi)
1630 : : dyew_rec(1,idir1,ipert1,idir2,ipert2)= dyew_rec(1,idir1,ipert1,idir2,ipert2)+&
1631 : : (rho_gerade1(2)*rho_ungerade2(1)+rho_ungerade1(2)*rho_gerade2(1)&
1632 : : -rho_gerade1(1)*rho_ungerade2(2)-rho_ungerade1(1)*rho_gerade2(2)) &
1633 : 12019824 : *mirror_paraperp/eta/dielt_perp*sin(phi)
1634 : : dyew_rec(2,idir1,ipert1,idir2,ipert2)= dyew_rec(2,idir1,ipert1,idir2,ipert2)-&
1635 : : (rho_gerade1(1)*rho_ungerade2(1)+rho_ungerade1(1)*rho_gerade2(1)&
1636 : : +rho_gerade1(2)*rho_ungerade2(2)+rho_ungerade1(2)*rho_gerade2(2)) &
1637 : 12019824 : *mirror_paraperp/eta/dielt_perp*sin(phi)
1638 : : dyew_rec(2,idir1,ipert1,idir2,ipert2)= dyew_rec(2,idir1,ipert1,idir2,ipert2)+&
1639 : : (rho_gerade1(2)*rho_ungerade2(1)+rho_ungerade1(2)*rho_gerade2(1)&
1640 : : -rho_gerade1(1)*rho_ungerade2(2)-rho_ungerade1(1)*rho_gerade2(2)) &
1641 : 12019824 : *mirror_paraperp/eta/dielt_perp*cos(phi)
1642 : : ! Finally, there is a term on the sum of charge, only for mirror charges
1643 : : ! Third, add the source charge (gerade ungerade)
1644 : : dyew_rec(1,idir1,ipert1,idir2,ipert2)= dyew_rec(1,idir1,ipert1,idir2,ipert2)-&
1645 : : (rho_gerade1(1)*rho_ungerade2(1)-rho_ungerade1(1)*rho_gerade2(1)&
1646 : : +rho_gerade1(2)*rho_ungerade2(2)-rho_ungerade1(2)*rho_gerade2(2)) &
1647 : 12019824 : *mirror_diff/eta/dielt_perp*cos(phi)
1648 : : dyew_rec(1,idir1,ipert1,idir2,ipert2)= dyew_rec(1,idir1,ipert1,idir2,ipert2)-&
1649 : : (rho_gerade1(2)*rho_ungerade2(1)-rho_ungerade1(2)*rho_gerade2(1)&
1650 : : -rho_gerade1(1)*rho_ungerade2(2)+rho_ungerade1(1)*rho_gerade2(2)) &
1651 : 12019824 : *mirror_diff/eta/dielt_perp*sin(phi)
1652 : : dyew_rec(2,idir1,ipert1,idir2,ipert2)= dyew_rec(2,idir1,ipert1,idir2,ipert2)+&
1653 : : (rho_gerade1(1)*rho_ungerade2(1)-rho_ungerade1(1)*rho_gerade2(1)&
1654 : : +rho_gerade1(2)*rho_ungerade2(2)-rho_ungerade1(2)*rho_gerade2(2)) &
1655 : 12019824 : *mirror_diff/eta/dielt_perp*sin(phi)
1656 : : dyew_rec(2,idir1,ipert1,idir2,ipert2)= dyew_rec(2,idir1,ipert1,idir2,ipert2)-&
1657 : : (rho_gerade1(2)*rho_ungerade2(1)-rho_ungerade1(2)*rho_gerade2(1)&
1658 : : +rho_gerade1(1)*rho_ungerade2(2)+rho_ungerade1(1)*rho_gerade2(2)) &
1659 : 12019824 : *mirror_diff/eta/dielt_perp*sin(phi)
1660 : : end if
1661 : : end do
1662 : : end do
1663 : : end do
1664 : : end do
1665 : : end if
1666 : : end do
1667 : : end do
1668 : 561 : ucsurf = rprimd_para(1,1)*rprimd_para(2,2)-rprimd_para(1,2)*rprimd_para(2,1)
1669 : : ! Renormalize by surface of periodic 2D lattice and out-of-plane dielectric constant
1670 [ + + + + : 208897 : dyew_rec = dyew_rec*(two_pi)/ucsurf
+ + + + +
+ ]
1671 : : ! Reciprocal summation completes. Remains some real-space contribution from first
1672 : : ! unit cells (ipert1 neq ipert2), impacting all IFCs the same way
1673 [ + + + + : 208897 : dyew_real = zero
+ + + + +
+ ]
1674 [ + + ]: 2425 : do ipert1=1,natom
1675 [ + + ]: 9121 : do ipert2=1,natom
1676 [ + + ]: 26784 : diff_xcart(:) =xcart(:,ipert2)-xcart(:,ipert1)
1677 [ + + + + ]: 87048 : rvec_dielt(:) = matmul(invdlt,diff_xcart)
1678 [ + + ]: 26784 : rvec_norm = dot_product(diff_xcart,rvec_dielt)
1679 : 6696 : sqrt_norm = dsqrt(rvec_norm)
1680 : 6696 : fac_erfc = sqrt_norm/dsqrt(two)/lambda
1681 : 6696 : fac_gauss = -rvec_norm/two/lambda**2
1682 [ + + ]: 28648 : do idir1=1,3
1683 [ + + ]: 87048 : do idir2=1,3
1684 [ + + ]: 80352 : if (ipert1 .NE. ipert2) then ! Only off-sites contributions
1685 : : ! First, contribution from eps^-1 (Rk'b-Rka) eps^-1
1686 : : fac_real = three*(one-erf(fac_erfc))/sqrt_norm**5+6*exp(fac_gauss)/rvec_norm**2/&
1687 : 43488 : sqrt(two_pi)/lambda+two*exp(fac_gauss)/rvec_norm/sqrt(two_pi)/lambda**3
1688 : : dyew_real(1,idir1,ipert1,idir2,ipert2) = dyew_real(1,idir1,ipert1,idir2,ipert2)+&
1689 : 43488 : fac_real*rvec_dielt(idir1)*rvec_dielt(idir2)
1690 : : ! Second contribution from esp^-1(alpha,beta)
1691 : 43488 : fac_real = (one-erf(fac_erfc))/sqrt_norm**3+exp(fac_gauss)/rvec_norm/sqrt(two_pi)/lambda
1692 : : dyew_real(1,idir1,ipert1,idir2,ipert2)=dyew_real(1,idir1,ipert1,idir2,ipert2)-&
1693 : 43488 : fac_real*dielt(idir1,idir2)
1694 : : end if
1695 : : end do
1696 : : end do
1697 : : end do
1698 : : end do
1699 [ + + + + : 208897 : dyew_real = dyew_real / dsqrt(detdlt)
+ + + + +
+ ]
1700 [ + + + + : 208897 : dyew = dyew_real + dyew_rec
+ + + + +
+ ]
1701 : 561 : end subroutine ewald9_2D
1702 : :
1703 : : end module m_ewald
1704 : : !!***
1705 : :
|
