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' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' include 'COMMON.LANGEVIN' 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 energia(0:n_ene),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, energia double precision stochforcvec(MAXRES6) common /stochcalc/ stochforcvec integer itime logical scale c scale=.true. icount_scale=0 if (lang.eq.1) then call sddir_precalc else if (lang.eq.2 .or. lang.eq.3) then call stochastic_force(stochforcvec) 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 stop endif c First step of the velocity Verlet algorithm if (lang.eq.2) then call sd_verlet1 else if (lang.eq.3) then call sd_verlet1_ciccotti else if (lang.eq.1) then call sddir_verlet1 else call verlet1 endif c Build the chain from the newly calculated coordinates call chainbuild_cart if (rattle) call rattle1 if (large.and. mod(itime,ntwe).eq.0) then write (iout,*) "Cartesian and internal coordinates: step 1" call cartprint call intout 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 tt0 = tcpu() c Calculate energy and forces call zerograd call etotal(energia) potE=energia(0)-energia(20) call cartgrad c Get the new accelerations call lagrangian t_enegrad=t_enegrad+tcpu()-tt0 c Determine maximum acceleration and scale down the timestep if needed call max_accel call predict_edrift(epdrift) if (amax.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 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,dimen do j=1,dimen 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,dimen do j=1,dimen 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,dimen stochforcvec(i)=fac_time*stochforcvec(i) enddo else if (lang.eq.1) then c Rescale the accelerations due to stochastic forces fac_time=1.0d0/dsqrt(fac_time) do i=1,dimen 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 call sd_verlet2 else if (lang.eq.3) then call sd_verlet2_ciccotti else if (lang.eq.1) then call sddir_verlet2 else call verlet2 endif if (rattle) call rattle2 totT=totT+d_time if (d_time.ne.d_time0) then d_time=d_time0 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,dimen do j=1,dimen 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 scale=.false. endif enddo c Calculate the kinetic and the total energy and the kinetic temperature 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/(dimen*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) d_t_old(j,i)=d_t(j,i) d_a_old(j,i)=d_a(j,i) enddo enddo if (mod(itime,ntwe).eq.0 .and. large) then 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 return end