diff --git a/CHANGELOG.rst b/CHANGELOG.rst index d3758ff..9a44eb5 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -3,6 +3,23 @@ Change Log ========== +0.9 +=== + +Added +----- + +* Implementation of revised damping parameters as desribed in + D. G. A. Smith, L. A. Burns, K. Patkowski, and C. D. Sherrill + J. Phys. Chem. Lett., 2016, 7, pp 2197–2203. + (Functionality should correspond to V3.2 Rev 0 of original dftd3 code.) + +Fixed +----- + +* Routine dftd3_pbc_dispersion delivers now correct stress tensor. + + 0.1 === diff --git a/README.rst b/README.rst index af737ea..f429c77 100644 --- a/README.rst +++ b/README.rst @@ -17,6 +17,9 @@ split into two parts: * Additional extensions which are necessary for the command line tool DFTD3 and the command line tool itself. +* Updated dftd3 code to include refitted/modified zero- and BJ-damped D3 + versions of Sherrill and coworkers (-bjm and -zerom) + (Functionality corresponds to V3.2 Rev0) Compilation =========== diff --git a/doc/man.pdf b/doc/man.pdf index 9f8abbe..da59a9b 100644 Binary files a/doc/man.pdf and b/doc/man.pdf differ diff --git a/lib/api.f90 b/lib/api.f90 index 235645e..c16d1ae 100644 --- a/lib/api.f90 +++ b/lib/api.f90 @@ -215,7 +215,7 @@ end subroutine dftd3_dispersion !! !! \param coords Coordinates of the atoms in atomic units. Shape: [3, nAtom]. !! \param izp Atomic number of each atom. Shape: [nAtom]. You can determine - !! the atomic number using the get_atomic_number() function. + !! the atomic number using the get_atomic_number() function. !! \param latvecs Lattice vectors in atomic units. Shape: [3, 3]. !! \param disp Calculated dispersion energy in atomic units. !! \param grads Calculated gradiens in atomic units, if present. @@ -271,6 +271,10 @@ subroutine dftd3_pbc_dispersion(this, coords, izp, latvecs, disp, grads, & & r2r4, this%r0ab, rcov, s6, s18, rs6, rs8, rs10, alp6, alp8, alp10, & & this%noabc, this%numgrad, this%version, grads, disp2, gnorm, & & stress, latvecs, rep_vdw, rep_cn, this%rthr, .false., this%cn_thr) + ! Note, the stress variable in pbcgdisp contains the *lattice derivatives* + ! on return, so it needs to be converted to obtain the stress tensor. + stress(:,:) = -matmul(stress, transpose(latvecs))& + & / abs(determinant(latvecs)) end subroutine dftd3_pbc_dispersion diff --git a/lib/core.f90 b/lib/core.f90 index 8acd473..c384f28 100644 --- a/lib/core.f90 +++ b/lib/core.f90 @@ -190,6 +190,93 @@ subroutine setfuncpar(func,version,TZ,s6,rs6,s18,rs18,alp) logical TZ ! double hybrid values revised according to procedure in the GMTKN30 pap + + if(version.eq.6)then + s6 =1.0d0 + alp =14.0d0 +! BJ damping with parameters from ... + select case (func) + case ("b2-plyp") + rs6 =0.486434 + s18 =0.672820 + rs18=3.656466 + s6 =0.640000 + case ("b3-lyp") + rs6 =0.278672 + s18 =1.466677 + rs18=4.606311 + case ("b97-d") + rs6 =0.240184 + s18 =1.206988 + rs18=3.864426 + case ("b-lyp") + rs6 =0.448486 + s18 =1.875007 + rs18=3.610679 + case ("b-p") + rs6 =0.821850 + s18 =3.140281 + rs18=2.728151 + case ("pbe") + rs6 =0.012092 + s18 =0.358940 + rs18=5.938951 + case ("pbe0") + rs6 =0.007912 + s18 =0.528823 + rs18=6.162326 + case ("lc-wpbe") + rs6 =0.563761 + s18 =0.906564 + rs18=3.593680 + case DEFAULT + call stoprun( 'functional name unknown' ) + end select + endif + + if(version.eq.5)then + s6 =1.0d0 + alp =14.0d0 +! zero damping with parameters from ... + select case (func) + case ("b2-plyp") + rs6 =1.313134 + s18 =0.717543 + rs18=0.016035 + s6 =0.640000 + case ("b3-lyp") + rs6 =1.338153 + s18 =1.532981 + rs18=0.013988 + case ("b97-d") + rs6 =1.151808 + s18 =1.020078 + rs18=0.035964 + case ("b-lyp") + rs6 =1.279637 + s18 =1.841686 + rs18=0.014370 + case ("b-p") + rs6 =1.233460 + s18 =1.945174 + rs18=0.000000 + case ("pbe") + rs6 =2.340218 + s18 =0.000000 + rs18=0.129434 + case ("pbe0") + rs6 =2.077949 + s18 =0.000081 + rs18=0.116755 + case ("lc-wpbe") + rs6 =1.366361 + s18 =1.280619 + rs18=0.003160 + case DEFAULT + call stoprun( 'functional name unknown' ) + end select + endif + ! DFT-D3 with Becke-Johnson finite-damping, variant 2 with their radii ! SE: Alp is only used in 3-body calculations if (version.eq.4)then @@ -746,12 +833,22 @@ subroutine edisp(max_elem,maxc,n,xyz,iz,c6ab,mxc,r2r4,r0ab,rcov, & !THR if (r2.gt.rthr) cycle r =sqrt(r2) - rr=r0ab(iz(jat),iz(iat))/r + tmp2=r0ab(iz(jat),iz(iat)) + rr=tmp2/r ! damping - tmp=rs6*rr - damp6 =1.d0/( 1.d0+6.d0*tmp**alp6 ) - tmp=rs8*rr - damp8 =1.d0/( 1.d0+6.d0*tmp**alp8 ) + if(version.eq.3)then + ! DFT-D3 zero-damp + tmp=rs6*rr + damp6 =1.d0/( 1.d0+6.d0*tmp**alp6 ) + tmp=rs8*rr + damp8 =1.d0/( 1.d0+6.d0*tmp**alp8 ) + else + ! DFT-D3M zero-damp + tmp=(r/(rs6*tmp2))+rs8*tmp2 + damp6 =1.d0/( 1.d0+6.d0*tmp**(-alp6) ) + tmp=(r/tmp2)+rs8*tmp2 + damp8 =1.d0/( 1.d0+6.d0*tmp**(-alp8) ) + endif ! get C6 call getc6(maxc,max_elem,c6ab,mxc,iz(iat),iz(jat), & & cn(iat),cn(jat),c6) @@ -761,13 +858,13 @@ subroutine edisp(max_elem,maxc,n,xyz,iz,c6ab,mxc,r2r4,r0ab,rcov, & ! r2r4 stored in main as sqrt c8 =3.0d0*c6*r2r4(iz(iat))*r2r4(iz(jat)) - ! DFT-D3 zero-damp - if (version.eq.3)then + ! DFT-D3 zero-damp or DFT-D3 M(zero) + if((version.eq.3).or.(version.eq.5))then e6=e6+c6*damp6/r6 e8=e8+c8*damp8/r8 end if - ! DFT-D3(BJ) - if (version.eq.4)then + ! DFT-D3(BJ) or DFT-D3M(BJ) + if((version.eq.4).or.(version.eq.6))then ! use BJ radius tmp=sqrt(c8/c6) e6=e6+ c6/(r6+(a1*tmp+a2)**6) @@ -971,7 +1068,7 @@ subroutine gdisp(max_elem,maxc,n,xyz,iz,c6ab,mxc,r2r4,r0ab,rcov, & ! 3333333333333333333333333333333333333333333333333333333333333333333333 ! zero damping - elseif (version.eq.3) then + elseif ((version.eq.3).or.(version.eq.5)) then if (echo)write(*,*) 'doing analytical gradient O(N^2) ...' call ncoord(n,rcov,iz,xyz,cn,cn_thr) @@ -1013,22 +1110,43 @@ subroutine gdisp(max_elem,maxc,n,xyz,iz,c6ab,mxc,r2r4,r0ab,rcov, & ! ! Calculates damping functions: - t6 = (r/(rs6*R0))**(-alp6) - damp6 =1.d0/( 1.d0+6.d0*t6 ) - t8 = (r/(rs8*R0))**(-alp8) - damp8 =1.d0/( 1.d0+6.d0*t8 ) + if (version.eq.3) then + t6 = (r/(rs6*R0))**(-alp6) + damp6 =1.d0/( 1.d0+6.d0*t6 ) + t8 = (r/(rs8*R0))**(-alp8) + damp8 =1.d0/( 1.d0+6.d0*t8 ) + + tmp1=s6*6.d0*damp6*C6/r7 + tmp2=s8*6.d0*C6*r42*damp8/r9 + ! d(r^(-6))/d(r_ij) + drij(linij)=drij(linij)-tmp1 & + & -4.d0*tmp2 + + + drij(linij)=drij(linij) & + & +tmp1*alp6*t6*damp6 & + & +3.d0*tmp2*alp8*t8*damp8 + !d(f_dmp)/d(r_ij) + + else ! version.eq.5 + t6 = (r/(rs6*R0)+R0*rs8)**(-alp6) + damp6 =1.d0/( 1.d0+6.d0*t6 ) + t8 = (r/R0+R0*rs8)**(-alp8) + damp8 =1.d0/( 1.d0+6.d0*t8 ) + + tmp1=s6*6.d0*damp6*C6/r7 + tmp2=s8*6.d0*C6*r42*damp8/r9 + ! d(r^(-6))/d(r_ij) + drij(linij)=drij(linij)-tmp1 & + & -4.d0*tmp2 + + + drij(linij)=drij(linij) & + & +tmp1*alp6*t6*damp6*r/(r+rs6*R0*R0*rs8) & + & +3.d0*tmp2*alp8*t8*damp8*r/(r+R0*R0*rs8) + !d(f_dmp)/d(r_ij) - tmp1=s6*6.d0*damp6*C6/r7 - tmp2=s8*6.d0*C6*r42*damp8/r9 - ! d(r^(-6))/d(r_ij) - drij(linij)=drij(linij)-tmp1 & - & -4.d0*tmp2 - - - drij(linij)=drij(linij) & - & +tmp1*alp6*t6*damp6 & - & +3.d0*tmp2*alp8*t8*damp8 - !d(f_dmp)/d(r_ij) + end if if ((.not.noabc).and.(r2.lt.abcthr)) then @@ -1063,11 +1181,11 @@ subroutine gdisp(max_elem,maxc,n,xyz,iz,c6ab,mxc,r2r4,r0ab,rcov, & - ! end if !version 3 + ! end if !version 3+5 ! BJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJBJ ! Becke-Johnson finite damping - elseif (version.eq.4) then + elseif ((version.eq.4).or.(version.eq.6)) then a1 =rs6 a2 =rs8 s8 =s18 @@ -2393,7 +2511,7 @@ subroutine pbcedisp(max_elem,maxc,n,xyz,iz,c6ab,mxc,r2r4,r0ab,& logical noabc integer iat,jat,kat - real(wp) r,r2,r6,r8,tmp,dx,dy,dz,c6,c8,c10,ang,rav + real(wp) r,r2,r6,r8,tmp,dx,dy,dz,c6,c8,c10,ang,rav,R0 real(wp) damp6,damp8,damp10,rr,thr,c9,r42,c12,r10,c14 real(wp) cn(n),rxyz(3),dxyz(3) real(wp) r2ab(n*n),cc6ab(n*n),dmp(n*n),d2(3),t1,t2,t3,tau(3) @@ -2468,7 +2586,7 @@ subroutine pbcedisp(max_elem,maxc,n,xyz,iz,c6ab,mxc,r2r4,r0ab,& - else if (version.eq.3) then + else if ((version.eq.3).or.(version.eq.5)) then ! DFT-D3(zero-damping) call pbcncoord(n,rcov,iz,xyz,cn,lat,rep_cn,crit_cn) @@ -2497,13 +2615,22 @@ subroutine pbcedisp(max_elem,maxc,n,xyz,iz,c6ab,mxc,r2r4,r0ab,& if (r2.gt.rthr) cycle r =sqrt(r2) - rr=r0ab(iz(jat),iz(iat))/r + R0=r0ab(iz(jat),iz(iat)) + rr=R0/r ! damping - tmp=rs6*rr - damp6 =1.d0/( 1.d0+6.d0*tmp**alp6 ) - tmp=rs8*rr - damp8 =1.d0/( 1.d0+6.d0*tmp**alp8 ) - + if(version.eq.3)then + ! DFT-D3 zero-damp + tmp=rs6*rr + damp6 =1.d0/( 1.d0+6.d0*tmp**alp6 ) + tmp=rs8*rr + damp8 =1.d0/( 1.d0+6.d0*tmp**alp8 ) + else + ! DFT-D3M zero-damp + tmp=(r/(rs6*R0))+rs8*R0 + damp6 =1.d0/( 1.d0+6.d0*tmp**(-alp6) ) + tmp=(r/R0)+rs8*R0 + damp8 =1.d0/( 1.d0+6.d0*tmp**(-alp8) ) + endif r6=r2**3 e6 =e6+damp6/r6* c6 @@ -2545,12 +2672,23 @@ subroutine pbcedisp(max_elem,maxc,n,xyz,iz,c6ab,mxc,r2r4,r0ab,& ! cutoff if (r2.gt.rthr) cycle r =sqrt(r2) - rr=r0ab(iz(jat),iz(iat))/r - ! damping - tmp=rs6*rr - damp6 =1.d0/( 1.d0+6.d0*tmp**alp6 ) - tmp=rs8*rr - damp8 =1.d0/( 1.d0+6.d0*tmp**alp8 ) + R0=r0ab(iz(jat),iz(iat)) + rr=R0/r + + ! damping + if(version.eq.3)then + ! DFT-D3 zero-damp + tmp=rs6*rr + damp6 =1.d0/( 1.d0+6.d0*tmp**alp6 ) + tmp=rs8*rr + damp8 =1.d0/( 1.d0+6.d0*tmp**alp8 ) + else + ! DFT-D3M zero-damp + tmp=(r/(rs6*R0))+rs8*R0 + damp6 =1.d0/( 1.d0+6.d0*tmp**(-alp6) ) + tmp=(r/R0)+rs8*R0 + damp8 =1.d0/( 1.d0+6.d0*tmp**(-alp8) ) + endif r6=r2**3 @@ -2569,7 +2707,7 @@ subroutine pbcedisp(max_elem,maxc,n,xyz,iz,c6ab,mxc,r2r4,r0ab,& end do end do ! write(*,*)'counter(edisp): ',counter - else if (version.eq.4) then + else if((version.eq.4).or.(version.eq.6)) then ! DFT-D3(BJ-damping) @@ -3257,7 +3395,7 @@ subroutine pbcgdisp(max_elem,maxc,n,xyz,iz,c6ab,mxc,r2r4,r0ab,& goto 999 end if - if (version.eq.3) then + if ((version.eq.3).or.(version.eq.5)) then !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! begin ZERO DAMPING GRADIENT @@ -3319,21 +3457,38 @@ subroutine pbcgdisp(max_elem,maxc,n,xyz,iz,c6ab,mxc,r2r4,r0ab,& ! ! Calculates damping functions: - t6 = (r/(rs6*R0))**(-alp6) - damp6 =1.d0/( 1.d0+6.d0*t6 ) - t8 = (r/(rs8*R0))**(-alp8) - damp8 =1.d0/( 1.d0+6.d0*t8 ) - - drij(tauz,tauy,taux,lin(iat,iat))=drij(tauz,tauy,taux,lin(iat,& - & iat))& - & +(-s6*(6.0/(r7)*C6*damp6)& - & -s8*(24.0/(r9)*C6*r42*damp8))*0.5d0 - - - drij(tauz,tauy,taux,lin(iat,iat))=drij(tauz,tauy,taux,lin(iat,& - & iat))& - & +(s6*C6/r7*6.d0*alp6*t6*damp6*damp6& - & +s8*C6*r42/r9*18.d0*alp8*t8*damp8*damp8)*0.5d0 + if (version.eq.3) then + t6 = (r/(rs6*R0))**(-alp6) + damp6 =1.d0/( 1.d0+6.d0*t6 ) + t8 = (r/(rs8*R0))**(-alp8) + damp8 =1.d0/( 1.d0+6.d0*t8 ) + + drij(tauz,tauy,taux,lin(iat,iat))=drij(tauz,tauy,taux,lin(iat,& + & iat))& + & +(-s6*(6.0/(r7)*C6*damp6)& + & -s8*(24.0/(r9)*C6*r42*damp8))*0.5d0 + + + drij(tauz,tauy,taux,lin(iat,iat))=drij(tauz,tauy,taux,lin(iat,& + & iat))& + & +(s6*C6/r7*6.d0*alp6*t6*damp6*damp6& + & +s8*C6*r42/r9*18.d0*alp8*t8*damp8*damp8)*0.5d0 + else !version.eq.5 + t6 = (r/(rs6*R0)+R0*rs8)**(-alp6) + damp6 =1.d0/( 1.d0+6.d0*t6 ) + t8 = (r/(R0)+R0*rs8)**(-alp8) + damp8 =1.d0/( 1.d0+6.d0*t8 ) + + tmp1=s6*6.d0*damp6*C6/r7 + tmp2=s8*6.d0*C6*r42*damp8/r9 + drij(tauz,tauy,taux,lin(iat,iat))=drij(tauz,tauy,taux,lin(iat, & + & iat)) - (tmp1 +4.d0*tmp2)*0.5d0 ! d(r^(-6))/d(r_ij) + + + drij(tauz,tauy,taux,lin(iat,iat))=drij(tauz,tauy,taux,lin(iat, & + & iat)) +(tmp1*alp6*t6*damp6*r/(r+rs6*R0*R0*rs8) & + & +3.d0*tmp2*alp8*t8*damp8*r/(r+R0*R0*rs8))*0.5d0 !d(f_dmp)/d(r_ij) + endif ! ! in dC6_rest all terms BUT C6-term is saved for the kat-loop ! @@ -3396,20 +3551,37 @@ subroutine pbcgdisp(max_elem,maxc,n,xyz,iz,c6ab,mxc,r2r4,r0ab,& ! ! Calculates damping functions: - t6 = (r/(rs6*R0))**(-alp6) - damp6 =1.d0/( 1.d0+6.d0*t6 ) - t8 = (r/(rs8*R0))**(-alp8) - damp8 =1.d0/( 1.d0+6.d0*t8 ) - - drij(tauz,tauy,taux,linij)=drij(tauz,tauy,taux,& - & linij)& - & -s6*(6.0/(r7)*C6*damp6)& - & -s8*(24.0/(r9)*C6*r42*damp8) - - drij(tauz,tauy,taux,linij)=drij(tauz,tauy,taux,& - & linij)& - & +s6*C6/r7*6.d0*alp6*t6*damp6*damp6& - & +s8*C6*r42/r9*18.d0*alp8*t8*damp8*damp8 + if (version.eq.3) then + t6 = (r/(rs6*R0))**(-alp6) + damp6 =1.d0/( 1.d0+6.d0*t6 ) + t8 = (r/(rs8*R0))**(-alp8) + damp8 =1.d0/( 1.d0+6.d0*t8 ) + + drij(tauz,tauy,taux,linij)=drij(tauz,tauy,taux,& + & linij)& + & -s6*(6.0/(r7)*C6*damp6)& + & -s8*(24.0/(r9)*C6*r42*damp8) + + drij(tauz,tauy,taux,linij)=drij(tauz,tauy,taux,& + & linij)& + & +s6*C6/r7*6.d0*alp6*t6*damp6*damp6& + & +s8*C6*r42/r9*18.d0*alp8*t8*damp8*damp8 + else !version.eq.5 + t6 = (r/(rs6*R0)+R0*rs8)**(-alp6) + damp6 =1.d0/( 1.d0+6.d0*t6 ) + t8 = (r/(R0)+R0*rs8)**(-alp8) + damp8 =1.d0/( 1.d0+6.d0*t8 ) + + tmp1=s6*6.d0*damp6*C6/r7 + tmp2=s8*6.d0*C6*r42*damp8/r9 + drij(tauz,tauy,taux,linij)=drij(tauz,tauy,taux, & + & linij) - (tmp1 +4.d0*tmp2) ! d(r^(-6))/d(r_ij) + + + drij(tauz,tauy,taux,linij)=drij(tauz,tauy,taux,linij) & + & +(tmp1*alp6*t6*damp6*r/(r+rs6*R0*R0*rs8) & + & +3.d0*tmp2*alp8*t8*damp8*r/(r+R0*R0*rs8)) !d(f_dmp)/d(r_ij) + endif ! ! in dC6_rest all terms BUT C6-term is saved for the kat-loop ! @@ -3434,7 +3606,7 @@ subroutine pbcgdisp(max_elem,maxc,n,xyz,iz,c6ab,mxc,r2r4,r0ab,& end do - elseif (version.eq.4) then + elseif ((version.eq.4).or.(version.eq.6)) then @@ -4785,16 +4957,25 @@ end function VECTORSIZE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + function determinant(aa) result(det) + real(wp), intent(in) :: aa(:,:) + real(wp) :: det + + det = aa(1,1) * (aa(2,2) * aa(3,3) - aa(3,2) * aa(2,3))& + & - aa(1,2) * (aa(2,1) * aa(3,3) - aa(3,1) * aa(2,3))& + & + aa(1,3) * (aa(2,1) * aa(3,2) - aa(3,1) * aa(2,2)) + + end function determinant + + subroutine inv_cell(x,a) real(wp), intent(in) :: x(3,3) real(wp), intent(out) :: a(3,3) integer i real(wp) det - a=0.0 - det=x(1,1)*x(2,2)*x(3,3)+x(1,2)*x(2,3)*x(3,1)+x(1,3)*x(2,1)*& - & x(3,2)-x(1,3)*x(2,2)*x(3,1)-x(1,2)*x(2,1)*x(3,3)-x(1,1)*& - & x(2,3)*x(3,2) + a = 0.0 + det = determinant(x) ! write(*,*)'Det:',det a(1,1)=x(2,2)*x(3,3)-x(2,3)*x(3,2) a(2,1)=x(2,3)*x(3,1)-x(2,1)*x(3,3) diff --git a/make.arch b/make.arch index dfbff81..892871d 100644 --- a/make.arch +++ b/make.arch @@ -1,7 +1,7 @@ # Make your choice -ARCH = x86_64-linux-gcc +ARCH = x86_64-linux-gnu -ifeq ($(ARCH),x86_64-linux-gcc) +ifeq ($(ARCH),x86_64-linux-gnu) FC = gfortran LN = gfortran FCFLAGS = diff --git a/prg/extras.f90 b/prg/extras.f90 index 866c609..5e9c0a1 100644 --- a/prg/extras.f90 +++ b/prg/extras.f90 @@ -55,6 +55,8 @@ subroutine printoptions write(*,*) '-old (DFT-D2)' write(*,*) '-zero (DFT-D3 original zero-damping)' write(*,*) '-bj (DFT-D3 with Becke-Johnson finite-damping)' + write(*,*) '-zerom (revised DFT-D3 original zero-damping)' + write(*,*) '-bjm (revised DFT-D3 with Becke-Johnson damping)' write(*,*) '-tz (use special parameters for TZ-type calculations)' write(*,*) 'variable parameters can be read from /.dftd3par.local' @@ -213,10 +215,19 @@ subroutine adisp(max_elem,maxc,n,xyz,iz,c6ab,mxc,r2r4,r0ab,rcov, & rr=R0/r r6=r2**3 - tmp=rs6*rr - damp6 =1.d0/( 1.d0+6.d0*tmp**alp6 ) - tmp=rs8*rr - damp8 =1.d0/( 1.d0+6.d0*tmp**alp8 ) + if(version.eq.3)then + ! DFT-D3 zero-damp + tmp=rs6*rr + damp6 =1.d0/( 1.d0+6.d0*tmp**alp6 ) + tmp=rs8*rr + damp8 =1.d0/( 1.d0+6.d0*tmp**alp8 ) + else + ! DFT-D3M zero-damp + tmp=(r/(rs6*R0))+rs8*R0 + damp6 =1.d0/( 1.d0+6.d0*tmp**(-alp6) ) + tmp=(r/R0)+rs8*R0 + damp8 =1.d0/( 1.d0+6.d0*tmp**(-alp8) ) + endif if (version.eq.2)then c6=c6ab(iz(jat),iz(iat),1,1,1) @@ -228,7 +239,7 @@ subroutine adisp(max_elem,maxc,n,xyz,iz,c6ab,mxc,r2r4,r0ab,rcov, & & cn(iat),cn(jat),c6) end if - if (version.eq.3)then + if((version.eq.3).or.(version.eq.5))then e6 =s6*autokcal*c6*damp6/r6 r8 =r6*r2 r42=r2r4(iz(iat))*r2r4(iz(jat)) @@ -236,7 +247,7 @@ subroutine adisp(max_elem,maxc,n,xyz,iz,c6ab,mxc,r2r4,r0ab,rcov, & e8 =s18*autokcal*c8*damp8/r8 end if - if (version.eq.4)then + if((version.eq.4).or.(version.eq.6))then r42=r2r4(iz(iat))*r2r4(iz(jat)) c8 =3.0d0*c6*r42 ! use BJ radius @@ -1650,10 +1661,17 @@ subroutine pbcadisp(max_elem,maxc,n,xyz,iz,c6ab,mxc,r2r4,r0ab,& rr=R0/r r6=r2**3 - tmp=rs6*rr - damp6 =1.d0/( 1.d0+6.d0*tmp**alp6 ) - tmp=rs8*rr - damp8 =1.d0/( 1.d0+6.d0*tmp**alp8 ) + if(version.eq.3)then + tmp=rs6*rr + damp6 =1.d0/( 1.d0+6.d0*tmp**alp6 ) + tmp=rs8*rr + damp8 =1.d0/( 1.d0+6.d0*tmp**alp8 ) + else + tmp=(r/(R0*rs6)+R0*rs8)**(-alp6) + damp6 =1.d0/( 1.d0+6.d0*tmp ) + tmp=(r/(R0)+R0*rs8)**(-alp8) + damp8 =1.d0/( 1.d0+6.d0*tmp ) + endif if (version.eq.2)then c6=c6ab(iz(jat),iz(iat),1,1,1) @@ -1669,7 +1687,7 @@ subroutine pbcadisp(max_elem,maxc,n,xyz,iz,c6ab,mxc,r2r4,r0ab,& & cn(iat),cn(jat),c6) end if - if (version.eq.3)then + if((version.eq.3).or.(version.eq.5))then r8 =r6*r2 r42=r2r4(iz(iat))*r2r4(iz(jat)) c8 =3.0d0*c6*r42 @@ -1682,7 +1700,7 @@ subroutine pbcadisp(max_elem,maxc,n,xyz,iz,c6ab,mxc,r2r4,r0ab,& end if end if - if (version.eq.4)then + if((version.eq.4).or.(version.eq.6))then r42=r2r4(iz(iat))*r2r4(iz(jat)) c8 =3.0d0*c6*r42 ! use BJ radius diff --git a/prg/main.f90 b/prg/main.f90 index 13880f9..bfc3e82 100644 --- a/prg/main.f90 +++ b/prg/main.f90 @@ -180,6 +180,8 @@ program dftd3_main if (index(ftmp,'-old') .ne.0) version=2 if (index(ftmp,'-zero') .ne.0) version=3 if (index(ftmp,'-bj') .ne.0) version=4 + if (index(ftmp,'-zerom') .ne.0) version=5 + if (index(ftmp,'-bjm') .ne.0) version=6 if (index(ftmp,'-min') .ne.0) then minc6=.true. j=0 @@ -331,13 +333,16 @@ program dftd3_main write(*,*)'J. Comput. Chem. 32 (2011), 1456-1465' write(*,*)'For DFT-D2 the reference is' write(*,*)'S. Grimme, J. Comput. Chem., 27 (2006), 1787-1799' + write(*,*)'For DFT-D3M or DFT-D3M(BJ) the reference is' + write(*,*)'D.G.A. Smith, L.A. Burns, K. Patkowski, and ' + write(*,*)'C.D. Sherrill, J. Phys. Chem. Lett. 7 (2016) 2197-2203' write(*,*) write(*,*)' files read : ' write(*,*)trim(etmp) if (.not.ex)write(*,*)trim(dtmp) end if - if (version.lt.2.or.version.gt.4)stop 'inacceptable version number' + if (version.lt.2.or.version.gt.6)stop 'inacceptable version number' !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ! all calculations start here @@ -400,7 +405,9 @@ program dftd3_main c8 =r2r4(iz(i))**2*3.0d0*c6 c10=(49.0d0/40.0d0)*c8**2/c6 dum=0.5*autoang*r0ab(z,z) - if (version.eq.4)dum=rs6*0.5*autoang*sqrt(c8/c6) + if((version.eq.4).or.(version.eq.6))then + dum=rs6*0.5*autoang*sqrt(c8/c6) + endif atmp=' ' if (fix(i)) then atmp='f' @@ -555,11 +562,16 @@ program dftd3_main ! output if (echo) then + if(version.lt.4)then write(*,'(/10x,'' DFT-D V'',i1)') version - else - write(*,'(/10x,'' DFT-D V3(BJ)'')') - end if + elseif(version.eq.4)then + write(*,'(/10x,'' DFT-D V3(BJ)'')') + elseif(version.eq.5)then + write(*,'(/10x,'' DFT-D V3 M'')') + elseif(version.eq.6)then + write(*,'(/10x,'' DFT-D V3 M(BJ)'')') + endif write(*,'('' DF '',a50)') func write(*,'('' parameters'')') if(version.eq.2)then @@ -575,13 +587,22 @@ program dftd3_main write(*,'('' alpha8 :'',f10.4)') alp8 write(*,'('' k1-k3 :'',3f10.4)') k1,k2,k3 end if - if(version.eq.4)then - write(*,'('' s6 :'',F23.15)') s6 - write(*,'('' s8 :'',F23.15)') s18 - write(*,'('' a1 :'',F23.15)') rs6 - write(*,'('' a2 :'',F23.15)') rs18 - write(*,'('' k1-k3 :'',3F23.15)') k1,k2,k3 + if((version.eq.4).or.(version.eq.6))then + write(*,'('' s6 :'',F10.4)') s6 + write(*,'('' s8 :'',F10.4)') s18 + write(*,'('' a1 :'',F10.4)') rs6 + write(*,'('' a2 :'',F10.4)') rs18 + write(*,'('' k1-k3 :'',3F10.4)') k1,k2,k3 end if + if(version.eq.5)then + write(*,'('' s6 :'',f10.4)') s6 + write(*,'('' s8 :'',f10.4)') s18 + write(*,'('' rs6 :'',f10.4)') rs6 + write(*,'('' beta :'',f10.4)') rs18 + write(*,'('' alpha6 :'',f10.4)') alp6 + write(*,'('' alpha8 :'',f10.4)') alp8 + write(*,'('' k1-k3 :'',3f10.4)') k1,k2,k3 + endif write(*,'('' Cutoff :'',f10.4,'' a.u.'')') sqrt(rthr) !*autoang write(*,'('' CN-Cutoff:'',f10.4,'' a.u.'')')sqrt(cn_thr)!*autoang ! if (pbc) then