TABLE OF CONTENTS


ABINIT/elt_ewald [ Functions ]

[ Top ] [ Functions ]

NAME

 elt_ewald

FUNCTION

 Compute 2nd derivatives of Ewald energy wrt strain for frozen wavefunction
 contributions to elastic tensor

COPYRIGHT

 Copyright (C) 1998-2018 ABINIT group (DRH, DCA, XG, 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 .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

INPUTS

 gmet(3,3)=metric tensor in reciprocal space (bohr^-2)
 gprimd(3,3)=dimensional primitive translations for reciprocal space (bohr^-1)
 mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
 comm_atom=--optional-- MPI communicator over atoms
 my_natom=number of atoms treated by current processor
 natom=number of atoms in unit cell
 ntypat=numbe of type of atoms
 rmet(3,3)=metric tensor in real space (bohr^2)
 rprimd(3,3)=dimensional primitive translation vectors (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

 elteew(6+3*natom,6)=2nd derivatives of Ewald energy wrt strain

PARENTS

      respfn

CHILDREN

      free_my_atmtab,get_my_atmtab,timab,wrtout,xmpi_sum

SOURCE

 44 #if defined HAVE_CONFIG_H
 45 #include "config.h"
 46 #endif
 47 
 48 #include "abi_common.h"
 49 
 50 
 51 subroutine elt_ewald(elteew,gmet,gprimd,my_natom,natom,ntypat,rmet,rprimd,&
 52 &                 typat,ucvol,xred,zion, &
 53 &                 mpi_atmtab,comm_atom) ! optional arguments (parallelism)
 54 
 55  use defs_basis
 56  use m_profiling_abi
 57  use m_xmpi
 58 
 59  use m_special_funcs,  only : abi_derfc
 60  use m_paral_atom,     only : get_my_atmtab, free_my_atmtab
 61 
 62 !This section has been created automatically by the script Abilint (TD).
 63 !Do not modify the following lines by hand.
 64 #undef ABI_FUNC
 65 #define ABI_FUNC 'elt_ewald'
 66  use interfaces_14_hidewrite
 67  use interfaces_18_timing
 68 !End of the abilint section
 69 
 70  implicit none
 71 
 72 !Arguments ------------------------------------
 73 !scalars
 74  integer,intent(in) :: my_natom,natom,ntypat
 75  real(dp),intent(in) :: ucvol
 76 !arrays
 77  integer,intent(in) :: typat(natom)
 78  integer,optional,intent(in) :: comm_atom
 79  integer,optional,target,intent(in) :: mpi_atmtab(:)
 80  real(dp),intent(in) :: gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3)
 81  real(dp),intent(in) :: xred(3,natom),zion(ntypat)
 82  real(dp),intent(out) :: elteew(6+3*natom,6)
 83 
 84 !Local variables-------------------------------
 85 !scalars
 86  integer :: ia,ia0,ib,ierr,ig1,ig2,ig3,ir1,ir2,ir3,is1,is2,jj,js,ka,kb,kd,kg,my_comm_atom,newg,newr,ng,nr
 87  logical :: my_atmtab_allocated,paral_atom
 88  real(dp) :: arg,ch,chsq,cos_term,d2derfc,d2gss,d2r,d2rs,dderfc,derfc_arg
 89  real(dp) :: dgss1,dgss2,direct,dr1,dr2,drs1,drs2,eew,eta,fac,fraca1,fraca2
 90  real(dp) :: fraca3,fracb1,fracb2,fracb3,gsq,gsum,r1,r2,r3,recip,reta
 91  real(dp) :: rmagn,rsq,sin_term,sumg,summi,summr,sumr,t1,term
 92  character(len=500) :: message
 93 !arrays
 94  integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/)
 95  integer,pointer :: my_atmtab(:)
 96  real(dp) :: d2gm(3,3,6,6),d2ris(3),d2rm(3,3,6,6),dgm(3,3,6),dris(3),drm(3,3,6)
 97  real(dp) :: t2(3),ts2(3),tsec(2),tt(3)
 98  real(dp),allocatable :: d2sumg(:,:),d2sumr(:,:),drhoisi(:,:),drhoisr(:,:)
 99  real(dp),allocatable :: mpibuf(:)
100 
101 ! *************************************************************************
102 
103 !DEBUG
104 !write(std_out,*)' elt_ewald : enter '
105 !stop
106 !ENDDEBUG
107 
108 !Compute 1st and 2nd derivatives of metric tensor wrt all strain components
109 !and store for use in inner loop below.
110 
111 !Set up parallelism over atoms
112  paral_atom=(present(comm_atom).and.(my_natom/=natom))
113  nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
114  my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
115  call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
116 
117 !Loop over 2nd strain index
118  do is2=1,6
119    kg=idx(2*is2-1);kd=idx(2*is2)
120    do jj = 1,3
121      drm(:,jj,is2)=rprimd(kg,:)*rprimd(kd,jj)+rprimd(kd,:)*rprimd(kg,jj)
122      dgm(:,jj,is2)=-(gprimd(kg,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(kg,jj))
123    end do
124 
125 !  Loop over 1st strain index, upper triangle only
126    do is1=1,is2
127 
128      ka=idx(2*is1-1);kb=idx(2*is1)
129      d2rm(:,:,is1,is2)=zero
130      d2gm(:,:,is1,is2)=zero
131      do jj = 1,3
132        if(ka==kg) d2rm(:,jj,is1,is2)=d2rm(:,jj,is1,is2)&
133 &       +rprimd(kb,:)*rprimd(kd,jj)+rprimd(kd,:)*rprimd(kb,jj)
134        if(ka==kd) d2rm(:,jj,is1,is2)=d2rm(:,jj,is1,is2)&
135 &       +rprimd(kb,:)*rprimd(kg,jj)+rprimd(kg,:)*rprimd(kb,jj)
136        if(kb==kg) d2rm(:,jj,is1,is2)=d2rm(:,jj,is1,is2)&
137 &       +rprimd(ka,:)*rprimd(kd,jj)+rprimd(kd,:)*rprimd(ka,jj)
138        if(kb==kd) d2rm(:,jj,is1,is2)=d2rm(:,jj,is1,is2)&
139 &       +rprimd(ka,:)*rprimd(kg,jj)+rprimd(kg,:)*rprimd(ka,jj)
140 
141        if(ka==kg) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
142 &       +gprimd(kb,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(kb,jj)
143        if(ka==kd) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
144 &       +gprimd(kb,:)*gprimd(kg,jj)+gprimd(kg,:)*gprimd(kb,jj)
145        if(kb==kg) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
146 &       +gprimd(ka,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(ka,jj)
147        if(kb==kd) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
148 &       +gprimd(ka,:)*gprimd(kg,jj)+gprimd(kg,:)*gprimd(ka,jj)
149      end do
150    end do !is1
151  end do !is2
152 
153 !Add up total charge and sum of $charge^2$ in cell
154  chsq=zero
155  ch=zero
156  do ia=1,natom
157    ch=ch+zion(typat(ia))
158    chsq=chsq+zion(typat(ia))**2
159  end do
160 
161 !Compute eta, the Ewald summation convergence parameter,
162 !for approximately optimized summations:
163  direct=rmet(1,1)+rmet(1,2)+rmet(1,3)+rmet(2,1)+&
164 & rmet(2,2)+rmet(2,3)+rmet(3,1)+rmet(3,2)+rmet(3,3)
165  recip=gmet(1,1)+gmet(1,2)+gmet(1,3)+gmet(2,1)+&
166 & gmet(2,2)+gmet(2,3)+gmet(3,1)+gmet(3,2)+gmet(3,3)
167 !Here, a bias is introduced, because G-space summation scales
168 !better than r space summation ! Note : debugging is the most
169 !easier at fixed eta.
170  eta=pi*200._dp/33.0_dp*sqrt(1.69_dp*recip/direct)
171 
172 !Conduct reciprocal space summations
173  fac=pi**2/eta ; gsum=zero
174  ABI_ALLOCATE(d2sumg,(6+3*natom,6))
175  ABI_ALLOCATE(drhoisr,(3,natom))
176  ABI_ALLOCATE(drhoisi,(3,natom))
177  d2sumg(:,:)=zero
178 
179 !Sum over G space, done shell after shell until all
180 !contributions are too small.
181  ng=0
182  do
183    ng=ng+1
184    newg=0
185 
186    do ig3=-ng,ng
187      do ig2=-ng,ng
188        do ig1=-ng,ng
189 
190 !        Exclude shells previously summed over
191          if(abs(ig1)==ng .or. abs(ig2)==ng .or. abs(ig3)==ng&
192 &         .or. ng==1 ) then
193 
194 !          gsq is G dot G = |G|^2
195            gsq=gmet(1,1)*dble(ig1*ig1)+gmet(2,2)*dble(ig2*ig2)+&
196 &           gmet(3,3)*dble(ig3*ig3)+2._dp*(gmet(2,1)*dble(ig1*ig2)+&
197 &           gmet(3,1)*dble(ig1*ig3)+gmet(3,2)*dble(ig3*ig2))
198 
199 !          Skip g=0:
200            if (gsq>1.0d-20) then
201              arg=fac*gsq
202 
203 !            Larger arg gives 0 contribution because of exp(-arg)
204              if (arg <= 80._dp) then
205 !              When any term contributes then include next shell
206                newg=1
207                term=exp(-arg)/gsq
208                summr = zero
209                summi = zero
210 !              Note that if reduced atomic coordinates xred drift outside
211 !              of unit cell (outside [0,1)) it is irrelevant in the following
212 !              term, which only computes a phase.
213                do ia=1,natom
214                  arg=two_pi*(ig1*xred(1,ia)+ig2*xred(2,ia)+ig3*xred(3,ia))
215 !                Sum real and imaginary parts (avoid complex variables)
216                  cos_term=cos(arg)
217                  sin_term=sin(arg)
218                  summr=summr+zion(typat(ia))*cos_term
219                  summi=summi+zion(typat(ia))*sin_term
220                  drhoisr(1,ia)=-two_pi*zion(typat(ia))*sin_term*dble(ig1)
221                  drhoisi(1,ia)= two_pi*zion(typat(ia))*cos_term*dble(ig1)
222                  drhoisr(2,ia)=-two_pi*zion(typat(ia))*sin_term*dble(ig2)
223                  drhoisi(2,ia)= two_pi*zion(typat(ia))*cos_term*dble(ig2)
224                  drhoisr(3,ia)=-two_pi*zion(typat(ia))*sin_term*dble(ig3)
225                  drhoisi(3,ia)= two_pi*zion(typat(ia))*cos_term*dble(ig3)
226                end do
227 
228 !              The following two checks avoid an annoying
229 !              underflow error message
230                if (abs(summr)<1.d-16) summr=zero
231                if (abs(summi)<1.d-16) summi=zero
232 
233 !              The product of term and summr**2 or summi**2 below
234 !              can underflow if not for checks above
235                t1=term*(summr*summr+summi*summi)
236                gsum=gsum+t1
237 !              Loop over 2nd strain index
238                do is2=1,6
239                  dgss2=dgm(1,1,is2)*dble(ig1*ig1)+dgm(2,2,is2)*dble(ig2*ig2)+&
240 &                 dgm(3,3,is2)*dble(ig3*ig3)+2._dp*(dgm(2,1,is2)*dble(ig1*ig2)+&
241 &                 dgm(3,1,is2)*dble(ig1*ig3)+dgm(3,2,is2)*dble(ig3*ig2))
242 !                Loop over 1st strain index, upper triangle only
243                  do is1=1,is2
244                    dgss1=dgm(1,1,is1)*dble(ig1*ig1)+dgm(2,2,is1)*dble(ig2*ig2)+&
245 &                   dgm(3,3,is1)*dble(ig3*ig3)+2._dp*(dgm(2,1,is1)*dble(ig1*ig2)+&
246 &                   dgm(3,1,is1)*dble(ig1*ig3)+dgm(3,2,is1)*dble(ig3*ig2))
247 
248                    d2gss=d2gm(1,1,is1,is2)*dble(ig1*ig1)+&
249 &                   d2gm(2,2,is1,is2)*dble(ig2*ig2)+&
250 &                   d2gm(3,3,is1,is2)*dble(ig3*ig3)+2._dp*&
251 &                   (d2gm(2,1,is1,is2)*dble(ig1*ig2)+&
252 &                   d2gm(3,1,is1,is2)*dble(ig1*ig3)+&
253 &                   d2gm(3,2,is1,is2)*dble(ig3*ig2))
254 
255                    d2sumg(is1,is2)=d2sumg(is1,is2)+&
256 &                   t1*((fac**2 + 2.0_dp*fac/gsq + 2.0_dp/(gsq**2))*dgss1*dgss2 -&
257 &                   0.5_dp*(fac + 1.0_dp/gsq)*d2gss)
258                    if(is1<=3) d2sumg(is1,is2)=d2sumg(is1,is2)+&
259 &                   t1*(fac + 1.0_dp/gsq)*dgss2
260                    if(is2<=3) d2sumg(is1,is2)=d2sumg(is1,is2)+&
261 &                   t1*(fac + 1.0_dp/gsq)*dgss1
262                    if(is1<=3 .and. is2<=3) d2sumg(is1,is2)=d2sumg(is1,is2)+t1
263 
264                  end do !is1
265 
266 !                Internal strain contributions
267                  do ia=1,natom
268                    js=7+3*(ia-1)
269                    t2(:)=2.0_dp*term*(summr*drhoisr(:,ia)+summi*drhoisi(:,ia))
270                    d2sumg(js:js+2,is2)=d2sumg(js:js+2,is2)-&
271 &                   (fac + 1.0_dp/gsq)*dgss2*t2(:)
272                    if(is2<=3) d2sumg(js:js+2,is2)=d2sumg(js:js+2,is2)-t2(:)
273                  end do
274                end do !is2
275 
276              end if ! End condition of not larger than 80.0
277            end if ! End skip g=0
278          end if ! End triple loop over G s and associated new shell condition
279        end do
280      end do
281    end do
282 
283 !  Check if new shell must be calculated
284    if (newg==0) exit
285  end do !  End the loop on ng (new shells). Note that there is one exit from this loop.
286 
287  sumg=gsum/(two_pi*ucvol)
288  d2sumg(:,:)=d2sumg(:,:)/(two_pi*ucvol)
289 
290  ABI_DEALLOCATE(drhoisr)
291  ABI_DEALLOCATE(drhoisi)
292 !Stress tensor is now computed elsewhere (ewald2) hence do not need
293 !length scale gradients (used to compute them here).
294 
295 !Conduct real space summations
296  reta=sqrt(eta)
297  fac=2._dp*sqrt(eta/pi)
298  ABI_ALLOCATE(d2sumr,(6+3*natom,6))
299  sumr=zero;d2sumr(:,:)=zero
300 
301 !In the following a summation is being conducted over all
302 !unit cells (ir1, ir2, ir3) so it is appropriate to map all
303 !reduced coordinates xred back into [0,1).
304 !
305 !Loop on shells in r-space as was done in g-space
306  nr=0
307  do
308    nr=nr+1
309    newr=0
310 !  
311    do ir3=-nr,nr
312      do ir2=-nr,nr
313        do ir1=-nr,nr
314          if( abs(ir3)==nr .or. abs(ir2)==nr .or. abs(ir1)==nr .or. nr==1 )then
315 
316            do ia0=1,my_natom
317              ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
318              js=7+3*(ia-1)
319 !            Map reduced coordinate xred(mu,ia) into [0,1)
320              fraca1=xred(1,ia)-aint(xred(1,ia))+0.5_dp-sign(0.5_dp,xred(1,ia))
321              fraca2=xred(2,ia)-aint(xred(2,ia))+0.5_dp-sign(0.5_dp,xred(2,ia))
322              fraca3=xred(3,ia)-aint(xred(3,ia))+0.5_dp-sign(0.5_dp,xred(3,ia))
323              do ib=1,natom
324                fracb1=xred(1,ib)-aint(xred(1,ib))+0.5_dp-sign(0.5_dp,xred(1,ib))
325                fracb2=xred(2,ib)-aint(xred(2,ib))+0.5_dp-sign(0.5_dp,xred(2,ib))
326                fracb3=xred(3,ib)-aint(xred(3,ib))+0.5_dp-sign(0.5_dp,xred(3,ib))
327                r1=dble(ir1)+fracb1-fraca1
328                r2=dble(ir2)+fracb2-fraca2
329                r3=dble(ir3)+fracb3-fraca3
330                rsq=rmet(1,1)*r1*r1+rmet(2,2)*r2*r2+rmet(3,3)*r3*r3+&
331 &               2.0_dp*(rmet(2,1)*r2*r1+rmet(3,2)*r3*r2+rmet(3,1)*r1*r3)
332 
333 !              Avoid zero denominators in 'term':
334                if (rsq>=1.0d-24) then
335 
336 !                Note: erfc(8) is about 1.1e-29,
337 !                so do not bother with larger arg.
338 !                Also: exp(-64) is about 1.6e-28,
339 !                so do not bother with larger arg**2 in exp.
340                  term=zero
341                  if (eta*rsq<64.0_dp) then
342                    newr=1
343                    rmagn=sqrt(rsq)
344                    arg=reta*rmagn
345 !                  derfc computes the complementary error function
346 !                  dderfc is the derivative of the complementary error function
347 !                  d2derfc is the 2nd derivative of the complementary error function
348                    dderfc=-fac*exp(-eta*rsq)
349                    d2derfc=-2._dp*eta*rmagn*dderfc
350                    derfc_arg = abi_derfc(arg)
351                    term=derfc_arg/rmagn
352                    sumr=sumr+zion(typat(ia))*zion(typat(ib))*term
353                    tt(:)=rmet(:,1)*r1+rmet(:,2)*r2+rmet(:,3)*r3
354                    dris(:)=tt(:)/rmagn
355 !                  Loop over 2nd strain index
356                    do is2=1,6
357                      drs2=drm(1,1,is2)*r1*r1+drm(2,2,is2)*r2*r2+&
358 &                     drm(3,3,is2)*r3*r3+&
359 &                     2.0_dp*(drm(2,1,is2)*r2*r1+drm(3,2,is2)*r3*r2+&
360 &                     drm(3,1,is2)*r1*r3)
361                      dr2=0.5_dp*drs2/rmagn
362 !                    Loop over 1st strain index, upper triangle only
363                      do is1=1,is2
364                        drs1=drm(1,1,is1)*r1*r1+drm(2,2,is1)*r2*r2+&
365 &                       drm(3,3,is1)*r3*r3+&
366 &                       2.0_dp*(drm(2,1,is1)*r2*r1+drm(3,2,is1)*r3*r2+&
367 &                       drm(3,1,is1)*r1*r3)
368                        dr1=0.5_dp*drs1/rmagn
369                        d2rs=d2rm(1,1,is1,is2)*r1*r1+d2rm(2,2,is1,is2)*r2*r2+&
370 &                       d2rm(3,3,is1,is2)*r3*r3+&
371 &                       2.0_dp*(d2rm(2,1,is1,is2)*r2*r1+d2rm(3,2,is1,is2)*r3*r2+&
372 &                       d2rm(3,1,is1,is2)*r1*r3)
373                        d2r=(0.25_dp*d2rs-dr1*dr2)/rmagn
374                        d2sumr(is1,is2)=d2sumr(is1,is2)+&
375 &                       zion(typat(ia))*zion(typat(ib))*&
376 &                       ((d2derfc-2.0_dp*dderfc/rmagn+2.0_dp*derfc_arg/rsq)*dr1*dr2+&
377 &                       (dderfc-derfc_arg/rmagn)*d2r)/rmagn
378                      end do !is1
379 !                    Internal strain contribution
380                      ts2(:)=drm(:,1,is2)*r1+drm(:,2,is2)*r2+drm(:,3,is2)*r3
381                      d2ris(:)=ts2(:)/rmagn-0.5_dp*drs2*tt(:)/(rsq*rmagn)
382 
383                      d2sumr(js:js+2,is2)=d2sumr(js:js+2,is2)-&
384 &                     2.0_dp*zion(typat(ia))*zion(typat(ib))*&
385 &                     ((d2derfc-2.0_dp*dderfc/rmagn+2.0_dp*derfc_arg/rsq)*dr2*dris(:)+&
386 &                     (dderfc-derfc_arg/rmagn)*d2ris(:))/rmagn
387                    end do !is2
388                  end if
389                end if ! End avoid zero denominators in'term'
390 
391              end do ! end loop over ib:
392            end do ! end loop over ia:
393          end if ! end triple loop over real space points and associated condition of new shell
394        end do
395      end do
396    end do
397 
398 !  Check if new shell must be calculated
399    if(newr==0) exit
400  end do !  End loop on nr (new shells). Note that there is an exit within the loop
401 
402 !In case of parallelism over atoms: communicate
403  if (paral_atom) then
404    call timab(48,1,tsec)
405    ABI_ALLOCATE(mpibuf,((6+3*natom)*6+1)) 
406    mpibuf(1:(6+3*natom)*6)=reshape(d2sumr(:,:),shape=(/((6+3*natom)*6)/))
407    mpibuf((6+3*natom)*6+1)=sumr
408    call xmpi_sum(mpibuf,my_comm_atom,ierr)
409    sumr=mpibuf((6+3*natom)*6+1)
410    d2sumr(:,:)=reshape(mpibuf(1:(6+3*natom)*6),shape=(/(6+3*natom),6/))
411    ABI_DEALLOCATE(mpibuf) 
412    call timab(48,2,tsec)
413  end if
414 
415  sumr=0.5_dp*sumr
416  d2sumr(:,:)=0.5_dp*d2sumr(:,:)
417  fac=pi*ch**2/(2.0_dp*eta*ucvol)
418 
419 !Finally assemble Ewald energy, eew
420  eew=sumg+sumr-chsq*reta/sqrt(pi)-fac
421 
422  elteew(:,:)=d2sumg(:,:)+d2sumr(:,:)
423 
424 !Additional term for all strains diagonal (from "fac" term in eew)
425  elteew(1:3,1:3)=elteew(1:3,1:3)-fac
426 
427 !Fill in lower triangle
428  do is2=2,6
429    do is1=1,is2-1
430      elteew(is2,is1)=elteew(is1,is2)
431    end do
432  end do
433 
434  ABI_DEALLOCATE(d2sumg)
435  ABI_DEALLOCATE(d2sumr)
436 
437 !Output the final values of ng and nr
438  write(message, '(a,i4,a,i4)' )' elt_ewald : nr and ng are ',nr,' and ',ng
439  call wrtout(std_out,message,'COLL')
440 
441 !Destroy atom table used for parallelism
442  call free_my_atmtab(my_atmtab,my_atmtab_allocated)
443 
444 end subroutine elt_ewald