zmiany w multichainie

diff --git a/source/unres/src_MD-M/chainbuild.F b/source/unres/src_MD-M/chainbuild.F
index c4970d3..766e14a 100644 (file)
--- a/source/unres/src_MD-M/chainbuild.F
+++ b/source/unres/src_MD-M/chainbuild.F
@@ -18,7 +18,8 @@ C Set lprn=.true. for debugging
       lprn = .false.
       perbox=.false.
       fail=.false.
-       if (perbox) then
+      print *, 'enter chainbuild' 
+      if (perbox) then
       cost=dcos(theta(3))
       sint=dsin(theta(3))
       print *,'before refsys'

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 574e19b..a080b08 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
@@ -627,6 +627,12 @@ c     if (icall.eq.0) lprn=.false.
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
+          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)
@@ -657,16 +663,53 @@ c            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
             alf1=alp(itypi)
             alf2=alp(itypj)
             alf12=0.5D0*(alf1+alf2)
-            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)
+          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
+      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)
             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
             rij=dsqrt(rrij)
             sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
-
+            sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj))
             if (sss.lt.1.0d0) then
 
 C Calculate angle-dependent terms of energy and contributions to their
@@ -717,6 +760,7 @@ C Calculate gradient components.
               fac=-expon*(e1+evdwij)*rij_shift
               sigder=fac*sigder
               fac=rij*fac
+            fac=fac+evdwij/(1.0-sss)*(-sssgrad)/sigmaii(itypi,itypj)*rij
 c              fac=0.0d0
 C Calculate the radial part of the gradient
               gg(1)=xj*fac
@@ -765,6 +809,12 @@ c     if (icall.eq.0) lprn=.false.
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
+          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)
@@ -795,16 +845,53 @@ c            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
             alf1=alp(itypi)
             alf2=alp(itypj)
             alf12=0.5D0*(alf1+alf2)
-            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)
+          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
+      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)
             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
             rij=dsqrt(rrij)
             sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
-
+            sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj))
             if (sss.gt.0.0d0) then
 
 C Calculate angle-dependent terms of energy and contributions to their
@@ -855,6 +942,7 @@ C Calculate gradient components.
               fac=-expon*(e1+evdwij)*rij_shift
               sigder=fac*sigder
               fac=rij*fac
+            fac=fac+evdwij/sss*sssgrad/sigmaii(itypi,itypj)*rij
 c              fac=0.0d0
 C Calculate the radial part of the gradient
               gg(1)=xj*fac
@@ -1278,7 +1366,10 @@ C Loop over i,i+2 and i,i+3 pairs of the peptide groups
 C
       do i=iturn3_start,iturn3_end
         if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1
-     &  .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1) cycle
+     &  .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1
+     &  .or. itype(i-1).eq.ntyp1
+     &  .or. itype(i+4).eq.ntyp1
+     &   ) cycle
         dxi=dc(1,i)
         dyi=dc(2,i)
         dzi=dc(3,i)
@@ -1288,6 +1379,12 @@ C
         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
         call eelecij_scale(i,i+2,ees,evdw1,eel_loc)
         if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
@@ -1296,7 +1393,10 @@ C
       do i=iturn4_start,iturn4_end
         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
      &    .or. itype(i+3).eq.ntyp1
-     &    .or. itype(i+4).eq.ntyp1) cycle
+     &    .or. itype(i+4).eq.ntyp1
+     &    .or. itype(i+5).eq.ntyp1
+     &    .or. itype(i-1).eq.ntyp1
+     &    ) cycle
         dxi=dc(1,i)
         dyi=dc(2,i)
         dzi=dc(3,i)
@@ -1316,7 +1416,10 @@ c
 c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
 c
       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
         dxi=dc(1,i)
         dyi=dc(2,i)
         dzi=dc(3,i)
@@ -1326,10 +1429,19 @@ c
         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
 c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
         num_conti=num_cont_hb(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
           call eelecij_scale(i,j,ees,evdw1,eel_loc)
         enddo ! j
         num_cont_hb(i)=num_conti
@@ -1398,16 +1510,54 @@ cd      write (iout,*) "eelecij",i,j
           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
           rrmij=1.0D0/rij
           rij=dsqrt(rij)
           rmij=1.0D0/rij
 c For extracting the short-range part of Evdwpp
           sss=sscale(rij/rpp(iteli,itelj))
-
+          sssgrad=sscagrad(rij/rpp(iteli,itelj))
           r3ij=rrmij*rmij
           r6ij=r3ij*r3ij  
           cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
@@ -1472,9 +1622,9 @@ cgrad            do l=1,3
 cgrad              gelc(l,k)=gelc(l,k)+ggg(l)
 cgrad            enddo
 cgrad          enddo
-          ggg(1)=facvdw*xj
-          ggg(2)=facvdw*yj
-          ggg(3)=facvdw*zj
+          ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
+          ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
+          ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
 c          do k=1,3
 c            ghalf=0.5D0*ggg(k)
 c            gvdwpp(k,i)=gvdwpp(k,i)+ghalf
@@ -1526,9 +1676,12 @@ cgrad              gelc(l,k)=gelc(l,k)+ggg(l)
 cgrad            enddo
 cgrad          enddo
 c 9/28/08 AL Gradient compotents will be summed only at the end
-          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
+          ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
+          ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
+          ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
           do k=1,3
             gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
             gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
@@ -2025,6 +2178,12 @@ c     & " iatel_e_vdw",iatel_e_vdw
         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_vdw(i),
 c     &   ' ielend',ielend_vdw(i)
@@ -2043,13 +2202,51 @@ c     &   ' ielend',ielend_vdw(i)
           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
           rrmij=1.0D0/rij
           rij=dsqrt(rij)
           sss=sscale(rij/rpp(iteli,itelj))
+            sssgrad=sscagrad(rij/rpp(iteli,itelj))
           if (sss.gt.0.0d0) then
             rmij=1.0D0/rij
             r3ij=rrmij*rmij
@@ -2067,9 +2264,12 @@ C
 C Calculate contributions to the Cartesian gradient.
 C
             facvdw=-6*rrmij*(ev1+evdwij)*sss
-            ggg(1)=facvdw*xj
-            ggg(2)=facvdw*yj
-            ggg(3)=facvdw*zj
+          ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
+          ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
+          ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
+C            ggg(1)=facvdw*xj
+C            ggg(2)=facvdw*yj
+C            ggg(3)=facvdw*zj
             do k=1,3
               gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
               gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
@@ -2100,7 +2300,7 @@ C
       dimension ggg(3)
       evdw2=0.0D0
       evdw2_14=0.0d0
-cd    print '(a)','Enter ESCP'
+CD        print '(a)','Enter ESCP KURWA'
 cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
       do i=iatscp_s,iatscp_e
         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
@@ -2108,7 +2308,12 @@ cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
         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))
-
+         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)
@@ -2119,15 +2324,46 @@ 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)
+      dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+      xj_safe=xj
+      yj_safe=yj
+      zj_safe=zj
+      subchap=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
+            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)
 
           sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
-
+          sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
           if (sss.lt.1.0d0) then
-
             fac=rrij**expon2
             e1=fac*fac*aad(itypj,iteli)
             e2=fac*bad(itypj,iteli)
@@ -2143,7 +2379,9 @@ C Uncomment following three lines for Ca-p interactions
 C
 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
 C
+             
             fac=-(evdwij+e1)*rrij*(1.0d0-sss)
+            fac=fac-(evdwij)*sssgrad*dsqrt(rrij)/rscp(itypj,iteli)
             ggg(1)=xj*fac
             ggg(2)=yj*fac
             ggg(3)=zj*fac
@@ -2209,6 +2447,12 @@ cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
         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))
+         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)
 
@@ -2220,13 +2464,43 @@ 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)
+      dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+      xj_safe=xj
+      yj_safe=yj
+      zj_safe=zj
+      subchap=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
+            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)
-
           sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
-
+          sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
           if (sss.gt.0.0d0) then
 
             fac=rrij**expon2
@@ -2245,6 +2519,7 @@ C
 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
 C
             fac=-(evdwij+e1)*rrij*sss
+            fac=fac+(evdwij)*sssgrad*dsqrt(rrij)/rscp(itypj,iteli)
             ggg(1)=xj*fac
             ggg(2)=yj*fac
             ggg(3)=zj*fac

diff --git a/source/unres/src_MD-M/intcartderiv.F b/source/unres/src_MD-M/intcartderiv.F
index f78cadd..0b75bf0 100644 (file)
--- a/source/unres/src_MD-M/intcartderiv.F
+++ b/source/unres/src_MD-M/intcartderiv.F
@@ -751,6 +751,7 @@ c        c(j,1)=c(j,1)
           c(j,i+nres)=c(j,i)+dc(j,i+nres)
         enddo
       enddo
+C      print *,'tutu'
 c      write (iout,*) "CHAINBUILD_CART"
 c      call cartprint
       call int_from_cart1(.false.)