Effective conservation of energy and momentum algorithm using switching
  potentials suitable for molecular dynamics simulation of thermodynamical
  systems
Effective conservation of energy and momentum
algorithm using switching potentials suitable for
molecular dynamics simulation of
thermodynamical systems
Christopher G. Jesudason ∗
Laboratory of Physics and Helsinki Institute of Physics,
P.O.Box 1100, FIN-02015 HUT, Finland.
Email: jesu@um.edu.my, chrysostomg@gmail.com
4 April, 2007
Abstract
During a crossover via a switching mechanism from one 2-body poten-
tial to another as might be applied in modeling (chemical) reactions in
the vicinity of bond formation, energy violations would occur due to finite
step size which determines the trajectory of the particles relative to the
potential interactions of the unbonded state by numerical (e.g. Verlet)
integration. This problem is overcome by an algorithm which preserves
the coordinates of the system for each move, but corrects for energy dis-
crepancies by ensuring both energy and momentum conservation in the
dynamics. The algorithm is tested for a hysteresis loop reaction model
with an without the implementation of the algorithm. The tests involve
checking the rate of energy flow out of the MD simulation box; in the
equilibrium state, no net rate of flows within experimental error should
be observed. The temperature and pressure of the box should also be
invariant within the range of fluctuation of these quantities. It is demon-
strated that the algorithm satisfies these criteria
AMS (MSC2000) Subject Classification. 00A71-2, 70H05, 80A20
1 PRELIMINARIES
The dimeric particle reaction simulated may be written
A2 (1.1)
∗on leave from Chemistry Department, University of Malaya, 50603 Kuala Lumpur,
Malaysia.
http://arxiv.org/abs/0704.0556v1
where k1 is the forward rate constant and k−1 is the backward rate constant.
The reaction simulation was conducted at extremely high temperatures which
are off-scale and not used in ordinary simulations of LJ (Lennard-Jones) fluids
where normally [1] the reduced temperatures T ∗ ranges ∼ 0.3 − 1.2, whereas
here, T ∗ ∼ 8.0 − 16.0, well above the supercritical regime of the LJ fluid At
these temperatures, the normal choices for time step increments do not obtain
without also taking into account energy-momentum conservation algorithms in
regions where there are abrupt changes of gradient [1, 2, 3]. The global literature
does not seem to cover such extreme conditions of simulation with discrete time
steps using the Verlet velocity algorithm. The units used here are reduced LJ
ones [1]. The simulation was at density ρ = 0.70 with 4096 atomic particls which
could react. The potentials used are as given in Fig. (1) where rb = 1.20 for the
vicinity where the bond of the dimer is broken and where 2 free particles emerge,
and rf = 0.85 is the point along the hysteresis potential curve where the dimer
is defined to exist for two previously free particles which collide. The reaction
proceeds as follows; all particles interact with the splined LJ pair potential uLJ
except for the dimeric pair (i, j) formed from particles i and j which interact
with a harmonic-like intermolecular potential modified by a switch u(r) given
u(r) = uvib(r)s(r) + uLJ [1− s(r)] (1.2)
where uvib(r) is the vibrational potential given by eq.(1.3) below
uvib(r) = u0 +
k(r − r0)
2 (1.3)
The switching function s(r) is defined as
s(r) =
)n (1.4)
where
s(r) → 1 if r < rsw
s(r) → 0 for r > rsw
The switching function becomes effective when the distance between the atoms
approach the value rsw (see Fig. (1) ). Some of the other parameters used in
the equations that follow include:
u0 = −10, r0 = 1.0, k ∼ 2446 (exact value is determined by the other input pa-
rameters), n = 100, rf = 0.85, rb = 1.20, and rsw = 1.11. Particles i and j above
also interact with all other particles not bonded to it via uLJ . Full simulation
details are given elsewhere [2]; suffice to say the activation energy at rf is ex-
tremely high at approximately 17.5. At rf , the molecular potential is turned on
where at this point there is actually a crossing of the potential curves although
the gradients of the molecular and free uLJ potentials are ”‘very close”’. On
the other hand, at rb , the switch forces the two curves to coalesce, but detailed
examination shows that there is an energy gap of about the same magnitude as
the cut-off point in a normal non-splined LJ potential (∼ 0.04 energy units),
meaning there is no crossing of the potentials. The current algorithm is applied
for both these cross-over regions with their different mechanisms of cross-over.
The MD cell is rectangular, with unit distance along the axis ( x direction)
of the cell length, whereas the breadth and height was both 1/16, implying a
thin pencil-like system where the thermostats were placed at the ends of the
MD cell, and the energy supplied per unit time step δt at both ends of the cell
(orthogonal to the x axis) in the vicinity of x = 0 and x = 1 maintained at
temperatures Th and Tl could be monitored, where this energy per unit step
time is respectively ǫh and ǫl. At equilibrium, (when Th = Tl), the net energy
supplied within statistical error (meaning 1-3 units of the standard error of the
ǫ distributions ) is zero, i.e. ǫl ≈ ǫh ≈ 0. The cell is divided up uniformly into
64 rectangular regions along the x axis and its thermodynamical properties of
temperature and pressure are probed. The resulting values of the ǫ’s and the
relative invariance of the pressure and temperature profiles would be a measure
of the accuracy of the algorithm from a thermodynamical point of view at the
steady state. For systems with a large number of particles such thermodynam-
ical criteria are appropriate. The synthetic thermostats now frequently used
in conjunction with ”‘non-Hamiltonian”’ MD [3] cannot be employed for this
type of study, where actual energy increments are sampled. The runs were for
4 million time steps, with averages taken over 100 dumps, where each variable
is sampled every 20 time steps. The final averages were over the 20-100 dump
values of averaged quantities.
0.8 1 1.2 1.4 1.6 1.8
r/LJ distance units
Potentials for simulation model
 intermolecular potential
s(r) switching function
atomic LJ potential
Figure 1: Potentials used for this work.
The temperature T and pressure p are computed by the equipartition and
0 20 40 60 80
x layer number
Figure 2: Temperature profile across the cell for different set conditions a−e for
temperature T ∗ and step time δt pairs (T ∗, δt) where a = (8.0, 2.0 ep− 3), b =
(8.0, 5.0 ep− 4), c = (8.0, 5.0 ep− 5), d = (12.0, 5.0 ep− 5), e = (16.0, 5.0 ep− 5).
The curves {l1, l3, t1, t2, t3} results with the application of the algorithm at rb
and rf with associated conditions l1 ⇔ a, l3 ⇔ b, t1 ⇔ c, t2 ⇔ d, t3 ⇔ e whilst
the curves {l2, l4, l5, l6, l7} are for the cases without implementing the algorithm
with the associated conditions l2 ⇔ a, l4 ⇔ b, l5 ⇔ c, l6 ⇔ d, l7 ⇔ e, where
ep x ≡ 10x.
0 20 40 60 80
x layer number
Figure 3: Pressure profile across the cell for different runs.The conditions of the
runs and the labeling of the curves are exactly as in Fig. (2).
Curve ǫh ǫl Mean Temperature
l1 -.2274E+00 ±0.19E-02 -.2295E+00 ±0.21E-02 0.9063E+01 ±0.62E-02
l2 -.5602E+00 ±0.22E-02 -.5596E+00 ±0.22E-02 0.1032E+02 ±0.63E-02
l3 -.4161E-01 ±0.14E-02 -.4089E-01 ±0.14E-02 0.8774E+01 ±0.79E-02
l4 -.5201E-01 ±0.16E-02 -.5103E-01 ±0.17E-02 0.8980E+01 ±0.98E-02
t1 -.5312E-03 ±0.92E-03 -.3334E-03 ±0.76E-03 0.8082E+01 ±0.49E-02
l5 0.1311E-02 ±0.82E-03 0.1147E-02 ±0.84E-03 0.7731E+01 ±0.97E-02
t2 -.6823E-03 ±0.12E-02 -.1507E-02 ±0.13E-02 0.1216E+02 ±0.17E-01
l6 0.7291E-02 ±0.12E-02 0.6343E-02 ±0.14E-02 0.1088E+02 ±0.15E-01
t3 -.9348E-03 ±0.18E-02 -.3379E-02 ±0.17E-02 0.1622E+02 ±0.18E-01
l7 0.1918E-01 ±0.14E-02 0.1938E-01 ±0.16E-02 0.1329E+02 ±0.20E-01
Table 1: Table with values for the mean heat supply per unit step and temper-
ature. The error is one unit of standard error for the quantities.
Virial expression given respectively by
pi.pi/mi
= 3NkBT andP = ρkBT +W/V
where W = − 1
w(rij) and the intermolecular pair Virial w(r) is given
by w(r) = r
dv(r)
with v being the potential.
2 ALGORITHMAND AND ANALYSIS OF NU-
MERICAL RESULTS
The velocity Verlet algorithm [4, p. 81]used here [1] and allied types generate
a trajectory at time nδt from that at (n− 1)δt with step increment δt through
a mapping Tm where (v(nδt), r(nδt)) = Tm(v((n − 1)δt), r((n − 1)δt)) which
does not scale linearly with δt. For a Hamiltonian H whose potential V is
dependent only on position r having momentum components pi, the system
without external perturbation has constant energy E, and the normal assump-
tion in MD (NAMD)is that for the nth step, ∆En = |H(nδt)− E| ≤ ǫ and also
i=1 ∆Ei ≤ ǫ
s for the specified ǫ′s. In the simulation under NAMD, the force
fields are constant and do not change for any one time step. In these cases, the
energy is a constant of the motion for any time interval δtT when no external
perturbations (e.g. due to thermostat interference) are impressed. When there
is a crossing of potentials at such a time interval interval from φb to φa at an
inter particle distance(icd ) rc (such as points rf and rb of Fig. (1)) of general
particle 1 and 2 (the (1, 2) particle pair) due to a reactive process (such as oc-
curs in either direction of (1.1)) a bifurcation occurs where the MD program
computes the next step coordinates as for the unreacted system (potential φb),
which needs to be corrected. Let the icd at time step i be ri (with φb potential)
and at step i + 1 after interval δt be rf = ri+1 where rf < rc < ri. Due to
this crossover, a different Hamiltonian H ′ is operative after point rc is crossed,
where under NAMD, the other coordinates not undergoing crossover are not
affected. For what follows, subscripts refer to the particle concerned. Let the
interparticle potential at rf be Ea = Ef = φa(rf ) and at rf be Eb = φb(rf ),
where ∆ = Eb − Ea. Then if rf be the final coordinate due to the φb potential
and force field, two questions may be asked: (i) Can the velocities of (1,2) be
scaled, so that there is no energy or momentum violation during the crossover
based on the φb trajectory calculation and (ii) Can a pseudo stochastic potential
be imposed from coordinates rc (at virtual time tc) to rf such that (i) above is
true? For (ii) we have
Theorem 2.1 A virtual potential which scales velocities to preserve momentum
and energy can be constructed about region rc.
proof The external work done δW on particles 1 and 2 over the time step is
proportional to the distance traveled since these forces are constant and so for
each of these particles i, Fext,i.∆ri = δWi where ∆ri is the distance increment
during at least part of the time step from rc to rf . For the non-reacting trajec-
tory over time λδt (λ ≤ 1) (virtual because it is not the correct path due to the
crossover at rc),
δW2 + δW1 − (φb(rf )− φb(rc)) = ∆
(K.E.) (2.1)
where ∆
(K.E.) is the change of kinetic energy for the (1, 2) pair from the
First Law between the end points rf , rc. Now over time interval tc to tf , for
the reactive trajectory, we introduce a ”‘virtual potential”’ V vir that will lead
to the same positional coordinates for the pair at the end of the time step with
different velocities than for the non-reactive transition leading to the transition
δW2 + δW1 − (V
vir(rf )− V
vir(rc)) = ∆
′(K.E.) (2.2)
where ∆
′(K.E.) is the change of kinetic energy for the pair with V vir turned
on and along this trajectory, the change of potential for V vir is equated to the
change in the K.E. of the pair as given in the results of theorem (2.2) for all
three orthogonal coordinates, i.e.
δV vir(r) − δφb(r) = δ
(K.E.x,y,z)−∆
′(K.E.)x,y,z
with momentum conservation, that is δV vir(ri) = δφa(ri) for the variation
along the ri coordinate, but δφa(ri) = −δK.E. along internuclear coordinate ri
whereas δV vir = −K.E. (scaled about all three axes). Continuity of potential
implies
φa(rf ) = V
vir(rf );φa(rc) = V
vir(rc);φb(rc) = V
vir(rc); (2.3)
Subtracting (2.1) from (2.2) and applying b.c.’s (2.6) leads to
∆ = φb(rf )− V
vir(rf ) = φb(rf )− φa(rf ) = Eb − Ea (2.4)
′(K.E.)−∆
(K.E.) (2.5)
The above shows that a conservative virtual potential could be said to be oper-
ating in the vicinity of the transition (from tc to ta) .•
Question (i) above leads to:
Theorem 2.2 Relative to the velocities at any rf due to the φb potential, the
rescaled velocities v ′ due to the potential difference ∆ leading to these final
velocities due to the virtual potential can have a form given by
′ = (1 + α)vi + β (2.6)
(where i = 1, 2) for a vector β.
proof From the v velocities at rf due to φb we compute the v
′ velocities at
rf due to the virtual potential. Since net change of momentum is due to the
external forces only, which is invariant for the (1, 2) pair, conservation of total
momentum relating v′ and v in (2.6) yields a definition of β ( summation from
1 to 2 for what follows, where the mass of particle i is mi)
β = −α
mivi/
mi (2.7)
Defining for any vector s2 = s.s,β2 = α2Q, where
(2.8)
then the rescaled velocities become from (2.6)
= (1 + α)2vi
2 + 2(1 + α)vi.β + β
2. (2.9)
With ∆ = Eb − Ea, Energy conservation implies
2 = ∆ (2.10)
The coupling of (2.9-2.10) leads, after several steps of algebra to
α2m1m2
2(m1 +m2)
2 + v2
2 − 2v1.v2
(2.11)
2αm1m2
2(m1 +m2)
2 + v2
2 − 2v1.v2
Defining a = (v1 − v2)
2, q = m2m1/[2(m1 +m2)], (q > 0, a ≥ 0), then the
above is equivalent to the quadratic equation
α2qa+ 2qaα−∆ = 0 (2.12)
and in simulations, only α is unknown and can be determined from (2.12) where
real solutions exist for ∆/qa ≥ −1. • The above Inequality leads to a certain
asymmetry concerning forward and backward reactions, even for reversible re-
actions where the region of formation and breakdown of molecules are located
in the same region with the reversal of the sign of approximate ∆. For this
simulation, a reaction in either direction (formation or breakdown of dimer )
proceeds if (??) is true; if not then the trajectory follows the one for the initial
trajectory without any reaction (i.e. no potential crossover).
Interpretation of results. Fig. (1) shows a rapidly changing potential curve
with several inflexion points used in the simulation at very high temperature
(as far as I know such ranges have not been reported in the literature for non-
synthetic methods) warranting smaller time steps; larger ones would introduce
errors due to the rapidly changing potential and high K.E.; thus, even with
the application of the algorithm between cordinates rf and rb, curves l1 and
l2 have too large a δt value to achieve equilibrium - meaning flat or invariant
- temperature (see Fig. (2) ) or pressure (see Fig. (3))or unit step thermostat
heat supply (see Table 1)(ǫh and ǫl) profiles where for these curves, the (ǫh, ǫl)
values show net heat absorption; the curve at t1 (with δt = 5.0 ep− 5 show flat
profiles (within statistical fluctuations and 2 standard errors of variation) for
temperature, pressure and net zero heat supply; and this choice of time step
interval was found adequate for runs at much higher temperatures (T = 12 and
T = 16) which was used to determine thermodynamical properties [2]. For this
δt value and all others, no reasonable stationary equilibrium conditions could
be obtained without the application of the algorithm (curves l2,l4,l5,l6 and l7).
The algorithm is seen to be effective over a wide temperature range for this
complex dimer reaction simulated under extreme values of thermodynamical
variables and the results here do not vary for longer runs and greater sampling
statistics (e.g. 6 or 10 million time steps). The thin, pencil-like geometry of the
rectangular cell with thermostats located at the ends would highlight the energy
non-conservation leading to a non-flat temperature distribution, as observed and
which was used to determine the regime of validity of the algorithm.
References
[1] J. M. Haile,Molecular Dynamics Simulation,JohnWiley & Sons,Inc.,New
York, 1992.
[2] C. G. Jesudason, Model hysteresis dimer molecule. I. Equilibrium prop-
erties. J. Math. Chem. JOMC, Accepted 2006.
[3] D. Frenkel and B. Smit, Understanding Molecular Simulations: From
Algorithms to Applications, Vol(1) of Computational Science Series, Aca-
demic Press, San Diego, Second Ed., 2002.
[4] M.P. Allen and D. J. Tildesley, Computer Simulation of Liquids,Oxford
Univ. Press, Oxford, 1992
	PRELIMINARIES
	ALGORITHM AND AND ANALYSIS OF NUMERICAL RESULTS
