changes needed for respa for evdw

diff --git a/source/unres/energy.f90 b/source/unres/energy.f90
index 1d58136..e360fd8 100644 (file)
--- a/source/unres/energy.f90
+++ b/source/unres/energy.f90
@@ -10999,6 +10999,20 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
       endif
       return
       end function sscale
+      real(kind=8) function sscale_grad(r)
+!      include "COMMON.SPLITELE"
+      real(kind=8) :: r,gamm
+      if(r.lt.r_cut-rlamb) then
+        sscale_grad=0.0d0
+      else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
+        gamm=(r-(r_cut-rlamb))/rlamb
+        sscale_grad=gamm*(6*gamm-6.0d0)/rlamb
+      else
+        sscale_grad=0d0
+      endif
+      return
+      end function sscale_grad
+
 !!!!!!!!!! PBCSCALE
       real(kind=8) function sscale_ele(r)
 !      include "COMMON.SPLITELE"
@@ -11648,7 +11662,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
 !el local variables
       integer :: iint,itypi,itypi1,itypj,subchap
       real(kind=8) :: rrij,xi,yi,zi,fac,sigm,epsi,sig,sig0ij,rij_shift
-      real(kind=8) :: sss,e1,e2,evdw
+      real(kind=8) :: sss,e1,e2,evdw,sss_grad
       real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
                     dist_temp, dist_init
 
@@ -11752,6 +11766,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
             sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
             sss_ele_cut=sscale_ele(1.0d0/(rij*sigma(itypi,itypj)))
             sss_ele_grad=sscagrad_ele(1.0d0/(rij*sigma(itypi,itypj)))
+            sss_grad=sscale_grad(1.0d0/(rij*sigmaii(itypi,itypj)))
             if (sss_ele_cut.le.0.0) cycle
             if (sss.lt.1.0d0) then
 
@@ -11805,8 +11820,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
               fac=-expon*(e1+evdwij)*rij_shift
               sigder=fac*sigder
               fac=rij*fac
-              fac=fac+evdwij*sss_ele_grad/sss_ele_cut&
-            /sigma(itypi,itypj)*rij
+              fac=fac+evdwij*(sss_ele_grad/sss_ele_cut&
+            /sigma(itypi,itypj)*rij-sss_grad/(1.0-sss)*rij  &
+            /sigmaii(itypi,itypj))
 !              fac=0.0d0
 ! Calculate the radial part of the gradient
               gg(1)=xj*fac
@@ -11845,7 +11861,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
 !el local variables
       integer :: iint,itypi,itypi1,itypj,subchap
       real(kind=8) :: rrij,xi,yi,zi,fac,sigm,epsi,sig0ij,sig
-      real(kind=8) :: sss,e1,e2,evdw,rij_shift
+      real(kind=8) :: sss,e1,e2,evdw,rij_shift,sss_grad
       real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
                     dist_temp, dist_init
       evdw=0.0D0
@@ -11955,6 +11971,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
             rij=dsqrt(rrij)
             sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
+            sss_grad=sscale_grad(1.0d0/(rij*sigmaii(itypi,itypj)))
             sss_ele_cut=sscale_ele(1.0d0/(rij*sigma(itypi,itypj)))
             sss_ele_grad=sscagrad_ele(1.0d0/(rij*sigma(itypi,itypj)))
             if (sss_ele_cut.le.0.0) cycle
@@ -11989,7 +12006,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
 !              write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
 !     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
               evdwij=evdwij*eps2rt*eps3rt
-              evdw=evdw+evdwij*sss
+              evdw=evdw+evdwij*sss*sss_ele_cut
               if (lprn) then
               sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
               epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
@@ -12011,8 +12028,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
               fac=-expon*(e1+evdwij)*rij_shift
               sigder=fac*sigder
               fac=rij*fac
-              fac=fac+evdwij*sss_ele_grad/sss_ele_cut&
-            /sigma(itypi,itypj)*rij
+              fac=fac+evdwij*(sss_ele_grad/sss_ele_cut&
+            /sigma(itypi,itypj)*rij+sss_grad/sss*rij  &
+            /sigmaii(itypi,itypj))
 
 !              fac=0.0d0
 ! Calculate the radial part of the gradient