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 d44e78f..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
@@ -12659,23 +12682,37 @@ C for vectorization reasons we will sumup at the end to avoid depenence of previ
        sc_aa_tube=sc_aa_tube_par(iti)
        sc_bb_tube=sc_bb_tube_par(iti)
        enetube(i+nres)=sc_aa_tube/rdiff6**2.0d0+sc_bb_tube/rdiff6
+C       enetube(i+nres)=0.0d0
 C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6
 C now we calculate gradient
        fac=-12.0d0*sc_aa_tube/rdiff6**2.0d0/rdiff-
      &       6.0d0*sc_bb_tube/rdiff6/rdiff
+C       fac=0.0
 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) go to 667
+         if (acavtub(iti).eq.0.0d0) then
+C go to 667
+         enecavtube(i+nres)=0.0
+         faccav=0.0
+         else
          denominator=(1.0+dcavtub(iti)*rdiff6*rdiff6)
-         enecavtube(i)=
-     &   acavtub(iti)*(rdiff+bcavtub(iti)*sqrt(rdiff)+ccavtub(iti))
+         enecavtube(i+nres)=
+     &   (bcavtub(iti)*rdiff+acavtub(iti)*sqrt(rdiff)+ccavtub(iti))
      &   /denominator
-         faccav=(acavtub(iti)*(1.0+bcavtub(iti)/2.0/sqrt(rdiff))
-     &   *denominator-acavtub(iti)*(rdiff+bcavtub(iti)*sqrt(rdiff)
-     &   +ccavtub(iti))*rdiff6**2.0d0/rdiff*dcavtub(iti))
+C         enecavtube(i)=0.0
+         faccav=((bcavtub(iti)*1.0d0+acavtub(iti)/2.0d0/sqrt(rdiff))
+     &   *denominator-(bcavtub(iti)*rdiff+acavtub(iti)*sqrt(rdiff)
+     &   +ccavtub(iti))*rdiff6**2.0d0/rdiff*dcavtub(iti)*12.0d0)
      &   /denominator**2.0d0
+C         faccav=0.0
          fac=fac+faccav
- 667     continue
+C 667     continue
+         endif
+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
@@ -12690,7 +12727,8 @@ C        if (acavtub(iti).eq.0.0) cycle
 
 
         do i=itube_start,itube_end
-          Etube=Etube+enetube(i)+enetube(i+nres)
+          Etube=Etube+enetube(i)+enetube(i+nres)+enecavtube(i)
+     & +enecavtube(i+nres)
         enddo
 C        print *,"ETUBE", etube
         return