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