introcuction of cavity term for pepgroup for nanosphere

[unres.git] / source / unres / src_MD-M / energy_p_new_barrier.F

diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F
index c037ff3..91a7c94 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new_barrier.F
@@ -12589,7 +12589,30 @@ C now we calculate gradient
      &       6.0d0*pep_bb_tube)/rdiff6/rdiff
 C       write(iout,'(a5,i4,f12.1,3f12.5)') "TU13",i,rdiff6,enetube(i),
 C     &rdiff,fac
-
+         if (acavtubpep.eq.0.0d0) then
+C go to 667
+         enecavtube(i)=0.0
+         faccav=0.0
+         else
+         denominator=(1.0+dcavtubpep*rdiff6*rdiff6)
+         enecavtube(i)=
+     &   (bcavtubpep*rdiff+acavtubpep*sqrt(rdiff)+ccavtubpep)
+     &   /denominator
+         enecavtube(i)=0.0
+         faccav=((bcavtubpep*1.0d0+acavtubpep/2.0d0/sqrt(rdiff))
+     &   *denominator-(bcavtubpep*rdiff+acavtubpep*sqrt(rdiff)
+     &   +ccavtubpep)*rdiff6**2.0d0/rdiff*dcavtubpep*12.0d0)
+     &   /denominator**2.0d0
+C         faccav=0.0
+C         fac=fac+faccav
+C 667     continue
+         endif
+C         print *,"TUT",i,iti,rdiff,rdiff6,acavtubpep,denominator,
+C     &   enecavtube(i),faccav
+C         print *,"licz=",
+C     & (bcavtub(iti)*rdiff+acavtub(iti)*sqrt(rdiff)+ccavtub(iti))
+CX         print *,"finene=",enetube(i+nres)+enecavtube(i)
+         
 C now direction of gg_tube vector
         do j=1,3
         gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac/2.0d0
@@ -12669,11 +12692,11 @@ C now direction of gg_tube vector
 C Now cavity term E=a(x+bsqrt(x)+c)/(1+dx^12)
          if (acavtub(iti).eq.0.0d0) then
 C go to 667
-         enecavtube(i)=0.0
+         enecavtube(i+nres)=0.0
          faccav=0.0
          else
          denominator=(1.0+dcavtub(iti)*rdiff6*rdiff6)
-         enecavtube(i)=
+         enecavtube(i+nres)=
      &   (bcavtub(iti)*rdiff+acavtub(iti)*sqrt(rdiff)+ccavtub(iti))
      &   /denominator
 C         enecavtube(i)=0.0
@@ -12685,11 +12708,11 @@ C         faccav=0.0
          fac=fac+faccav
 C 667     continue
          endif
-         print *,"TUT",i,iti,rdiff,rdiff6,acavtub(iti),denominator,
-     &   enecavtube(i),faccav
-         print *,"licz=",
-     & (bcavtub(iti)*rdiff+acavtub(iti)*sqrt(rdiff)+ccavtub(iti))
-         print *,"finene=",enetube(i+nres)+enecavtube(i)
+C         print *,"TUT",i,iti,rdiff,rdiff6,acavtub(iti),denominator,
+C     &   enecavtube(i),faccav
+C         print *,"licz=",
+C     & (bcavtub(iti)*rdiff+acavtub(iti)*sqrt(rdiff)+ccavtub(iti))
+C         print *,"finene=",enetube(i+nres)+enecavtube(i)
          do j=1,3
           gg_tube_SC(j,i)=gg_tube_SC(j,i)+vectube(j)*fac
           gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac
@@ -12705,6 +12728,7 @@ C        if (acavtub(iti).eq.0.0) cycle
 
         do i=itube_start,itube_end
           Etube=Etube+enetube(i)+enetube(i+nres)+enecavtube(i)
+     & +enecavtube(i+nres)
         enddo
 C        print *,"ETUBE", etube
         return