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
if(tnp .or. tnp1) then
#ifndef G77
- write (iout,'(a3,7f)') "TTT",EK,s_np,potE,pi_np,Csplit,
+ write (iout,'(a3,7f20.10)') "TTT",EK,s_np,potE,pi_np,Csplit,
& E_long,energia_short(0)
#else
write (iout,'(a3,7f20.10)') "TTT",EK,s_np,potE,pi_np,Csplit,
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
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"
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
call chainbuild_cart
call kinetic(EK)
if (tbf) then
- call verlet_bath(EK)
+ call verlet_bath
endif
kinetic_T=2.0d0/(dimen3*Rb)*EK
if(me.eq.king.or..not.out1file)then
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
t_enegrad=t_enegrad+tcpu()-tt0
#endif
endif
+
+
+
return
end
c-----------------------------------------------------------
c & " d_t_work_new",d_t_work_new(ii)
enddo
enddo
+ call flush(iout)
c diagnostics
c Ek1=0.0d0
c ii=0
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
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