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)
11 integer ierror,ierrcode
13 include 'COMMON.SETUP'
14 include 'COMMON.CONTROL'
18 include 'COMMON.LANGEVIN'
20 include 'COMMON.LANGEVIN.lang0'
22 include 'COMMON.CHAIN'
23 include 'COMMON.DERIV'
25 include 'COMMON.LOCAL'
26 include 'COMMON.INTERACT'
27 include 'COMMON.IOUNITS'
28 include 'COMMON.NAMES'
29 include 'COMMON.TIME1'
31 double precision vcm(3),incr(3)
32 double precision cm(3),L(3)
33 integer ilen,count,rstcount
36 integer maxcount_scale /20/
38 double precision stochforcvec(MAXRES6)
39 common /stochcalc/ stochforcvec
47 else if (lang.eq.2 .or. lang.eq.3) then
49 call stochastic_force(stochforcvec)
52 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
54 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
61 icount_scale=icount_scale+1
62 if (icount_scale.gt.maxcount_scale) then
64 & "ERROR: too many attempts at scaling down the time step. ",
65 & "amax=",amax,"epdrift=",epdrift,
66 & "damax=",damax,"edriftmax=",edriftmax,
70 call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
74 c First step of the velocity Verlet algorithm
79 else if (lang.eq.3) then
81 call sd_verlet1_ciccotti
83 else if (lang.eq.1) then
88 c Build the chain from the newly calculated coordinates
90 if (rattle) call rattle1
92 if (large.and. mod(itime,ntwe).eq.0) then
93 write (iout,*) "Cartesian and internal coordinates: step 1"
98 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
99 & (dc(j,i+nres),j=1,3)
101 write (iout,*) "Accelerations"
103 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
104 & (d_a(j,i+nres),j=1,3)
106 write (iout,*) "Velocities, step 1"
108 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
109 & (d_t(j,i+nres),j=1,3)
118 c Calculate energy and forces
120 call etotal(potEcomp)
123 t_etotal=t_etotal+MPI_Wtime()-tt0
125 t_etotal=t_etotal+tcpu()-tt0
128 potE=potEcomp(0)-potEcomp(20)
130 c Get the new accelerations
133 t_enegrad=t_enegrad+MPI_Wtime()-tt0
135 t_enegrad=t_enegrad+tcpu()-tt0
137 c Determine maximum acceleration and scale down the timestep if needed
139 amax=amax/(itime_scal+1)**2
140 call predict_edrift(epdrift)
141 if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
142 c Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
144 ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax))
146 itime_scal=itime_scal+ifac_time
147 c fac_time=dmin1(damax/amax,0.5d0)
148 fac_time=0.5d0**ifac_time
149 d_time=d_time*fac_time
150 if (lang.eq.2 .or. lang.eq.3) then
152 c write (iout,*) "Calling sd_verlet_setup: 1"
153 c Rescale the stochastic forces and recalculate or restore
154 c the matrices of tinker integrator
155 if (itime_scal.gt.maxflag_stoch) then
156 if (large) write (iout,'(a,i5,a)')
157 & "Calculate matrices for stochastic step;",
158 & " itime_scal ",itime_scal
160 call sd_verlet_p_setup
162 call sd_verlet_ciccotti_setup
164 write (iout,'(2a,i3,a,i3,1h.)')
165 & "Warning: cannot store matrices for stochastic",
166 & " integration because the index",itime_scal,
167 & " is greater than",maxflag_stoch
168 write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",
169 & " integration Langevin algorithm for better efficiency."
170 else if (flag_stoch(itime_scal)) then
171 if (large) write (iout,'(a,i5,a,l1)')
172 & "Restore matrices for stochastic step; itime_scal ",
173 & itime_scal," flag ",flag_stoch(itime_scal)
176 pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
177 afric_mat(i,j)=afric0_mat(i,j,itime_scal)
178 vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
179 prand_mat(i,j)=prand0_mat(i,j,itime_scal)
180 vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
181 vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
185 if (large) write (iout,'(2a,i5,a,l1)')
186 & "Calculate & store matrices for stochastic step;",
187 & " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
189 call sd_verlet_p_setup
191 call sd_verlet_ciccotti_setup
193 flag_stoch(ifac_time)=.true.
196 pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
197 afric0_mat(i,j,itime_scal)=afric_mat(i,j)
198 vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
199 prand0_mat(i,j,itime_scal)=prand_mat(i,j)
200 vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
201 vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
205 fac_time=1.0d0/dsqrt(fac_time)
207 stochforcvec(i)=fac_time*stochforcvec(i)
210 else if (lang.eq.1) then
211 c Rescale the accelerations due to stochastic forces
212 fac_time=1.0d0/dsqrt(fac_time)
214 d_as_work(i)=d_as_work(i)*fac_time
217 if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)')
218 & "itime",itime," Timestep scaled down to ",
219 & d_time," ifac_time",ifac_time," itime_scal",itime_scal
221 c Second step of the velocity Verlet algorithm
226 else if (lang.eq.3) then
228 call sd_verlet2_ciccotti
230 else if (lang.eq.1) then
235 if (rattle) call rattle2
237 if (d_time.ne.d_time0) then
240 if (lang.eq.2 .or. lang.eq.3) then
241 if (large) write (iout,'(a)')
242 & "Restore original matrices for stochastic step"
243 c write (iout,*) "Calling sd_verlet_setup: 2"
244 c Restore the matrices of tinker integrator if the time step has been restored
247 pfric_mat(i,j)=pfric0_mat(i,j,0)
248 afric_mat(i,j)=afric0_mat(i,j,0)
249 vfric_mat(i,j)=vfric0_mat(i,j,0)
250 prand_mat(i,j)=prand0_mat(i,j,0)
251 vrand_mat1(i,j)=vrand0_mat1(i,j,0)
252 vrand_mat2(i,j)=vrand0_mat2(i,j,0)
261 c Calculate the kinetic and the total energy and the kinetic temperature
266 c write (iout,*) "step",itime," EK",EK," EK1",EK1
268 c Couple the system to Berendsen bath if needed
269 if (tbf .and. lang.eq.0) then
272 kinetic_T=2.0d0/(dimen3*Rb)*EK
273 c Backup the coordinates, velocities, and accelerations
277 d_t_old(j,i)=d_t(j,i)
278 d_a_old(j,i)=d_a(j,i)
282 if (mod(itime,ntwe).eq.0 .and. large) then
283 write (iout,*) "Velocities, step 2"
285 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
286 & (d_t(j,i+nres),j=1,3)