Adam's changes
[unres.git] / source / unres / src-HCD-5D / MD_A-MTS.F
index ca52aaa..fcef69e 100644 (file)
@@ -166,14 +166,14 @@ c   Entering the MD loop
      &      .and. mod(itime,count_reset_vel).eq.0) then
           call random_vel
           write(iout,'(a,f20.2)') 
-     &     "Velocities reset to random values, time",totT      
+     &     "Velocities reset to random values, time",totT
           do i=0,2*nres
             do j=1,3
               d_t_old(j,i)=d_t(j,i)
             enddo
           enddo
         endif
-               if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
+        if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
           call inertia_tensor  
           call vcm_vel(vcm)
           do j=1,3
@@ -182,7 +182,7 @@ c   Entering the MD loop
           call kinetic(EK)
           kinetic_T=2.0d0/(dimen3*Rb)*EK
           scalfac=dsqrt(T_bath/kinetic_T)
-          write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT      
+          write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
           do i=0,2*nres
             do j=1,3
               d_t_old(j,i)=scalfac*d_t(j,i)
@@ -209,6 +209,7 @@ c Variable time step algorithm.
           stop
 #endif
         endif
+        itime_mat=itime
         if (ntwe.ne.0) then
          if (mod(itime,ntwe).eq.0) then
            call statout(itime)
@@ -263,10 +264,10 @@ C          call check_ecartint
         if (rstcount.eq.1000.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)
@@ -1256,7 +1257,6 @@ c      write (iout,*) "friction accelerations"
       call fivediaginv_mult(dimen,fric_work, d_af_work)
 c      write (iout,*) "stochastic acceleratios"
       call fivediaginv_mult(dimen,stochforcvec, d_as_work)
-c      write (iout,*) "Leaving sddir_precalc"
 #else
       call ginv_mult(fric_work, d_af_work)
       call ginv_mult(stochforcvec, d_as_work)
@@ -1266,6 +1266,7 @@ c      write (iout,*) "Leaving sddir_precalc"
       write (iout,'(3f10.5)') (d_af_work(i),i=1,dimen3)
       write (iout,*) "d_as_work"
       write (iout,'(3f10.5)') (d_as_work(i),i=1,dimen3)
+      write (iout,*) "Leaving sddir_precalc"
 #endif
       return
       end
@@ -1670,6 +1671,7 @@ 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)
       write (iout,*) "init_MD INDPDB",indpdb
       d_time0=d_time
 c      write(iout,*) "d_time", d_time
@@ -1757,10 +1759,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)
@@ -1811,10 +1813,11 @@ c      rest2name = prefix(:ilen(prefix))//'.rst'
          write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
      &   (d_t(j,i+nres),j=1,3)
        enddo
+      endif
+      if (lang.eq.0) then
 c  Zeroing the total angular momentum of the system
        write(iout,*) "Calling the zero-angular 
      & momentum subroutine"
-      endif
       call inertia_tensor  
 c  Getting the potential energy and forces and velocities and accelerations
       call vcm_vel(vcm)
@@ -1829,9 +1832,26 @@ c Removing the velocity of the center of mass
       if(me.eq.king.or..not.out1file)then
        write (iout,*) "vcm right after adjustment:"
        write (iout,*) (vcm(j),j=1,3) 
+       write (iout,*) "Initial velocities after adjustment"
+       do i=0,nres
+         write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
+     &   (d_t(j,i+nres),j=1,3)
+       enddo
+       call flush(iout)
       endif
-      call flush(iout)
-      write (iout,*) "init_MD before initial structure REST ",rest
+      endif
+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
@@ -1912,9 +1932,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
@@ -1922,7 +1945,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)
@@ -1957,12 +1982,12 @@ c 8/22/17 AL Loop to produce a low-energy random conformation
                   cycle
                 endif
               endif
-#ifdef SEARCHSC
               if (me.eq.king.or..not.out1file) then
               write (iout,*) "Energies after removing overlaps"
               call etotal(energia(0))
               call enerprint(energia(0))
               endif
+#ifdef SEARCHSC
 ! Search for better SC rotamers if requested
               if (searchsc) then
                 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
@@ -1973,12 +1998,31 @@ 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
+! Remove SC overlaps if requested
+              if (overlapsc) then
+                write (iout,*) 'Calling OVERLAP_SC'
+                call overlap_sc(fail)
+                if (fail) then 
+                  write (iout,*) 
+     &            "Failed to remove overlap"
+                endif
+              endif
+              if (me.eq.king.or..not.out1file) then
+              write (iout,*) "Energies after removing overlaps"
+              call etotal(energia(0))
+              call enerprint(energia(0))
+              endif
           endif
 C 8/22/17 AL Minimize initial structure
           if (dccart) then
@@ -2004,6 +2048,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)
@@ -2011,6 +2066,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
@@ -2050,6 +2106,7 @@ c       write(iout,*) (potEcomp(i),i=0,n_ene)
 c      write (iout,*) "PotE-homology",potE
       totE=EK+potE
       itime=0
+      itime_mat=itime
       if (ntwe.ne.0) call statout(itime)
       if(me.eq.king.or..not.out1file)
      &  write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
@@ -2187,7 +2244,7 @@ c-----------------------------------------------------------
       double precision xv,sigv,lowb,highb,vec_afm(3),Ek1,Ek2,Ek3,aux
       integer i,ii,j,k,l,ind
       double precision anorm_distr
-      logical lprn /.true./
+      logical lprn /.false./
 #ifdef FIVEDIAG
       integer ichain,n,innt,inct,ibeg,ierr
       double precision work(8*maxres6)
@@ -2351,6 +2408,9 @@ c        write (iout,*) "ichain",ichain," innt",innt," inct",inct
 !
 #define WLOS
 #ifdef WLOS
+      if (nnt.eq.1) then
+        d_t(:,0)=d_t(:,1)
+      endif
       do i=1,nres
         if (itype(i).eq.10 .or. itype(i).eq.ntyp1) then
           do j=1,3
@@ -2363,11 +2423,17 @@ c        write (iout,*) "ichain",ichain," innt",innt," inct",inct
           enddo
         end if
       enddo
+      d_t(:,nres)=0.0d0
       d_t(:,nct)=0.0d0
+      d_t(:,2*nres)=0.0d0
+      if (nnt.gt.1) then
+        d_t(:,0)=d_t(:,1)
+        d_t(:,1)=0.0d0
+      endif
 c      d_a(:,0)=d_a(:,1)
 c      d_a(:,1)=0.0d0
 c      write (iout,*) "Shifting accelerations"
-      do ichain=1,nchain
+      do ichain=2,nchain
 c        write (iout,*) "ichain",chain_border1(1,ichain)-1,
 c     &     chain_border1(1,ichain)
         d_t(:,chain_border1(1,ichain)-1)=d_t(:,chain_border1(1,ichain))
@@ -2797,8 +2863,8 @@ c      enddo
         dc(j,0)=dc_work(j)
         d_t(j,0)=d_t_work(j)
       enddo
-      ind=3    
-      do i=nnt,nct-1   
+      ind=3
+      do i=nnt,nct-1
         do j=1,3
           dc(j,i)=dc_work(ind+j)
           d_t(j,i)=d_t_work(ind+j)