correction in restart - adding dc0 -critical bug
[unres.git] / source / unres / src_MD-M / MREMD.F
index a3ec5a7..7a72b17 100644 (file)
@@ -572,10 +572,10 @@ C            call enerprint(remd_ene(0,i))
            endif
            open(irest2,file=rest2name,status='unknown')
            write(irest2,*) totT,EK,potE,totE,t_bath
-           do i=1,2*nres
+           do i=0,2*nres
             write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
            enddo
-           do i=1,2*nres
+           do i=0,2*nres
             write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
            enddo
            if(usampl) then
@@ -1232,8 +1232,8 @@ c-----------------------------------------------------------------------
       include 'COMMON.SBRIDGE'
       include 'COMMON.INTERACT'
                
-      real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
-     &     d_restart2(3,2*maxres*maxprocs)
+      real d_restart1(3,2*(maxres+1)*maxprocs),r_d(3,2*(maxres+1)),
+     &     d_restart2(3,2*(maxres+1)*maxprocs)
       real t5_restart1(5)
       integer iret,itmp
       integer*2 i_index
@@ -1250,23 +1250,23 @@ c-----------------------------------------------------------------------
      &      t_restart1,5,mpi_real,king,CG_COMM,ierr)
 
 
-       do i=1,2*nres
+       do i=0,2*nres
          do j=1,3
            r_d(j,i)=d_t(j,i)
          enddo
        enddo
-       call mpi_gather(r_d,3*2*nres,mpi_real,
-     &           d_restart1,3*2*nres,mpi_real,king,
+       call mpi_gather(r_d,3*2*(nres+1),mpi_real,
+     &           d_restart1,3*2*(nres+1),mpi_real,king,
      &           CG_COMM,ierr)
 
 
-       do i=1,2*nres
+       do i=0,2*nres
          do j=1,3
            r_d(j,i)=dc(j,i)
          enddo
        enddo
-       call mpi_gather(r_d,3*2*nres,mpi_real,
-     &           d_restart2,3*2*nres,mpi_real,king,
+       call mpi_gather(r_d,3*2*(nres+1),mpi_real,
+     &           d_restart2,3*2*(nres+1),mpi_real,king,
      &           CG_COMM,ierr)
 
        if(me.eq.king) then
@@ -1295,16 +1295,18 @@ c-----------------------------------------------------------------------
          enddo
 
          do il=0,nodes-1
-           do i=1,2*nres
+           do i=0,2*nres
             do j=1,3
-             call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
+             call 
+     &      xdrffloat_(ixdrf, d_restart1(j,i+2*(nres+1)*il), iret)
             enddo
            enddo
          enddo
          do il=0,nodes-1
-           do i=1,2*nres
+           do i=0,2*nres
             do j=1,3
-             call xdrffloat_(ixdrf, d_restart2(j,i+2*nres*il), iret)
+             call 
+     & xdrffloat_(ixdrf, d_restart2(j,i+2*(nres+1)*il), iret)
             enddo
            enddo
          enddo
@@ -1355,16 +1357,18 @@ c-----------------------------------------------------------------------
          enddo
 
          do il=0,nodes-1
-           do i=1,2*nres
+           do i=0,2*nres
             do j=1,3
-             call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
+             call 
+     &    xdrffloat(ixdrf, d_restart1(j,i+2*(nres+1)*il), iret)
             enddo
            enddo
          enddo
          do il=0,nodes-1
-           do i=1,2*nres
+           do i=0,2*nres
             do j=1,3
-             call xdrffloat(ixdrf, d_restart2(j,i+2*nres*il), iret)
+             call
+     &  xdrffloat(ixdrf, d_restart2(j,i+2*(nres+1)*il), iret)
             enddo
            enddo
          enddo
@@ -1631,7 +1635,7 @@ c end debugging
       include 'COMMON.CHAIN'
       include 'COMMON.SBRIDGE'
       include 'COMMON.INTERACT'
-      real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
+      real d_restart1(3,2*(maxres+1)*maxprocs),r_d(3,2*(maxres+1)),
      &                 t5_restart1(5)
       integer*2 i_index
      &            (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
@@ -1707,45 +1711,47 @@ c end debugging
 
          if(me.eq.king)then
               do il=0,nodes-1
-               do i=1,2*nres
+               do i=0,2*nres
 c                read(irest2,'(3e15.5)') 
 c     &                (d_restart1(j,i+2*nres*il),j=1,3)
             do j=1,3
 #ifdef AIX
-             call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
+             call xdrffloat_(ixdrf, d_restart1(j,i+2*(nres+1)*il), iret)
 #else
-             call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
+             call xdrffloat(ixdrf, d_restart1(j,i+2*(nres+1)*il), iret)
 #endif
             enddo
                enddo
               enddo
          endif
-         call mpi_scatter(d_restart1,3*2*nres,mpi_real,
-     &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
+         call mpi_scatter(d_restart1,3*2*(nres+1),mpi_real,
+     &           r_d,3*2*(nres+1),mpi_real,king,CG_COMM,ierr)
 
-         do i=1,2*nres
+         do i=0,2*nres
            do j=1,3
             d_t(j,i)=r_d(j,i)
            enddo
          enddo
          if(me.eq.king)then 
               do il=0,nodes-1
-               do i=1,2*nres
+               do i=0,2*nres
 c                read(irest2,'(3e15.5)') 
 c     &                (d_restart1(j,i+2*nres*il),j=1,3)
             do j=1,3
 #ifdef AIX
-             call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
+             call 
+     &   xdrffloat_(ixdrf, d_restart1(j,i+2*(nres+1)*il), iret)
 #else
-             call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
+             call 
+     &   xdrffloat(ixdrf, d_restart1(j,i+2*(nres+1)*il), iret)
 #endif
             enddo
                enddo
               enddo
          endif
-         call mpi_scatter(d_restart1,3*2*nres,mpi_real,
-     &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
-         do i=1,2*nres
+         call mpi_scatter(d_restart1,3*2*(nres+1),mpi_real,
+     &           r_d,3*2*(nres+1),mpi_real,king,CG_COMM,ierr)
+         do i=0,2*nres
            do j=1,3
             dc(j,i)=r_d(j,i)
            enddo
@@ -1841,14 +1847,14 @@ c     &                (d_restart1(j,i+2*nres*il),j=1,3)
 
          if(me.eq.king)then
               do il=0,nodes-1
-               do i=1,2*nres
+               do i=0,2*nres
                 read(irest2,'(3e15.5)') 
-     &                (d_restart1(j,i+2*nres*il),j=1,3)
+     &                (d_restart1(j,i+2*(nres+1)*il),j=1,3)
                enddo
               enddo
          endif
-         call mpi_scatter(d_restart1,3*2*nres,mpi_real,
-     &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
+         call mpi_scatter(d_restart1,3*2*(nres+1),mpi_real,
+     &           r_d,3*2*(nres+1),mpi_real,king,CG_COMM,ierr)
 
          do i=1,2*nres
            do j=1,3
@@ -1857,14 +1863,14 @@ c     &                (d_restart1(j,i+2*nres*il),j=1,3)
          enddo
          if(me.eq.king)then 
               do il=0,nodes-1
-               do i=1,2*nres
+               do i=0,2*nres
                 read(irest2,'(3e15.5)') 
-     &                (d_restart1(j,i+2*nres*il),j=1,3)
+     &                (d_restart1(j,i+2*(nres+1)*il),j=1,3)
                enddo
               enddo
          endif
-         call mpi_scatter(d_restart1,3*2*nres,mpi_real,
-     &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
+         call mpi_scatter(d_restart1,3*2*(nres+1),mpi_real,
+     &           r_d,3*2*(nres+1),mpi_real,king,CG_COMM,ierr)
          do i=1,2*nres
            do j=1,3
             dc(j,i)=r_d(j,i)
@@ -1874,227 +1880,3 @@ c     &                (d_restart1(j,i+2*nres*il),j=1,3)
         return
         end
 c------------------------------------------
-      subroutine returnbox
-      include 'DIMENSIONS'
-      include 'mpif.h'
-      include 'COMMON.CONTROL'
-      include 'COMMON.VAR'
-      include 'COMMON.MD'
-#ifndef LANG0
-      include 'COMMON.LANGEVIN'
-#else
-      include 'COMMON.LANGEVIN.lang0'
-#endif
-      include 'COMMON.CHAIN'
-      include 'COMMON.DERIV'
-      include 'COMMON.GEO'
-      include 'COMMON.LOCAL'
-      include 'COMMON.INTERACT'
-      include 'COMMON.IOUNITS'
-      include 'COMMON.NAMES'
-      include 'COMMON.TIME1'
-      include 'COMMON.REMD'
-      include 'COMMON.SETUP'
-      include 'COMMON.MUCA'
-      include 'COMMON.HAIRPIN'
-        j=1
-        chain_beg=1
-C        do i=1,nres
-C       write(*,*) 'initial', i,j,c(j,i)
-C        enddo
-        do i=1,nres-1
-         if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
-          chain_end=i
-          if (allareout.eq.1) then
-            ireturnval=int(c(j,i)/boxxsize)
-            if (c(j,i).le.0) ireturnval=ireturnval-1
-            do k=chain_beg,chain_end
-              c(j,k)=c(j,k)-ireturnval*boxxsize
-              c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
-            enddo
-           endif
-           chain_beg=i+1
-           allareout=1
-         else
-          if (int(c(j,i)/boxxsize).eq.0) allareout=0
-         endif
-        enddo
-         if (allareout.eq.1) then
-            ireturnval=int(c(j,i)/boxxsize)
-            if (c(j,i).le.0) ireturnval=ireturnval-1
-            do k=chain_beg,nres
-              c(j,k)=c(j,k)-ireturnval*boxxsize
-              c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
-            enddo
-          endif
-C NO JUMP 
-C        do i=1,nres
-C        write(*,*) 'befor no jump', i,j,c(j,i)
-C        enddo
-        nojumpval=0
-        do i=2,nres
-           if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
-             difference=abs(c(j,i-1)-c(j,i))
-C             print *,'diff', difference
-             if (difference.gt.boxxsize/2.0) then
-                if (c(j,i-1).gt.c(j,i)) then
-                  nojumpval=1
-                 else
-                   nojumpval=-1
-                 endif
-              else
-              nojumpval=0
-              endif
-              endif
-              c(j,i)=c(j,i)+nojumpval*boxxsize
-              c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
-         enddo
-       nojumpval=0
-        do i=2,nres
-           if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
-             difference=abs(c(j,i-1)-c(j,i))
-             if (difference.gt.boxxsize/2.0) then
-                if (c(j,i-1).gt.c(j,i)) then
-                  nojumpval=1
-                 else
-                   nojumpval=-1
-                 endif
-              else
-              nojumpval=0
-              endif
-             endif
-              c(j,i)=c(j,i)+nojumpval*boxxsize
-              c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
-         enddo
-
-C        do i=1,nres
-C        write(*,*) 'after no jump', i,j,c(j,i)
-C        enddo
-
-C NOW Y dimension
-        j=2
-        chain_beg=1
-        do i=1,nres-1
-         if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
-          chain_end=i
-          if (allareout.eq.1) then
-            ireturnval=int(c(j,i)/boxysize)
-            if (c(j,i).le.0) ireturnval=ireturnval-1
-            do k=chain_beg,chain_end
-              c(j,k)=c(j,k)-ireturnval*boxysize
-              c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
-            enddo
-           endif
-           chain_beg=i+1
-           allareout=1
-         else
-          if (int(c(j,i)/boxysize).eq.0) allareout=0
-         endif
-        enddo
-         if (allareout.eq.1) then
-            ireturnval=int(c(j,i)/boxysize)
-            if (c(j,i).le.0) ireturnval=ireturnval-1
-            do k=chain_beg,nres
-              c(j,k)=c(j,k)-ireturnval*boxysize
-              c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
-            enddo
-          endif
-        nojumpval=0
-        do i=2,nres
-           if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
-             difference=abs(c(j,i-1)-c(j,i))
-             if (difference.gt.boxysize/2.0) then
-                if (c(j,i-1).gt.c(j,i)) then
-                  nojumpval=1
-                 else
-                   nojumpval=-1
-                 endif
-              else
-              nojumpval=0
-              endif
-           endif
-              c(j,i)=c(j,i)+nojumpval*boxysize
-              c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize
-         enddo
-      nojumpval=0
-        do i=2,nres
-           if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
-             difference=abs(c(j,i-1)-c(j,i))
-             if (difference.gt.boxysize/2.0) then
-                if (c(j,i-1).gt.c(j,i)) then
-                  nojumpval=1
-                 else
-                   nojumpval=-1
-                 endif
-              else
-              nojumpval=0
-              endif
-            endif
-              c(j,i)=c(j,i)+nojumpval*boxysize
-              c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize
-         enddo
-
-        j=3
-        chain_beg=1
-        do i=1,nres-1
-         if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
-          chain_end=i
-          if (allareout.eq.1) then
-            ireturnval=int(c(j,i)/boxysize)
-            if (c(j,i).le.0) ireturnval=ireturnval-1
-            do k=chain_beg,chain_end
-              c(j,k)=c(j,k)-ireturnval*boxzsize
-              c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
-            enddo
-           endif
-           chain_beg=i+1
-           allareout=1
-         else
-          if (int(c(j,i)/boxzsize).eq.0) allareout=0
-         endif
-        enddo
-         if (allareout.eq.1) then
-            ireturnval=int(c(j,i)/boxzsize)
-            if (c(j,i).le.0) ireturnval=ireturnval-1
-            do k=chain_beg,nres
-              c(j,k)=c(j,k)-ireturnval*boxzsize
-              c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
-            enddo
-          endif
-        nojumpval=0
-        do i=2,nres
-           if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
-             difference=abs(c(j,i-1)-c(j,i))
-             if (difference.gt.(boxzsize/2.0)) then
-                if (c(j,i-1).gt.c(j,i)) then
-                  nojumpval=1
-                 else
-                   nojumpval=-1
-                 endif
-              else
-              nojumpval=0
-              endif
-            endif
-              c(j,i)=c(j,i)+nojumpval*boxzsize
-              c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize
-         enddo
-       nojumpval=0
-        do i=2,nres
-           if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
-             difference=abs(c(j,i-1)-c(j,i))
-             if (difference.gt.boxzsize/2.0) then
-                if (c(j,i-1).gt.c(j,i)) then
-                  nojumpval=1
-                 else
-                   nojumpval=-1
-                 endif
-              else
-              nojumpval=0
-              endif
-            endif
-              c(j,i)=c(j,i)+nojumpval*boxzsize
-              c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize
-         enddo
-
-        return
-        end