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 7644 : subroutine ewald(eew,gmet,grewtn,gsqcut,icutcoul,natom,ngfft,nkpt,ntypat,rcut,&
76 7644 : 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 7644 : real(dp),allocatable :: gcutoff(:)
99 : ! *************************************************************************
100 :
101 : !This is the minimum argument of an exponential, with some safety
102 7644 : minexparg=log(tiny(0._dp))+five
103 :
104 : !Add up total charge and sum of $charge^2$ in cell
105 :
106 7644 : chsq=0._dp
107 7644 : ch=0._dp
108 29994 : do ia=1,natom
109 22350 : ch=ch+zion(typat(ia))
110 29994 : 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 7644 : 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 7644 : 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 7644 : 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 7644 : eta=pi*200.0_dp/33.0_dp*sqrt(1.69_dp*recip/direct)
139 : end if
140 :
141 : !Conduct reciprocal space summations
142 7644 : fac=pi**2/eta
143 7644 : gsum=0._dp
144 97044 : grewtn(:,:)=0.0_dp
145 :
146 : !Initialize Gcut-off array from m_gtermcutoff
147 : !ABI_MALLOC(gcutoff,(ngfft(1)*ngfft(2)*ngfft(3)))
148 7644 : 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 7644 : ng=0
154 : do
155 318596 : ng=ng+1
156 318596 : 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 318596 : ii=1
165 15902136 : do ig3=-ng,ng
166 1241482172 : do ig2=-ng,ng
167 >18129*10^7 : do ig1=-ng,ng
168 : ! Exclude shells previously summed over
169 >18127*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 7169034148 : & gmet(3,1)*dble(ig1*ig3)+gmet(3,2)*dble(ig3*ig2))
175 :
176 : ! Skip g=0:
177 7169034148 : if (gsq>1.0d-20) then
178 7169026504 : arg=fac*gsq
179 :
180 : ! Larger arg gives 0 contribution because of exp(-arg)
181 7169026504 : if (arg <= -minexparg ) then
182 : ! When any term contributes then include next shell
183 1843245502 : newg=1
184 :
185 : if((abs(ig1).lt.ngfft(1)).and.&
186 1843245502 : &(abs(ig2).lt.ngfft(2)).and.&
187 : &(abs(ig3).lt.ngfft(3))) then
188 566536172 : ig23=ngfft(1)*(abs(ig2)+ngfft(2)*(abs(ig3)))
189 566536172 : 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 566536172 : term=exp(-arg)/gsq*gcutoff(ii)
194 1276709330 : else if (icutcoul.ne.3) then
195 : term=zero !exp(-arg)/gsq
196 : else
197 1276487284 : term=exp(-arg)/gsq
198 : endif
199 :
200 1843245502 : summr = 0.0_dp
201 1843245502 : 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 7593912754 : do ia=1,natom
211 5750667252 : 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 5750667252 : summr=summr+zion(typat(ia))*cos(arg)
214 7593912754 : summi=summi+zion(typat(ia))*sin(arg)
215 : end do
216 :
217 : ! The following two checks avoid an annoying underflow error msg
218 1843245502 : if (abs(summr)<1.d-16) summr=0.0_dp
219 1843245502 : 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 1843245502 : t1=term*(summr*summr+summi*summi)
224 1843245502 : gsum=gsum+t1
225 :
226 7593912754 : do ia=1,natom
227 : ! Again only phase is computed so xred may fall outside [0,1).
228 5750667252 : arg=two_pi*(ig1*xred(1,ia)+ig2*xred(2,ia)+ig3*xred(3,ia))
229 5750667252 : phr= cos(arg)
230 5750667252 : phi=-sin(arg)
231 : ! (note: do not need real part, commented out)
232 : ! c1r=(phr*summr-phi*summi)*(term*zion(typat(ia)))
233 5750667252 : c1i=(phi*summr+phr*summi)*(term*zion(typat(ia)))
234 : ! compute coordinate gradients
235 5750667252 : grewtn(1,ia)=grewtn(1,ia)-c1i*ig1
236 5750667252 : grewtn(2,ia)=grewtn(2,ia)-c1i*ig2
237 7593912754 : 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 318596 : 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 7644 : 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 7644 : term=-2._dp/ucvol
261 97044 : grewtn(:,:)=grewtn(:,:)*term
262 : !call DSCAL(3*natom,term,grewtn,1)
263 :
264 : !Conduct real space summations
265 7644 : reta=sqrt(eta)
266 7644 : fac=2._dp*sqrt(eta/pi)
267 7644 : 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 7644 : nr=0
275 : do
276 22241 : nr=nr+1
277 22241 : 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 134942 : do ir3=-nr,nr
287 782687 : do ir2=-nr,nr
288 4912683 : do ir1=-nr,nr
289 4799982 : if( abs(ir3)==nr .or. abs(ir2)==nr .or. abs(ir1)==nr .or. nr==1 )then
290 :
291 11479280 : do ia=1,natom
292 : ! Map reduced coordinate xred(mu,ia) into [0,1)
293 8759650 : fraca1=xred(1,ia)-aint(xred(1,ia))+0.5_dp-sign(0.5_dp,xred(1,ia))
294 8759650 : fraca2=xred(2,ia)-aint(xred(2,ia))+0.5_dp-sign(0.5_dp,xred(2,ia))
295 8759650 : fraca3=xred(3,ia)-aint(xred(3,ia))+0.5_dp-sign(0.5_dp,xred(3,ia))
296 8759650 : drdta1=0.0_dp
297 8759650 : drdta2=0.0_dp
298 8759650 : drdta3=0.0_dp
299 :
300 103188764 : 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 94429114 : fracb1=xred(1,ib)-aint(xred(1,ib))+0.5_dp-sign(0.5_dp,xred(1,ib))
304 94429114 : fracb2=xred(2,ib)-aint(xred(2,ib))+0.5_dp-sign(0.5_dp,xred(2,ib))
305 94429114 : fracb3=xred(3,ib)-aint(xred(3,ib))+0.5_dp-sign(0.5_dp,xred(3,ib))
306 94429114 : r1=dble(ir1)+fracb1-fraca1
307 94429114 : r2=dble(ir2)+fracb2-fraca2
308 94429114 : r3=dble(ir3)+fracb3-fraca3
309 : rsq=rmet(1,1)*r1*r1+rmet(2,2)*r2*r2+rmet(3,3)*r3*r3+&
310 94429114 : & 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 103188764 : 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 94406764 : term=0._dp
318 94406764 : if (eta*rsq<64.0_dp) then
319 5549124 : newr=1
320 5549124 : rmagn=sqrt(rsq)
321 5549124 : arg=reta*rmagn
322 : ! derfc is the real(dp) complementary error function
323 5549124 : derfc_arg = abi_derfc(arg)
324 5549124 : term=derfc_arg/rmagn
325 5549124 : sumr=sumr+zion(typat(ia))*zion(typat(ib))*term
326 : term=zion(typat(ia))*zion(typat(ib))*&
327 5549124 : & (term+fac*exp(-eta*rsq))/rsq
328 : ! Length scale grads now handled with stress tensor in ewald2
329 5549124 : r1a1d=rmet(1,1)*r1+rmet(1,2)*r2+rmet(1,3)*r3
330 5549124 : r2a2d=rmet(2,1)*r1+rmet(2,2)*r2+rmet(2,3)*r3
331 5549124 : r3a3d=rmet(3,1)*r1+rmet(3,2)*r2+rmet(3,3)*r3
332 : ! Compute terms related to coordinate gradients
333 5549124 : drdta1=drdta1+term*r1a1d
334 5549124 : drdta2=drdta2+term*r2a2d
335 5549124 : 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 8759650 : grewtn(1,ia)=grewtn(1,ia)+drdta1
341 8759650 : grewtn(2,ia)=grewtn(2,ia)+drdta2
342 12911887 : 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 22241 : if(newr==0) exit
351 : end do ! End loop on nr (new shells). Note that there is an exit within the loop
352 : !
353 7644 : sumr=0.5_dp*sumr
354 7644 : fac=pi*ch**2.0_dp/(2.0_dp*eta*ucvol)
355 :
356 : !Finally assemble Ewald energy, eew
357 7644 : 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 7634 : eew=sumg+sumr-chsq*reta/sqrt(pi)-fac
362 : end if
363 :
364 7644 : 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 7644 : 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 6302 : 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 6302 : call matr3inv(rprimd,gprimd)
439 :
440 : !This is the minimum argument of an exponential, with some safety
441 6302 : minexparg=log(tiny(0._dp))+five
442 :
443 : !Add up total charge and sum of charge^2 in cell
444 6302 : ch=0._dp
445 24956 : do ia=1,natom
446 24956 : 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 6302 : & 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 6302 : & 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 6302 : eta=pi*200.0_dp/33.0_dp*sqrt(1.69_dp*recip/direct)
458 :
459 6302 : fac=pi**2/eta
460 :
461 : !Conduct reciprocal space summations
462 6302 : strg(1:6)=0.0_dp
463 :
464 : !Sum over G space, done shell after shell until all
465 : !contributions are too small
466 6302 : ng=0
467 : do
468 270096 : ng=ng+1
469 270096 : newg=0
470 :
471 13808420 : do ig3=-ng,ng
472 1110349444 : do ig2=-ng,ng
473 >17019*10^7 : do ig1=-ng,ng
474 :
475 : ! Exclude shells previously summed over
476 >17017*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 6418953326 : g1=gprimd(1,1)*ig1+gprimd(1,2)*ig2+gprimd(1,3)*ig3
481 6418953326 : g2=gprimd(2,1)*ig1+gprimd(2,2)*ig2+gprimd(2,3)*ig3
482 6418953326 : g3=gprimd(3,1)*ig1+gprimd(3,2)*ig2+gprimd(3,3)*ig3
483 : ! Compute |G|^2 (no pi factors)
484 6418953326 : gsq=(g1**2+g2**2+g3**2)
485 :
486 : ! skip g=0:
487 6418953326 : if (gsq>1.0d-20) then
488 6418947024 : arg1=fac*gsq
489 :
490 : ! larger arg1 gives 0 contribution because of exp(-arg1)
491 6418947024 : if (arg1<= -minexparg) then
492 : ! When any term contributes then include next shell
493 1644036704 : newg=1
494 1644036704 : term1=exp(-arg1)/arg1
495 1644036704 : summr = 0.0_dp
496 1644036704 : summi = 0.0_dp
497 6750263198 : do ia=1,natom
498 5106226494 : 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 5106226494 : summr=summr+zion(typat(ia))*cos(arg2)
501 6750263198 : summi=summi+zion(typat(ia))*sin(arg2)
502 : end do
503 :
504 : ! Avoid underflow error messages
505 1644036704 : if (abs(summr)<1.d-16) summr=0.0_dp
506 1644036704 : if (abs(summi)<1.d-16) summi=0.0_dp
507 :
508 1644036704 : term2=(2._dp/gsq)*(1._dp+arg1)
509 1644036704 : t1=term2*g1*g1-1._dp
510 1644036704 : t2=term2*g2*g2-1._dp
511 1644036704 : t3=term2*g3*g3-1._dp
512 1644036704 : t4=term2*g2*g3
513 1644036704 : t5=term2*g1*g3
514 1644036704 : t6=term2*g1*g2
515 1644036704 : term3=term1*(summr*summr+summi*summi)
516 1644036704 : strg(1)=strg(1)+t1*term3
517 1644036704 : strg(2)=strg(2)+t2*term3
518 1644036704 : strg(3)=strg(3)+t3*term3
519 1644036704 : strg(4)=strg(4)+t4*term3
520 1644036704 : strg(5)=strg(5)+t5*term3
521 1644036704 : 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 270096 : 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 6302 : reta=sqrt(eta)
538 6302 : strr(1:6)=0.0_dp
539 :
540 : !Loop on shells in r-space as was done in g-space
541 6302 : nr=0
542 : do
543 18130 : nr=nr+1
544 18130 : newr=0
545 :
546 109696 : do ir3=-nr,nr
547 635250 : do ir2=-nr,nr
548 3993998 : do ir1=-nr,nr
549 3902432 : if( abs(ir3)==nr .or. abs(ir2)==nr .or. abs(ir1)==nr .or. nr==1 )then
550 :
551 9405204 : do ia=1,natom
552 : ! Convert reduced atomic coordinates to [0,1)
553 7199330 : fraca1=xred(1,ia)-aint(xred(1,ia))+0.5_dp-sign(0.5_dp,xred(1,ia))
554 7199330 : fraca2=xred(2,ia)-aint(xred(2,ia))+0.5_dp-sign(0.5_dp,xred(2,ia))
555 7199330 : fraca3=xred(3,ia)-aint(xred(3,ia))+0.5_dp-sign(0.5_dp,xred(3,ia))
556 92731002 : do ib=1,natom
557 82154794 : fracb1=xred(1,ib)-aint(xred(1,ib))+0.5_dp-sign(0.5_dp,xred(1,ib))
558 82154794 : fracb2=xred(2,ib)-aint(xred(2,ib))+0.5_dp-sign(0.5_dp,xred(2,ib))
559 82154794 : fracb3=xred(3,ib)-aint(xred(3,ib))+0.5_dp-sign(0.5_dp,xred(3,ib))
560 82154794 : r1=ir1+fracb1-fraca1
561 82154794 : r2=ir2+fracb2-fraca2
562 82154794 : r3=ir3+fracb3-fraca3
563 : ! Convert from reduced to cartesian coordinates
564 82154794 : r1c=rprimd(1,1)*r1+rprimd(1,2)*r2+rprimd(1,3)*r3
565 82154794 : r2c=rprimd(2,1)*r1+rprimd(2,2)*r2+rprimd(2,3)*r3
566 82154794 : r3c=rprimd(3,1)*r1+rprimd(3,2)*r2+rprimd(3,3)*r3
567 : ! Compute |r|^2
568 82154794 : rsq=r1c**2+r2c**2+r3c**2
569 82154794 : rmagn=sqrt(rsq)
570 :
571 : ! Avoid zero denominators in 'term':
572 89354124 : 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 82136140 : arg3=reta*rmagn
577 82136140 : if (arg3<8.0_dp) then
578 4500952 : newr=1
579 : ! derfc computes the complementary error function
580 : ! dderfc is the derivative of the complementary error function
581 4500952 : dderfc=(-2/sqrt(pi))*exp(-eta*rsq)
582 4500952 : derfc_arg = abi_derfc(arg3)
583 4500952 : term3=dderfc-derfc_arg/arg3
584 4500952 : term4=zion(typat(ia))*zion(typat(ib))*term3
585 4500952 : strr(1)=strr(1)+term4*r1c*r1c/rsq
586 4500952 : strr(2)=strr(2)+term4*r2c*r2c/rsq
587 4500952 : strr(3)=strr(3)+term4*r3c*r3c/rsq
588 4500952 : strr(4)=strr(4)+term4*r2c*r3c/rsq
589 4500952 : strr(5)=strr(5)+term4*r1c*r3c/rsq
590 4500952 : 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 18130 : 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 6302 : fac = pi/(2._dp*ucvol*eta)
611 6302 : stress(1)=(0.5_dp*reta*strr(1)+fac*(strg(1)+(ch**2)))/ucvol
612 6302 : stress(2)=(0.5_dp*reta*strr(2)+fac*(strg(2)+(ch**2)))/ucvol
613 6302 : stress(3)=(0.5_dp*reta*strr(3)+fac*(strg(3)+(ch**2)))/ucvol
614 6302 : stress(4)=(0.5_dp*reta*strr(4)+fac*strg(4))/ucvol
615 6302 : stress(5)=(0.5_dp*reta*strr(5)+fac*strg(5))/ucvol
616 6302 : stress(6)=(0.5_dp*reta*strr(6)+fac*strg(6))/ucvol
617 :
618 6302 : 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 181081 : 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 362162 : real(dp) :: c3r(2*mr+1),cosqxred(natom),wdielt(3,3),eig_dielt(3),gpq(3),gpqfac(3,3),gpqgpq(3,3)
710 362162 : real(dp) :: invdlt(3,3),ircar(3),ircax(3),rr(3),sinqxred(natom)
711 362162 : real(dp) :: xredcar(3,natom),xredcax(3,natom),xredicar(3),xredicax(3),xx(3)
712 : real(dp) :: gprimbyacell(3,3) !,tsec(2)
713 181081 : real(dp),allocatable :: dyddt(:,:,:,:,:), dydqt(:,:,:,:,:,:), dyqqt(:,:,:,:,:,:,:)
714 181081 : real(dp),allocatable :: work(:)
715 181081 : complex(dp) :: exp2piqx(natom)
716 181081 : 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 16422890 : if (all(zeff == zero).and.all(qdrp_cart == zero)) then
722 3128282 : dyew = zero; return
723 : end if
724 13103851 : 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 141171 : dipquad_=0; if(present(dipquad)) dipquad_=dipquad
731 141171 : quadquad_=0; if(present(quadquad)) quadquad_=quadquad
732 :
733 : ! Deactivate real space sums for quadrupolar fields or for dipdip = -1
734 141171 : ewald_option = 0; if (present(option)) ewald_option = option
735 141171 : 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 141027 : minexparg=log(tiny(0._dp))+five
740 141171 : if (ewald_option == 1) minexparg=-20.0_dp
741 :
742 : ! initialize complex phase factors
743 491992 : do ia = 1, natom
744 : arga = two_pi*( (qphon(1))*xred(1,ia)&
745 : +(qphon(2))*xred(2,ia)&
746 350821 : +(qphon(3))*xred(3,ia) )
747 491992 : exp2piqx(ia) = exp(arga*j_dpc)
748 : end do
749 141171 : ng_expxq = 1000
750 423513 : ABI_MALLOC(expx1, (-ng_expxq:ng_expxq, natom))
751 282342 : ABI_MALLOC(expx2, (-ng_expxq:ng_expxq, natom))
752 282342 : ABI_MALLOC(expx3, (-ng_expxq:ng_expxq, natom))
753 491992 : do ia = 1, natom
754 702484813 : do ig1 = -ng_expxq, ng_expxq
755 701992821 : expx1(ig1, ia) = exp(ig1*two_pi*xred(1,ia)*j_dpc)
756 701992821 : expx2(ig1, ia) = exp(ig1*two_pi*xred(2,ia)*j_dpc)
757 702343642 : expx3(ig1, ia) = exp(ig1*two_pi*xred(3,ia)*j_dpc)
758 : end do
759 : end do
760 :
761 141171 : gprimbyacell = gprim
762 564684 : gprimbyacell(:,1) = gprimbyacell(:,1) / acell(1)
763 564684 : gprimbyacell(:,2) = gprimbyacell(:,2) / acell(2)
764 564684 : 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 141171 : 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 141171 : gmet(2,2)+gmet(2,3)+gmet(3,1)+gmet(3,2)+gmet(3,3)
771 :
772 141171 : 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 141171 : 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 141171 : inv4eta = one / four / eta
814 :
815 705855 : ABI_MALLOC(dyddt,(2,3,natom,3,natom))
816 705855 : ABI_MALLOC(dydqt,(2,3,natom,3,natom,3))
817 705855 : ABI_MALLOC(dyqqt,(2,3,natom,3,natom,3,3))
818 :
819 37471225 : dyddt = zero
820 112554846 : dydqt = zero
821 337805709 : dyqqt = zero
822 :
823 : !Sum terms over g space:
824 : ng=0
825 : do
826 1898918 : ng=ng+1
827 :
828 : ! if needed, update the complex phases for larger G vectors
829 1898918 : 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 1898918 : newg=0
851 34005968 : do ig3=-ng,ng
852 802421358 : do ig2=-ng,ng
853 27092131162 : do ig1=-ng,ng
854 27060024112 : 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 4240540255 : gprimbyacell(1,2)+(ig3+qphon(3))*gprimbyacell(1,3)
858 : gpq(2)=(ig1+qphon(1))*gprimbyacell(2,1)+(ig2+qphon(2))*&
859 4240540255 : gprimbyacell(2,2)+(ig3+qphon(3))*gprimbyacell(2,3)
860 : gpq(3)=(ig1+qphon(1))*gprimbyacell(3,1)+(ig2+qphon(2))*&
861 4240540255 : gprimbyacell(3,2)+(ig3+qphon(3))*gprimbyacell(3,3)
862 4240540255 : gsq=zero
863 16962161020 : do jj=1,3
864 55127023315 : do ii=1,3
865 38164862295 : gpqgpq(ii,jj)=gpq(ii)*gpq(jj)
866 50886483060 : gsq=gsq+gpqgpq(ii,jj)*dielt(ii,jj)
867 : end do
868 : end do
869 :
870 : ! Skip q=0:
871 4240540255 : 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 4240520764 : arg1=(two_pi**2)*gsq* inv4eta
883 :
884 : ! Larger arg gives 0 contribution:
885 4240520764 : if (arg1<= -minexparg ) then
886 1118988781 : newg=1
887 :
888 : ! Here calculate the term
889 1118988781 : term1=exp(-arg1)/gsq
890 4475955124 : do jj=1,3
891 14546854153 : do ii=1,3
892 13427865372 : 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 4594649808 : do ia=1,natom
899 3475661027 : cosqxred(ia)= real(exp2piqx(ia)*expx1(ig1, ia)*expx2(ig2, ia)*expx3(ig3, ia))
900 4594649808 : sinqxred(ia)=aimag(exp2piqx(ia)*expx1(ig1, ia)*expx2(ig2, ia)*expx3(ig3, ia))
901 : end do
902 :
903 : ! First, the diagonal terms
904 4475955124 : do nu=1,3
905 14902938205 : do ia=1,natom
906 34637915586 : do mu=nu,3
907 31280949243 : 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 3475661027 : do ib=2,natom
914 9936789130 : do ia=1,ib-1
915 : ! phase factor dipole-dipole
916 6461128103 : cddr=cosqxred(ia)*cosqxred(ib)+sinqxred(ia)*sinqxred(ib)
917 6461128103 : cddi=sinqxred(ia)*cosqxred(ib)-cosqxred(ia)*sinqxred(ib)
918 :
919 : ! Dipole-dipole contribution
920 28201184658 : do nu=1,3
921 64611281030 : do mu=nu,3
922 38766768618 : dyddt(1,mu,ia,nu,ib)=dyddt(1,mu,ia,nu,ib)+gpqfac(mu,nu)*cddr
923 58150152927 : 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 1118988781 : if (do_quadrupole) then
930 415900706 : do ib=1,natom
931 997315061 : do ia=1,natom
932 :
933 : ! phase factor for dipole-quadrupole
934 581414355 : cqdr=cosqxred(ia)*sinqxred(ib)-sinqxred(ia)*cosqxred(ib)
935 581414355 : cqdi=cosqxred(ia)*cosqxred(ib)+sinqxred(ia)*sinqxred(ib)
936 :
937 : ! phase factor quadrupole-quadrupole
938 581414355 : cqqr=cosqxred(ia)*cosqxred(ib)+sinqxred(ia)*sinqxred(ib)
939 581414355 : cqqi=sinqxred(ia)*cosqxred(ib)-cosqxred(ia)*sinqxred(ib)
940 :
941 : ! Dipole-quadrupole contribution
942 2325657420 : do ii=1,3
943 7558386615 : do jj=1,3
944 22675159845 : do kk=1,3
945 15698187585 : g3=gpq(ii)*gpq(jj)*gpq(kk)
946 15698187585 : dydqt(1,ii,ia,jj,ib,kk)=dydqt(1,ii,ia,jj,ib,kk)+g3*term1*cqdr
947 20930916780 : 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 2605368201 : do ii=1,3
954 7558386615 : do jj=1,3
955 22675159845 : do kk=1,3
956 68025479535 : do ll=1,3
957 47094562755 : g4 = gpq(ii)*gpq(jj)*gpq(kk)*gpq(ll)
958 47094562755 : dyqqt(1,ii,ia,jj,ib,kk,ll)=dyqqt(1,ii,ia,jj,ib,kk,ll)+g4*term1*cqqr
959 62792750340 : 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 1898918 : if(newg==0)exit
977 : end do
978 :
979 : ! Multiplies by common factor
980 141171 : fact1=4.0_dp*pi/ucvol
981 491992 : do ib=1,natom
982 1266182 : do ia=1,ib
983 3447581 : do nu=1,3
984 7741900 : do mu=nu,3
985 4645140 : dyddt(1,mu,ia,nu,ib)=dyddt(1,mu,ia,nu,ib)*fact1
986 6967710 : 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 141171 : if (do_quadrupole) then
992 6039836 : dydqt=dydqt*fact1/two * two_pi
993 18131300 : dyqqt=dyqqt*fact1/four * two_pi ** 2
994 : end if
995 :
996 141171 : reta=sqrt(eta)
997 141171 : reta3=-eta*reta
998 :
999 : !Calculating the inverse (transpose) of the dielectric tensor
1000 141171 : 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 141171 : dielt(1,2)*dielt(2,1)*dielt(3,3)
1007 :
1008 141171 : 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 141171 : inv_detdlt = one / sqrt(detdlt)
1020 141171 : fact3=reta3 * inv_detdlt
1021 :
1022 141171 : if (ewald_option /= 1) then
1023 : ! Preparing the loop on real space
1024 491557 : do ia=1,natom
1025 1543150 : 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 1402124 : xred(3,ia)*acell(3)*rprim(ii,3) )*reta
1029 : end do
1030 : end do
1031 491557 : do ia=1,natom
1032 1543150 : do ii=1,3
1033 : xredcax(ii,ia)= invdlt(1,ii)*xredcar(ii,ia)+&
1034 : invdlt(2,ii)*xredcar(ii,ia)+&
1035 1402124 : invdlt(3,ii)*xredcar(ii,ia)
1036 : end do
1037 : end do
1038 :
1039 : ! Prepare the evaluation of exp(iq*R)
1040 2820802052 : do ir=-mr,mr
1041 2820661026 : arg1=-two_pi*qphon(1)*ir
1042 2820661026 : arg2=-two_pi*qphon(2)*ir
1043 2820661026 : arg3=-two_pi*qphon(3)*ir
1044 2820661026 : c1r(ir+mr+1)=cos(arg1)
1045 2820661026 : c1i(ir+mr+1)=sin(arg1)
1046 2820661026 : c2r(ir+mr+1)=cos(arg2)
1047 2820661026 : c2i(ir+mr+1)=sin(arg2)
1048 2820661026 : c3r(ir+mr+1)=cos(arg3)
1049 2820802052 : c3i(ir+mr+1)=sin(arg3)
1050 : end do
1051 :
1052 1216218 : do nr=1,mr
1053 1216218 : newr=0
1054 :
1055 : ! Begin big loop on real space vectors
1056 15491768 : do ir3=-nr,nr
1057 242017034 : do ir2=-nr,nr
1058 :
1059 : ! Here, construct the cosine and sine of q*R for components 2 and 3
1060 226525266 : c23r = c2r(ir2+mr+1) * c3r(ir3+mr+1) - c2i(ir2+mr+1) * c3i(ir3+mr+1)
1061 226525266 : 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 226525266 : c23r=c23r * fact3
1065 226525266 : c23i=c23i * fact3
1066 :
1067 4669927318 : do ir1=-nr,nr
1068 4655651768 : 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 1197715766 : c123r = c1r(ir1+mr+1) * c23r - c1i(ir1+mr+1) * c23i
1072 1197715766 : c123i = c1i(ir1+mr+1) * c23r + c1r(ir1+mr+1) * c23i
1073 :
1074 4790863064 : do ii=1,3
1075 : ircar(ii)= ( ir1*acell(1)*rprim(ii,1)+&
1076 : ir2*acell(2)*rprim(ii,2)+&
1077 4790863064 : ir3*acell(3)*rprim(ii,3) ) * reta
1078 : end do
1079 4790863064 : do ii=1,3
1080 : ircax(ii)= invdlt(1,ii)*ircar(ii)+&
1081 : invdlt(2,ii)*ircar(ii)+&
1082 4790863064 : invdlt(3,ii)*ircar(ii)
1083 : end do
1084 :
1085 : ! Here loops on atoms
1086 4111740513 : do ib=1,natom
1087 11656098988 : do ii=1,3
1088 8742074241 : xredicar(ii)=ircar(ii)-xredcar(ii,ib)
1089 11656098988 : xredicax(ii)=ircax(ii)-xredcax(ii,ib)
1090 : end do
1091 13018818028 : do ia=1,ib
1092 22702667116 : do ii=1,3
1093 17027000337 : rr(ii)=xredicar(ii)+xredcar(ii,ia)
1094 22702667116 : xx(ii)=xredicax(ii)+xredcax(ii,ia)
1095 : end do
1096 :
1097 5675666779 : 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 8589691526 : 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 1307663559 : if (y2 >= y2min) then
1106 1307313028 : newr=1
1107 1307313028 : yy=sqrt(y2)
1108 1307313028 : invy=1.0_dp/yy
1109 1307313028 : invy2=invy**2
1110 1307313028 : derfc_yy = abi_derfc(yy)
1111 1307313028 : term2=derfc_yy*invy*invy2
1112 1307313028 : term3=fact2*exp(-y2)*invy2
1113 1307313028 : term4=-(term2+term3)
1114 1307313028 : term5=(3.0_dp*term2+term3*(3.0_dp+2.0_dp*y2))*invy2
1115 5229252112 : do nu=1,3
1116 13073130280 : do mu=nu,3
1117 7843878168 : dyddt(1,mu,ia,nu,ib)=dyddt(1,mu,ia,nu,ib)+c123r*(xx(nu)*xx(mu)*term5+term4*invdlt(nu,mu))
1118 11765817252 : 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 350531 : 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 1402124 : do nu=1,3
1132 4556903 : do mu=1,3
1133 : dyddt(1,mu,ia,nu,ib)=dyddt(1,mu,ia,nu,ib)+&
1134 4206372 : 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 1216218 : if(newr==0)exit
1149 1216218 : 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 350821 : do ib=1,natom-1
1155 979771 : do nu=1,3
1156 2108707 : do ia=ib+1,natom
1157 4439271 : do mu=nu,3
1158 2540214 : dyddt(1,mu,ia,nu,ib)= dyddt(1,mu,ib,nu,ia)
1159 3810321 : 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 491992 : do ib=1,natom
1166 1193634 : do nu=2,3
1167 3447581 : do ia=1,natom
1168 6689437 : do mu=1,nu-1
1169 3592677 : dyddt(1,mu,ia,nu,ib)=dyddt(1,nu,ia,mu,ib)
1170 5987795 : 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 37471225 : dyew = zero
1179 491992 : do ib=1,natom
1180 1544455 : do nu=1,3
1181 4995961 : do ia=1,natom
1182 15423171 : do mu=1,3
1183 46704801 : do ii=1,3
1184 140114403 : do jj=1,3
1185 : ! dipole-dipole correction
1186 : dyew(1,mu,ia,nu,ib)=dyew(1,mu,ia,nu,ib) + &
1187 97002279 : 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 97002279 : zeff(ii,mu,ia)*zeff(jj,nu,ib)*dyddt(2,ii,ia,jj,ib)
1190 129336372 : if (do_quadrupole) then
1191 20417832 : do kk=1,3
1192 15313374 : 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 20417832 : 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 141171 : ABI_FREE(expx1)
1222 141171 : ABI_FREE(expx2)
1223 141171 : ABI_FREE(expx3)
1224 141171 : ABI_FREE(dyddt)
1225 141171 : ABI_FREE(dydqt)
1226 141171 : 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 :
|
