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 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 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 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,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 #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,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 #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 call verlet2 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,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 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/(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) 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 .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 endif return end