X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD%2FMD_A-MTS.F;h=29e4fcb1f2761e5cf611e91830d1515e964a8ebc;hb=c51f4eb43d54ac93f92168b225453afc56383ba2;hp=042fcc0b97aa4d134fe71c27114eb1d609a2bfef;hpb=a45ce81aae2c857c30b821e39e45a7a08356de4d;p=unres.git

diff --git a/source/unres/src_MD/MD_A-MTS.F b/source/unres/src_MD/MD_A-MTS.F
index 042fcc0..29e4fcb 100644
--- 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
 
 
@@ -1888,6 +1892,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 +2034,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-----------------------------------------------------------
@@ -3379,3 +3397,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