Adam's unres update

[unres.git] / source / unres / src-HCD-5D / MD_A-MTS.F

diff --git a/source/unres/src-HCD-5D/MD_A-MTS.F b/source/unres/src-HCD-5D/MD_A-MTS.F
index d82cf17..8cebc8f 100644 (file)
--- a/source/unres/src-HCD-5D/MD_A-MTS.F
+++ b/source/unres/src-HCD-5D/MD_A-MTS.F
@@ -65,6 +65,7 @@ c
       t_enegrad=0.0d0
       t_sdsetup=0.0d0
       write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
+      write (iout,*) "Restart frequency:",irest_freq
 #ifdef MPI
       tt0=MPI_Wtime()
 #else
@@ -261,13 +262,13 @@ C          call check_ecartint
              call cartout(totT)
           endif
         endif
-        if (rstcount.eq.1000.or.itime.eq.n_timestep) then
+        if (rstcount.eq.irest_freq.or.itime.eq.n_timestep) then
            open(irest2,file=rest2name,status='unknown')
            write(irest2,*) totT,EK,potE,totE,t_bath
-           do i=1,2*nres
+           do i=0,2*nres-1
             write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
            enddo
-           do i=1,2*nres
+           do i=0,2*nres-1
             write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
            enddo
           close(irest2)
@@ -290,7 +291,7 @@ C          call check_ecartint
      & '  End of MD calculation  '
 #ifdef TIMING_ENE
       write (iout,*) "time for etotal",t_etotal," elong",t_elong,
-     &  " eshort",t_eshort
+     &  " eshort",t_eshort," list",time_list
       write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,
      & " time_fricmatmult",time_fricmatmult," time_fsample ",
      & time_fsample
@@ -429,7 +430,7 @@ c Calculate energy and forces
         call zerograd
         call etotal(potEcomp)
 ! AL 4/17/17: Reduce the steps if NaNs occurred.
-        if (potEcomp(0).gt.0.99e20 .or. isnan(potEcomp(0)).gt.0) then
+        if (potEcomp(0).gt.0.99e20 .or. isnan(potEcomp(0))) then
           d_time=d_time/2
           cycle
         endif
@@ -1671,6 +1672,8 @@ c  Set up the initial conditions of a MD simulation
       integer iran_num
       double precision etot
       logical fail
+      integer i_start_models(0:nodes-1)
+      double precision potEcomp_long(0:Max_Ene)
       write (iout,*) "init_MD INDPDB",indpdb
       d_time0=d_time
 c      write(iout,*) "d_time", d_time
@@ -1758,10 +1761,10 @@ c        statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
          if (me.eq.king)
      &     inquire(file=mremd_rst_name,exist=file_exist)
 #ifdef MPI
-           write (*,*) me," Before broadcast: file_exist",file_exist
+c           write (*,*) me," Before broadcast: file_exist",file_exist
          call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
      &          IERR)
-         write (*,*) me," After broadcast: file_exist",file_exist
+c         write (*,*) me," After broadcast: file_exist",file_exist
 c        inquire(file=mremd_rst_name,exist=file_exist)
 #endif
         if(me.eq.king.or..not.out1file)
@@ -1809,7 +1812,7 @@ c      rest2name = prefix(:ilen(prefix))//'.rst'
       if(me.eq.king.or..not.out1file)then
        write (iout,*) "Initial velocities"
        do i=0,nres
-         write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
+         write (iout,'(i7,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
      &   (d_t(j,i+nres),j=1,3)
        enddo
       endif
@@ -1839,7 +1842,18 @@ c Removing the velocity of the center of mass
        call flush(iout)
       endif
       endif
-      write (iout,*) "init_MD before initial structure REST ",rest
+c      write (iout,*) "init_MD before initial structure REST ",rest
+      if(start_from_model .and. (me.eq.king .or. .not. out1file))
+     &    write(iout,*) 'START_FROM_MODELS is ON'
+      if(start_from_model .and. rest) then 
+        if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
+         write(iout,*) 
+     &   'START_FROM_MODELS is OFF because the run is restarted'
+         write(iout,*) 'Remove restart keyword from input'
+        endif
+      endif
+c      write (iout,*) "rest ",rest," start_from_model",start_from_model,
+c     &  " nmodel_start",nmodel_start," preminim",preminim
       if (.not.rest) then
   122   continue
         if (iranconf.ne.0) then
@@ -1920,9 +1934,12 @@ c 8/22/17 AL Loop to produce a low-energy random conformation
         else if (preminim) then
           if (start_from_model) then
             n_model_try=0
-            do while (fail .and. n_model_try.lt.constr_homology)
+            fail=.true.
+            list_model_try=0
+            do while (fail .and. n_model_try.lt.nmodel_start)
+              write (iout,*) "n_model_try",n_model_try
               do
-                i_model=iran_num(1,constr_homology)
+                i_model=iran_num(1,nmodel_start)
                 do k=1,n_model_try
                   if (i_model.eq.list_model_try(k)) exit
                 enddo
@@ -1930,7 +1947,9 @@ c 8/22/17 AL Loop to produce a low-energy random conformation
               enddo
               n_model_try=n_model_try+1
               list_model_try(n_model_try)=i_model
-              write (iout,*) 'starting from model ',i_model
+              if (me.eq.king .or. .not. out1file)
+     &         write (iout,*) 'Trying to start from model ',
+     &         pdbfiles_chomo(i_model)(:ilen(pdbfiles_chomo(i_model)))
               do i=1,2*nres
                 do j=1,3
                   c(j,i)=chomo(j,i,i_model)
@@ -1938,6 +1957,7 @@ c 8/22/17 AL Loop to produce a low-energy random conformation
               enddo
               call int_from_cart(.true.,.false.)
               call sc_loc_geom(.false.)
+              dc(:,0)=c(:,1)
               do i=1,nres-1
                 do j=1,3
                   dc(j,i)=c(j,i+1)-c(j,i)
@@ -1981,10 +2001,14 @@ c 8/22/17 AL Loop to produce a low-energy random conformation
               call etotal(energia(0))
 #endif
             enddo
-            if (n_model_try.gt.constr_homology) then
+            call MPI_Gather(i_model,1,MPI_INTEGER,i_start_models(0),
+     &        1,MPI_INTEGER,king,CG_COMM,IERROR)
+            if (n_model_try.gt.nmodel_start .and. 
+     &         (me.eq.king .or. out1file)) then
               write (iout,*) 
      &    "All models have irreparable overlaps. Trying randoms starts."
               iranconf=1
+              i_model=nmodel_start+1
               goto 122
             endif
           else
@@ -2027,6 +2051,17 @@ C 8/22/17 AL Minimize initial structure
 #endif
           endif
         endif
+        if (nmodel_start.gt.0 .and. me.eq.king) then
+          write (iout,'(a)') "Task  Starting model"
+          do i=0,nodes-1
+            if (i_start_models(i).gt.nmodel_start) then
+              write (iout,'(i4,2x,a)') i,"RANDOM STRUCTURE"
+            else
+              write(iout,'(i4,2x,a)')i,pdbfiles_chomo(i_start_models(i))
+     &          (:ilen(pdbfiles_chomo(i_start_models(i))))
+            endif
+          enddo
+        endif
       endif ! .not. rest
       call chainbuild_cart
       call kinetic(EK)
@@ -2034,6 +2069,7 @@ C 8/22/17 AL Minimize initial structure
         call verlet_bath
       endif
       kinetic_T=2.0d0/(dimen3*Rb)*EK
+      write (iout,*) "Initial kinetic energy",EK," kinetic T",kinetic_T
       if(me.eq.king.or..not.out1file)then
        call cartprint
        call intout
@@ -2046,6 +2082,7 @@ C 8/22/17 AL Minimize initial structure
       call zerograd
       call etotal(potEcomp)
       if (large) call enerprint(potEcomp)
+      if (RESPA) call etotal_long(potEcomp_long)
 #ifdef TIMING_ENE
 #ifdef MPI
       t_etotal=t_etotal+MPI_Wtime()-tt0
@@ -2058,7 +2095,7 @@ c      write (iout,*) "PotE-homology",potE-potEcomp(27)
       call cartgrad
       call lagrangian
       call max_accel
-      if (amax*d_time .gt. dvmax) then
+      if (amax*d_time .gt. dvmax .and. .not. respa) then
         d_time=d_time*dvmax/amax
         if(me.eq.king.or..not.out1file) write (iout,*) 
      &   "Time step reduced to",d_time,