1c
2c
3c     ###################################################
4c     ##  COPYRIGHT (C)  2011  by  Jay William Ponder  ##
5c     ##              All Rights Reserved              ##
6c     ###################################################
7c
8c     #############################################################
9c     ##                                                         ##
10c     ##  subroutine respa  --  r-RESPA molecular dynamics step  ##
11c     ##                                                         ##
12c     #############################################################
13c
14c
15c     "respa" performs a single multiple time step molecular dynamics
16c     step using the reversible reference system propagation algorithm
17c     (r-RESPA) via a Verlet core with the potential split into fast-
18c     and slow-evolving portions
19c
20c     literature references:
21c
22c     D. D. Humphreys, R. A. Friesner and B. J. Berne, "A Multiple-
23c     Time-Step Molecular Dynamics Algorithm for Macromolecules",
24c     Journal of Physical Chemistry, 98, 6885-6892 (1994)
25c
26c     X. Qian and T. Schlick, "Efficient Multiple-Time-Step Integrators
27c     with Distance-Based Force Splitting for Particle-Mesh-Ewald
28c     Molecular Dynamics Simulations", Journal of Chemical Physics,
29c     115, 4019-4029 (2001)
30c
31c
32      subroutine respa (istep,dt)
33      use atomid
34      use atoms
35      use freeze
36      use ielscf
37      use mdstuf
38      use moldyn
39      use polar
40      use units
41      use usage
42      use virial
43      implicit none
44      integer i,j,k,m
45      integer istep
46      integer nalt
47      real*8 dt,dt_2
48      real*8 dta,dta_2
49      real*8 epot,etot
50      real*8 eksum,eps
51      real*8 temp,pres
52      real*8 ealt,dalt
53      real*8 term
54      real*8 ekin(3,3)
55      real*8 stress(3,3)
56      real*8 viralt(3,3)
57      real*8, allocatable :: xold(:)
58      real*8, allocatable :: yold(:)
59      real*8, allocatable :: zold(:)
60      real*8, allocatable :: derivs(:,:)
61c
62c
63c     set some time values for the dynamics integration
64c
65      eps =  0.00000001d0
66      dalt = arespa
67      nalt = int(dt/(dalt+eps)) + 1
68      dalt = dble(nalt)
69      dt_2 = 0.5d0 * dt
70      dta = dt / dalt
71      dta_2 = 0.5d0 * dta
72c
73c     store the current atom positions, then find half-step
74c     velocities via velocity Verlet recursion
75c
76      do i = 1, nuse
77         m = iuse(i)
78         do j = 1, 3
79            v(j,m) = v(j,m) + a(j,m)*dt_2
80         end do
81      end do
82c
83c     initialize virial from fast-evolving potential energy terms
84c
85      do i = 1, 3
86         do j = 1, 3
87            viralt(j,i) = 0.0d0
88         end do
89      end do
90c
91c     perform dynamic allocation of some local arrays
92c
93      allocate (xold(n))
94      allocate (yold(n))
95      allocate (zold(n))
96      allocate (derivs(3,n))
97c
98c     find fast-evolving velocities and positions via Verlet recursion
99c
100      do k = 1, nalt
101         do i = 1, nuse
102            m = iuse(i)
103            do j = 1, 3
104               v(j,m) = v(j,m) + aalt(j,m)*dta_2
105            end do
106            xold(m) = x(m)
107            yold(m) = y(m)
108            zold(m) = z(m)
109            x(m) = x(m) + v(1,m)*dta
110            y(m) = y(m) + v(2,m)*dta
111            z(m) = z(m) + v(3,m)*dta
112         end do
113         if (use_rattle)  call rattle (dta,xold,yold,zold)
114c
115c     get the fast-evolving potential energy and atomic forces
116c
117         call gradfast (ealt,derivs)
118c
119c     use Newton's second law to get fast-evolving accelerations;
120c     update fast-evolving velocities using the Verlet recursion
121c
122         do i = 1, nuse
123            m = iuse(i)
124            do j = 1, 3
125               aalt(j,m) = -ekcal * derivs(j,m) / mass(m)
126               v(j,m) = v(j,m) + aalt(j,m)*dta_2
127            end do
128         end do
129         if (use_rattle)  call rattle2 (dta)
130c
131c     find average virial from fast-evolving potential terms
132c
133         do i = 1, 3
134            do j = 1, 3
135               viralt(j,i) = viralt(j,i) + vir(j,i)/dalt
136            end do
137         end do
138      end do
139c
140c     apply Verlet half-step updates for any auxiliary dipoles
141c
142      if (use_ielscf) then
143         do i = 1, nuse
144            m = iuse(i)
145            do j = 1, 3
146               vaux(j,m) = vaux(j,m) + aaux(j,m)*dt_2
147               vpaux(j,m) = vpaux(j,m) + apaux(j,m)*dt_2
148               uaux(j,m) = uaux(j,m) + vaux(j,m)*dt
149               upaux(j,m) = upaux(j,m) + vpaux(j,m)*dt
150            end do
151         end do
152      end if
153c
154c     get the slow-evolving potential energy and atomic forces
155c
156      call gradslow (epot,derivs)
157      epot = epot + ealt
158c
159c     make half-step temperature and pressure corrections
160c
161      call temper2 (dt,temp)
162      call pressure2 (epot,temp)
163c
164c     use Newton's second law to get the slow accelerations;
165c     find full-step velocities using velocity Verlet recursion
166c
167      do i = 1, nuse
168         m = iuse(i)
169         do j = 1, 3
170            a(j,m) = -ekcal * derivs(j,m) / mass(m)
171            v(j,m) = v(j,m) + a(j,m)*dt_2
172         end do
173      end do
174c
175c     apply Verlet full-step updates for any auxiliary dipoles
176c
177      if (use_ielscf) then
178         term = 2.0d0 / (dt*dt)
179         do i = 1, nuse
180            m = iuse(i)
181            do j = 1, 3
182               aaux(j,m) = term * (uind(j,m)-uaux(j,m))
183               apaux(j,m) = term * (uinp(j,m)-upaux(j,m))
184               vaux(j,m) = vaux(j,m) + aaux(j,m)*dt_2
185               vpaux(j,m) = vpaux(j,m) + apaux(j,m)*dt_2
186            end do
187         end do
188      end if
189c
190c     perform deallocation of some local arrays
191c
192      deallocate (xold)
193      deallocate (yold)
194      deallocate (zold)
195      deallocate (derivs)
196c
197c     find the constraint-corrected full-step velocities
198c
199      if (use_rattle)  call rattle2 (dt)
200c
201c     increment total virial from sum of fast and slow parts
202c
203      do i = 1, 3
204         do j = 1, 3
205            vir(j,i) = vir(j,i) + viralt(j,i)
206         end do
207      end do
208c
209c     make full-step temperature and pressure corrections
210c
211      call temper (dt,eksum,ekin,temp)
212      call pressure (dt,epot,ekin,temp,pres,stress)
213c
214c     total energy is sum of kinetic and potential energies
215c
216      etot = eksum + epot
217c
218c     compute statistics and save trajectory for this step
219c
220      call mdstat (istep,dt,etot,epot,eksum,temp,pres)
221      call mdsave (istep,dt,epot,eksum)
222      call mdrest (istep)
223      return
224      end
225c
226c
227c     ##################################################################
228c     ##                                                              ##
229c     ##  subroutine gradfast  --  fast energy & gradient components  ##
230c     ##                                                              ##
231c     ##################################################################
232c
233c
234c     "gradfast" calculates the potential energy and first derivatives
235c     for the fast-evolving local valence potential energy terms
236c
237c
238      subroutine gradfast (energy,derivs)
239      use limits
240      use potent
241      implicit none
242      real*8 energy
243      real*8 derivs(3,*)
244      logical save_vdw,save_repuls
245      logical save_disp,save_charge
246      logical save_chgdpl,save_dipole
247      logical save_mpole,save_polar
248      logical save_chgtrn,save_rxnfld
249      logical save_solv,save_list
250c
251c
252c     save the original state of slow-evolving potentials
253c
254      save_vdw = use_vdw
255      save_repuls = use_repuls
256      save_disp = use_disp
257      save_charge = use_charge
258      save_chgdpl = use_chgdpl
259      save_dipole = use_dipole
260      save_mpole = use_mpole
261      save_polar = use_polar
262      save_chgtrn = use_chgtrn
263      save_rxnfld = use_rxnfld
264      save_solv = use_solv
265      save_list = use_list
266c
267c     turn off slow-evolving nonbonded potential energy terms
268c
269      use_vdw = .false.
270      use_repuls = .false.
271      use_disp = .false.
272      use_charge = .false.
273      use_chgdpl = .false.
274      use_dipole = .false.
275      use_mpole = .false.
276      use_polar = .false.
277      use_chgtrn = .false.
278      use_rxnfld = .false.
279      use_solv = .false.
280      use_list = .false.
281c
282c     get energy and gradient for fast-evolving potential terms
283c
284      call gradient (energy,derivs)
285c
286c     restore the original state of slow-evolving potentials
287c
288      use_vdw = save_vdw
289      use_repuls = save_repuls
290      use_disp = save_disp
291      use_charge = save_charge
292      use_chgdpl = save_chgdpl
293      use_dipole = save_dipole
294      use_mpole = save_mpole
295      use_polar = save_polar
296      use_chgtrn = save_chgtrn
297      use_rxnfld = save_rxnfld
298      use_solv = save_solv
299      use_list = save_list
300      return
301      end
302c
303c
304c     ##################################################################
305c     ##                                                              ##
306c     ##  subroutine gradslow  --  slow energy & gradient components  ##
307c     ##                                                              ##
308c     ##################################################################
309c
310c
311c     "gradslow" calculates the potential energy and first derivatives
312c     for the slow-evolving nonbonded potential energy terms
313c
314c
315      subroutine gradslow (energy,derivs)
316      use potent
317      implicit none
318      real*8 energy
319      real*8 derivs(3,*)
320      logical save_bond,save_angle
321      logical save_strbnd,save_urey
322      logical save_angang,save_opbend
323      logical save_opdist,save_improp
324      logical save_imptor,save_tors
325      logical save_pitors,save_strtor
326      logical save_angtor,save_tortor
327      logical save_geom,save_metal
328      logical save_extra
329c
330c
331c     save the original state of fast-evolving potentials
332c
333      save_bond = use_bond
334      save_angle = use_angle
335      save_strbnd = use_strbnd
336      save_urey = use_urey
337      save_angang = use_angang
338      save_opbend = use_opbend
339      save_opdist = use_opdist
340      save_improp = use_improp
341      save_imptor = use_imptor
342      save_tors = use_tors
343      save_pitors = use_pitors
344      save_strtor = use_strtor
345      save_angtor = use_angtor
346      save_tortor = use_tortor
347      save_geom = use_geom
348      save_metal = use_metal
349      save_extra = use_extra
350c
351c     turn off fast-evolving valence potential energy terms
352c
353      use_bond = .false.
354      use_angle = .false.
355      use_strbnd = .false.
356      use_urey = .false.
357      use_angang = .false.
358      use_opbend = .false.
359      use_opdist = .false.
360      use_improp = .false.
361      use_imptor = .false.
362      use_tors = .false.
363      use_pitors = .false.
364      use_strtor = .false.
365      use_angtor = .false.
366      use_tortor = .false.
367      use_geom = .false.
368      use_metal = .false.
369      use_extra = .false.
370c
371c     get energy and gradient for slow-evolving potential terms
372c
373      call gradient (energy,derivs)
374c
375c     restore the original state of fast-evolving potentials
376c
377      use_bond = save_bond
378      use_angle = save_angle
379      use_strbnd = save_strbnd
380      use_urey = save_urey
381      use_angang = save_angang
382      use_opbend = save_opbend
383      use_opdist = save_opdist
384      use_improp = save_improp
385      use_imptor = save_imptor
386      use_tors = save_tors
387      use_pitors = save_pitors
388      use_strtor = save_strtor
389      use_angtor = save_angtor
390      use_tortor = save_tortor
391      use_geom = save_geom
392      use_metal = save_metal
393      use_extra = save_extra
394      return
395      end
396