projects / unres.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
introduce lipid good grad wrong symplex
[unres.git] / source / unres / src_MD-M / energy_p_new-sep_barrier.F
diff --git a/source/unres/src_MD-M/energy_p_new-sep_barrier.F b/source/unres/src_MD-M/energy_p_new-sep_barrier.F
index bba2354..c0b4c84 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new-sep_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new-sep_barrier.F
@@ -104,8 +104,8 @@ cd   &                  'iend=',iend(i,iint)
               rrij=1.0D0/rij
               eps0ij=eps(itypi,itypj)
               fac=rrij**expon2
-              e1=fac*fac*aa(itypi,itypj)
-              e2=fac*bb(itypi,itypj)
+              e1=fac*fac*aa
+              e2=fac*bb
               evdwij=e1+e2
               evdw=evdw+(1.0d0-sss)*evdwij
 C 
@@ -193,8 +193,8 @@ C Change 12/1/95 to calculate four-body interactions
               rrij=1.0D0/rij
               eps0ij=eps(itypi,itypj)
               fac=rrij**expon2
-              e1=fac*fac*aa(itypi,itypj)
-              e2=fac*bb(itypi,itypj)
+              e1=fac*fac*aa
+              e2=fac*bb
               evdwij=e1+e2
               evdw=evdw+sss*evdwij
 C 
@@ -277,8 +277,8 @@ C
             if (sss.lt.1.0d0) then
               r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
               fac=r_shift_inv**expon
-              e1=fac*fac*aa(itypi,itypj)
-              e2=fac*bb(itypi,itypj)
+              e1=fac*fac*aa
+              e2=fac*bb
               evdwij=e_augm+e1+e2
 cd            sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
 cd            epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
@@ -360,8 +360,8 @@ C
             if (sss.gt.0.0d0) then
               r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
               fac=r_shift_inv**expon
-              e1=fac*fac*aa(itypi,itypj)
-              e2=fac*bb(itypi,itypj)
+              e1=fac*fac*aa
+              e2=fac*bb
               evdwij=e_augm+e1+e2
 cd            sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
 cd            epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
@@ -474,16 +474,16 @@ C Calculate the angle-dependent terms of energy & contributions to derivatives.
 C Calculate whole angle-dependent part of epsilon and contributions
 C to its derivatives
               fac=(rrij*sigsq)**expon2
-              e1=fac*fac*aa(itypi,itypj)
-              e2=fac*bb(itypi,itypj)
+              e1=fac*fac*aa
+              e2=fac*bb
               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
               eps2der=evdwij*eps3rt
               eps3der=evdwij*eps2rt
               evdwij=evdwij*eps2rt*eps3rt
               evdw=evdw+evdwij*(1.0d0-sss)
               if (lprn) then
-              sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
-              epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+              sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+              epsi=bb**2/aa
 cd              write (iout,'(2(a3,i3,2x),15(0pf7.3))')
 cd     &          restyp(itypi),i,restyp(itypj),j,
 cd     &          epsi,sigm,chi1,chi2,chip1,chip2,
@@ -587,16 +587,16 @@ C Calculate the angle-dependent terms of energy & contributions to derivatives.
 C Calculate whole angle-dependent part of epsilon and contributions
 C to its derivatives
               fac=(rrij*sigsq)**expon2
-              e1=fac*fac*aa(itypi,itypj)
-              e2=fac*bb(itypi,itypj)
+              e1=fac*fac*aa
+              e2=fac*bb
               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
               eps2der=evdwij*eps3rt
               eps3der=evdwij*eps2rt
               evdwij=evdwij*eps2rt*eps3rt
               evdw=evdw+evdwij*sss
               if (lprn) then
-              sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
-              epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+              sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+              epsi=bb**2/aa
 cd              write (iout,'(2(a3,i3,2x),15(0pf7.3))')
 cd     &          restyp(itypi),i,restyp(itypj),j,
 cd     &          epsi,sigm,chi1,chi2,chip1,chip2,
@@ -701,6 +701,33 @@ c            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
           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-
+     &        ((positi-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+         sslipj=sscalelip(fracinbuf)
+         ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
+        elseif (zi.gt.bufliptop) then
+         fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
+         sslipj=sscalelip(fracinbuf)
+         ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
+        else
+         sslipj=1.0d0
+         ssgradlipj=0.0
+        endif
+       else
+         sslipj=0.0d0
+         ssgradlipj=0.0
+       endif
+      aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+     &  +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
+      bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+     &  +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
+
       dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
       xj_safe=xj
       yj_safe=yj
@@ -761,8 +788,8 @@ cd     &          rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
 c---------------------------------------------------------------
               rij_shift=1.0D0/rij_shift 
               fac=rij_shift**expon
-              e1=fac*fac*aa(itypi,itypj)
-              e2=fac*bb(itypi,itypj)
+              e1=fac*fac*aa
+              e2=fac*bb
               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
               eps2der=evdwij*eps3rt
               eps3der=evdwij*eps2rt
@@ -771,8 +798,8 @@ c     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
               evdwij=evdwij*eps2rt*eps3rt
               evdw=evdw+evdwij*(1.0d0-sss)
               if (lprn) then
-              sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
-              epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+              sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+              epsi=bb**2/aa
               write (iout,'(2(a3,i3,2x),17(0pf7.3))')
      &          restyp(itypi),i,restyp(itypj),j,
      &          epsi,sigm,chi1,chi2,chip1,chip2,
@@ -795,6 +822,8 @@ C Calculate the radial part of the gradient
               gg(1)=xj*fac
               gg(2)=yj*fac
               gg(3)=zj*fac
+              gg_lipi(3)=ssgradlipi*evdwij
+              gg_lipj(3)=ssgradlipj*evdwij
 C Calculate angular part of the gradient.
               call sc_grad_scale(1.0d0-sss)
             endif
@@ -883,6 +912,32 @@ c            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
           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-
+     &        ((positi-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+         sslipj=sscalelip(fracinbuf)
+         ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
+        elseif (zi.gt.bufliptop) then
+         fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
+         sslipj=sscalelip(fracinbuf)
+         ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
+        else
+         sslipj=1.0d0
+         ssgradlipj=0.0
+        endif
+       else
+         sslipj=0.0d0
+         ssgradlipj=0.0
+       endif
+      aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+     &  +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
+      bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+     &  +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
       dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
       xj_safe=xj
       yj_safe=yj
@@ -943,8 +998,8 @@ cd     &          rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
 c---------------------------------------------------------------
               rij_shift=1.0D0/rij_shift 
               fac=rij_shift**expon
-              e1=fac*fac*aa(itypi,itypj)
-              e2=fac*bb(itypi,itypj)
+              e1=fac*fac*aa
+              e2=fac*bb
               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
               eps2der=evdwij*eps3rt
               eps3der=evdwij*eps2rt
@@ -953,8 +1008,8 @@ c     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
               evdwij=evdwij*eps2rt*eps3rt
               evdw=evdw+evdwij*sss
               if (lprn) then
-              sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
-              epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+              sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+              epsi=bb**2/aa
               write (iout,'(2(a3,i3,2x),17(0pf7.3))')
      &          restyp(itypi),i,restyp(itypj),j,
      &          epsi,sigm,chi1,chi2,chip1,chip2,
@@ -977,6 +1032,8 @@ C Calculate the radial part of the gradient
               gg(1)=xj*fac
               gg(2)=yj*fac
               gg(3)=zj*fac
+              gg_lipi(3)=ssgradlipi*evdwij
+              gg_lipj(3)=ssgradlipj*evdwij
 C Calculate angular part of the gradient.
               call sc_grad_scale(sss)
             endif
@@ -1073,8 +1130,8 @@ C I hate to put IF's in the loops, but here don't have another choice!!!!
 c---------------------------------------------------------------
               rij_shift=1.0D0/rij_shift 
               fac=rij_shift**expon
-              e1=fac*fac*aa(itypi,itypj)
-              e2=fac*bb(itypi,itypj)
+              e1=fac*fac*aa
+              e2=fac*bb
               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
               eps2der=evdwij*eps3rt
               eps3der=evdwij*eps2rt
@@ -1083,8 +1140,8 @@ c---------------------------------------------------------------
               evdwij=evdwij*eps2rt*eps3rt
               evdw=evdw+(evdwij+e_augm)*(1.0d0-sss)
               if (lprn) then
-              sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
-              epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+              sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+              epsi=bb**2/aa
               write (iout,'(2(a3,i3,2x),17(0pf7.3))')
      &          restyp(itypi),i,restyp(itypj),j,
      &          epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
@@ -1195,8 +1252,8 @@ C I hate to put IF's in the loops, but here don't have another choice!!!!
 c---------------------------------------------------------------
               rij_shift=1.0D0/rij_shift 
               fac=rij_shift**expon
-              e1=fac*fac*aa(itypi,itypj)
-              e2=fac*bb(itypi,itypj)
+              e1=fac*fac*aa
+              e2=fac*bb
               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
               eps2der=evdwij*eps3rt
               eps3der=evdwij*eps2rt
@@ -1205,8 +1262,8 @@ c---------------------------------------------------------------
               evdwij=evdwij*eps2rt*eps3rt
               evdw=evdw+(evdwij+e_augm)*sss
               if (lprn) then
-              sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
-              epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+              sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+              epsi=bb**2/aa
               write (iout,'(2(a3,i3,2x),17(0pf7.3))')
      &          restyp(itypi),i,restyp(itypj),j,
      &          epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
@@ -1263,10 +1320,10 @@ c      write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12
       enddo 
 c      write (iout,*) "gg",(gg(k),k=1,3)
       do k=1,3
-        gvdwx(k,i)=gvdwx(k,i)-gg(k)
+        gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
      &            +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
      &          +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*scalfac
-        gvdwx(k,j)=gvdwx(k,j)+gg(k)
+        gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)
      &            +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
      &          +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*scalfac
 c        write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
@@ -1278,8 +1335,8 @@ C
 C Calculate the components of the gradient in DC and X
 C
       do l=1,3
-        gvdwc(l,i)=gvdwc(l,i)-gg(l)
-        gvdwc(l,j)=gvdwc(l,j)+gg(l)
+        gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(k)
+        gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(k)
       enddo
       return
       end
UNRESPACK 3.x