update HCD-5D

[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 08852ba..4498b4a 100644 (file)
--- a/source/unres/src-HCD-5D/MD_A-MTS.F
+++ b/source/unres/src-HCD-5D/MD_A-MTS.F
@@ -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)
@@ -1811,10 +1812,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,8 +1831,14 @@ 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
       endif
-      call flush(iout)
       write (iout,*) "init_MD before initial structure REST ",rest
       if (.not.rest) then
   122   continue
@@ -2050,6 +2058,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:",
@@ -2351,6 +2360,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 +2375,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))