Practical Construction of Modified Hamiltonians.
One of the most fruitful ways to analyze the effects of discretization error in the numerical solution of a system of differential equations is to examine the "modified equations," which are equations that are exactly satisfied by the (approximate) discrete solution.  These do not actually exist in general but rather are defined by an asymptotic expansion in powers of the discretization parameter. Nonetheless, if the expansion is suitably truncated, the resulting modified equations have a solution which is remarkably close to the discrete solution. In the case of a Hamiltonian system of ordinary differential equations, the modified equations are also Hamiltonian if and only if the integrator is symplectic. Evidence for the existence of a Hamiltonian for a particular calculation is obtained by calculating modified Hamiltonians and monitoring how well they are conserved. Also, energy drifts caused by numerical instability are better revealed by evaluating modified Hamiltonians. Doing this calculation would normally be complicated and highly dependent on the details of the method, even if differences are used to approximate derivatives. A relatively simple procedure is presented here, nearly independent of the internal structure of the integrator, for obtaining highly accurate estimates for modified Hamiltonians. As a bonus of the method of construction, the modified Hamiltonians are exactly conserved by a numerical solution in the case of a quadratic Hamiltonian.
Introduction
One of the most fruitful ways to analyze the effects of discretization error in the numerical solution
of differential equations is to examine the "modified equations," which are the equations that are
exactly satisfied by the (approximate) discrete solution. These do not actually exist (in general),
but rather are defined by an asymptotic expansion in powers of the discretization parameter.
Nonetheless, if the expansion is suitably truncated, the resulting modified equations have a solution
which is remarkably close to the discrete solution [9]. In the case of a Hamiltonian system of
The work of R. D. Skeel was supported in part by NSF Grants DMS-9971830, DBI-9974555 and NIH Grant
P41RR05969 and completed while visiting the Mathematics Department, University of California, San Diego.
y
z
ordinary differential equations, the modified equations are also Hamiltonian if and only if the
integrator is symplectic. The existence of a modified, or "shadow" [4], Hamiltonian is an indicator
of the validity of statistical estimates calculated from long time integration of chaotic Hamiltonian
systems [18]. In addition, the modified Hamiltonian is a more sensitive indicator than is the
original Hamiltonian of drift in the energy (caused by instability). Evidence for the existence of a
Hamiltonian for a particular calculation can be obtained by calculating modified Hamiltonians and
monitoring how well they are conserved. Doing this calculation would normally be complicated and
highly dependent on the details of the method, even if differences are used to approximate higher
derivatives. Presented here is a relatively simple procedure, nearly independent of the internal
structure of the integrator, for obtaining highly accurate estimates for modified Hamiltonians.
Consider a step by step numerical integrator x evolves an approximate
solution x n - x(nh) for a system of ordinary differential equations -
f(x). For such discrete
solutions there exists modified equations -
defined by an asymptotic expansion such
that formally the numerical solution x (nh). The modified right-hand-side function f h is
defined uniquely by postulating an asymptotic expansion f in powers of h,
substituting this into the equations for the numerical solution, expanding in powers of h, and
equating coefficients [26, 6, 22]. The asymptotic expansion does not generally converge except for
(reasonable integrators applied to) linear differential equations.
A Hamiltonian system is of the form
I
for some Hamiltonian H(x), . The modified equation for an integrator \Phi h applied to
this system is Hamiltonian, i.e., f modified Hamiltonian H h (x), if and only if
the integrator is symplectic [23, 20]. The integrator is symplectic if \Phi h;x (x)J T \Phi h;x (x) j J . There
is theoretical [2, 8, 18] and empirical evidence that
very small error
for a very long time where x h is the solution for a suitably truncated Hamiltonian H h . In what
follows we assume that H h is such a Hamiltonian and we neglect the very small error.
If we plot total energy as a function of time for a numerical integrator such as leapfrog/St-ormer/
Verlet applied to a molecular dynamics simulation, we get a graph like Fig. 3. What we observe
are large fluctuations in the original Hamiltonian, as the trajectory moves on a hypersurface of
constant modified Hamiltonian. A small drift or jump in the energy would be obscured by the
fluctuations. A plot of a modified Hamiltonian might be more revealing. As an example, the plots
of modified Hamiltonians in Fig. 4 show a clear rise in energy already in a 400-step simulation.
This indicates that plots of suitable modified Hamiltonians can make it easier to test integration
algorithms for instability and programming bugs. Details of this and other numerical tests are
given in section 2. Before continuing, it is worth emphasizing that the concern of this paper is
stability monitoring-not the monitoring and enhancement of accuracy, as in [4] and [15].
The goal is to construct an approximate modified Hamiltonian
that can be conveniently assembled from quantities, such as forces and energies, already available
from the numerical integration. We consider the special separable Hamiltonian H(q;
which the system is of the form
A "brute force" approach would be to determine an asymptotic expansion for H h and of the
quantities available for making an approximation and then to form a suitable linear combination
of the latter. By such a matching of asymptotic expansions one could derive the following modified
Hamiltonians for the leapfrog method:
Here a superscript n denotes evaluation at q n , the centered difference operator is defined by ffi w
the averaging operator is defined by -w
are defined in terms of q n , p n by the leapfrog method.
An easier and more elegant construction is presented in Secs. 3-5. The technique is developed
only for splitting methods. It is likely that a similar construction is also possible for symplectic
implicit Runge-Kutta methods. The idea is to add a new position and conjugate momentum variable
to get an extended Hamiltonian -
H h (y) which is homogeneous of order 2. For such a Hamiltonian
Jy h (t). Thus the problem is reduced to that of forming an approximation for
using the numerical solution of an extended Hamiltonian system. It is plausible that such a
construction might be useful theoretically due to the existence of robust approximation techniques.
Eq. (1) for H [2] contains an h 2 term which is not needed for achieving 2nd order accuracy. It
is present because the truncations H [2k] are designed to exactly conserve energy for the numerical
solution when H is quadratic. (See Secs. 4.1 and 5.1.) This is a very useful property because typical
applications, including molecular dynamics, are dominated by harmonic motion. The existence of
a modified Hamiltonian that is exactly conserved for a quadratic Hamiltonian is noted in [17,
Eq. (4.7b)], and the search for similar methods having this property was central to the results of
this paper. For a quadratic Hamiltonian the modified Hamiltonian H h actually exists (if h is not
too large), but H [2k] 6= H h . simple derivation of H h for the one-dimensional case is given in [22].)
Also, it should be noted that the Hamiltonians H [2k] will not detect numerical instability in the
case of quadratic Hamiltonians H .
The modified Hamiltonians H [2k] (x), defined by Eqs. (15), (8), (14), (7) are computed
and plotted as functions of time for numerical solutions generated by the leapfrog method, given
energy
time (fs)
8th order
6th order
4th order
energy
2nd order

Figure

1: Energy and various truncations of modified Hamiltonians for decalanine.
by (9). The unmodified Hamiltonians are those of classical molecular dynamics. The testing was
done with a molecular dynamics program written by the second author, which is compatible with
NAMD [11, 16] but limited in features to facilitate algorithm testing.
The first couple of experiments demonstrate the quality of the modified Hamiltonians. The test
problem is a 66-atom peptide, decalanine, in a vacuum [1]. The force field parameters are those of
CHARMM 22 for proteins [13, 14] without cutoffs for nonbonded forces.

Figure

1 shows a plot of the Hamiltonian and 2nd, 4th, 6th, and 8th order modified Hamiltonian
approximations vs. time for 100000 fs (femtoseconds) for a step size with the energy
sampled every 8th step. The level graph at the top is the 8th order truncation, the one just barely
beneath it is 6th order, and the one under that is 4th order. The greatly fluctuating graph is the
energy itself and the undulating one well below it is the 2nd order truncation. Note how well the
asymptotic theory holds for the higher order truncations-one could not obtain such flat plots by
simply smoothing the original Hamiltonian.

Figure

2 expands the vertical scale to show fluctuations in the 8th, 6th, and 4th order truncations
of modified Hamiltonians.
An explanation is in order concerning the initial drop in energy. Because a symplectic method
preserves volume in phase space and because there is less phase space volume at lower energies,
it can be inferred that the first part of the trajectory is simply the second half of a very unusual
fluctuation. In other words the initial conditions are atypical, i.e., not properly equilibrated (with
respect to the original Hamiltonian). This is particularly well revealed by the plot of the 2nd order
truncation.
The remaining experiments demonstrate the ability of modified Hamiltonians to detect instabil-
ity. The test problem is a set of 125 water molecules harmonically restrained to a 10 - A-radius sphere.
The water is based on the TIP3P model [25] without cutoffs and with flexibility incorporated by
adding bond stretching and angle bending harmonic terms (cf. Ref. [12]).
energy
time (fs)
8th order
6th order
4th order

Figure

2: Closer look at higher order truncations for decalanine.
energy
time (fs)
energy

Figure

3: Energy for flexible water with step size 2.5 fs.

Figure

3 shows a plot of the energy vs. time for 1 000 fs for a step size with the energy
sampled every step. Note that the large fluctuations make it is difficult to determine whether or
not there is energy drift.

Figure

4 shows a plot of the 6th and 8th order modified Hamiltonians for the same step size
energy
time (fs)
8th order
6th order

Figure

4: 6th and 8th order truncations with step size 2.5 fs.
fs. An upward energy drift is now obvious. The 2nd and 4th order approximations are
not shown because neither of them were as flat. Normal mode analysis for this system [10] shows
that the 250 fastest frequencies have periods in the range 9.8-10.2 fs and use of the formula in [22,
p. 131] shows that a 2.5 fs step size is 30% of the effective period for discrete leapfrog dynamics. It
is remarkable that the 8th order approximation is the flattest, even for such a large step size.

Figure

5, shows a plot of the 6th and 8th order modified Hamiltonians for step size
There is no apparent upward drift of the energy. Theoretically instability due to 4:1 resonance [21]
should occur for the leapfrog method at h \Delta angular frequency =
2, which is in the range 2.2-2.3 fs
for flexible water.
3 Augmenting the Integrator
We assume that one step of size h for the given method applied to a system with Hamiltonian
H is the composition of exact h-flows for Hamiltonian systems with Hamiltonians H 1 , H 2 , . ,
HL . Each H l (x) is assumed to be sufficiently smooth on some domain containing the infinite time
trajectory. For example,
1. the leapfrog method for separable Hamiltonian systems H(q;
2. the Rowlands method [19] for special separable Hamiltonian systems uses H 1
3. double time-stepping [7, 24] uses
4 U fast (q), H 2
fast (q),
energy
time (fs)
8th order
6th order

Figure

5: 6th and 8th order truncations with step size 2.15 fs.
4. Molly [5] does the same as double time-stepping except for the substitution of U slow (A(q))
for U slow (q) where A(q) is a local temporal averaging of q over vibrational motion.
We defined the homogeneous extension of a Hamiltonian by
If H is quadratic, then -
H is homogeneous of order 2:
. The extended Hamiltonian yields the augmented system
With initial condition and the system simplifies to
For
p+U(q), the extended Hamiltonian is -
and the simplified augmented system is
Remark. The association of ff with q rather than fi is of practical importance in that we want
to get values of -
p is calculated.
The following proposition shows that the value of the extended Hamiltonian can be calculated
knowing just the solution:
H(y) be the homogeneous extension of a given Hamiltonian H(x), and let y(t)
be a solution of the extended Hamiltonian system with ff initially 1. Then
J is the matrix
I
of augmented dimension.
Proof. Differentiating Eq. (4) with respect to oe gives -
Because -
H is a homogeneous extension of H , the solution of -
H "includes" that of H and we have
Of course, the goal is not to calculate the original Hamiltonian, for which we know a formula
but not the solution; rather, it is to calculate a modified Hamiltonian, for which we know the
solution (at grid points) but not a formula. Therefore, we must augment the integrator so that
its solution at grid points is that of the homogeneous extension of the modified Hamiltonian.
For an integrator that is a composition of Hamiltonian flows, this is accomplished by using the
homogeneous extension of each of the constituent Hamiltonians. More specifically, we define the
augmented method y
H to be the composition of exact flows for systems with
Hamiltonians -
HL where
H l (q; ff;
Lemma 1 (Commutativity) The method \Psi h defined above has a modified Hamiltonian
where H h (q; p) is the modified Hamiltonian of the original method \Phi h , i.e., the following diagram
commutes:
homogeneous extension
discretization # discretization
homogeneous extension
Proof. The modified Hamiltonian -
H h for method \Psi h can be expressed as an asymptotic expansion
using the Baker-Campbell-Hausdorf formula [20]. This formula combines Hamiltonians using the
operations of scalar multiplication, addition, and the Poisson bracket fH;
x JN x . It is thus
sufficient to show that each of these three commute with the operation of forming the homogeneous
extension. We show this only for the last of these. The homogeneous extension of the Poisson
bracket is
This is exactly the same as the (extended) Poisson bracket of ff
Remark. The aim is to discretize the extended Hamiltonian so that this commutativity property
holds. Extension of this technique to implicit Runge-Kutta methods would require an augmentation
of the method so that commutativity holds.
The following corollary allows the value of the Hamiltonian to be approximated from known
values of y h (t) at grid points:
Proposition 2 Let x h (t) and y h (t) be the solutions for modified Hamiltonians H h and -
tively. Then
Jy h (t):
Proof. Similar to that of Prop. 1. 2
4 Using Full Step Values
This section presents the construction of H [2k] for even values of k.
Let y h (t) be the solution of the modified extended Hamiltonian system with initial condition
y. It has values y h
(t) be the degree k polynomial
interpolant of these values. (For large k it may be preferable, instead, to use trigonometric
interpolation suitably modified [3].)
From Prop. 2, 1-
jh
Z jh=2
\Gammaj h=22 -
The interpolant - k where the error
(t) and
with the brackets denoting a 1)th divided difference. Noting
that -
jh
Z jh=2
\Gammaj h=22 -
Je(t)dt \Gammajh
Z jh=2
\Gammaj h=22 y h (t) T -
Z jh=2
\Gammaj h=22 -
Je(t)dt +O(h 2k+2 )
Z jh=2
\Gammaj h=2
where the second equation is obtained by integrating by parts and where fl(t) def
Jy [k+1]
(t).
This can be expressed as an expansion
By forming a suitable linear combination of the values H k;j , it is expected that
one can get -
H h with the first k=2 \Gamma 1 leading error terms eliminated:
linear combination of the H
Note. The value -
contains a leading term that is only O(h k ), so it is not useful for
eliminating error terms.
The case is the 4th order accurate formula
For the case
\Gammah
and
Z 2h
\Gamma2h
and hence,
Below are given formulas for H [8] and for H [4] in terms of values of y h (t) at grid points. Let a j
be the jth centered difference of y h (t) at
where the centered difference operator is defined by ffi and the
averaging operator is defined by
The 4th degree interpolant in divided difference form is
Hence,
6 a 3 s(s
and
6 a 4 s(s
and we have2
powers of s:
Averaging over
and averaging over \Gamma2 - s - 2 yields
90 A 14 \Gamma 107
Therefore,
For a 2nd degree interpolant it follows from Eq. (5) that
An implementation of these formulas might calculate H [2k] consecutive values of n in
terms of quantities A n
defined in terms of centered differences of y n which can
be obtained from the x n . (Only 1st and higher differences of fi n are needed.)
Example 1. To make this concrete, we calculate H [4] (x) for the leapfrog method, as given by
Eq. (2). The leapfrog method advances one step by
We have
Suppressing the n in the superscript,
y
whence
not needed7 7
From
so
Therefore,
4.1 The case of a quadratic Hamiltonian
The following result implies that, in the case where H(x) is quadratic, the numerical solution
exactly conserves an approximate modified Hamiltonian which is a linear functional of 1-
where -(t) is a linear combination of numerical solution values.
Proposition 3 Assume that \Phi h is the composition of flows for systems with quadratic Hamiltonians
and that \Psi h is constructed as in Proposition 1. Then the quantity
a i;j \Psi i
where the sum is taken over a finite set of pairs of integers, is exactly conserved by method \Psi h .
Proof. The mapping \Psi h is the composition of
flows for systems with homogeneous quadratic Hamiltonians. Then
a i;j (S i Sy) T -
a
a i;j \Psi i
(y):5 Using Intermediate Values
This section presents the construction of H [2k] for odd values of k.
For most numerical integrators one can define "sensible" mid-step values,
and these can be used instead of full step values to get an estimate accurate up to O(h 2k ). We
assume that \Psi
\Gammah=2
\Psi h=2 is a composition of exact flows of homogeneously
extended Hamiltonians.
Remark. It is not necessary that the mid-step values be approximations to y(t) at midpoints
nor that \Psi h be time symmetric (\Psi \Gamma1
we need is that \Psi where each of
\Psi 1;h , \Psi 2;h is a composition of exact flows of homogeneously extended Hamiltonians.
For example, the leapfrog method separates into half steps -
\Gammah=2
\Psi h=2 with -
\Psi h=2 as follows:
Remark. For the leapfrog method the estimate over an interval from (n\Gamma 1
2 k)h to (n+ 1
k is odd actually uses only values of energy and forces from the shorter interval from (n \Gamma 1k
to (n
The mid-step values are values at midpoints of some function z h (t) which can be used to
construct the Hamiltonian:
Proposition 4 Let z h
z h (t) T -
Jz h (t):
Proof. For any real s we define \Psi s
sh-flow for -
. Because - h is symplectic, z h Hamiltonian system
with Hamiltonian -
. Also, - h is the composition of flows of
Hamiltonians that are 2nd order homogeneous, and hence -
is homogeneous of order 2.
z h (t) T -
(t) be the degree k polynomial interpolant of the values z h
(y),
As before, let
jh
Z jh=2
\Gammaj h=22 -
The interpolant - k where the error
h (t) and
z [k+1]
Similar to before, we get
Z jh=2
\Gammaj h=2
z h (t) T -
Jz [k+1]
h (t). This can be expressed as an expansion
Again, it is expected that a suitable linear combination of the different values of H k;j
yields -
H h with the first (k leading error terms eliminated:
linear combination of the H
Note. It seems not to be possible to combine values obtained from full steps with those from
half steps to further increase the order of accuracy, because the error expansions for the two kinds
of averages do not have terms in common that can cancel.
For we have the 2nd order formula
For
Z h=2
\Gammah=2
and
Z 3h=2
\Gamma3h=2
e
Let b j be the jth centered difference of z h (t) at using mid-step values:
The 3rd degree interpolant is
Hence,
and
and we
powers of s:
Averaging
and averaging over \Gamma 3
yields
Therefore,
For a 1st degree interpolant it follows from Eq. (13) that
Example 2. We calculate H [2] (x) for the leapfrog method, as given by Eq. (1). We have
J(-y n
Suppressing the n in the superscript,
y
e
Therefore,
Example 3. We calculate H [6] (x) for the leapfrog method, as given by Eq. (3). From Eqs. (16)
we have
We have h-ffiF from Eq. (11b), and ffi 2 (q T F ) is given by Eq. (12). Then
and, therefore,
5.1 The case of a quadratic Hamiltonian
Proposition 5 Assume that \Phi h is the composition of flows for systems with quadratic Hamiltoni-
ans, that \Psi h is constructed as in Proposition 1, and that -
\Psi h=2 is as assumed at the beginning of
this section. Then the quantity
a i;j
where the sum is taken over a finite set of pairs of integers, is exactly conserved by method \Psi h .
Proof. Because -
are the compositions of flows for systems with homogeneous quadratic
Hamiltonians, the mappings -
a
a
a i;j
Acknowledgment

The authors are grateful for the assistance of Justin Wozniak, who did preliminary tests of the
second order truncation for the H'enon-Heiles Hamiltonian and for decalanine.



--R

http://www.
On the Hamiltonian interpolation of near to the identity symplectic mappings with application to symplectic integration algorithms.
Residual acceleration as a measure of the accuracy of molecular dynamics simulations.
Shadow mass and the relationship between velocity and momentum in symplectic numerical integration.

On the scope of the method of modified equations.
Generalized Verlet algorithm for efficient molecular dynamics simulations with long-range interactions
The life-span of backward error analysis for numerical integrators
Asymptotic expansions and backward analysis for numerical inte- grators
Longer time steps for molecular dynamics.
NAMD2: Greater scalability for parallel molecular dynamics.
Molecular Modelling: Principles and Applications.


Common molecular dynamics algorithms revisited: Accuracy and optimal time steps of St-ormer-leapfrog

An analysis of the accuracy of Langevin and molecular dynamics algorithm.
Backward error analysis for numerical integrators.
A numerical algorithm for Hamiltonian systems.
Numerical Hamiltonian Problems.
Nonlinear resonance artifacts in molecular dynamics simulations.
Integration schemes for molecular dynamics and related applications.
Some Geometric and Numerical Methods for Perturbed Integrable Systems.
Reversible multiple time scale molecular dynamics.

The modified equation approach to the stability and accuracy analysis of finite difference methods.
--TR

--CTR
Robert D. Engle , Robert D. Skeel , Matthew Drees, Monitoring energy drift with shadow Hamiltonians, Journal of Computational Physics, v.206 n.2, p.432-452, 1 July 2005
Jess A. Izaguirre , Scott S. Hampton, Shadow hybrid Monte Carlo: an efficient propagator in phase space of macromolecules, Journal of Computational Physics, v.200 n.2, p.581-604, November 2004
integrator;symplectic;modified equation;backward error;hamiltonian;numerical