Merge branch 'devel' into feature-ga

[unres.git] / source / unres / src_MD / MD_A-MTS.F

diff --git a/source/unres/src_MD/MD_A-MTS.F b/source/unres/src_MD/MD_A-MTS.F
index 042fcc0..acbffa9 100644 (file)
--- a/source/unres/src_MD/MD_A-MTS.F
+++ b/source/unres/src_MD/MD_A-MTS.F
@@ -592,13 +592,15 @@ c Backup the coordinates, velocities, and accelerations
         H=(HNose1-H0)*s_np
 cd        write (iout,'(a,10f)') "hhh",EK,s_np,potE,pi_np,H0
 cd     &   ,EK+potE+pi_np**2/(2*Q_np)+dimen3*0.001986d0*t_bath*log(s_np)
-        write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
+cd        write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
+          hhh=h
        endif
 
        if(tnh) then
         HNose1=Hnose_nh(EK,potE)
         H=HNose1-H0
-        write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
+        hhh=h
+cd        write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
        endif
 
        if (large) then
@@ -1104,14 +1106,16 @@ c Backup the coordinates, velocities, and accelerations
         H=(HNose1-H0)*s_np
 cd        write (iout,'(a,10f)') "hhh",EK,s_np,potE,pi_np,H0
 cd     &   ,EK+potE+pi_np**2/(2*Q_np)+dimen3*0.001986d0*t_bath*log(s_np)
-        write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
+cd        write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
+          hhh=h
 cd        write (iout,'(a,3f)') "EE2 NP S, pi",totT, s_np, pi_np
        endif
 
        if(tnh) then
         HNose1=Hnose_nh(EK,potE)
         H=HNose1-H0
-        write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
+cd        write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
+        hhh=h
        endif
 
 
@@ -1802,6 +1806,7 @@ 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                    
+       call flush(iout)
 c  Zeroing the total angular momentum of the system
        write(iout,*) "Calling the zero-angular 
      & momentum subroutine"
@@ -1819,6 +1824,7 @@ 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) 
+       call flush(iout)
       endif
       if (.not.rest) then              
          call chainbuild
@@ -1888,6 +1894,17 @@ c Removing the velocity of the center of mass
        write (iout,*) 'H0= ',H0
       endif
 
+      if (hmc.gt.0) then
+         hmc_acc=0
+         hmc_etot=potE+EK
+          if(me.eq.king.or..not.out1file)
+     &       write(iout,*) 'HMC',hmc_etot,potE,EK
+         do i=1,2*nres
+           do j=1,3
+            dc_hmc(j,i)=dc(j,i)
+           enddo
+         enddo
+      endif
 
       call cartgrad
       call lagrangian
@@ -2019,6 +2036,9 @@ C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
         t_enegrad=t_enegrad+tcpu()-tt0
 #endif
       endif
+
+
+
       return
       end
 c-----------------------------------------------------------
@@ -2060,6 +2080,7 @@ c          write (iout,*) "i",i," ii",ii," geigen",geigen(i),
 c     &      " d_t_work_new",d_t_work_new(ii)
         enddo
       enddo
+      call flush(iout)
 c diagnostics
 c      Ek1=0.0d0
 c      ii=0
@@ -2098,6 +2119,7 @@ c Transfer to the d_t vector
 c      do i=0,nres-1
 c        write (iout,*) "d_t",i,(d_t(j,i),j=1,3)
 c      enddo
+c      call flush(iout)
       do i=nnt,nct
         if (itype(i).ne.10) then
           do j=1,3
@@ -3379,3 +3401,61 @@ cd      write(iout,'(a,4f)') 'mmm',EK,potE,HNose1,pi_np
 
       return
       end
+
+      subroutine hmc_test(itime)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.CONTROL'
+      include 'COMMON.MD'
+      include 'COMMON.CHAIN'
+
+           hmc_acc=hmc_acc+1
+           delta=-(potE+EK-hmc_etot)/(Rb*t_bath)
+           if (delta .lt. -50.0d0) then
+                delta=0.0d0
+           else
+                delta=dexp(delta)
+           endif
+           xxx=ran_number(0.0d0,1.0d0)
+
+           if (me.eq.king .or. .not. out1file)
+     &       write(iout,'(a8,i5,6f10.4)') 
+     &        'HMC',itime,potE+EK,potE,EK,hmc_etot,delta,xxx
+
+           if (delta .le. xxx) then
+            do i=1,2*nres
+             do j=1,3
+              dc(j,i)=dc_hmc(j,i)
+             enddo
+            enddo
+            itime=itime-hmc
+            totT=totThmc
+           else
+            if (me.eq.king .or. .not. out1file)
+     &       write(iout,*) 'HMC accepting new'
+            totThmc=totT
+            do i=1,2*nres
+             do j=1,3
+              dc_hmc(j,i)=dc(j,i)
+             enddo
+            enddo
+           endif
+
+           call chainbuild_cart
+           call random_vel
+           do i=0,2*nres
+            do j=1,3
+              d_t_old(j,i)=d_t(j,i)
+            enddo
+           enddo
+           call kinetic(EK)
+           kinetic_T=2.0d0/(dimen3*Rb)*EK
+           call etotal(potEcomp)
+           potE=potEcomp(0)
+           hmc_etot=potE+EK
+           if (me.eq.king .or. .not. out1file)
+     &      write(iout,'(a8,i5,3f10.4)')'HMC new',itime,potE+EK,potE,EK
+
+
+      return
+      end