corrections
[unres.git] / source / unres / src-HCD-5D / MREMD.F
index 78a7404..f22e2f6 100644 (file)
@@ -57,11 +57,12 @@ cold      integer nup(0:maxprocs),ndown(0:maxprocs)
       integer rep2i(0:maxprocs),ireqi(maxprocs)
       integer icache_all(maxprocs)
       integer status(MPI_STATUS_SIZE),statusi(MPI_STATUS_SIZE,maxprocs)
-      logical synflag,end_of_run,file_exist /.false./,ovrtim
+      logical synflag,end_of_run,file_exist /.false./,ovrtim,first_pass
       double precision t_bath_temp,delta,ene_iex_iex,ene_i_i,ene_iex_i,
      & ene_i_iex,xxx,tmp,econstr_temp_iex,econstr_temp_i
       integer iran_num
       double precision ran_number
+      integer i_econstr/20/
 
 cdeb      imin_itime_old=0
       ntwx_cache=0
@@ -144,8 +145,8 @@ cold       endif
 cd       print '(i4,a6,100i4)',me,' ndown',(ndown(i),i=0,ndown(0))
 
        
-       write (*,*) "Processor",me," rest",rest,"
-     &   restart1fie",restart1file
+c       write (*,*) "Processor",me," rest",rest,"
+c     &   restart1fie",restart1file
        if(rest.and.restart1file) then 
            if (me.eq.king)
      &     inquire(file=mremd_rst_name,exist=file_exist)
@@ -175,6 +176,20 @@ cd           write (*,*) me," After broadcast: file_exist",file_exist
                enddo
               enddo
            endif 
+           stdfp=dsqrt(2*Rb*t_bath/d_time)
+           do i=1,ntyp
+             stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
+           enddo
+           if (lang.gt.0 .and. .not.surfarea) then
+             do i=nnt,nct-1
+               stdforcp(i)=stdfp*dsqrt(gamp)
+             enddo
+             do i=nnt,nct
+              if (itype(i).ne.ntyp1) stdforcsc(i)=stdfsc(iabs(itype(i)))
+     &                   *dsqrt(gamsc(iabs(itype(i))))
+              enddo
+           endif
+
        endif
 
        if(me.eq.king) then
@@ -237,6 +252,20 @@ cd           write (*,*) me," After broadcast: file_exist",file_exist
               write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",
      &                    (ndowna(i,il),i=1,ndowna(0,il))
              enddo
+
+             stdfp=dsqrt(2*Rb*t_bath/d_time)
+             do i=1,ntyp
+               stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
+             enddo
+             if (lang.gt.0 .and. .not.surfarea) then
+             do i=nnt,nct-1
+               stdforcp(i)=stdfp*dsqrt(gamp)
+             enddo
+             do i=nnt,nct
+              if (itype(i).ne.ntyp1) stdforcsc(i)=stdfsc(iabs(itype(i)))
+     &                   *dsqrt(gamsc(iabs(itype(i))))
+             enddo
+         endif
         ELSE IF (.not.(rest.and.file_exist)) THEN
          do il=1,remd_m(1)
           ifirst(il)=il
@@ -292,7 +321,7 @@ cd       print *,'ttt',me,remd_tlist,(remd_t(i),i=1,nrep)
           iset=i2set(me)
           if(me.eq.king.or..not.out1file) 
      &     write(iout,*) me,"iset=",iset,"t_bath=",t_bath
-          call flush(iout)
+c          call flush(iout)
        endif        
 c
        stdfp=dsqrt(2*Rb*t_bath/d_time)
@@ -395,6 +424,8 @@ c   Entering the MD loop
 #endif
       itime=0
       end_of_run=.false.
+      first_pass=.not.rest
+c      write (iout,*) "first_pass",first_pass
       do while(.not.end_of_run)
         itime=itime+1
         if(itime.eq.n_timestep.and.me.eq.king) end_of_run=.true.
@@ -694,7 +725,7 @@ ctime           call flush(iout)
                write(iout,*) "MIN ii_write=",ii_write
              endif
             endif
-           call flush(iout)
+c           call flush(iout)
         endif
         if (synflag) then
 c Update the time safety limiy
@@ -828,7 +859,7 @@ cd end
           time05=MPI_WTIME()
           if (me.eq.king .or. .not. out1file) then
             write(iout,*) 'REMD writing traj time=',time05-time04
-            call flush(iout)
+c            call flush(iout)
           endif
 
 
@@ -913,7 +944,10 @@ c               write (iout,*) "ene_i_i",remd_ene(0,i)
 c               call flush(iout)
 c               write (iout,*) "rescaling weights with temperature",
 c     &          remd_t_bath(iex)
-               if (real(ene_i_i).ne.real(remd_ene(0,i))) then
+c               write (iout,*) "first_pass",first_pass
+               if (.not.first_pass.and.
+     &           real(ene_i_i).ne.real(remd_ene(0,i))) 
+     &         then
                 write (iout,*) "ERROR: inconsistent energies:",i,
      &            ene_i_i,remd_ene(0,i)
                endif
@@ -931,7 +965,8 @@ c               write (iout,*) "0,iex",remd_t_bath(iex)
 c               call enerprint(remd_ene(0,iex))
 
                call sum_energy(remd_ene(0,iex),.false.)
-               if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
+               if (.not.first_pass.and.
+     &           real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
                 write (iout,*) "ERROR: inconsistent energies:",iex,
      &            ene_iex_iex,remd_ene(0,iex)
                endif
@@ -1036,6 +1071,7 @@ c              call flush(iout)
             endif
            enddo
            enddo
+           first_pass=.false.
 cd           write (iout,*) "exchange completed"
 cd           call flush(iout) 
         ELSE
@@ -1062,12 +1098,12 @@ cd            write(iout,*) "i_dir=",i_dir
                iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
 c 9/1/17 AL: Correction; otherwise the restraint energies are mis-assigned 
 c       on failed replica exchange attempt
-               econstr_temp_i=remd_ene(20,i)
-               econstr_temp_iex=remd_ene(20,iex)
+               econstr_temp_i=remd_ene(i_econstr,i)
+               econstr_temp_iex=remd_ene(i_econstr,iex)
 c 9/11/17 AL: Adaptive sampling (temperature dependent restraints potentials)
                if (adaptive) then
-                 remd_ene(20,i)=remd_ene(n_ene+5,i)
-                 remd_ene(20,iex)=remd_ene(n_ene+6,iex)
+                 remd_ene(i_econstr,i)=remd_ene(n_ene+5,i)
+                 remd_ene(i_econstr,iex)=remd_ene(n_ene+6,iex)
                endif
             elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then
 
@@ -1076,10 +1112,10 @@ c 9/11/17 AL: Adaptive sampling (temperature dependent restraints potentials)
                i_iset1=i_iset+1
                i_mset1=iran_num(1,mset(i_iset1))
                iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
-               econstr_temp_i=remd_ene(20,i)
-               econstr_temp_iex=remd_ene(20,iex)
-               remd_ene(20,i)=remd_ene(n_ene+3,i)
-               remd_ene(20,iex)=remd_ene(n_ene+4,iex)
+               econstr_temp_i=remd_ene(i_econstr,i)
+               econstr_temp_iex=remd_ene(i_econstr,iex)
+               remd_ene(i_econstr,i)=remd_ene(n_ene+3,i)
+               remd_ene(i_econstr,iex)=remd_ene(n_ene+4,iex)
 
             elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then
 
@@ -1088,14 +1124,14 @@ c 9/11/17 AL: Adaptive sampling (temperature dependent restraints potentials)
                i_iset1=i_iset+1
                i_mset1=iran_num(1,mset(i_iset1))
                iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
-               econstr_temp_i=remd_ene(20,i)
-               econstr_temp_iex=remd_ene(20,iex)
+               econstr_temp_i=remd_ene(i_econstr,i)
+               econstr_temp_iex=remd_ene(i_econstr,iex)
                if (adaptive) then
-                 remd_ene(20,i)=remd_ene(n_ene+7,i)
-                 remd_ene(20,iex)=remd_ene(n_ene+8,iex)
+                 remd_ene(i_econstr,i)=remd_ene(n_ene+7,i)
+                 remd_ene(i_econstr,iex)=remd_ene(n_ene+8,iex)
                else
-                 remd_ene(20,i)=remd_ene(n_ene+3,i)
-                 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
+                 remd_ene(i_econstr,i)=remd_ene(n_ene+3,i)
+                 remd_ene(i_econstr,iex)=remd_ene(n_ene+4,iex)
                endif
             else
                goto 444 
@@ -1195,8 +1231,8 @@ cd              write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
               else
                remd_ene(0,iex)=ene_iex_iex
                remd_ene(0,i)=ene_i_i
-               remd_ene(20,iex)=econstr_temp_iex
-               remd_ene(20,i)=econstr_temp_i
+               remd_ene(i_econstr,iex)=econstr_temp_iex
+               remd_ene(i_econstr,i)=econstr_temp_i
               endif
 
 cd      do il=1,nset
@@ -1232,7 +1268,7 @@ c-------------------------------------
               enddo
              endif
 
-             call flush(iout)
+c             call flush(iout)
 
 cd              write (iout,'(a6,100i4)') "ifirst",
 cd     &                    (ifirst(i),i=1,remd_m(1))
@@ -1278,12 +1314,22 @@ co     &    " rescaling weights with temperature",t_bath
          do i=1,ntyp
            stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
          enddo
-
-cde         write(iout,*) 'REMD after',me,t_bath
+c Compute the standard deviations of stochastic forces for Langevin dynamics
+c if the friction coefficients do not depend on surface area
+         if (lang.gt.0 .and. .not.surfarea) then
+           do i=nnt,nct-1
+             stdforcp(i)=stdfp*dsqrt(gamp)
+           enddo
+           do i=nnt,nct
+             if (itype(i).ne.ntyp1) stdforcsc(i)=stdfsc(iabs(itype(i)))
+     &                   *dsqrt(gamsc(iabs(itype(i))))
+           enddo
+         endif
+cde           write(iout,*) 'REMD after',me,t_bath
            time08=MPI_WTIME()
            if (me.eq.king .or. .not. out1file) then
             write(iout,*) 'REMD exchange time=',time08-time00
-            call flush(iout)
+c            call flush(iout)
            endif
         endif
       enddo
@@ -1344,7 +1390,7 @@ c-----------------------------------------------------------------------
       include 'COMMON.SBRIDGE'
       include 'COMMON.INTERACT'
                
-      real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
+      real d_restart1(3,2*maxres*maxprocs),r_d(3,0:2*maxres-1),
      &     d_restart2(3,2*maxres*maxprocs)
       real t5_restart1(5)
       integer iret,itmp
@@ -1364,7 +1410,7 @@ c-----------------------------------------------------------------------
      &      t_restart1,5,mpi_real,king,CG_COMM,ierr)
 
 
-       do i=1,2*nres
+       do i=0,2*nres-1
          do j=1,3
            r_d(j,i)=d_t(j,i)
          enddo
@@ -1374,7 +1420,7 @@ c-----------------------------------------------------------------------
      &           CG_COMM,ierr)
 
 
-       do i=1,2*nres
+       do i=0,2*nres-1
          do j=1,3
            r_d(j,i)=dc(j,i)
          enddo
@@ -1755,14 +1801,14 @@ 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*maxprocs),r_d(3,0:2*maxres-1),
      &                 t5_restart1(5)
       integer*2 i_index
      &            (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
       common /przechowalnia/ d_restart1
       integer i,j,il,il1,ixdrf,iret,itmp
       integer ierr
-      write (*,*) "Processor",me," called read1restart"
+c      write (*,*) "Processor",me," called read1restart"
 
          if(me.eq.king)then
               open(irest2,file=mremd_rst_name,status='unknown')
@@ -1849,7 +1895,7 @@ c     &                (d_restart1(j,i+2*nres*il),j=1,3)
          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
+         do i=0,2*nres-1
            do j=1,3
             d_t(j,i)=r_d(j,i)
            enddo
@@ -1871,7 +1917,7 @@ c     &                (d_restart1(j,i+2*nres*il),j=1,3)
          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
+         do i=0,2*nres-1
            do j=1,3
             dc(j,i)=r_d(j,i)
            enddo
@@ -1952,7 +1998,7 @@ c     &           CG_COMM,ierr)
       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*maxprocs),r_d(3,0:2*maxres-1),
      &                 t5_restart1(5)
       common /przechowalnia/ d_restart1
       integer i,j,il,itmp
@@ -1988,7 +2034,7 @@ c     &           CG_COMM,ierr)
          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
+         do i=0,2*nres-1
            do j=1,3
             d_t(j,i)=r_d(j,i)
            enddo
@@ -2003,7 +2049,7 @@ c     &           CG_COMM,ierr)
          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
+         do i=0,2*nres-1
            do j=1,3
             dc(j,i)=r_d(j,i)
            enddo