projects / unres.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
dodanie parametrow do nanotube
[unres.git] / source / wham / src-M / energy_p_new.F
diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F
index a57f604..b5f5705 100644 (file)
--- a/source/wham/src-M/energy_p_new.F
+++ b/source/wham/src-M/energy_p_new.F
@@ -250,11 +250,28 @@ C
      &                wturn6*fact(5)*gcorr6_turn(j,i)+
      &                wsccor*fact(2)*gsccorc(j,i)
      &               +wliptran*gliptranc(j,i)
+     &                 +welec*gshieldc(j,i)
+     &                 +welec*gshieldc_loc(j,i)
+     &                 +wcorr*gshieldc_ec(j,i)
+     &                 +wcorr*gshieldc_loc_ec(j,i)
+     &                 +wturn3*gshieldc_t3(j,i)
+     &                 +wturn3*gshieldc_loc_t3(j,i)
+     &                 +wturn4*gshieldc_t4(j,i)
+     &                 +wturn4*gshieldc_loc_t4(j,i)
+     &                 +wel_loc*gshieldc_ll(j,i)
+     &                 +wel_loc*gshieldc_loc_ll(j,i)
+
           gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+
      &                  wbond*gradbx(j,i)+
      &                  wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+
      &                  wsccor*fact(2)*gsccorx(j,i)
      &                 +wliptran*gliptranx(j,i)
+     &                 +welec*gshieldx(j,i)
+     &                 +wcorr*gshieldx_ec(j,i)
+     &                 +wturn3*gshieldx_t3(j,i)
+     &                 +wturn4*gshieldx_t4(j,i)
+     &                 +wel_loc*gshieldx_ll(j,i)
+
         else
           gradc(j,i,icg)=fact(1)*wsc*gvdwc(j,i)
      &                +fact(1)*wscp*gvdwc_scp(j,i)+
@@ -312,11 +329,28 @@ C
      &                wturn6*fact(5)*gcorr6_turn(j,i)+
      &                wsccor*fact(2)*gsccorc(j,i)
      &               +wliptran*gliptranc(j,i)
+     &                 +welec*gshieldc(j,i)
+     &                 +welec*gshieldc_loc(j,i)
+     &                 +wcorr*gshieldc_ec(j,i)
+     &                 +wcorr*gshieldc_loc_ec(j,i)
+     &                 +wturn3*gshieldc_t3(j,i)
+     &                 +wturn3*gshieldc_loc_t3(j,i)
+     &                 +wturn4*gshieldc_t4(j,i)
+     &                 +wturn4*gshieldc_loc_t4(j,i)
+     &                 +wel_loc*gshieldc_ll(j,i)
+     &                 +wel_loc*gshieldc_loc_ll(j,i)
+
           gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+
      &                  wbond*gradbx(j,i)+
      &                  wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+
      &                  wsccor*fact(1)*gsccorx(j,i)
      &                 +wliptran*gliptranx(j,i)
+     &                 +welec*gshieldx(j,i)
+     &                 +wcorr*gshieldx_ec(j,i)
+     &                 +wturn3*gshieldx_t3(j,i)
+     &                 +wturn4*gshieldx_t4(j,i)
+     &                 +wel_loc*gshieldx_ll(j,i)
+
               else
           gradc(j,i,icg)=fact(1)*wsc*gvdwc(j,i)+
      &                   fact(1)*wscp*gvdwc_scp(j,i)+
@@ -2171,6 +2205,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)
@@ -2216,6 +2275,28 @@ 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
@@ -2722,6 +2803,7 @@ 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)
 c          write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
 C          write (iout,'(a6,2i5,0pf7.3)')
 C     &            'eelloc',i,j,eel_loc_ij
@@ -2779,11 +2861,13 @@ C Partial derivatives in virtual-bond dihedral angles gamma
      &            (a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
      &           +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j))
      &    *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
 
           gel_loc_loc(j-1)=gel_loc_loc(j-1)+ 
      &            (a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
      &           +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j))
      &    *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
 
 cd          call checkint3(i,j,mu1,mu2,a22,a23,a32,a33,acipa,eel_loc_ij)
 cd          write(iout,*) 'agg  ',agg
@@ -2797,6 +2881,7 @@ C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
             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
@@ -2809,18 +2894,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
@@ -3083,6 +3172,37 @@ C Third- and fourth-order contributions from turns
       double precision agg(3,4),aggi(3,4),aggi1(3,4),
      &    aggj(3,4),aggj1(3,4),a_temp(2,2)
       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,j1,j2
+          zj=(c(3,j)+c(3,j+1))/2.0d0
+C          xj=mod(xj,boxxsize)
+C          if (xj.lt.0) xj=xj+boxxsize
+C          yj=mod(yj,boxysize)
+C          if (yj.lt.0) yj=yj+boxysize
+          zj=mod(zj,boxzsize)
+          if (zj.lt.0) zj=zj+boxzsize
+C          if ((zj.lt.0).or.(xj.lt.0).or.(yj.lt.0)) write (*,*) "CHUJ"
+       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
+
       if (j.eq.i+2) then
       if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
 C changes suggested by Ana to avoid out of bounds
@@ -3120,8 +3240,11 @@ C        fac_shield(j)=0.6
 
         eello_turn3=eello_turn3+0.5d0*(pizda(1,1)+pizda(2,2))
      &  *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
         eello_t3=0.5d0*(pizda(1,1)+pizda(2,2))
      &  *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
 
 cd        write (2,*) 'i,',i,' j',j,'eello_turn3',
 cd     &    0.5d0*(pizda(1,1)+pizda(2,2)),
@@ -3176,6 +3299,8 @@ C Derivatives in gamma(i)
         call matmat2(a_temp(1,1),pizda(1,1),pizda(1,1))
         gel_loc_turn3(i)=gel_loc_turn3(i)+0.5d0*(pizda(1,1)+pizda(2,2))
      &   *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
 C Derivatives in gamma(i+1)
         call matmat2(EUg(1,1,i+1),EUgder(1,1,i+2),auxmat2(1,1))
         call transpose2(auxmat2(1,1),pizda(1,1))
@@ -3183,6 +3308,7 @@ C Derivatives in gamma(i+1)
         gel_loc_turn3(i+1)=gel_loc_turn3(i+1)
      &    +0.5d0*(pizda(1,1)+pizda(2,2))
      &   *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
 
 C Cartesian derivatives
         do l=1,3
@@ -3194,6 +3320,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+1.0d0)
 
           a_temp(1,1)=aggi1(l,1)
           a_temp(1,2)=aggi1(l,2)
@@ -3203,6 +3330,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+1.0d0)
 
           a_temp(1,1)=aggj(l,1)
           a_temp(1,2)=aggj(l,2)
@@ -3212,6 +3340,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+1.0d0)
 
           a_temp(1,1)=aggj1(l,1)
           a_temp(1,2)=aggj1(l,2)
@@ -3221,6 +3350,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+1.0d0)
 
         enddo
         endif
@@ -3330,6 +3460,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+1.0d0)
 
 C Derivatives in gamma(i+1)
         call transpose2(EUgder(1,1,i+2),e2tder(1,1))
@@ -3340,6 +3471,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+1.0d0)
 
 C Derivatives in gamma(i+2)
         call transpose2(EUgder(1,1,i+3),e3tder(1,1))
@@ -3353,6 +3485,7 @@ 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+1.0d0)
 
 C Cartesian derivatives
 
@@ -3375,6 +3508,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+1.0d0)
 
           enddo
         endif
@@ -3395,6 +3529,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+1.0d0)
 
           a_temp(1,1)=aggi1(l,1)
           a_temp(1,2)=aggi1(l,2)
@@ -3411,6 +3546,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+1.0d0)
 
           a_temp(1,1)=aggj(l,1)
           a_temp(1,2)=aggj(l,2)
@@ -3427,6 +3563,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+1.0d0)
 
           a_temp(1,1)=aggj1(l,1)
           a_temp(1,2)=aggj1(l,2)
@@ -3443,8 +3580,17 @@ 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+1.0d0)
 
         enddo
+         gshieldc_t4(3,i)=gshieldc_t4(3,i)+
+     &     ssgradlipi*eello_t4/4.0d0*lipscale
+         gshieldc_t4(3,j)=gshieldc_t4(3,j)+
+     &     ssgradlipj*eello_t4/4.0d0*lipscale
+         gshieldc_t4(3,i-1)=gshieldc_t4(3,i-1)+
+     &     ssgradlipi*eello_t4/4.0d0*lipscale
+         gshieldc_t4(3,j-1)=gshieldc_t4(3,j-1)+
+     &     ssgradlipj*eello_t4/4.0d0*lipscale
         endif
  178  continue
       endif          
UNRESPACK 3.x