small chanegs in nanotube + working wham for lipid

diff --git a/source/unres/src_MD-M/COMMON.TORSION b/source/unres/src_MD-M/COMMON.TORSION
index e7c54a0..20cf153 100644 (file)
--- a/source/unres/src_MD-M/COMMON.TORSION
+++ b/source/unres/src_MD-M/COMMON.TORSION
@@ -40,14 +40,15 @@ C           surface
       double precision b1,b2,cc,dd,ee,ctilde,dtilde,b2tilde,b1tilde,
      & b,bnew1,bnew2,eenew,gtb1,gtb2,eeold,gtee
       integer nloctyp,iloctyp(-ntyp1:ntyp1),itype2loc(-ntyp1:ntyp1)
-      common/fourier/ b1(2,maxres),b2(2,maxres),b(13,-ntyp:ntyp),
+      common/fourier/ b1(2,-maxres:maxres),b2(2,-maxres:maxres),
+     &     b(13,-ntyp:ntyp),
      &    bnew1(3,2,-ntyp:ntyp),bnew2(3,2,-ntyp:ntyp),
      &    cc(2,2,-ntyp:ntyp),
      &    dd(2,2,-ntyp:ntyp),eeold(2,2,-ntyp:ntyp),
      &    eenew(2,-ntyp:ntyp),
      &    ee(2,2,maxres),
      &    ctilde(2,2,-ntyp:ntyp),
-     &    dtilde(2,2,-ntyp:ntyp),b1tilde(2,maxres),
-     &    b2tilde(2,maxres),
+     &    dtilde(2,2,-ntyp:ntyp),b1tilde(2,-maxres:maxres),
+     &    b2tilde(2,-maxres:maxres),
      &    gtb1(2,maxres),gtb2(2,maxres),gtEE(2,2,maxres),
      &    nloctyp,iloctyp,itype2loc

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 1390fe3..e3a17dd 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new_barrier.F
@@ -4647,8 +4647,8 @@ c     &   a33*gmuji2(4)
 #endif
 cd          write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
 
-          if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
-     &            'eelloc',i,j,eel_loc_ij
+          if (energy_dec) write (iout,'(a6,2i5,0pf7.3,2f7.3)')
+     &            'eelloc',i,j,eel_loc_ij,a22*muij(1),a23*muij(2)
 c           if (eel_loc_ij.ne.0)
 c     &      write (iout,'(a4,2i4,8f9.5)')'chuj',
 c     &     i,j,a22,muij(1),a23,muij(2),a32,muij(3),a33,muij(4)
@@ -12611,17 +12611,17 @@ C now calculate distance from center of tube and direction vectors
       zmin=boxzsize
 
         do j=-1,1
-         vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize)
+         vectube(1)=dmod((c(1,i)+c(1,i+1))/2.0d0,boxxsize)
          vectube(1)=vectube(1)+boxxsize*j
-         vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxysize)
+         vectube(2)=dmod((c(2,i)+c(2,i+1))/2.0d0,boxysize)
          vectube(2)=vectube(2)+boxysize*j
-         vectube(3)=mod((c(3,i)+c(3,i+1))/2.0d0,boxzsize)
+         vectube(3)=dmod((c(3,i)+c(3,i+1))/2.0d0,boxzsize)
          vectube(3)=vectube(3)+boxzsize*j
 
 
-         xminact=abs(vectube(1)-tubecenter(1))
-         yminact=abs(vectube(2)-tubecenter(2))
-         zminact=abs(vectube(3)-tubecenter(3))
+         xminact=dabs(vectube(1)-tubecenter(1))
+         yminact=dabs(vectube(2)-tubecenter(2))
+         zminact=dabs(vectube(3)-tubecenter(3))
 
            if (xmin.gt.xminact) then
             xmin=xminact
@@ -12674,13 +12674,13 @@ C go to 667
          enecavtube(i)=0.0
          faccav=0.0
          else
-         denominator=(1.0+dcavtubpep*rdiff6*rdiff6)
+         denominator=(1.0d0+dcavtubpep*rdiff6*rdiff6)
          enecavtube(i)=
-     &   (bcavtubpep*rdiff+acavtubpep*sqrt(rdiff)+ccavtubpep)
+     &   (bcavtubpep*rdiff+acavtubpep*dsqrt(rdiff)+ccavtubpep)
      &   /denominator
          enecavtube(i)=0.0
-         faccav=((bcavtubpep*1.0d0+acavtubpep/2.0d0/sqrt(rdiff))
-     &   *denominator-(bcavtubpep*rdiff+acavtubpep*sqrt(rdiff)
+         faccav=((bcavtubpep*1.0d0+acavtubpep/2.0d0/dsqrt(rdiff))
+     &   *denominator-(bcavtubpep*rdiff+acavtubpep*dsqrt(rdiff)
      &   +ccavtubpep)*rdiff6**2.0d0/rdiff*dcavtubpep*12.0d0)
      &   /denominator**2.0d0
 C         faccav=0.0
@@ -12701,7 +12701,7 @@ C now direction of gg_tube vector
         enddo
 
        do i=itube_start,itube_end
-        enecavtube(i)=0.0 
+        enecavtube(i)=0.0d0
 C Lets not jump over memory as we use many times iti
          iti=itype(i)
 C lets ommit dummy atoms for now
@@ -12713,17 +12713,17 @@ C      .or.(iti.eq.10)
       ymin=boxysize
       zmin=boxzsize
         do j=-1,1
-         vectube(1)=mod((c(1,i+nres)),boxxsize)
+         vectube(1)=dmod((c(1,i+nres)),boxxsize)
          vectube(1)=vectube(1)+boxxsize*j
-         vectube(2)=mod((c(2,i+nres)),boxysize)
+         vectube(2)=dmod((c(2,i+nres)),boxysize)
          vectube(2)=vectube(2)+boxysize*j
-         vectube(3)=mod((c(3,i+nres)),boxzsize)
+         vectube(3)=dmod((c(3,i+nres)),boxzsize)
          vectube(3)=vectube(3)+boxzsize*j
 
 
-         xminact=abs(vectube(1)-tubecenter(1))
-         yminact=abs(vectube(2)-tubecenter(2))
-         zminact=abs(vectube(3)-tubecenter(3))
+         xminact=dabs(vectube(1)-tubecenter(1))
+         yminact=dabs(vectube(2)-tubecenter(2))
+         zminact=dabs(vectube(3)-tubecenter(3))
 
            if (xmin.gt.xminact) then
             xmin=xminact
@@ -12772,16 +12772,16 @@ 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+nres)=0.0
-         faccav=0.0
+         enecavtube(i+nres)=0.0d0
+         faccav=0.0d0
          else
-         denominator=(1.0+dcavtub(iti)*rdiff6*rdiff6)
+         denominator=(1.0d0+dcavtub(iti)*rdiff6*rdiff6)
          enecavtube(i+nres)=
-     &   (bcavtub(iti)*rdiff+acavtub(iti)*sqrt(rdiff)+ccavtub(iti))
+     &   (bcavtub(iti)*rdiff+acavtub(iti)*dsqrt(rdiff)+ccavtub(iti))
      &   /denominator
 C         enecavtube(i)=0.0
-         faccav=((bcavtub(iti)*1.0d0+acavtub(iti)/2.0d0/sqrt(rdiff))
-     &   *denominator-(bcavtub(iti)*rdiff+acavtub(iti)*sqrt(rdiff)
+         faccav=((bcavtub(iti)*1.0d0+acavtub(iti)/2.0d0/dsqrt(rdiff))
+     &   *denominator-(bcavtub(iti)*rdiff+acavtub(iti)*dsqrt(rdiff)
      &   +ccavtub(iti))*rdiff6**2.0d0/rdiff*dcavtub(iti)*12.0d0)
      &   /denominator**2.0d0
 C         faccav=0.0

diff --git a/source/wham/src-M/COMMON.SHIELD b/source/wham/src-M/COMMON.SHIELD
index 1f96c94..ac7ac98 100644 (file)
--- a/source/wham/src-M/COMMON.SHIELD
+++ b/source/wham/src-M/COMMON.SHIELD
@@ -1,9 +1,11 @@
        double precision VSolvSphere,VSolvSphere_div,long_r_sidechain,
      & short_r_sidechain,fac_shield,grad_shield_side,grad_shield,
-     & buff_shield,wshield,grad_shield_loc            
+     & buff_shield,wshield,grad_shield_loc,lipscale,sslipi,sslipj,
+     & ssgradlipi,ssgradlipj
        integer  ishield_list,shield_list,ees0plist
        common /shield/ VSolvSphere,VSolvSphere_div,buff_shield,
-     & long_r_sidechain(ntyp),
+     & long_r_sidechain(ntyp),sslipi,sslipj,
+     & ssgradlipi,ssgradlipj,lipscale,
      & short_r_sidechain(ntyp),fac_shield(maxres),wshield,
      & grad_shield_side(3,maxcont,-1:maxres),grad_shield(3,-1:maxres),
      &  grad_shield_loc(3,maxcont,-1:maxres),

diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F
index 378eaae..a36e177 100644 (file)
--- a/source/wham/src-M/energy_p_new.F
+++ b/source/wham/src-M/energy_p_new.F
@@ -2398,8 +2398,12 @@ C          fac_shield(j)=0.6
           fac_shield(j)=1.0
           eesij=(el1+el2)
           ees=ees+eesij
+     &*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+
           endif
           evdw1=evdw1+evdwij*sss
+     &*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+
 c             write (iout,'(a6,2i5,0pf7.3,2i5,2e11.3)') 
 c     &'evdw1',i,j,evdwij
 c     &,iteli,itelj,aaa,evdw1
@@ -2414,7 +2418,11 @@ C Calculate contributions to the Cartesian gradient.
 C
 #ifdef SPLITELE
           facvdw=-6*rrmij*(ev1+evdwij)*sss
+     &    *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+
           facel=-3*rrmij*(el1+eesij)
+     &*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+
           fac1=fac
           erij(1)=xj*rmij
           erij(2)=yj*rmij
@@ -2498,8 +2506,13 @@ C          ggg(2)=facvdw*yj
 C          ggg(3)=facvdw*zj
           if (sss.gt.0.0) then
           ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj
+     &    *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
           ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj
+     &    *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+
           ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj
+     &    *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+
           else
           ggg(1)=0.0
           ggg(2)=0.0
@@ -2510,6 +2523,11 @@ C          ggg(3)=facvdw*zj
             gvdwpp(k,i)=gvdwpp(k,i)+ghalf
             gvdwpp(k,j)=gvdwpp(k,j)+ghalf
           enddo
+           gvdwpp(3,j)=gvdwpp(3,j)+
+     &     sss*ssgradlipj*evdwij/2.0d0*lipscale**2
+           gvdwpp(3,i)=gvdwpp(3,i)+
+     &     sss*ssgradlipi*evdwij/2.0d0*lipscale**2
+
 *
 * Loop over residues i+1 thru j-1.
 *
@@ -2520,6 +2538,8 @@ C          ggg(3)=facvdw*zj
           enddo
 #else
           facvdw=(ev1+evdwij)*sss
+     &    *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+
           facel=el1+eesij  
           fac1=fac
           fac=-3*rrmij*(facvdw+facvdw+facel)
@@ -2564,6 +2584,8 @@ cd   &          (dcosg(k),k=1,3)
           do k=1,3
             ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k) 
      &      *fac_shield(i)**2*fac_shield(j)**2
+     &      *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+
           enddo
           do k=1,3
             ghalf=0.5D0*ggg(k)
@@ -2571,11 +2593,15 @@ cd   &          (dcosg(k),k=1,3)
      &               +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
      &               + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
      &           *fac_shield(i)**2*fac_shield(j)**2
+     &      *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+
 
             gelc(k,j)=gelc(k,j)+ghalf
      &               +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
      &               + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
      &           *fac_shield(i)**2*fac_shield(j)**2
+     &      *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+
           enddo
           do k=i+1,j-1
             do l=1,3
@@ -2831,9 +2857,11 @@ C           fac_shield(j)=0.6
           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
+C          write (iout,'(a3,i4,a3,i4,a8,4f8.3)') 
+C     & 'i',i,' j',j,' eel_loc_ij',eel_loc_ij,sslipi,
+C     & sslipj,lipscale
+C          write (iout,'(a6,2i5,0pf7.3,2f7.3)')
+C     &            'eelloc',i,j,eel_loc_ij,a22*muij(1),a23*muij(2)
 C          write(iout,*) 'muije=',i,j,muij(1),muij(2),muij(3),muij(4)
 c          write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3)
 C          eel_loc=eel_loc+eel_loc_ij
@@ -3272,7 +3300,9 @@ C        fac_shield(j)=0.6
         eello_t3=0.5d0*(pizda(1,1)+pizda(2,2))
      &  *fac_shield(i)*fac_shield(j)
      &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
-
+          write (iout,'(a3,i4,a3,i4,a8,4f8.3)') 
+     & 'i',i,' j',j,' eturn3',eello_t3,sslipi,
+     & sslipj,lipscale
 cd        write (2,*) 'i,',i,' j',j,'eello_turn3',
 cd     &    0.5d0*(pizda(1,1)+pizda(2,2)),
 cd     &    ' eello_turn3_num',4*eello_turn3_num
@@ -3431,6 +3461,8 @@ C        fac_shield(j)=0.6
 
         eello_turn4=eello_turn4-(s1+s2+s3)
      &  *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
         eello_t4=-(s1+s2+s3)
      &  *fac_shield(i)*fac_shield(j)
 

diff --git a/source/wham/src-M/parmread.F b/source/wham/src-M/parmread.F
index 2fc37d3..ce1bf31 100644 (file)
--- a/source/wham/src-M/parmread.F
+++ b/source/wham/src-M/parmread.F
@@ -71,6 +71,7 @@ c If reading not own parameters, skip assignment
       call reada(controlcard,"ATRISS",atriss,0.3D0)
       call reada(controlcard,"BTRISS",btriss,0.02D0)
       call reada(controlcard,"CTRISS",ctriss,1.0D0)
+      call reada(controlcard,"LIPSCALE",lipscale,1.3D0)
       dyn_ss=(index(controlcard,'DYN_SS').gt.0)
       write(iout,*) "ATRISS",atriss
       write(iout,*) "BTRISS",btriss