TABLE OF CONTENTS


ABINIT/m_ewald [ Modules ]

[ Top ] [ Modules ]

NAME

  m_ewald

FUNCTION

  This module gathers routines to compute the Ewald energy and its derivatives

COPYRIGHT

  Copyright (C) 2014-2018 ABINIT group (DCA, XG, JJC, GMR)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_ewald
28 
29  use defs_basis
30  use m_abicore
31  use m_errors
32  use m_splines
33 
34  use m_special_funcs,  only : abi_derfc
35  use m_symtk,          only : matr3inv
36 
37  implicit none
38 
39  private
40 
41  public :: ewald    ! Compute Ewald energy and derivatives with respect to xred
42  public :: ewald2   ! Derivative of the Ewald energy with respect to strain.
43  public :: ewald9   ! Compute ewald contribution to the dynamical matrix, at a given
44                     ! q wavevector, including anisotropic dielectric tensor and effective charges
45 
46 contains

m_ewald/ewald [ Functions ]

[ Top ] [ m_ewald ] [ Functions ]

NAME

 ewald

FUNCTION

 Compute Ewald energy and derivatives with respect to dimensionless
 reduced atom coordinates xred.

INPUTS

 gmet(3,3)=metric tensor in reciprocal space (bohr^-2)
 natom=number of atoms in unit cell
 ntypat=numbe of type of atoms
 rmet(3,3)=metric tensor in real space (bohr^2)
 typat(natom)=integer label of each type of atom (1,2,...)
 ucvol=unit cell volume (bohr^3)
 xred(3,natom)=relative coords of atoms in unit cell (dimensionless)
 zion(ntypat)=charge on each type of atom (real number)

OUTPUT

 eew=final ewald energy in hartrees
 grewtn(3,natom)=grads of eew wrt xred(3,natom), hartrees.

PARENTS

      setvtr

CHILDREN

      matr3inv,spline

SOURCE

 80 subroutine ewald(eew,gmet,grewtn,natom,ntypat,rmet,typat,ucvol,xred,zion)
 81 
 82 
 83 !This section has been created automatically by the script Abilint (TD).
 84 !Do not modify the following lines by hand.
 85 #undef ABI_FUNC
 86 #define ABI_FUNC 'ewald'
 87 !End of the abilint section
 88 
 89  implicit none
 90 
 91 !Arguments ------------------------------------
 92 !scalars
 93  integer,intent(in) :: natom,ntypat
 94  real(dp),intent(in) :: ucvol
 95  real(dp),intent(out) :: eew
 96 !arrays
 97  integer,intent(in) :: typat(natom)
 98  real(dp),intent(in) :: gmet(3,3),rmet(3,3),xred(3,natom),zion(ntypat)
 99  real(dp),intent(out) :: grewtn(3,natom)
100 
101 !Local variables-------------------------------
102 !scalars
103  integer :: ia,ib,ig1,ig2,ig3,ir1,ir2,ir3,newg,newr,ng,nr
104  real(dp) :: arg,c1i,ch,chsq,derfc_arg,direct,drdta1,drdta2,drdta3,eta,fac
105  real(dp) :: fraca1,fraca2,fraca3,fracb1,fracb2,fracb3,gsq,gsum,phi,phr,r1
106  real(dp) :: minexparg
107  real(dp) :: r1a1d,r2,r2a2d,r3,r3a3d,recip,reta,rmagn,rsq,sumg,summi,summr,sumr
108  real(dp) :: t1,term
109  character(len=500) :: message
110 
111 ! *************************************************************************
112 
113 !This is the minimum argument of an exponential, with some safety
114  minexparg=log(tiny(0._dp))+five
115 
116 !Add up total charge and sum of $charge^2$ in cell
117 
118  chsq=0._dp
119  ch=0._dp
120  do ia=1,natom
121    ch=ch+zion(typat(ia))
122    chsq=chsq+zion(typat(ia))**2
123  end do
124 
125 !Compute eta, the Ewald summation convergence parameter,
126 !for approximately optimized summations:
127  direct=rmet(1,1)+rmet(1,2)+rmet(1,3)+rmet(2,1)+&
128 & rmet(2,2)+rmet(2,3)+rmet(3,1)+rmet(3,2)+rmet(3,3)
129  recip=gmet(1,1)+gmet(1,2)+gmet(1,3)+gmet(2,1)+&
130 & gmet(2,2)+gmet(2,3)+gmet(3,1)+gmet(3,2)+gmet(3,3)
131 !A bias is introduced, because G-space summation scales
132 !better than r space summation ! Note : debugging is the most
133 !easier at fixed eta.
134  eta=pi*200.0_dp/33.0_dp*sqrt(1.69_dp*recip/direct)
135 
136 !Conduct reciprocal space summations
137  fac=pi**2/eta
138  gsum=0._dp
139  grewtn(:,:)=0.0_dp
140 
141 !Sum over G space, done shell after shell until all
142 !contributions are too small.
143  ng=0
144  do
145    ng=ng+1
146    newg=0
147 !   Instead of this warning that most normal users do not understand (because they are doing GS calculations, and not RF calculations), 
148 !   one should optimize this routine. But usually this is a very small fraction of any ABINIT run.
149 !   if (ng > 20 .and. mod(ng,10)==0) then
150 !      write (message,'(3a,I10)') "Very large box of G neighbors in ewald: you probably do not want to do this.", ch10,&
151 !&       " If you have a metal consider setting dipdip 0.  ng = ", ng
152 !      MSG_WARNING(message)
153 !   end if
154 
155    do ig3=-ng,ng
156      do ig2=-ng,ng
157        do ig1=-ng,ng
158 
159 !        Exclude shells previously summed over
160          if(abs(ig1)==ng .or. abs(ig2)==ng .or. abs(ig3)==ng .or. ng==1 ) then
161 
162 !          gsq is G dot G = |G|^2
163            gsq=gmet(1,1)*dble(ig1*ig1)+gmet(2,2)*dble(ig2*ig2)+&
164 &           gmet(3,3)*dble(ig3*ig3)+2._dp*(gmet(2,1)*dble(ig1*ig2)+&
165 &           gmet(3,1)*dble(ig1*ig3)+gmet(3,2)*dble(ig3*ig2))
166 
167 !          Skip g=0:
168            if (gsq>1.0d-20) then
169              arg=fac*gsq
170 
171 !            Larger arg gives 0 contribution because of exp(-arg)
172              if (arg <= -minexparg ) then
173 !              When any term contributes then include next shell
174                newg=1
175                term=exp(-arg)/gsq
176                summr = 0.0_dp
177                summi = 0.0_dp
178 
179 !              XG 20180531  : the two do-loops on ia should be merged, in order to spare
180 !              the waste of computing twice the sin and cos.
181 
182 !              Note that if reduced atomic coordinates xred drift outside
183 !              of unit cell (outside [0,1)) it is irrelevant in the following
184 !              term, which only computes a phase.
185                do ia=1,natom
186                  arg=two_pi*(ig1*xred(1,ia)+ig2*xred(2,ia)+ig3*xred(3,ia))
187 !                Sum real and imaginary parts (avoid complex variables)
188                  summr=summr+zion(typat(ia))*cos(arg)
189                  summi=summi+zion(typat(ia))*sin(arg)
190                end do
191 
192 !              The following two checks avoid an annoying underflow error message
193                if (abs(summr)<1.d-16) summr=0.0_dp
194                if (abs(summi)<1.d-16) summi=0.0_dp
195 
196 !              The product of term and summr**2 or summi**2 below
197 !              can underflow if not for checks above
198                t1=term*(summr*summr+summi*summi)
199                gsum=gsum+t1
200 
201                do ia=1,natom
202 !                Again only phase is computed so xred may fall outside [0,1).
203                  arg=two_pi*(ig1*xred(1,ia)+ig2*xred(2,ia)+ig3*xred(3,ia))
204                  phr= cos(arg)
205                  phi=-sin(arg)
206 !                (note: do not need real part, commented out)
207 !                c1r=(phr*summr-phi*summi)*(term*zion(typat(ia)))
208                  c1i=(phi*summr+phr*summi)*(term*zion(typat(ia)))
209 !                compute coordinate gradients
210                  grewtn(1,ia)=grewtn(1,ia)-c1i*ig1
211                  grewtn(2,ia)=grewtn(2,ia)-c1i*ig2
212                  grewtn(3,ia)=grewtn(3,ia)-c1i*ig3
213                end do
214 
215              end if ! End condition of not larger than -minexparg
216            end if ! End skip g=0
217          end if ! End triple loop over G s and associated new shell condition
218 
219        end do
220      end do
221    end do
222 
223 !  Check if new shell must be calculated
224    if (newg==0) exit
225 
226  end do !  End the loop on ng (new shells). Note that there is one exit from this loop.
227 
228  sumg=gsum/(two_pi*ucvol)
229 
230 !Stress tensor is now computed elsewhere (ewald2) hence do not need
231 !length scale gradients (used to compute them here).
232 
233 !normalize coordinate gradients by unit cell volume ucvol
234  term=-2._dp/ucvol
235  grewtn(:,:)=grewtn(:,:)*term
236 !call DSCAL(3*natom,term,grewtn,1)
237 
238 !Conduct real space summations
239  reta=sqrt(eta)
240  fac=2._dp*sqrt(eta/pi)
241  sumr=0.0_dp
242 
243 !In the following a summation is being conducted over all
244 !unit cells (ir1, ir2, ir3) so it is appropriate to map all
245 !reduced coordinates xred back into [0,1).
246 !
247 !Loop on shells in r-space as was done in g-space
248  nr=0
249  do
250    nr=nr+1
251    newr=0
252 !   Instead of this warning that most normal users do not understand (because they are doing GS calculations, and not RF calculations),
253 !   one should optimize this routine. But usually this is a very small fraction of any ABINIT run.
254 !   if (nr > 20 .and. mod(nr,10)==0) then
255 !      write (message,'(3a,I10)') "Very large box of R neighbors in ewald: you probably do not want to do this.", ch10,&
256 !&       " If you have a metal consider setting dipdip 0.  nr = ", nr
257 !      MSG_WARNING(message)
258 !   end if
259 !
260    do ir3=-nr,nr
261      do ir2=-nr,nr
262        do ir1=-nr,nr
263          if( abs(ir3)==nr .or. abs(ir2)==nr .or. abs(ir1)==nr .or. nr==1 )then
264 
265            do ia=1,natom
266 !            Map reduced coordinate xred(mu,ia) into [0,1)
267              fraca1=xred(1,ia)-aint(xred(1,ia))+0.5_dp-sign(0.5_dp,xred(1,ia))
268              fraca2=xred(2,ia)-aint(xred(2,ia))+0.5_dp-sign(0.5_dp,xred(2,ia))
269              fraca3=xred(3,ia)-aint(xred(3,ia))+0.5_dp-sign(0.5_dp,xred(3,ia))
270              drdta1=0.0_dp
271              drdta2=0.0_dp
272              drdta3=0.0_dp
273 
274              do ib=1,natom
275 !              fraca and fracb should be precomputedi and become arrays with natom dimension. 
276 !              Also the combination with dble(ir1), dble(ir2), dble(ir3) or fraca should be done outside of the ib loop.
277                fracb1=xred(1,ib)-aint(xred(1,ib))+0.5_dp-sign(0.5_dp,xred(1,ib))
278                fracb2=xred(2,ib)-aint(xred(2,ib))+0.5_dp-sign(0.5_dp,xred(2,ib))
279                fracb3=xred(3,ib)-aint(xred(3,ib))+0.5_dp-sign(0.5_dp,xred(3,ib))
280                r1=dble(ir1)+fracb1-fraca1
281                r2=dble(ir2)+fracb2-fraca2
282                r3=dble(ir3)+fracb3-fraca3
283                rsq=rmet(1,1)*r1*r1+rmet(2,2)*r2*r2+rmet(3,3)*r3*r3+&
284 &               2.0_dp*(rmet(2,1)*r2*r1+rmet(3,2)*r3*r2+rmet(3,1)*r1*r3)
285 
286 !              Avoid zero denominators in 'term':
287                if (rsq>=1.0d-24) then
288 
289 !                Note: erfc(8) is about 1.1e-29, so do not bother with larger arg.
290 !                Also: exp(-64) is about 1.6e-28, so do not bother with larger arg**2 in exp.
291                  term=0._dp
292                  if (eta*rsq<64.0_dp) then
293                    newr=1
294                    rmagn=sqrt(rsq)
295                    arg=reta*rmagn
296 !                  derfc is the real(dp) complementary error function
297                    derfc_arg = abi_derfc(arg)
298                    term=derfc_arg/rmagn
299                    sumr=sumr+zion(typat(ia))*zion(typat(ib))*term
300                    term=zion(typat(ia))*zion(typat(ib))*&
301 &                   (term+fac*exp(-eta*rsq))/rsq
302 !                  Length scale grads now handled with stress tensor in ewald2
303                    r1a1d=rmet(1,1)*r1+rmet(1,2)*r2+rmet(1,3)*r3
304                    r2a2d=rmet(2,1)*r1+rmet(2,2)*r2+rmet(2,3)*r3
305                    r3a3d=rmet(3,1)*r1+rmet(3,2)*r2+rmet(3,3)*r3
306 !                  Compute terms related to coordinate gradients
307                    drdta1=drdta1+term*r1a1d
308                    drdta2=drdta2+term*r2a2d
309                    drdta3=drdta3+term*r3a3d
310                  end if
311                end if ! End avoid zero denominators in'term'
312              end do ! end loop over ib:
313 
314              grewtn(1,ia)=grewtn(1,ia)+drdta1
315              grewtn(2,ia)=grewtn(2,ia)+drdta2
316              grewtn(3,ia)=grewtn(3,ia)+drdta3
317            end do ! end loop over ia:
318          end if
319        end do ! end triple loop over real space points and associated condition of new shell
320      end do
321    end do
322 
323 !  Check if new shell must be calculated
324    if(newr==0) exit
325  end do ! End loop on nr (new shells). Note that there is an exit within the loop
326 !
327  sumr=0.5_dp*sumr
328  fac=pi*ch**2/(2.0_dp*eta*ucvol)
329 
330 !Finally assemble Ewald energy, eew
331  eew=sumg+sumr-chsq*reta/sqrt(pi)-fac
332 
333 !DEBUG
334 !write(std_out,*)'eew=sumg+sumr-chsq*reta/sqrt(pi)-fac'
335 !write(std_out,*)eew,sumg,sumr,chsq*reta/sqrt(pi),fac
336 !ENDDEBUG
337 
338 !Length scale grads handled with stress tensor, ewald2
339 
340 !Output the final values of ng and nr
341 ! write(message, '(a,a,i4,a,i4)' )ch10,' ewald : nr and ng are ',nr,' and ',ng
342 ! call wrtout(std_out,message,'COLL')
343 
344 end subroutine ewald

m_ewald/ewald2 [ Functions ]

[ Top ] [ m_ewald ] [ Functions ]

NAME

 ewald2

FUNCTION

 Compute the part of the stress tensor coming from the Ewald energy
 which is calculated by derivating the Ewald energy with respect to strain.
 See Nielsen and Martin, Phys. Rev. B 32, 3792 (1985) [[cite:Nielsen1985a]].
 Definition of stress tensor is $(1/ucvol)*d(Etot)/d(strain(a,b))$.

INPUTS

 gmet(3,3)=metric tensor in reciprocal space (bohr^-2)
 natom=number of atoms in umit cell
 ntypat=number of type of atoms
 rmet(3,3)=metric tensor in real space (bohr^2) (inverse transpose of gmet)
 rprimd(3,3)=dimensional primitive translations in real space (bohr)
 typat(natom)=integer label of each type of atom (1,2,...)
 ucvol=unit cell volume (bohr^3)
 xred(3,natom)=relative coords of atoms in unit cell (dimensionless)
 zion(ntypat)=charge on each type of atom (real number)

OUTPUT

 $stress(6)=(1/ucvol)*gradient$ of Ewald energy with respect to strain,
      in hartrees/bohr^3
 Cartesian components of stress are provided for this symmetric
 tensor in the order 11 22 33 32 31 21.

PARENTS

      stress

CHILDREN

      matr3inv,spline

SOURCE

385 subroutine ewald2(gmet,natom,ntypat,rmet,rprimd,stress,typat,ucvol,xred,zion)
386 
387 
388 !This section has been created automatically by the script Abilint (TD).
389 !Do not modify the following lines by hand.
390 #undef ABI_FUNC
391 #define ABI_FUNC 'ewald2'
392 !End of the abilint section
393 
394  implicit none
395 
396 !Arguments ------------------------------------
397 !scalars
398  integer,intent(in) :: natom,ntypat
399  real(dp),intent(in) :: ucvol
400 !arrays
401  integer,intent(in) :: typat(natom)
402  real(dp),intent(in) :: gmet(3,3),rmet(3,3),rprimd(3,3),xred(3,natom)
403  real(dp),intent(in) :: zion(ntypat)
404  real(dp),intent(out) :: stress(6)
405 
406 !Local variables-------------------------------
407 !scalars
408  integer :: ia,ib,ig1,ig2,ig3,ir1,ir2,ir3,newg,newr,ng,nr
409  real(dp) :: arg1,arg2,arg3,ch,dderfc,derfc_arg,direct,eta,fac,fraca1
410  real(dp) :: fraca2,fraca3,fracb1,fracb2,fracb3,g1,g2,g3,gsq,r1,r1c,r2,r2c
411  real(dp) :: minexparg
412  real(dp) :: r3,r3c,recip,reta,rmagn,rsq,summi,summr,t1,t2,t3,t4,t5,t6,term1
413  real(dp) :: term2,term3,term4
414 !arrays
415  real(dp) :: gprimd(3,3),strg(6),strr(6)
416 
417 ! *************************************************************************
418 
419 !Define dimensional reciprocal space primitive translations gprimd
420 !(inverse transpose of rprimd)
421  call matr3inv(rprimd,gprimd)
422 
423 !This is the minimum argument of an exponential, with some safety
424  minexparg=log(tiny(0._dp))+five
425 
426 !Add up total charge and sum of charge^2 in cell
427  ch=0._dp
428  do ia=1,natom
429    ch=ch+zion(typat(ia))
430  end do
431 
432 !Compute eta, the Ewald summation convergence parameter,
433 !for approximately optimized summations:
434  direct=rmet(1,1)+rmet(1,2)+rmet(1,3)+rmet(2,1)+&
435 & rmet(2,2)+rmet(2,3)+rmet(3,1)+rmet(3,2)+rmet(3,3)
436  recip=gmet(1,1)+gmet(1,2)+gmet(1,3)+gmet(2,1)+&
437 & gmet(2,2)+gmet(2,3)+gmet(3,1)+gmet(3,2)+gmet(3,3)
438 !Here, a bias is introduced, because G-space summation scales
439 !better than r space summation !
440  eta=pi*200.0_dp/33.0_dp*sqrt(1.69_dp*recip/direct)
441 
442  fac=pi**2/eta
443 
444 !Conduct reciprocal space summations
445  strg(1:6)=0.0_dp
446 
447 !Sum over G space, done shell after shell until all
448 !contributions are too small
449  ng=0
450  do
451    ng=ng+1
452    newg=0
453 
454    do ig3=-ng,ng
455      do ig2=-ng,ng
456        do ig1=-ng,ng
457 
458 !        Exclude shells previously summed over
459          if(abs(ig1)==ng .or. abs(ig2)==ng .or. abs(ig3)==ng .or. ng==1 ) then
460 
461 !          Compute Cartesian components of each G
462            g1=gprimd(1,1)*ig1+gprimd(1,2)*ig2+gprimd(1,3)*ig3
463            g2=gprimd(2,1)*ig1+gprimd(2,2)*ig2+gprimd(2,3)*ig3
464            g3=gprimd(3,1)*ig1+gprimd(3,2)*ig2+gprimd(3,3)*ig3
465 !          Compute |G|^2 (no pi factors)
466            gsq=(g1**2+g2**2+g3**2)
467 
468 !          skip g=0:
469            if (gsq>1.0d-20) then
470              arg1=fac*gsq
471 
472 !            larger arg1 gives 0 contribution because of exp(-arg1)
473              if (arg1<= -minexparg) then
474 !              When any term contributes then include next shell
475                newg=1
476                term1=exp(-arg1)/arg1
477                summr = 0.0_dp
478                summi = 0.0_dp
479                do ia=1,natom
480                  arg2=two_pi*(ig1*xred(1,ia)+ig2*xred(2,ia)+ig3*xred(3,ia))
481 !                Sum real and imaginary parts (avoid complex variables)
482                  summr=summr+zion(typat(ia))*cos(arg2)
483                  summi=summi+zion(typat(ia))*sin(arg2)
484                end do
485 
486 !              Avoid underflow error messages
487                if (abs(summr)<1.d-16) summr=0.0_dp
488                if (abs(summi)<1.d-16) summi=0.0_dp
489 
490                term2=(2._dp/gsq)*(1._dp+arg1)
491                t1=term2*g1*g1-1._dp
492                t2=term2*g2*g2-1._dp
493                t3=term2*g3*g3-1._dp
494                t4=term2*g2*g3
495                t5=term2*g1*g3
496                t6=term2*g1*g2
497                term3=term1*(summr*summr+summi*summi)
498                strg(1)=strg(1)+t1*term3
499                strg(2)=strg(2)+t2*term3
500                strg(3)=strg(3)+t3*term3
501                strg(4)=strg(4)+t4*term3
502                strg(5)=strg(5)+t5*term3
503                strg(6)=strg(6)+t6*term3
504 
505              end if ! End condition not being larger than -minexparg
506            end if ! End skip g=0
507 
508          end if ! End triple loop and condition of new shell
509        end do
510      end do
511    end do
512 
513 !  Check if new shell must be calculated
514    if (newg==0) exit
515  end do ! End loop on new shell. Note that there is an "exit" instruction within the loop
516 
517 
518 !Conduct real space summations
519  reta=sqrt(eta)
520  strr(1:6)=0.0_dp
521 
522 !Loop on shells in r-space as was done in g-space
523  nr=0
524  do
525    nr=nr+1
526    newr=0
527 
528    do ir3=-nr,nr
529      do ir2=-nr,nr
530        do ir1=-nr,nr
531          if( abs(ir3)==nr .or. abs(ir2)==nr .or. abs(ir1)==nr .or. nr==1 )then
532 
533            do ia=1,natom
534 !            Convert reduced atomic coordinates to [0,1)
535              fraca1=xred(1,ia)-aint(xred(1,ia))+0.5_dp-sign(0.5_dp,xred(1,ia))
536              fraca2=xred(2,ia)-aint(xred(2,ia))+0.5_dp-sign(0.5_dp,xred(2,ia))
537              fraca3=xred(3,ia)-aint(xred(3,ia))+0.5_dp-sign(0.5_dp,xred(3,ia))
538              do ib=1,natom
539                fracb1=xred(1,ib)-aint(xred(1,ib))+0.5_dp-sign(0.5_dp,xred(1,ib))
540                fracb2=xred(2,ib)-aint(xred(2,ib))+0.5_dp-sign(0.5_dp,xred(2,ib))
541                fracb3=xred(3,ib)-aint(xred(3,ib))+0.5_dp-sign(0.5_dp,xred(3,ib))
542                r1=ir1+fracb1-fraca1
543                r2=ir2+fracb2-fraca2
544                r3=ir3+fracb3-fraca3
545 !              Convert from reduced to cartesian coordinates
546                r1c=rprimd(1,1)*r1+rprimd(1,2)*r2+rprimd(1,3)*r3
547                r2c=rprimd(2,1)*r1+rprimd(2,2)*r2+rprimd(2,3)*r3
548                r3c=rprimd(3,1)*r1+rprimd(3,2)*r2+rprimd(3,3)*r3
549 !              Compute |r|^2
550                rsq=r1c**2+r2c**2+r3c**2
551                rmagn=sqrt(rsq)
552 
553 !              Avoid zero denominators in 'term':
554                if (rmagn>=1.0d-12) then
555 
556 !                Note: erfc(8) is about 1.1e-29, so do not bother with larger arg.
557 !                Also: exp(-64) is about 1.6e-28, so do not bother with larger arg**2 in exp.
558                  arg3=reta*rmagn
559                  if (arg3<8.0_dp) then
560                    newr=1
561 !                  derfc computes the complementary error function
562 !                  dderfc is the derivative of the complementary error function
563                    dderfc=(-2/sqrt(pi))*exp(-eta*rsq)
564                    derfc_arg = abi_derfc(arg3)
565                    term3=dderfc-derfc_arg/arg3
566                    term4=zion(typat(ia))*zion(typat(ib))*term3
567                    strr(1)=strr(1)+term4*r1c*r1c/rsq
568                    strr(2)=strr(2)+term4*r2c*r2c/rsq
569                    strr(3)=strr(3)+term4*r3c*r3c/rsq
570                    strr(4)=strr(4)+term4*r2c*r3c/rsq
571                    strr(5)=strr(5)+term4*r1c*r3c/rsq
572                    strr(6)=strr(6)+term4*r1c*r2c/rsq
573                  end if ! End the condition of not being to large
574                end if ! End avoid zero denominator
575 
576              end do ! End loop over ib:
577            end do  ! End loop over ia:
578 
579          end if ! End triple loop overs real space points, and associated new shell condition
580        end do
581      end do
582    end do
583 
584 !  Check if new shell must be calculated
585    if(newr==0) exit
586  end do ! End loop on new shells
587 
588 !Finally assemble stress tensor coming from Ewald energy, stress
589 !(note division by unit cell volume in accordance with definition
590 !found in Nielsen and Martin, Phys. Rev. B 32, 3792 (1985) [[cite:Nielsen1985a]]
591 
592  fac = pi/(2._dp*ucvol*eta)
593  stress(1)=(0.5_dp*reta*strr(1)+fac*(strg(1)+(ch**2)))/ucvol
594  stress(2)=(0.5_dp*reta*strr(2)+fac*(strg(2)+(ch**2)))/ucvol
595  stress(3)=(0.5_dp*reta*strr(3)+fac*(strg(3)+(ch**2)))/ucvol
596  stress(4)=(0.5_dp*reta*strr(4)+fac*strg(4))/ucvol
597  stress(5)=(0.5_dp*reta*strr(5)+fac*strg(5))/ucvol
598  stress(6)=(0.5_dp*reta*strr(6)+fac*strg(6))/ucvol
599 
600 end subroutine ewald2

m_ewald/ewald9 [ Functions ]

[ Top ] [ m_ewald ] [ Functions ]

NAME

 ewald9

FUNCTION

 Compute ewald contribution to the dynamical matrix, at a given
 q wavevector, including anisotropic dielectric tensor and effective charges
 See Phys. Rev. B 55, 10355 (1997) [[cite:Gonze1997a]], equations (71) to (75).

INPUTS

 acell = lengths by which lattice vectors are multiplied
 dielt(3,3)=dielectric tensor
 gmet(3,3) = metric in reciprocal space.
 gprim(3,3)=dimensionless primitive translations in reciprocal space
 natom=number of atoms in unit cell
 qphon(3)=phonon wavevector (same system of coordinates as the
  reciprocal lattice vectors)
 rmet = metric in real space
 rprim(3,3)=dimensionless primitive translations in real space
 sumg0: if=1, the sum in reciprocal space must include g=0,
  if=0, this contribution must be skipped (q=0 singularity)
 ucvol=unit cell volume in (whatever length scale units)**3
 xred(3,natom)=relative coords of atoms in unit cell (dimensionless)
 zeff(3,3,natom)=effective charge on each atom, versus electric
  field and atomic displacement

OUTPUT

 dyew(2,3,natom,3,natom)= Ewald part of the dynamical matrix,
  second energy derivative wrt xred(3,natom) in Hartrees
 Set to zero if all(zeff == zero)

NOTES

 1. The q=0 part should be subtracted, by another call to
 the present routine, with q=0. The present routine correspond
 to the quantity written A-bar in the explanatory notes.
 If q=0 is asked, sumg0 should be put to 0. Otherwise, it should be put to 1.
 2. Because this routine can be used many times in the
 evaluation of phonons in ppddb9, it has been
 optimized carefully. There is still possibility
 for improvement, by using bloking on G and R!
 3. There can be small numerical variations due to the
 fact that the input dielectric tensor is usually
 not perfectly symmetric ....

PARENTS

      ddb_hybrid,m_dynmat,m_effective_potential,m_ifc

CHILDREN

      matr3inv,spline

SOURCE

 656 subroutine ewald9(acell,dielt,dyew,gmet,gprim,natom,qphon,rmet,rprim,sumg0,ucvol,xred,zeff)
 657 
 658 
 659 !This section has been created automatically by the script Abilint (TD).
 660 !Do not modify the following lines by hand.
 661 #undef ABI_FUNC
 662 #define ABI_FUNC 'ewald9'
 663 !End of the abilint section
 664 
 665  implicit none
 666 
 667 !Arguments -------------------------------
 668 !scalars
 669  integer,intent(in) :: natom,sumg0
 670  real(dp),intent(in) :: ucvol
 671 !arrays
 672  real(dp),intent(in) :: acell(3),dielt(3,3),gmet(3,3),gprim(3,3),qphon(3)
 673  real(dp),intent(in) :: rmet(3,3),rprim(3,3),xred(3,natom),zeff(3,3,natom)
 674  real(dp),intent(out) :: dyew(2,3,natom,3,natom)
 675 
 676 !Local variables -------------------------
 677 !scalars
 678  integer,parameter :: mr=10000,ny2_spline=1024*10
 679  integer :: i2,ia,ib,ig1,ig2,ig3,ii,ir,ir1,ir2,ir3,jj,mu,newg,newr,ng,nr,nu,ng_expxq
 680  real(dp),parameter :: fac=4.0_dp/3.0_dp/sqrt(pi)
 681  real(dp),parameter :: fact2=2.0_dp/sqrt(pi)
 682  real(dp),parameter :: y2max=64.0_dp, y2min=1.0d-24
 683  real(dp) :: arg1,arg2,arg3,arga,c123i,c123r,c23i,c23r,detdlt,inv_detdlt
 684  real(dp) :: direct,eta,fact1,fact3,gsq,recip,reta,reta3,inv4eta
 685  real(dp) :: minexparg
 686  real(dp) :: term1,term2,term3,term4,term5,y2,yy,invy,invy2,derfc_yy
 687  character(len=500) :: message
 688 !arrays
 689  real(dp) :: c1i(2*mr+1),c1r(2*mr+1),c2i(2*mr+1),c2r(2*mr+1),c3i(2*mr+1)
 690  real(dp) :: c3r(2*mr+1),cosqxred(natom),gpq(3),gpqfac(3,3),gpqgpq(3,3)
 691  real(dp) :: invdlt(3,3),ircar(3),ircax(3),rr(3),sinqxred(natom)
 692  real(dp) :: xredcar(3,natom),xredcax(3,natom),xredicar(3),xredicax(3),xx(3)
 693  real(dp) :: gprimbyacell(3,3)
 694  real(dp),allocatable :: dyewt(:,:,:,:,:)
 695  complex(dpc) :: exp2piqx(natom)
 696  complex(dpc),allocatable :: expx1(:,:), expx2(:,:), expx3(:,:)
 697 
 698 !#define DEV_USESPLINE
 699 #ifdef DEV_USESPLINE
 700  integer :: jspl
 701  real(dp) :: aa,bb,cc,dd,step,stepm1,step2div6
 702  logical,save :: first_call=.True.
 703  real(dp),save :: t4spl(ny2_spline,2),t5spl(ny2_spline,2),y2vals(ny2_spline)
 704 #endif
 705 
 706 ! *********************************************************************
 707 
 708  ! This routine is expensive so skip the calculation and return zeros if zeff == zero.
 709  ! Typically this happens when the DDB file does not contains zeff but dipdip = 1 is used (default).
 710  if (all(zeff == zero)) then
 711    dyew = zero; return
 712  end if
 713 
 714 !This is the minimum argument of an exponential, with some safety
 715  minexparg=log(tiny(0._dp))+five
 716 
 717 #ifdef DEV_USESPLINE
 718  step = (0.1_dp + y2max - y2min) / (ny2_spline - 1)
 719  stepm1 = one / step; step2div6 = step**2/six
 720  if (first_call) then
 721    first_call = .False.
 722    do ii=1,ny2_spline
 723      y2 = y2min + (ii-1) * step
 724      y2vals(ii) = y2
 725      yy=sqrt(y2)
 726      invy=1.0_dp/yy
 727      invy2=invy**2
 728      derfc_yy = abi_derfc(yy)
 729      term2=derfc_yy*invy*invy2
 730      term3=fact2*exp(-y2)*invy2
 731      term4=-(term2+term3)
 732      term5=(3.0_dp*term2+term3*(3.0_dp+2.0_dp*y2))*invy2
 733      t4spl(ii,1) = term4
 734      t5spl(ii,1) = term5
 735    end do
 736 
 737    call spline(y2vals, t4spl(:,1), ny2_spline, zero, zero, t4spl(:,2))
 738    call spline(y2vals, t5spl(:,1), ny2_spline, zero, zero, t5spl(:,2))
 739    !do ii=1,ny2_spline
 740    !  write(345,*)y2vals(ii),t4spl(ii,:),t5spl(ii,:)
 741    !end do
 742    !write(std_out,*)"spline tables created"; stop
 743  end if
 744 #endif
 745 
 746  ABI_ALLOCATE(dyewt,(2,3,natom,3,natom))
 747 
 748 ! initialize complex phase factors
 749  do ia = 1, natom
 750    arga = two_pi*( (qphon(1))*xred(1,ia)&
 751 &                 +(qphon(2))*xred(2,ia)&
 752 &                 +(qphon(3))*xred(3,ia) )
 753    exp2piqx(ia) = exp(arga*j_dpc)
 754  end do
 755  ng_expxq = 1000
 756  ABI_ALLOCATE(expx1, (-ng_expxq:ng_expxq, natom))
 757  ABI_ALLOCATE(expx2, (-ng_expxq:ng_expxq, natom))
 758  ABI_ALLOCATE(expx3, (-ng_expxq:ng_expxq, natom))
 759  do ia = 1, natom
 760    do ig1 = -ng_expxq, ng_expxq
 761      expx1(ig1, ia) = exp(ig1*two_pi*xred(1,ia)*j_dpc)
 762      expx2(ig1, ia) = exp(ig1*two_pi*xred(2,ia)*j_dpc)
 763      expx3(ig1, ia) = exp(ig1*two_pi*xred(3,ia)*j_dpc)
 764    end do
 765  end do
 766 
 767  gprimbyacell = gprim
 768  gprimbyacell(:,1) = gprimbyacell(:,1) / acell(1)
 769  gprimbyacell(:,2) = gprimbyacell(:,2) / acell(2)
 770  gprimbyacell(:,3) = gprimbyacell(:,3) / acell(3)
 771 
 772 !compute eta for approximately optimized summations:
 773  direct=rmet(1,1)+rmet(1,2)+rmet(1,3)+rmet(2,1)+&
 774 & rmet(2,2)+rmet(2,3)+rmet(3,1)+rmet(3,2)+rmet(3,3)
 775  recip=gmet(1,1)+gmet(1,2)+gmet(1,3)+gmet(2,1)+&
 776 & gmet(2,2)+gmet(2,3)+gmet(3,1)+gmet(3,2)+gmet(3,3)
 777  eta=pi*100.0_dp/33.0_dp*sqrt(1.69_dp*recip/direct)
 778  inv4eta = one / four / eta
 779 
 780  dyew = zero
 781  dyewt = zero
 782 
 783 !Sum terms over g space:
 784  ng=0
 785  do
 786    ng=ng+1
 787 
 788 ! if needed, update the complex phases for larger G vectors
 789    if (ng > ng_expxq) then
 790      !write(std_out,*)"have to realloc"
 791      ABI_DEALLOCATE(expx1)
 792      ABI_DEALLOCATE(expx2)
 793      ABI_DEALLOCATE(expx3)
 794 
 795      ng_expxq = ng_expxq*2
 796 ! TODO: half of this space is not needed, as it contains the complex conjugate of the other half.
 797 ! present duplication avoids if statements inside the loop, however
 798      ABI_ALLOCATE(expx1, (-ng_expxq:ng_expxq, natom))
 799      ABI_ALLOCATE(expx2, (-ng_expxq:ng_expxq, natom))
 800      ABI_ALLOCATE(expx3, (-ng_expxq:ng_expxq, natom))
 801      do ia = 1, natom
 802        do ig1 = -ng_expxq, ng_expxq
 803          expx1(ig1, ia) = exp(ig1*two_pi*xred(1,ia)*j_dpc)
 804          expx2(ig1, ia) = exp(ig1*two_pi*xred(2,ia)*j_dpc)
 805          expx3(ig1, ia) = exp(ig1*two_pi*xred(3,ia)*j_dpc)
 806        end do
 807      end do
 808    end if
 809 
 810    newg=0
 811    do ig3=-ng,ng
 812      do ig2=-ng,ng
 813        do ig1=-ng,ng
 814          if(abs(ig1)==ng .or. abs(ig2)==ng .or. abs(ig3)==ng .or. ng==1 )then
 815 
 816            gpq(1)=(ig1+qphon(1))*gprimbyacell(1,1)+(ig2+qphon(2))*&
 817 &           gprimbyacell(1,2)+(ig3+qphon(3))*gprimbyacell(1,3)
 818            gpq(2)=(ig1+qphon(1))*gprimbyacell(2,1)+(ig2+qphon(2))*&
 819 &           gprimbyacell(2,2)+(ig3+qphon(3))*gprimbyacell(2,3)
 820            gpq(3)=(ig1+qphon(1))*gprimbyacell(3,1)+(ig2+qphon(2))*&
 821 &           gprimbyacell(3,2)+(ig3+qphon(3))*gprimbyacell(3,3)
 822            gsq=zero
 823            do jj=1,3
 824              do ii=1,3
 825                gpqgpq(ii,jj)=gpq(ii)*gpq(jj)
 826                gsq=gsq+gpqgpq(ii,jj)*dielt(ii,jj)
 827              end do
 828            end do
 829 
 830 !          Skip q=0:
 831            if (gsq<1.0d-20) then
 832              if (sumg0==1) then
 833                write(message,'(a,a,a,a,a)' )&
 834 &               'The phonon wavelength should not be zero :',ch10,&
 835 &               'there are non-analytical terms that cannot be treated.',ch10,&
 836 &               'Action: subtract this wavelength from the input file.'
 837                MSG_ERROR(message)
 838              end if
 839 
 840            else
 841 
 842              arg1=(two_pi**2)*gsq* inv4eta
 843 
 844 !            Larger arg gives 0 contribution:
 845              if (arg1<= -minexparg ) then
 846                newg=1
 847 
 848 !              Here calculate the term
 849                term1=exp(-arg1)/gsq
 850                do jj=1,3
 851                  do ii=1,3
 852                    gpqfac(ii,jj)=gpqgpq(ii,jj)*term1
 853                  end do
 854                end do
 855 
 856 ! MJV: replaced old calls to cos and sin. Checked for 10 tests in v2 that max error is about 6.e-15, usually < 2.e-15
 857                do ia=1,natom
 858                  cosqxred(ia)=real(exp2piqx(ia)*expx1(ig1, ia)*expx2(ig2, ia)*expx3(ig3, ia))
 859                  sinqxred(ia)=aimag(exp2piqx(ia)*expx1(ig1, ia)*expx2(ig2, ia)*expx3(ig3, ia))
 860                end do
 861 
 862 !              First, the diagonal terms
 863                do nu=1,3
 864                  do ia=1,natom
 865                    do mu=nu,3
 866                      dyewt(1,mu,ia,nu,ia)=dyewt(1,mu,ia,nu,ia)+gpqfac(mu,nu)
 867                    end do
 868                  end do
 869                end do
 870 
 871 !              Then, the non-diagonal ones
 872                do ib=2,natom
 873                  do ia=1,ib-1
 874                    c123r=cosqxred(ia)*cosqxred(ib)+sinqxred(ia)*sinqxred(ib)
 875                    c123i=sinqxred(ia)*cosqxred(ib)-cosqxred(ia)*sinqxred(ib)
 876 !                  The most inner loop
 877                    do nu=1,3
 878                      do mu=nu,3
 879                        dyewt(1,mu,ia,nu,ib)=dyewt(1,mu,ia,nu,ib)+gpqfac(mu,nu)*c123r
 880                        dyewt(2,mu,ia,nu,ib)=dyewt(2,mu,ia,nu,ib)+gpqfac(mu,nu)*c123i
 881                      end do
 882                    end do
 883                  end do
 884                end do
 885 
 886              end if ! endif exp() argument is smaller than -minexparg
 887            end if ! Endif g/=0 :
 888          end if ! End triple summation over Gs:
 889        end do
 890      end do
 891    end do
 892 
 893 !  Check if new shell must be calculated
 894    if(newg==0)exit
 895  end do
 896 
 897 !Multiplies by common factor
 898  fact1=4.0_dp*pi/ucvol
 899  do ib=1,natom
 900    do ia=1,ib
 901      do nu=1,3
 902        do mu=nu,3
 903          dyewt(1,mu,ia,nu,ib)=dyewt(1,mu,ia,nu,ib)*fact1
 904          dyewt(2,mu,ia,nu,ib)=dyewt(2,mu,ia,nu,ib)*fact1
 905        end do
 906      end do
 907    end do
 908  end do
 909 
 910  reta=sqrt(eta)
 911  reta3=-eta*reta
 912 
 913 !Calculating the inverse (transpose) of the dielectric tensor
 914  call matr3inv(dielt,invdlt)
 915 !Calculating the determinant of the dielectric tensor
 916  detdlt=dielt(1,1)*dielt(2,2)*dielt(3,3)+dielt(1,3)*dielt(2,1)*&
 917 & dielt(3,2)+dielt(1,2)*dielt(2,3)*dielt(3,1)-dielt(1,3)*&
 918 & dielt(2,2)*dielt(3,1)-dielt(1,1)*dielt(2,3)*dielt(3,2)-&
 919 & dielt(1,2)*dielt(2,1)*dielt(3,3)
 920 
 921  if(detdlt<tol6)then
 922    write(message, '(a,es16.6,11a)' )&
 923 &   'The determinant of the dielectrix matrix, detdlt=',detdlt,' is smaller than 1.0d-6.',ch10,&
 924 &   'The use of the dipole-dipole model for interatomic force constants is not possible.',ch10,&
 925 &   'It is likely that you have not treated the electric field perturbations,',ch10,&
 926 &   'because you not are dealing with an insulator, so that',ch10,&
 927 &   'your dielectric matrix was simply set to zero in the Derivative DataBase.',ch10,&
 928 &   'Action: set the input variable dipdip to 0 .'
 929    MSG_ERROR(message)
 930  end if
 931 
 932  inv_detdlt = one / sqrt(detdlt)
 933  fact3=reta3 * inv_detdlt
 934 
 935 !Preparing the loop on real space
 936  do ia=1,natom
 937    do ii=1,3
 938      xredcar(ii,ia)=(xred(1,ia)*acell(1)*rprim(ii,1)+&
 939 &     xred(2,ia)*acell(2)*rprim(ii,2)+&
 940 &     xred(3,ia)*acell(3)*rprim(ii,3) )*reta
 941    end do
 942  end do
 943  do ia=1,natom
 944    do ii=1,3
 945      xredcax(ii,ia)= invdlt(1,ii)*xredcar(ii,ia)+&
 946 &     invdlt(2,ii)*xredcar(ii,ia)+&
 947 &     invdlt(3,ii)*xredcar(ii,ia)
 948    end do
 949  end do
 950 
 951 !Prepare the evaluation of exp(iq*R)
 952  do ir=-mr,mr
 953    arg1=-two_pi*qphon(1)*ir
 954    arg2=-two_pi*qphon(2)*ir
 955    arg3=-two_pi*qphon(3)*ir
 956    c1r(ir+mr+1)=cos(arg1)
 957    c1i(ir+mr+1)=sin(arg1)
 958    c2r(ir+mr+1)=cos(arg2)
 959    c2i(ir+mr+1)=sin(arg2)
 960    c3r(ir+mr+1)=cos(arg3)
 961    c3i(ir+mr+1)=sin(arg3)
 962  end do
 963 
 964  do nr=1,mr
 965    newr=0
 966 
 967 !  Begin big loop on real space vectors
 968    do ir3=-nr,nr
 969      do ir2=-nr,nr
 970 
 971 !      Here, construct the cosine and sine of q*R for components 2 and 3
 972        c23r = c2r(ir2+mr+1) * c3r(ir3+mr+1) - c2i(ir2+mr+1) * c3i(ir3+mr+1)
 973        c23i = c2i(ir2+mr+1) * c3r(ir3+mr+1) + c2r(ir2+mr+1) * c3i(ir3+mr+1)
 974 
 975 !      Also multiplies by fact3, because it is a rather economical place to do so
 976        c23r=c23r * fact3
 977        c23i=c23i * fact3
 978 
 979        do ir1=-nr,nr
 980          if( abs(ir3)==nr .or. abs(ir2)==nr .or. abs(ir1)==nr .or. nr==1 )then
 981 
 982 !          This is the real part and imaginary part of the phase factor exp(iq*R)
 983            c123r = c1r(ir1+mr+1) * c23r - c1i(ir1+mr+1) * c23i
 984            c123i = c1i(ir1+mr+1) * c23r + c1r(ir1+mr+1) * c23i
 985 
 986            do ii=1,3
 987              ircar(ii)= ( ir1*acell(1)*rprim(ii,1)+&
 988 &             ir2*acell(2)*rprim(ii,2)+&
 989 &             ir3*acell(3)*rprim(ii,3) ) * reta
 990            end do
 991            do ii=1,3
 992              ircax(ii)=   invdlt(1,ii)*ircar(ii)+&
 993 &             invdlt(2,ii)*ircar(ii)+&
 994 &             invdlt(3,ii)*ircar(ii)
 995            end do
 996 
 997 !          Here loops on atoms
 998            do ib=1,natom
 999              do ii=1,3
1000                xredicar(ii)=ircar(ii)-xredcar(ii,ib)
1001                xredicax(ii)=ircax(ii)-xredcax(ii,ib)
1002              end do
1003              do ia=1,ib
1004                do ii=1,3
1005                  rr(ii)=xredicar(ii)+xredcar(ii,ia)
1006                  xx(ii)=xredicax(ii)+xredcax(ii,ia)
1007                end do
1008 
1009                y2=rr(1)*xx(1)+rr(2)*xx(2)+rr(3)*xx(3)
1010 
1011 !              The atoms should not be too far of each other
1012                if (y2 < y2max) then
1013 !                Note: erfc(8) is about 1.1e-29, so dont bother with larger y.
1014 !                Also: exp(-64) is about 1.6e-28, do dont bother with larger y**2 in exp.
1015 
1016 !                Avoid zero denominators in term:
1017                  if (y2 >= y2min) then
1018                    newr=1
1019 ! TODO : find a workaround here - the sqrt, erf and exp functions are slow
1020 !        and for dense q meshes this takes forever
1021 !        could tabulate and spline the full function of yy on a very fine grid
1022 !        and look it up there...
1023 
1024 #ifdef DEV_USESPLINE
1025                    if (y2 > 1.0_dp) then
1026                      ! Spline function
1027                      jspl = 1 + int((y2 - y2min) * stepm1); dd = y2 - y2vals(jspl)
1028                      bb = dd * stepm1
1029                      aa = one - bb
1030                      cc = aa*(aa**2-one) * step2div6
1031                      dd = bb*(bb**2-one) * step2div6
1032 
1033                      term4 = aa*t4spl(jspl,1) + bb*t4spl(jspl+1,1) + cc*t4spl(jspl,2) + dd*t4spl(jspl+1,2)
1034                      term5 = aa*t5spl(jspl,1) + bb*t5spl(jspl+1,1) + cc*t5spl(jspl,2) + dd*t5spl(jspl+1,2)
1035                    else
1036                      ! Evaluate values (the function increases quickly for x --> 0, the spline
1037                      ! is unstable and this leads to SIGFPE.
1038                      yy=sqrt(y2)
1039                      invy=1.0_dp/yy
1040                      invy2=invy**2
1041                      derfc_yy = abi_derfc(yy)
1042                      term2=derfc_yy*invy*invy2
1043                      term3=fact2*exp(-y2)*invy2
1044                      term4=-(term2+term3)
1045                      term5=(3.0_dp*term2+term3*(3.0_dp+2.0_dp*y2))*invy2
1046                    end if
1047 
1048 #else
1049                    yy=sqrt(y2)
1050                    invy=1.0_dp/yy
1051                    invy2=invy**2
1052                    derfc_yy = abi_derfc(yy)
1053                    term2=derfc_yy*invy*invy2
1054                    term3=fact2*exp(-y2)*invy2
1055                    term4=-(term2+term3)
1056                    term5=(3.0_dp*term2+term3*(3.0_dp+2.0_dp*y2))*invy2
1057 #endif
1058                    do nu=1,3
1059                      do mu=nu,3
1060                        dyewt(1,mu,ia,nu,ib)=dyewt(1,mu,ia,nu,ib)+&
1061 &                       c123r*(xx(nu)*xx(mu)*term5+term4*invdlt(nu,mu))
1062                        dyewt(2,mu,ia,nu,ib)=dyewt(2,mu,ia,nu,ib)+&
1063 &                       c123i*(xx(nu)*xx(mu)*term5+term4*invdlt(nu,mu))
1064                      end do
1065                    end do
1066                  else
1067 !                  If zero denominator, the atoms should be identical
1068                    if (ia/=ib)then
1069                      write(message, '(5a,i0,a,i0,a)' )&
1070 &                     'The distance between two atoms seem to vanish.',ch10,&
1071 &                     'This is not allowed.',ch10,&
1072 &                     'Action: check the input for the atoms number',ia,' and',ib,'.'
1073                      MSG_ERROR(message)
1074                    else
1075 !                    This is the correction when the atoms are identical
1076                      do nu=1,3
1077                        do mu=1,3
1078                          dyewt(1,mu,ia,nu,ib)=dyewt(1,mu,ia,nu,ib)+&
1079 &                         fac*reta3*invdlt(nu,mu) * inv_detdlt
1080                        end do
1081                      end do
1082                    end if
1083                  end if ! End the condition for avoiding zero denominators
1084                end if ! End the condition of too large distance between atoms
1085              end do
1086            end do ! End loop over ia and ib :
1087          end if ! End triple loop over real space points:
1088        end do ! ir1
1089      end do ! ir2
1090    end do ! ir3
1091 
1092 !  Check if new shell must be calculated
1093    if(newr==0)exit
1094    if(newr==1 .and. nr==mr) MSG_BUG('mr is too small')
1095  end do
1096 
1097 !Now, symmetrizes
1098  do ib=1,natom-1
1099    do nu=1,3
1100      do ia=ib+1,natom
1101        do mu=nu,3
1102          dyewt(1,mu,ia,nu,ib)=dyewt(1,mu,ib,nu,ia)
1103          dyewt(2,mu,ia,nu,ib)=-dyewt(2,mu,ib,nu,ia)
1104        end do
1105      end do
1106    end do
1107  end do
1108 
1109  do ib=1,natom
1110    do nu=2,3
1111      do ia=1,natom
1112        do mu=1,nu-1
1113          dyewt(1,mu,ia,nu,ib)=dyewt(1,nu,ia,mu,ib)
1114          dyewt(2,mu,ia,nu,ib)=dyewt(2,nu,ia,mu,ib)
1115        end do
1116      end do
1117    end do
1118  end do
1119 
1120 !Tests
1121 !write(std_out,*)' ewald9 : take into account the effective charges '
1122  do ib=1,natom
1123    do nu=1,3
1124      do ia=1,natom
1125        do mu=1,3
1126          do i2=1,2
1127            dyew(i2,mu,ia,nu,ib)=zero
1128            do jj=1,3
1129              do ii=1,3
1130                dyew(i2,mu,ia,nu,ib)=dyew(i2,mu,ia,nu,ib) + &
1131                 zeff(ii,mu,ia)*zeff(jj,nu,ib)*dyewt(i2,ii,ia,jj,ib)
1132              end do
1133            end do
1134          end do
1135        end do
1136      end do
1137    end do
1138  end do
1139 
1140  ABI_DEALLOCATE(expx1)
1141  ABI_DEALLOCATE(expx2)
1142  ABI_DEALLOCATE(expx3)
1143  ABI_DEALLOCATE(dyewt)
1144 
1145 end subroutine ewald9