zmiany w WHAM PBC przed sprawdzeniem differow

diff --git a/.gitignore b/.gitignore
index 2e4973f..d0edd88 100644 (file)
--- a/.gitignore
+++ b/.gitignore
@@ -31,3 +31,5 @@ period_build
 build_devel
 build_theta
 build_new
+build_prere/
+period_build2

diff --git a/source/wham/src-M/COMMON.CHAIN b/source/wham/src-M/COMMON.CHAIN
index 8dcdd98..b147f08 100644 (file)
--- a/source/wham/src-M/COMMON.CHAIN
+++ b/source/wham/src-M/COMMON.CHAIN
@@ -12,3 +12,7 @@
      &nsup,nstart_sup,anatemp,
      &nend_sup,chain_length,tabperm(maxperm,maxsym),nperm,
      & nstart_seq,ishift_pdb
+      double precision boxxsize,boxysize,boxzsize,enecut,sscut,sss,
+     &sssgrad
+      common /box/  boxxsize,boxysize,boxzsize,enecut,sscut,sss,sssgrad
+

diff --git a/source/wham/src-M/COMMON.SPLITELE b/source/wham/src-M/COMMON.SPLITELE
new file mode 100644 (file)
index 0000000..a2f0447
--- /dev/null
+++ b/source/wham/src-M/COMMON.SPLITELE
@@ -0,0 +1,2 @@
+      double precision r_cut,rlamb
+      common /splitele/ r_cut,rlamb

diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F
index 8f932cc..a4f5fb4 100644 (file)
--- a/source/wham/src-M/energy_p_new.F
+++ b/source/wham/src-M/energy_p_new.F
@@ -797,6 +797,14 @@ c      if (icall.gt.0) lprn=.true.
         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
+
         dxi=dc_norm(1,nres+i)
         dyi=dc_norm(2,nres+i)
         dzi=dc_norm(3,nres+i)
@@ -830,15 +838,59 @@ c           chip12=0.0D0
 c           alf1=0.0D0
 c           alf2=0.0D0
 c           alf12=0.0D0
-            xj=c(1,nres+j)-xi
-            yj=c(2,nres+j)-yi
-            zj=c(3,nres+j)-zi
+            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
+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
+
             dxj=dc_norm(1,nres+j)
             dyj=dc_norm(2,nres+j)
             dzj=dc_norm(3,nres+j)
 c            write (iout,*) i,j,xj,yj,zj
             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
             rij=dsqrt(rrij)
+            sss=sscale((1.0d0/rij)/sigma(itypi,itypj))
+            sssgrad=sscagrad((1.0d0/rij)/sigma(itypi,itypj))
+            if (sss.le.0.0) cycle
 C Calculate angle-dependent terms of energy and contributions to their
 C derivatives.
             call sc_angular
@@ -861,9 +913,9 @@ c---------------------------------------------------------------
             eps3der=evdwij*eps2rt
             evdwij=evdwij*eps2rt*eps3rt
             if (bb(itypi,itypj).gt.0) then
-              evdw=evdw+evdwij
+              evdw=evdw+evdwij*sss
             else
-              evdw_t=evdw_t+evdwij
+              evdw_t=evdw_t+evdwij*sss
             endif
             ij=icant(itypi,itypj)
             aux=eps1*eps2rt**2*eps3rt**2
@@ -892,6 +944,7 @@ C Calculate gradient components.
             fac=-expon*(e1+evdwij)*rij_shift
             sigder=fac*sigder
             fac=rij*fac
+            fac=fac+evdwij/sss*sssgrad/sigma(itypi,itypj)*rij
 C Calculate the radial part of the gradient
             gg(1)=xj*fac
             gg(2)=yj*fac
@@ -1852,7 +1905,10 @@ cd      write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
         gcorr_loc(i)=0.0d0
       enddo
       do i=iatel_s,iatel_e
-        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
+     &  .or. itype(i+2).eq.ntyp1
+     &  .or. itype(i-1).eq.ntyp1
+     &) cycle
         if (itel(i).eq.0) goto 1215
         dxi=dc(1,i)
         dyi=dc(2,i)
@@ -1863,10 +1919,21 @@ cd      write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
         xmedi=c(1,i)+0.5d0*dxi
         ymedi=c(2,i)+0.5d0*dyi
         zmedi=c(3,i)+0.5d0*dzi
+          xmedi=mod(xmedi,boxxsize)
+          if (xmedi.lt.0) xmedi=xmedi+boxxsize
+          ymedi=mod(ymedi,boxysize)
+          if (ymedi.lt.0) ymedi=ymedi+boxysize
+          zmedi=mod(zmedi,boxzsize)
+          if (zmedi.lt.0) zmedi=zmedi+boxzsize
         num_conti=0
 c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
         do j=ielstart(i),ielend(i)
-          if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle
+          if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
+     & .or.itype(j+2).eq.ntyp1
+     & .or.itype(j-1).eq.ntyp1
+     &) cycle
+C
+C) cycle
           if (itel(j).eq.0) goto 1216
           ind=ind+1
           iteli=itel(i)
@@ -1888,10 +1955,50 @@ C End diagnostics
           dx_normj=dc_norm(1,j)
           dy_normj=dc_norm(2,j)
           dz_normj=dc_norm(3,j)
-          xj=c(1,j)+0.5D0*dxj-xmedi
-          yj=c(2,j)+0.5D0*dyj-ymedi
-          zj=c(3,j)+0.5D0*dzj-zmedi
+          xj=c(1,j)+0.5D0*dxj
+          yj=c(2,j)+0.5D0*dyj
+          zj=c(3,j)+0.5D0*dzj
+         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
+      dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+      xj_safe=xj
+      yj_safe=yj
+      zj_safe=zj
+      isubchap=0
+      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
+            isubchap=1
+          endif
+       enddo
+       enddo
+       enddo
+       if (isubchap.eq.1) then
+          xj=xj_temp-xmedi
+          yj=yj_temp-ymedi
+          zj=zj_temp-zmedi
+       else
+          xj=xj_safe-xmedi
+          yj=yj_safe-ymedi
+          zj=zj_safe-zmedi
+       endif
+
           rij=xj*xj+yj*yj+zj*zj
+            sss=sscale(sqrt(rij))
+            sssgrad=sscagrad(sqrt(rij))
           rrmij=1.0D0/rij
           rij=dsqrt(rij)
           rmij=1.0D0/rij
@@ -1915,7 +2022,7 @@ c          write (iout,*) "i",i,iteli," j",j,itelj," eesij",eesij
 C 12/26/95 - for the evaluation of multi-body H-bonding interactions
           ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
           ees=ees+eesij
-          evdw1=evdw1+evdwij
+          evdw1=evdw1+evdwij*sss
 c             write (iout,'(a6,2i5,0pf7.3,2i5,2e11.3)') 
 c     &'evdw1',i,j,evdwij
 c     &,iteli,itelj,aaa,evdw1
@@ -1929,7 +2036,7 @@ C
 C Calculate contributions to the Cartesian gradient.
 C
 #ifdef SPLITELE
-          facvdw=-6*rrmij*(ev1+evdwij) 
+          facvdw=-6*rrmij*(ev1+evdwij)*sss
           facel=-3*rrmij*(el1+eesij)
           fac1=fac
           erij(1)=xj*rmij
@@ -1955,9 +2062,18 @@ C
               gelc(l,k)=gelc(l,k)+ggg(l)
             enddo
           enddo
-          ggg(1)=facvdw*xj
-          ggg(2)=facvdw*yj
-          ggg(3)=facvdw*zj
+C          ggg(1)=facvdw*xj
+C          ggg(2)=facvdw*yj
+C          ggg(3)=facvdw*zj
+          if (sss.gt.0.0) then
+          ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj
+          ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj
+          ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj
+          else
+          ggg(1)=0.0
+          ggg(2)=0.0
+          ggg(3)=0.0
+          endif
           do k=1,3
             ghalf=0.5D0*ggg(k)
             gvdwpp(k,i)=gvdwpp(k,i)+ghalf
@@ -1972,7 +2088,7 @@ C
             enddo
           enddo
 #else
-          facvdw=ev1+evdwij 
+          facvdw=(ev1+evdwij)*sss
           facel=el1+eesij  
           fac1=fac
           fac=-3*rrmij*(facvdw+facvdw+facel)
@@ -2818,7 +2934,13 @@ c     &   " iscp",(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i))
         xi=0.5D0*(c(1,i)+c(1,i+1))
         yi=0.5D0*(c(2,i)+c(2,i+1))
         zi=0.5D0*(c(3,i)+c(3,i+1))
-
+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
         do iint=1,nscp_gr(i)
 
         do j=iscpstart(i,iint),iscpend(i,iint)
@@ -2829,28 +2951,73 @@ c         xj=c(1,nres+j)-xi
 c         yj=c(2,nres+j)-yi
 c         zj=c(3,nres+j)-zi
 C Uncomment following three lines for Ca-p interactions
-          xj=c(1,j)-xi
-          yj=c(2,j)-yi
-          zj=c(3,j)-zi
+          xj=c(1,j)
+          yj=c(2,j)
+          zj=c(3,j)
+C returning the 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
+      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 jth atom
+      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
           rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+C sss is scaling function for smoothing the cutoff gradient otherwise
+C the gradient would not be continuouse
+          sss=sscale(1.0d0/(dsqrt(rrij)))
+          if (sss.le.0.0d0) cycle
+          sssgrad=sscagrad(1.0d0/(dsqrt(rrij)))
           fac=rrij**expon2
           e1=fac*fac*aad(itypj,iteli)
           e2=fac*bad(itypj,iteli)
           if (iabs(j-i) .le. 2) then
             e1=scal14*e1
             e2=scal14*e2
-            evdw2_14=evdw2_14+e1+e2
+            evdw2_14=evdw2_14+(e1+e2)*sss
           endif
           evdwij=e1+e2
 c          write (iout,'(a6,2i5,0pf7.3,2i3,3e11.3)')
 c     &        'evdw2',i,j,evdwij,iteli,itypj,fac,aad(itypj,iteli),
 c     &       bad(itypj,iteli)
-          evdw2=evdw2+evdwij
+          evdw2=evdw2+evdwij*sss
           if (calc_grad) then
 C
 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
 C
-          fac=-(evdwij+e1)*rrij
+          fac=-(evdwij+e1)*rrij*sss
+          fac=fac+(evdwij)*sssgrad*dsqrt(rrij)/expon
           ggg(1)=xj*fac
           ggg(2)=yj*fac
           ggg(3)=zj*fac
@@ -3085,16 +3252,19 @@ c
       estr1=0.0d0
 c      write (iout,*) "distchainmax",distchainmax
       do i=nnt+1,nct
-        if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then
-          estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax)
-          do j=1,3
-          gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax)
-     &      *dc(j,i-1)/vbld(i)
-          enddo
-          if (energy_dec) write(iout,*)
-     &       "estr1",i,vbld(i),distchainmax,
-     &       gnmr1(vbld(i),-1.0d0,distchainmax)
-        else
+        if (itype(i-1).eq.ntyp1 .and. itype(i).eq.ntyp1) cycle
+C          estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax)
+C          do j=1,3
+C          gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax)
+C     &      *dc(j,i-1)/vbld(i)
+C          enddo
+C          if (energy_dec) write(iout,*)
+C     &       "estr1",i,vbld(i),distchainmax,
+C     &       gnmr1(vbld(i),-1.0d0,distchainmax)
+C        else
+         if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then
+        diff = vbld(i)-vbldpDUM
+         else
           diff = vbld(i)-vbldp0
 c          write (iout,*) i,vbld(i),vbldp0,diff,AKP*diff*diff
           estr=estr+diff*diff
@@ -3186,7 +3356,9 @@ c      write (iout,*) "nres",nres
 c     write (*,'(a,i2)') 'EBEND ICG=',icg
 c      write (iout,*) ithet_start,ithet_end
       do i=ithet_start,ithet_end
-        if (itype(i-1).eq.ntyp1) cycle
+C        if (itype(i-1).eq.ntyp1) cycle
+        if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1
+     &  .or.itype(i).eq.ntyp1) cycle
 C Zero the energy function and its derivative at 0 or pi.
         call splinthet(theta(i),0.5d0*delta,ss,ssd)
         it=itype(i-1)
@@ -3203,7 +3375,7 @@ C Zero the energy function and its derivative at 0 or pi.
           ichir22=isign(1,itype(i))
          endif
 
-        if (i.gt.3 .and. itype(i-2).ne.ntyp1) then
+        if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
 #ifdef OSF
           phii=phi(i)
           icrc=0
@@ -3218,7 +3390,7 @@ C Zero the energy function and its derivative at 0 or pi.
           y(1)=0.0D0
           y(2)=0.0D0
         endif
-        if (i.lt.nres .and. itype(i).ne.ntyp1) then
+        if (i.lt.nres .and. itype(i+1).ne.ntyp1) then
 #ifdef OSF
           phii1=phi(i+1)
           icrc=0
@@ -3433,7 +3605,9 @@ C
       etheta=0.0D0
 c      write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
       do i=ithet_start,ithet_end
-        if (itype(i-1).eq.ntyp1) cycle
+C        if (itype(i-1).eq.ntyp1) cycle
+        if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1
+     &  .or.itype(i).eq.ntyp1) cycle
         if (iabs(itype(i+1)).eq.20) iblock=2
         if (iabs(itype(i+1)).ne.20) iblock=1
         dethetai=0.0d0
@@ -3445,7 +3619,7 @@ c      write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
           coskt(k)=dcos(k*theti2)
           sinkt(k)=dsin(k*theti2)
         enddo
-        if (i.gt.3 .and. itype(i-2).ne.ntyp1) then
+        if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
 #ifdef OSF
           phii=phi(i)
           if (phii.ne.phii) phii=150.0
@@ -3465,7 +3639,7 @@ c      write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
             sinph1(k)=0.0d0
           enddo 
         endif
-        if (i.lt.nres .and. itype(i).ne.ntyp1) then
+        if (i.lt.nres .and. itype(i+1).ne.ntyp1) then
 #ifdef OSF
           phii1=phi(i+1)
           if (phii1.ne.phii1) phii1=150.0
@@ -4401,8 +4575,10 @@ C Set lprn=.true. for debugging
 c      lprn=.true.
       etors=0.0D0
       do i=iphi_start,iphi_end
-        if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1
-     &       .or. itype(i).eq.ntyp1) cycle
+        if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1
+     &      .or. itype(i).eq.ntyp1 .or. itype(i-3).eq.ntyp1) cycle
+C        if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1
+C     &       .or. itype(i).eq.ntyp1) cycle
         if (itel(i-2).eq.0 .or. itel(i-1).eq.0) goto 1215
          if (iabs(itype(i)).eq.20) then
          iblock=2
@@ -4502,8 +4678,11 @@ C Set lprn=.true. for debugging
 c     lprn=.true.
       etors_d=0.0D0
       do i=iphi_start,iphi_end-1
-        if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1
-     &      .or. itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+C        if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1
+C     &      .or. itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+         if ((itype(i-2).eq.ntyp1).or.itype(i-3).eq.ntyp1.or.
+     &  (itype(i-1).eq.ntyp1).or.(itype(i).eq.ntyp1).or.
+     &  (itype(i+1).eq.ntyp1)) cycle
         if (itel(i-2).eq.0 .or. itel(i-1).eq.0 .or. itel(i).eq.0) 
      &     goto 1215
         itori=itortyp(itype(i-2))
@@ -7512,4 +7691,34 @@ C-----------------------------------------------------------------------------
       scalar=sc
       return
       end
+C-----------------------------------------------------------------------
+      double precision function sscale(r)
+      double precision r,gamm
+      include "COMMON.SPLITELE"
+      if(r.lt.r_cut-rlamb) then
+        sscale=1.0d0
+      else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
+        gamm=(r-(r_cut-rlamb))/rlamb
+        sscale=1.0d0+gamm*gamm*(2*gamm-3.0d0)
+      else
+        sscale=0d0
+      endif
+      return
+      end
+C-----------------------------------------------------------------------
+C-----------------------------------------------------------------------
+      double precision function sscagrad(r)
+      double precision r,gamm
+      include "COMMON.SPLITELE"
+      if(r.lt.r_cut-rlamb) then
+        sscagrad=0.0d0
+      else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
+        gamm=(r-(r_cut-rlamb))/rlamb
+        sscagrad=gamm*(6*gamm-6.0d0)/rlamb
+      else
+        sscagrad=0.0d0
+      endif
+      return
+      end
+C-----------------------------------------------------------------------
 

diff --git a/source/wham/src-M/parmread.F b/source/wham/src-M/parmread.F
index c9ac2a3..ea6cb14 100644 (file)
--- a/source/wham/src-M/parmread.F
+++ b/source/wham/src-M/parmread.F
@@ -150,7 +150,7 @@ c Read the virtual-bond parameters, masses, and moments of inertia
 c and Stokes' radii of the peptide group and side chains
 c
 #ifdef CRYST_BOND
-      read (ibond,*) vbldp0,akp
+      read (ibond,*) vbldp0,vbldpdum,akp
       do i=1,ntyp
         nbondterm(i)=1
         read (ibond,*) vbldsc0(1,i),aksc(1,i)
@@ -162,7 +162,7 @@ c
         endif
       enddo
 #else
-      read (ibond,*) ijunk,vbldp0,akp,rjunk
+      read (ibond,*) ijunk,vbldp0,vbldpdum,akp,rjunk
       do i=1,ntyp
         read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
      &   j=1,nbondterm(i))

diff --git a/source/wham/src-M/readpdb.f b/source/wham/src-M/readpdb.f
index a65709e..daddd52 100644 (file)
--- a/source/wham/src-M/readpdb.f
+++ b/source/wham/src-M/readpdb.f
@@ -26,7 +26,9 @@ C geometry.
           goto 10
         else if (card(:3).eq.'TER') then
 C End current chain
-          ires_old=ires+1 
+c          ires_old=ires+1 
+          ires_old=ires+2
+          itype(ires_old-1)=ntyp1 
           itype(ires_old)=ntyp1
           ibeg=2
 c          write (iout,*) "Chain ended",ires,ishift,ires_old
@@ -85,14 +87,51 @@ C system
       nres=ires
       do i=2,nres-1
 c        write (iout,*) i,itype(i)
+
         if (itype(i).eq.ntyp1) then
-c          write (iout,*) "dummy",i,itype(i)
-          do j=1,3
-            c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
-c            c(j,i)=(c(j,i-1)+c(j,i+1))/2
-            dc(j,i)=c(j,i)
-          enddo
-        endif
+         if (itype(i+1).eq.ntyp1) then
+C 16/01/2014 by Adasko: Adding to dummy atoms in the chain
+C first is connected prevous chain (itype(i+1).eq.ntyp1)=true
+C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
+C           if (unres_pdb) then
+C 2/15/2013 by Adam: corrected insertion of the last dummy residue
+C            call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
+C            if (fail) then
+C              e2(1)=0.0d0
+C              e2(2)=1.0d0
+C              e2(3)=0.0d0
+C            endif !fail
+C            do j=1,3
+C             c(j,i)=c(j,i-1)-1.9d0*e2(j)
+C            enddo
+C           else   !unres_pdb
+           do j=1,3
+             dcj=(c(j,i-2)-c(j,i-3))/2.0
+             c(j,i)=c(j,i-1)+dcj
+             c(j,nres+i)=c(j,i)
+           enddo     
+C          endif   !unres_pdb
+         else     !itype(i+1).eq.ntyp1
+C          if (unres_pdb) then
+C 2/15/2013 by Adam: corrected insertion of the first dummy residue
+C            call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
+C            if (fail) then
+C              e2(1)=0.0d0
+C              e2(2)=1.0d0
+C              e2(3)=0.0d0
+C            endif
+C            do j=1,3
+C              c(j,i)=c(j,i+1)-1.9d0*e2(j)
+C            enddo
+C          else !unres_pdb
+           do j=1,3
+            dcj=(c(j,i+3)-c(j,i+2))/2.0
+            c(j,i)=c(j,i+1)-dcj
+            c(j,nres+i)=c(j,i)
+           enddo
+C          endif !unres_pdb
+         endif !itype(i+1).eq.ntyp1
+        endif  !itype.eq.ntyp1
       enddo
 C Calculate the CM of the last side chain.
       call sccenter(ires,iii,sccor)
@@ -102,7 +141,7 @@ C Calculate the CM of the last side chain.
         nres=nres+1
         itype(nres)=ntyp1
         do j=1,3
-          dcj=c(j,nres-2)-c(j,nres-3)
+          dcj=(c(j,nres-2)-c(j,nres-3))/2.0
           c(j,nres)=c(j,nres-1)+dcj
           c(j,2*nres)=c(j,nres)
         enddo
@@ -120,7 +159,7 @@ C Calculate the CM of the last side chain.
         nsup=nsup-1
         nstart_sup=2
         do j=1,3
-          dcj=c(j,4)-c(j,3)
+          dcj=(c(j,4)-c(j,3))/2.0
           c(j,1)=c(j,2)-dcj
           c(j,nres+1)=c(j,1)
         enddo

diff --git a/source/wham/src-M/readrtns.F b/source/wham/src-M/readrtns.F
index 36c13b1..9023170 100644 (file)
--- a/source/wham/src-M/readrtns.F
+++ b/source/wham/src-M/readrtns.F
@@ -17,6 +17,7 @@
       include "COMMON.FREE"
       include "COMMON.CONTROL"
       include "COMMON.ENERGIES"
+      include "COMMON.SPLITELE"
       character*800 controlcard
       integer i,j,k,ii,n_ene_found
       integer ind,itype1,itype2,itypf,itypsc,itypp
@@ -74,6 +75,12 @@
       call readi(controlcard,"RESCALE",rescale_mode,1)
       check_conf=index(controlcard,"NO_CHECK_CONF").eq.0
       call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
+       call reada(controlcard,'BOXX',boxxsize,100.0d0)
+       call reada(controlcard,'BOXY',boxysize,100.0d0)
+       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
+c Cutoff range for interactions
+      call reada(controlcard,"R_CUT",r_cut,15.0d0)
+      call reada(controlcard,"LAMBDA",rlamb,0.3d0)
       call readi(controlcard,'SYM',symetr,1)
       write (iout,*) "DISTCHAINMAX",distchainmax
       write (iout,*) "delta",delta