--- /dev/null
+ subroutine MD
+c------------------------------------------------
+c The driver for molecular dynamics subroutines
+c------------------------------------------------
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+#ifdef MPI
+ include "mpif.h"
+ integer IERROR,ERRCODE
+#endif
+ include 'COMMON.SETUP'
+ include 'COMMON.CONTROL'
+ include 'COMMON.VAR'
+ include 'COMMON.MD'
+#ifndef LANG0
+ include 'COMMON.LANGEVIN'
+#else
+ include 'COMMON.LANGEVIN.lang0'
+#endif
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.NAMES'
+ include 'COMMON.TIME1'
+ double precision cm(3),L(3),vcm(3)
+#ifdef VOUT
+ double precision v_work(maxres6),v_transf(maxres6)
+#endif
+ integer ilen,rstcount
+ external ilen
+ character*50 tytul
+ common /gucio/ cm
+ integer itime
+c
+#ifdef MPI
+ if (ilen(tmpdir).gt.0)
+ & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"
+ & //liczba(:ilen(liczba))//'.rst')
+#else
+ if (ilen(tmpdir).gt.0)
+ & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_"//'.rst')
+#endif
+ t_MDsetup=0.0d0
+ t_langsetup=0.0d0
+ t_MD=0.0d0
+ t_enegrad=0.0d0
+ t_sdsetup=0.0d0
+ write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
+#ifdef MPI
+ tt0=MPI_Wtime()
+#else
+ tt0 = tcpu()
+#endif
+c Determine the inverse of the inertia matrix.
+ call setup_MD_matrices
+c Initialize MD
+ call init_MD
+#ifdef MPI
+ t_MDsetup = MPI_Wtime()-tt0
+#else
+ t_MDsetup = tcpu()-tt0
+#endif
+ rstcount=0
+c Entering the MD loop
+#ifdef MPI
+ tt0 = MPI_Wtime()
+#else
+ tt0 = tcpu()
+#endif
+ if (lang.eq.2 .or. lang.eq.3) then
+#ifndef LANG0
+ call setup_fricmat
+ if (lang.eq.2) then
+ call sd_verlet_p_setup
+ else
+ call sd_verlet_ciccotti_setup
+ endif
+ do i=1,dimen3
+ do j=1,dimen3
+ pfric0_mat(i,j,0)=pfric_mat(i,j)
+ afric0_mat(i,j,0)=afric_mat(i,j)
+ vfric0_mat(i,j,0)=vfric_mat(i,j)
+ prand0_mat(i,j,0)=prand_mat(i,j)
+ vrand0_mat1(i,j,0)=vrand_mat1(i,j)
+ vrand0_mat2(i,j,0)=vrand_mat2(i,j)
+ enddo
+ enddo
+ flag_stoch(0)=.true.
+ do i=1,maxflag_stoch
+ flag_stoch(i)=.false.
+ enddo
+#else
+ write (iout,*)
+ & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
+#ifdef MPI
+ call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
+#endif
+ stop
+#endif
+ else if (lang.eq.1 .or. lang.eq.4) then
+ call setup_fricmat
+ endif
+#ifdef MPI
+ t_langsetup=MPI_Wtime()-tt0
+ tt0=MPI_Wtime()
+#else
+ t_langsetup=tcpu()-tt0
+ tt0=tcpu()
+#endif
+ do itime=1,n_timestep
+ rstcount=rstcount+1
+ if (lang.gt.0 .and. surfarea .and.
+ & mod(itime,reset_fricmat).eq.0) then
+ if (lang.eq.2 .or. lang.eq.3) then
+#ifndef LANG0
+ call setup_fricmat
+ if (lang.eq.2) then
+ call sd_verlet_p_setup
+ else
+ call sd_verlet_ciccotti_setup
+ endif
+ do i=1,dimen3
+ do j=1,dimen3
+ pfric0_mat(i,j,0)=pfric_mat(i,j)
+ afric0_mat(i,j,0)=afric_mat(i,j)
+ vfric0_mat(i,j,0)=vfric_mat(i,j)
+ prand0_mat(i,j,0)=prand_mat(i,j)
+ vrand0_mat1(i,j,0)=vrand_mat1(i,j)
+ vrand0_mat2(i,j,0)=vrand_mat2(i,j)
+ enddo
+ enddo
+ flag_stoch(0)=.true.
+ do i=1,maxflag_stoch
+ flag_stoch(i)=.false.
+ enddo
+#endif
+ else if (lang.eq.1 .or. lang.eq.4) then
+ call setup_fricmat
+ endif
+ write (iout,'(a,i10)')
+ & "Friction matrix reset based on surface area, itime",itime
+ endif
+ if (reset_vel .and. tbf .and. lang.eq.0
+ & .and. mod(itime,count_reset_vel).eq.0) then
+ call random_vel
+ write(iout,'(a,f20.2)')
+ & "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
+ call inertia_tensor
+ call vcm_vel(vcm)
+ do j=1,3
+ d_t(j,0)=d_t(j,0)-vcm(j)
+ enddo
+ 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
+ do i=0,2*nres
+ do j=1,3
+ d_t_old(j,i)=scalfac*d_t(j,i)
+ enddo
+ enddo
+ endif
+ if (lang.ne.4) then
+ if (RESPA) then
+c Time-reversible RESPA algorithm
+c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
+ call RESPA_step(itime)
+ else
+c Variable time step algorithm.
+ call velverlet_step(itime)
+ endif
+ else
+#ifdef BROWN
+ call brown_step(itime)
+#else
+ print *,"Brown dynamics not here!"
+#ifdef MPI
+ call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
+#endif
+ stop
+#endif
+ endif
+ if (ntwe.ne.0) then
+ if (mod(itime,ntwe).eq.0) call statout(itime)
+#ifdef VOUT
+ do j=1,3
+ v_work(j)=d_t(j,0)
+ enddo
+ ind=3
+ do i=nnt,nct-1
+ do j=1,3
+ ind=ind+1
+ v_work(ind)=d_t(j,i)
+ enddo
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ do j=1,3
+ ind=ind+1
+ v_work(ind)=d_t(j,i+nres)
+ enddo
+ endif
+ enddo
+
+ write (66,'(80f10.5)')
+ & ((d_t(j,i),j=1,3),i=0,nres-1),((d_t(j,i+nres),j=1,3),i=1,nres)
+ do i=1,ind
+ v_transf(i)=0.0d0
+ do j=1,ind
+ v_transf(i)=v_transf(i)+gvec(j,i)*v_work(j)
+ enddo
+ v_transf(i)= v_transf(i)*dsqrt(geigen(i))
+ enddo
+ write (67,'(80f10.5)') (v_transf(i),i=1,ind)
+#endif
+ endif
+ if (mod(itime,ntwx).eq.0) then
+ write (tytul,'("time",f8.2)') totT
+ if(mdpdb) then
+ call pdbout(potE,tytul,ipdb)
+ else
+ call cartout(totT)
+ endif
+ endif
+ 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
+ write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
+ enddo
+ do i=1,2*nres
+ write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
+ enddo
+ close(irest2)
+ rstcount=0
+ endif
+ enddo
+#ifdef MPI
+ t_MD=MPI_Wtime()-tt0
+#else
+ t_MD=tcpu()-tt0
+#endif
+ write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
+ & ' Timing ',
+ & 'MD calculations setup:',t_MDsetup,
+ & 'Energy & gradient evaluation:',t_enegrad,
+ & 'Stochastic MD setup:',t_langsetup,
+ & 'Stochastic MD step setup:',t_sdsetup,
+ & 'MD steps:',t_MD
+ write (iout,'(/28(1h=),a25,27(1h=))')
+ & ' End of MD calculation '
+#ifdef TIMING_ENE
+ write (iout,*) "time for etotal",t_etotal," elong",t_elong,
+ & " eshort",t_eshort
+ write (iout,*) "time_fric",time_fric," time_stoch",time_stoch,
+ & " time_fricmatmult",time_fricmatmult," time_fsample ",
+ & time_fsample
+#endif
+ return
+ end
+c-------------------------------------------------------------------------------
+ subroutine velverlet_step(itime)
+c-------------------------------------------------------------------------------
+c Perform a single velocity Verlet step; the time step can be rescaled if
+c increments in accelerations exceed the threshold
+c-------------------------------------------------------------------------------
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+#ifdef MPI
+ include 'mpif.h'
+ integer ierror,ierrcode
+#endif
+ include 'COMMON.SETUP'
+ include 'COMMON.CONTROL'
+ include 'COMMON.VAR'
+ include 'COMMON.MD'
+#ifndef LANG0
+ include 'COMMON.LANGEVIN'
+#else
+ include 'COMMON.LANGEVIN.lang0'
+#endif
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.NAMES'
+ include 'COMMON.TIME1'
+ include 'COMMON.MUCA'
+ double precision vcm(3),incr(3)
+ double precision cm(3),L(3)
+ integer ilen,count,rstcount
+ external ilen
+ character*50 tytul
+ integer maxcount_scale /20/
+ common /gucio/ cm
+ double precision stochforcvec(MAXRES6)
+ common /stochcalc/ stochforcvec
+ integer itime
+ logical scale
+ double precision HNose1,HNose,HNose_nh,H,vtnp(maxres6)
+ double precision vtnp_(maxres6),vtnp_a(maxres6)
+c
+ scale=.true.
+ icount_scale=0
+ if (lang.eq.1) then
+ call sddir_precalc
+ else if (lang.eq.2 .or. lang.eq.3) then
+#ifndef LANG0
+ call stochastic_force(stochforcvec)
+#else
+ write (iout,*)
+ & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
+#ifdef MPI
+ call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
+#endif
+ stop
+#endif
+ endif
+ itime_scal=0
+ do while (scale)
+ icount_scale=icount_scale+1
+ if (icount_scale.gt.maxcount_scale) then
+ write (iout,*)
+ & "ERROR: too many attempts at scaling down the time step. ",
+ & "amax=",amax,"epdrift=",epdrift,
+ & "damax=",damax,"edriftmax=",edriftmax,
+ & "d_time=",d_time
+ call flush(iout)
+#ifdef MPI
+ call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
+#endif
+ stop
+ endif
+c First step of the velocity Verlet algorithm
+ if (lang.eq.2) then
+#ifndef LANG0
+ call sd_verlet1
+#endif
+ else if (lang.eq.3) then
+#ifndef LANG0
+ call sd_verlet1_ciccotti
+#endif
+ else if (lang.eq.1) then
+ call sddir_verlet1
+ else if (tnp1) then
+ call tnp1_step1
+ else if (tnp) then
+ call tnp_step1
+ else
+ if (tnh) then
+ call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8)
+ do i=0,2*nres
+ do j=1,3
+ d_t_old(j,i)=d_t_old(j,i)*scale_nh
+ enddo
+ enddo
+ endif
+ call verlet1
+ endif
+c Build the chain from the newly calculated coordinates
+ call chainbuild_cart
+ if (rattle) call rattle1
+ if (ntwe.ne.0) then
+ if (large.and. mod(itime,ntwe).eq.0) then
+ write (iout,*) "Cartesian and internal coordinates: step 1"
+ call cartprint
+ call intout
+ write (iout,*) "dC"
+ do i=0,nres
+ write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
+ & (dc(j,i+nres),j=1,3)
+ enddo
+ write (iout,*) "Accelerations"
+ do i=0,nres
+ write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
+ & (d_a(j,i+nres),j=1,3)
+ enddo
+ write (iout,*) "Velocities, step 1"
+ 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
+ endif
+ endif
+#ifdef MPI
+ tt0 = MPI_Wtime()
+#else
+ tt0 = tcpu()
+#endif
+c Calculate energy and forces
+ call zerograd
+ call etotal(potEcomp)
+#ifdef TIMING_ENE
+#ifdef MPI
+ t_etotal=t_etotal+MPI_Wtime()-tt0
+#else
+ t_etotal=t_etotal+tcpu()-tt0
+#endif
+#endif
+ E_old=potE
+ potE=potEcomp(0)-potEcomp(20)
+ call cartgrad
+c Get the new accelerations
+ call lagrangian
+#ifdef MPI
+ t_enegrad=t_enegrad+MPI_Wtime()-tt0
+#else
+ t_enegrad=t_enegrad+tcpu()-tt0
+#endif
+c Determine maximum acceleration and scale down the timestep if needed
+ call max_accel
+ amax=amax/(itime_scal+1)**2
+ call predict_edrift(epdrift)
+ if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
+c Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
+ scale=.true.
+ ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax))
+ & /dlog(2.0d0)+1
+ itime_scal=itime_scal+ifac_time
+c fac_time=dmin1(damax/amax,0.5d0)
+ fac_time=0.5d0**ifac_time
+ d_time=d_time*fac_time
+ if (lang.eq.2 .or. lang.eq.3) then
+#ifndef LANG0
+c write (iout,*) "Calling sd_verlet_setup: 1"
+c Rescale the stochastic forces and recalculate or restore
+c the matrices of tinker integrator
+ if (itime_scal.gt.maxflag_stoch) then
+ if (large) write (iout,'(a,i5,a)')
+ & "Calculate matrices for stochastic step;",
+ & " itime_scal ",itime_scal
+ if (lang.eq.2) then
+ call sd_verlet_p_setup
+ else
+ call sd_verlet_ciccotti_setup
+ endif
+ write (iout,'(2a,i3,a,i3,1h.)')
+ & "Warning: cannot store matrices for stochastic",
+ & " integration because the index",itime_scal,
+ & " is greater than",maxflag_stoch
+ write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",
+ & " integration Langevin algorithm for better efficiency."
+ else if (flag_stoch(itime_scal)) then
+ if (large) write (iout,'(a,i5,a,l1)')
+ & "Restore matrices for stochastic step; itime_scal ",
+ & itime_scal," flag ",flag_stoch(itime_scal)
+ do i=1,dimen3
+ do j=1,dimen3
+ pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
+ afric_mat(i,j)=afric0_mat(i,j,itime_scal)
+ vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
+ prand_mat(i,j)=prand0_mat(i,j,itime_scal)
+ vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
+ vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
+ enddo
+ enddo
+ else
+ if (large) write (iout,'(2a,i5,a,l1)')
+ & "Calculate & store matrices for stochastic step;",
+ & " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
+ if (lang.eq.2) then
+ call sd_verlet_p_setup
+ else
+ call sd_verlet_ciccotti_setup
+ endif
+ flag_stoch(ifac_time)=.true.
+ do i=1,dimen3
+ do j=1,dimen3
+ pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
+ afric0_mat(i,j,itime_scal)=afric_mat(i,j)
+ vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
+ prand0_mat(i,j,itime_scal)=prand_mat(i,j)
+ vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
+ vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
+ enddo
+ enddo
+ endif
+ fac_time=1.0d0/dsqrt(fac_time)
+ do i=1,dimen3
+ stochforcvec(i)=fac_time*stochforcvec(i)
+ enddo
+#endif
+ else if (lang.eq.1) then
+c Rescale the accelerations due to stochastic forces
+ fac_time=1.0d0/dsqrt(fac_time)
+ do i=1,dimen3
+ d_as_work(i)=d_as_work(i)*fac_time
+ enddo
+ endif
+ if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)')
+ & "itime",itime," Timestep scaled down to ",
+ & d_time," ifac_time",ifac_time," itime_scal",itime_scal
+ else
+c Second step of the velocity Verlet algorithm
+ if (lang.eq.2) then
+#ifndef LANG0
+ call sd_verlet2
+#endif
+ else if (lang.eq.3) then
+#ifndef LANG0
+ call sd_verlet2_ciccotti
+#endif
+ else if (lang.eq.1) then
+ call sddir_verlet2
+ else if (tnp1) then
+ call tnp1_step2
+ else if (tnp) then
+ call tnp_step2
+ else
+ call verlet2
+ if (tnh) then
+ call kinetic(EK)
+ call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8)
+ do i=0,2*nres
+ do j=1,3
+ d_t(j,i)=d_t(j,i)*scale_nh
+ enddo
+ enddo
+ endif
+ endif
+ if (rattle) call rattle2
+ totT=totT+d_time
+ if (d_time.ne.d_time0) then
+ d_time=d_time0
+#ifndef LANG0
+ if (lang.eq.2 .or. lang.eq.3) then
+ if (large) write (iout,'(a)')
+ & "Restore original matrices for stochastic step"
+c write (iout,*) "Calling sd_verlet_setup: 2"
+c Restore the matrices of tinker integrator if the time step has been restored
+ do i=1,dimen3
+ do j=1,dimen3
+ pfric_mat(i,j)=pfric0_mat(i,j,0)
+ afric_mat(i,j)=afric0_mat(i,j,0)
+ vfric_mat(i,j)=vfric0_mat(i,j,0)
+ prand_mat(i,j)=prand0_mat(i,j,0)
+ vrand_mat1(i,j)=vrand0_mat1(i,j,0)
+ vrand_mat2(i,j)=vrand0_mat2(i,j,0)
+ enddo
+ enddo
+ endif
+#endif
+ endif
+ scale=.false.
+ endif
+ enddo
+c Calculate the kinetic and the total energy and the kinetic temperature
+ if (tnp .or. tnp1) then
+ do i=0,2*nres
+ do j=1,3
+ d_t_old(j,i)=d_t(j,i)
+ d_t(j,i)=d_t(j,i)/s_np
+ enddo
+ enddo
+ endif
+ call kinetic(EK)
+ totE=EK+potE
+c diagnostics
+c call kinetic1(EK1)
+c write (iout,*) "step",itime," EK",EK," EK1",EK1
+c end diagnostics
+c Couple the system to Berendsen bath if needed
+ if (tbf .and. lang.eq.0) then
+ call verlet_bath
+ endif
+ kinetic_T=2.0d0/(dimen3*Rb)*EK
+c Backup the coordinates, velocities, and accelerations
+ do i=0,2*nres
+ do j=1,3
+ dc_old(j,i)=dc(j,i)
+ if(.not.(tnp .or. tnp1)) d_t_old(j,i)=d_t(j,i)
+ d_a_old(j,i)=d_a(j,i)
+ enddo
+ enddo
+ if (ntwe.ne.0) then
+ if (mod(itime,ntwe).eq.0) then
+
+ if(tnp .or. tnp1) then
+ HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen3)
+ 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)
+cd write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
+ hhh=h
+ endif
+
+ if(tnh) then
+ HNose1=Hnose_nh(EK,potE)
+ H=HNose1-H0
+ hhh=h
+cd write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
+ endif
+
+ if (large) then
+ itnp=0
+ do j=1,3
+ itnp=itnp+1
+ vtnp(itnp)=d_t(j,0)
+ enddo
+ do i=nnt,nct-1
+ do j=1,3
+ itnp=itnp+1
+ vtnp(itnp)=d_t(j,i)
+ enddo
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ inres=i+nres
+ do j=1,3
+ itnp=itnp+1
+ vtnp(itnp)=d_t(j,inres)
+ enddo
+ endif
+ enddo
+
+c Transform velocities from UNRES coordinate space to cartesian and Gvec
+c eigenvector space
+
+ do i=1,dimen3
+ vtnp_(i)=0.0d0
+ vtnp_a(i)=0.0d0
+ do j=1,dimen3
+ vtnp_(i)=vtnp_(i)+Gvec(j,i)*vtnp(j)
+ vtnp_a(i)=vtnp_a(i)+A(i,j)*vtnp(j)
+ enddo
+ vtnp_(i)=vtnp_(i)*dsqrt(geigen(i))
+ enddo
+
+ do i=1,dimen3
+ write (iout,'("WWW",i3,3f10.5)') i,vtnp(i),vtnp_(i),vtnp_a(i)
+ enddo
+
+ write (iout,*) "Velocities, step 2"
+ 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
+ endif
+ endif
+ endif
+ return
+ end
+c-------------------------------------------------------------------------------
+ subroutine RESPA_step(itime)
+c-------------------------------------------------------------------------------
+c Perform a single RESPA step.
+c-------------------------------------------------------------------------------
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+#ifdef MPI
+ include 'mpif.h'
+ integer IERROR,ERRCODE
+#endif
+ include 'COMMON.SETUP'
+ include 'COMMON.CONTROL'
+ include 'COMMON.VAR'
+ include 'COMMON.MD'
+#ifndef LANG0
+ include 'COMMON.LANGEVIN'
+#else
+ include 'COMMON.LANGEVIN.lang0'
+#endif
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.NAMES'
+ include 'COMMON.TIME1'
+ double precision energia_short(0:n_ene),
+ & energia_long(0:n_ene)
+ double precision cm(3),L(3),vcm(3),incr(3)
+ double precision dc_old0(3,0:maxres2),d_t_old0(3,0:maxres2),
+ & d_a_old0(3,0:maxres2)
+ integer ilen,count,rstcount
+ external ilen
+ character*50 tytul
+ integer maxcount_scale /10/
+ common /gucio/ cm,energia_short
+ double precision stochforcvec(MAXRES6)
+ common /stochcalc/ stochforcvec
+ integer itime
+ logical scale
+ double precision vtnp(maxres6), vtnp_(maxres6), vtnp_a(maxres6)
+ common /cipiszcze/ itt
+ itt=itime
+ if (ntwe.ne.0) then
+ if (large.and. mod(itime,ntwe).eq.0) then
+ write (iout,*) "***************** RESPA itime",itime
+ write (iout,*) "Cartesian and internal coordinates: step 0"
+c call cartprint
+ call pdbout(0.0d0,"cipiszcze",iout)
+ call intout
+ write (iout,*) "Accelerations from long-range forces"
+ do i=0,nres
+ write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
+ & (d_a(j,i+nres),j=1,3)
+ enddo
+ write (iout,*) "Velocities, step 0"
+ 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
+ endif
+ endif
+c
+c Perform the initial RESPA step (increment velocities)
+c write (iout,*) "*********************** RESPA ini"
+ if (tnp1) then
+ call tnp_respa_step1
+ else if (tnp) then
+ call tnp_respa_step1
+ else
+ if (tnh.and..not.xiresp) then
+ call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8)
+ do i=0,2*nres
+ do j=1,3
+ d_t(j,i)=d_t(j,i)*scale_nh
+ enddo
+ enddo
+ endif
+ call RESPA_vel
+ endif
+
+cd if(tnp .or. tnp1) then
+cd write (iout,'(a,3f)') "EE1 NP S, pi",totT, s_np, pi_np
+cd endif
+
+ if (ntwe.ne.0) then
+ if (mod(itime,ntwe).eq.0 .and. large) then
+ write (iout,*) "Velocities, end"
+ 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
+ endif
+ endif
+c Compute the short-range forces
+#ifdef MPI
+ tt0 =MPI_Wtime()
+#else
+ tt0 = tcpu()
+#endif
+C 7/2/2009 commented out
+c call zerograd
+c call etotal_short(energia_short)
+ if (tnp.or.tnp1) potE=energia_short(0)
+c call cartgrad
+c call lagrangian
+C 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
+ do i=0,2*nres
+ do j=1,3
+ d_a(j,i)=d_a_short(j,i)
+ enddo
+ enddo
+ if (ntwe.ne.0) then
+ if (large.and. mod(itime,ntwe).eq.0) then
+ write (iout,*) "energia_short",energia_short(0)
+ write (iout,*) "Accelerations from short-range forces"
+ do i=0,nres
+ write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
+ & (d_a(j,i+nres),j=1,3)
+ enddo
+ endif
+ endif
+#ifdef MPI
+ t_enegrad=t_enegrad+MPI_Wtime()-tt0
+#else
+ t_enegrad=t_enegrad+tcpu()-tt0
+#endif
+ do i=0,2*nres
+ do j=1,3
+ dc_old(j,i)=dc(j,i)
+ if(.not.(tnp .or. tnp1)) d_t_old(j,i)=d_t(j,i)
+ d_a_old(j,i)=d_a(j,i)
+ enddo
+ enddo
+c 6/30/08 A-MTS: attempt at increasing the split number
+ do i=0,2*nres
+ do j=1,3
+ dc_old0(j,i)=dc_old(j,i)
+ d_t_old0(j,i)=d_t_old(j,i)
+ d_a_old0(j,i)=d_a_old(j,i)
+ enddo
+ enddo
+ if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
+ if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
+c
+ scale=.true.
+ d_time0=d_time
+ do while (scale)
+
+ scale=.false.
+c write (iout,*) "itime",itime," ntime_split",ntime_split
+c Split the time step
+ d_time=d_time0/ntime_split
+c Perform the short-range RESPA steps (velocity Verlet increments of
+c positions and velocities using short-range forces)
+c write (iout,*) "*********************** RESPA split"
+ do itsplit=1,ntime_split
+ if (lang.eq.1) then
+ call sddir_precalc
+ else if (lang.eq.2 .or. lang.eq.3) then
+#ifndef LANG0
+ call stochastic_force(stochforcvec)
+#else
+ write (iout,*)
+ & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
+#ifdef MPI
+ call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
+#endif
+ stop
+#endif
+ endif
+c First step of the velocity Verlet algorithm
+ if (lang.eq.2) then
+#ifndef LANG0
+ call sd_verlet1
+#endif
+ else if (lang.eq.3) then
+#ifndef LANG0
+ call sd_verlet1_ciccotti
+#endif
+ else if (lang.eq.1) then
+ call sddir_verlet1
+ else if (tnp1) then
+ call tnp1_respa_i_step1
+ else if (tnp) then
+ call tnp_respa_i_step1
+ else
+ if (tnh.and.xiresp) then
+ call kinetic(EK)
+ call nhcint(EK,scale_nh,wdtii,wdtii2,wdtii4,wdtii8)
+ do i=0,2*nres
+ do j=1,3
+ d_t_old(j,i)=d_t_old(j,i)*scale_nh
+ enddo
+ enddo
+cd write(iout,*) "SSS1",itsplit,EK,scale_nh
+ endif
+ call verlet1
+ endif
+c Build the chain from the newly calculated coordinates
+ call chainbuild_cart
+ if (rattle) call rattle1
+ if (ntwe.ne.0) then
+ if (large.and. mod(itime,ntwe).eq.0) then
+ write (iout,*) "***** ITSPLIT",itsplit
+ write (iout,*) "Cartesian and internal coordinates: step 1"
+ call pdbout(0.0d0,"cipiszcze",iout)
+c call cartprint
+ call intout
+ write (iout,*) "Velocities, step 1"
+ 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
+ endif
+ endif
+#ifdef MPI
+ tt0 = MPI_Wtime()
+#else
+ tt0 = tcpu()
+#endif
+c Calculate energy and forces
+ call zerograd
+ call etotal_short(energia_short)
+ E_old=potE
+ potE=energia_short(0)
+#ifdef TIMING_ENE
+#ifdef MPI
+ t_eshort=t_eshort+MPI_Wtime()-tt0
+#else
+ t_eshort=t_eshort+tcpu()-tt0
+#endif
+#endif
+ call cartgrad
+c Get the new accelerations
+ call lagrangian
+C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
+ do i=0,2*nres
+ do j=1,3
+ d_a_short(j,i)=d_a(j,i)
+ enddo
+ enddo
+ if (ntwe.ne.0) then
+ if (large.and. mod(itime,ntwe).eq.0) then
+ write (iout,*)"energia_short",energia_short(0)
+ write (iout,*) "Accelerations from short-range forces"
+ do i=0,nres
+ write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
+ & (d_a(j,i+nres),j=1,3)
+ enddo
+ endif
+ endif
+c 6/30/08 A-MTS
+c Determine maximum acceleration and scale down the timestep if needed
+ call max_accel
+ amax=amax/ntime_split**2
+ call predict_edrift(epdrift)
+ if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0)
+ & write (iout,*) "amax",amax," damax",damax,
+ & " epdrift",epdrift," epdriftmax",epdriftmax
+c Exit loop and try with increased split number if the change of
+c acceleration is too big
+ if (amax.gt.damax .or. epdrift.gt.edriftmax) then
+ if (ntime_split.lt.maxtime_split) then
+ scale=.true.
+ ntime_split=ntime_split*2
+ do i=0,2*nres
+ do j=1,3
+ dc_old(j,i)=dc_old0(j,i)
+ d_t_old(j,i)=d_t_old0(j,i)
+ d_a_old(j,i)=d_a_old0(j,i)
+ enddo
+ enddo
+ write (iout,*) "acceleration/energy drift too large",amax,
+ & epdrift," split increased to ",ntime_split," itime",itime,
+ & " itsplit",itsplit
+ exit
+ else
+ write (iout,*)
+ & "Uh-hu. Bumpy landscape. Maximum splitting number",
+ & maxtime_split,
+ & " already reached!!! Trying to carry on!"
+ endif
+ endif
+#ifdef MPI
+ t_enegrad=t_enegrad+MPI_Wtime()-tt0
+#else
+ t_enegrad=t_enegrad+tcpu()-tt0
+#endif
+c Second step of the velocity Verlet algorithm
+ if (lang.eq.2) then
+#ifndef LANG0
+ call sd_verlet2
+#endif
+ else if (lang.eq.3) then
+#ifndef LANG0
+ call sd_verlet2_ciccotti
+#endif
+ else if (lang.eq.1) then
+ call sddir_verlet2
+ else if (tnp1) then
+ call tnp1_respa_i_step2
+ else if (tnp) then
+ call tnp_respa_i_step2
+ else
+ call verlet2
+ if (tnh.and.xiresp) then
+ call kinetic(EK)
+ call nhcint(EK,scale_nh,wdtii,wdtii2,wdtii4,wdtii8)
+ do i=0,2*nres
+ do j=1,3
+ d_t(j,i)=d_t(j,i)*scale_nh
+ enddo
+ enddo
+cd write(iout,*) "SSS2",itsplit,EK,scale_nh
+ endif
+ endif
+ if (rattle) call rattle2
+c Backup the coordinates, velocities, and accelerations
+ if (tnp .or. tnp1) then
+ do i=0,2*nres
+ do j=1,3
+ d_t_old(j,i)=d_t(j,i)
+ if (tnp) d_t(j,i)=d_t(j,i)/s_np
+ if (tnp1) d_t(j,i)=d_t(j,i)/s_np
+ enddo
+ enddo
+ endif
+
+ do i=0,2*nres
+ do j=1,3
+ dc_old(j,i)=dc(j,i)
+ if(.not.(tnp .or. tnp1)) d_t_old(j,i)=d_t(j,i)
+ d_a_old(j,i)=d_a(j,i)
+ enddo
+ enddo
+ enddo
+
+ enddo ! while scale
+
+c Restore the time step
+ d_time=d_time0
+c Compute long-range forces
+#ifdef MPI
+ tt0 =MPI_Wtime()
+#else
+ tt0 = tcpu()
+#endif
+ call zerograd
+ call etotal_long(energia_long)
+ E_long=energia_long(0)
+ potE=energia_short(0)+energia_long(0)
+#ifdef TIMING_ENE
+#ifdef MPI
+ t_elong=t_elong+MPI_Wtime()-tt0
+#else
+ t_elong=t_elong+tcpu()-tt0
+#endif
+#endif
+ call cartgrad
+ call lagrangian
+#ifdef MPI
+ t_enegrad=t_enegrad+MPI_Wtime()-tt0
+#else
+ t_enegrad=t_enegrad+tcpu()-tt0
+#endif
+c Compute accelerations from long-range forces
+ if (ntwe.ne.0) then
+ if (large.and. mod(itime,ntwe).eq.0) then
+ write (iout,*) "energia_long",energia_long(0)
+ write (iout,*) "Cartesian and internal coordinates: step 2"
+c call cartprint
+ call pdbout(0.0d0,"cipiszcze",iout)
+ call intout
+ write (iout,*) "Accelerations from long-range forces"
+ do i=0,nres
+ write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
+ & (d_a(j,i+nres),j=1,3)
+ enddo
+ write (iout,*) "Velocities, step 2"
+ 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
+ endif
+ endif
+c Compute the final RESPA step (increment velocities)
+c write (iout,*) "*********************** RESPA fin"
+ if (tnp1) then
+ call tnp_respa_step2
+ else if (tnp) then
+ call tnp_respa_step2
+ else
+ call RESPA_vel
+ if (tnh.and..not.xiresp) then
+ call kinetic(EK)
+ call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8)
+ do i=0,2*nres
+ do j=1,3
+ d_t(j,i)=d_t(j,i)*scale_nh
+ enddo
+ enddo
+ endif
+ endif
+
+ if (tnp .or. tnp1) then
+ do i=0,2*nres
+ do j=1,3
+ d_t(j,i)=d_t_old(j,i)/s_np
+ enddo
+ enddo
+ endif
+
+c Compute the complete potential energy
+ do i=0,n_ene
+ potEcomp(i)=energia_short(i)+energia_long(i)
+ enddo
+ potE=potEcomp(0)-potEcomp(20)
+c potE=energia_short(0)+energia_long(0)
+ totT=totT+d_time
+c Calculate the kinetic and the total energy and the kinetic temperature
+ call kinetic(EK)
+ totE=EK+potE
+c Couple the system to Berendsen bath if needed
+ if (tbf .and. lang.eq.0) then
+ call verlet_bath
+ endif
+ kinetic_T=2.0d0/(dimen3*Rb)*EK
+c Backup the coordinates, velocities, and accelerations
+ if (ntwe.ne.0) then
+ if (mod(itime,ntwe).eq.0 .and. large) then
+ write (iout,*) "Velocities, end"
+ 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
+ endif
+
+ if (mod(itime,ntwe).eq.0) then
+
+ if(tnp .or. tnp1) then
+#ifndef G77
+ 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,
+ & E_long,energia_short(0)
+#endif
+ HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen3)
+ 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)
+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
+cd write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
+ hhh=h
+ endif
+
+
+ if (large) then
+ itnp=0
+ do j=1,3
+ itnp=itnp+1
+ vtnp(itnp)=d_t(j,0)
+ enddo
+ do i=nnt,nct-1
+ do j=1,3
+ itnp=itnp+1
+ vtnp(itnp)=d_t(j,i)
+ enddo
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ inres=i+nres
+ do j=1,3
+ itnp=itnp+1
+ vtnp(itnp)=d_t(j,inres)
+ enddo
+ endif
+ enddo
+
+c Transform velocities from UNRES coordinate space to cartesian and Gvec
+c eigenvector space
+
+ do i=1,dimen3
+ vtnp_(i)=0.0d0
+ vtnp_a(i)=0.0d0
+ do j=1,dimen3
+ vtnp_(i)=vtnp_(i)+Gvec(j,i)*vtnp(j)
+ vtnp_a(i)=vtnp_a(i)+A(i,j)*vtnp(j)
+ enddo
+ vtnp_(i)=vtnp_(i)*dsqrt(geigen(i))
+ enddo
+
+ do i=1,dimen3
+ write (iout,'("WWW",i3,3f10.5)') i,vtnp(i),vtnp_(i),vtnp_a(i)
+ enddo
+
+ endif
+ endif
+ endif
+
+
+ return
+ end
+c---------------------------------------------------------------------
+ subroutine RESPA_vel
+c First and last RESPA step (incrementing velocities using long-range
+c forces).
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CONTROL'
+ include 'COMMON.VAR'
+ include 'COMMON.MD'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.NAMES'
+ do j=1,3
+ d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time
+ enddo
+ do i=nnt,nct-1
+ do j=1,3
+ d_t(j,i)=d_t(j,i)+0.5d0*d_a(j,i)*d_time
+ enddo
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ inres=i+nres
+ do j=1,3
+ d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
+ enddo
+ endif
+ enddo
+ return
+ end
+c-----------------------------------------------------------------
+ subroutine verlet1
+c Applying velocity Verlet algorithm - step 1 to coordinates
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CONTROL'
+ include 'COMMON.VAR'
+ include 'COMMON.MD'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.NAMES'
+ double precision adt,adt2
+
+#ifdef DEBUG
+ write (iout,*) "VELVERLET1 START: DC"
+ do i=0,nres
+ write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
+ & (dc(j,i+nres),j=1,3)
+ enddo
+#endif
+ do j=1,3
+ adt=d_a_old(j,0)*d_time
+ adt2=0.5d0*adt
+ dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
+ d_t_new(j,0)=d_t_old(j,0)+adt2
+ d_t(j,0)=d_t_old(j,0)+adt
+ enddo
+ do i=nnt,nct-1
+ do j=1,3
+ adt=d_a_old(j,i)*d_time
+ adt2=0.5d0*adt
+ dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
+ d_t_new(j,i)=d_t_old(j,i)+adt2
+ d_t(j,i)=d_t_old(j,i)+adt
+ enddo
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ inres=i+nres
+ do j=1,3
+ adt=d_a_old(j,inres)*d_time
+ adt2=0.5d0*adt
+ dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
+ d_t_new(j,inres)=d_t_old(j,inres)+adt2
+ d_t(j,inres)=d_t_old(j,inres)+adt
+ enddo
+ endif
+ enddo
+#ifdef DEBUG
+ write (iout,*) "VELVERLET1 END: DC"
+ do i=0,nres
+ write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
+ & (dc(j,i+nres),j=1,3)
+ enddo
+#endif
+ return
+ end
+c---------------------------------------------------------------------
+ subroutine verlet2
+c Step 2 of the velocity Verlet algorithm: update velocities
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CONTROL'
+ include 'COMMON.VAR'
+ include 'COMMON.MD'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.NAMES'
+ do j=1,3
+ d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time
+ enddo
+ do i=nnt,nct-1
+ do j=1,3
+ d_t(j,i)=d_t_new(j,i)+0.5d0*d_a(j,i)*d_time
+ enddo
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ inres=i+nres
+ do j=1,3
+ d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
+ enddo
+ endif
+ enddo
+ return
+ end
+c-----------------------------------------------------------------
+ subroutine sddir_precalc
+c Applying velocity Verlet algorithm - step 1 to coordinates
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+#ifdef MPI
+ include 'mpif.h'
+#endif
+ include 'COMMON.CONTROL'
+ include 'COMMON.VAR'
+ include 'COMMON.MD'
+#ifndef LANG0
+ include 'COMMON.LANGEVIN'
+#else
+ include 'COMMON.LANGEVIN.lang0'
+#endif
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.NAMES'
+ include 'COMMON.TIME1'
+ double precision stochforcvec(MAXRES6)
+ common /stochcalc/ stochforcvec
+c
+c Compute friction and stochastic forces
+c
+#ifdef MPI
+ time00=MPI_Wtime()
+#else
+ time00=tcpu()
+#endif
+ call friction_force
+#ifdef MPI
+ time_fric=time_fric+MPI_Wtime()-time00
+ time00=MPI_Wtime()
+#else
+ time_fric=time_fric+tcpu()-time00
+ time00=tcpu()
+#endif
+ call stochastic_force(stochforcvec)
+#ifdef MPI
+ time_stoch=time_stoch+MPI_Wtime()-time00
+#else
+ time_stoch=time_stoch+tcpu()-time00
+#endif
+c
+c Compute the acceleration due to friction forces (d_af_work) and stochastic
+c forces (d_as_work)
+c
+ call ginv_mult(fric_work, d_af_work)
+ call ginv_mult(stochforcvec, d_as_work)
+ return
+ end
+c---------------------------------------------------------------------
+ subroutine sddir_verlet1
+c Applying velocity Verlet algorithm - step 1 to velocities
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CONTROL'
+ include 'COMMON.VAR'
+ include 'COMMON.MD'
+#ifndef LANG0
+ include 'COMMON.LANGEVIN'
+#else
+ include 'COMMON.LANGEVIN.lang0'
+#endif
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.NAMES'
+c Revised 3/31/05 AL: correlation between random contributions to
+c position and velocity increments included.
+ double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3)
+ double precision adt,adt2
+c
+c Add the contribution from BOTH friction and stochastic force to the
+c coordinates, but ONLY the contribution from the friction forces to velocities
+c
+ do j=1,3
+ adt=(d_a_old(j,0)+d_af_work(j))*d_time
+ adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
+ dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
+ d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
+ d_t(j,0)=d_t_old(j,0)+adt
+ enddo
+ ind=3
+ do i=nnt,nct-1
+ do j=1,3
+ adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
+ adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
+ dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
+ d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
+ d_t(j,i)=d_t_old(j,i)+adt
+ enddo
+ ind=ind+3
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ inres=i+nres
+ do j=1,3
+ adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
+ adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
+ dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
+ d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
+ d_t(j,inres)=d_t_old(j,inres)+adt
+ enddo
+ ind=ind+3
+ endif
+ enddo
+ return
+ end
+c---------------------------------------------------------------------
+ subroutine sddir_verlet2
+c Calculating the adjusted velocities for accelerations
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CONTROL'
+ include 'COMMON.VAR'
+ include 'COMMON.MD'
+#ifndef LANG0
+ include 'COMMON.LANGEVIN'
+#else
+ include 'COMMON.LANGEVIN.lang0'
+#endif
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.NAMES'
+ double precision stochforcvec(MAXRES6),d_as_work1(MAXRES6)
+ double precision cos60 /0.5d0/, sin60 /0.86602540378443864676d0/
+c Revised 3/31/05 AL: correlation between random contributions to
+c position and velocity increments included.
+c The correlation coefficients are calculated at low-friction limit.
+c Also, friction forces are now not calculated with new velocities.
+
+c call friction_force
+ call stochastic_force(stochforcvec)
+c
+c Compute the acceleration due to friction forces (d_af_work) and stochastic
+c forces (d_as_work)
+c
+ call ginv_mult(stochforcvec, d_as_work1)
+
+c
+c Update velocities
+c
+ do j=1,3
+ d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j))
+ & +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
+ enddo
+ ind=3
+ do i=nnt,nct-1
+ do j=1,3
+ d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j))
+ & +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
+ enddo
+ ind=ind+3
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ inres=i+nres
+ do j=1,3
+ d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres)
+ & +d_af_work(ind+j))+sin60*d_as_work(ind+j)
+ & +cos60*d_as_work1(ind+j))*d_time
+ enddo
+ ind=ind+3
+ endif
+ enddo
+ return
+ end
+c---------------------------------------------------------------------
+ subroutine max_accel
+c
+c Find the maximum difference in the accelerations of the the sites
+c at the beginning and the end of the time step.
+c
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CONTROL'
+ include 'COMMON.VAR'
+ include 'COMMON.MD'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ double precision aux(3),accel(3),accel_old(3),dacc
+ do j=1,3
+c aux(j)=d_a(j,0)-d_a_old(j,0)
+ accel_old(j)=d_a_old(j,0)
+ accel(j)=d_a(j,0)
+ enddo
+ amax=0.0d0
+ do i=nnt,nct
+c Backbone
+ if (i.lt.nct) then
+c 7/3/08 changed to asymmetric difference
+ do j=1,3
+c accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
+ accel_old(j)=accel_old(j)+0.5d0*d_a_old(j,i)
+ accel(j)=accel(j)+0.5d0*d_a(j,i)
+c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
+ if (dabs(accel(j)).gt.dabs(accel_old(j))) then
+ dacc=dabs(accel(j)-accel_old(j))
+ if (dacc.gt.amax) amax=dacc
+ endif
+ enddo
+ endif
+ enddo
+c Side chains
+ do j=1,3
+c accel(j)=aux(j)
+ accel_old(j)=d_a_old(j,0)
+ accel(j)=d_a(j,0)
+ enddo
+ if (nnt.eq.2) then
+ do j=1,3
+ accel_old(j)=accel_old(j)+d_a_old(j,1)
+ accel(j)=accel(j)+d_a(j,1)
+ enddo
+ endif
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ do j=1,3
+c accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
+ accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
+ accel(j)=accel(j)+d_a(j,i+nres)
+ enddo
+ endif
+ do j=1,3
+c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
+ if (dabs(accel(j)).gt.dabs(accel_old(j))) then
+ dacc=dabs(accel(j)-accel_old(j))
+ if (dacc.gt.amax) amax=dacc
+ endif
+ enddo
+ do j=1,3
+ accel_old(j)=accel_old(j)+d_a_old(j,i)
+ accel(j)=accel(j)+d_a(j,i)
+c aux(j)=aux(j)+d_a(j,i)-d_a_old(j,i)
+ enddo
+ enddo
+ return
+ end
+c---------------------------------------------------------------------
+ subroutine predict_edrift(epdrift)
+c
+c Predict the drift of the potential energy
+c
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CONTROL'
+ include 'COMMON.VAR'
+ include 'COMMON.MD'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.MUCA'
+ double precision epdrift,epdriftij
+c Drift of the potential energy
+ epdrift=0.0d0
+ do i=nnt,nct
+c Backbone
+ if (i.lt.nct) then
+ do j=1,3
+ epdriftij=dabs((d_a(j,i)-d_a_old(j,i))*gcart(j,i))
+ if (lmuca) epdriftij=epdriftij*factor
+c write (iout,*) "back",i,j,epdriftij
+ if (epdriftij.gt.epdrift) epdrift=epdriftij
+ enddo
+ endif
+c Side chains
+ if (itype(i).ne.10) then
+ do j=1,3
+ epdriftij=
+ & dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
+ if (lmuca) epdriftij=epdriftij*factor
+c write (iout,*) "side",i,j,epdriftij
+ if (epdriftij.gt.epdrift) epdrift=epdriftij
+ enddo
+ endif
+ enddo
+ epdrift=0.5d0*epdrift*d_time*d_time
+c write (iout,*) "epdrift",epdrift
+ return
+ end
+c-----------------------------------------------------------------------
+ subroutine verlet_bath
+c
+c Coupling to the thermostat by using the Berendsen algorithm
+c
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CONTROL'
+ include 'COMMON.VAR'
+ include 'COMMON.MD'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.NAMES'
+ double precision T_half,fact
+c
+ T_half=2.0d0/(dimen3*Rb)*EK
+ fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0))
+c write(iout,*) "T_half", T_half
+c write(iout,*) "EK", EK
+c write(iout,*) "fact", fact
+ do j=1,3
+ d_t(j,0)=fact*d_t(j,0)
+ enddo
+ do i=nnt,nct-1
+ do j=1,3
+ d_t(j,i)=fact*d_t(j,i)
+ enddo
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ inres=i+nres
+ do j=1,3
+ d_t(j,inres)=fact*d_t(j,inres)
+ enddo
+ endif
+ enddo
+ return
+ end
+c---------------------------------------------------------
+ subroutine init_MD
+c Set up the initial conditions of a MD simulation
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+#ifdef MP
+ include 'mpif.h'
+ character*16 form
+ integer IERROR,ERRCODE
+#endif
+ include 'COMMON.SETUP'
+ include 'COMMON.CONTROL'
+ include 'COMMON.VAR'
+ include 'COMMON.MD'
+#ifndef LANG0
+ include 'COMMON.LANGEVIN'
+#else
+ include 'COMMON.LANGEVIN.lang0'
+#endif
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.NAMES'
+ include 'COMMON.REMD'
+ real*8 energia_long(0:n_ene),
+ & energia_short(0:n_ene),vcm(3),incr(3),E_short
+ double precision cm(3),L(3),xv,sigv,lowb,highb
+ double precision varia(maxvar)
+ character*256 qstr
+ integer ilen
+ external ilen
+ character*50 tytul
+ logical file_exist
+ common /gucio/ cm
+ d_time0=d_time
+c write(iout,*) "d_time", d_time
+c Compute the standard deviations of stochastic forces for Langevin dynamics
+c if the friction coefficients do not depend on surface area
+ if (lang.gt.0 .and. .not.surfarea) then
+ do i=nnt,nct-1
+ stdforcp(i)=stdfp*dsqrt(gamp)
+ enddo
+ do i=nnt,nct
+ stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(itype(i)))
+ enddo
+ endif
+c Open the pdb file for snapshotshots
+#ifdef MPI
+ if(mdpdb) then
+ if (ilen(tmpdir).gt.0)
+ & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
+ & liczba(:ilen(liczba))//".pdb")
+ open(ipdb,
+ & file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
+ & //".pdb")
+ else
+#ifdef NOXDR
+ if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
+ & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
+ & liczba(:ilen(liczba))//".x")
+ cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
+ & //".x"
+#else
+ if (ilen(tmpdir).gt.0 .and. (me.eq.king .or. .not.traj1file))
+ & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD"//
+ & liczba(:ilen(liczba))//".cx")
+ cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
+ & //".cx"
+#endif
+ endif
+#else
+ if(mdpdb) then
+ if (ilen(tmpdir).gt.0)
+ & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.pdb")
+ open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
+ else
+ if (ilen(tmpdir).gt.0)
+ & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_MD.cx")
+ cartname=prefix(:ilen(prefix))//"_MD.cx"
+ endif
+#endif
+ if (usampl) then
+ write (qstr,'(256(1h ))')
+ ipos=1
+ do i=1,nfrag
+ iq = qinfrag(i,iset)*10
+ iw = wfrag(i,iset)/100
+ if (iw.gt.0) then
+ if(me.eq.king.or..not.out1file)
+ & write (iout,*) "Frag",qinfrag(i,iset),wfrag(i,iset),iq,iw
+ write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
+ ipos=ipos+7
+ endif
+ enddo
+ do i=1,npair
+ iq = qinpair(i,iset)*10
+ iw = wpair(i,iset)/100
+ if (iw.gt.0) then
+ if(me.eq.king.or..not.out1file)
+ & write (iout,*) "Pair",i,qinpair(i,iset),wpair(i,iset),iq,iw
+ write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
+ ipos=ipos+7
+ endif
+ enddo
+c pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
+#ifdef NOXDR
+c cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
+#else
+c cartname=cartname(:ilen(cartname)-3)//qstr(:ipos-1)//'.cx'
+#endif
+c statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
+ endif
+ icg=1
+ if (rest) then
+ if (restart1file) then
+ if (me.eq.king)
+ & inquire(file=mremd_rst_name,exist=file_exist)
+ write (*,*) me," Before broadcast: file_exist",file_exist
+#ifdef MPI
+ call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
+ & IERR)
+ write (*,*) me," After broadcast: file_exist",file_exist
+#endif
+c inquire(file=mremd_rst_name,exist=file_exist)
+ if(me.eq.king.or..not.out1file)
+ & write(iout,*) "Initial state read by master and distributed"
+ else
+ if (ilen(tmpdir).gt.0)
+ & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'
+ & //liczba(:ilen(liczba))//'.rst')
+ inquire(file=rest2name,exist=file_exist)
+ endif
+ if(file_exist) then
+ if(.not.restart1file) then
+ if(me.eq.king.or..not.out1file)
+ & write(iout,*) "Initial state will be read from file ",
+ & rest2name(:ilen(rest2name))
+ call readrst
+ endif
+ call rescale_weights(t_bath)
+ else
+ if(me.eq.king.or..not.out1file)then
+ if (restart1file) then
+ write(iout,*) "File ",mremd_rst_name(:ilen(mremd_rst_name)),
+ & " does not exist"
+ else
+ write(iout,*) "File ",rest2name(:ilen(rest2name)),
+ & " does not exist"
+ endif
+ write(iout,*) "Initial velocities randomly generated"
+ endif
+ call random_vel
+ totT=0.0d0
+ endif
+ else
+c Generate initial velocities
+ if(me.eq.king.or..not.out1file)
+ & write(iout,*) "Initial velocities randomly generated"
+ call random_vel
+ totT=0.0d0
+ endif
+c rest2name = prefix(:ilen(prefix))//'.rst'
+ if(me.eq.king.or..not.out1file)then
+ write (iout,*) "Initial velocities"
+ 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)
+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)
+c write (iout,*) "velocity of the center of the mass:"
+c write (iout,*) (vcm(j),j=1,3)
+ do j=1,3
+ d_t(j,0)=d_t(j,0)-vcm(j)
+ enddo
+c Removing the velocity of the center of mass
+ call vcm_vel(vcm)
+ 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
+ if(iranconf.ne.0) then
+ if (overlapsc) then
+ print *, 'Calling OVERLAP_SC'
+ call overlap_sc(fail)
+ endif
+
+ if (searchsc) then
+ call sc_move(2,nres-1,10,1d10,nft_sc,etot)
+ print *,'SC_move',nft_sc,etot
+ if(me.eq.king.or..not.out1file)
+ & write(iout,*) 'SC_move',nft_sc,etot
+ endif
+
+ if(dccart)then
+ print *, 'Calling MINIM_DC'
+ call minim_dc(etot,iretcode,nfun)
+ else
+ call geom_to_var(nvar,varia)
+ print *,'Calling MINIMIZE.'
+ call minimize(etot,varia,iretcode,nfun)
+ call var_to_geom(nvar,varia)
+ endif
+ if(me.eq.king.or..not.out1file)
+ & write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
+ endif
+ endif
+ call chainbuild_cart
+ call kinetic(EK)
+ if (tbf) then
+ call verlet_bath
+ endif
+ kinetic_T=2.0d0/(dimen3*Rb)*EK
+ if(me.eq.king.or..not.out1file)then
+ call cartprint
+ call intout
+ endif
+#ifdef MPI
+ tt0=MPI_Wtime()
+#else
+ tt0=tcpu()
+#endif
+ call zerograd
+ call etotal(potEcomp)
+#ifdef TIMING_ENE
+#ifdef MPI
+ t_etotal=t_etotal+MPI_Wtime()-tt0
+#else
+ t_etotal=t_etotal+tcpu()-tt0
+#endif
+#endif
+ potE=potEcomp(0)
+
+ if(tnp .or. tnp1) then
+ s_np=1.0
+ pi_np=0.0
+ HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen3)
+ H0=Hnose1
+ write(iout,*) 'H0= ',H0
+ endif
+
+ if(tnh) then
+ HNose1=Hnose_nh(EK,potE)
+ H0=HNose1
+ 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
+ call max_accel
+ if (amax*d_time .gt. dvmax) then
+ d_time=d_time*dvmax/amax
+ if(me.eq.king.or..not.out1file) write (iout,*)
+ & "Time step reduced to",d_time,
+ & " because of too large initial acceleration."
+ endif
+ if(me.eq.king.or..not.out1file)then
+ write(iout,*) "Potential energy and its components"
+ call enerprint(potEcomp)
+c write(iout,*) (potEcomp(i),i=0,n_ene)
+ endif
+ potE=potEcomp(0)-potEcomp(20)
+ totE=EK+potE
+ itime=0
+ if (ntwe.ne.0) call statout(itime)
+ if(me.eq.king.or..not.out1file)
+ & write (iout,'(/a/3(a25,1pe14.5/))') "Initial:",
+ & " Kinetic energy",EK," potential energy",potE,
+ & " total energy",totE," maximum acceleration ",
+ & amax
+ if (large) then
+ write (iout,*) "Initial coordinates"
+ do i=1,nres
+ write (iout,'(i3,3f10.5,3x,3f10.5)') i,(c(j,i),j=1,3),
+ & (c(j,i+nres),j=1,3)
+ enddo
+ write (iout,*) "Initial dC"
+ do i=0,nres
+ write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
+ & (dc(j,i+nres),j=1,3)
+ enddo
+ write (iout,*) "Initial velocities"
+ write (iout,"(13x,' backbone ',23x,' side chain')")
+ 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
+ write (iout,*) "Initial accelerations"
+ do i=0,nres
+c write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
+ write (iout,'(i3,3f15.10,3x,3f15.10)') i,(d_a(j,i),j=1,3),
+ & (d_a(j,i+nres),j=1,3)
+ enddo
+ endif
+ do i=0,2*nres
+ do j=1,3
+ dc_old(j,i)=dc(j,i)
+ d_t_old(j,i)=d_t(j,i)
+ d_a_old(j,i)=d_a(j,i)
+ enddo
+c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
+ enddo
+ if (RESPA) then
+#ifdef MPI
+ tt0 =MPI_Wtime()
+#else
+ tt0 = tcpu()
+#endif
+ call zerograd
+ call etotal_short(energia_short)
+#ifdef TIMING_ENE
+#ifdef MPI
+ t_eshort=t_eshort+MPI_Wtime()-tt0
+#else
+ t_eshort=t_eshort+tcpu()-tt0
+#endif
+#endif
+
+ if(tnp .or. tnp1) then
+ E_short=energia_short(0)
+ HNose1=Hnose(EK,s_np,E_short,pi_np,Q_np,t_bath,dimen3)
+ Csplit=Hnose1
+c Csplit =110
+c_new_var_csplit Csplit=H0-E_long
+c Csplit = H0-energia_short(0)
+ write(iout,*) 'Csplit= ',Csplit
+ endif
+
+
+ call cartgrad
+ call lagrangian
+ if(.not.out1file .and. large) then
+ write (iout,*) "energia_long",energia_long(0),
+ & " energia_short",energia_short(0),
+ & " total",energia_long(0)+energia_short(0)
+ write (iout,*) "Initial fast-force accelerations"
+ do i=0,nres
+ write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
+ & (d_a(j,i+nres),j=1,3)
+ enddo
+ endif
+C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
+ do i=0,2*nres
+ do j=1,3
+ d_a_short(j,i)=d_a(j,i)
+ enddo
+ enddo
+#ifdef MPI
+ tt0=MPI_Wtime()
+#else
+ tt0=tcpu()
+#endif
+ call zerograd
+ call etotal_long(energia_long)
+#ifdef TIMING_ENE
+#ifdef MPI
+ t_elong=t_elong+MPI_Wtime()-tt0
+#else
+ t_elong=t_elong+tcpu()-tt0
+#endif
+#endif
+ call cartgrad
+ call lagrangian
+ if(.not.out1file .and. large) then
+ write (iout,*) "energia_long",energia_long(0)
+ write (iout,*) "Initial slow-force accelerations"
+ do i=0,nres
+ write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
+ & (d_a(j,i+nres),j=1,3)
+ enddo
+ endif
+#ifdef MPI
+ t_enegrad=t_enegrad+MPI_Wtime()-tt0
+#else
+ t_enegrad=t_enegrad+tcpu()-tt0
+#endif
+ endif
+
+
+
+ return
+ end
+c-----------------------------------------------------------
+ subroutine random_vel
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CONTROL'
+ include 'COMMON.VAR'
+ include 'COMMON.MD'
+#ifndef LANG0
+ include 'COMMON.LANGEVIN'
+#else
+ include 'COMMON.LANGEVIN.lang0'
+#endif
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.NAMES'
+ include 'COMMON.TIME1'
+ double precision xv,sigv,lowb,highb
+c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m
+c First generate velocities in the eigenspace of the G matrix
+c write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
+c call flush(iout)
+c write (iout,*) "RANDOM_VEL dimen",dimen
+ xv=0.0d0
+ ii=0
+ do i=1,dimen
+ do k=1,3
+ ii=ii+1
+ sigv=dsqrt((Rb*t_bath)/geigen(i))
+ lowb=-5*sigv
+ highb=5*sigv
+ d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
+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
+c do i=1,dimen
+c do k=1,3
+c ii=ii+1
+c Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2
+c enddo
+c enddo
+c write (iout,*) "Ek from eigenvectors",Ek1
+c end diagnostics
+c Transform velocities to UNRES coordinate space
+ do k=0,2
+ do i=1,dimen
+ ind=(i-1)*3+k+1
+ d_t_work(ind)=0.0d0
+ do j=1,dimen
+ d_t_work(ind)=d_t_work(ind)
+ & +Gvec(i,j)*d_t_work_new((j-1)*3+k+1)
+ enddo
+c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
+c call flush(iout)
+ enddo
+ enddo
+c Transfer to the d_t vector
+ do j=1,3
+ d_t(j,0)=d_t_work(j)
+ enddo
+ ind=3
+ do i=nnt,nct-1
+ do j=1,3
+ ind=ind+1
+ d_t(j,i)=d_t_work(ind)
+ enddo
+ enddo
+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
+ ind=ind+1
+ d_t(j,i+nres)=d_t_work(ind)
+ enddo
+ endif
+ enddo
+c call kinetic(EK)
+c write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",
+c & 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1
+c call flush(iout)
+ return
+ end
+#ifndef LANG0
+c-----------------------------------------------------------
+ subroutine sd_verlet_p_setup
+c Sets up the parameters of stochastic Verlet algorithm
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+#ifdef MPI
+ include 'mpif.h'
+#endif
+ include 'COMMON.CONTROL'
+ include 'COMMON.VAR'
+ include 'COMMON.MD'
+#ifndef LANG0
+ include 'COMMON.LANGEVIN'
+#else
+ include 'COMMON.LANGEVIN.lang0'
+#endif
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.NAMES'
+ include 'COMMON.TIME1'
+ double precision emgdt(MAXRES6),
+ & pterm,vterm,rho,rhoc,vsig,
+ & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
+ & afric_vec(MAXRES6),prand_vec(MAXRES6),
+ & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
+ logical lprn /.false./
+ double precision zero /1.0d-8/, gdt_radius /0.05d0/
+ double precision ktm
+#ifdef MPI
+ tt0 = MPI_Wtime()
+#else
+ tt0 = tcpu()
+#endif
+c
+c AL 8/17/04 Code adapted from tinker
+c
+c Get the frictional and random terms for stochastic dynamics in the
+c eigenspace of mass-scaled UNRES friction matrix
+c
+ do i = 1, dimen
+ gdt = fricgam(i) * d_time
+c
+c Stochastic dynamics reduces to simple MD for zero friction
+c
+ if (gdt .le. zero) then
+ pfric_vec(i) = 1.0d0
+ vfric_vec(i) = d_time
+ afric_vec(i) = 0.5d0 * d_time * d_time
+ prand_vec(i) = 0.0d0
+ vrand_vec1(i) = 0.0d0
+ vrand_vec2(i) = 0.0d0
+c
+c Analytical expressions when friction coefficient is large
+c
+ else
+ if (gdt .ge. gdt_radius) then
+ egdt = dexp(-gdt)
+ pfric_vec(i) = egdt
+ vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
+ afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
+ pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
+ vterm = 1.0d0 - egdt**2
+ rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
+c
+c Use series expansions when friction coefficient is small
+c
+ else
+ gdt2 = gdt * gdt
+ gdt3 = gdt * gdt2
+ gdt4 = gdt2 * gdt2
+ gdt5 = gdt2 * gdt3
+ gdt6 = gdt3 * gdt3
+ gdt7 = gdt3 * gdt4
+ gdt8 = gdt4 * gdt4
+ gdt9 = gdt4 * gdt5
+ afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
+ & - gdt5/120.0d0 + gdt6/720.0d0
+ & - gdt7/5040.0d0 + gdt8/40320.0d0
+ & - gdt9/362880.0d0) / fricgam(i)**2
+ vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
+ pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
+ pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
+ & + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
+ & + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
+ & + 127.0d0*gdt9/90720.0d0
+ vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
+ & - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
+ & - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
+ & - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
+ rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
+ & - 17.0d0*gdt2/1280.0d0
+ & + 17.0d0*gdt3/6144.0d0
+ & + 40967.0d0*gdt4/34406400.0d0
+ & - 57203.0d0*gdt5/275251200.0d0
+ & - 1429487.0d0*gdt6/13212057600.0d0)
+ end if
+c
+c Compute the scaling factors of random terms for the nonzero friction case
+c
+ ktm = 0.5d0*d_time/fricgam(i)
+ psig = dsqrt(ktm*pterm) / fricgam(i)
+ vsig = dsqrt(ktm*vterm)
+ rhoc = dsqrt(1.0d0 - rho*rho)
+ prand_vec(i) = psig
+ vrand_vec1(i) = vsig * rho
+ vrand_vec2(i) = vsig * rhoc
+ end if
+ end do
+ if (lprn) then
+ write (iout,*)
+ & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
+ & " vrand_vec2"
+ do i=1,dimen
+ write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
+ & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
+ enddo
+ endif
+c
+c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
+c
+#ifndef LANG0
+ call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
+ call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
+ call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
+ call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
+ call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec1,vrand_mat1)
+ call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
+#endif
+#ifdef MPI
+ t_sdsetup=t_sdsetup+MPI_Wtime()
+#else
+ t_sdsetup=t_sdsetup+tcpu()-tt0
+#endif
+ return
+ end
+c-------------------------------------------------------------
+ subroutine eigtransf1(n,ndim,ab,d,c)
+ implicit none
+ integer n,ndim
+ double precision ab(ndim,ndim,n),c(ndim,n),d(ndim)
+ integer i,j,k
+ do i=1,n
+ do j=1,n
+ c(i,j)=0.0d0
+ do k=1,n
+ c(i,j)=c(i,j)+ab(k,j,i)*d(k)
+ enddo
+ enddo
+ enddo
+ return
+ end
+c-------------------------------------------------------------
+ subroutine eigtransf(n,ndim,a,b,d,c)
+ implicit none
+ integer n,ndim
+ double precision a(ndim,n),b(ndim,n),c(ndim,n),d(ndim)
+ integer i,j,k
+ do i=1,n
+ do j=1,n
+ c(i,j)=0.0d0
+ do k=1,n
+ c(i,j)=c(i,j)+a(i,k)*b(k,j)*d(k)
+ enddo
+ enddo
+ enddo
+ return
+ end
+c-------------------------------------------------------------
+ subroutine sd_verlet1
+c Applying stochastic velocity Verlet algorithm - step 1 to velocities
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CONTROL'
+ include 'COMMON.VAR'
+ include 'COMMON.MD'
+#ifndef LANG0
+ include 'COMMON.LANGEVIN'
+#else
+ include 'COMMON.LANGEVIN.lang0'
+#endif
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.NAMES'
+ double precision stochforcvec(MAXRES6)
+ common /stochcalc/ stochforcvec
+ logical lprn /.false./
+
+c write (iout,*) "dc_old"
+c do i=0,nres
+c write (iout,'(i5,3f10.5,5x,3f10.5)')
+c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
+c enddo
+ do j=1,3
+ dc_work(j)=dc_old(j,0)
+ d_t_work(j)=d_t_old(j,0)
+ d_a_work(j)=d_a_old(j,0)
+ enddo
+ ind=3
+ do i=nnt,nct-1
+ do j=1,3
+ dc_work(ind+j)=dc_old(j,i)
+ d_t_work(ind+j)=d_t_old(j,i)
+ d_a_work(ind+j)=d_a_old(j,i)
+ enddo
+ ind=ind+3
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ do j=1,3
+ dc_work(ind+j)=dc_old(j,i+nres)
+ d_t_work(ind+j)=d_t_old(j,i+nres)
+ d_a_work(ind+j)=d_a_old(j,i+nres)
+ enddo
+ ind=ind+3
+ endif
+ enddo
+#ifndef LANG0
+ if (lprn) then
+ write (iout,*)
+ & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
+ & " vrand_mat2"
+ do i=1,dimen
+ do j=1,dimen
+ write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
+ & vfric_mat(i,j),afric_mat(i,j),
+ & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
+ enddo
+ enddo
+ endif
+ do i=1,dimen
+ ddt1=0.0d0
+ ddt2=0.0d0
+ do j=1,dimen
+ dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
+ & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
+ ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
+ ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
+ enddo
+ d_t_work_new(i)=ddt1+0.5d0*ddt2
+ d_t_work(i)=ddt1+ddt2
+ enddo
+#endif
+ do j=1,3
+ dc(j,0)=dc_work(j)
+ d_t(j,0)=d_t_work(j)
+ enddo
+ 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)
+ enddo
+ ind=ind+3
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ inres=i+nres
+ do j=1,3
+ dc(j,inres)=dc_work(ind+j)
+ d_t(j,inres)=d_t_work(ind+j)
+ enddo
+ ind=ind+3
+ endif
+ enddo
+ return
+ end
+c--------------------------------------------------------------------------
+ subroutine sd_verlet2
+c Calculating the adjusted velocities for accelerations
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CONTROL'
+ include 'COMMON.VAR'
+ include 'COMMON.MD'
+#ifndef LANG0
+ include 'COMMON.LANGEVIN'
+#else
+ include 'COMMON.LANGEVIN.lang0'
+#endif
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.NAMES'
+ double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
+ common /stochcalc/ stochforcvec
+c
+c Compute the stochastic forces which contribute to velocity change
+c
+ call stochastic_force(stochforcvecV)
+
+#ifndef LANG0
+ do i=1,dimen
+ ddt1=0.0d0
+ ddt2=0.0d0
+ do j=1,dimen
+ ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
+ ddt2=ddt2+vrand_mat1(i,j)*stochforcvec(j)+
+ & vrand_mat2(i,j)*stochforcvecV(j)
+ enddo
+ d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
+ enddo
+#endif
+ do j=1,3
+ d_t(j,0)=d_t_work(j)
+ enddo
+ ind=3
+ do i=nnt,nct-1
+ do j=1,3
+ d_t(j,i)=d_t_work(ind+j)
+ enddo
+ ind=ind+3
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ inres=i+nres
+ do j=1,3
+ d_t(j,inres)=d_t_work(ind+j)
+ enddo
+ ind=ind+3
+ endif
+ enddo
+ return
+ end
+c-----------------------------------------------------------
+ subroutine sd_verlet_ciccotti_setup
+c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
+c version
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+#ifdef MPI
+ include 'mpif.h'
+#endif
+ include 'COMMON.CONTROL'
+ include 'COMMON.VAR'
+ include 'COMMON.MD'
+#ifndef LANG0
+ include 'COMMON.LANGEVIN'
+#else
+ include 'COMMON.LANGEVIN.lang0'
+#endif
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.NAMES'
+ include 'COMMON.TIME1'
+ double precision emgdt(MAXRES6),
+ & pterm,vterm,rho,rhoc,vsig,
+ & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
+ & afric_vec(MAXRES6),prand_vec(MAXRES6),
+ & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
+ logical lprn /.false./
+ double precision zero /1.0d-8/, gdt_radius /0.05d0/
+ double precision ktm
+#ifdef MPI
+ tt0 = MPI_Wtime()
+#else
+ tt0 = tcpu()
+#endif
+c
+c AL 8/17/04 Code adapted from tinker
+c
+c Get the frictional and random terms for stochastic dynamics in the
+c eigenspace of mass-scaled UNRES friction matrix
+c
+ do i = 1, dimen
+ write (iout,*) "i",i," fricgam",fricgam(i)
+ gdt = fricgam(i) * d_time
+c
+c Stochastic dynamics reduces to simple MD for zero friction
+c
+ if (gdt .le. zero) then
+ pfric_vec(i) = 1.0d0
+ vfric_vec(i) = d_time
+ afric_vec(i) = 0.5d0*d_time*d_time
+ prand_vec(i) = afric_vec(i)
+ vrand_vec2(i) = vfric_vec(i)
+c
+c Analytical expressions when friction coefficient is large
+c
+ else
+ egdt = dexp(-gdt)
+ pfric_vec(i) = egdt
+ vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
+ afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
+ prand_vec(i) = afric_vec(i)
+ vrand_vec2(i) = vfric_vec(i)
+c
+c Compute the scaling factors of random terms for the nonzero friction case
+c
+c ktm = 0.5d0*d_time/fricgam(i)
+c psig = dsqrt(ktm*pterm) / fricgam(i)
+c vsig = dsqrt(ktm*vterm)
+c prand_vec(i) = psig*afric_vec(i)
+c vrand_vec2(i) = vsig*vfric_vec(i)
+ end if
+ end do
+ if (lprn) then
+ write (iout,*)
+ & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
+ & " vrand_vec2"
+ do i=1,dimen
+ write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
+ & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
+ enddo
+ endif
+c
+c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
+c
+ call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
+ call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
+ call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
+ call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
+ call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
+#ifdef MPI
+ t_sdsetup=t_sdsetup+MPI_Wtime()
+#else
+ t_sdsetup=t_sdsetup+tcpu()-tt0
+#endif
+ return
+ end
+c-------------------------------------------------------------
+ subroutine sd_verlet1_ciccotti
+c Applying stochastic velocity Verlet algorithm - step 1 to velocities
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+#ifdef MPI
+ include 'mpif.h'
+#endif
+ include 'COMMON.CONTROL'
+ include 'COMMON.VAR'
+ include 'COMMON.MD'
+#ifndef LANG0
+ include 'COMMON.LANGEVIN'
+#else
+ include 'COMMON.LANGEVIN.lang0'
+#endif
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.NAMES'
+ double precision stochforcvec(MAXRES6)
+ common /stochcalc/ stochforcvec
+ logical lprn /.false./
+
+c write (iout,*) "dc_old"
+c do i=0,nres
+c write (iout,'(i5,3f10.5,5x,3f10.5)')
+c & i,(dc_old(j,i),j=1,3),(dc_old(j,i+nres),j=1,3)
+c enddo
+ do j=1,3
+ dc_work(j)=dc_old(j,0)
+ d_t_work(j)=d_t_old(j,0)
+ d_a_work(j)=d_a_old(j,0)
+ enddo
+ ind=3
+ do i=nnt,nct-1
+ do j=1,3
+ dc_work(ind+j)=dc_old(j,i)
+ d_t_work(ind+j)=d_t_old(j,i)
+ d_a_work(ind+j)=d_a_old(j,i)
+ enddo
+ ind=ind+3
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ do j=1,3
+ dc_work(ind+j)=dc_old(j,i+nres)
+ d_t_work(ind+j)=d_t_old(j,i+nres)
+ d_a_work(ind+j)=d_a_old(j,i+nres)
+ enddo
+ ind=ind+3
+ endif
+ enddo
+
+#ifndef LANG0
+ if (lprn) then
+ write (iout,*)
+ & "pfric_mat, vfric_mat, afric_mat, prand_mat, vrand_mat1,",
+ & " vrand_mat2"
+ do i=1,dimen
+ do j=1,dimen
+ write (iout,'(2i5,6e15.5)') i,j,pfric_mat(i,j),
+ & vfric_mat(i,j),afric_mat(i,j),
+ & prand_mat(i,j),vrand_mat1(i,j),vrand_mat2(i,j)
+ enddo
+ enddo
+ endif
+ do i=1,dimen
+ ddt1=0.0d0
+ ddt2=0.0d0
+ do j=1,dimen
+ dc_work(i)=dc_work(i)+vfric_mat(i,j)*d_t_work(j)
+ & +afric_mat(i,j)*d_a_work(j)+prand_mat(i,j)*stochforcvec(j)
+ ddt1=ddt1+pfric_mat(i,j)*d_t_work(j)
+ ddt2=ddt2+vfric_mat(i,j)*d_a_work(j)
+ enddo
+ d_t_work_new(i)=ddt1+0.5d0*ddt2
+ d_t_work(i)=ddt1+ddt2
+ enddo
+#endif
+ do j=1,3
+ dc(j,0)=dc_work(j)
+ d_t(j,0)=d_t_work(j)
+ enddo
+ 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)
+ enddo
+ ind=ind+3
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ inres=i+nres
+ do j=1,3
+ dc(j,inres)=dc_work(ind+j)
+ d_t(j,inres)=d_t_work(ind+j)
+ enddo
+ ind=ind+3
+ endif
+ enddo
+ return
+ end
+c--------------------------------------------------------------------------
+ subroutine sd_verlet2_ciccotti
+c Calculating the adjusted velocities for accelerations
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CONTROL'
+ include 'COMMON.VAR'
+ include 'COMMON.MD'
+#ifndef LANG0
+ include 'COMMON.LANGEVIN'
+#else
+ include 'COMMON.LANGEVIN.lang0'
+#endif
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.NAMES'
+ double precision stochforcvec(MAXRES6),stochforcvecV(MAXRES6)
+ common /stochcalc/ stochforcvec
+c
+c Compute the stochastic forces which contribute to velocity change
+c
+ call stochastic_force(stochforcvecV)
+#ifndef LANG0
+ do i=1,dimen
+ ddt1=0.0d0
+ ddt2=0.0d0
+ do j=1,dimen
+
+ ddt1=ddt1+vfric_mat(i,j)*d_a_work(j)
+c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
+ ddt2=ddt2+vrand_mat2(i,j)*stochforcvec(j)
+ enddo
+ d_t_work(i)=d_t_work_new(i)+0.5d0*ddt1+ddt2
+ enddo
+#endif
+ do j=1,3
+ d_t(j,0)=d_t_work(j)
+ enddo
+ ind=3
+ do i=nnt,nct-1
+ do j=1,3
+ d_t(j,i)=d_t_work(ind+j)
+ enddo
+ ind=ind+3
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ inres=i+nres
+ do j=1,3
+ d_t(j,inres)=d_t_work(ind+j)
+ enddo
+ ind=ind+3
+ endif
+ enddo
+ return
+ end
+#endif
+c------------------------------------------------------
+ double precision function HNose(ek,s,e,pi,Q,t_bath,dimenl)
+ implicit none
+ double precision ek,s,e,pi,Q,t_bath,Rb
+ integer dimenl
+ Rb=0.001986d0
+ HNose=ek+e+pi**2/(2*Q)+dimenl*Rb*t_bath*log(s)
+c print '(6f15.5,i5,a2,2f15.5)',ek,s,e,pi,Q,t_bath,dimenl,"--",
+c & pi**2/(2*Q),dimenl*Rb*t_bath*log(s)
+ return
+ end
+c-----------------------------------------------------------------
+ double precision function HNose_nh(eki,e)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.MD'
+ HNose_nh=eki+e+dimen3*Rb*t_bath*xlogs(1)+qmass(1)*vlogs(1)**2/2
+ do i=2,nnos
+ HNose_nh=HNose_nh+qmass(i)*vlogs(i)**2/2+Rb*t_bath*xlogs(i)
+ enddo
+c write(4,'(5e15.5)')
+c & vlogs(1),xlogs(1),HNose,eki,e
+ return
+ end
+c-----------------------------------------------------------------
+ SUBROUTINE NHCINT(akin,scale,wdti,wdti2,wdti4,wdti8)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.MD'
+ double precision akin,gnkt,dt,aa,gkt,scale
+ double precision wdti(maxyosh),wdti2(maxyosh),
+ & wdti4(maxyosh),wdti8(maxyosh)
+ integer i,iresn,iyosh,inos,nnos1
+
+ dt=d_time
+ nnos1=nnos+1
+ GKT = Rb*t_bath
+ GNKT = dimen3*GKT
+ akin=akin*2
+
+
+C THIS ROUTINE DOES THE NOSE-HOOVER PART OF THE
+C INTEGRATION FROM t=0 TO t=DT/2
+C GET THE TOTAL KINETIC ENERGY
+ SCALE = 1.D0
+c CALL GETKINP(MASS,VX,VY,VZ,AKIN)
+C UPDATE THE FORCES
+ GLOGS(1) = (AKIN - GNKT)/QMASS(1)
+C START THE MULTIPLE TIME STEP PROCEDURE
+ DO IRESN = 1,NRESN
+ DO IYOSH = 1,NYOSH
+C UPDATE THE THERMOSTAT VELOCITIES
+ VLOGS(NNOS) = VLOGS(NNOS) + GLOGS(NNOS)*WDTI4(IYOSH)
+ DO INOS = 1,NNOS-1
+ AA = EXP(-WDTI8(IYOSH)*VLOGS(NNOS1-INOS) )
+ VLOGS(NNOS-INOS) = VLOGS(NNOS-INOS)*AA*AA
+ & + WDTI4(IYOSH)*GLOGS(NNOS-INOS)*AA
+ ENDDO
+C UPDATE THE PARTICLE VELOCITIES
+ AA = EXP(-WDTI2(IYOSH)*VLOGS(1) )
+ SCALE = SCALE*AA
+C UPDATE THE FORCES
+ GLOGS(1) = (SCALE*SCALE*AKIN - GNKT)/QMASS(1)
+C UPDATE THE THERMOSTAT POSITIONS
+ DO INOS = 1,NNOS
+ XLOGS(INOS) = XLOGS(INOS) + VLOGS(INOS)*WDTI2(IYOSH)
+ ENDDO
+C UPDATE THE THERMOSTAT VELOCITIES
+ DO INOS = 1,NNOS-1
+ AA = EXP(-WDTI8(IYOSH)*VLOGS(INOS+1) )
+ VLOGS(INOS) = VLOGS(INOS)*AA*AA
+ & + WDTI4(IYOSH)*GLOGS(INOS)*AA
+ GLOGS(INOS+1) = (QMASS(INOS)*VLOGS(INOS)*VLOGS(INOS)
+ & -GKT)/QMASS(INOS+1)
+ ENDDO
+ VLOGS(NNOS) = VLOGS(NNOS) + GLOGS(NNOS)*WDTI4(IYOSH)
+ ENDDO
+ ENDDO
+C UPDATE THE PARTICLE VELOCITIES
+c outside of this subroutine
+c DO I = 1,N
+c VX(I) = VX(I)*SCALE
+c VY(I) = VY(I)*SCALE
+c VZ(I) = VZ(I)*SCALE
+c ENDDO
+ RETURN
+ END
+c-----------------------------------------------------------------
+ subroutine tnp1_respa_i_step1
+c Applying Nose-Poincare algorithm - step 1 to coordinates
+c JPSJ 70 75 (2001) S. Nose
+c
+c d_t is not updated here
+c
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CONTROL'
+ include 'COMMON.VAR'
+ include 'COMMON.MD'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.NAMES'
+ double precision adt,adt2,tmp
+
+ tmp=1+pi_np/(2*Q_np)*0.5*d_time
+ s12_np=s_np*tmp**2
+ pistar=pi_np/tmp
+ s12_dt=d_time/s12_np
+ d_time_s12=d_time*0.5*s12_np
+
+ do j=1,3
+ d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s12
+ dc(j,0)=dc_old(j,0)+d_t_new(j,0)*s12_dt
+ enddo
+ do i=nnt,nct-1
+ do j=1,3
+ d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s12
+ dc(j,i)=dc_old(j,i)+d_t_new(j,i)*s12_dt
+ enddo
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ inres=i+nres
+ do j=1,3
+ d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s12
+ dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*s12_dt
+ enddo
+ endif
+ enddo
+ return
+ end
+c---------------------------------------------------------------------
+ subroutine tnp1_respa_i_step2
+c Step 2 of the velocity Verlet algorithm: update velocities
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CONTROL'
+ include 'COMMON.VAR'
+ include 'COMMON.MD'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.NAMES'
+
+ double precision d_time_s12
+
+ do i=0,2*nres
+ do j=1,3
+ d_t(j,i)=d_t_new(j,i)
+ enddo
+ enddo
+
+ call kinetic(EK)
+ EK=EK/s12_np**2
+
+ d_time_s12=0.5d0*s12_np*d_time
+
+ do j=1,3
+ d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s12
+ enddo
+ do i=nnt,nct-1
+ do j=1,3
+ d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s12
+ enddo
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ inres=i+nres
+ do j=1,3
+ d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s12
+ enddo
+ endif
+ enddo
+
+ pistar=pistar+(EK-0.5*(E_old+potE)
+ & -dimen3*Rb*t_bath*log(s12_np)+Csplit-dimen3*Rb*t_bath)*d_time
+ tmp=1+pistar/(2*Q_np)*0.5*d_time
+ s_np=s12_np*tmp**2
+ pi_np=pistar/tmp
+
+ return
+ end
+c-------------------------------------------------------
+
+ subroutine tnp1_step1
+c Applying Nose-Poincare algorithm - step 1 to coordinates
+c JPSJ 70 75 (2001) S. Nose
+c
+c d_t is not updated here
+c
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CONTROL'
+ include 'COMMON.VAR'
+ include 'COMMON.MD'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.NAMES'
+ double precision adt,adt2,tmp
+
+ tmp=1+pi_np/(2*Q_np)*0.5*d_time
+ s12_np=s_np*tmp**2
+ pistar=pi_np/tmp
+ s12_dt=d_time/s12_np
+ d_time_s12=d_time*0.5*s12_np
+
+ do j=1,3
+ d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s12
+ dc(j,0)=dc_old(j,0)+d_t_new(j,0)*s12_dt
+ enddo
+ do i=nnt,nct-1
+ do j=1,3
+ d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s12
+ dc(j,i)=dc_old(j,i)+d_t_new(j,i)*s12_dt
+ enddo
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ inres=i+nres
+ do j=1,3
+ d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s12
+ dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*s12_dt
+ enddo
+ endif
+ enddo
+ return
+ end
+c---------------------------------------------------------------------
+ subroutine tnp1_step2
+c Step 2 of the velocity Verlet algorithm: update velocities
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CONTROL'
+ include 'COMMON.VAR'
+ include 'COMMON.MD'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.NAMES'
+
+ double precision d_time_s12
+
+ do i=0,2*nres
+ do j=1,3
+ d_t(j,i)=d_t_new(j,i)
+ enddo
+ enddo
+
+ call kinetic(EK)
+ EK=EK/s12_np**2
+
+ d_time_s12=0.5d0*s12_np*d_time
+
+ do j=1,3
+ d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s12
+ enddo
+ do i=nnt,nct-1
+ do j=1,3
+ d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s12
+ enddo
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ inres=i+nres
+ do j=1,3
+ d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s12
+ enddo
+ endif
+ enddo
+
+cd write(iout,*) 'pistar',pistar,EK,E_old,potE,s12_np
+ pistar=pistar+(EK-0.5*(E_old+potE)
+ & -dimen3*Rb*t_bath*log(s12_np)+H0-dimen3*Rb*t_bath)*d_time
+ tmp=1+pistar/(2*Q_np)*0.5*d_time
+ s_np=s12_np*tmp**2
+ pi_np=pistar/tmp
+
+ return
+ end
+
+c-----------------------------------------------------------------
+ subroutine tnp_respa_i_step1
+c Applying Nose-Poincare algorithm - step 1 to coordinates
+c J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird
+c
+c d_t is not updated here, it is destroyed
+c
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CONTROL'
+ include 'COMMON.VAR'
+ include 'COMMON.MD'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.NAMES'
+ double precision C_np,d_time_s,tmp,d_time_ss
+
+ d_time_s=d_time*0.5*s_np
+ct2 d_time_s=d_time*0.5*s12_np
+
+ do j=1,3
+ d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s
+ enddo
+ do i=nnt,nct-1
+ do j=1,3
+ d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s
+ enddo
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ inres=i+nres
+ do j=1,3
+ d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s
+ enddo
+ endif
+ enddo
+
+ do i=0,2*nres
+ do j=1,3
+ d_t(j,i)=d_t_new(j,i)
+ enddo
+ enddo
+
+ call kinetic(EK)
+ EK=EK/s_np**2
+
+ C_np=0.5*d_time*(dimen3*Rb*t_bath*(1.0+log(s_np))-EK+potE-Csplit)
+ & -pi_np
+
+ pistar=-2.0*C_np/(1.0+sqrt(1.0-C_np*d_time/Q_np))
+ tmp=0.5*d_time*pistar/Q_np
+ s12_np=s_np*(1.0+tmp)/(1.0-tmp)
+
+ d_time_ss=0.5*d_time*(1.0/s12_np+1.0/s_np)
+ct2 d_time_ss=d_time/s12_np
+c d_time_ss=0.5*d_time*(1.0/sold_np+1.0/s_np)
+
+ do j=1,3
+ dc(j,0)=dc_old(j,0)+d_t_new(j,0)*d_time_ss
+ enddo
+ do i=nnt,nct-1
+ do j=1,3
+ dc(j,i)=dc_old(j,i)+d_t_new(j,i)*d_time_ss
+ enddo
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ inres=i+nres
+ do j=1,3
+ dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*d_time_ss
+ enddo
+ endif
+ enddo
+
+ return
+ end
+c---------------------------------------------------------------------
+
+ subroutine tnp_respa_i_step2
+c Step 2 of the velocity Verlet algorithm: update velocities
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CONTROL'
+ include 'COMMON.VAR'
+ include 'COMMON.MD'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.NAMES'
+
+ double precision d_time_s
+
+ EK=EK*(s_np/s12_np)**2
+ HNose1=Hnose(EK,s12_np,potE,pistar,Q_np,t_bath,dimen3)
+ pi_np=pistar+0.5*d_time*(2*EK-dimen3*Rb*t_bath
+ & -HNose1+Csplit)
+
+cr print '(a,5f)','i_step2',EK,potE,HNose1,pi_np,E_long
+ d_time_s=d_time*0.5*s12_np
+c d_time_s=d_time*0.5*s_np
+
+ do j=1,3
+ d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s
+ enddo
+ do i=nnt,nct-1
+ do j=1,3
+ d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s
+ enddo
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ inres=i+nres
+ do j=1,3
+ d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s
+ enddo
+ endif
+ enddo
+
+ s_np=s12_np
+
+ return
+ end
+c-----------------------------------------------------------------
+ subroutine tnp_respa_step1
+c Applying Nose-Poincare algorithm - step 1 to vel for RESPA
+c J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird
+c
+c d_t is not updated here, it is destroyed
+c
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CONTROL'
+ include 'COMMON.VAR'
+ include 'COMMON.MD'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.NAMES'
+ double precision C_np,d_time_s,tmp,d_time_ss
+ double precision energia(0:n_ene)
+
+ d_time_s=d_time*0.5*s_np
+
+ do j=1,3
+ d_t_old(j,0)=d_t_old(j,0)+d_a(j,0)*d_time_s
+ enddo
+ do i=nnt,nct-1
+ do j=1,3
+ d_t_old(j,i)=d_t_old(j,i)+d_a(j,i)*d_time_s
+ enddo
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ inres=i+nres
+ do j=1,3
+ d_t_old(j,inres)=d_t_old(j,inres)+d_a(j,inres)*d_time_s
+ enddo
+ endif
+ enddo
+
+
+c C_np=0.5*d_time*(dimen3*Rb*t_bath*(1.0+log(s_np))-EK+potE-H0)
+c & -pi_np
+c
+c pistar=-2.0*C_np/(1.0+sqrt(1.0-C_np*d_time/Q_np))
+c tmp=0.5*d_time*pistar/Q_np
+c s12_np=s_np*(1.0+tmp)/(1.0-tmp)
+c write(iout,*) 'tnp_respa_step1',s_np,s12_np,EK,potE,C_np,pistar,tmp
+
+ct1 pi_np=pistar
+c sold_np=s_np
+c s_np=s12_np
+
+c-------------------------------------
+c test of reviewer's comment
+ pi_np=pi_np-0.5*d_time*(E_long+Csplit-H0)
+cr print '(a,3f)','1 pi_np,s_np',pi_np,s_np,E_long
+c-------------------------------------
+
+ return
+ end
+c---------------------------------------------------------------------
+ subroutine tnp_respa_step2
+c Step 2 of the velocity Verlet algorithm: update velocities for RESPA
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CONTROL'
+ include 'COMMON.VAR'
+ include 'COMMON.MD'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.NAMES'
+
+ double precision d_time_s
+
+ct1 s12_np=s_np
+ct2 pistar=pi_np
+
+ct call kinetic(EK)
+ct HNose1=Hnose(EK,s12_np,potE,pistar,Q_np,t_bath,dimen3)
+ct pi_np=pistar+0.5*d_time*(2*EK-dimen3*Rb*t_bath)
+ct & -0.5*d_time*(HNose1-H0)
+
+c-------------------------------------
+c test of reviewer's comment
+ pi_np=pi_np-0.5*d_time*(E_long+Csplit-H0)
+cr print '(a,3f)','2 pi_np,s_np',pi_np,s_np,E_long
+c-------------------------------------
+ d_time_s=d_time*0.5*s_np
+
+ do j=1,3
+ d_t_old(j,0)=d_t_old(j,0)+d_a(j,0)*d_time_s
+ enddo
+ do i=nnt,nct-1
+ do j=1,3
+ d_t_old(j,i)=d_t_old(j,i)+d_a(j,i)*d_time_s
+ enddo
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ inres=i+nres
+ do j=1,3
+ d_t_old(j,inres)=d_t_old(j,inres)+d_a(j,inres)*d_time_s
+ enddo
+ endif
+ enddo
+
+cd s_np=s12_np
+
+ return
+ end
+c---------------------------------------------------------------------
+ subroutine tnp_step1
+c Applying Nose-Poincare algorithm - step 1 to coordinates
+c J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird
+c
+c d_t is not updated here, it is destroyed
+c
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CONTROL'
+ include 'COMMON.VAR'
+ include 'COMMON.MD'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.NAMES'
+ double precision C_np,d_time_s,tmp,d_time_ss
+
+ d_time_s=d_time*0.5*s_np
+
+ do j=1,3
+ d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s
+ enddo
+ do i=nnt,nct-1
+ do j=1,3
+ d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s
+ enddo
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ inres=i+nres
+ do j=1,3
+ d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s
+ enddo
+ endif
+ enddo
+
+ do i=0,2*nres
+ do j=1,3
+ d_t(j,i)=d_t_new(j,i)
+ enddo
+ enddo
+
+ call kinetic(EK)
+ EK=EK/s_np**2
+
+ C_np=0.5*d_time*(dimen3*Rb*t_bath*(1.0+log(s_np))-EK+potE-H0)
+ & -pi_np
+
+ pistar=-2.0*C_np/(1.0+sqrt(1.0-C_np*d_time/Q_np))
+ tmp=0.5*d_time*pistar/Q_np
+ s12_np=s_np*(1.0+tmp)/(1.0-tmp)
+c write(iout,*) 'tnp_step1',s_np,s12_np,EK,potE,C_np,pistar,tmp
+
+ d_time_ss=0.5*d_time*(1.0/s12_np+1.0/s_np)
+
+ do j=1,3
+ dc(j,0)=dc_old(j,0)+d_t_new(j,0)*d_time_ss
+ enddo
+ do i=nnt,nct-1
+ do j=1,3
+ dc(j,i)=dc_old(j,i)+d_t_new(j,i)*d_time_ss
+ enddo
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ inres=i+nres
+ do j=1,3
+ dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*d_time_ss
+ enddo
+ endif
+ enddo
+
+ return
+ end
+c-----------------------------------------------------------------
+ subroutine tnp_step2
+c Step 2 of the velocity Verlet algorithm: update velocities
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CONTROL'
+ include 'COMMON.VAR'
+ include 'COMMON.MD'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.GEO'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.NAMES'
+
+ double precision d_time_s
+
+ EK=EK*(s_np/s12_np)**2
+ HNose1=Hnose(EK,s12_np,potE,pistar,Q_np,t_bath,dimen3)
+ pi_np=pistar+0.5*d_time*(2*EK-dimen3*Rb*t_bath)
+ & -0.5*d_time*(HNose1-H0)
+
+cd write(iout,'(a,4f)') 'mmm',EK,potE,HNose1,pi_np
+ d_time_s=d_time*0.5*s12_np
+
+ do j=1,3
+ d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s
+ enddo
+ do i=nnt,nct-1
+ do j=1,3
+ d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s
+ enddo
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.10) then
+ inres=i+nres
+ do j=1,3
+ d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s
+ enddo
+ endif
+ enddo
+
+ s_np=s12_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
projects / unres.git / blobdiff