TABLE OF CONTENTS


ABINIT/dfpt_ewald [ Functions ]

[ Top ] [ Functions ]

NAME

 dfpt_ewald

FUNCTION

 Compute ewald contribution to the dynamical matrix, at a given q wavevector.
 Note: the q=0 part should be subtracted, by another call to
 the present routine, with q=0. The present routine correspond
 to the quantity C_bar defined in Eq.(24) or (27) in Phys. Rev. B 55, 10355 (1997). 
 The two calls correspond to Eq.(23) of the same paper.
 If q=0 is asked, sumg0 should be put to 0. Otherwise, it should be put to 1.

COPYRIGHT

 Copyright (C) 1999-2018 ABINIT group (XG)
 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 (length units **-2)
 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
 qphon(3)=phonon wavevector (same system of coordinates as the
          reciprocal lattice vectors)
 rmet(3,3)=metric tensor in real space (length units squared)
 sumg0: if=1, the sum in reciprocal space must include g=0,
   if=0, this contribution must be skipped (q=0 singularity)
 typat(natom)=integer label of each type of atom (1,2,...)
 ucvol=unit cell volume in (whatever length scale units)**3
 xred(3,natom)=relative coords of atoms in unit cell (dimensionless)
 zion(ntypat)=charge on each type of atom (real number)

OUTPUT

 dyew(2,3,natom,3,natom)= Ewald part of the dynamical matrix,
    second energy derivative wrt xred(3,natom), Hartrees.

PARENTS

      respfn

CHILDREN

      free_my_atmtab,get_my_atmtab,timab,xmpi_sum

SOURCE

 50 #if defined HAVE_CONFIG_H
 51 #include "config.h"
 52 #endif
 53 
 54 #include "abi_common.h"
 55 
 56 
 57 subroutine dfpt_ewald(dyew,gmet,my_natom,natom,qphon,rmet,sumg0,typat,ucvol,xred,zion, &
 58 &                 mpi_atmtab,comm_atom ) ! optional arguments (parallelism))
 59 
 60  use defs_basis
 61  use m_profiling_abi
 62  use m_errors
 63  use m_xmpi
 64 
 65  use m_special_funcs,  only : abi_derfc
 66  use m_paral_atom,     only : get_my_atmtab, free_my_atmtab
 67 
 68 !This section has been created automatically by the script Abilint (TD).
 69 !Do not modify the following lines by hand.
 70 #undef ABI_FUNC
 71 #define ABI_FUNC 'dfpt_ewald'
 72  use interfaces_18_timing
 73 !End of the abilint section
 74 
 75  implicit none
 76 
 77 !Arguments -------------------------------
 78 !scalars
 79  integer,intent(in) :: my_natom,natom,sumg0
 80  real(dp),intent(in) :: ucvol
 81 !arrays
 82  integer,intent(in) :: typat(natom)
 83  integer,optional,intent(in) :: comm_atom
 84  integer,optional,target,intent(in) :: mpi_atmtab(:)
 85  real(dp),intent(in) :: gmet(3,3),qphon(3),rmet(3,3),xred(3,natom),zion(*)
 86  real(dp),intent(out) :: dyew(2,3,natom,3,natom)
 87 
 88 !Local variables -------------------------
 89 !nr, ng affect convergence of sums (nr=3,ng=5 is not good enough):
 90 !scalars
 91  integer,parameter :: im=2,ng=10,nr=6,re=1
 92  integer :: ia,ia0,ib,ierr,ig1,ig2,ig3,ii,ir1,ir2,ir3,mu,my_comm_atom,nu
 93  logical :: my_atmtab_allocated,paral_atom
 94  real(dp) :: arg,arga,argb,c1i,c1r,da1,da2,da3,derfc_arg
 95  real(dp) :: direct,dot1,dot2,dot3,dotr1,dotr2,dotr3
 96  real(dp) :: eta,fac,gdot12,gdot13,gdot23,gsq,gsum,norm1
 97  real(dp) :: r1,r2,r3,rdot12,rdot13,rdot23,recip,reta
 98  real(dp) :: reta3m,rmagn,rsq,term,term1,term2
 99  real(dp) :: term3
100  character(len=500) :: message
101 !arrays
102  real(dp) :: tsec(2)
103  integer,pointer :: my_atmtab(:)
104  real(dp) :: gpq(3),rq(3)
105 
106 ! *************************************************************************
107 
108 !Set up parallelism over atoms
109  paral_atom=(present(comm_atom).and.(my_natom/=natom))
110  nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
111  my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
112  call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
113 
114 !Compute eta for approximately optimized summations:
115  direct=rmet(1,1)+rmet(1,2)+rmet(1,3)+rmet(2,1)+&
116 & 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 & gmet(2,2)+gmet(2,3)+gmet(3,1)+gmet(3,2)+gmet(3,3)
119  eta=pi*(dble(ng)/dble(nr))*sqrt(1.69_dp*recip/direct)
120 
121 !Test Ewald s summation
122 !eta=1.2_dp*eta
123 
124 !Sum terms over g space:
125  fac=pi**2/eta
126  gsum=zero
127  da1=zero
128  da2=zero
129  da3=zero
130  dyew(:,:,:,:,:)=zero
131  do ig3=-ng,ng
132    do ig2=-ng,ng
133      do ig1=-ng,ng
134        gpq(1)=dble(ig1)+qphon(1)
135        gpq(2)=dble(ig2)+qphon(2)
136        gpq(3)=dble(ig3)+qphon(3)
137        gdot12=gmet(2,1)*gpq(1)*gpq(2)
138        gdot13=gmet(3,1)*gpq(1)*gpq(3)
139        gdot23=gmet(3,2)*gpq(2)*gpq(3)
140        dot1=gmet(1,1)*gpq(1)**2+gdot12+gdot13
141        dot2=gmet(2,2)*gpq(2)**2+gdot12+gdot23
142        dot3=gmet(3,3)*gpq(3)**2+gdot13+gdot23
143        gsq=dot1+dot2+dot3
144 !      Skip q=0:
145        if (gsq<1.0d-20) then
146          if (sumg0==1) then
147            write(message,'(5a)')&
148 &           'The phonon wavelength should not be zero : ',ch10,&
149 &           'there are non-analytical terms that the code cannot handle.',ch10,&
150 &           'Action : subtract this wavelength from the input.'
151            MSG_ERROR(message)
152          end if
153        else
154          arg=fac*gsq
155 !        Larger arg gives 0 contribution:
156          if (arg <= 80._dp) then
157            term=exp(-arg)/gsq
158            do ia0=1,my_natom
159              ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
160              arga=two_pi*(gpq(1)*xred(1,ia)+gpq(2)*xred(2,ia)+gpq(3)*xred(3,ia))
161              do ib=1,ia
162                argb=two_pi*(gpq(1)*xred(1,ib)+gpq(2)*xred(2,ib)+gpq(3)*xred(3,ib))
163                arg=arga-argb
164                c1r=cos(arg)*term
165                c1i=sin(arg)*term
166 
167                do mu=1,3
168                  do nu=1,mu
169                    dyew(re,mu,ia,nu,ib)=dyew(re,mu,ia,nu,ib)+gpq(mu)*gpq(nu)*c1r
170                    dyew(im,mu,ia,nu,ib)=dyew(im,mu,ia,nu,ib)+gpq(mu)*gpq(nu)*c1i
171                  end do
172                end do
173 
174              end do
175            end do
176          end if
177 !        Endif g/=0 :
178        end if
179 !      End triple loop over G s:
180      end do
181    end do
182  end do
183 
184 !End G summation by accounting for some common factors.
185 !(for the charges:see end of routine)
186  norm1=4.0_dp*pi/ucvol
187  do ia0=1,my_natom
188    ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
189    do ib=1,ia
190      do mu=1,3
191        do nu=1,mu
192          dyew(:,mu,ia,nu,ib)=dyew(:,mu,ia,nu,ib)*norm1
193        end do
194      end do
195    end do
196  end do
197 
198 
199 !Do sums over real space:
200  reta=sqrt(eta)
201  reta3m=-eta*reta
202  fac=4._dp/3.0_dp/sqrt(pi)
203  do ir3=-nr,nr
204    do ir2=-nr,nr
205      do ir1=-nr,nr
206        arg=-two_pi*(qphon(1)*ir1+qphon(2)*ir2+qphon(3)*ir3)
207        c1r=cos(arg)*reta3m
208        c1i=sin(arg)*reta3m
209        do ia0=1,my_natom 
210          ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
211          do ib=1,ia
212            r1=dble(ir1)+xred(1,ia)-xred(1,ib)
213            r2=dble(ir2)+xred(2,ia)-xred(2,ib)
214            r3=dble(ir3)+xred(3,ia)-xred(3,ib)
215            rdot12=rmet(2,1)*r1*r2
216            rdot13=rmet(3,1)*r1*r3
217            rdot23=rmet(3,2)*r2*r3
218            dotr1=rmet(1,1)*r1**2+rdot12+rdot13
219            dotr2=rmet(2,2)*r2**2+rdot12+rdot23
220            dotr3=rmet(3,3)*r3**2+rdot13+rdot23
221            rsq=dotr1+dotr2+dotr3
222            rmagn=sqrt(rsq)
223 !          Avoid zero denominators in term :
224            if (rmagn>=1.0d-12) then
225              arg=reta*rmagn
226              term=zero
227              if (arg<8.0_dp) then
228 !              Note: erfc(8) is about 1.1e-29,
229 !              so don t bother with larger arg.
230 !              Also: exp(-64) is about 1.6e-28,
231 !              so don t bother with larger arg**2 in exp.
232                derfc_arg = abi_derfc(arg)
233                term=derfc_arg/arg**3
234                term1=2.0_dp/sqrt(pi)*exp(-arg**2)/arg**2
235                term2=-(term+term1)
236                term3=(3*term+term1*(3.0_dp+2.0_dp*arg**2))/rsq
237                rq(1)=rmet(1,1)*r1+rmet(1,2)*r2+rmet(1,3)*r3
238                rq(2)=rmet(2,1)*r1+rmet(2,2)*r2+rmet(2,3)*r3
239                rq(3)=rmet(3,1)*r1+rmet(3,2)*r2+rmet(3,3)*r3
240                do mu=1,3
241                  do nu=1,mu
242                    dyew(re,mu,ia,nu,ib)=dyew(re,mu,ia,nu,ib)+&
243 &                   c1r*(rq(mu)*rq(nu)*term3+rmet(mu,nu)*term2)
244                    dyew(im,mu,ia,nu,ib)=dyew(im,mu,ia,nu,ib)+&
245 &                   c1i*(rq(mu)*rq(nu)*term3+rmet(mu,nu)*term2)
246                  end do
247                end do
248              end if
249            else
250              if (ia/=ib)then
251                write(message,'(a,a,a,a,a,i5,a,i5,a)')&
252 &               'The distance between two atoms vanishes.',ch10,&
253 &               'This is not allowed.',ch10,&
254 &               'Action: check the input for the atoms number',ia,' and',ib,'.'
255                MSG_ERROR(message)
256              else
257                do mu=1,3
258                  do nu=1,mu
259                    dyew(re,mu,ia,nu,ib)=dyew(re,mu,ia,nu,ib)+&
260 &                   fac*reta3m*rmet(mu,nu)
261                  end do
262                end do
263              end if
264            end if
265 
266          end do ! End loop over ib:
267        end do ! End loop over ia:
268      end do ! End triple loop over real space points:
269    end do
270  end do
271 
272 !Take account of the charges
273 !write(std_out,*)' '
274  do ia0=1,my_natom     
275    ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
276    do ib=1,ia
277      do mu=1,3
278        do nu=1,mu
279          do ii=1,2
280 !          write(std_out,*)dyew(ii,mu,ia,nu,ib)
281            dyew(ii,mu,ia,nu,ib)=dyew(ii,mu,ia,nu,ib)*&
282 &           zion(typat(ia))*zion(typat(ib))
283          end do
284        end do
285      end do
286    end do
287  end do
288 
289 !Symmetrize with respect to the directions
290  do ia0=1,my_natom
291    ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
292    do ib=1,ia
293      do mu=1,3
294        do nu=1,mu
295          dyew(re,nu,ia,mu,ib)=dyew(re,mu,ia,nu,ib)
296          dyew(im,nu,ia,mu,ib)=dyew(im,mu,ia,nu,ib)
297        end do
298      end do
299    end do
300  end do
301 
302 !In case of parallelism over atoms: communicate
303  if (paral_atom) then
304    call timab(48,1,tsec)
305    call xmpi_sum(dyew,my_comm_atom,ierr)
306    call timab(48,2,tsec)
307  end if
308 
309 !Fill the upper part of the matrix, with the hermitian conjugate
310  do ia=1,natom
311    do ib=1,ia
312      do nu=1,3
313        do mu=1,3
314          dyew(re,mu,ib,nu,ia)=dyew(re,mu,ia,nu,ib)
315          dyew(im,mu,ib,nu,ia)=-dyew(im,mu,ia,nu,ib)
316        end do
317      end do
318    end do
319  end do
320 
321 !Destroy atom table used for parallelism
322  call free_my_atmtab(my_atmtab,my_atmtab_allocated)
323 
324 end subroutine dfpt_ewald