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 : 7665 : subroutine ewald(eew,gmet,grewtn,gsqcut,icutcoul,natom,ngfft,nkpt,ntypat,rcut,&
76 : 7665 : 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 : 7665 : real(dp),allocatable :: gcutoff(:)
99 : : ! *************************************************************************
100 : :
101 : : !This is the minimum argument of an exponential, with some safety
102 : 7665 : minexparg=log(tiny(0._dp))+five
103 : :
104 : : !Add up total charge and sum of $charge^2$ in cell
105 : :
106 : 7665 : chsq=0._dp
107 : 7665 : ch=0._dp
108 [ + + ]: 30053 : do ia=1,natom
109 : 22388 : ch=ch+zion(typat(ia))
110 : 30053 : 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 : 7665 : 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 : 7665 : 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 [ - + ]: 7665 : 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 : 7665 : eta=pi*200.0_dp/33.0_dp*sqrt(1.69_dp*recip/direct)
139 : : end if
140 : :
141 : : !Conduct reciprocal space summations
142 : 7665 : fac=pi**2/eta
143 : 7665 : gsum=0._dp
144 [ + + + + ]: 97217 : grewtn(:,:)=0.0_dp
145 : :
146 : : !Initialize Gcut-off array from m_gtermcutoff
147 : : !ABI_MALLOC(gcutoff,(ngfft(1)*ngfft(2)*ngfft(3)))
148 : 7665 : 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 : 7665 : ng=0
154 : : do
155 : 319278 : ng=ng+1
156 : 319278 : 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 : 319278 : ii=1
165 [ + + ]: 15926878 : do ig3=-ng,ng
166 [ + + ]: 1242636236 : do ig2=-ng,ng
167 [ + + ]: >18135*10^7 : do ig1=-ng,ng
168 : : ! Exclude shells previously summed over
169 [ + + + + : >18133*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 : 7175526837 : & gmet(3,1)*dble(ig1*ig3)+gmet(3,2)*dble(ig3*ig2))
175 : :
176 : : ! Skip g=0:
177 [ + + ]: 7175526837 : if (gsq>1.0d-20) then
178 : 7175519172 : arg=fac*gsq
179 : :
180 : : ! Larger arg gives 0 contribution because of exp(-arg)
181 [ + + ]: 7175519172 : if (arg <= -minexparg ) then
182 : : ! When any term contributes then include next shell
183 : 1845811084 : newg=1
184 : :
185 : : if((abs(ig1).lt.ngfft(1)).and.&
186 [ + + + + : 1845811084 : &(abs(ig2).lt.ngfft(2)).and.&
+ + ]
187 : : &(abs(ig3).lt.ngfft(3))) then
188 : 568867324 : ig23=ngfft(1)*(abs(ig2)+ngfft(2)*(abs(ig3)))
189 : 568867324 : 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 : 568867324 : term=exp(-arg)/gsq*gcutoff(ii)
194 [ + + ]: 1276943760 : else if (icutcoul.ne.3) then
195 : : term=zero !exp(-arg)/gsq
196 : : else
197 : 1276721714 : term=exp(-arg)/gsq
198 : : endif
199 : :
200 : 1845811084 : summr = 0.0_dp
201 : 1845811084 : 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 [ + + ]: 7600368796 : do ia=1,natom
211 : 5754557712 : 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 : 5754557712 : summr=summr+zion(typat(ia))*cos(arg)
214 : 7600368796 : summi=summi+zion(typat(ia))*sin(arg)
215 : : end do
216 : :
217 : : ! The following two checks avoid an annoying underflow error msg
218 [ + + ]: 1845811084 : if (abs(summr)<1.d-16) summr=0.0_dp
219 [ + + ]: 1845811084 : 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 : 1845811084 : t1=term*(summr*summr+summi*summi)
224 : 1845811084 : gsum=gsum+t1
225 : :
226 [ + + ]: 7600368796 : do ia=1,natom
227 : : ! Again only phase is computed so xred may fall outside [0,1).
228 : 5754557712 : arg=two_pi*(ig1*xred(1,ia)+ig2*xred(2,ia)+ig3*xred(3,ia))
229 : 5754557712 : phr= cos(arg)
230 : 5754557712 : phi=-sin(arg)
231 : : ! (note: do not need real part, commented out)
232 : : ! c1r=(phr*summr-phi*summi)*(term*zion(typat(ia)))
233 : 5754557712 : c1i=(phi*summr+phr*summi)*(term*zion(typat(ia)))
234 : : ! compute coordinate gradients
235 : 5754557712 : grewtn(1,ia)=grewtn(1,ia)-c1i*ig1
236 : 5754557712 : grewtn(2,ia)=grewtn(2,ia)-c1i*ig2
237 : 7600368796 : 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 [ + + ]: 319278 : 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 : 7665 : 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 : 7665 : term=-2._dp/ucvol
261 [ + + + + ]: 97217 : grewtn(:,:)=grewtn(:,:)*term
262 : : !call DSCAL(3*natom,term,grewtn,1)
263 : :
264 : : !Conduct real space summations
265 : 7665 : reta=sqrt(eta)
266 : 7665 : fac=2._dp*sqrt(eta/pi)
267 : 7665 : 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 : 7665 : nr=0
275 : : do
276 : 22300 : nr=nr+1
277 : 22300 : 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 [ + + ]: 135288 : do ir3=-nr,nr
287 [ + + ]: 784580 : do ir2=-nr,nr
288 [ + + ]: 4923540 : do ir1=-nr,nr
289 [ + + + + : 4810552 : if( abs(ir3)==nr .or. abs(ir2)==nr .or. abs(ir1)==nr .or. nr==1 )then
+ + + + ]
290 : :
291 [ + + ]: 11497773 : do ia=1,natom
292 : : ! Map reduced coordinate xred(mu,ia) into [0,1)
293 : 8771812 : fraca1=xred(1,ia)-aint(xred(1,ia))+0.5_dp-sign(0.5_dp,xred(1,ia))
294 : 8771812 : fraca2=xred(2,ia)-aint(xred(2,ia))+0.5_dp-sign(0.5_dp,xred(2,ia))
295 : 8771812 : fraca3=xred(3,ia)-aint(xred(3,ia))+0.5_dp-sign(0.5_dp,xred(3,ia))
296 : 8771812 : drdta1=0.0_dp
297 : 8771812 : drdta2=0.0_dp
298 : 8771812 : drdta3=0.0_dp
299 : :
300 [ + + ]: 103224750 : 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 : 94452938 : fracb1=xred(1,ib)-aint(xred(1,ib))+0.5_dp-sign(0.5_dp,xred(1,ib))
304 : 94452938 : fracb2=xred(2,ib)-aint(xred(2,ib))+0.5_dp-sign(0.5_dp,xred(2,ib))
305 : 94452938 : fracb3=xred(3,ib)-aint(xred(3,ib))+0.5_dp-sign(0.5_dp,xred(3,ib))
306 : 94452938 : r1=dble(ir1)+fracb1-fraca1
307 : 94452938 : r2=dble(ir2)+fracb2-fraca2
308 : 94452938 : r3=dble(ir3)+fracb3-fraca3
309 : : rsq=rmet(1,1)*r1*r1+rmet(2,2)*r2*r2+rmet(3,3)*r3*r3+&
310 : 94452938 : & 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 [ + + ]: 103224750 : 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 : 94430550 : term=0._dp
318 [ + + ]: 94430550 : if (eta*rsq<64.0_dp) then
319 : 5554160 : newr=1
320 : 5554160 : rmagn=sqrt(rsq)
321 : 5554160 : arg=reta*rmagn
322 : : ! derfc is the real(dp) complementary error function
323 : 5554160 : derfc_arg = abi_derfc(arg)
324 : 5554160 : term=derfc_arg/rmagn
325 : 5554160 : sumr=sumr+zion(typat(ia))*zion(typat(ib))*term
326 : : term=zion(typat(ia))*zion(typat(ib))*&
327 : 5554160 : & (term+fac*exp(-eta*rsq))/rsq
328 : : ! Length scale grads now handled with stress tensor in ewald2
329 : 5554160 : r1a1d=rmet(1,1)*r1+rmet(1,2)*r2+rmet(1,3)*r3
330 : 5554160 : r2a2d=rmet(2,1)*r1+rmet(2,2)*r2+rmet(2,3)*r3
331 : 5554160 : r3a3d=rmet(3,1)*r1+rmet(3,2)*r2+rmet(3,3)*r3
332 : : ! Compute terms related to coordinate gradients
333 : 5554160 : drdta1=drdta1+term*r1a1d
334 : 5554160 : drdta2=drdta2+term*r2a2d
335 : 5554160 : 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 : 8771812 : grewtn(1,ia)=grewtn(1,ia)+drdta1
341 : 8771812 : grewtn(2,ia)=grewtn(2,ia)+drdta2
342 : 12933072 : 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 [ + + ]: 22300 : if(newr==0) exit
351 : : end do ! End loop on nr (new shells). Note that there is an exit within the loop
352 : : !
353 : 7665 : sumr=0.5_dp*sumr
354 : 7665 : fac=pi*ch**2.0_dp/(2.0_dp*eta*ucvol)
355 : :
356 : : !Finally assemble Ewald energy, eew
357 [ + + ]: 7665 : 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 : 7655 : eew=sumg+sumr-chsq*reta/sqrt(pi)-fac
362 : : end if
363 : :
364 [ - + ]: 7665 : 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 : 7665 : 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 : 6319 : 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 : 6319 : call matr3inv(rprimd,gprimd)
439 : :
440 : : !This is the minimum argument of an exponential, with some safety
441 : 6319 : minexparg=log(tiny(0._dp))+five
442 : :
443 : : !Add up total charge and sum of charge^2 in cell
444 : 6319 : ch=0._dp
445 [ + + ]: 25007 : do ia=1,natom
446 : 25007 : 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 : 6319 : & 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 : 6319 : & 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 : 6319 : eta=pi*200.0_dp/33.0_dp*sqrt(1.69_dp*recip/direct)
458 : :
459 : 6319 : fac=pi**2/eta
460 : :
461 : : !Conduct reciprocal space summations
462 : 6319 : strg(1:6)=0.0_dp
463 : :
464 : : !Sum over G space, done shell after shell until all
465 : : !contributions are too small
466 : 6319 : ng=0
467 : : do
468 : 270606 : ng=ng+1
469 : 270606 : newg=0
470 : :
471 [ + + ]: 13825250 : do ig3=-ng,ng
472 [ + + ]: 1111041344 : do ig2=-ng,ng
473 [ + + ]: >17022*10^7 : do ig1=-ng,ng
474 : :
475 : : ! Exclude shells previously summed over
476 [ + + + + : >17021*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 : 6422812003 : g1=gprimd(1,1)*ig1+gprimd(1,2)*ig2+gprimd(1,3)*ig3
481 : 6422812003 : g2=gprimd(2,1)*ig1+gprimd(2,2)*ig2+gprimd(2,3)*ig3
482 : 6422812003 : g3=gprimd(3,1)*ig1+gprimd(3,2)*ig2+gprimd(3,3)*ig3
483 : : ! Compute |G|^2 (no pi factors)
484 : 6422812003 : gsq=(g1**2+g2**2+g3**2)
485 : :
486 : : ! skip g=0:
487 [ + + ]: 6422812003 : if (gsq>1.0d-20) then
488 : 6422805684 : arg1=fac*gsq
489 : :
490 : : ! larger arg1 gives 0 contribution because of exp(-arg1)
491 [ + + ]: 6422805684 : if (arg1<= -minexparg) then
492 : : ! When any term contributes then include next shell
493 : 1645361582 : newg=1
494 : 1645361582 : term1=exp(-arg1)/arg1
495 : 1645361582 : summr = 0.0_dp
496 : 1645361582 : summi = 0.0_dp
497 [ + + ]: 6754237832 : do ia=1,natom
498 : 5108876250 : 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 : 5108876250 : summr=summr+zion(typat(ia))*cos(arg2)
501 : 6754237832 : summi=summi+zion(typat(ia))*sin(arg2)
502 : : end do
503 : :
504 : : ! Avoid underflow error messages
505 [ + + ]: 1645361582 : if (abs(summr)<1.d-16) summr=0.0_dp
506 [ + + ]: 1645361582 : if (abs(summi)<1.d-16) summi=0.0_dp
507 : :
508 : 1645361582 : term2=(2._dp/gsq)*(1._dp+arg1)
509 : 1645361582 : t1=term2*g1*g1-1._dp
510 : 1645361582 : t2=term2*g2*g2-1._dp
511 : 1645361582 : t3=term2*g3*g3-1._dp
512 : 1645361582 : t4=term2*g2*g3
513 : 1645361582 : t5=term2*g1*g3
514 : 1645361582 : t6=term2*g1*g2
515 : 1645361582 : term3=term1*(summr*summr+summi*summi)
516 : 1645361582 : strg(1)=strg(1)+t1*term3
517 : 1645361582 : strg(2)=strg(2)+t2*term3
518 : 1645361582 : strg(3)=strg(3)+t3*term3
519 : 1645361582 : strg(4)=strg(4)+t4*term3
520 : 1645361582 : strg(5)=strg(5)+t5*term3
521 : 1645361582 : 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 [ + + ]: 270606 : 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 : 6319 : reta=sqrt(eta)
538 : 6319 : strr(1:6)=0.0_dp
539 : :
540 : : !Loop on shells in r-space as was done in g-space
541 : 6319 : nr=0
542 : : do
543 : 18181 : nr=nr+1
544 : 18181 : newr=0
545 : :
546 [ + + ]: 110002 : do ir3=-nr,nr
547 [ + + ]: 636967 : do ir2=-nr,nr
548 [ + + ]: 4004079 : do ir1=-nr,nr
549 [ + + + + : 3912258 : if( abs(ir3)==nr .or. abs(ir2)==nr .or. abs(ir1)==nr .or. nr==1 )then
+ + + + ]
550 : :
551 [ + + ]: 9422697 : do ia=1,natom
552 : : ! Convert reduced atomic coordinates to [0,1)
553 : 7210992 : fraca1=xred(1,ia)-aint(xred(1,ia))+0.5_dp-sign(0.5_dp,xred(1,ia))
554 : 7210992 : fraca2=xred(2,ia)-aint(xred(2,ia))+0.5_dp-sign(0.5_dp,xred(2,ia))
555 : 7210992 : fraca3=xred(3,ia)-aint(xred(3,ia))+0.5_dp-sign(0.5_dp,xred(3,ia))
556 [ + + ]: 92774403 : do ib=1,natom
557 : 82178118 : fracb1=xred(1,ib)-aint(xred(1,ib))+0.5_dp-sign(0.5_dp,xred(1,ib))
558 : 82178118 : fracb2=xred(2,ib)-aint(xred(2,ib))+0.5_dp-sign(0.5_dp,xred(2,ib))
559 : 82178118 : fracb3=xred(3,ib)-aint(xred(3,ib))+0.5_dp-sign(0.5_dp,xred(3,ib))
560 : 82178118 : r1=ir1+fracb1-fraca1
561 : 82178118 : r2=ir2+fracb2-fraca2
562 : 82178118 : r3=ir3+fracb3-fraca3
563 : : ! Convert from reduced to cartesian coordinates
564 : 82178118 : r1c=rprimd(1,1)*r1+rprimd(1,2)*r2+rprimd(1,3)*r3
565 : 82178118 : r2c=rprimd(2,1)*r1+rprimd(2,2)*r2+rprimd(2,3)*r3
566 : 82178118 : r3c=rprimd(3,1)*r1+rprimd(3,2)*r2+rprimd(3,3)*r3
567 : : ! Compute |r|^2
568 : 82178118 : rsq=r1c**2+r2c**2+r3c**2
569 : 82178118 : rmagn=sqrt(rsq)
570 : :
571 : : ! Avoid zero denominators in 'term':
572 [ + + ]: 89389110 : 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 : 82159430 : arg3=reta*rmagn
577 [ + + ]: 82159430 : if (arg3<8.0_dp) then
578 : 4505916 : newr=1
579 : : ! derfc computes the complementary error function
580 : : ! dderfc is the derivative of the complementary error function
581 : 4505916 : dderfc=(-2/sqrt(pi))*exp(-eta*rsq)
582 : 4505916 : derfc_arg = abi_derfc(arg3)
583 : 4505916 : term3=dderfc-derfc_arg/arg3
584 : 4505916 : term4=zion(typat(ia))*zion(typat(ib))*term3
585 : 4505916 : strr(1)=strr(1)+term4*r1c*r1c/rsq
586 : 4505916 : strr(2)=strr(2)+term4*r2c*r2c/rsq
587 : 4505916 : strr(3)=strr(3)+term4*r3c*r3c/rsq
588 : 4505916 : strr(4)=strr(4)+term4*r2c*r3c/rsq
589 : 4505916 : strr(5)=strr(5)+term4*r1c*r3c/rsq
590 : 4505916 : 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 [ + + ]: 18181 : 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 : 6319 : fac = pi/(2._dp*ucvol*eta)
611 : 6319 : stress(1)=(0.5_dp*reta*strr(1)+fac*(strg(1)+(ch**2)))/ucvol
612 : 6319 : stress(2)=(0.5_dp*reta*strr(2)+fac*(strg(2)+(ch**2)))/ucvol
613 : 6319 : stress(3)=(0.5_dp*reta*strr(3)+fac*(strg(3)+(ch**2)))/ucvol
614 : 6319 : stress(4)=(0.5_dp*reta*strr(4)+fac*strg(4))/ucvol
615 : 6319 : stress(5)=(0.5_dp*reta*strr(5)+fac*strg(5))/ucvol
616 : 6319 : stress(6)=(0.5_dp*reta*strr(6)+fac*strg(6))/ucvol
617 : :
618 : 6319 : 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 : : !! dielt(3,3)=dielectric tensor
634 : : !! gmet(3,3) = metric in reciprocal space.
635 : : !! gprimd(3,3)=dimensional primitive translations in reciprocal space (bohr^-1)
636 : : !! natom=number of atoms in unit cell
637 : : !! qphon(3)=phonon wavevector (same system of coordinates as the reciprocal lattice vectors)
638 : : !! rmet = metric in real space
639 : : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
640 : : !! sumg0: if=1, the sum in reciprocal space must include g=0,
641 : : !! if=0, this contribution must be skipped (q=0 singularity)
642 : : !! ucvol=unit cell volume in (whatever length scale units)**3
643 : : !! xred(3,natom)=relative coords of atoms in unit cell (dimensionless)
644 : : !! zeff(3,3,natom)=effective charge on each atom, versus electric
645 : : !! field and atomic displacement
646 : : !! qdrp_cart(3,3,3,natom)=Quadrupole tensor on each atom in cartesian cordinates
647 : : !! option= 0: use old implementation;
648 : : !! 1: reduce the smalest argument of the exponentials to be evaluated,
649 : : !! set eta to 1 and skip real space sum, leads to a significant speedup
650 : : !! [dipquad] = if 1, atmfrc has been build without dipole-quadrupole part
651 : : !! [quadquad] = if 1, atmfrc has been build without quadrupole-quadrupole part
652 : : !!
653 : : !! OUTPUT
654 : : !! dyew(2,3,natom,3,natom)= Ewald part of the dynamical matrix,
655 : : !! second energy derivative wrt xred(3,natom) in Hartrees
656 : : !! Set to zero if all(zeff == zero)
657 : : !! eta: parameter used to split R and G-space summation
658 : : !!
659 : : !! NOTES
660 : : !! 1. The q=0 part should be subtracted, by another call to
661 : : !! the present routine, with q=0. The present routine correspond
662 : : !! to the quantity written A-bar in the explanatory notes.
663 : : !! If q=0 is asked, sumg0 should be put to 0. Otherwise, it should be put to 1.
664 : : !! 2. Because this routine can be used many times in the
665 : : !! evaluation of phonons in ppddb9, it has been
666 : : !! optimized carefully. There is still possibility
667 : : !! for improvement, by using bloking on G and R!
668 : : !! 3. There can be small numerical variations due to the
669 : : !! fact that the input dielectric tensor is usually
670 : : !! not perfectly symmetric.
671 : : !!
672 : : !! SOURCE
673 : :
674 : 176318 : subroutine ewald9(dielt,dyew,gmet,gprimd,natom,qphon,rmet,rprimd,sumg0,ucvol,xred,zeff, qdrp_cart, eta, &
675 : : option, dipquad, quadquad) ! optional
676 : :
677 : : !Arguments -------------------------------
678 : : !scalars
679 : : integer,intent(in) :: natom,sumg0
680 : : integer,optional,intent(in) :: option, dipquad, quadquad
681 : : real(dp),intent(in) :: ucvol
682 : : !arrays
683 : : real(dp),intent(in) :: dielt(3,3),gmet(3,3),gprimd(3,3),qphon(3)
684 : : real(dp),intent(in) :: rmet(3,3),rprimd(3,3),xred(3,natom),zeff(3,3,natom)
685 : : real(dp),intent(in) :: qdrp_cart(3,3,3,natom)
686 : : real(dp),intent(out) :: dyew(2,3,natom,3,natom), eta
687 : :
688 : : !Local variables -------------------------
689 : : !scalars
690 : : integer,parameter :: mr=10000
691 : : integer :: ia,ib,ig1,ig2,ig3,ii,ll,kk,ir,ir1,ir2,ir3,jj
692 : : integer :: info,lwork,mu,newg,newr,ng,nr,nu,ng_expxq
693 : : integer :: ewald_option
694 : : integer :: dipquad_,quadquad_
695 : : logical :: do_quadrupole
696 : : logical, save :: firstcall = .TRUE.
697 : : real(dp),parameter :: fac=4.0_dp/3.0_dp/sqrt(pi)
698 : : real(dp),parameter :: fact2=2.0_dp/sqrt(pi)
699 : : real(dp),parameter :: y2max=64.0_dp, y2min=1.0d-24
700 : : real(dp) :: cddi,cddr,cqdi,cqdr,cqqi,cqqr,g3,g4
701 : : real(dp) :: arg1,arg2,arg3,arga,c123r,c123i,c23i,c23r,detdlt,inv_detdlt
702 : : real(dp) :: direct,fact1,fact3,gsq,recip,reta,reta3,inv4eta
703 : : real(dp) :: minexparg,sigma_max
704 : : real(dp) :: term1,term2,term3,term4,term5,y2,yy,invy,invy2,derfc_yy
705 : : character(len=700) :: msg
706 : : !arrays
707 : : real(dp) :: c1i(2*mr+1),c1r(2*mr+1),c2i(2*mr+1),c2r(2*mr+1),c3i(2*mr+1)
708 : 352636 : real(dp) :: c3r(2*mr+1),cosqxred(natom),wdielt(3,3),eig_dielt(3),gpq(3),gpqfac(3,3),gpqgpq(3,3)
709 : 352636 : real(dp) :: invdlt(3,3),ircar(3),ircax(3),rr(3),sinqxred(natom)
710 : 352636 : real(dp) :: xredcar(3,natom),xredcax(3,natom),xredicar(3),xredicax(3),xx(3)
711 : 176318 : real(dp),allocatable :: dyddt(:,:,:,:,:), dydqt(:,:,:,:,:,:), dyqqt(:,:,:,:,:,:,:)
712 : 176318 : real(dp),allocatable :: work(:)
713 [ + + ]: 176318 : complex(dp) :: exp2piqx(natom)
714 : 176318 : complex(dp),allocatable :: expx1(:,:), expx2(:,:), expx3(:,:)
715 : : ! *********************************************************************
716 : :
717 : : ! This routine is expensive so skip the calculation and return zeros if zeff == zero.
718 : : ! Typically this happens when the DDB file does not contains zeff but dipdip = 1 is used (default).
719 [ + + + + : 16413364 : if (all(zeff == zero).and.all(qdrp_cart == zero)) then
+ + + + +
+ + + + +
+ + + + +
+ ]
720 [ + + + + : 3128282 : dyew = zero; return
+ + + + +
+ ]
721 : : end if
722 [ + + + + : 13099088 : do_quadrupole = any(qdrp_cart /= zero)
+ + + + +
+ ]
723 : :
724 : : ! Keep track of total time spent.
725 : : !call timab(1749, 1, tsec)
726 : :
727 : : ! Initialize dipquad and quadquad options
728 [ + + ]: 136408 : dipquad_=0; if(present(dipquad)) dipquad_=dipquad
729 [ + + ]: 136408 : quadquad_=0; if(present(quadquad)) quadquad_=quadquad
730 : :
731 : : ! Deactivate real space sums for quadrupolar fields or for dipdip = -1
732 [ + + ]: 136408 : ewald_option = 0; if (present(option)) ewald_option = option
733 [ + + + + ]: 136408 : if (do_quadrupole.and.(dipquad_==1.or.quadquad_==1)) ewald_option = 1
734 : : !ewald_option = 0
735 : :
736 : : !This is the minimum argument of an exponential, with some safety
737 : 136264 : minexparg=log(tiny(0._dp))+five
738 [ + + ]: 136408 : if (ewald_option == 1) minexparg=-20.0_dp
739 : :
740 : : ! initialize complex phase factors
741 [ + + ]: 475400 : do ia = 1, natom
742 : : arga = two_pi*( (qphon(1))*xred(1,ia)&
743 : : +(qphon(2))*xred(2,ia)&
744 : 338992 : +(qphon(3))*xred(3,ia) )
745 : 475400 : exp2piqx(ia) = exp(arga*j_dpc)
746 : : end do
747 : 136408 : ng_expxq = 1000
748 [ - + + - : 409224 : ABI_MALLOC(expx1, (-ng_expxq:ng_expxq, natom))
+ - + - -
+ - + ]
749 [ - + + - : 272816 : ABI_MALLOC(expx2, (-ng_expxq:ng_expxq, natom))
- + - + ]
750 [ - + + - : 272816 : ABI_MALLOC(expx3, (-ng_expxq:ng_expxq, natom))
- + - + ]
751 [ + + ]: 475400 : do ia = 1, natom
752 [ + + ]: 678798392 : do ig1 = -ng_expxq, ng_expxq
753 : 678322992 : expx1(ig1, ia) = exp(ig1*two_pi*xred(1,ia)*j_dpc)
754 : 678322992 : expx2(ig1, ia) = exp(ig1*two_pi*xred(2,ia)*j_dpc)
755 : 678661984 : expx3(ig1, ia) = exp(ig1*two_pi*xred(3,ia)*j_dpc)
756 : : end do
757 : : end do
758 : :
759 : : ! compute eta for approximately optimized summations:
760 : : direct=rmet(1,1)+rmet(1,2)+rmet(1,3)+rmet(2,1)+&
761 : 136408 : rmet(2,2)+rmet(2,3)+rmet(3,1)+rmet(3,2)+rmet(3,3)
762 : : recip=gmet(1,1)+gmet(1,2)+gmet(1,3)+gmet(2,1)+&
763 : 136408 : gmet(2,2)+gmet(2,3)+gmet(3,1)+gmet(3,2)+gmet(3,3)
764 : :
765 : 136408 : eta=pi*100.0_dp/33.0_dp*sqrt(1.69_dp*recip/direct)
766 : :
767 : : ! Compute a material-dependent width for the Gaussians that hopefully
768 : : ! will make the Ewald real-space summation unnecessary.
769 [ + + ]: 136408 : if (ewald_option == 1) then
770 : :
771 : 145 : wdielt(:,:)=dielt(:,:)
772 : :
773 : : ! Diagonalize dielectric matrix
774 : 145 : lwork=-1
775 [ - + ]: 145 : ABI_MALLOC(work,(10))
776 : 145 : call dsyev('N','U',3, wdielt, 3, eig_dielt, work, lwork,info)
777 : 145 : lwork=nint(work(1))
778 : 145 : ABI_FREE(work)
779 : :
780 [ - + + - : 435 : ABI_MALLOC(work,(lwork))
+ - + - -
+ - + ]
781 : 145 : call dsyev('V','U',3, wdielt, 3, eig_dielt, work, lwork,info)
782 [ - + ]: 145 : ABI_FREE(work)
783 : :
784 : : ! This is a tentative maximum value for the gaussian width in real space
785 : 145 : sigma_max=three
786 : :
787 : : ! Set eta taking into account that the eps_inf is used as a metric in reciprocal space
788 [ - + + - : 725 : eta=sqrt(maxval(eig_dielt))/sigma_max
+ + ]
789 : :
790 [ + + ]: 145 : if (firstcall) then
791 : 3 : firstcall = .FALSE.
792 : 3 : write(msg, '(4a,f9.4,9a)' ) ch10,&
793 : 3 : ' Warning : due to the use of quadrupolar fields, the width of the reciprocal space gaussians', ch10, &
794 : 3 : ' in ewald9 has been set to eta= ', eta, ' 1/bohr and the real-space sums have been neglected.', ch10, &
795 : 3 : ' One should check whether this choice leads to correct results for the specific system under study', &
796 : 3 : ' and q-point grid.',ch10, &
797 : 3 : ' It is recommended to check that calculations with dipdip=1 and -1 (both with dipquad=0 and quadquad=0)', ch10, &
798 : 6 : ' lead to identical results. Otherwise increase the resolution of the q-point grid and repeat this test.', ch10
799 [ + + ]: 9 : call wrtout([ab_out,std_out], msg)
800 : : end if
801 : :
802 : : !Internally eta is the square of the gaussians width
803 : 145 : eta=eta*eta
804 : : end if
805 : :
806 : 136408 : inv4eta = one / four / eta
807 : :
808 [ - + + - : 682040 : ABI_MALLOC(dyddt,(2,3,natom,3,natom))
+ - - + +
- + - + -
- + - + ]
809 [ - + + - : 682040 : ABI_MALLOC(dydqt,(2,3,natom,3,natom,3))
- + + - +
- + - + -
- + - + ]
810 [ - + + - : 682040 : ABI_MALLOC(dyqqt,(2,3,natom,3,natom,3,3))
- + + - +
- + - + -
- + - + ]
811 : :
812 [ + + + + : 36502136 : dyddt = zero
+ + + + +
+ ]
813 [ + + + + : 109642816 : dydqt = zero
+ + + + +
+ + + ]
814 [ + + + + : 329064856 : dyqqt = zero
+ + + + +
+ + + +
+ ]
815 : :
816 : : !Sum terms over g space:
817 : : ng=0
818 : : do
819 : 1842318 : ng=ng+1
820 : :
821 : : ! if needed, update the complex phases for larger G vectors
822 [ - + ]: 1842318 : if (ng > ng_expxq) then
823 : : !write(std_out,*)"have to realloc"
824 [ # # ]: 0 : ABI_FREE(expx1)
825 [ # # ]: 0 : ABI_FREE(expx2)
826 [ # # ]: 0 : ABI_FREE(expx3)
827 : :
828 : 0 : ng_expxq = ng_expxq*2
829 : : ! TODO: half of this space is not needed, as it contains the complex conjugate of the other half.
830 : : ! present duplication avoids if statements inside the loop, however
831 [ # # # # : 0 : ABI_MALLOC(expx1, (-ng_expxq:ng_expxq, natom))
# # # # #
# # # # #
# # ]
832 [ # # # # : 0 : ABI_MALLOC(expx2, (-ng_expxq:ng_expxq, natom))
# # # # #
# # # ]
833 [ # # # # : 0 : ABI_MALLOC(expx3, (-ng_expxq:ng_expxq, natom))
# # # # #
# # # ]
834 [ # # ]: 0 : do ia = 1, natom
835 [ # # ]: 0 : do ig1 = -ng_expxq, ng_expxq
836 : 0 : expx1(ig1, ia) = exp(ig1*two_pi*xred(1,ia)*j_dpc)
837 : 0 : expx2(ig1, ia) = exp(ig1*two_pi*xred(2,ia)*j_dpc)
838 : 0 : expx3(ig1, ia) = exp(ig1*two_pi*xred(3,ia)*j_dpc)
839 : : end do
840 : : end do
841 : : end if
842 : :
843 : 1842318 : newg=0
844 [ + + ]: 33141558 : do ig3=-ng,ng
845 [ + + ]: 786789988 : do ig2=-ng,ng
846 [ + + ]: 26764192702 : do ig1=-ng,ng
847 [ + + + + : 26732893462 : if(abs(ig1)==ng .or. abs(ig2)==ng .or. abs(ig3)==ng .or. ng==1 )then
+ + + + ]
848 : :
849 : : gpq(1)=(ig1+qphon(1))*gprimd(1,1)+(ig2+qphon(2))*&
850 : 4161174652 : gprimd(1,2)+(ig3+qphon(3))*gprimd(1,3)
851 : : gpq(2)=(ig1+qphon(1))*gprimd(2,1)+(ig2+qphon(2))*&
852 : 4161174652 : gprimd(2,2)+(ig3+qphon(3))*gprimd(2,3)
853 : : gpq(3)=(ig1+qphon(1))*gprimd(3,1)+(ig2+qphon(2))*&
854 : 4161174652 : gprimd(3,2)+(ig3+qphon(3))*gprimd(3,3)
855 : 4161174652 : gsq=zero
856 [ + + ]: 16644698608 : do jj=1,3
857 [ + + ]: 54095270476 : do ii=1,3
858 : 37450571868 : gpqgpq(ii,jj)=gpq(ii)*gpq(jj)
859 : 49934095824 : gsq=gsq+gpqgpq(ii,jj)*dielt(ii,jj)
860 : : end do
861 : : end do
862 : :
863 : : ! Skip q=0:
864 [ + + ]: 4161174652 : if (gsq<1.0d-20) then
865 [ - + ]: 19501 : if (sumg0==1) then
866 : : write(msg,'(5a)' )&
867 : 0 : 'The phonon wavelength should not be zero :',ch10,&
868 : 0 : 'there are non-analytical terms that cannot be treated.',ch10,&
869 : 0 : 'Action: subtract this wavelength from the input file.'
870 : 0 : ABI_ERROR(msg)
871 : : end if
872 : :
873 : : else
874 : :
875 : 4161155151 : arg1=(two_pi**2)*gsq* inv4eta
876 : :
877 : : ! Larger arg gives 0 contribution:
878 [ + + ]: 4161155151 : if (arg1<= -minexparg ) then
879 : 1095877376 : newg=1
880 : :
881 : : ! Here calculate the term
882 : 1095877376 : term1=exp(-arg1)/gsq
883 [ + + ]: 4383509504 : do jj=1,3
884 [ + + ]: 14246405888 : do ii=1,3
885 : 13150528512 : gpqfac(ii,jj)=gpqgpq(ii,jj)*term1
886 : : end do
887 : : end do
888 : :
889 : : ! MJV: replaced old calls to cos and sin.
890 : : ! Checked for 10 tests in v2 that max error is about 6.e-15, usually < 2.e-15
891 [ + + ]: 4519989265 : do ia=1,natom
892 : 3424111889 : cosqxred(ia)= real(exp2piqx(ia)*expx1(ig1, ia)*expx2(ig2, ia)*expx3(ig3, ia))
893 : 4519989265 : sinqxred(ia)=aimag(exp2piqx(ia)*expx1(ig1, ia)*expx2(ig2, ia)*expx3(ig3, ia))
894 : : end do
895 : :
896 : : ! First, the diagonal terms
897 [ + + ]: 4383509504 : do nu=1,3
898 [ + + ]: 14655845171 : do ia=1,natom
899 [ + + ]: 34104639129 : do mu=nu,3
900 : 30817007001 : dyddt(1,mu,ia,nu,ia)=dyddt(1,mu,ia,nu,ia)+gpqfac(mu,nu)
901 : : end do
902 : : end do
903 : : end do
904 : :
905 : : ! Then, the non-diagonal ones
906 [ + + ]: 3424111889 : do ib=2,natom
907 [ + + ]: 9851475931 : do ia=1,ib-1
908 : : ! phase factor dipole-dipole
909 : 6427364042 : cddr=cosqxred(ia)*cosqxred(ib)+sinqxred(ia)*sinqxred(ib)
910 : 6427364042 : cddi=sinqxred(ia)*cosqxred(ib)-cosqxred(ia)*sinqxred(ib)
911 : :
912 : : ! Dipole-dipole contribution
913 [ + + ]: 28037690681 : do nu=1,3
914 [ + + ]: 64273640420 : do mu=nu,3
915 : 38564184252 : dyddt(1,mu,ia,nu,ib)=dyddt(1,mu,ia,nu,ib)+gpqfac(mu,nu)*cddr
916 : 57846276378 : dyddt(2,mu,ia,nu,ib)=dyddt(2,mu,ia,nu,ib)+gpqfac(mu,nu)*cddi
917 : : end do
918 : : end do
919 : : end do
920 : : end do
921 : :
922 [ + + ]: 1095877376 : if (do_quadrupole) then
923 [ + + ]: 341240163 : do ib=1,natom
924 [ + + ]: 803577258 : do ia=1,natom
925 : :
926 : : ! phase factor for dipole-quadrupole
927 : 462337095 : cqdr=cosqxred(ia)*sinqxred(ib)-sinqxred(ia)*cosqxred(ib)
928 : 462337095 : cqdi=cosqxred(ia)*cosqxred(ib)+sinqxred(ia)*sinqxred(ib)
929 : :
930 : : ! phase factor quadrupole-quadrupole
931 : 462337095 : cqqr=cosqxred(ia)*cosqxred(ib)+sinqxred(ia)*sinqxred(ib)
932 : 462337095 : cqqi=sinqxred(ia)*cosqxred(ib)-cosqxred(ia)*sinqxred(ib)
933 : :
934 : : ! Dipole-quadrupole contribution
935 [ + + ]: 1849348380 : do ii=1,3
936 [ + + ]: 6010382235 : do jj=1,3
937 [ + + ]: 18031146705 : do kk=1,3
938 : 12483101565 : g3=gpq(ii)*gpq(jj)*gpq(kk)
939 : 12483101565 : dydqt(1,ii,ia,jj,ib,kk)=dydqt(1,ii,ia,jj,ib,kk)+g3*term1*cqdr
940 : 16644135420 : dydqt(2,ii,ia,jj,ib,kk)=dydqt(2,ii,ia,jj,ib,kk)+g3*term1*cqdi
941 : : end do ! kk
942 : : end do ! jj
943 : : end do ! ii
944 : :
945 : : ! Quadrupole-quadrupole contribution
946 [ + + ]: 2077510023 : do ii=1,3
947 [ + + ]: 6010382235 : do jj=1,3
948 [ + + ]: 18031146705 : do kk=1,3
949 [ + + ]: 54093440115 : do ll=1,3
950 : 37449304695 : g4 = gpq(ii)*gpq(jj)*gpq(kk)*gpq(ll)
951 : 37449304695 : dyqqt(1,ii,ia,jj,ib,kk,ll)=dyqqt(1,ii,ia,jj,ib,kk,ll)+g4*term1*cqqr
952 : 49932406260 : dyqqt(2,ii,ia,jj,ib,kk,ll)=dyqqt(2,ii,ia,jj,ib,kk,ll)+g4*term1*cqqi
953 : : end do
954 : : end do ! kk
955 : : end do ! jj
956 : : end do ! ii
957 : : end do ! ia
958 : : end do ! ib
959 : : end if
960 : :
961 : : end if ! endif exp() argument is smaller than -minexparg
962 : : end if ! Endif g/=0 :
963 : : end if ! End triple summation over Gs:
964 : : end do
965 : : end do
966 : : end do
967 : :
968 : : ! Check if new shell must be calculated
969 [ + + ]: 1842318 : if(newg==0)exit
970 : : end do
971 : :
972 : : ! Multiplies by common factor
973 : 136408 : fact1=4.0_dp*pi/ucvol
974 [ + + ]: 475400 : do ib=1,natom
975 [ + + ]: 1228392 : do ia=1,ib
976 [ + + ]: 3350960 : do nu=1,3
977 [ + + ]: 7529920 : do mu=nu,3
978 : 4517952 : dyddt(1,mu,ia,nu,ib)=dyddt(1,mu,ia,nu,ib)*fact1
979 : 6776928 : dyddt(2,mu,ia,nu,ib)=dyddt(2,mu,ia,nu,ib)*fact1
980 : : end do
981 : : end do
982 : : end do
983 : : end do
984 [ + + ]: 136408 : if (do_quadrupole) then
985 [ + + + + : 3127806 : dydqt=dydqt*fact1/two * two_pi
+ + + + +
+ + + ]
986 [ + + + + : 9390447 : dyqqt=dyqqt*fact1/four * two_pi ** 2
+ + + + +
+ + + +
+ ]
987 : : end if
988 : :
989 : 136408 : reta=sqrt(eta)
990 : 136408 : reta3=-eta*reta
991 : :
992 : : !Calculating the inverse (transpose) of the dielectric tensor
993 : 136408 : call matr3inv(dielt,invdlt)
994 : :
995 : : !Calculating the determinant of the dielectric tensor
996 : : detdlt=dielt(1,1)*dielt(2,2)*dielt(3,3)+dielt(1,3)*dielt(2,1)*&
997 : : dielt(3,2)+dielt(1,2)*dielt(2,3)*dielt(3,1)-dielt(1,3)*&
998 : : dielt(2,2)*dielt(3,1)-dielt(1,1)*dielt(2,3)*dielt(3,2)-&
999 : 136408 : dielt(1,2)*dielt(2,1)*dielt(3,3)
1000 : :
1001 [ - + ]: 136408 : if(detdlt<tol6)then
1002 : : write(msg, '(a,es16.6,11a)' )&
1003 : 0 : 'The determinant of the dielectrix matrix, detdlt=',detdlt,' is smaller than 1.0d-6.',ch10,&
1004 : 0 : 'The use of the dipole-dipole model for interatomic force constants is not possible.',ch10,&
1005 : 0 : 'It is likely that you have not treated the electric field perturbations,',ch10,&
1006 : 0 : 'because you not are dealing with an insulator, so that',ch10,&
1007 : 0 : 'your dielectric matrix was simply set to zero in the Derivative DataBase.',ch10,&
1008 : 0 : 'Action: set the input variable dipdip to 0 .'
1009 : 0 : ABI_ERROR(msg)
1010 : : end if
1011 : :
1012 : 136408 : inv_detdlt = one / sqrt(detdlt)
1013 : 136408 : fact3=reta3 * inv_detdlt
1014 : :
1015 [ + + ]: 136408 : if (ewald_option /= 1) then
1016 : : ! Preparing the loop on real space
1017 [ + + ]: 474965 : do ia=1,natom
1018 [ + + ]: 1491071 : do ii=1,3
1019 : : xredcar(ii,ia)=(xred(1,ia)*rprimd(ii,1)+&
1020 : : xred(2,ia)*rprimd(ii,2)+&
1021 : 1354808 : xred(3,ia)*rprimd(ii,3) )*reta
1022 : : end do
1023 : : end do
1024 [ + + ]: 474965 : do ia=1,natom
1025 [ + + ]: 1491071 : do ii=1,3
1026 : : xredcax(ii,ia)= invdlt(1,ii)*xredcar(ii,ia)+&
1027 : : invdlt(2,ii)*xredcar(ii,ia)+&
1028 : 1354808 : invdlt(3,ii)*xredcar(ii,ia)
1029 : : end do
1030 : : end do
1031 : :
1032 : : ! Prepare the evaluation of exp(iq*R)
1033 [ + + ]: 2725532526 : do ir=-mr,mr
1034 : 2725396263 : arg1=-two_pi*qphon(1)*ir
1035 : 2725396263 : arg2=-two_pi*qphon(2)*ir
1036 : 2725396263 : arg3=-two_pi*qphon(3)*ir
1037 : 2725396263 : c1r(ir+mr+1)=cos(arg1)
1038 : 2725396263 : c1i(ir+mr+1)=sin(arg1)
1039 : 2725396263 : c2r(ir+mr+1)=cos(arg2)
1040 : 2725396263 : c2i(ir+mr+1)=sin(arg2)
1041 : 2725396263 : c3r(ir+mr+1)=cos(arg3)
1042 : 2725532526 : c3i(ir+mr+1)=sin(arg3)
1043 : : end do
1044 : :
1045 [ + - ]: 1174417 : do nr=1,mr
1046 : 1174417 : newr=0
1047 : :
1048 : : ! Begin big loop on real space vectors
1049 [ + + ]: 14990442 : do ir3=-nr,nr
1050 [ + + ]: 235265067 : do ir2=-nr,nr
1051 : :
1052 : : ! Here, construct the cosine and sine of q*R for components 2 and 3
1053 : 220274625 : c23r = c2r(ir2+mr+1) * c3r(ir3+mr+1) - c2i(ir2+mr+1) * c3i(ir3+mr+1)
1054 : 220274625 : c23i = c2i(ir2+mr+1) * c3r(ir3+mr+1) + c2r(ir2+mr+1) * c3i(ir3+mr+1)
1055 : :
1056 : : ! Also multiplies by fact3, because it is a rather economical place to do so
1057 : 220274625 : c23r=c23r * fact3
1058 : 220274625 : c23i=c23i * fact3
1059 : :
1060 [ + + ]: 4567160227 : do ir1=-nr,nr
1061 [ + + + + : 4553344202 : if( abs(ir3)==nr .or. abs(ir2)==nr .or. abs(ir1)==nr .or. nr==1 )then
+ + + + ]
1062 : :
1063 : : ! This is the real part and imaginary part of the phase factor exp(iq*R)
1064 : 1165387049 : c123r = c1r(ir1+mr+1) * c23r - c1i(ir1+mr+1) * c23i
1065 : 1165387049 : c123i = c1i(ir1+mr+1) * c23r + c1r(ir1+mr+1) * c23i
1066 : :
1067 [ + + ]: 4661548196 : do ii=1,3
1068 : : ircar(ii)= ( ir1*rprimd(ii,1)+&
1069 : : ir2*rprimd(ii,2)+&
1070 : 4661548196 : ir3*rprimd(ii,3) ) * reta
1071 : : end do
1072 [ + + ]: 4661548196 : do ii=1,3
1073 : : ircax(ii)= invdlt(1,ii)*ircar(ii)+&
1074 : : invdlt(2,ii)*ircar(ii)+&
1075 : 4661548196 : invdlt(3,ii)*ircar(ii)
1076 : : end do
1077 : :
1078 : : ! Here loops on atoms
1079 [ + + ]: 3993426279 : do ib=1,natom
1080 [ + + ]: 11312156920 : do ii=1,3
1081 : 8484117690 : xredicar(ii)=ircar(ii)-xredcar(ii,ib)
1082 : 11312156920 : xredicax(ii)=ircax(ii)-xredcax(ii,ib)
1083 : : end do
1084 [ + + ]: 12675805186 : do ia=1,ib
1085 [ + + ]: 22058785516 : do ii=1,3
1086 : 16544089137 : rr(ii)=xredicar(ii)+xredcar(ii,ia)
1087 : 22058785516 : xx(ii)=xredicax(ii)+xredcax(ii,ia)
1088 : : end do
1089 : :
1090 : 5514696379 : y2=rr(1)*xx(1)+rr(2)*xx(2)+rr(3)*xx(3)
1091 : :
1092 : : ! The atoms should not be too far of each other
1093 [ + + ]: 8342735609 : if (y2 < y2max) then
1094 : : ! Note: erfc(8) is about 1.1e-29, so dont bother with larger y.
1095 : : ! Also: exp(-64) is about 1.6e-28, do dont bother with larger y**2 in exp.
1096 : :
1097 : : ! Avoid zero denominators in term:
1098 [ + + ]: 1268356686 : if (y2 >= y2min) then
1099 : 1268017984 : newr=1
1100 : 1268017984 : yy=sqrt(y2)
1101 : 1268017984 : invy=1.0_dp/yy
1102 : 1268017984 : invy2=invy**2
1103 : 1268017984 : derfc_yy = abi_derfc(yy)
1104 : 1268017984 : term2=derfc_yy*invy*invy2
1105 : 1268017984 : term3=fact2*exp(-y2)*invy2
1106 : 1268017984 : term4=-(term2+term3)
1107 : 1268017984 : term5=(3.0_dp*term2+term3*(3.0_dp+2.0_dp*y2))*invy2
1108 [ + + ]: 5072071936 : do nu=1,3
1109 [ + + ]: 12680179840 : do mu=nu,3
1110 : 7608107904 : dyddt(1,mu,ia,nu,ib)=dyddt(1,mu,ia,nu,ib)+c123r*(xx(nu)*xx(mu)*term5+term4*invdlt(nu,mu))
1111 : 11412161856 : dyddt(2,mu,ia,nu,ib)=dyddt(2,mu,ia,nu,ib)+c123i*(xx(nu)*xx(mu)*term5+term4*invdlt(nu,mu))
1112 : : end do
1113 : : end do
1114 : : else
1115 : : ! If zero denominator, the atoms should be identical
1116 [ - + ]: 338702 : if (ia/=ib)then
1117 : : write(msg, '(5a,i0,a,i0,a)' )&
1118 : 0 : 'The distance between two atoms seem to vanish.',ch10,&
1119 : 0 : 'This is not allowed.',ch10,&
1120 : 0 : 'Action: check the input for the atoms number',ia,' and',ib,'.'
1121 : 0 : ABI_ERROR(msg)
1122 : : else
1123 : : ! This is the correction when the atoms are identical
1124 [ + + ]: 1354808 : do nu=1,3
1125 [ + + ]: 4403126 : do mu=1,3
1126 : : dyddt(1,mu,ia,nu,ib)=dyddt(1,mu,ia,nu,ib)+&
1127 : 4064424 : fac*reta3*invdlt(nu,mu) * inv_detdlt
1128 : : end do
1129 : : end do
1130 : : end if
1131 : : end if ! End the condition for avoiding zero denominators
1132 : : end if ! End the condition of too large distance between atoms
1133 : : end do
1134 : : end do ! End loop over ia and ib :
1135 : : end if ! End triple loop over real space points:
1136 : : end do ! ir1
1137 : : end do ! ir2
1138 : : end do ! ir3
1139 : :
1140 : : ! Check if new shell must be calculated
1141 [ + + ]: 1174417 : if(newr==0)exit
1142 [ - + ]: 1174417 : if(newr==1 .and. nr==mr) ABI_BUG('mr is too small')
1143 : : end do
1144 : : end if ! check if should compute real part
1145 : :
1146 : : ! Now, symmetrizes
1147 [ + + ]: 338992 : do ib=1,natom-1
1148 [ + + ]: 946744 : do nu=1,3
1149 [ + + ]: 2052336 : do ia=ib+1,natom
1150 [ + + ]: 4333752 : do mu=nu,3
1151 : 2484000 : dyddt(1,mu,ia,nu,ib)= dyddt(1,mu,ib,nu,ia)
1152 : 3726000 : dyddt(2,mu,ia,nu,ib)=-dyddt(2,mu,ib,nu,ia)
1153 : : end do
1154 : : end do
1155 : : end do
1156 : : end do
1157 : :
1158 [ + + ]: 475400 : do ib=1,natom
1159 [ + + ]: 1153384 : do nu=2,3
1160 [ + + ]: 3350960 : do ia=1,natom
1161 [ + + ]: 6512944 : do mu=1,nu-1
1162 : 3500976 : dyddt(1,mu,ia,nu,ib)=dyddt(1,nu,ia,mu,ib)
1163 : 5834960 : dyddt(2,mu,ia,nu,ib)=dyddt(2,nu,ia,mu,ib)
1164 : : end do
1165 : : end do
1166 : : end do
1167 : : end do
1168 : :
1169 : : !Tests
1170 : : !write(std_out,*)' ewald9 : take into account the effective charges '
1171 [ + + + + : 36502136 : dyew = zero
+ + + + +
+ ]
1172 [ + + ]: 475400 : do ib=1,natom
1173 [ + + ]: 1492376 : do nu=1,3
1174 [ + + ]: 4856944 : do ia=1,natom
1175 [ + + ]: 15020880 : do mu=1,3
1176 [ + + ]: 45512688 : do ii=1,3
1177 [ + + ]: 136538064 : do jj=1,3
1178 : : ! dipole-dipole correction
1179 : : dyew(1,mu,ia,nu,ib)=dyew(1,mu,ia,nu,ib) + &
1180 : 94526352 : zeff(ii,mu,ia)*zeff(jj,nu,ib)*dyddt(1,ii,ia,jj,ib)
1181 : : dyew(2,mu,ia,nu,ib)=dyew(2,mu,ia,nu,ib) + &
1182 : 94526352 : zeff(ii,mu,ia)*zeff(jj,nu,ib)*dyddt(2,ii,ia,jj,ib)
1183 [ + + ]: 126035136 : if (do_quadrupole) then
1184 [ + + ]: 10514124 : do kk=1,3
1185 [ + + ]: 7885593 : if (dipquad_==1) then
1186 : : ! dipole-quadrupole correction
1187 : : dyew(1,mu,ia,nu,ib)=dyew(1,mu,ia,nu,ib) + &
1188 : : (zeff(ii,nu,ib)*qdrp_cart(kk,jj,mu,ia) - &
1189 : 139968 : zeff(ii,mu,ia)*qdrp_cart(kk,jj,nu,ib)) * dydqt(1,ii,ia,jj,ib,kk)
1190 : : dyew(2,mu,ia,nu,ib)=dyew(2,mu,ia,nu,ib) + &
1191 : : (zeff(ii,nu,ib)*qdrp_cart(kk,jj,mu,ia) - &
1192 : 139968 : zeff(ii,mu,ia)*qdrp_cart(kk,jj,nu,ib)) * dydqt(2,ii,ia,jj,ib,kk)
1193 : : end if
1194 : :
1195 : : ! quadrupole-quadrupole correction
1196 [ + + ]: 10514124 : if (quadquad_==1) then
1197 [ + + ]: 559872 : do ll=1,3
1198 : : dyew(1,mu,ia,nu,ib)=dyew(1,mu,ia,nu,ib) + &
1199 : 419904 : (qdrp_cart(ll,ii,mu,ia)*qdrp_cart(kk,jj,nu,ib)) * dyqqt(1,ii,ia,jj,ib,kk,ll)
1200 : : dyew(2,mu,ia,nu,ib)=dyew(2,mu,ia,nu,ib) + &
1201 : 559872 : (qdrp_cart(ll,ii,mu,ia)*qdrp_cart(kk,jj,nu,ib)) * dyqqt(2,ii,ia,jj,ib,kk,ll)
1202 : : end do
1203 : : end if
1204 : : end do
1205 : : end if
1206 : :
1207 : : end do
1208 : : end do
1209 : : end do
1210 : : end do
1211 : : end do
1212 : : end do
1213 : :
1214 [ - + ]: 136408 : ABI_FREE(expx1)
1215 [ - + ]: 136408 : ABI_FREE(expx2)
1216 [ - + ]: 136408 : ABI_FREE(expx3)
1217 [ - + ]: 136408 : ABI_FREE(dyddt)
1218 [ - + ]: 136408 : ABI_FREE(dydqt)
1219 [ - + ]: 136408 : ABI_FREE(dyqqt)
1220 : : !call timab(1749, 2, tsec)
1221 : :
1222 : : end subroutine ewald9
1223 : : !!***
1224 : :
1225 : : !!****f* m_ewald/ewald9_2D
1226 : : !!
1227 : : !! NAME
1228 : : !! ewald9_2D
1229 : : !!
1230 : : !! FUNCTION
1231 : : !! Compute the long-range electrostatics contribution to interatomic force constants
1232 : : !! in the bi-dimensional (2D) case, considering the 2D is embedded in a dielectric environment
1233 : : !! and has a given dielectric thickness. The singularity of the Coulomb potential is treated
1234 : : !! using the Ewald summation approach. It is possible to input a more complicate model
1235 : : !! with two consecutive dielectric slabs.
1236 : : !!
1237 : : !! INPUTS
1238 : : !! natom=number of atoms in unit cell
1239 : : !! rprimd(3,3)=dimensional primitive translations for real space (bohr)
1240 : : !! xred(3,natom)=reduced coordinates of the atoms
1241 : : !! dielt(3,3)=dielectric tensor of the 2D
1242 : : !! dyew(2,3,natom,3,natom)=long-range electrostatics IFCs following Ewald
1243 : : !! qphon(3)=phonon wavevector in reduced coordinates
1244 : : !! zeff(3,3,natom)=Born effective charge tensor
1245 : : !! qdrp_cart(3,3,3,natom)=Dynamical quadrupoles
1246 : : !! dielt_env=dielectric constant of the embedding environment (1 in vacuum)
1247 : : !! thick(2)=dielectric thicknesses of the slab, first value correspond to the outer
1248 : : !! dielectric, second to the inner dielectric slab (if any)
1249 : : !! sys_dim= system dimensionality (indicates axis without periodicity)
1250 : : !!
1251 : : !! OUTPUT
1252 : : !! dyew(2,3,natom,3,natom)=long-range electrostatics IFCs following Ewald
1253 : : !!
1254 : : !! SOURCE
1255 : :
1256 : 561 : subroutine ewald9_2D(natom,rprimd,xred,dielt,dyew,qphon,zeff,qdrp_cart,dielt_env,thick,sys_dim)
1257 : :
1258 : : !Arguments -------------------------------
1259 : : !scalars
1260 : : real(dp), intent(in) :: dielt_env
1261 : : integer :: natom, sys_dim
1262 : : !arrays
1263 : : real(dp),intent(in) :: rprimd(3,3),thick(2),xred(3,natom),dielt(3,3),qphon(3)
1264 : : real(dp),intent(in) :: zeff(3,3,natom),qdrp_cart(3,3,3,natom)
1265 : : real(dp),intent(out) :: dyew(2,3,natom,3,natom)
1266 : : character(len=700) :: msg
1267 : :
1268 : : !Local variables -------------------------
1269 : : !scalars
1270 : : integer :: gmax,idir1,idir2,ibz1,ibz2,ipert1,ipert2,inner_thick,ndir,mdir
1271 : : real(dp) :: detdlt, delta_perp, lambda, dielt_perp,dielt_perp1,dielt_perp2, eta,eta1,xi, dielt_eff
1272 : : real(dp) :: dielt_eff1,dielt_eff2, norm_kvec, norm_kvec0, phi
1273 : : real(dp) :: ewald_fun, ewald_fun1, ewald_fun2, rflct_coeff, out_thick
1274 : : real(dp) :: rflct_coeff1, rflct_coeff2, rprimd_perp
1275 : : real(dp) :: fac_erfc, fac_ewald1, fac_ewald2, fac_ewald2b,fac_exp
1276 : : real(dp) :: fac_exp1, trans_fun, fac_gauss, fac_mirror, fac_mirror1, fac_real
1277 : : real(dp) :: mean2_perp, mirror_diff, mirror_parapara, mirror_paraperp, mirror_perpperp
1278 : : real(dp) :: rvec_norm, sqrt_norm, ucsurf, xmean
1279 : : logical, save :: firstcall = .TRUE.
1280 : : !arrays
1281 : : integer :: periodic_dir(3)
1282 : 1122 : real(dp) :: dyew_real(2,3,natom,3,natom),dyew_rec(2,3,natom,3,natom),kvec(2),invdlt_para(2,2)
1283 : : real(dp) :: rprimd_para(2,2), gprimd_para(2,2), gvec(2)
1284 : : real(dp) :: norm_dielt(2), dielt_para(2,2), invdlt(3,3), qvec(3), qvec_para(2)
1285 : 1122 : real(dp) :: diff_xcart(3), xcart_para(2,natom),xcart_perp(natom)
1286 : 1122 : real(dp) :: kvec_para(2),zeff_para(2,3,natom), zeff_perp(3,natom)
1287 : 1122 : real(dp) :: xcart(3,natom), gprimd(3,3), rvec_dielt(3)
1288 : 1122 : real(dp) :: qdrp_parapara(2,2,3,natom),qdrp_perpperp(3,natom), qdrp_paraperp(2,3,natom)
1289 : : real(dp) :: rho_gerade1(2),rho_gerade2(2), rho_ungerade1(2), rho_ungerade2(2)
1290 : : ! *************************************************************************
1291 : :
1292 : 561 : periodic_dir(:) = 0
1293 [ - + ]: 561 : if (sys_dim ==2) then ! 2D along x
1294 : 0 : periodic_dir(2) =1 ; periodic_dir(3) = 1
1295 [ - + ]: 561 : elseif (sys_dim==3) then ! 2D along y
1296 : 0 : periodic_dir(1) = 1 ; periodic_dir(3) = 1
1297 [ + - ]: 561 : elseif (sys_dim==4) then ! 2D along z
1298 : 561 : periodic_dir(1) = 1 ; periodic_dir(2) = 1
1299 : : end if
1300 : :
1301 [ + + + + ]: 8017 : xcart=zero
1302 [ + + ]: 2425 : do ipert1=1,natom
1303 [ + + + + : 30385 : xcart(:,ipert1)=matmul(rprimd, xred(:,ipert1))
+ + ]
1304 : : end do
1305 : :
1306 : 561 : call matr3inv(rprimd,gprimd)
1307 [ + + + + : 26657 : zeff_para=zero ; zeff_perp=zero
+ + + + +
+ ]
1308 [ + + ]: 2244 : qvec(:) = qphon(1)*gprimd(:,1)+qphon(2)*gprimd(:,2)
1309 : : ndir=0
1310 [ + + ]: 2244 : do idir1=1,3
1311 [ + + ]: 2244 : if (periodic_dir(idir1)==1) then
1312 : 1122 : ndir=ndir+1
1313 : 1122 : qvec_para(ndir)=qvec(idir1)
1314 [ + + + + ]: 16034 : zeff_para(ndir,:,:) = zeff(idir1,:,:)
1315 [ + + ]: 4850 : xcart_para(ndir,:) = xcart(idir1,:)
1316 : : else
1317 [ + + + + ]: 8017 : zeff_perp(:,:) = zeff(idir1,:,:)
1318 [ + + ]: 2425 : xcart_perp(:) = xcart(idir1,:)
1319 [ + - ]: 561 : if (qvec(idir1)>tol6) then !Check if phonon mode is not out-of-plane
1320 : : write(msg, '(a,es16.6,5a)')&
1321 : 0 : 'The phonon wavevector along the confined direction is',qvec(idir1),' 1/Bohr >1.0d-6',ch10,&
1322 : 0 : 'The phonon wavevector should be purely along the periodic direction', ch10, &
1323 : 0 : 'when using Ewald summation in 2D. Please check your input file and structure'
1324 : 0 : ABI_ERROR(msg)
1325 : : end if
1326 : : end if
1327 : : end do
1328 [ + + ]: 2425 : xmean = sum(xcart_perp)/natom
1329 [ + + ]: 2425 : xcart_perp(:)=xcart_perp(:)-xmean
1330 : :
1331 : : !Calculating the inverse (transpose) of the dielectric tensor
1332 : 561 : call matr3inv(dielt,invdlt)
1333 : : !Calculating the determinant of the dielectric tensor
1334 : : detdlt=dielt(1,1)*dielt(2,2)*dielt(3,3)+dielt(1,3)*dielt(2,1)*&
1335 : : & dielt(3,2)+dielt(1,2)*dielt(2,3)*dielt(3,1)-dielt(1,3)*&
1336 : : & dielt(2,2)*dielt(3,1)-dielt(1,1)*dielt(2,3)*dielt(3,2)-&
1337 : 561 : & dielt(1,2)*dielt(2,1)*dielt(3,3)
1338 : :
1339 [ - + ]: 561 : if(detdlt<tol6)then
1340 : : write(msg, '(a,es16.6,11a)' )&
1341 : 0 : 'The determinant of the dielectrix matrix, detdlt=',detdlt,' is smaller than 1.0d-6.',ch10,&
1342 : 0 : 'The use of the dipole-dipole model for interatomic force constants is not possible.',ch10,&
1343 : 0 : 'It is likely that you have not treated the electric field perturbations,',ch10,&
1344 : 0 : 'because you not are dealing with an insulator, so that',ch10,&
1345 : 0 : 'your dielectric matrix was simply set to zero in the Derivative DataBase.',ch10,&
1346 : 0 : 'Action: set the input variable dipdip to 0 .'
1347 : 0 : ABI_ERROR(msg)
1348 : : end if
1349 : :
1350 : : ! The dielectric tensor must be diagonal in the confined direction
1351 : 561 : rprimd_para = zero ; rprimd_perp = zero
1352 : 561 : dielt_para = zero ; dielt_perp = zero
1353 : 561 : invdlt_para = zero ; gprimd_para = zero
1354 : 561 : ndir=0
1355 [ + + ]: 2244 : do idir1=1,3
1356 [ + + ]: 1683 : if (periodic_dir(idir1)==1) then
1357 : 1122 : ndir=ndir+1
1358 : : end if
1359 : 1683 : mdir=0
1360 [ + + ]: 7293 : do idir2=1,3
1361 [ + + ]: 5049 : if (periodic_dir(idir2)==1) then
1362 : 3366 : mdir=mdir+1
1363 : : end if
1364 [ + + + + : 5049 : if ((periodic_dir(idir1)==1 .and. periodic_dir(idir2)==0) .or. &
+ + + + ]
1365 : 1683 : (periodic_dir(idir1)==0 .and. periodic_dir(idir2)==1)) then
1366 [ + - - + ]: 2244 : if (abs(dielt(idir1,idir2))>tol6 .or. abs(rprimd(idir1,idir2))>tol6) then !No cross in-plane out-of-plane are allowed
1367 : : write(msg, '(7a)' )&
1368 : 0 : 'The dielectric matrix shows off-diagonal components in the confined direction larger than 1d-6',ch10,&
1369 : 0 : 'This is forbidden when considering the Ewald summation for 2D systems. Please check if your', ch10, &
1370 : 0 : 'confined direction is correctly specified in the anaddb input or if the vacuum size if sufficiently', ch10, &
1371 : 0 : 'large in the confined direction to avoid spurious interactions between unit cells'
1372 : 0 : ABI_ERROR(msg)
1373 : : end if
1374 [ + + + + ]: 32068 : qdrp_paraperp(ndir,:,:) = qdrp_cart(idir1,idir2,:,:)
1375 [ + + + - ]: 2805 : elseif (periodic_dir(idir1)==0 .and. periodic_dir(idir2)==0) then
1376 : 561 : dielt_perp = dielt(idir1,idir2)
1377 : 561 : rprimd_perp = rprimd(idir1,idir2)
1378 [ + + + + ]: 8017 : qdrp_perpperp(:,:) = qdrp_cart(idir1,idir2,:,:)
1379 : : else
1380 : 2244 : dielt_para(ndir,mdir) = dielt(idir1,idir2)
1381 : 2244 : rprimd_para(ndir,mdir)=rprimd(idir1,idir2)
1382 : 2244 : gprimd_para(ndir,mdir) = gprimd(idir1,idir2)
1383 : 2244 : invdlt_para(ndir,mdir)=invdlt(idir1,idir2)
1384 [ + + + + ]: 32068 : qdrp_parapara(ndir,mdir,:,:) = qdrp_cart(idir1,idir2,:,:)
1385 : : end if
1386 : : end do
1387 : : end do
1388 : :
1389 : :
1390 : : ! Consistency check for dielectric thicknesses
1391 : 561 : inner_thick = thick(2)
1392 : 561 : out_thick = thick(1)
1393 : :
1394 [ + + ]: 2425 : do ipert1=1,natom
1395 [ + + ]: 9121 : do ipert2=1,natom
1396 [ + - - + ]: 8560 : if (xcart_perp(ipert1)-xcart_perp(ipert2)>out_thick .or. xcart_perp(ipert1)+xcart_perp(ipert2)>out_thick) then
1397 : : write(msg, '(5a)' )&
1398 : 0 : 'Some atoms seem to be located away from the dielectric slab.',ch10,&
1399 : 0 : 'The present model only allows for this specific scenario. Please increases slightly',ch10,&
1400 : 0 : 'the dielectric slab thickness'
1401 : 0 : ABI_ERROR(msg)
1402 : : end if
1403 : : end do
1404 : : end do
1405 [ + + ]: 561 : if (inner_thick>zero) then
1406 [ - + ]: 190 : if (inner_thick>out_thick) then
1407 : : write(msg, '(3a)' )&
1408 : 0 : 'When considering consecutive dielectric slab models, inner thickness',ch10,&
1409 : 0 : 'should be smaller than the outer thickness. Please check your input file'
1410 : 0 : ABI_ERROR(msg)
1411 : : else
1412 [ + + ]: 570 : do ipert1=1,natom
1413 [ + + ]: 1330 : do ipert2=1,natom
1414 [ + - - + ]: 1140 : if (xcart_perp(ipert1)-xcart_perp(ipert2)>inner_thick .or. xcart_perp(ipert1)+xcart_perp(ipert2)>inner_thick) then
1415 : : write(msg, '(5a)' )&
1416 : 0 : 'Some atoms seem to be located away from the inner dielectric slab.',ch10,&
1417 : 0 : 'The present model only allows for this specific scenario. Please increases slightly',ch10,&
1418 : 0 : 'the inner dielectric slab thickness'
1419 : 0 : ABI_ERROR(msg)
1420 : : end if
1421 : : end do
1422 : : end do
1423 : : end if
1424 : : end if
1425 : :
1426 : : ! First needs to determine the Gaussian broadening intrinsic to the Ewald summation. In 2D, both the real and
1427 : : ! reciprocal summation are related to the complementary error function. We want to restrict the real-part to
1428 : : ! the first Wigner cell. We use the fact that sqrt(1-e^{-x^2}) < erf(x) < sqrt(1-e^{-4x^2/pi})
1429 : : ! and invert those relationships to estimate the broadening required to restrict the real-part summation of
1430 : : !the Ewald summation; here fixes the threshold to 1e-9 for contribution from later unit cells
1431 : : rvec_dielt = zero
1432 : 561 : ndir=0
1433 [ + + ]: 1683 : do idir1=1,2
1434 [ + + + + : 10659 : norm_dielt(idir1) = dot_product(rprimd_para(idir1,:),matmul(invdlt_para(:,:),rprimd_para(idir1,:)))
+ + ]
1435 : : end do
1436 : :
1437 [ - + + - : 2244 : lambda = dsqrt(maxval(norm_dielt))/dsqrt(-two*dlog(one-(one-tol9)**2))
+ + ]
1438 : :
1439 : : ! Now that we have computed the value of lambda, we can compute the Ewald summation, starting from the
1440 : : ! reciprocal sum. We need first to estimate the max. number of reciprocal vectors we need to consider
1441 : : ! for the Ewald summation. We use here a stricter tolerance than for the determination of lambda
1442 : : ! Note that the worst case scenario is always when considering Rka=Rk'b in this case
1443 : :
1444 [ - + + - : 2244 : gmax = int(dsqrt(-two*dlog(one-(one-tol12)**2))/lambda/(two_pi/dsqrt(minval(norm_dielt))))
+ + ]
1445 [ + + ]: 561 : if (firstcall) then
1446 : 2 : firstcall = .FALSE.
1447 : 2 : write(msg, '(6a,f9.4,3a,i3,1a)' ) ch10,&
1448 : 2 : ' Ewald treatment of 2D long-range electrostatics interatomic force constants', ch10, &
1449 : 2 : ' To restrict the real-part summation of Ewald to the first unit cell, the Gaussian broadening', ch10, &
1450 : 2 : ' has been set to ', lambda, ' 1/Bohr. For the reciprocal sum, this corresponds to max.', ch10, &
1451 : 4 : ' ',2*gmax-1, ' Brillouin zone repetitions in either in-plane directions'
1452 [ + + ]: 6 : call wrtout([ab_out,std_out], msg)
1453 : : end if
1454 : :
1455 [ + + + + : 208897 : dyew_rec = zero
+ + + + +
+ ]
1456 : : ! If one dielectric slab model, same dielectric for both regions
1457 : : ! Otherwise, inner dielectric ~1 and the other has been computed
1458 : : ! accordingly in the anaddb driver
1459 [ + + ]: 561 : if (inner_thick> zero) then
1460 : : dielt_perp1 = one
1461 : : else
1462 : 371 : dielt_perp1=dielt_perp
1463 : : end if
1464 : 561 : dielt_perp2 = dielt_perp
1465 [ + + ]: 8216 : do ibz1 = -gmax,gmax
1466 [ + + ]: 114681 : do ibz2 = -gmax,gmax
1467 [ + + ]: 319395 : gvec(:) = ibz1*gprimd_para(:,1)+ibz2*gprimd_para(:,2)
1468 [ + + ]: 319395 : kvec(:) = gvec(:) + qvec_para(:)
1469 [ + + ]: 319395 : kvec(:) = kvec(:)*two_pi
1470 [ + + + + ]: 745255 : kvec_para(:) = matmul(dielt_para,kvec)
1471 [ + + ]: 319395 : norm_kvec0 = dot_product(kvec,kvec)
1472 [ + + ]: 319395 : norm_kvec = dot_product(kvec,kvec_para)
1473 [ + + ]: 114120 : if (abs(norm_kvec)>tol6) then !Remove G=q=0 case
1474 : 106456 : eta = dsqrt(norm_kvec/dielt_perp)
1475 : 106456 : eta1 = dsqrt(norm_kvec/dielt_perp1)
1476 : 106456 : xi = dsqrt(norm_kvec/dielt_perp2)
1477 : : ! Effective dielectric constants (depends on direction)
1478 : 106456 : dielt_eff = dsqrt(norm_kvec*dielt_perp/norm_kvec0)
1479 : 106456 : dielt_eff1 = dsqrt(norm_kvec*dielt_perp1/norm_kvec0)
1480 : 106456 : dielt_eff2 = dsqrt(norm_kvec*dielt_perp2/norm_kvec0)
1481 : : ! Reflection coefficient at the dielectric interfaces
1482 : 106456 : rflct_coeff = (dielt_eff-dielt_env)/(dielt_eff+dielt_env)
1483 : 106456 : rflct_coeff2 = (dielt_eff2-dielt_env)/(dielt_eff2+dielt_env)
1484 : 106456 : trans_fun = (one+rflct_coeff2*dexp(-xi*(out_thick-inner_thick)))
1485 : 106456 : trans_fun = trans_fun/(one-rflct_coeff2*dexp(-xi*(out_thick-inner_thick)))
1486 : 106456 : rflct_coeff1 = (dielt_eff1*trans_fun-dielt_eff2)/(dielt_eff1*trans_fun+dielt_eff2)
1487 : : ! Dipole-dipole charges prefactors
1488 : 106456 : fac_exp = rflct_coeff*dexp(-eta*out_thick)
1489 : 106456 : fac_exp1 = rflct_coeff1*dexp(-eta1*inner_thick)
1490 : 106456 : fac_mirror = two*fac_exp/(one-fac_exp**2)
1491 : 106456 : fac_mirror1 = two*fac_exp1/(one-fac_exp1**2)
1492 : : ! Ewald factor in error function
1493 : 106456 : fac_ewald1= eta*dsqrt(dielt_perp)*lambda/dsqrt(two)
1494 [ + + ]: 486310 : do ipert1=1,natom
1495 [ + + ]: 1913786 : do ipert2=1,natom
1496 : 1427476 : delta_perp = (xcart_perp(ipert2)-xcart_perp(ipert1))
1497 : 1427476 : mean2_perp = (xcart_perp(ipert2)+xcart_perp(ipert1))
1498 : 1427476 : fac_ewald2= delta_perp/lambda/sqrt(two*dielt_perp)
1499 : 1427476 : fac_ewald2b= delta_perp/lambda/sqrt(two*dielt_perp1)
1500 : : ! Ewald function and derivatives (eta factorized)
1501 : : ewald_fun = half*(dexp(-eta*delta_perp)*(one-erf(fac_ewald1-fac_ewald2)))+ &
1502 : 1427476 : half*(dexp(eta*delta_perp)*(one-erf(fac_ewald1+fac_ewald2)))
1503 : : ewald_fun1 = half*(-dexp(-eta*delta_perp)*(one-erf(fac_ewald1-fac_ewald2)))+ &
1504 : 1427476 : half*(dexp(eta*delta_perp)*(one-erf(fac_ewald1+fac_ewald2)))
1505 : : ! For second derivative, there is in principle a Gaussian term as well
1506 : : ! However, by an appropriate choice of the electrostatic gauge (mean average
1507 : : ! potential), we can neglect it. This approximation has been validated
1508 : : ! with respect to real-space dipoles and exact calculated points
1509 : : ewald_fun2 = half*(dexp(-eta1*delta_perp)*(one-erf(fac_ewald1-fac_ewald2b)))+ &
1510 : 1427476 : half*(dexp(eta1*delta_perp)*(one-erf(fac_ewald1+fac_ewald2b)))
1511 : 1427476 : ewald_fun2=-ewald_fun2*eta1**2
1512 : :
1513 : : ! Phase factor and mirror terms
1514 [ + + ]: 4282428 : phi = dot_product(kvec,xcart_para(:,ipert1)-xcart_para(:,ipert2))
1515 : 1427476 : mirror_parapara = fac_mirror*(dcosh(eta*mean2_perp)+fac_exp*dcosh(eta*delta_perp))
1516 [ + + ]: 1427476 : if (inner_thick >tol6) then
1517 : 91940 : mirror_perpperp = fac_mirror1*(dcosh(eta1*mean2_perp)-fac_exp1*dcosh(eta1*delta_perp))
1518 : : else
1519 : 1335536 : mirror_perpperp = fac_mirror*(dcosh(eta*mean2_perp)-fac_exp1*dcosh(eta*delta_perp))
1520 : : end if
1521 : 1427476 : mirror_paraperp = fac_mirror*fac_exp*dsinh(eta*delta_perp)
1522 : 1427476 : mirror_diff = fac_mirror*dsinh(eta*mean2_perp)
1523 : : ! Then compute the charge prefactor
1524 [ + + ]: 6089758 : do idir1=1,3
1525 [ + + ]: 18557188 : do idir2=1,3
1526 [ + + + + : 115625556 : rho_gerade1(1) = -half*dot_product(kvec,matmul(qdrp_parapara(:,:,idir1,ipert1),kvec))
+ + ]
1527 [ + + + + : 115625556 : rho_gerade2(1) = -half*dot_product(kvec,matmul(qdrp_parapara(:,:,idir2,ipert2),kvec))
+ + ]
1528 : 12847284 : rho_gerade1(1) = rho_gerade1(1)+half*eta1**2*qdrp_perpperp(idir1,ipert1)
1529 : 12847284 : rho_gerade2(1) = rho_gerade2(1)+half*eta1**2*qdrp_perpperp(idir2,ipert2)
1530 [ + + ]: 38541852 : rho_gerade1(2) = -dot_product(kvec,zeff_para(:,idir1,ipert1))
1531 [ + + ]: 38541852 : rho_gerade2(2) = -dot_product(kvec,zeff_para(:,idir2,ipert2))
1532 : 12847284 : rho_ungerade1(1) = -zeff_perp(idir1,ipert1)
1533 [ + + ]: 38541852 : rho_ungerade1(2) = eta1*dot_product(kvec,qdrp_paraperp(:,idir1,ipert1))
1534 : 12847284 : rho_ungerade2(1) = -zeff_perp(idir2,ipert2)
1535 [ + + ]: 38541852 : rho_ungerade2(2) = eta1*dot_product(kvec,qdrp_paraperp(:,idir2,ipert2))
1536 : : ! First, add the source charges (gerade gerade)
1537 : : dyew_rec(1,idir1,ipert1,idir2,ipert2)= dyew_rec(1,idir1,ipert1,idir2,ipert2)+&
1538 : : (rho_gerade1(1)*rho_gerade2(1)+rho_gerade1(2)*rho_gerade2(2))&
1539 : 12847284 : *ewald_fun/eta/dielt_perp*cos(phi)
1540 : : dyew_rec(1,idir1,ipert1,idir2,ipert2)= dyew_rec(1,idir1,ipert1,idir2,ipert2)-&
1541 : : (rho_gerade1(2)*rho_gerade2(1)-rho_gerade1(1)*rho_gerade2(2))&
1542 : 12847284 : *ewald_fun/eta/dielt_perp*sin(phi)
1543 : : dyew_rec(2,idir1,ipert1,idir2,ipert2)= dyew_rec(2,idir1,ipert1,idir2,ipert2)+&
1544 : : (rho_gerade1(1)*rho_gerade2(1)+rho_gerade1(2)*rho_gerade2(2))&
1545 : 12847284 : *ewald_fun/eta/dielt_perp*sin(phi)
1546 : : dyew_rec(2,idir1,ipert1,idir2,ipert2)= dyew_rec(2,idir1,ipert1,idir2,ipert2)+&
1547 : : (rho_gerade1(2)*rho_gerade2(1)-rho_gerade1(1)*rho_gerade2(2))&
1548 : 12847284 : *ewald_fun/eta/dielt_perp*cos(phi)
1549 : : ! Second, add the source charges (ungerade ungerade)
1550 : : dyew_rec(1,idir1,ipert1,idir2,ipert2)= dyew_rec(1,idir1,ipert1,idir2,ipert2)+&
1551 : : (rho_ungerade1(1)*rho_ungerade2(1)+rho_ungerade1(2)*rho_ungerade2(2))&
1552 : 12847284 : *ewald_fun2/eta1/dielt_perp1*cos(phi)
1553 : : dyew_rec(1,idir1,ipert1,idir2,ipert2)= dyew_rec(1,idir1,ipert1,idir2,ipert2)-&
1554 : : (rho_ungerade1(2)*rho_ungerade2(1)-rho_ungerade1(1)*rho_ungerade2(2))&
1555 : 12847284 : *ewald_fun2/eta1/dielt_perp1*sin(phi)
1556 : : dyew_rec(2,idir1,ipert1,idir2,ipert2)= dyew_rec(2,idir1,ipert1,idir2,ipert2)+&
1557 : : (rho_ungerade1(1)*rho_ungerade2(1)+rho_ungerade1(2)*rho_ungerade2(2))&
1558 : 12847284 : *ewald_fun2/eta1/dielt_perp1*sin(phi)
1559 : : dyew_rec(2,idir1,ipert1,idir2,ipert2)= dyew_rec(2,idir1,ipert1,idir2,ipert2)+&
1560 : : (rho_ungerade1(2)*rho_ungerade2(1)-rho_ungerade1(1)*rho_ungerade2(2))&
1561 : 12847284 : *ewald_fun2/eta1/dielt_perp1*cos(phi)
1562 : : ! Third, add the source charge (gerade ungerade)
1563 : : ! For sake of consistenty, only used when there is only one dielectric thickness
1564 [ + + ]: 12847284 : if (inner_thick <tol6) then
1565 : : dyew_rec(1,idir1,ipert1,idir2,ipert2)= dyew_rec(1,idir1,ipert1,idir2,ipert2)-&
1566 : : (rho_gerade1(1)*rho_ungerade2(1)+rho_ungerade1(1)*rho_gerade2(1)&
1567 : : +rho_gerade1(2)*rho_ungerade2(2)+rho_ungerade1(2)*rho_gerade2(2)) &
1568 : 12019824 : *ewald_fun1/eta/dielt_perp*cos(phi)
1569 : : dyew_rec(1,idir1,ipert1,idir2,ipert2)= dyew_rec(1,idir1,ipert1,idir2,ipert2)+&
1570 : : (rho_gerade1(2)*rho_ungerade2(1)+rho_ungerade1(2)*rho_gerade2(1)&
1571 : : -rho_gerade1(1)*rho_ungerade2(2)-rho_ungerade1(1)*rho_gerade2(2)) &
1572 : 12019824 : *ewald_fun1/eta/dielt_perp*sin(phi)
1573 : : dyew_rec(2,idir1,ipert1,idir2,ipert2)= dyew_rec(2,idir1,ipert1,idir2,ipert2)-&
1574 : : (rho_gerade1(1)*rho_ungerade2(1)+rho_ungerade1(1)*rho_gerade2(1)&
1575 : : +rho_gerade1(2)*rho_ungerade2(2)+rho_ungerade1(2)*rho_gerade2(2)) &
1576 : 12019824 : *ewald_fun1/eta/dielt_perp*sin(phi)
1577 : : dyew_rec(2,idir1,ipert1,idir2,ipert2)= dyew_rec(2,idir1,ipert1,idir2,ipert2)-&
1578 : : (rho_gerade1(2)*rho_ungerade2(1)+rho_ungerade1(2)*rho_gerade2(1)&
1579 : : -rho_gerade1(1)*rho_ungerade2(2)-rho_ungerade1(1)*rho_gerade2(2)) &
1580 : 12019824 : *ewald_fun1/eta/dielt_perp*cos(phi)
1581 : : end if
1582 : : ! Now add the interactions with the mirror charges... para para
1583 : : dyew_rec(1,idir1,ipert1,idir2,ipert2)= dyew_rec(1,idir1,ipert1,idir2,ipert2)+&
1584 : : (rho_gerade1(1)*rho_gerade2(1)+rho_gerade1(2)*rho_gerade2(2))&
1585 : 12847284 : *mirror_parapara/eta/dielt_perp*cos(phi)
1586 : : dyew_rec(1,idir1,ipert1,idir2,ipert2)= dyew_rec(1,idir1,ipert1,idir2,ipert2)-&
1587 : : (rho_gerade1(2)*rho_gerade2(1)-rho_gerade1(1)*rho_gerade2(2))&
1588 : 12847284 : *mirror_parapara/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_gerade2(1)+rho_gerade1(2)*rho_gerade2(2))&
1591 : 12847284 : *mirror_parapara/eta/dielt_perp*sin(phi)
1592 : : dyew_rec(2,idir1,ipert1,idir2,ipert2)= dyew_rec(2,idir1,ipert1,idir2,ipert2)+&
1593 : : (rho_gerade1(2)*rho_gerade2(1)-rho_gerade1(1)*rho_gerade2(2))&
1594 : 12847284 : *mirror_parapara/eta/dielt_perp*cos(phi)
1595 : : ! Now with perp perp
1596 : : dyew_rec(1,idir1,ipert1,idir2,ipert2)= dyew_rec(1,idir1,ipert1,idir2,ipert2)+&
1597 : : (rho_ungerade1(1)*rho_ungerade2(1)+rho_ungerade1(2)*rho_ungerade2(2))&
1598 : 12847284 : *mirror_perpperp*eta1/dielt_perp1*cos(phi)
1599 : : dyew_rec(1,idir1,ipert1,idir2,ipert2)= dyew_rec(1,idir1,ipert1,idir2,ipert2)-&
1600 : : (rho_ungerade1(2)*rho_ungerade2(1)-rho_ungerade1(1)*rho_ungerade2(2))&
1601 : 12847284 : *mirror_perpperp*eta1/dielt_perp1*sin(phi)
1602 : : dyew_rec(2,idir1,ipert1,idir2,ipert2)= dyew_rec(2,idir1,ipert1,idir2,ipert2)+&
1603 : : (rho_ungerade1(1)*rho_ungerade2(1)+rho_ungerade1(2)*rho_ungerade2(2))&
1604 : 12847284 : *mirror_perpperp*eta1/dielt_perp1*sin(phi)
1605 : : dyew_rec(2,idir1,ipert1,idir2,ipert2)= dyew_rec(2,idir1,ipert1,idir2,ipert2)+&
1606 : : (rho_ungerade1(2)*rho_ungerade2(1)-rho_ungerade1(1)*rho_ungerade2(2))&
1607 : 12847284 : *mirror_perpperp*eta1/dielt_perp1*cos(phi)
1608 : : ! Third, add the source charge (gerade ungerade)
1609 [ + + ]: 17129712 : if (inner_thick <tol6) then
1610 : : dyew_rec(1,idir1,ipert1,idir2,ipert2)= dyew_rec(1,idir1,ipert1,idir2,ipert2)+&
1611 : : (rho_gerade1(1)*rho_ungerade2(1)+rho_ungerade1(1)*rho_gerade2(1)&
1612 : : +rho_gerade1(2)*rho_ungerade2(2)+rho_ungerade1(2)*rho_gerade2(2)) &
1613 : 12019824 : *mirror_paraperp/eta/dielt_perp*cos(phi)
1614 : : dyew_rec(1,idir1,ipert1,idir2,ipert2)= dyew_rec(1,idir1,ipert1,idir2,ipert2)+&
1615 : : (rho_gerade1(2)*rho_ungerade2(1)+rho_ungerade1(2)*rho_gerade2(1)&
1616 : : -rho_gerade1(1)*rho_ungerade2(2)-rho_ungerade1(1)*rho_gerade2(2)) &
1617 : 12019824 : *mirror_paraperp/eta/dielt_perp*sin(phi)
1618 : : dyew_rec(2,idir1,ipert1,idir2,ipert2)= dyew_rec(2,idir1,ipert1,idir2,ipert2)-&
1619 : : (rho_gerade1(1)*rho_ungerade2(1)+rho_ungerade1(1)*rho_gerade2(1)&
1620 : : +rho_gerade1(2)*rho_ungerade2(2)+rho_ungerade1(2)*rho_gerade2(2)) &
1621 : 12019824 : *mirror_paraperp/eta/dielt_perp*sin(phi)
1622 : : dyew_rec(2,idir1,ipert1,idir2,ipert2)= dyew_rec(2,idir1,ipert1,idir2,ipert2)+&
1623 : : (rho_gerade1(2)*rho_ungerade2(1)+rho_ungerade1(2)*rho_gerade2(1)&
1624 : : -rho_gerade1(1)*rho_ungerade2(2)-rho_ungerade1(1)*rho_gerade2(2)) &
1625 : 12019824 : *mirror_paraperp/eta/dielt_perp*cos(phi)
1626 : : ! Finally, there is a term on the sum of charge, only for mirror charges
1627 : : ! Third, add the source charge (gerade ungerade)
1628 : : dyew_rec(1,idir1,ipert1,idir2,ipert2)= dyew_rec(1,idir1,ipert1,idir2,ipert2)-&
1629 : : (rho_gerade1(1)*rho_ungerade2(1)-rho_ungerade1(1)*rho_gerade2(1)&
1630 : : +rho_gerade1(2)*rho_ungerade2(2)-rho_ungerade1(2)*rho_gerade2(2)) &
1631 : 12019824 : *mirror_diff/eta/dielt_perp*cos(phi)
1632 : : dyew_rec(1,idir1,ipert1,idir2,ipert2)= dyew_rec(1,idir1,ipert1,idir2,ipert2)-&
1633 : : (rho_gerade1(2)*rho_ungerade2(1)-rho_ungerade1(2)*rho_gerade2(1)&
1634 : : -rho_gerade1(1)*rho_ungerade2(2)+rho_ungerade1(1)*rho_gerade2(2)) &
1635 : 12019824 : *mirror_diff/eta/dielt_perp*sin(phi)
1636 : : dyew_rec(2,idir1,ipert1,idir2,ipert2)= dyew_rec(2,idir1,ipert1,idir2,ipert2)+&
1637 : : (rho_gerade1(1)*rho_ungerade2(1)-rho_ungerade1(1)*rho_gerade2(1)&
1638 : : +rho_gerade1(2)*rho_ungerade2(2)-rho_ungerade1(2)*rho_gerade2(2)) &
1639 : 12019824 : *mirror_diff/eta/dielt_perp*sin(phi)
1640 : : dyew_rec(2,idir1,ipert1,idir2,ipert2)= dyew_rec(2,idir1,ipert1,idir2,ipert2)-&
1641 : : (rho_gerade1(2)*rho_ungerade2(1)-rho_ungerade1(2)*rho_gerade2(1)&
1642 : : +rho_gerade1(1)*rho_ungerade2(2)+rho_ungerade1(1)*rho_gerade2(2)) &
1643 : 12019824 : *mirror_diff/eta/dielt_perp*sin(phi)
1644 : : end if
1645 : : end do
1646 : : end do
1647 : : end do
1648 : : end do
1649 : : end if
1650 : : end do
1651 : : end do
1652 : 561 : ucsurf = rprimd_para(1,1)*rprimd_para(2,2)-rprimd_para(1,2)*rprimd_para(2,1)
1653 : : ! Renormalize by surface of periodic 2D lattice and out-of-plane dielectric constant
1654 [ + + + + : 208897 : dyew_rec = dyew_rec*(two_pi)/ucsurf
+ + + + +
+ ]
1655 : : ! Reciprocal summation completes. Remains some real-space contribution from first
1656 : : ! unit cells (ipert1 neq ipert2), impacting all IFCs the same way
1657 [ + + + + : 208897 : dyew_real = zero
+ + + + +
+ ]
1658 [ + + ]: 2425 : do ipert1=1,natom
1659 [ + + ]: 9121 : do ipert2=1,natom
1660 [ + + ]: 26784 : diff_xcart(:) =xcart(:,ipert2)-xcart(:,ipert1)
1661 [ + + + + ]: 87048 : rvec_dielt(:) = matmul(invdlt,diff_xcart)
1662 [ + + ]: 26784 : rvec_norm = dot_product(diff_xcart,rvec_dielt)
1663 : 6696 : sqrt_norm = dsqrt(rvec_norm)
1664 : 6696 : fac_erfc = sqrt_norm/dsqrt(two)/lambda
1665 : 6696 : fac_gauss = -rvec_norm/two/lambda**2
1666 [ + + ]: 28648 : do idir1=1,3
1667 [ + + ]: 87048 : do idir2=1,3
1668 [ + + ]: 80352 : if (ipert1 .NE. ipert2) then ! Only off-sites contributions
1669 : : ! First, contribution from eps^-1 (Rk'b-Rka) eps^-1
1670 : : fac_real = three*(one-erf(fac_erfc))/sqrt_norm**5+6*exp(fac_gauss)/rvec_norm**2/&
1671 : 43488 : sqrt(two_pi)/lambda+two*exp(fac_gauss)/rvec_norm/sqrt(two_pi)/lambda**3
1672 : : dyew_real(1,idir1,ipert1,idir2,ipert2) = dyew_real(1,idir1,ipert1,idir2,ipert2)+&
1673 : 43488 : fac_real*rvec_dielt(idir1)*rvec_dielt(idir2)
1674 : : ! Second contribution from esp^-1(alpha,beta)
1675 : 43488 : fac_real = (one-erf(fac_erfc))/sqrt_norm**3+exp(fac_gauss)/rvec_norm/sqrt(two_pi)/lambda
1676 : : dyew_real(1,idir1,ipert1,idir2,ipert2)=dyew_real(1,idir1,ipert1,idir2,ipert2)-&
1677 : 43488 : fac_real*dielt(idir1,idir2)
1678 : : end if
1679 : : end do
1680 : : end do
1681 : : end do
1682 : : end do
1683 [ + + + + : 208897 : dyew_real = dyew_real / dsqrt(detdlt)
+ + + + +
+ ]
1684 [ + + + + : 208897 : dyew = dyew_real + dyew_rec
+ + + + +
+ ]
1685 : 561 : end subroutine ewald9_2D
1686 : :
1687 : : end module m_ewald
1688 : : !!***
1689 : :
|
