1 c-------------------------------------------------------------------------------
2 subroutine RESPA_step(itime)
3 c-------------------------------------------------------------------------------
4 c Perform a single RESPA step.
5 c-------------------------------------------------------------------------------
6 implicit real*8 (a-h,o-z)
10 integer IERROR,ERRCODE
12 include 'COMMON.SETUP'
13 include 'COMMON.CONTROL'
17 include 'COMMON.LANGEVIN'
19 include 'COMMON.LANGEVIN.lang0'
21 include 'COMMON.CHAIN'
22 include 'COMMON.DERIV'
24 include 'COMMON.LOCAL'
25 include 'COMMON.INTERACT'
26 include 'COMMON.IOUNITS'
27 include 'COMMON.NAMES'
28 include 'COMMON.TIME1'
29 double precision energia_short(0:n_ene),
30 & energia_long(0:n_ene)
31 double precision cm(3),L(3),vcm(3),incr(3)
32 double precision dc_old0(3,0:maxres2),d_t_old0(3,0:maxres2),
33 & d_a_old0(3,0:maxres2)
34 integer ilen,count,rstcount
37 integer maxcount_scale /10/
39 double precision stochforcvec(MAXRES6)
40 common /stochcalc/ stochforcvec
43 common /cipiszcze/ itt
46 if (large.and. mod(itime,ntwe).eq.0) then
47 write (iout,*) "***************** RESPA itime",itime
48 write (iout,*) "Cartesian and internal coordinates: step 0"
50 call pdbout(0.0d0,"cipiszcze",iout)
52 write (iout,*) "Accelerations from long-range forces"
54 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
55 & (d_a(j,i+nres),j=1,3)
57 write (iout,*) "Velocities, step 0"
59 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
60 & (d_t(j,i+nres),j=1,3)
65 c Perform the initial RESPA step (increment velocities)
66 c write (iout,*) "*********************** RESPA ini"
69 if (mod(itime,ntwe).eq.0 .and. large) then
70 write (iout,*) "Velocities, end"
72 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
73 & (d_t(j,i+nres),j=1,3)
77 c Compute the short-range forces
83 C 7/2/2009 commented out
85 c call etotal_short(energia_short)
88 C 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
91 d_a(j,i)=d_a_short(j,i)
95 if (large.and. mod(itime,ntwe).eq.0) then
96 write (iout,*) "energia_short",energia_short(0)
97 write (iout,*) "Accelerations from short-range forces"
99 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
100 & (d_a(j,i+nres),j=1,3)
105 t_enegrad=t_enegrad+MPI_Wtime()-tt0
107 t_enegrad=t_enegrad+tcpu()-tt0
112 d_t_old(j,i)=d_t(j,i)
113 d_a_old(j,i)=d_a(j,i)
116 c 6/30/08 A-MTS: attempt at increasing the split number
119 dc_old0(j,i)=dc_old(j,i)
120 d_t_old0(j,i)=d_t_old(j,i)
121 d_a_old0(j,i)=d_a_old(j,i)
124 if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
125 if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
132 c write (iout,*) "itime",itime," ntime_split",ntime_split
133 c Split the time step
134 d_time=d_time0/ntime_split
135 c Perform the short-range RESPA steps (velocity Verlet increments of
136 c positions and velocities using short-range forces)
137 c write (iout,*) "*********************** RESPA split"
138 do itsplit=1,ntime_split
141 else if (lang.eq.2 .or. lang.eq.3) then
143 call stochastic_force(stochforcvec)
146 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
148 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
153 c First step of the velocity Verlet algorithm
158 else if (lang.eq.3) then
160 call sd_verlet1_ciccotti
162 else if (lang.eq.1) then
167 c Build the chain from the newly calculated coordinates
169 if (rattle) call rattle1
171 if (large.and. mod(itime,ntwe).eq.0) then
172 write (iout,*) "***** ITSPLIT",itsplit
173 write (iout,*) "Cartesian and internal coordinates: step 1"
174 call pdbout(0.0d0,"cipiszcze",iout)
177 write (iout,*) "Velocities, step 1"
179 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
180 & (d_t(j,i+nres),j=1,3)
189 c Calculate energy and forces
191 call etotal_short(energia_short)
194 t_eshort=t_eshort+MPI_Wtime()-tt0
196 t_eshort=t_eshort+tcpu()-tt0
200 c Get the new accelerations
202 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
205 d_a_short(j,i)=d_a(j,i)
209 if (large.and. mod(itime,ntwe).eq.0) then
210 write (iout,*)"energia_short",energia_short(0)
211 write (iout,*) "Accelerations from short-range forces"
213 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
214 & (d_a(j,i+nres),j=1,3)
219 c Determine maximum acceleration and scale down the timestep if needed
221 amax=amax/ntime_split**2
222 call predict_edrift(epdrift)
223 if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0)
224 & write (iout,*) "amax",amax," damax",damax,
225 & " epdrift",epdrift," epdriftmax",epdriftmax
226 c Exit loop and try with increased split number if the change of
227 c acceleration is too big
228 if (amax.gt.damax .or. epdrift.gt.edriftmax) then
229 if (ntime_split.lt.maxtime_split) then
231 ntime_split=ntime_split*2
234 dc_old(j,i)=dc_old0(j,i)
235 d_t_old(j,i)=d_t_old0(j,i)
236 d_a_old(j,i)=d_a_old0(j,i)
239 write (iout,*) "acceleration/energy drift too large",amax,
240 & epdrift," split increased to ",ntime_split," itime",itime,
245 & "Uh-hu. Bumpy landscape. Maximum splitting number",
247 & " already reached!!! Trying to carry on!"
251 t_enegrad=t_enegrad+MPI_Wtime()-tt0
253 t_enegrad=t_enegrad+tcpu()-tt0
255 c Second step of the velocity Verlet algorithm
260 else if (lang.eq.3) then
262 call sd_verlet2_ciccotti
264 else if (lang.eq.1) then
269 if (rattle) call rattle2
270 c Backup the coordinates, velocities, and accelerations
274 d_t_old(j,i)=d_t(j,i)
275 d_a_old(j,i)=d_a(j,i)
282 c Restore the time step
284 c Compute long-range forces
291 call etotal_long(energia_long)
294 t_elong=t_elong+MPI_Wtime()-tt0
296 t_elong=t_elong+tcpu()-tt0
302 t_enegrad=t_enegrad+MPI_Wtime()-tt0
304 t_enegrad=t_enegrad+tcpu()-tt0
306 c Compute accelerations from long-range forces
308 if (large.and. mod(itime,ntwe).eq.0) then
309 write (iout,*) "energia_long",energia_long(0)
310 write (iout,*) "Cartesian and internal coordinates: step 2"
312 call pdbout(0.0d0,"cipiszcze",iout)
314 write (iout,*) "Accelerations from long-range forces"
316 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
317 & (d_a(j,i+nres),j=1,3)
319 write (iout,*) "Velocities, step 2"
321 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
322 & (d_t(j,i+nres),j=1,3)
326 c Compute the final RESPA step (increment velocities)
327 c write (iout,*) "*********************** RESPA fin"
329 c Compute the complete potential energy
331 potEcomp(i)=energia_short(i)+energia_long(i)
333 potE=potEcomp(0)-potEcomp(20)
334 c potE=energia_short(0)+energia_long(0)
336 c Calculate the kinetic and the total energy and the kinetic temperature
339 c Couple the system to Berendsen bath if needed
340 if (tbf .and. lang.eq.0) then
343 kinetic_T=2.0d0/(dimen3*Rb)*EK
344 c Backup the coordinates, velocities, and accelerations
346 if (mod(itime,ntwe).eq.0 .and. large) then
347 write (iout,*) "Velocities, end"
349 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
350 & (d_t(j,i+nres),j=1,3)