After long DEBUG of energy function

diff --git a/source/unres/src_MD-M/elecont.f b/source/unres/src_MD-M/elecont.f
index 4faf3e5..73325f2 100644 (file)
--- a/source/unres/src_MD-M/elecont.f
+++ b/source/unres/src_MD-M/elecont.f
@@ -52,6 +52,12 @@ c      data epp    / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/
         xmedi=xi+0.5*dxi
         ymedi=yi+0.5*dyi
         zmedi=zi+0.5*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
         do 4 j=i+2,nct-1
           if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) goto 4
           ind=ind+1
@@ -66,9 +72,49 @@ c      data epp    / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/
           dxj=c(1,j+1)-c(1,j)
           dyj=c(2,j+1)-c(2,j)
           dzj=c(3,j+1)-c(3,j)
-          xj=c(1,j)+0.5*dxj-xmedi
-          yj=c(2,j)+0.5*dyj-ymedi
-          zj=c(3,j)+0.5*dzj-zmedi
+          xj=c(1,j)+0.5*dxj
+          yj=c(2,j)+0.5*dyj
+          zj=c(3,j)+0.5*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.0/(xj*xj+yj*yj+zj*zj)
           rmij=sqrt(rrmij)
           r3ij=rrmij*rmij
@@ -94,7 +140,7 @@ c      data epp    / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/
             econt(ncont)=eesij
           endif
           ees=ees+eesij
-          evdw=evdw+evdwij
+          evdw=evdw+evdwij*sss
     4   continue
     1 continue
       if (lprint) then

diff --git a/source/unres/src_MD-M/energy_p_new-sep.F b/source/unres/src_MD-M/energy_p_new-sep.F
index 85f3e1f..de386c6 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new-sep.F
+++ b/source/unres/src_MD-M/energy_p_new-sep.F
@@ -2128,6 +2128,12 @@ c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
         xmedi=c(1,i)+0.5d0*dxi
         ymedi=c(2,i)+0.5d0*dyi
         zmedi=c(3,i)+0.5d0*dzi
+          xmedi=mod(xmedi,boxysize)
+          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)
@@ -2143,13 +2149,50 @@ c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(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
-          rij=xj*xj+yj*yj+zj*zj
+          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)
+      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(sqrt(rij))
           if (sss.lt.1.0d0) then
             rmij=1.0D0/rij
             r3ij=rrmij*rmij
@@ -2168,9 +2211,9 @@ C
 C Calculate contributions to the Cartesian gradient.
 C
             facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)
-            ggg(1)=facvdw*xj
-            ggg(2)=facvdw*yj
-            ggg(3)=facvdw*zj
+            ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj
+            ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj
+            ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj
 
             do k=1,3
               ghalf=0.5D0*ggg(k)
@@ -2227,6 +2270,12 @@ c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
         xmedi=c(1,i)+0.5d0*dxi
         ymedi=c(2,i)+0.5d0*dyi
         zmedi=c(3,i)+0.5d0*dzi
+          xmedi=mod(xmedi,boxysize)
+          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)
@@ -2242,13 +2291,50 @@ c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(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)
+     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(sqrt(rij))
           if (sss.gt.0.0d0) then
             rmij=1.0D0/rij
             r3ij=rrmij*rmij
@@ -2267,9 +2353,9 @@ 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
+            ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj
+            ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj
 
             do k=1,3
               ghalf=0.5D0*ggg(k)

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 a080b08..7d2d27f 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
@@ -1311,6 +1311,7 @@ cd        write(iout,*) 'EE',EE(:,:,i)
 cd      enddo
 cd      call check_vecgrad
 cd      stop
+C      print *,"WCHODZE3"
       if (icheckgrad.eq.1) then
         do i=1,nres-1
           fac=1.0d0/dsqrt(scalar(dc(1,i),dc(1,i)))
@@ -1406,6 +1407,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=num_cont_hb(i)
         call eelecij_scale(i,i+3,ees,evdw1,eel_loc)
         if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1) 
@@ -1496,6 +1503,7 @@ C 13-go grudnia roku pamietnego...
      &                   0.0d0,0.0d0,1.0d0/
 c          time00=MPI_Wtime()
 cd      write (iout,*) "eelecij",i,j
+C      print *,"WCHODZE2"
           ind=ind+1
           iteli=itel(i)
           itelj=itel(j)
@@ -1519,7 +1527,7 @@ cd      write (iout,*) "eelecij",i,j
           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
+      dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
       xj_safe=xj
       yj_safe=yj
       zj_safe=zj
@@ -1530,7 +1538,7 @@ cd      write (iout,*) "eelecij",i,j
           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
+          dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
           if(dist_temp.lt.dist_init) then
             dist_init=dist_temp
             xj_temp=xj
@@ -2164,6 +2172,7 @@ c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
       double precision scal_el /0.5d0/
 #endif
       evdw1=0.0D0
+C      print *,"WCHODZE"
 c      write (iout,*) "iatel_s_vdw",iatel_s_vdw,
 c     & " iatel_e_vdw",iatel_e_vdw
       call flush(iout)
@@ -2211,7 +2220,7 @@ c     &   ' ielend',ielend_vdw(i)
           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
+      dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
       xj_safe=xj
       yj_safe=yj
       zj_safe=zj
@@ -2222,7 +2231,7 @@ c     &   ' ielend',ielend_vdw(i)
           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
+          dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
           if(dist_temp.lt.dist_init) then
             dist_init=dist_temp
             xj_temp=xj

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 23aafae..ad328c0 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new_barrier.F
@@ -2010,7 +2010,7 @@ C
       include 'COMMON.VECTORS'
       include 'COMMON.FFIELD'
       dimension ggg(3)
-cd      write(iout,*) 'In EELEC_soft_sphere'
+C      write(iout,*) 'In EELEC_soft_sphere'
       ees=0.0D0
       evdw1=0.0D0
       eel_loc=0.0d0 
@@ -2025,6 +2025,12 @@ cd      write(iout,*) 'In EELEC_soft_sphere'
         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)
@@ -2038,10 +2044,49 @@ c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
           dxj=dc(1,j)
           dyj=dc(2,j)
           dzj=dc(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-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**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))
           if (rij.lt.r0ijsq) then
             evdw1ij=0.25d0*(rij-r0ijsq)**2
             fac=rij-r0ijsq
@@ -2049,13 +2094,13 @@ c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
             evdw1ij=0.0d0
             fac=0.0d0
           endif
-          evdw1=evdw1+evdw1ij
+          evdw1=evdw1+evdw1ij*sss
 C
 C Calculate contributions to the Cartesian gradient.
 C
-          ggg(1)=fac*xj
-          ggg(2)=fac*yj
-          ggg(3)=fac*zj
+          ggg(1)=fac*xj*sssgrad
+          ggg(2)=fac*yj*sssgrad
+          ggg(3)=fac*zj*sssgrad
           do k=1,3
             gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
             gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
@@ -2928,31 +2973,6 @@ C 14/01/2014 TURN3,TUNR4 does no go under periodic boundry condition
         xmedi=c(1,i)+0.5d0*dxi
         ymedi=c(2,i)+0.5d0*dyi
         zmedi=c(3,i)+0.5d0*dzi
-C Return atom into box, boxxsize is size of box in x dimension
-c  184   continue
-c        if (xmedi.gt.((0.5d0)*boxxsize)) xmedi=xmedi-boxxsize
-c        if (xmedi.lt.((-0.5d0)*boxxsize)) xmedi=xmedi+boxxsize
-C Condition for being inside the proper box
-c        if ((xmedi.gt.((0.5d0)*boxxsize)).or.
-c     &       (xmedi.lt.((-0.5d0)*boxxsize))) then
-c        go to 184
-c        endif
-c  185   continue
-c        if (ymedi.gt.((0.5d0)*boxysize)) ymedi=ymedi-boxysize
-c        if (ymedi.lt.((-0.5d0)*boxysize)) ymedi=ymedi+boxysize
-cC Condition for being inside the proper box
-c        if ((ymedi.gt.((0.5d0)*boxysize)).or.
-c     &       (ymedi.lt.((-0.5d0)*boxysize))) then
-c        go to 185
-c        endif
-c  186   continue
-c        if (zmedi.gt.((0.5d0)*boxzsize)) zmedi=zmedi-boxzsize
-c        if (zmedi.lt.((-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize
-cC Condition for being inside the proper box
-c        if ((zmedi.gt.((0.5d0)*boxzsize)).or.
-c     &       (zmedi.lt.((-0.5d0)*boxzsize))) then
-c        go to 186
-c        endif
           xmedi=mod(xmedi,boxxsize)
           if (xmedi.lt.0) xmedi=xmedi+boxxsize
           ymedi=mod(ymedi,boxysize)
@@ -3169,7 +3189,8 @@ C          zj=c(3,j)+0.5D0*dzj-zmedi
           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
+          if ((zj.lt.0).or.(xj.lt.0).or.(yj.lt.0)) write (*,*) "CHUJ"
+      dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
       xj_safe=xj
       yj_safe=yj
       zj_safe=zj
@@ -3180,7 +3201,7 @@ C          zj=c(3,j)+0.5D0*dzj-zmedi
           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
+          dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
           if(dist_temp.lt.dist_init) then
             dist_init=dist_temp
             xj_temp=xj