1c 2c $Id$ 3c 4 subroutine xc_vdw_der(s6,s8,sr6,sr8,a1,a2,n,x,z,force) 5c 6c S. Grimme J Comp Chem 25, 1463 (2004) 7c U. Zimmerli, M Parrinello and P. Koumoutsakos, JCP. 120, 2693 (2004) 8c Q. Wu and W. Yang, JCP. 116, 515 (2002) 9c 10 implicit none 11#include "mafdecls.fh" 12#include "errquit.fh" 13#include "xc_vdw.fh" 14#include "global.fh" 15 double precision s6,s8,sr6,sr8,a1,a2 16 integer n 17 double precision x(3,n),force(3,n) 18 integer z(n) 19c 20 integer i,j,k,A,l_cnij,k_cnij,l_cnijk,k_cnijk 21 double precision c6ij_sk 22 external c6ij_sk 23 double precision drajdxa 24 double precision ff1,rr,ff 25 double precision fdmp,f1dmp,cnA,cnj 26 external c6cn,crd_nr 27 double precision c6cn,crd_nr 28 double precision fac6,fac8,fdmp6,fdmp6a,fdmp8,fdmp8a,Qfac 29 double precision rAj,rAk,rjk,r0aj,r0ak,r0jk,c6Aj,grad_c6(3) 30 double precision dxAj,dyAj,dzAj,dxAk,dyAk,dzAk,dxjk,dyjk,dzjk 31 double precision tmp6,tmp6a,tmp8,tmp8a 32 double precision xc_fdmpbj, xc_fdmpbj_d1 33 external xc_fdmpbj, xc_fdmpbj_d1 34c 35c Derivatives of Grimme dispersion term 36c 37c DFT-D1 / DFT-D2 38c 39 if (ivdw.le.2) then 40 do A=1,n 41 force(1,A)=0d0 42 force(2,A)=0d0 43 force(3,A)=0d0 44 if (Z(A).ne.0) then 45 do j=1,n 46 if(A.ne.j) then 47 rAj=sqrt( 48 + (x(1,A)-x(1,j))**2 + 49 + (x(2,A)-x(2,j))**2 + 50 + (x(3,A)-x(3,j))**2) 51c protect from NaNs caused by bqs 52 if(raj.lt.1d30.or.abs(raj).lt.1d-8) then 53 r0aj=r0(z(A))+r0(z(j)) 54 ff= fdmp(rAj,r0aj) 55 ff1= f1dmp(rAj,r0aj,ff) 56 rr=c6ij_sk(A,j,z)/(rAj**6)* 57 * ((-6d0*ff/rAj)+ff1) 58 do i=1,3 59 drAjdxa=(x(i,A)-x(i,j))/rAj 60 force(i,A)=force(i,A)-rr*drAjdxa 61 enddo 62 endif 63 endif 64 enddo 65 endif 66 enddo 67 if(abs(s6-1d0).gt.1d-9) 68 F call dscal(3*n,s6,force,1) 69c 70c DFT-D3 71c 72 else if (ivdw.eq.3) then 73c 74c Precompute coordinate derivatives C6 dependency 75c 76 if (.not.ma_push_get(mt_dbl,3*n,'xcvdw cnij',l_cnij,k_cnij)) 77 & call errquit('xcvdw cnij: cannot allocate cnij',0, MA_ERR) 78 if (.not.ma_push_get(mt_dbl,3*n*n,'vdw cnijk',l_cnijk,k_cnijk)) 79 & call errquit('vdw cnijk: cannot allocate cnijk',0, MA_ERR) 80c 81 call crd_nr_der(n,x,z,dbl_mb(k_cnij),dbl_mb(k_cnijk)) 82c 83 do A=1,n 84 force(1,A)=0.0d0 85 force(2,A)=0.0d0 86 force(3,A)=0.0d0 87 if (Z(A).ne.0) then 88 do j=1,n 89 if(A.ne.j) then 90 dxAj=x(1,A)-x(1,j) 91 dyAj=x(2,A)-x(2,j) 92 dzAj=x(3,A)-x(3,j) 93 rAj=dxAj**2+dyAj**2+dzAj**2 94c 95c Two center derivatives. Grimme uses screening to reduce 96c computational work 97c 98c Screening r^2 distance vs threshold of 20000.0 99c 100 if (rAj.gt.20000.d0) goto 901 101c 102c Factors 103c 104 r0aj=r0AB(z(A),z(j)) 105 Qfac=Qatom(z(A))*Qatom(z(j)) 106 fac6=(dsqrt(rAj)/(sr6*r0aj))**(-alpha) 107 fac8=(dsqrt(rAj)/(sr8*r0aj))**(-(alpha+2.0d0)) 108 fdmp6=1.0d0/(1.0d0+6.0d0*fac6) 109 fdmp8=1.0d0/(1.0d0+6.0d0*fac8) 110c 111c Coordination dependent C6_AB value 112c 113 cnA=crd_nr(A,n,x,z) 114 cnj=crd_nr(j,n,x,z) 115 c6Aj=c6cn(z(A),z(j),cnA,cnj) 116c 117c Get gradient for coordination number dependent C6 118c 119 call c6_grad(grad_c6,A,j,A,x,z,n, 120 & dbl_mb(k_cnij),dbl_mb(k_cnijk)) 121c 122 tmp6=6.0d0*fdmp6*s6*c6Aj/(rAj**4.0d0) 123 tmp8=6.0d0*fdmp8*s8*c6Aj*Qfac/(rAj**5.0d0) 124c 125c dx contribution to A 126c 127 tmp6a=tmp6*dxAj 128 tmp8a=tmp8*dxAj 129 force(1,A)=force(1,A) 130 $ +(1.0d0-fdmp6*fac6*alpha)*tmp6a 131 $ -fdmp6*s6*grad_c6(1)/(rAj**3.0d0) 132 $ +(4.0d0-3.0d0*fdmp8*fac8*(alpha+2.0d0))*tmp8a 133 $ -3.0d0*fdmp8*s8*grad_c6(1)*Qfac/(rAj**4.0d0) 134c 135c dy contribution to A 136c 137 tmp6a=tmp6*dyAj 138 tmp8a=tmp8*dyAj 139 force(2,A)=force(2,A) 140 $ +(1.0d0-fdmp6*fac6*alpha)*tmp6a 141 $ -fdmp6*s6*grad_c6(2)/(rAj**3.0d0) 142 $ +(4.0d0-3.0d0*fdmp8*fac8*(alpha+2.0d0))*tmp8a 143 $ -3.0d0*fdmp8*s8*grad_c6(2)*Qfac/(rAj**4.0d0) 144c 145c dz contribution to A 146c 147 tmp6a=tmp6*dzAj 148 tmp8a=tmp8*dzAj 149 force(3,A)=force(3,A) 150 $ +(1.0d0-fdmp6*fac6*alpha)*tmp6a 151 $ -fdmp6*s6*grad_c6(3)/(rAj**3.0d0) 152 $ +(4.0d0-3.0d0*fdmp8*fac8*(alpha+2.0d0))*tmp8a 153 $ -3.0d0*fdmp8*s8*grad_c6(3)*Qfac/(rAj**4.0d0) 154 901 continue 155 endif 156 enddo 157c 158c Three center derivatives. Grimme uses aggressive screening 159c to get this N^3 contribution back to N^2 160c 161 do j=2,n 162 if(A.ne.j) then 163 rAj=sqrt( 164 + (x(1,A)-x(1,j))**2 + 165 + (x(2,A)-x(2,j))**2 + 166 + (x(3,A)-x(3,j))**2) 167 r0aj=r0AB(z(A),z(j)) 168c 169c Screening per Grimme 170c 171 if (rAj.gt.1600d0*r0aj/r0AB(1,1)) goto 910 172c 173c Third center involved 174c 175 do k=1,j-1 176 if(A.ne.k) then 177 dxAk=x(1,A)-x(1,k) 178 dyAk=x(2,A)-x(2,k) 179 dzAk=x(3,A)-x(3,k) 180 rAk=dxAk**2+dyAk**2+dzAk**2 181 r0ak=r0AB(z(A),z(k)) 182 dxjk=x(1,j)-x(1,k) 183 dyjk=x(2,j)-x(2,k) 184 dzjk=x(3,j)-x(3,k) 185 rjk=dxjk**2+dyjk**2+dzjk**2 186 r0jk=r0AB(z(j),z(k)) 187c 188c Screening r^2 distance vs threshold of 1600.0*(radii Ak) 189c 190 if ((rAk.gt.1600.0d0*r0ak/r0AB(1,1)).or. 191 $ (rjk.gt.1600.0d0*r0jk/r0AB(1,1))) goto 911 192c 193c Get gradient for coordination number dependent C6 for three centers 194c 195 call c6_grad(grad_c6,j,k,A,x,z,n, 196 & dbl_mb(k_cnij),dbl_mb(k_cnijk)) 197 fac6=(sr6*r0jk/dsqrt(rjk))**(alpha) 198 fac8=(sr8*r0jk/dsqrt(rjk))**(alpha+2.0d0) 199 fdmp6=1.0d0/(1.0d0+6.0d0*fac6) 200 fdmp8=1.0d0/(1.0d0+6.0d0*fac8) 201c 202c dx, dy, and dz contribution to A 203c 204 Qfac=Qatom(z(j))*Qatom(z(k)) 205 force(1,A)=force(1,A) 206 $ -fdmp6*s6*grad_c6(1)/(rjk**3.0d0) 207 $ -3.0d0*fdmp8*s8*grad_c6(1)*Qfac/(rjk**4.0d0) 208 force(2,A)=force(2,A) 209 $ -fdmp6*s6*grad_c6(2)/(rjk**3.0d0) 210 $ -3.0d0*fdmp8*s8*grad_c6(2)*Qfac/(rjk**4.0d0) 211 force(3,A)=force(3,A) 212 $ -fdmp6*s6*grad_c6(3)/(rjk**3.0d0) 213 $ -3.0d0*fdmp8*s8*grad_c6(3)*Qfac/(rjk**4.0d0) 214 endif 215 911 continue 216 enddo 217 910 continue 218 endif 219 enddo 220 endif 221 enddo 222 if (.not.ma_pop_stack(l_cnijk)) 223 $ call errquit('xcvdw cnijk: cannot pop cnijk',4, MA_ERR) 224 if (.not.ma_pop_stack(l_cnij)) 225 $ call errquit('xcvdw cnij: cannot pop cnij',4, MA_ERR) 226c 227c DFT-D3BJ 228c 229 else if (ivdw.eq.4) then 230c 231c Precompute coordinate derivatives C6 dependency 232c 233 if (.not.ma_push_get(mt_dbl,3*n,'xcvdw cnij',l_cnij,k_cnij)) 234 & call errquit('xcvdw cnij: cannot allocate cnij',0, MA_ERR) 235 if (.not.ma_push_get(mt_dbl,3*n*n,'vdw cnijk',l_cnijk,k_cnijk)) 236 & call errquit('vdw cnijk: cannot allocate cnijk',0, MA_ERR) 237c 238 call crd_nr_der(n,x,z,dbl_mb(k_cnij),dbl_mb(k_cnijk)) 239c 240 do A=1,n 241 force(1,A)=0.0d0 242 force(2,A)=0.0d0 243 force(3,A)=0.0d0 244 if (Z(A).ne.0) then 245 do j=1,n 246 if(A.ne.j.and.Z(j).ne.0) then 247 dxAj=x(1,A)-x(1,j) 248 dyAj=x(2,A)-x(2,j) 249 dzAj=x(3,A)-x(3,j) 250 rAj=dxAj**2+dyAj**2+dzAj**2 251c 252c Two center derivatives. Grimme uses screening to reduce 253c computational work 254c 255c Screening r^2 distance vs threshold of 20000.0 256c 257 if (rAj.gt.20000.d0) goto 941 258c 259c Factors 260c 261 rAj = dsqrt(rAj) 262 Qfac=Qatom(z(A))*Qatom(z(j)) 263c 264c Coordination dependent C6_AB value 265c 266 cnA=crd_nr(A,n,x,z) 267 cnj=crd_nr(j,n,x,z) 268 c6Aj=c6cn(z(A),z(j),cnA,cnj) 269 c8=3.0d0*c6Aj*Qfac 270 r0Aj=dsqrt(3.0d0*Qfac) 271c 272c Get gradient for coordination number dependent C6 273c 274 call c6_grad(grad_c6,A,j,A,x,z,n, 275 & dbl_mb(k_cnij),dbl_mb(k_cnijk)) 276c 277c dx contribution to A 278c 279 force(1,A)=force(1,A) 280 $ -s6*c6Aj*xc_fdmpbj_d1(rAj,r0Aj,a1,a2,6) 281 $ *dxAj/rAj 282 $ -s8*c8*xc_fdmpbj_d1(rAj,r0Aj,a1,a2,8) 283 $ *dxAj/rAj 284 $ -s6*grad_c6(1)*xc_fdmpbj(rAj,r0Aj,a1,a2,6) 285 $ -s8*3.0d0*grad_c6(1)*Qfac*xc_fdmpbj(rAj,r0Aj, 286 $ a1,a2,8) 287c 288c dy contribution to A 289c 290 force(2,A)=force(2,A) 291 $ -s6*c6Aj*xc_fdmpbj_d1(rAj,r0Aj,a1,a2,6) 292 $ *dyAj/rAj 293 $ -s8*c8*xc_fdmpbj_d1(rAj,r0Aj,a1,a2,8) 294 $ *dyAj/rAj 295 $ -s6*grad_c6(2)*xc_fdmpbj(rAj,r0Aj,a1,a2,6) 296 $ -s8*3.0d0*grad_c6(2)*Qfac*xc_fdmpbj(rAj,r0Aj, 297 $ a1,a2,8) 298c 299c dz contribution to A 300c 301 force(3,A)=force(3,A) 302 $ -s6*c6Aj*xc_fdmpbj_d1(rAj,r0Aj,a1,a2,6) 303 $ *dzAj/rAj 304 $ -s8*c8*xc_fdmpbj_d1(rAj,r0Aj,a1,a2,8) 305 $ *dzAj/rAj 306 $ -s6*grad_c6(3)*xc_fdmpbj(rAj,r0Aj,a1,a2,6) 307 $ -s8*3.0d0*grad_c6(3)*Qfac*xc_fdmpbj(rAj,r0Aj, 308 $ a1,a2,8) 309 941 continue 310 endif 311 enddo 312c 313c Three center derivatives. Grimme uses aggressive screening 314c to get this N^3 contribution back to N^2 315c 316 do j=1,n 317 if(A.ne.j.and.z(j).ne.0) then 318 rAj=sqrt( 319 + (x(1,A)-x(1,j))**2 + 320 + (x(2,A)-x(2,j))**2 + 321 + (x(3,A)-x(3,j))**2) 322 r0aj=r0AB(z(A),z(j)) 323c 324c Screening per Grimme 325c 326 if (rAj.gt.1600d0*r0aj/r0AB(1,1)) goto 950 327c 328c Third center involved 329c 330 do k=1,n 331 if(A.ne.k.and.k.ne.j.and.z(k).ne.0) then 332 dxAk=x(1,A)-x(1,k) 333 dyAk=x(2,A)-x(2,k) 334 dzAk=x(3,A)-x(3,k) 335 rAk=dxAk**2+dyAk**2+dzAk**2 336 r0ak=r0AB(z(A),z(k)) 337 dxjk=x(1,j)-x(1,k) 338 dyjk=x(2,j)-x(2,k) 339 dzjk=x(3,j)-x(3,k) 340 rjk=dxjk**2+dyjk**2+dzjk**2 341 r0jk=r0AB(z(j),z(k)) 342c 343c Screening r^2 distance vs threshold of 1600.0*(radii Ak) 344c 345 if ((rAk.gt.1600.0d0*r0ak/r0AB(1,1)).or. 346 $ (rjk.gt.1600.0d0*r0jk/r0AB(1,1))) goto 951 347c 348c Get gradient for coordination number dependent C6 for three centers 349c 350 call c6_grad(grad_c6,j,k,A,x,z,n, 351 & dbl_mb(k_cnij),dbl_mb(k_cnijk)) 352 rjk=dsqrt(rjk) 353c 354c dx, dy, and dz contribution to A 355c 356 Qfac=Qatom(z(j))*Qatom(z(k)) 357 fdmp6=xc_fdmpbj(rjk,dsqrt(3.0d0*Qfac),a1,a2,6) 358 fdmp8=xc_fdmpbj(rjk,dsqrt(3.0d0*Qfac),a1,a2,8) 359 force(1,A)=force(1,A) 360 $ -fdmp6*s6*grad_c6(1) 361 $ -3.0d0*fdmp8*s8*grad_c6(1)*Qfac 362 force(2,A)=force(2,A) 363 $ -fdmp6*s6*grad_c6(2) 364 $ -3.0d0*fdmp8*s8*grad_c6(2)*Qfac 365 force(3,A)=force(3,A) 366 $ -fdmp6*s6*grad_c6(3) 367 $ -3.0d0*fdmp8*s8*grad_c6(3)*Qfac 368 endif 369 951 continue 370 enddo 371 950 continue 372 endif 373 enddo 374 endif 375 enddo 376 if (.not.ma_pop_stack(l_cnijk)) 377 $ call errquit('xcvdw cnijk: cannot pop cnijk',4, MA_ERR) 378 if (.not.ma_pop_stack(l_cnij)) 379 $ call errquit('xcvdw cnij: cannot pop cnij',4, MA_ERR) 380 endif 381c 382#ifdef DEBUG 383 write(6,*) ' gradient vdw called' 384#endif 385 return 386 end 387C> 388C> \brief Evaluate the dispersion energy 389C> 390C> This function evaluates the dispersion energy based on an empirical 391C> expression. This routine supports multiple expressions commonly 392C> used in density functional theory. 393C> 394C> The DFT-D3 correction with BJ damping [4,5] is given by 395C> \f{eqnarray*}{ 396C> E_{\mathrm{disp}}^{D3(BJ)} &=& -\frac{1}{2}\sum_{A\ne B} 397C> s_6\frac{C_6^{AB}}{R_{AB}^6+[f(R^0_{AB})]^6} + 398C> s_8\frac{C_8^{AB}}{R_{AB}^8+[f(R^0_{AB})]^8} \\\\ 399C> f(R^0_{AB}) &=& a_1R^0_{AB}+a_2 \\\\ 400C> R^0_{AB} &=& \sqrt{\frac{C_8^{AB}}{C_6^{AB}}} 401C> \f} 402C> 403C> \return The dispersion energy \f$E_{\mathrm{disp}}\f$. 404C> 405C> ### References ### 406C> 407C> [1] S. Grimme, 408C> "Accurate description of van der Waals complexes by density 409C> functional theory including empirical corrections", 410C> J. Comp. Chem. (2004) <b>25</b>, 1463-1473, DOI: 411C> <a href="https://doi.org/10.1002/jcc.20078"> 412C> 10.1002/jcc.20078</a>. 413C> 414C> [2] U. Zimmerli, M. Parrinello, P. Koumoutsakos, 415C> "Dispersion corrections to density functionals for water 416C> aromatic interactions", 417C> J. Chem. Phys. (2004) <b>120</b>, 2693, DOI: 418C> <a href="https://doi.org/10.1063/1.1637034"> 419C> 10.1063/1.1637034</a>. 420C> 421C> [3] Q. Wu, W. Yang, 422C> "Empirical correction to density functional theory for van der 423C> Waals interactions", 424C> J. Chem. Phys. (2002) <b>116</b>, 515, DOI: 425C> <a href="https://doi.org/10.1063/1.1424928"> 426C> 10.1063/1.1424928</a>. 427C> 428C> [4] A.D. Becke, E.R. Johnson, 429C> "A unified density-functional treatment of dynamical, 430C> nondynamical and dispersion correlations", 431C> J. Chem. Phys. (2007) <b>127</b> 124108, DOI: 432C> <a href="https://doi.org/10.1063/1.2768530"> 433C> 10.1063/1.2768530</a> (See appendix C). 434C> 435C> [5] S. Grimme, S. Ehrlich, L. Goerigk, 436C> "Effect of the damping function in dispersion corrected 437C> density functional theory", J. Comput. Chem. (2011) 438C> <b>32</b>, pp. 1456-1465, DOI: 439C> <a href="https://doi.org/10.1002/jcc.21759"> 440C> 10.1002/jcc.21759</a> (See Eqs.(5-6)). 441C> 442 double precision function xc_vdw_e(s6,s8,sr6,sr8,a1,a2,n,x,z) 443c 444c S. Grimme J Comp Chem 25, 1463 (2004) 445c U. Zimmerli, M Parrinello and P. Koumoutsakos, JCP. 120, 2693 (2004) 446c Q. Wu and W. Yang, JCP. 116, 515 (2002) 447c 448 implicit none 449#include "errquit.fh" 450#include "xc_vdw.fh" 451 double precision s6 !< [Input] The \f$s_6\f$ coefficient 452 double precision s8 !< [Input] The \f$s_8\f$ coefficient 453 double precision sr6 !< [Input] The \f$s_{r,6}\f$ coefficient 454 double precision sr8 !< [Input] The \f$s_{r,8}\f$ coefficient 455 double precision a1 !< [Input] The \f$a_1\f$ coefficient 456 double precision a2 !< [Input] The \f$a_2\f$ coefficient 457 integer n !< [Input] The number of atoms 458 double precision x(3,n) !< [Input] The atomic coordinates 459 integer z(n) !< [Input] The atomic numbers of the atoms 460c 461 integer i,j 462 integer i6, i8 463 parameter(i6 = 6) 464 parameter(i8 = 8) 465 double precision fdmp, fdmp3, cni, cnj,c6d3,xc_fdmpbj 466 double precision c6ij_sk,rij,c6cn,crd_nr,e6,e8,r0bj 467 external c6ij_sk,c6cn,nxtask,crd_nr,fdmp3,xc_fdmpbj 468c 469 xc_vdw_e=0d0 470 e6=0.0d0 471 e8=0.0d0 472c 473c DFT-D1 / DFT-D2 474c 475 if (ivdw.le.2) then 476 do i=1,n-1 477 if (Z(i).ne.0) then 478 if (r0(z(i)).le.0.0d0) then 479 call errquit("xc_vdw_e: no Grimme parameters for element", 480 + int(z(i)),UERR) 481 endif 482 do j=i+1,n 483 if (Z(j).ne.0) then 484 rij=dsqrt( 485 + (x(1,i)-x(1,j))**2 + 486 + (x(2,i)-x(2,j))**2 + 487 + (x(3,i)-x(3,j))**2) 488c protect from NaNs caused by bqs 489 if(rij.gt.1d30.or.abs(rij).lt.1d-4) then 490 xc_vdw_e=0d0 491 else 492 xc_vdw_e=xc_vdw_e-c6ij_sk(i,j,z)* 493 * fdmp(rij,r0(z(i))+r0(z(j)))* 494 * (rij)**(-6.0d0) 495 endif 496 endif 497 enddo 498 endif 499 enddo 500 xc_vdw_e=xc_vdw_e*s6 501c 502c DFT-D3 503c 504c As off August, 2011 Grimme states: "Adding three-body corrections is 505c currently not recommended, as very little is known about the three- 506c body behaviour of common DFs in overlapping density regions." 507c http://toc.uni-muenster.de/DFTD3/data/man.pdf, section 1.3. 508c Hence the three-body terms have not been implemented. 509c 510c The reference to three-center derivatives in the gradient code 511c refers to contributions that come from differentiating the 512c coordination dependent dispersion coefficients. 513c 514 else if (ivdw.eq.3) then 515 do i=1,n-1 516 if (Z(i).ne.0) then 517 do j=i+1,n 518 if (Z(j).ne.0) then 519 rij=dsqrt( 520 + (x(1,i)-x(1,j))**2 + 521 + (x(2,i)-x(2,j))**2 + 522 + (x(3,i)-x(3,j))**2) 523c protect from NaNs caused by bqs 524 if(rij.lt.1d30.or.abs(rij).lt.1d-8) then 525 cni=crd_nr(i,n,x,z) 526 cnj=crd_nr(j,n,x,z) 527 c6d3=c6cn(z(i),z(j),cni,cnj) 528 c8=3.0d0*c6d3*Qatom(z(i))*Qatom(z(j)) 529 e6=e6-c6d3*fdmp3(rij,r0AB(z(i),z(j))*sr6,alpha)* 530 * (rij)**(-6.0d0) 531 e8=e8-c8*fdmp3(rij,r0AB(z(i),z(j))*sr8,alpha+2.0d0)* 532 * (rij)**(-8.0d0) 533 endif 534 endif 535 enddo 536 endif 537 enddo 538 xc_vdw_e=e6*s6+e8*s8 539 else if (ivdw.eq.4) then 540 do i=1,n-1 541 if (Z(i).ne.0) then 542 do j=i+1,n 543 if (Z(j).ne.0) then 544 rij=dsqrt( 545 + (x(1,i)-x(1,j))**2 + 546 + (x(2,i)-x(2,j))**2 + 547 + (x(3,i)-x(3,j))**2) 548c protect from NaNs caused by bqs 549 if(rij.lt.1d30.or.abs(rij).lt.1d-8) then 550 cni=crd_nr(i,n,x,z) 551 cnj=crd_nr(j,n,x,z) 552 c6d3=c6cn(z(i),z(j),cni,cnj) 553 c8=3.0d0*c6d3*Qatom(z(i))*Qatom(z(j)) 554 r0bj=dsqrt(c8/c6d3) 555 e6=e6-c6d3*xc_fdmpbj(rij,r0bj,a1,a2,i6) 556 e8=e8-c8*xc_fdmpbj(rij,r0bj,a1,a2,i8) 557 endif 558 endif 559 enddo 560 endif 561 enddo 562 xc_vdw_e=e6*s6+e8*s8 563 endif 564 return 565 end 566c 567 subroutine xc_vdw(rtdb,geom,exc,force,what) 568 implicit none 569 character *(*) what 570 integer geom,rtdb 571 double precision exc,force(*),s6,s8,sr6,sr8,a1,a2 572#include "geom.fh" 573#include "mafdecls.fh" 574#include "errquit.fh" 575#include "util.fh" 576#include "stdio.fh" 577#include "global.fh" 578#include "rtdb.fh" 579#include "xc_vdw.fh" 580c 581 integer n 582 integer itags,ltags,i_xyz,l_xyz,icharge,lcharge, 583 I l_fvdw,i_fvdw, i_xyz2,l_xyz2,i_iz2,l_iz2 584 external xc_vdw_e 585 double precision xc_vdw_e,evdw,scalea 586 integer iz,lz,i 587 logical xc_vdw_init 588 external xc_vdw_init 589 logical oprint,oprinth 590 logical stat 591 logical use_nwxc_disp,out1 592c 593 double precision delta,delta_default 594 double precision s6_in 595c 596 oprint = util_print('vdw', print_medium) 597 oprinth = util_print('vdw high', print_high) 598c 599c Allocate memory blocks 600c 601 if (.not. geom_ncent(geom, n)) 602 & call errquit('xcvdw: geom_ncent failed',geom, GEOM_ERR) 603 if (.not.MA_push_get(MT_Dbl,n*3,'xyz',l_xyz,i_xyz)) 604 & call errquit('xcvdw: cannot allocate xyz',0, MA_ERR) 605 if (.not.MA_Push_Get(MT_int,n,'atns',lz,iz)) 606 & call errquit('xcvdw: cannot allocate atns',0, MA_ERR) 607 if (.not.MA_Push_Get(MT_Byte,n*16,'tags',ltags,itags)) 608 & call errquit('xcvdw: cannot allocate tags',0, MA_ERR) 609 if (.not.MA_Push_Get(MT_Dbl,n,'charge',lcharge,icharge)) 610 & call errquit('xcvdw: cannot allocate charge',0, MA_ERR) 611 if (.not. geom_cart_get2(geom, n, Byte_MB(itags), 612 & Dbl_MB(i_xyz), Dbl_MB(icharge), int_mb(iz))) 613 & call errquit('xcvdw: geom_cart_get failed',74, GEOM_ERR) 614 if (.not.ma_pop_stack(lcharge)) 615 & call errquit('xcvdw: cannot pop stack',3, MA_ERR) 616c 617c Which Grimme dispersion version 618c 619 if (.not.rtdb_get(rtdb,'dft:ivdw',MT_INT,1,ivdw)) 620 & ivdw = 2 621c 622c get rid of bqs 623c 624 if (.not.MA_push_get(mt_dbl,n*3,'xyz',l_xyz2,i_xyz2)) 625 & call errquit('xcvdw: cannot allocate xyz',0, MA_ERR) 626 if (.not.MA_push_get(mt_int,n,'iz',l_iz2,i_iz2)) 627 & call errquit('xcvdw: cannot allocate xyz',0, MA_ERR) 628 call xc_vdw_nobqs(geom, Byte_MB(itags), 629 N n, dbl_mb(i_xyz), int_mb(iz), 630 T dbl_mb(i_xyz2),int_mb(i_iz2)) 631 if (.not.ma_chop_stack(l_xyz2)) 632 C call errquit('xcvdw: cannot pop stack',14, MA_ERR) 633c 634 use_nwxc_disp = .false. 635 if(util_module_avail("nwxc")) then 636 call nwxc_getvals("nwxc_has_disp",out1) 637 use_nwxc_disp = out1 638 if (use_nwxc_disp) then 639c get ivdw here as this data is needed for xc_vdw_init 640 call nwxc_get_disp(ivdw,s6,s8,sr6,sr8,a1,a2,alpha) 641 endif 642 endif 643c 644c conversion factor angs 2 au 645c 646 if(.not.geom_get_ang2au(geom, scalea)) call 647 S errquit('xcvdw: gang2au failed',0,0) 648c 649c Initialize are variables 650c 651 if(.not.xc_vdw_init(scalea)) 652 & call errquit('xcvdw: vwdinit failed',0, 0) 653c 654c Read in some user defined parameters 655c 656 if (.not.rtdb_get(rtdb,'dft:vdwalpha',MT_DBL,1,alpha)) 657 & alpha = 20.0d0 658 if (ivdw.eq.3) alpha = 14.0d0 659 if (ivdw.eq.4) alpha = 14.0d0 660c 661c Get proper scaling factors depending on Grimme dispersion version 662c 663 if (.not.use_nwxc_disp) then 664 if(.not.rtdb_get(rtdb,'dft:vdw_s6',mt_dbl,1,s6)) s6=0d0 665 if(.not.rtdb_get(rtdb,'dft:vdw_s8',mt_dbl,1,s8)) s8=0d0 666 if(.not.rtdb_get(rtdb,'dft:vdw_sr6',mt_dbl,1,sr6)) sr6=0d0 667 if(.not.rtdb_get(rtdb,'dft:vdw_sr8',mt_dbl,1,sr8)) sr8=0d0 668 if(.not.rtdb_get(rtdb,'dft:vdw_a1',mt_dbl,1,a1)) a1=0d0 669 if(.not.rtdb_get(rtdb,'dft:vdw_a2',mt_dbl,1,a2)) a2=0d0 670 call get_scaling_fac(s6,s8,sr6,sr8,a1,a2) 671 if(rtdb_get(rtdb, 'dft:vdw', mt_dbl, 1, s6_in)) then 672 s6=s6_in 673 if(ga_nodeid().eq.0)write(6,*) ' WARNING: vdw s6 set = ',s6 674 endif 675 else 676c get the scaling factors here (just being paranoid) 677 call nwxc_get_disp(ivdw,s6,s8,sr6,sr8,a1,a2,alpha) 678 endif 679c 680 if(what.eq.'energy') then 681c 682c Compute energy contribution 683c 684 if(oprinth.and.ga_nodeid().eq.0) then 685 write(luout,*) ' s6 =',s6 686 write(luout,*) ' s8 =',s8 687 write(luout,*) ' sr6 =',sr6 688 write(luout,*) ' sr8 =',sr8 689 write(luout,*) ' alpha =',alpha 690 write(luout,*) ' ivdw =',ivdw 691 endif 692c 693 evdw=xc_vdw_e(s6,s8,sr6,sr8,a1,a2,n,dbl_mb(i_xyz),int_mb(iz)) 694c 695 if(oprint.and.ga_nodeid().eq.0) then 696 write(luout,*) 697 D ' Dispersion Parameters' 698 write(luout,*) 699 D ' ---------------------' 700 if (ivdw.eq.1.or.ivdw.eq.2) then 701 write(luout,222) s6, evdw 702 222 format(/ 703 & ' s6 scale factor :', f22.12/ 704 & ' vdW contrib :', f22.12/) 705 endif 706 if (ivdw.eq.3) then 707 write(luout,223) s6, s8, sr6, sr8, evdw 708 223 format(/ 709 & ' DFT-D3 Model ', / 710 & ' s8 scale factor :', f22.12/ 711 & ' sr6 scale factor :', f22.12/ 712 & ' sr8 scale factor :', f22.12/ 713 & ' vdW contrib :', f22.12/) 714 endif 715 if (ivdw.eq.4) then 716 write(luout,224) s6, s8, a1, a2, evdw 717 224 format(/ 718 & ' DFT-D3BJ Model ', / 719 & ' s6 scale factor :', f22.12/ 720 & ' s8 scale factor :', f22.12/ 721 & ' a1 parameter :', f22.12/ 722 & ' a2 parameter :', f22.12/ 723 & ' vdW contrib :', f22.12/) 724 endif 725 endif 726c 727c Add contribution to Exc 728c 729 Exc=Exc+evdw 730c 731 elseif(what.eq.'forces') then 732c 733c Gradient calculation 734c 735 if (.not.MA_push_get(MT_Dbl,n*3,'xyz',l_fvdw,i_fvdw)) 736 & call errquit('xcvdw: cannot allocate forcev',0, MA_ERR) 737c 738 call xc_vdw_der(s6,s8,sr6,sr8,a1,a2,n,dbl_mb(i_xyz),int_mb(iz), 739 D dbl_mb(i_fvdw)) 740c 741 if(oprinth.and.ga_nodeid().eq.0) then 742 write(luout,*) ' vdW contrib for S6=',s6 743 do i=1,n 744 write(luout,'(I2,3F10.7," F = ",3(1PE13.5))') 745 Z int_mb(iz+i-1), 746 X dbl_mb(i_xyz+3*(i-1)), 747 Y dbl_mb(i_xyz+3*(i-1)+1), 748 Z dbl_mb(i_xyz+3*(i-1)+2), 749 X dbl_mb(i_fvdw+3*(i-1)), 750 Y dbl_mb(i_fvdw+3*(i-1)+1), 751 Z dbl_mb(i_fvdw+3*(i-1)+2) 752 enddo 753 write(luout,*) ' before vdw contr @@@@@' 754 do i=1,n 755 write(luout,'(I2,3F10.7," F = ",3(1PE13.5))') 756 Z int_mb(iz+i-1), 757 X dbl_mb(i_xyz+3*(i-1)), 758 Y dbl_mb(i_xyz+3*(i-1)+1), 759 Z dbl_mb(i_xyz+3*(i-1)+2), 760 X force(1+3*(i-1)), 761 Y force(1+3*(i-1)+1), 762 Z force(1+3*(i-1)+2) 763 enddo 764 765 endif 766c 767c Add to force matrix 768c 769 call daxpy(3*n,1d0,dbl_mb(i_fvdw),1,force,1) 770c 771 if(oprinth.and.ga_nodeid().eq.0) then 772 write(luout,*) ' after vdw contr @@@@@' 773 do i=1,n 774 write(luout,'(I2,3F10.7," F = ",3(1PE13.5))') 775 Z int_mb(iz+i-1), 776 X dbl_mb(i_xyz+3*(i-1)), 777 Y dbl_mb(i_xyz+3*(i-1)+1), 778 Z dbl_mb(i_xyz+3*(i-1)+2), 779 X force(1+3*(i-1)), 780 Y force(1+3*(i-1)+1), 781 Z force(1+3*(i-1)+2) 782 enddo 783 endif 784c 785 elseif(what.eq.'hessian') then 786c 787c Hessian calculation, numerical from analytical gradients 788c 789c Get delta as used in a numerical hessian DFT calculation 790c 791 delta_default = 0.01d0 792 if (.not.rtdb_get(rtdb,'stpr_gen:delta',MT_DBL,1,delta)) 793 & delta = delta_default 794c 795 call xc_vdw_hess(delta,s6,s8,sr6,sr8,a1,a2,n, 796 & dbl_mb(i_xyz),int_mb(iz)) 797c 798 if(oprint.and.ga_nodeid().eq.0) then 799 write(luout,*) ' s6 = ',s6 800 write(luout,*) ' vdw to hessian contribution is done' 801 endif 802 endif 803c 804c Clean up 805c 806 if (.not.ma_chop_stack(l_xyz)) 807 C call errquit('xcvdw: cannot pop stack',4, MA_ERR) 808c 809 return 810 end 811c 812 subroutine xc_vdw_hess(delta,s6,s8,sr6,sr8,a1,a2,n,x,z) 813 implicit none 814c This function makes vdw empirical correction to the hessian 815c must be called before thermochemical data and vibrational 816c analysis are generated. 817#include "inp.fh" 818#include "util.fh" 819#include "mafdecls.fh" 820#include "stdio.fh" 821#include "errquit.fh" 822#include "global.fh" 823 double precision s6,s8,sr6,sr8,a1,a2 824 integer n 825 double precision x(3,n) 826 integer z(n) 827 double precision l_force(3,n) 828 double precision r_force(3,n) 829 double precision hessvdw(3,n,3,n) 830 double precision dval0 831 double precision dval1 832 double precision delta 833 integer i, j, A, B 834 integer n3xyz 835 integer nat2 836 integer nhesst 837 integer lenhess 838 integer l_exybs,k_exybs 839 integer l_exybt,k_exybt 840c 841 character*(nw_max_path_len) filehess 842c 843c dispersion contribution to hessian 844c 845c in principle this task is very fast, 846c so only master node works, and read-write to disk. 847c 848 if (ga_nodeid().eq.0) then 849c 850 call util_file_name('hess', .false., .false.,filehess) 851 852 lenhess=inp_strlen(filehess) 853 n3xyz=3*n 854 nhesst=n3xyz*(n3xyz+1)/2 855 nat2=n3xyz*n3xyz 856c 857 if (.not.ma_push_get(mt_dbl,nat2,'xcvdwhess exybs ', 858 * l_exybs,k_exybs)) 859 & call errquit('xcvdwhess: cannot allocate exybs',0, MA_ERR) 860 call dfill(nat2,0.0d00,dbl_mb(k_exybs),1) 861c 862 if (.not.ma_push_get(mt_dbl,nhesst,'xcvdwhess exybt ', 863 * l_exybt,k_exybt)) 864 & call errquit('xcvdwhess: cannot allocate exybt',0, MA_ERR) 865 call dfill(nhesst,0.0d00,dbl_mb(k_exybt),1) 866c 867 write(LuOut,* ) 'Read old hessian file : ', filehess 868c 869 call dfill(nhesst,0.d0,dbl_mb(k_exybt),1) 870 !write(6,* ) 'leee ' 871 open(unit=69,file=filehess,form='formatted',status='old', 872 & err=99900,access='sequential') 873 do j = 0,(nhesst-1) 874 read(69,*,err=99901,end=99902) dval0 875 dbl_mb(k_exybt+j) = dval0 876 !write(6,* ) 'dval ', dval0 877 enddo 878 close(unit=69,status='keep') 879c 880c complete the square matrix from triangle values 881c 882 call trin2squa(dbl_mb(k_exybs),dbl_mb(k_exybt),n3xyz) 883c 884 write(LuOut,* ) 'vdW contribution to hessian ' 885c 886 call output(dbl_mb(k_exybs), 1, n3xyz, 1, n3xyz, n3xyz, n3xyz, 1) 887c 888 do A = 1, n 889 do i = 1, 3 890 do B=1, n 891 do j=1, 3 892 r_force(j,B)=0.0d0 893 l_force(j,B)=0.0d0 894 enddo 895 enddo 896c right displacement 897 x(i,A) = x(i,A) + delta 898 call xc_vdw_der(s6,s8,sr6,sr8,a1,a2,n,x,z,r_force) 899c left displacement 900 x(i,A) = x(i,A) - 2.0d0*delta 901 call xc_vdw_der(s6,s8,sr6,sr8,a1,a2,n,x,z,l_force) 902c back to original position 903 x(i,A) = x(i,A) + delta 904c 905 do B=1,n 906 do j=1, 3 907 dval1 = (r_force(j,B)-l_force(j,B)) / (2.0d00*delta) 908 hessvdw(i,A,j,B)=dval1 909 enddo 910 enddo 911 enddo 912 enddo 913c 914 call output(hessvdw,1,n3xyz,1,n3xyz, 915 & n3xyz,n3xyz,1) 916c 917 call daxpy(nat2,1d0,hessvdw,1,dbl_mb(k_exybs),1) 918c 919c:write the final hessian 920c 921 call output(dbl_mb(k_exybs),1,n3xyz,1,n3xyz, 922 & n3xyz,n3xyz,1) 923 call stpr_wrt_fd_from_sq(dbl_mb(k_exybs),n3xyz,filehess) 924 write(LuOut,* ) 'New hessian file vdw corrected has been 925 . written:', filehess 926 if (.not.ma_chop_stack(l_exybs)) 927 C call errquit('xcvdwhess: cannot pop stack exybs',4, MA_ERR) 928 929 endif 930c 931 call ga_sync 932c 933 return 93499900 continue 935 write(luout,*)'hess_file => ',filehess 936 call errquit('xc_vdw_hess 99900', 911, DISK_ERR) 93799901 continue 938 write(luout,*)'hess_file => ',filehess 939 call errquit('xc_vdw_hess 99901', 911, DISK_ERR) 94099902 continue 941 write(luout,*)'hess_file => ',filehess 942 call errquit('xc_vdw_hess 99902', 911, DISK_ERR) 943 end 944c 945 subroutine xc_vdw_to_hessian(rtdb) 946c 947#include "mafdecls.fh" 948#include "errquit.fh" 949#include "rtdb.fh" 950#include "geom.fh" 951 integer rtdb 952 integer geom 953 double precision dum 954 double precision dum1(1) 955 character*255 name 956c 957 if (.not. rtdb_cget(rtdb,'geometry', 1, name)) 958 $ name = 'geometry' 959c 960 if (.not. geom_create(geom, name)) 961 $ call errquit('xc_vdw_to_hessian: geom_create failed !', 962 $ 0,GEOM_ERR) 963c 964 if (.not. geom_rtdb_load(rtdb, geom, name)) 965 $ call errquit('xc_vdw_to_hessian: no geometry load form rtdb', 0, 966 $ GEOM_ERR) 967c 968 call xc_vdw(rtdb, geom, dum, dum1, 'hessian') 969c 970 if(.not. geom_destroy(geom)) 971 $ call errquit('xc_vdw_to_hessian: geom_create failed !', 972 $ 0,GEOM_ERR) 973c 974 return 975 end 976 subroutine xc_vdw_nobqs(geom, tags, 977 N n, xyz, z, 978 T tempxyz, tempz) 979 implicit none 980#include "geom.fh" 981#include "inp.fh" 982 double precision xyz(3,*),tempxyz(3,*) 983 character*16 tags(*) 984 integer n 985 integer z(*),tempz(*) 986 integer geom 987c 988 integer i,n_left 989c 990 do i=1,n 991 tempz(i)=z(i) 992 tempxyz(1,i)=xyz(1,i) 993 tempxyz(2,i)=xyz(2,i) 994 tempxyz(3,i)=xyz(3,i) 995 enddo 996 n_left=0 997 do i=1,n 998 n_left=n_left+1 999 z(n_left)=tempz(i) 1000 if (inp_compare(.false.,tags(i)(1:2),'bq')) then 1001c move bq at a crazy large distance to kill its contribution 1002 xyz(1,n_left)=1d30 1003 xyz(2,n_left)=1d30 1004 xyz(3,n_left)=1d30 1005 else 1006 xyz(1,n_left)=tempxyz(1,i) 1007 xyz(2,n_left)=tempxyz(2,i) 1008 xyz(3,n_left)=tempxyz(3,i) 1009 endif 1010 enddo 1011 n=n_left 1012 return 1013 end 1014