1 c-------------------------------------------------------------------------------
2 subroutine velverlet_step(itime)
3 c-------------------------------------------------------------------------------
4 c Perform a single velocity Verlet step; the time step can be rescaled if
5 c increments in accelerations exceed the threshold
6 c-------------------------------------------------------------------------------
7 implicit real*8 (a-h,o-z)
9 include 'COMMON.CONTROL'
12 include 'COMMON.LANGEVIN'
13 include 'COMMON.CHAIN'
14 include 'COMMON.DERIV'
16 include 'COMMON.LOCAL'
17 include 'COMMON.INTERACT'
18 include 'COMMON.IOUNITS'
19 include 'COMMON.NAMES'
20 include 'COMMON.TIME1'
22 double precision energia(0:n_ene),vcm(3),incr(3)
23 double precision cm(3),L(3)
24 integer ilen,count,rstcount
27 integer maxcount_scale /20/
28 common /gucio/ cm, energia
29 double precision stochforcvec(MAXRES6)
30 common /stochcalc/ stochforcvec
38 else if (lang.eq.2 .or. lang.eq.3) then
39 call stochastic_force(stochforcvec)
43 icount_scale=icount_scale+1
44 if (icount_scale.gt.maxcount_scale) then
46 & "ERROR: too many attempts at scaling down the time step. ",
47 & "amax=",amax,"epdrift=",epdrift,
48 & "damax=",damax,"edriftmax=",edriftmax,
52 c First step of the velocity Verlet algorithm
55 else if (lang.eq.3) then
56 call sd_verlet1_ciccotti
57 else if (lang.eq.1) then
62 c Build the chain from the newly calculated coordinates
64 if (rattle) call rattle1
65 if (large.and. mod(itime,ntwe).eq.0) then
66 write (iout,*) "Cartesian and internal coordinates: step 1"
69 write (iout,*) "Accelerations"
71 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
72 & (d_a(j,i+nres),j=1,3)
74 write (iout,*) "Velocities, step 1"
76 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
77 & (d_t(j,i+nres),j=1,3)
81 c Calculate energy and forces
84 potE=energia(0)-energia(20)
86 c Get the new accelerations
88 t_enegrad=t_enegrad+tcpu()-tt0
89 c Determine maximum acceleration and scale down the timestep if needed
91 call predict_edrift(epdrift)
92 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
93 c Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
95 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax))
97 itime_scal=itime_scal+ifac_time
98 c fac_time=dmin1(damax/amax,0.5d0)
99 fac_time=0.5d0**ifac_time
100 d_time=d_time*fac_time
101 if (lang.eq.2 .or. lang.eq.3) then
102 c write (iout,*) "Calling sd_verlet_setup: 1"
103 c Rescale the stochastic forces and recalculate or restore
104 c the matrices of tinker integrator
105 if (itime_scal.gt.maxflag_stoch) then
106 if (large) write (iout,'(a,i5,a)')
107 & "Calculate matrices for stochastic step;",
108 & " itime_scal ",itime_scal
110 call sd_verlet_p_setup
112 call sd_verlet_ciccotti_setup
114 write (iout,'(2a,i3,a,i3,1h.)')
115 & "Warning: cannot store matrices for stochastic",
116 & " integration because the index",itime_scal,
117 & " is greater than",maxflag_stoch
118 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",
119 & " integration Langevin algorithm for better efficiency."
120 else if (flag_stoch(itime_scal)) then
121 if (large) write (iout,'(a,i5,a,l1)')
122 & "Restore matrices for stochastic step; itime_scal ",
123 & itime_scal," flag ",flag_stoch(itime_scal)
126 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
127 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
128 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
129 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
130 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
131 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
135 if (large) write (iout,'(2a,i5,a,l1)')
136 & "Calculate & store matrices for stochastic step;",
137 & " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
139 call sd_verlet_p_setup
141 call sd_verlet_ciccotti_setup
143 flag_stoch(ifac_time)=.true.
146 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
147 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
148 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
149 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
150 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
151 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
155 fac_time=1.0d0/dsqrt(fac_time)
157 stochforcvec(i)=fac_time*stochforcvec(i)
159 else if (lang.eq.1) then
160 c Rescale the accelerations due to stochastic forces
161 fac_time=1.0d0/dsqrt(fac_time)
163 d_as_work(i)=d_as_work(i)*fac_time
166 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)')
167 & "itime",itime," Timestep scaled down to ",
168 & d_time," ifac_time",ifac_time," itime_scal",itime_scal
170 c Second step of the velocity Verlet algorithm
173 else if (lang.eq.3) then
174 call sd_verlet2_ciccotti
175 else if (lang.eq.1) then
180 if (rattle) call rattle2
182 if (d_time.ne.d_time0) then
184 if (lang.eq.2 .or. lang.eq.3) then
185 if (large) write (iout,'(a)')
186 & "Restore original matrices for stochastic step"
187 c write (iout,*) "Calling sd_verlet_setup: 2"
188 c Restore the matrices of tinker integrator if the time step has been restored
191 pfric_mat(i,j)=pfric0_mat(i,j,0)
192 afric_mat(i,j)=afric0_mat(i,j,0)
193 vfric_mat(i,j)=vfric0_mat(i,j,0)
194 prand_mat(i,j)=prand0_mat(i,j,0)
195 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
196 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
204 c Calculate the kinetic and the total energy and the kinetic temperature
209 c write (iout,*) "step",itime," EK",EK," EK1",EK1
211 c Couple the system to Berendsen bath if needed
212 if (tbf .and. lang.eq.0) then
215 kinetic_T=2.0d0/(dimen*Rb)*EK
216 c Backup the coordinates, velocities, and accelerations
220 d_t_old(j,i)=d_t(j,i)
221 d_a_old(j,i)=d_a(j,i)
224 if (mod(itime,ntwe).eq.0 .and. large) then
225 write (iout,*) "Velocities, step 2"
227 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
228 & (d_t(j,i+nres),j=1,3)