critical change for outer field potential -debug off
[unres.git] / source / unres / src_MD-M / MREMD.F
index b7d2d1f..fd83ca0 100644 (file)
@@ -304,8 +304,11 @@ c------------------------------------------------
 #endif
 c Determine the inverse of the inertia matrix.
       call setup_MD_matrices
+      write(iout,*), "TU", dc(1,0)
 c Initialize MD
       call init_MD
+      write(iout,*), "TU2", dc(1,0)
+
       if (rest) then  
        if (me.eq.king .or. .not. out1file)
      &  write (iout,'(a60,f10.5)') "REMD restart Temperature:",t_bath
@@ -525,7 +528,9 @@ c Variable time step algorithm.
               ugamma_cache(i,ntwx_cache)=ugamma(i)
               uscdiff_cache(i,ntwx_cache)=uscdiff(i)
             enddo
-
+C            print *,'przed returnbox'
+            call returnbox
+C            call enerprint(remd_ene(0,i))
             do i=1,nres*2
              do j=1,3
               c_cache(j,i,ntwx_cache)=c(j,i)
@@ -570,10 +575,10 @@ c Variable time step algorithm.
            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
@@ -831,7 +836,7 @@ c     &          remd_t_bath(iex)
                call rescale_weights(remd_t_bath(iex))
 
 c               write (iout,*) "0,i",remd_t_bath(iex)
-c               call enerprint(remd_ene(0,i))
+               call enerprint(remd_ene(0,i))
 
                call sum_energy(remd_ene(0,i),.false.)
 c               write (iout,*) "ene_i_iex",remd_ene(0,i)
@@ -959,7 +964,7 @@ cd            write(iout,*) "########",ii
 
 cd            write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
 
-            i_dir=iran_num(1,3)
+             i_dir=iran_num(1,3)
 cd            write(iout,*) "i_dir=",i_dir
 
             if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then            
@@ -1170,7 +1175,7 @@ co     &    " rescaling weights with temperature",t_bath
          stdfp=dsqrt(2*Rb*t_bath/d_time)
          do i=1,ntyp
            stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
-         enddo 
+         enddo
 
 cde         write(iout,*) 'REMD after',me,t_bath
            time08=MPI_WTIME()
@@ -1230,8 +1235,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,0:2*(maxres)*maxprocs),r_d(3,0:2*maxres),
+     &     d_restart2(3,0:2*(maxres)*maxprocs)
       real t5_restart1(5)
       integer iret,itmp
       integer*2 i_index
@@ -1248,23 +1253,26 @@ 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
+C       write (*,*) dc(j,0),"TU3"
+       do j=1,3
+       dc(j,0)=c(j,1)
+       enddo
+       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
@@ -1293,16 +1301,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
@@ -1353,16 +1363,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
@@ -1629,7 +1641,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,0:2*(maxres)*maxprocs),r_d(3,0:2*(maxres)),
      &                 t5_restart1(5)
       integer*2 i_index
      &            (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
@@ -1705,45 +1717,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
@@ -1814,7 +1828,7 @@ c     &                (d_restart1(j,i+2*nres*il),j=1,3)
       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,0:2*maxres*maxprocs),r_d(3,0:2*maxres),
      &                 t5_restart1(5)
       common /przechowalnia/ d_restart1
          if(me.eq.king)then
@@ -1839,31 +1853,31 @@ 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 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
                 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)
-         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
@@ -1871,4 +1885,4 @@ c     &                (d_restart1(j,i+2*nres*il),j=1,3)
         if(me.eq.king) close(irest2)
         return
         end
-          
+c------------------------------------------