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)
8 include 'COMMON.CONTROL'
11 include 'COMMON.LANGEVIN'
12 include 'COMMON.CHAIN'
13 include 'COMMON.DERIV'
15 include 'COMMON.LOCAL'
16 include 'COMMON.INTERACT'
17 include 'COMMON.IOUNITS'
18 include 'COMMON.NAMES'
19 include 'COMMON.TIME1'
20 double precision energia(0:n_ene),energia_short(0:n_ene),
21 & energia_long(0:n_ene)
22 double precision cm(3),L(3),vcm(3),incr(3)
23 integer ilen,count,rstcount
26 integer maxcount_scale /10/
28 double precision stochforcvec(MAXRES6)
29 common /stochcalc/ stochforcvec
32 if (large.and. mod(itime,ntwe).eq.0) then
33 write (iout,*) "***************** RESPA itime",itime
34 write (iout,*) "Cartesian and internal coordinates: step 0"
37 write (iout,*) "Accelerations"
39 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
40 & (d_a(j,i+nres),j=1,3)
42 write (iout,*) "Velocities, step 0"
44 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
45 & (d_t(j,i+nres),j=1,3)
49 c Perform the initial RESPA step (increment velocities)
50 c write (iout,*) "*********************** RESPA ini"
52 if (mod(itime,ntwe).eq.0 .and. large) then
53 write (iout,*) "Velocities, end"
55 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
56 & (d_t(j,i+nres),j=1,3)
59 c Compute the short-range forces
61 call etotal_short(energia_short)
73 d_time=d_time/ntime_split
74 c Perform the short-range RESPSA steps (velocity Verlet increments of
75 c positions and velocities using short-range forces)
76 c write (iout,*) "*********************** RESPA split"
77 do itsplit=1,ntime_split
80 else if (lang.eq.2 .or. lang.eq.3) then
81 call stochastic_force(stochforcvec)
83 c First step of the velocity Verlet algorithm
86 else if (lang.eq.3) then
87 call sd_verlet1_ciccotti
88 else if (lang.eq.1) then
93 c Build the chain from the newly calculated coordinates
95 if (rattle) call rattle1
96 if (large.and. mod(itime,ntwe).eq.0) then
97 write (iout,*) "Cartesian and internal coordinates: step 1"
100 write (iout,*) "Accelerations"
102 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
103 & (d_a(j,i+nres),j=1,3)
105 write (iout,*) "Velocities, step 1"
107 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
108 & (d_t(j,i+nres),j=1,3)
112 c Calculate energy and forces
114 call etotal_short(energia_short)
116 c Get the new accelerations
118 t_enegrad=t_enegrad+tcpu()-tt0
119 c Second step of the velocity Verlet algorithm
122 else if (lang.eq.3) then
123 call sd_verlet2_ciccotti
124 else if (lang.eq.1) then
129 if (rattle) call rattle2
130 c Backup the coordinates, velocities, and accelerations
134 d_t_old(j,i)=d_t(j,i)
135 d_a_old(j,i)=d_a(j,i)
139 c Restore the time step
141 c Compute long-range forces
143 call etotal_long(energia_long)
145 c Compute accelerations from long-range forces
147 if (large.and. mod(itime,ntwe).eq.0) then
148 write (iout,*) "Cartesian and internal coordinates: step 2"
151 write (iout,*) "Accelerations"
153 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
154 & (d_a(j,i+nres),j=1,3)
156 write (iout,*) "Velocities, step 2"
158 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
159 & (d_t(j,i+nres),j=1,3)
162 c Compute the final RESPA step (increment velocities)
163 c write (iout,*) "*********************** RESPA fin"
165 c Compute the complete potential energy
166 potE=energia_short(0)+energia_long(0)
168 c Calculate the kinetic and the total energy and the kinetic temperature
171 c Couple the system to Berendsen bath if needed
172 if (tbf .and. lang.eq.0) then
175 kinetic_T=2.0d0/(dimen*Rb)*EK
176 c Backup the coordinates, velocities, and accelerations
177 if (mod(itime,ntwe).eq.0 .and. large) then
178 write (iout,*) "Velocities, end"
180 write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
181 & (d_t(j,i+nres),j=1,3)