small fix for Gly as 1st residue and in dynamic disulfides

diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F
index cba6b5e..be6ed2a 100644 (file)
--- a/source/wham/src-M/energy_p_new.F
+++ b/source/wham/src-M/energy_p_new.F
@@ -23,7 +23,7 @@ cMS$ATTRIBUTES C ::  proc_proc
       include 'COMMON.SBRIDGE'
       include 'COMMON.CHAIN'
       double precision fact(6)
-cd      write(iout, '(a,i2)')'Calling etotal ipot=',ipot
+Cd      write(iout, '(a,i2)')'Calling etotal ipot=',ipot
 cd    print *,'nnt=',nnt,' nct=',nct
 C
 C Compute the side-chain and electrostatic interaction energy
@@ -825,6 +825,7 @@ c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
 c      if (icall.gt.0) lprn=.true.
       ind=0
       do i=iatsc_s,iatsc_e
+C        write(iout,*) i,"i",iatsc_s,iatsc_e
         itypi=iabs(itype(i))
         if (itypi.eq.ntyp1) cycle
         itypi1=iabs(itype(i+1))
@@ -891,6 +892,8 @@ C     &                        'evdw',i,j,evdwij,'tss',evdw,evdw_t
              enddo! k
             ELSE
             ind=ind+1
+C            write(iout,*) j,"j",istart(i,iint),iend(i,iint)
+
             itypj=iabs(itype(j))
             if (itypj.eq.ntyp1) cycle
             dscj_inv=vbld_inv(j+nres)
@@ -1955,7 +1958,7 @@ C 13-go grudnia roku pamietnego...
       double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
      &                   0.0d0,1.0d0,0.0d0,
      &                   0.0d0,0.0d0,1.0d0/
-cd      write(iout,*) 'In EELEC'
+      write(iout,*) 'In EELEC'
 cd      do i=1,nloctyp
 cd        write(iout,*) 'Type',i
 cd        write(iout,*) 'B1',B1(:,i)
@@ -2013,15 +2016,16 @@ C      write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
         gcorr_loc(i)=0.0d0
       enddo
       do i=iatel_s,iatel_e
-          if (i.eq.1) then 
-           if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1
-     &  .or. itype(i+2).eq.ntyp1) cycle
-          else
+          write (iout,*) i,"i2",itype(i)
+          if (i.eq.1) cycle 
+C           if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1
+C     &  .or. itype(i+2).eq.ntyp1) cycle
+C          else
         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
      &  .or. itype(i+2).eq.ntyp1
      &  .or. itype(i-1).eq.ntyp1
      &) cycle
-         endif
+C         endif
         if (itel(i).eq.0) goto 1215
         dxi=dc(1,i)
         dyi=dc(2,i)
@@ -2041,6 +2045,8 @@ C      write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
         num_conti=0
 C        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
         do j=ielstart(i),ielend(i)
+C           write(iout,*) j,"j2"
+
           if (j.eq.1) then
            if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
      & .or.itype(j+2).eq.ntyp1
@@ -2051,6 +2057,7 @@ C        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
      & .or.itype(j-1).eq.ntyp1
      &) cycle
          endif
+C           write(iout,*) j,"j2"
 C
 C) cycle
           if (itel(j).eq.0) goto 1216
@@ -2844,7 +2851,7 @@ C Cartesian derivatives
      &      +0.5d0*(pizda(1,1)+pizda(2,2))
         enddo
         endif
-      else if (j.eq.i+3 .and. itype(i+2).ne.ntyp1) then
+      else if (j.eq.i+3 .and. itype(i+2).ne.ntyp1.and.(i.gt.1)) then
       if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
 C changes suggested by Ana to avoid out of bounds
      & .or.((i+5).gt.nres)
@@ -3682,7 +3689,7 @@ c     &    rad2deg*phii,rad2deg*phii1,ethetai
 c 1215   continue
       enddo
       ethetacnstr=0.0d0
-C      print *,ithetaconstr_start,ithetaconstr_end,"TU"
+      print *,ntheta_constr,"TU"
       do i=1,ntheta_constr
         itheta=itheta_constr(i)
         thetiii=theta(itheta)
@@ -3706,7 +3713,7 @@ C     &    i,itheta,rad2deg*thetiii,
 C     &    rad2deg*theta_constr0(i),  rad2deg*theta_drange(i),
 C     &    rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4,
 C     &    gloc(itheta+nphi-2,icg)
-        endif
+C        endif
       enddo
 C Ufff.... We've done all this!!! 
       return

diff --git a/source/wham/src-M/ssMD.F b/source/wham/src-M/ssMD.F
index 5080b18..bbf3206 100644 (file)
--- a/source/wham/src-M/ssMD.F
+++ b/source/wham/src-M/ssMD.F
@@ -150,11 +150,113 @@ c-------END TESTING CODE
       dyi=dc_norm(2,nres+i)
       dzi=dc_norm(3,nres+i)
       dsci_inv=vbld_inv(i+nres)
+        xi=c(1,nres+i)
+        yi=c(2,nres+i)
+        zi=c(3,nres+i)
+C returning the ith atom to box
+          xi=mod(xi,boxxsize)
+          if (xi.lt.0) xi=xi+boxxsize
+          yi=mod(yi,boxysize)
+          if (yi.lt.0) yi=yi+boxysize
+          zi=mod(zi,boxzsize)
+          if (zi.lt.0) zi=zi+boxzsize
+       if ((zi.gt.bordlipbot)
+     &.and.(zi.lt.bordliptop)) then
+C the energy transfer exist
+        if (zi.lt.buflipbot) then
+C what fraction I am in
+         fracinbuf=1.0d0-
+     &        ((zi-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+         sslipi=sscalelip(fracinbuf)
+         ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
+        elseif (zi.gt.bufliptop) then
+         fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
+         sslipi=sscalelip(fracinbuf)
+         ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
+        else
+         sslipi=1.0d0
+         ssgradlipi=0.0
+        endif
+       else
+         sslipi=0.0d0
+         ssgradlipi=0.0
+       endif
+            xj=c(1,nres+j)
+            yj=c(2,nres+j)
+            zj=c(3,nres+j)
+C returning jth atom to box
+          xj=mod(xj,boxxsize)
+          if (xj.lt.0) xj=xj+boxxsize
+          yj=mod(yj,boxysize)
+          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
+      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
 
       itypj=itype(j)
-      xj=c(1,nres+j)-c(1,nres+i)
-      yj=c(2,nres+j)-c(2,nres+i)
-      zj=c(3,nres+j)-c(3,nres+i)
+C checking the distance
+      dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+      xj_safe=xj
+      yj_safe=yj
+      zj_safe=zj
+      subchap=0
+C finding the closest
+      do xshift=-1,1
+      do yshift=-1,1
+      do zshift=-1,1
+          xj=xj_safe+xshift*boxxsize
+          yj=yj_safe+yshift*boxysize
+          zj=zj_safe+zshift*boxzsize
+          dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+          if(dist_temp.lt.dist_init) then
+            dist_init=dist_temp
+            xj_temp=xj
+            yj_temp=yj
+            zj_temp=zj
+            subchap=1
+          endif
+       enddo
+       enddo
+       enddo
+       if (subchap.eq.1) then
+          xj=xj_temp-xi
+          yj=yj_temp-yi
+          zj=zj_temp-zi
+       else
+          xj=xj_safe-xi
+          yj=yj_safe-yi
+          zj=zj_safe-zi
+       endif
+
+C      xj=c(1,nres+j)-c(1,nres+i)
+C      yj=c(2,nres+j)-c(2,nres+i)
+C      zj=c(3,nres+j)-c(3,nres+i)
       dxj=dc_norm(1,nres+j)
       dyj=dc_norm(2,nres+j)
       dzj=dc_norm(3,nres+j)
@@ -187,9 +289,9 @@ c      om12=dxi*dxj+dyi*dyj+dzi*dzj
 
       ljXs=sig-sig0ij
       ljA=eps1*eps2rt**2*eps3rt**2
-      ljB=ljA*bb(itypi,itypj)
-      ljA=ljA*aa(itypi,itypj)
-      ljxm=ljXs+(-2.0D0*aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
+      ljB=ljA*bb
+      ljA=ljA*aa
+      ljxm=ljXs+(-2.0D0*aa/bb)**(1.0D0/6.0D0)
 
       ssXs=d0cm
       deltat1=1.0d0-om1
@@ -223,7 +325,7 @@ c-------TESTING CODE
 c     Stop and plot energy and derivative as a function of distance
       if (checkstop) then
         ssm=ssC-0.25D0*ssB*ssB/ssA
-        ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
+        ljm=-0.25D0*ljB*bb/aa
         if (ssm.lt.ljm .and.
      &       dabs(rij-0.5d0*(ssxm+ljxm)).lt.0.35d0*(ljxm-ssxm)) then
           nicheck=1000
@@ -248,8 +350,8 @@ c-------END TESTING CODE
         havebond=.false.
         ljd=rij-ljXs
         fac=(1.0D0/ljd)**expon
-        e1=fac*fac*aa(itypi,itypj)
-        e2=fac*bb(itypi,itypj)
+        e1=fac*fac*aa
+        e2=fac*bb
         eij=eps1*eps2rt*eps3rt*(e1+e2)
 C        write(iout,*) eij,'TU?1'
         eps2der=eij*eps3rt
@@ -316,8 +418,8 @@ C         write(iout,*) eij,'TU?3'
           eom12=fac1*d_ssxm(3)+fac2*d_xm(3)+h1*d_ssm(3)
         else
           havebond=.false.
-          ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
-          d_ljm(1)=-0.5D0*bb(itypi,itypj)/aa(itypi,itypj)*ljB
+          ljm=-0.25D0*ljB*bb/aa
+          d_ljm(1)=-0.5D0*bb/aa*ljB
           d_ljm(2)=d_ljm(1)*(0.5D0*eps2rt_om2/eps2rt+alf2/eps3rt)
           d_ljm(3)=d_ljm(1)*(0.5D0*eps1_om12+0.5D0*eps2rt_om12/eps2rt-
      +         alf12/eps3rt)