TABLE OF CONTENTS
ABINIT/dfpt_ewald [ 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