corrected cluster for lipid

[unres.git] / source / cluster / wham / src-M / energy_p_new.F

diff --git a/source/cluster/wham/src-M/energy_p_new.F b/source/cluster/wham/src-M/energy_p_new.F
index dfdaba6..206d327 100644 (file)
--- a/source/cluster/wham/src-M/energy_p_new.F
+++ b/source/cluster/wham/src-M/energy_p_new.F
@@ -2237,6 +2237,31 @@ C         endif
           if (ymedi.lt.0) ymedi=ymedi+boxysize
           zmedi=mod(zmedi,boxzsize)
           if (zmedi.lt.0) zmedi=zmedi+boxzsize
+          zmedi2=mod(zmedi,boxzsize)
+          if (zmedi2.lt.0) zmedi2=zmedi2+boxzsize
+       if ((zmedi2.gt.bordlipbot)
+     &.and.(zmedi2.lt.bordliptop)) then
+C the energy transfer exist
+        if (zmedi2.lt.buflipbot) then
+C what fraction I am in
+         fracinbuf=1.0d0-
+     &        ((zmedi2-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+         sslipi=sscalelip(fracinbuf)
+         ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
+        elseif (zmedi2.gt.bufliptop) then
+         fracinbuf=1.0d0-((bordliptop-zmedi2)/lipbufthick)
+         sslipi=sscalelip(fracinbuf)
+         ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
+        else
+         sslipi=1.0d0
+         ssgradlipi=0.0d0
+        endif
+       else
+         sslipi=0.0d0
+         ssgradlipi=0.0d0
+       endif
+
         num_conti=0
 c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
         do j=ielstart(i),ielend(i)
@@ -2280,6 +2305,29 @@ C End diagnostics
           if (yj.lt.0) yj=yj+boxysize
           zj=mod(zj,boxzsize)
           if (zj.lt.0) zj=zj+boxzsize
+       if ((zj.gt.bordlipbot)
+     &.and.(zj.lt.bordliptop)) then
+C the energy transfer exist
+        if (zj.lt.buflipbot) then
+C what fraction I am in
+         fracinbuf=1.0d0-
+     &        ((zj-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+         sslipj=sscalelip(fracinbuf)
+         ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
+        elseif (zj.gt.bufliptop) then
+         fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick)
+         sslipj=sscalelip(fracinbuf)
+         ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
+        else
+         sslipj=1.0d0
+         ssgradlipj=0.0
+        endif
+       else
+         sslipj=0.0d0
+         ssgradlipj=0.0
+       endif
+
       dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
       xj_safe=xj
       yj_safe=yj
@@ -2349,15 +2397,18 @@ C#undef DEBUG
           el1=el1*fac_shield(i)**2*fac_shield(j)**2
           el2=el2*fac_shield(i)**2*fac_shield(j)**2
           eesij=(el1+el2)
+C     &*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
           ees=ees+eesij
           else
           fac_shield(i)=1.0
           fac_shield(j)=1.0
           eesij=(el1+el2)
+     &*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
           ees=ees+eesij
           endif
 C          ees=ees+eesij
           evdw1=evdw1+evdwij*sss
+     &*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
 cd          write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
 cd     &      iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
 cd     &      1.0D0/dsqrt(rrmij),evdwij,eesij,
@@ -2788,6 +2839,8 @@ C           fac_shield(j)=0.6
           endif
           eel_loc_ij=eel_loc_ij
      &    *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
           eel_loc=eel_loc+eel_loc_ij
 C Partial derivatives in virtual-bond dihedral angles gamma
           if (calc_grad) then
@@ -2852,9 +2905,10 @@ cd          write(iout,*) 'aggj1',aggj1
 
 C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
           do l=1,3
-            ggg(l)=agg(l,1)*muij(1)+
-     &          agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4)
+            ggg(l)=(agg(l,1)*muij(1)+
+     &          agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4))
      &    *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
 
           enddo
           do k=i+2,j2
@@ -2867,18 +2921,22 @@ C Remaining derivatives of eello
             gel_loc(l,i)=gel_loc(l,i)+aggi(l,1)*muij(1)+
      &          aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4)
      &    *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
 
             gel_loc(l,i+1)=gel_loc(l,i+1)+aggi1(l,1)*muij(1)+
      &          aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4)
      &    *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
 
             gel_loc(l,j)=gel_loc(l,j)+aggj(l,1)*muij(1)+
      &          aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4)
      &    *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
 
             gel_loc(l,j1)=gel_loc(l,j1)+aggj1(l,1)*muij(1)+
      &          aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4)
      &    *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
 
           enddo
           endif
@@ -3178,8 +3236,11 @@ C        fac_shield(j)=0.6
         endif
         eello_turn3=eello_turn3+0.5d0*(pizda(1,1)+pizda(2,2))
      &  *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+
         eello_t3=0.5d0*(pizda(1,1)+pizda(2,2))
      &  *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
 
 cd        write (2,*) 'i,',i,' j',j,'eello_turn3',
 cd     &    0.5d0*(pizda(1,1)+pizda(2,2)),
@@ -3253,6 +3314,7 @@ C Cartesian derivatives
           gcorr3_turn(l,i)=gcorr3_turn(l,i)
      &      +0.5d0*(pizda(1,1)+pizda(2,2))
      &   *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
 
           a_temp(1,1)=aggi1(l,1)
           a_temp(1,2)=aggi1(l,2)
@@ -3262,6 +3324,7 @@ C Cartesian derivatives
           gcorr3_turn(l,i+1)=gcorr3_turn(l,i+1)
      &      +0.5d0*(pizda(1,1)+pizda(2,2))
      &   *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
 
           a_temp(1,1)=aggj(l,1)
           a_temp(1,2)=aggj(l,2)
@@ -3271,6 +3334,7 @@ C Cartesian derivatives
           gcorr3_turn(l,j)=gcorr3_turn(l,j)
      &      +0.5d0*(pizda(1,1)+pizda(2,2))
      &   *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
 
           a_temp(1,1)=aggj1(l,1)
           a_temp(1,2)=aggj1(l,2)
@@ -3280,6 +3344,7 @@ C Cartesian derivatives
           gcorr3_turn(l,j1)=gcorr3_turn(l,j1)
      &      +0.5d0*(pizda(1,1)+pizda(2,2))
      &   *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
 
         enddo
         endif
@@ -3333,8 +3398,11 @@ C        fac_shield(j)=0.6
         endif
         eello_turn4=eello_turn4-(s1+s2+s3)
      &  *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+
         eello_t4=-(s1+s2+s3)
      &  *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
 
 cd        write (2,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3),
 cd     &    ' eello_turn4_num',8*eello_turn4_num
@@ -3390,6 +3458,7 @@ C     &     *2.0
         s3=0.5d0*(pizda(1,1)+pizda(2,2))
         gel_loc_turn4(i)=gel_loc_turn4(i)-(s1+s3)
      &  *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
 
 C Derivatives in gamma(i+1)
         call transpose2(EUgder(1,1,i+2),e2tder(1,1))
@@ -3400,6 +3469,7 @@ C Derivatives in gamma(i+1)
         s3=0.5d0*(pizda(1,1)+pizda(2,2))
         gel_loc_turn4(i+1)=gel_loc_turn4(i+1)-(s2+s3)
      &  *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
 
 C Derivatives in gamma(i+2)
         call transpose2(EUgder(1,1,i+3),e3tder(1,1))
@@ -3413,7 +3483,8 @@ C Derivatives in gamma(i+2)
         s3=0.5d0*(pizda(1,1)+pizda(2,2))
         gel_loc_turn4(i+2)=gel_loc_turn4(i+2)-(s1+s2+s3)
      &  *fac_shield(i)*fac_shield(j)
-
+     &*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+ 
 C Cartesian derivatives
 C Derivatives of this turn contributions in DC(i+2)
         if (j.lt.nres-1) then
@@ -3434,6 +3505,7 @@ C Derivatives of this turn contributions in DC(i+2)
             ggg(l)=-(s1+s2+s3)
             gcorr4_turn(l,i+2)=gcorr4_turn(l,i+2)-(s1+s2+s3)
      &  *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
 
           enddo
         endif
@@ -3454,6 +3526,7 @@ C Remaining derivatives of this turn contribution
           s3=0.5d0*(pizda(1,1)+pizda(2,2))
           gcorr4_turn(l,i)=gcorr4_turn(l,i)-(s1+s2+s3)
      &  *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
 
           a_temp(1,1)=aggi1(l,1)
           a_temp(1,2)=aggi1(l,2)
@@ -3470,6 +3543,7 @@ C Remaining derivatives of this turn contribution
           s3=0.5d0*(pizda(1,1)+pizda(2,2))
           gcorr4_turn(l,i+1)=gcorr4_turn(l,i+1)-(s1+s2+s3)
      &  *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
 
           a_temp(1,1)=aggj(l,1)
           a_temp(1,2)=aggj(l,2)
@@ -3486,6 +3560,7 @@ C Remaining derivatives of this turn contribution
           s3=0.5d0*(pizda(1,1)+pizda(2,2))
           gcorr4_turn(l,j)=gcorr4_turn(l,j)-(s1+s2+s3)
      &  *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
 
           a_temp(1,1)=aggj1(l,1)
           a_temp(1,2)=aggj1(l,2)
@@ -3502,6 +3577,7 @@ C Remaining derivatives of this turn contribution
           s3=0.5d0*(pizda(1,1)+pizda(2,2))
           gcorr4_turn(l,j1)=gcorr4_turn(l,j1)-(s1+s2+s3)
      &  *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
 
         enddo
         endif