+++ /dev/null
-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
projects / unres.git / blobdiff