TABLE OF CONTENTS
ABINIT/m_ewald [ Modules ]
NAME
m_ewald
FUNCTION
This module gathers routines to compute the Ewald energy and its derivatives
COPYRIGHT
Copyright (C) 2014-2024 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 .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_ewald 23 24 use defs_basis 25 use m_abicore 26 use m_errors 27 use m_splines 28 use m_time 29 use m_xmpi 30 31 use m_gtermcutoff, only : termcutoff 32 use m_special_funcs, only : abi_derfc 33 use m_symtk, only : matr3inv 34 35 implicit none 36 37 private 38 39 public :: ewald ! Compute Ewald energy and derivatives with respect to xred 40 public :: ewald2 ! Derivative of the Ewald energy with respect to strain. 41 public :: ewald9 ! Compute ewald contribution to the dynamical matrix, at a given 42 ! q wavevector, including anisotropic dielectric tensor and effective charges 43 44 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.
SOURCE
72 subroutine ewald(eew,gmet,grewtn,gsqcut,icutcoul,natom,ngfft,nkpt,ntypat,rcut,rmet,rprimd,typat,ucvol,vcutgeo,xred,zion) 73 74 !Arguments ------------------------------------ 75 !scalars 76 integer,intent(in) :: icutcoul,natom,nkpt,ntypat 77 real(dp),intent(in) :: gsqcut,rcut,ucvol 78 real(dp),intent(out) :: eew 79 !arrays 80 integer,intent(in) :: ngfft(18),typat(natom) 81 real(dp),intent(in) :: gmet(3,3),rmet(3,3),rprimd(3,3),xred(3,natom),vcutgeo(3),zion(ntypat) 82 real(dp),intent(out) :: grewtn(3,natom) 83 84 !Local variables------------------------------- 85 !scalars 86 integer :: ia,ib,ig1,ig2,ig3,ig23,ii,ir1,ir2,ir3,newg,newr,ng,nr 87 real(dp) :: arg,c1i,ch,chsq,derfc_arg,direct,drdta1,drdta2,drdta3,eta,fac 88 real(dp) :: fraca1,fraca2,fraca3,fracb1,fracb2,fracb3,gsq,gsum,phi,phr,r1 89 real(dp) :: minexparg 90 real(dp) :: r1a1d,r2,r2a2d,r3,r3a3d,recip,reta,rmagn,rsq,sumg,summi,summr,sumr 91 real(dp) :: t1,term ,zcut !, gcart_para, gcart_perp 92 !character(len=500) :: msg 93 !arrays 94 real(dp),allocatable :: gcutoff(:) 95 96 ! ************************************************************************* 97 98 !This is the minimum argument of an exponential, with some safety 99 minexparg=log(tiny(0._dp))+five 100 101 !Add up total charge and sum of $charge^2$ in cell 102 103 chsq=0._dp 104 ch=0._dp 105 do ia=1,natom 106 ch=ch+zion(typat(ia)) 107 chsq=chsq+zion(typat(ia))**2 108 end do 109 110 !Compute eta, the Ewald summation convergence parameter, 111 !for approximately optimized summations: 112 direct=rmet(1,1)+rmet(1,2)+rmet(1,3)+rmet(2,1)+& 113 & rmet(2,2)+rmet(2,3)+rmet(3,1)+rmet(3,2)+rmet(3,3) 114 recip=gmet(1,1)+gmet(1,2)+gmet(1,3)+gmet(2,1)+& 115 & gmet(2,2)+gmet(2,3)+gmet(3,1)+gmet(3,2)+gmet(3,3) 116 !A bias is introduced, because G-space summation scales 117 !better than r space summation ! Note : debugging is the most 118 !easier at fixed eta. 119 zcut=SQRT(DOT_PRODUCT(rprimd(:,3),rprimd(:,3)))/2.0_dp 120 if(icutcoul.eq.1) then 121 eta=SQRT(16.0_dp/SQRT(DOT_PRODUCT(rprimd(:,1),rprimd(:,1)))) 122 ! else if (icutcoul.eq.2) then 123 ! zcut=SQRT(DOT_PRODUCT(rprimd(:,3),rprimd(:,3)))/2.0_dp 124 ! eta=217.6_dp/zcut**2.0_dp 125 ! eta=1.0_dp/zcut**2.0_dp 126 ! eta=SQRT(16.0_dp/SQRT(DOT_PRODUCT(rprimd(:,1),rprimd(:,1)))) 127 ! else if (icutcoul.eq.2) then 128 ! zcut=SQRT(DOT_PRODUCT(rprimd(:,3),rprimd(:,3)))/2.0_dp 129 ! eta=SQRT(8.0_dp/zcut) 130 ! eta=SQRT(16.0_dp/SQRT(DOT_PRODUCT(rprimd(:,1),rprimd(:,1)))) 131 ! else if (icutcoul.eq.2) then 132 ! zcut=SQRT(DOT_PRODUCT(rprimd(:,3),rprimd(:,3)))/2.0_dp 133 ! eta=SQRT(8.0_dp/zcut) 134 else 135 eta=pi*200.0_dp/33.0_dp*sqrt(1.69_dp*recip/direct) 136 end if 137 138 !Conduct reciprocal space summations 139 fac=pi**2/eta 140 gsum=0._dp 141 grewtn(:,:)=0.0_dp 142 143 !Initialize Gcut-off array from m_gtermcutoff 144 !ABI_MALLOC(gcutoff,(ngfft(1)*ngfft(2)*ngfft(3))) 145 call termcutoff(gcutoff,gsqcut,icutcoul,ngfft,nkpt,rcut,rprimd,vcutgeo) 146 147 !if (icutcoul.eq.3) then 148 !Sum over G space, done shell after shell until all 149 !contributions are too small. 150 ng=0 151 do 152 ng=ng+1 153 newg=0 154 ! Instead of this warning that most normal users do not understand (because they are doing GS calculations, and not RF calculations), 155 ! one should optimize this routine. But usually this is a very small fraction of any ABINIT run. 156 ! if (ng > 20 .and. mod(ng,10)==0) then 157 ! write (msg,'(3a,I10)') "Very large box of G neighbors in ewald: you probably do not want to do this.", ch10,& 158 !& " If you have a metal consider setting dipdip 0. ng = ", ng 159 ! ABI_WARNING(msg) 160 ! end if 161 ii=1 162 do ig3=-ng,ng 163 do ig2=-ng,ng 164 do ig1=-ng,ng 165 ! Exclude shells previously summed over 166 if(abs(ig1)==ng .or. abs(ig2)==ng .or. abs(ig3)==ng .or. ng==1 ) then 167 168 ! gsq is G dot G = |G|^2 169 gsq=gmet(1,1)*dble(ig1*ig1)+gmet(2,2)*dble(ig2*ig2)+& 170 & gmet(3,3)*dble(ig3*ig3)+2._dp*(gmet(2,1)*dble(ig1*ig2)+& 171 & gmet(3,1)*dble(ig1*ig3)+gmet(3,2)*dble(ig3*ig2)) 172 173 ! Skip g=0: 174 if (gsq>1.0d-20) then 175 arg=fac*gsq 176 177 ! Larger arg gives 0 contribution because of exp(-arg) 178 if (arg <= -minexparg ) then 179 ! When any term contributes then include next shell 180 newg=1 181 182 if((abs(ig1).lt.ngfft(1)).and.& 183 &(abs(ig2).lt.ngfft(2)).and.& 184 &(abs(ig3).lt.ngfft(3))) then 185 ig23=ngfft(1)*(abs(ig2)+ngfft(2)*(abs(ig3))) 186 ii=abs(ig1)+ig23+1 187 !term= ( exp(-arg) + gcutoff(ii) - 1.0_dp )/gsq 188 !term=exp(-arg)/gsq*gcutoff(ii) 189 !term= ( exp(-arg) + gcutoff(ii) - 1.0_dp)/gsq 190 term=exp(-arg)/gsq*gcutoff(ii) 191 else if (icutcoul.ne.3) then 192 term=zero !exp(-arg)/gsq 193 else 194 term=exp(-arg)/gsq 195 endif 196 197 summr = 0.0_dp 198 summi = 0.0_dp 199 200 201 ! XG 20180531 : the two do-loops on ia should be merged, in order to spare 202 ! the waste of computing twice the sin and cos. 203 204 ! Note that if reduced atomic coordinates xred drift outside 205 ! of unit cell (outside [0,1)) it is irrelevant in the following 206 ! term, which only computes a phase. 207 do ia=1,natom 208 arg=two_pi*(ig1*xred(1,ia)+ig2*xred(2,ia)+ig3*xred(3,ia)) 209 ! Sum real and imaginary parts (avoid complex variables) 210 summr=summr+zion(typat(ia))*cos(arg) 211 summi=summi+zion(typat(ia))*sin(arg) 212 end do 213 214 ! The following two checks avoid an annoying underflow error msg 215 if (abs(summr)<1.d-16) summr=0.0_dp 216 if (abs(summi)<1.d-16) summi=0.0_dp 217 218 ! The product of term and summr**2 or summi**2 below 219 ! can underflow if not for checks above 220 t1=term*(summr*summr+summi*summi) 221 gsum=gsum+t1 222 223 do ia=1,natom 224 ! Again only phase is computed so xred may fall outside [0,1). 225 arg=two_pi*(ig1*xred(1,ia)+ig2*xred(2,ia)+ig3*xred(3,ia)) 226 phr= cos(arg) 227 phi=-sin(arg) 228 ! (note: do not need real part, commented out) 229 ! c1r=(phr*summr-phi*summi)*(term*zion(typat(ia))) 230 c1i=(phi*summr+phr*summi)*(term*zion(typat(ia))) 231 ! compute coordinate gradients 232 grewtn(1,ia)=grewtn(1,ia)-c1i*ig1 233 grewtn(2,ia)=grewtn(2,ia)-c1i*ig2 234 grewtn(3,ia)=grewtn(3,ia)-c1i*ig3 235 end do 236 237 end if ! End condition of not larger than -minexparg 238 end if ! End skip g=0 239 end if ! End triple loop over G s and associated new shell condition 240 241 end do 242 end do 243 end do 244 245 ! Check if new shell must be calculated 246 if (newg==0) exit 247 248 end do ! End the loop on ng (new shells). Note that there is one exit from this loop. 249 !endif 250 251 sumg=gsum/(two_pi*ucvol) 252 253 !Stress tensor is now computed elsewhere (ewald2) hence do not need 254 !length scale gradients (used to compute them here). 255 256 !normalize coordinate gradients by unit cell volume ucvol 257 term=-2._dp/ucvol 258 grewtn(:,:)=grewtn(:,:)*term 259 !call DSCAL(3*natom,term,grewtn,1) 260 261 !Conduct real space summations 262 reta=sqrt(eta) 263 fac=2._dp*sqrt(eta/pi) 264 sumr=0.0_dp 265 266 !In the following a summation is being conducted over all 267 !unit cells (ir1, ir2, ir3) so it is appropriate to map all 268 !reduced coordinates xred back into [0,1). 269 ! 270 !Loop on shells in r-space as was done in g-space 271 nr=0 272 do 273 nr=nr+1 274 newr=0 275 ! Instead of this warning that most normal users do not understand (because they are doing GS calculations, and not RF calculations), 276 ! one should optimize this routine. But usually this is a very small fraction of any ABINIT run. 277 ! if (nr > 20 .and. mod(nr,10)==0) then 278 ! write (msg,'(3a,I10)') "Very large box of R neighbors in ewald: you probably do not want to do this.", ch10,& 279 !& " If you have a metal consider setting dipdip 0. nr = ", nr 280 ! ABI_WARNING(msg) 281 ! end if 282 ! 283 do ir3=-nr,nr 284 do ir2=-nr,nr 285 do ir1=-nr,nr 286 if( abs(ir3)==nr .or. abs(ir2)==nr .or. abs(ir1)==nr .or. nr==1 )then 287 288 do ia=1,natom 289 ! Map reduced coordinate xred(mu,ia) into [0,1) 290 fraca1=xred(1,ia)-aint(xred(1,ia))+0.5_dp-sign(0.5_dp,xred(1,ia)) 291 fraca2=xred(2,ia)-aint(xred(2,ia))+0.5_dp-sign(0.5_dp,xred(2,ia)) 292 fraca3=xred(3,ia)-aint(xred(3,ia))+0.5_dp-sign(0.5_dp,xred(3,ia)) 293 drdta1=0.0_dp 294 drdta2=0.0_dp 295 drdta3=0.0_dp 296 297 do ib=1,natom 298 ! fraca and fracb should be precomputedi and become arrays with natom dimension. 299 ! Also the combination with dble(ir1), dble(ir2), dble(ir3) or fraca should be done outside of the ib loop. 300 fracb1=xred(1,ib)-aint(xred(1,ib))+0.5_dp-sign(0.5_dp,xred(1,ib)) 301 fracb2=xred(2,ib)-aint(xred(2,ib))+0.5_dp-sign(0.5_dp,xred(2,ib)) 302 fracb3=xred(3,ib)-aint(xred(3,ib))+0.5_dp-sign(0.5_dp,xred(3,ib)) 303 r1=dble(ir1)+fracb1-fraca1 304 r2=dble(ir2)+fracb2-fraca2 305 r3=dble(ir3)+fracb3-fraca3 306 rsq=rmet(1,1)*r1*r1+rmet(2,2)*r2*r2+rmet(3,3)*r3*r3+& 307 & 2.0_dp*(rmet(2,1)*r2*r1+rmet(3,2)*r3*r2+rmet(3,1)*r1*r3) 308 309 ! Avoid zero denominators in 'term': 310 if (rsq>=1.0d-24) then 311 312 ! Note: erfc(8) is about 1.1e-29, so do not bother with larger arg. 313 ! Also: exp(-64) is about 1.6e-28, so do not bother with larger arg**2 in exp. 314 term=0._dp 315 if (eta*rsq<64.0_dp) then 316 newr=1 317 rmagn=sqrt(rsq) 318 arg=reta*rmagn 319 ! derfc is the real(dp) complementary error function 320 derfc_arg = abi_derfc(arg) 321 term=derfc_arg/rmagn 322 sumr=sumr+zion(typat(ia))*zion(typat(ib))*term 323 term=zion(typat(ia))*zion(typat(ib))*& 324 & (term+fac*exp(-eta*rsq))/rsq 325 ! Length scale grads now handled with stress tensor in ewald2 326 r1a1d=rmet(1,1)*r1+rmet(1,2)*r2+rmet(1,3)*r3 327 r2a2d=rmet(2,1)*r1+rmet(2,2)*r2+rmet(2,3)*r3 328 r3a3d=rmet(3,1)*r1+rmet(3,2)*r2+rmet(3,3)*r3 329 ! Compute terms related to coordinate gradients 330 drdta1=drdta1+term*r1a1d 331 drdta2=drdta2+term*r2a2d 332 drdta3=drdta3+term*r3a3d 333 end if 334 end if ! End avoid zero denominators in'term' 335 end do ! end loop over ib: 336 337 grewtn(1,ia)=grewtn(1,ia)+drdta1 338 grewtn(2,ia)=grewtn(2,ia)+drdta2 339 grewtn(3,ia)=grewtn(3,ia)+drdta3 340 end do ! end loop over ia: 341 end if 342 end do ! end triple loop over real space points and associated condition of new shell 343 end do 344 end do 345 346 ! Check if new shell must be calculated 347 if(newr==0) exit 348 end do ! End loop on nr (new shells). Note that there is an exit within the loop 349 ! 350 sumr=0.5_dp*sumr 351 fac=pi*ch**2.0_dp/(2.0_dp*eta*ucvol) 352 353 !Finally assemble Ewald energy, eew 354 if(icutcoul.ne.3) then 355 !eew=sumg+sumr-chsq*reta/sqrt(pi)-fac 356 eew=sumg+sumr-chsq*reta/sqrt(pi) 357 else 358 eew=sumg+sumr-chsq*reta/sqrt(pi)-fac 359 end if 360 361 ABI_FREE(gcutoff) 362 363 !DEBUG 364 !write(std_out,*)'eew=sumg+sumr-chsq*reta/sqrt(pi)-fac' 365 !write(std_out,*)eew,sumg,sumr,chsq*reta/sqrt(pi),fac 366 !ENDDEBUG 367 368 !Length scale grads handled with stress tensor, ewald2 369 370 !Output the final values of ng and nr 371 ! write(msg, '(a,a,i4,a,i4)' )ch10,' ewald : nr and ng are ',nr,' and ',ng 372 ! call wrtout(std_out,msg,'COLL') 373 374 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.
SOURCE
409 subroutine ewald2(gmet,natom,ntypat,rmet,rprimd,stress,typat,ucvol,xred,zion) 410 411 !Arguments ------------------------------------ 412 !scalars 413 integer,intent(in) :: natom,ntypat 414 real(dp),intent(in) :: ucvol 415 !arrays 416 integer,intent(in) :: typat(natom) 417 real(dp),intent(in) :: gmet(3,3),rmet(3,3),rprimd(3,3),xred(3,natom) 418 real(dp),intent(in) :: zion(ntypat) 419 real(dp),intent(out) :: stress(6) 420 421 !Local variables------------------------------- 422 !scalars 423 integer :: ia,ib,ig1,ig2,ig3,ir1,ir2,ir3,newg,newr,ng,nr 424 real(dp) :: arg1,arg2,arg3,ch,dderfc,derfc_arg,direct,eta,fac,fraca1 425 real(dp) :: fraca2,fraca3,fracb1,fracb2,fracb3,g1,g2,g3,gsq,r1,r1c,r2,r2c 426 real(dp) :: minexparg 427 real(dp) :: r3,r3c,recip,reta,rmagn,rsq,summi,summr,t1,t2,t3,t4,t5,t6,term1 428 real(dp) :: term2,term3,term4 429 !arrays 430 real(dp) :: gprimd(3,3),strg(6),strr(6) 431 432 ! ************************************************************************* 433 434 !Define dimensional reciprocal space primitive translations gprimd 435 !(inverse transpose of rprimd) 436 call matr3inv(rprimd,gprimd) 437 438 !This is the minimum argument of an exponential, with some safety 439 minexparg=log(tiny(0._dp))+five 440 441 !Add up total charge and sum of charge^2 in cell 442 ch=0._dp 443 do ia=1,natom 444 ch=ch+zion(typat(ia)) 445 end do 446 447 !Compute eta, the Ewald summation convergence parameter, 448 !for approximately optimized summations: 449 direct=rmet(1,1)+rmet(1,2)+rmet(1,3)+rmet(2,1)+& 450 & rmet(2,2)+rmet(2,3)+rmet(3,1)+rmet(3,2)+rmet(3,3) 451 recip=gmet(1,1)+gmet(1,2)+gmet(1,3)+gmet(2,1)+& 452 & gmet(2,2)+gmet(2,3)+gmet(3,1)+gmet(3,2)+gmet(3,3) 453 !Here, a bias is introduced, because G-space summation scales 454 !better than r space summation ! 455 eta=pi*200.0_dp/33.0_dp*sqrt(1.69_dp*recip/direct) 456 457 fac=pi**2/eta 458 459 !Conduct reciprocal space summations 460 strg(1:6)=0.0_dp 461 462 !Sum over G space, done shell after shell until all 463 !contributions are too small 464 ng=0 465 do 466 ng=ng+1 467 newg=0 468 469 do ig3=-ng,ng 470 do ig2=-ng,ng 471 do ig1=-ng,ng 472 473 ! Exclude shells previously summed over 474 if(abs(ig1)==ng .or. abs(ig2)==ng .or. abs(ig3)==ng .or. ng==1 ) then 475 476 ! Compute Cartesian components of each G 477 g1=gprimd(1,1)*ig1+gprimd(1,2)*ig2+gprimd(1,3)*ig3 478 g2=gprimd(2,1)*ig1+gprimd(2,2)*ig2+gprimd(2,3)*ig3 479 g3=gprimd(3,1)*ig1+gprimd(3,2)*ig2+gprimd(3,3)*ig3 480 ! Compute |G|^2 (no pi factors) 481 gsq=(g1**2+g2**2+g3**2) 482 483 ! skip g=0: 484 if (gsq>1.0d-20) then 485 arg1=fac*gsq 486 487 ! larger arg1 gives 0 contribution because of exp(-arg1) 488 if (arg1<= -minexparg) then 489 ! When any term contributes then include next shell 490 newg=1 491 term1=exp(-arg1)/arg1 492 summr = 0.0_dp 493 summi = 0.0_dp 494 do ia=1,natom 495 arg2=two_pi*(ig1*xred(1,ia)+ig2*xred(2,ia)+ig3*xred(3,ia)) 496 ! Sum real and imaginary parts (avoid complex variables) 497 summr=summr+zion(typat(ia))*cos(arg2) 498 summi=summi+zion(typat(ia))*sin(arg2) 499 end do 500 501 ! Avoid underflow error messages 502 if (abs(summr)<1.d-16) summr=0.0_dp 503 if (abs(summi)<1.d-16) summi=0.0_dp 504 505 term2=(2._dp/gsq)*(1._dp+arg1) 506 t1=term2*g1*g1-1._dp 507 t2=term2*g2*g2-1._dp 508 t3=term2*g3*g3-1._dp 509 t4=term2*g2*g3 510 t5=term2*g1*g3 511 t6=term2*g1*g2 512 term3=term1*(summr*summr+summi*summi) 513 strg(1)=strg(1)+t1*term3 514 strg(2)=strg(2)+t2*term3 515 strg(3)=strg(3)+t3*term3 516 strg(4)=strg(4)+t4*term3 517 strg(5)=strg(5)+t5*term3 518 strg(6)=strg(6)+t6*term3 519 520 end if ! End condition not being larger than -minexparg 521 end if ! End skip g=0 522 523 end if ! End triple loop and condition of new shell 524 end do 525 end do 526 end do 527 528 ! Check if new shell must be calculated 529 if (newg==0) exit 530 end do ! End loop on new shell. Note that there is an "exit" instruction within the loop 531 532 533 !Conduct real space summations 534 reta=sqrt(eta) 535 strr(1:6)=0.0_dp 536 537 !Loop on shells in r-space as was done in g-space 538 nr=0 539 do 540 nr=nr+1 541 newr=0 542 543 do ir3=-nr,nr 544 do ir2=-nr,nr 545 do ir1=-nr,nr 546 if( abs(ir3)==nr .or. abs(ir2)==nr .or. abs(ir1)==nr .or. nr==1 )then 547 548 do ia=1,natom 549 ! Convert reduced atomic coordinates to [0,1) 550 fraca1=xred(1,ia)-aint(xred(1,ia))+0.5_dp-sign(0.5_dp,xred(1,ia)) 551 fraca2=xred(2,ia)-aint(xred(2,ia))+0.5_dp-sign(0.5_dp,xred(2,ia)) 552 fraca3=xred(3,ia)-aint(xred(3,ia))+0.5_dp-sign(0.5_dp,xred(3,ia)) 553 do ib=1,natom 554 fracb1=xred(1,ib)-aint(xred(1,ib))+0.5_dp-sign(0.5_dp,xred(1,ib)) 555 fracb2=xred(2,ib)-aint(xred(2,ib))+0.5_dp-sign(0.5_dp,xred(2,ib)) 556 fracb3=xred(3,ib)-aint(xred(3,ib))+0.5_dp-sign(0.5_dp,xred(3,ib)) 557 r1=ir1+fracb1-fraca1 558 r2=ir2+fracb2-fraca2 559 r3=ir3+fracb3-fraca3 560 ! Convert from reduced to cartesian coordinates 561 r1c=rprimd(1,1)*r1+rprimd(1,2)*r2+rprimd(1,3)*r3 562 r2c=rprimd(2,1)*r1+rprimd(2,2)*r2+rprimd(2,3)*r3 563 r3c=rprimd(3,1)*r1+rprimd(3,2)*r2+rprimd(3,3)*r3 564 ! Compute |r|^2 565 rsq=r1c**2+r2c**2+r3c**2 566 rmagn=sqrt(rsq) 567 568 ! Avoid zero denominators in 'term': 569 if (rmagn>=1.0d-12) then 570 571 ! Note: erfc(8) is about 1.1e-29, so do not bother with larger arg. 572 ! Also: exp(-64) is about 1.6e-28, so do not bother with larger arg**2 in exp. 573 arg3=reta*rmagn 574 if (arg3<8.0_dp) then 575 newr=1 576 ! derfc computes the complementary error function 577 ! dderfc is the derivative of the complementary error function 578 dderfc=(-2/sqrt(pi))*exp(-eta*rsq) 579 derfc_arg = abi_derfc(arg3) 580 term3=dderfc-derfc_arg/arg3 581 term4=zion(typat(ia))*zion(typat(ib))*term3 582 strr(1)=strr(1)+term4*r1c*r1c/rsq 583 strr(2)=strr(2)+term4*r2c*r2c/rsq 584 strr(3)=strr(3)+term4*r3c*r3c/rsq 585 strr(4)=strr(4)+term4*r2c*r3c/rsq 586 strr(5)=strr(5)+term4*r1c*r3c/rsq 587 strr(6)=strr(6)+term4*r1c*r2c/rsq 588 end if ! End the condition of not being to large 589 end if ! End avoid zero denominator 590 591 end do ! End loop over ib: 592 end do ! End loop over ia: 593 594 end if ! End triple loop overs real space points, and associated new shell condition 595 end do 596 end do 597 end do 598 599 ! Check if new shell must be calculated 600 if(newr==0) exit 601 end do ! End loop on new shells 602 603 !Finally assemble stress tensor coming from Ewald energy, stress 604 !(note division by unit cell volume in accordance with definition 605 !found in Nielsen and Martin, Phys. Rev. B 32, 3792 (1985) [[cite:Nielsen1985a]] 606 607 fac = pi/(2._dp*ucvol*eta) 608 stress(1)=(0.5_dp*reta*strr(1)+fac*(strg(1)+(ch**2)))/ucvol 609 stress(2)=(0.5_dp*reta*strr(2)+fac*(strg(2)+(ch**2)))/ucvol 610 stress(3)=(0.5_dp*reta*strr(3)+fac*(strg(3)+(ch**2)))/ucvol 611 stress(4)=(0.5_dp*reta*strr(4)+fac*strg(4))/ucvol 612 stress(5)=(0.5_dp*reta*strr(5)+fac*strg(5))/ucvol 613 stress(6)=(0.5_dp*reta*strr(6)+fac*strg(6))/ucvol 614 615 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 (72) to (75). This has been generalized to quadrupoles. Delivers the left hand side of Eq.(72), possibly generalized.
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 qdrp_cart(3,3,3,natom)=Quadrupole tensor on each atom in cartesian cordinates option= 0: use old implementation; 1: reduce the smalest argument of the exponentials to be evaluated, set eta to 1 and skip real space sum, leads to a significant speedup [dipquad] = if 1, atmfrc has been build without dipole-quadrupole part [quadquad] = if 1, atmfrc has been build without quadrupole-quadrupole part
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 ....
SOURCE
671 subroutine ewald9(acell,dielt,dyew,gmet,gprim,natom,qphon,rmet,rprim,sumg0,ucvol,xred,zeff, qdrp_cart, & 672 option, dipquad, quadquad) ! optional 673 674 !Arguments ------------------------------- 675 !scalars 676 integer,intent(in) :: natom,sumg0 677 integer,optional,intent(in) :: option, dipquad, quadquad 678 real(dp),intent(in) :: ucvol 679 !arrays 680 real(dp),intent(in) :: acell(3),dielt(3,3),gmet(3,3),gprim(3,3),qphon(3) 681 real(dp),intent(in) :: rmet(3,3),rprim(3,3),xred(3,natom),zeff(3,3,natom) 682 real(dp),intent(in) :: qdrp_cart(3,3,3,natom) 683 real(dp),intent(out) :: dyew(2,3,natom,3,natom) 684 685 !Local variables ------------------------- 686 !scalars 687 integer,parameter :: mr=10000 688 integer :: ia,ib,ig1,ig2,ig3,ii,ll,kk,ir,ir1,ir2,ir3,jj 689 integer :: info,lwork,mu,newg,newr,ng,nr,nu,ng_expxq 690 integer :: ewald_option 691 integer :: dipquad_,quadquad_ 692 logical :: do_quadrupole 693 logical, save :: firstcall = .TRUE. 694 real(dp),parameter :: fac=4.0_dp/3.0_dp/sqrt(pi) 695 real(dp),parameter :: fact2=2.0_dp/sqrt(pi) 696 real(dp),parameter :: y2max=64.0_dp, y2min=1.0d-24 697 real(dp) :: cddi,cddr,cqdi,cqdr,cqqi,cqqr,g3,g4 698 real(dp) :: arg1,arg2,arg3,arga,c123r,c123i,c23i,c23r,detdlt,inv_detdlt 699 real(dp) :: direct,eta,fact1,fact3,gsq,recip,reta,reta3,inv4eta 700 real(dp) :: minexparg,sigma_max 701 real(dp) :: term1,term2,term3,term4,term5,y2,yy,invy,invy2,derfc_yy 702 character(len=700) :: msg 703 !arrays 704 real(dp) :: c1i(2*mr+1),c1r(2*mr+1),c2i(2*mr+1),c2r(2*mr+1),c3i(2*mr+1) 705 real(dp) :: c3r(2*mr+1),cosqxred(natom),wdielt(3,3),eig_dielt(3),gpq(3),gpqfac(3,3),gpqgpq(3,3) 706 real(dp) :: invdlt(3,3),ircar(3),ircax(3),rr(3),sinqxred(natom) 707 real(dp) :: xredcar(3,natom),xredcax(3,natom),xredicar(3),xredicax(3),xx(3) 708 real(dp) :: gprimbyacell(3,3) !,tsec(2) 709 real(dp),allocatable :: dyddt(:,:,:,:,:), dydqt(:,:,:,:,:,:), dyqqt(:,:,:,:,:,:,:) 710 real(dp),allocatable :: work(:) 711 complex(dpc) :: exp2piqx(natom) 712 complex(dpc),allocatable :: expx1(:,:), expx2(:,:), expx3(:,:) 713 714 ! ********************************************************************* 715 716 ! This routine is expensive so skip the calculation and return zeros if zeff == zero. 717 ! Typically this happens when the DDB file does not contains zeff but dipdip = 1 is used (default). 718 if (all(zeff == zero).and.all(qdrp_cart == zero)) then 719 dyew = zero; return 720 end if 721 do_quadrupole = any(qdrp_cart /= zero) 722 723 ! Keep track of total time spent. 724 !call timab(1749, 1, tsec) 725 726 ! Initialize dipquad and quadquad options 727 dipquad_=0; if(present(dipquad)) dipquad_=dipquad 728 quadquad_=0; if(present(quadquad)) quadquad_=quadquad 729 730 ! Deactivate real space sums for quadrupolar fields or for dipdip = -1 731 ewald_option = 0; if (present(option)) ewald_option = option 732 if (do_quadrupole.and.(dipquad_==1.or.quadquad_==1)) ewald_option = 1 733 !ewald_option = 0 734 735 !This is the minimum argument of an exponential, with some safety 736 minexparg=log(tiny(0._dp))+five 737 if (ewald_option == 1) minexparg=-20.0_dp 738 739 ! initialize complex phase factors 740 do ia = 1, natom 741 arga = two_pi*( (qphon(1))*xred(1,ia)& 742 +(qphon(2))*xred(2,ia)& 743 +(qphon(3))*xred(3,ia) ) 744 exp2piqx(ia) = exp(arga*j_dpc) 745 end do 746 ng_expxq = 1000 747 ABI_MALLOC(expx1, (-ng_expxq:ng_expxq, natom)) 748 ABI_MALLOC(expx2, (-ng_expxq:ng_expxq, natom)) 749 ABI_MALLOC(expx3, (-ng_expxq:ng_expxq, natom)) 750 do ia = 1, natom 751 do ig1 = -ng_expxq, ng_expxq 752 expx1(ig1, ia) = exp(ig1*two_pi*xred(1,ia)*j_dpc) 753 expx2(ig1, ia) = exp(ig1*two_pi*xred(2,ia)*j_dpc) 754 expx3(ig1, ia) = exp(ig1*two_pi*xred(3,ia)*j_dpc) 755 end do 756 end do 757 758 gprimbyacell = gprim 759 gprimbyacell(:,1) = gprimbyacell(:,1) / acell(1) 760 gprimbyacell(:,2) = gprimbyacell(:,2) / acell(2) 761 gprimbyacell(:,3) = gprimbyacell(:,3) / acell(3) 762 763 !compute eta for approximately optimized summations: 764 direct=rmet(1,1)+rmet(1,2)+rmet(1,3)+rmet(2,1)+& 765 & rmet(2,2)+rmet(2,3)+rmet(3,1)+rmet(3,2)+rmet(3,3) 766 recip=gmet(1,1)+gmet(1,2)+gmet(1,3)+gmet(2,1)+& 767 & gmet(2,2)+gmet(2,3)+gmet(3,1)+gmet(3,2)+gmet(3,3) 768 eta=pi*100.0_dp/33.0_dp*sqrt(1.69_dp*recip/direct) 769 770 ! Compute a material-dependent width for the Gaussians that hopefully 771 ! will make the Ewald real-space summation unnecessary. 772 if (ewald_option == 1) then 773 774 wdielt(:,:)=dielt(:,:) 775 776 !Diagonalize dielectric matrix 777 lwork=-1 778 ABI_MALLOC(work,(10)) 779 call dsyev('N','U',3, wdielt, 3, eig_dielt, work, lwork,info) 780 lwork=nint(work(1)) 781 ABI_FREE(work) 782 783 ABI_MALLOC(work,(lwork)) 784 call dsyev('V','U',3, wdielt, 3, eig_dielt, work, lwork,info) 785 ABI_FREE(work) 786 787 !This is a tentative maximum value for the gaussian width in real space 788 sigma_max=three 789 790 !Set eta taking into account that the eps_inf is used as a metric in 791 !reciprocal space 792 eta=sqrt(maxval(eig_dielt))/sigma_max 793 794 if (firstcall) then 795 firstcall = .FALSE. 796 write(msg, '(4a,f9.4,9a)' ) ch10,& 797 ' Warning : due to the use of quadrupolar fields, the width of the reciprocal space gaussians', ch10, & 798 ' in ewald9 has been set to eta= ', eta, ' 1/bohr and the real-space sums have been neglected.', ch10, & 799 ' One should check whether this choice leads to correct results for the specific system under study', & 800 ' and q-point grid.',ch10, & 801 ' It is recommended to check that calculations with dipdip=1 and -1 (both with dipquad=0 and quadquad=0)', ch10, & 802 ' lead to identical results. Otherwise increase the resolution of the q-point grid and repeat this test.', ch10 803 call wrtout([ab_out,std_out], msg) 804 end if 805 806 !Internally eta is the square of the gaussians width 807 eta=eta*eta 808 end if 809 810 inv4eta = one / four / eta 811 812 ABI_MALLOC(dyddt,(2,3,natom,3,natom)) 813 ABI_MALLOC(dydqt,(2,3,natom,3,natom,3)) 814 ABI_MALLOC(dyqqt,(2,3,natom,3,natom,3,3)) 815 816 dyddt = zero 817 dydqt = zero 818 dyqqt = zero 819 820 !Sum terms over g space: 821 ng=0 822 do 823 ng=ng+1 824 825 ! if needed, update the complex phases for larger G vectors 826 if (ng > ng_expxq) then 827 !write(std_out,*)"have to realloc" 828 ABI_FREE(expx1) 829 ABI_FREE(expx2) 830 ABI_FREE(expx3) 831 832 ng_expxq = ng_expxq*2 833 ! TODO: half of this space is not needed, as it contains the complex conjugate of the other half. 834 ! present duplication avoids if statements inside the loop, however 835 ABI_MALLOC(expx1, (-ng_expxq:ng_expxq, natom)) 836 ABI_MALLOC(expx2, (-ng_expxq:ng_expxq, natom)) 837 ABI_MALLOC(expx3, (-ng_expxq:ng_expxq, natom)) 838 do ia = 1, natom 839 do ig1 = -ng_expxq, ng_expxq 840 expx1(ig1, ia) = exp(ig1*two_pi*xred(1,ia)*j_dpc) 841 expx2(ig1, ia) = exp(ig1*two_pi*xred(2,ia)*j_dpc) 842 expx3(ig1, ia) = exp(ig1*two_pi*xred(3,ia)*j_dpc) 843 end do 844 end do 845 end if 846 847 newg=0 848 do ig3=-ng,ng 849 do ig2=-ng,ng 850 do ig1=-ng,ng 851 if(abs(ig1)==ng .or. abs(ig2)==ng .or. abs(ig3)==ng .or. ng==1 )then 852 853 gpq(1)=(ig1+qphon(1))*gprimbyacell(1,1)+(ig2+qphon(2))*& 854 & gprimbyacell(1,2)+(ig3+qphon(3))*gprimbyacell(1,3) 855 gpq(2)=(ig1+qphon(1))*gprimbyacell(2,1)+(ig2+qphon(2))*& 856 & gprimbyacell(2,2)+(ig3+qphon(3))*gprimbyacell(2,3) 857 gpq(3)=(ig1+qphon(1))*gprimbyacell(3,1)+(ig2+qphon(2))*& 858 & gprimbyacell(3,2)+(ig3+qphon(3))*gprimbyacell(3,3) 859 gsq=zero 860 do jj=1,3 861 do ii=1,3 862 gpqgpq(ii,jj)=gpq(ii)*gpq(jj) 863 gsq=gsq+gpqgpq(ii,jj)*dielt(ii,jj) 864 end do 865 end do 866 867 ! Skip q=0: 868 if (gsq<1.0d-20) then 869 if (sumg0==1) then 870 write(msg,'(a,a,a,a,a)' )& 871 'The phonon wavelength should not be zero :',ch10,& 872 'there are non-analytical terms that cannot be treated.',ch10,& 873 'Action: subtract this wavelength from the input file.' 874 ABI_ERROR(msg) 875 end if 876 877 else 878 879 arg1=(two_pi**2)*gsq* inv4eta 880 881 ! Larger arg gives 0 contribution: 882 if (arg1<= -minexparg ) then 883 newg=1 884 885 ! Here calculate the term 886 term1=exp(-arg1)/gsq 887 do jj=1,3 888 do ii=1,3 889 gpqfac(ii,jj)=gpqgpq(ii,jj)*term1 890 end do 891 end do 892 893 ! 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 894 do ia=1,natom 895 cosqxred(ia)= real(exp2piqx(ia)*expx1(ig1, ia)*expx2(ig2, ia)*expx3(ig3, ia)) 896 sinqxred(ia)=aimag(exp2piqx(ia)*expx1(ig1, ia)*expx2(ig2, ia)*expx3(ig3, ia)) 897 end do 898 899 ! First, the diagonal terms 900 do nu=1,3 901 do ia=1,natom 902 do mu=nu,3 903 dyddt(1,mu,ia,nu,ia)=dyddt(1,mu,ia,nu,ia)+gpqfac(mu,nu) 904 end do 905 end do 906 end do 907 908 ! Then, the non-diagonal ones 909 do ib=2,natom 910 do ia=1,ib-1 911 ! phase factor dipole-dipole 912 cddr=cosqxred(ia)*cosqxred(ib)+sinqxred(ia)*sinqxred(ib) 913 cddi=sinqxred(ia)*cosqxred(ib)-cosqxred(ia)*sinqxred(ib) 914 915 ! Dipole-dipole contribution 916 do nu=1,3 917 do mu=nu,3 918 dyddt(1,mu,ia,nu,ib)=dyddt(1,mu,ia,nu,ib)+gpqfac(mu,nu)*cddr 919 dyddt(2,mu,ia,nu,ib)=dyddt(2,mu,ia,nu,ib)+gpqfac(mu,nu)*cddi 920 end do 921 end do 922 end do 923 end do 924 925 if (do_quadrupole) then 926 do ib=1,natom 927 do ia=1,natom 928 929 ! phase factor for dipole-quadrupole 930 cqdr=cosqxred(ia)*sinqxred(ib)-sinqxred(ia)*cosqxred(ib) 931 cqdi=cosqxred(ia)*cosqxred(ib)+sinqxred(ia)*sinqxred(ib) 932 933 ! phase factor quadrupole-quadrupole 934 cqqr=cosqxred(ia)*cosqxred(ib)+sinqxred(ia)*sinqxred(ib) 935 cqqi=sinqxred(ia)*cosqxred(ib)-cosqxred(ia)*sinqxred(ib) 936 937 ! Dipole-quadrupole contribution 938 do ii=1,3 939 do jj=1,3 940 do kk=1,3 941 g3=gpq(ii)*gpq(jj)*gpq(kk) 942 dydqt(1,ii,ia,jj,ib,kk)=dydqt(1,ii,ia,jj,ib,kk)+g3*term1*cqdr 943 dydqt(2,ii,ia,jj,ib,kk)=dydqt(2,ii,ia,jj,ib,kk)+g3*term1*cqdi 944 end do ! kk 945 end do ! jj 946 end do ! ii 947 948 ! Quadrupole-quadrupole contribution 949 do ii=1,3 950 do jj=1,3 951 do kk=1,3 952 do ll=1,3 953 g4 = gpq(ii)*gpq(jj)*gpq(kk)*gpq(ll) 954 dyqqt(1,ii,ia,jj,ib,kk,ll)=dyqqt(1,ii,ia,jj,ib,kk,ll)+g4*term1*cqqr 955 dyqqt(2,ii,ia,jj,ib,kk,ll)=dyqqt(2,ii,ia,jj,ib,kk,ll)+g4*term1*cqqi 956 end do 957 end do ! kk 958 end do ! jj 959 end do ! ii 960 end do ! ia 961 end do ! ib 962 end if 963 964 end if ! endif exp() argument is smaller than -minexparg 965 end if ! Endif g/=0 : 966 end if ! End triple summation over Gs: 967 end do 968 end do 969 end do 970 971 ! Check if new shell must be calculated 972 if(newg==0)exit 973 end do 974 975 !Multiplies by common factor 976 fact1=4.0_dp*pi/ucvol 977 do ib=1,natom 978 do ia=1,ib 979 do nu=1,3 980 do mu=nu,3 981 dyddt(1,mu,ia,nu,ib)=dyddt(1,mu,ia,nu,ib)*fact1 982 dyddt(2,mu,ia,nu,ib)=dyddt(2,mu,ia,nu,ib)*fact1 983 end do 984 end do 985 end do 986 end do 987 if (do_quadrupole) then 988 dydqt=dydqt*fact1/two * two_pi 989 dyqqt=dyqqt*fact1/four * two_pi ** 2 990 end if 991 992 reta=sqrt(eta) 993 reta3=-eta*reta 994 995 !Calculating the inverse (transpose) of the dielectric tensor 996 call matr3inv(dielt,invdlt) 997 !Calculating the determinant of the dielectric tensor 998 detdlt=dielt(1,1)*dielt(2,2)*dielt(3,3)+dielt(1,3)*dielt(2,1)*& 999 & dielt(3,2)+dielt(1,2)*dielt(2,3)*dielt(3,1)-dielt(1,3)*& 1000 & dielt(2,2)*dielt(3,1)-dielt(1,1)*dielt(2,3)*dielt(3,2)-& 1001 & dielt(1,2)*dielt(2,1)*dielt(3,3) 1002 1003 if(detdlt<tol6)then 1004 write(msg, '(a,es16.6,11a)' )& 1005 'The determinant of the dielectrix matrix, detdlt=',detdlt,' is smaller than 1.0d-6.',ch10,& 1006 'The use of the dipole-dipole model for interatomic force constants is not possible.',ch10,& 1007 'It is likely that you have not treated the electric field perturbations,',ch10,& 1008 'because you not are dealing with an insulator, so that',ch10,& 1009 'your dielectric matrix was simply set to zero in the Derivative DataBase.',ch10,& 1010 'Action: set the input variable dipdip to 0 .' 1011 ABI_ERROR(msg) 1012 end if 1013 1014 inv_detdlt = one / sqrt(detdlt) 1015 fact3=reta3 * inv_detdlt 1016 1017 if (ewald_option /= 1) then 1018 ! Preparing the loop on real space 1019 do ia=1,natom 1020 do ii=1,3 1021 xredcar(ii,ia)=(xred(1,ia)*acell(1)*rprim(ii,1)+& 1022 xred(2,ia)*acell(2)*rprim(ii,2)+& 1023 xred(3,ia)*acell(3)*rprim(ii,3) )*reta 1024 end do 1025 end do 1026 do ia=1,natom 1027 do ii=1,3 1028 xredcax(ii,ia)= invdlt(1,ii)*xredcar(ii,ia)+& 1029 invdlt(2,ii)*xredcar(ii,ia)+& 1030 invdlt(3,ii)*xredcar(ii,ia) 1031 end do 1032 end do 1033 1034 ! Prepare the evaluation of exp(iq*R) 1035 do ir=-mr,mr 1036 arg1=-two_pi*qphon(1)*ir 1037 arg2=-two_pi*qphon(2)*ir 1038 arg3=-two_pi*qphon(3)*ir 1039 c1r(ir+mr+1)=cos(arg1) 1040 c1i(ir+mr+1)=sin(arg1) 1041 c2r(ir+mr+1)=cos(arg2) 1042 c2i(ir+mr+1)=sin(arg2) 1043 c3r(ir+mr+1)=cos(arg3) 1044 c3i(ir+mr+1)=sin(arg3) 1045 end do 1046 1047 do nr=1,mr 1048 newr=0 1049 1050 ! Begin big loop on real space vectors 1051 do ir3=-nr,nr 1052 do ir2=-nr,nr 1053 1054 ! Here, construct the cosine and sine of q*R for components 2 and 3 1055 c23r = c2r(ir2+mr+1) * c3r(ir3+mr+1) - c2i(ir2+mr+1) * c3i(ir3+mr+1) 1056 c23i = c2i(ir2+mr+1) * c3r(ir3+mr+1) + c2r(ir2+mr+1) * c3i(ir3+mr+1) 1057 1058 ! Also multiplies by fact3, because it is a rather economical place to do so 1059 c23r=c23r * fact3 1060 c23i=c23i * fact3 1061 1062 do ir1=-nr,nr 1063 if( abs(ir3)==nr .or. abs(ir2)==nr .or. abs(ir1)==nr .or. nr==1 )then 1064 1065 ! This is the real part and imaginary part of the phase factor exp(iq*R) 1066 c123r = c1r(ir1+mr+1) * c23r - c1i(ir1+mr+1) * c23i 1067 c123i = c1i(ir1+mr+1) * c23r + c1r(ir1+mr+1) * c23i 1068 1069 do ii=1,3 1070 ircar(ii)= ( ir1*acell(1)*rprim(ii,1)+& 1071 ir2*acell(2)*rprim(ii,2)+& 1072 ir3*acell(3)*rprim(ii,3) ) * reta 1073 end do 1074 do ii=1,3 1075 ircax(ii)= invdlt(1,ii)*ircar(ii)+& 1076 invdlt(2,ii)*ircar(ii)+& 1077 invdlt(3,ii)*ircar(ii) 1078 end do 1079 1080 ! Here loops on atoms 1081 do ib=1,natom 1082 do ii=1,3 1083 xredicar(ii)=ircar(ii)-xredcar(ii,ib) 1084 xredicax(ii)=ircax(ii)-xredcax(ii,ib) 1085 end do 1086 do ia=1,ib 1087 do ii=1,3 1088 rr(ii)=xredicar(ii)+xredcar(ii,ia) 1089 xx(ii)=xredicax(ii)+xredcax(ii,ia) 1090 end do 1091 1092 y2=rr(1)*xx(1)+rr(2)*xx(2)+rr(3)*xx(3) 1093 1094 ! The atoms should not be too far of each other 1095 if (y2 < y2max) then 1096 ! Note: erfc(8) is about 1.1e-29, so dont bother with larger y. 1097 ! Also: exp(-64) is about 1.6e-28, do dont bother with larger y**2 in exp. 1098 1099 ! Avoid zero denominators in term: 1100 if (y2 >= y2min) then 1101 newr=1 1102 yy=sqrt(y2) 1103 invy=1.0_dp/yy 1104 invy2=invy**2 1105 derfc_yy = abi_derfc(yy) 1106 term2=derfc_yy*invy*invy2 1107 term3=fact2*exp(-y2)*invy2 1108 term4=-(term2+term3) 1109 term5=(3.0_dp*term2+term3*(3.0_dp+2.0_dp*y2))*invy2 1110 do nu=1,3 1111 do mu=nu,3 1112 dyddt(1,mu,ia,nu,ib)=dyddt(1,mu,ia,nu,ib)+c123r*(xx(nu)*xx(mu)*term5+term4*invdlt(nu,mu)) 1113 dyddt(2,mu,ia,nu,ib)=dyddt(2,mu,ia,nu,ib)+c123i*(xx(nu)*xx(mu)*term5+term4*invdlt(nu,mu)) 1114 end do 1115 end do 1116 else 1117 ! If zero denominator, the atoms should be identical 1118 if (ia/=ib)then 1119 write(msg, '(5a,i0,a,i0,a)' )& 1120 'The distance between two atoms seem to vanish.',ch10,& 1121 'This is not allowed.',ch10,& 1122 'Action: check the input for the atoms number',ia,' and',ib,'.' 1123 ABI_ERROR(msg) 1124 else 1125 ! This is the correction when the atoms are identical 1126 do nu=1,3 1127 do mu=1,3 1128 dyddt(1,mu,ia,nu,ib)=dyddt(1,mu,ia,nu,ib)+& 1129 fac*reta3*invdlt(nu,mu) * inv_detdlt 1130 end do 1131 end do 1132 end if 1133 end if ! End the condition for avoiding zero denominators 1134 end if ! End the condition of too large distance between atoms 1135 end do 1136 end do ! End loop over ia and ib : 1137 end if ! End triple loop over real space points: 1138 end do ! ir1 1139 end do ! ir2 1140 end do ! ir3 1141 1142 ! Check if new shell must be calculated 1143 if(newr==0)exit 1144 if(newr==1 .and. nr==mr) ABI_BUG('mr is too small') 1145 end do 1146 end if ! check if should compute real part 1147 1148 !Now, symmetrizes 1149 do ib=1,natom-1 1150 do nu=1,3 1151 do ia=ib+1,natom 1152 do mu=nu,3 1153 dyddt(1,mu,ia,nu,ib)= dyddt(1,mu,ib,nu,ia) 1154 dyddt(2,mu,ia,nu,ib)=-dyddt(2,mu,ib,nu,ia) 1155 end do 1156 end do 1157 end do 1158 end do 1159 1160 do ib=1,natom 1161 do nu=2,3 1162 do ia=1,natom 1163 do mu=1,nu-1 1164 dyddt(1,mu,ia,nu,ib)=dyddt(1,nu,ia,mu,ib) 1165 dyddt(2,mu,ia,nu,ib)=dyddt(2,nu,ia,mu,ib) 1166 end do 1167 end do 1168 end do 1169 end do 1170 1171 !Tests 1172 !write(std_out,*)' ewald9 : take into account the effective charges ' 1173 dyew = zero 1174 do ib=1,natom 1175 do nu=1,3 1176 do ia=1,natom 1177 do mu=1,3 1178 do ii=1,3 1179 do jj=1,3 1180 ! dipole-dipole correction 1181 dyew(1,mu,ia,nu,ib)=dyew(1,mu,ia,nu,ib) + & 1182 zeff(ii,mu,ia)*zeff(jj,nu,ib)*dyddt(1,ii,ia,jj,ib) 1183 dyew(2,mu,ia,nu,ib)=dyew(2,mu,ia,nu,ib) + & 1184 zeff(ii,mu,ia)*zeff(jj,nu,ib)*dyddt(2,ii,ia,jj,ib) 1185 if (do_quadrupole) then 1186 do kk=1,3 1187 if (dipquad_==1) then 1188 ! dipole-quadrupole correction 1189 dyew(1,mu,ia,nu,ib)=dyew(1,mu,ia,nu,ib) + & 1190 (zeff(ii,nu,ib)*qdrp_cart(kk,jj,mu,ia) - & 1191 zeff(ii,mu,ia)*qdrp_cart(kk,jj,nu,ib)) * dydqt(1,ii,ia,jj,ib,kk) 1192 dyew(2,mu,ia,nu,ib)=dyew(2,mu,ia,nu,ib) + & 1193 (zeff(ii,nu,ib)*qdrp_cart(kk,jj,mu,ia) - & 1194 zeff(ii,mu,ia)*qdrp_cart(kk,jj,nu,ib)) * dydqt(2,ii,ia,jj,ib,kk) 1195 end if 1196 1197 ! quadrupole-quadrupole correction 1198 if (quadquad_==1) then 1199 do ll=1,3 1200 dyew(1,mu,ia,nu,ib)=dyew(1,mu,ia,nu,ib) + & 1201 (qdrp_cart(ll,ii,mu,ia)*qdrp_cart(kk,jj,nu,ib)) * dyqqt(1,ii,ia,jj,ib,kk,ll) 1202 dyew(2,mu,ia,nu,ib)=dyew(2,mu,ia,nu,ib) + & 1203 (qdrp_cart(ll,ii,mu,ia)*qdrp_cart(kk,jj,nu,ib)) * dyqqt(2,ii,ia,jj,ib,kk,ll) 1204 end do 1205 end if 1206 end do 1207 end if 1208 1209 end do 1210 end do 1211 end do 1212 end do 1213 end do 1214 end do 1215 1216 ABI_FREE(expx1) 1217 ABI_FREE(expx2) 1218 ABI_FREE(expx3) 1219 ABI_FREE(dyddt) 1220 ABI_FREE(dydqt) 1221 ABI_FREE(dyqqt) 1222 1223 !call timab(1749, 2, tsec) 1224 1225 end subroutine ewald9