Accurate energy conservation in molecular dynamics simulation Текст научной статьи по специальности «Физика»

ACCURATE ENERGY CONSERVATION IN MOLECULAR DYNAMICS SIMULATION

O.A. Zolotov, V. E. Zalizniak Siberian Federal University, 79 Svobodny Prospect, Krasnoyarsk 660041, Russia ozolot_@mail.ru, vzalizniak@sfU-kras.ru

PACS 02.70.Ns, 45.20.dh, 45.50.Jf

In molecular dynamics, Hamiltonian systems of differential equations are numerically integrated using some sym-plectic method. Symplectic integrators are simple algorithms that appear to be well-suited for large scale simulations. One feature of these simulations is that there is an unphysical drift in the energy of the system over long integration periods. A drift in the energy is more obvious when a relatively long time step is used. In this article, a special approach, based on symplectic discretization and momenta corrections, is presented. The proposed method conserves the total energy of the system over the interval of simulation for any acceptable time step. A new approach to perform a constant-temperature molecular dynamics simulation is also presented. Numerical experiments illustrating these approaches are described.

Keywords: Hamiltonian systems, symplectic numerical methods, energy conservation, molecular dynamics. 1. Introduction

In the field of molecular dynamics (MD), simulations of systems modeling materials or molecules at the microscopic scale are performed. There is a particular interest in large scale MD simulations, involving perhaps as many as several million atoms, over very long time scales. Such large simulations are of interest not only in standard macromolecular modeling but also in the modeling of various nanotechnology ideas. Molecular dynamics simulation first involves setting up the initial positions and velocities of the particles, and then the numerical integration of the classical equations of motion.

In a classical mechanics model, each atom is a particle with position rn, momentum pn and mass mn. Positions and momenta are determined by the equations of motion:

dpn _ f (t) drn _ Pn

dt n ' dt mn' rn(0), pn(0) are given, n _ 1,..., N,

where fn(t) is the force acting on the particle at time t and N is number of particles. Let us denote r _ (ri,..., rN) and p _ (pi,..., pN). Forces are determined from the system potential energy U(r):

dU (r)

fn(t)_ --u(-)_-Vr„U (r). O-n

These equations are a special case of a Hamiltonian system of differential equations:

dpn dH (r, p) drn dH (r, p)

dt d rn dt dpn

with Hamiltonian H _ K + U where K is the kinetic energy:

N n2

H (r p) = £ + U (r). (1)

Hamiltonian (1) is not explicitly time-dependent and its value is the total energy: H (r, p) = E. Classical dynamics is time reversible and the total energy of an isolated system is conserved. If the external forces on the individual particles sum to zero, the total momentum,

N

P = 5] Pn,

Jn 1

n=1

is conserved. If the external force acting on a system is zero, then the total angular momentum of the system,

N

L 'y ] rn X pn,

n=l

is constant.

Numerical integration is the most computationally-intensive part of an MD simulation. The desirable qualities for a simulation algorithm might be as follows:

(a) It should be fast, and require little memory. Numerical integration algorithms require a number of force evaluations in order to advance a trajectory by a given time step. The computational time is proportional to the number of force evaluations taken. Therefore, algorithms with one force evaluations per time step are preferable;

(b) It should allow the use of a long time step;

(c) It should reproduce the classical trajectory as closely as possible;

(d) It should satisfy the known conservation laws for energy, total momentum and angular momentum, and be time-reversible.

Symplectic schemes have been found to be effective for MD simulations with classical dynamics, especially for long time integration (see, for example, [1-3] and references therein). The most commonly used symplectic algorithm is the Verlet algorithm [4] which appears in many different formulations. In what follows, we will use a second order method, first disclosed by Rowlands [5]:

k n

rk+l/2 = rk + & rn rn + Tk 2m.

Pn+1 = Pn - Tk (Vr„U(r))

pk+i

rk+1 = rk+1/2 + pn

rn rn +Tk 2mn

n = 1,... ,N; k = 0,1,.

k + 1/2 ;

(2)

For symplectic algorithms, there exists a modified Hamiltonian Hm, which can be obtained from an asymptotic expansion of the original analytic Hamiltonian H. In the case of second order scheme, the general expression for the modified Hamiltonian is given in [6] and has the form:

Hm (r, p) = H (r, p) + T2g (r, p) + O (t4) .

r

n

n

Symplectic integrators do not in fact conserve the value of true Hamiltonian H; instead, they reproduce more closely the value of modified Hamiltonian Hm. The accuracy of the energy conservation for the true Hamiltonian is dependent on the time step.

In numerical solutions of classical dynamics, all variables are discretized, and this could destroy the energy conservation. In MD simulations, interactions beyond a cutoff distance rc are usually ignored. This means that the forces are not exactly correct. This also affects the energy conservation in MD simulations [7].

2. Constant energy simulation

Let us consider a system of N particles and U(r) is the potential energy of the system. Let E0 be the initial total energy of the system. Given coordinates and momenta at some time tk-1, (rk-1, pk-1), the one step given by equation (2) generates (rk, pk). Then, the total energy of the system E (rk, pk) = K (pk) + U (rk) = E (tk) and because the total energy is not conserved: E (tk) = E0. To achieve the energy conservation, we introduce corrections Apk to momenta pk with the proviso that

1 y ÍMA 2

2 ^ mn

n=1

— min.

Using the method of Lagrange multipliers, we introduce functional

W (Apk■ *) = 1 E ^ + * (1t ^ + U (rk) - E0

n=1 \ n=1

When variation SW = 0 for any arbitrary variations S (Apk) and SÀ, the functional W reaches its minimal value. Then, the condition SW = 0 leads to the following system of equations:

—Apn + A— (pn + Apn

mn mn

N („k )2 (3)

1 f ip4Ap^+U (rk )=E,.

2 m,

n=1

Solving the first of equations (3) we obtain:

Apn = -y+apn ^ pn + Apn = yt+aPn = apn.

The second equation (3) gives us the value of coefficient a:

1 f MiMl! + U (rk) = a2K (pk) + U (rk) = E0 a = ./E- Uf >.

2 ^ mn v ; VF ; V ; 0 y K (pk)

0

In summary, the above results lead to the following algorithm (ECI1 - energy conserving integrator 1):

„k+1/2

rn + Tk pn

2mn '

Pn

pn - Tk (Vr„U (r))

k+1/2 '

pfc+1 n

„k+1

a pn,

r

n = 1,

k+1/2 n

ak = N ;

+ Tk

2mn

(4)

'Ep - U (rk+1) K (p*) : k = 0,1,....

In algorithm (4) the difference of two terms E0 — U ~ K. In the case of a few body problems one can meet with the loss of numerical accuracy when kinetic energy K is small. Then we introduce corrections Apk to momenta pk and corrections Ark to positions rk with the proviso that

N (ApJ

— -i-

2

1

2

n=1

mn

2 k ^

+ -

,k)2

n=1

It is assumed that every particle n is connected to position rn through the mediation of virtual spring with the stiffness A particle may not be easily displaced from the position rn, which means that one should choose:

ks ~ max

n

d 2U

d r2

Using the method of Lagrange multipliers, we introduce functional

N

(Ap.

k 2 k N nk

s

n=1

m«

2

n=1

+ (Arn)2 + A( 2E

N ^ + Apn

+ U (rk + Ark) - E,

n=1

mn

The condition SW = 0, for any arbitrary variations S (Apk), S (Ark) and SA, leads to the following system of equations:

—Apn + A— (pn + Apn

nn mn mn

ksArn + A (VrnU(r)) )2

(5)

N ^n + Apn

E

n=1

mn

+ U (rk + Ark) = Ep

Let us use the following approximations:

(VrnU(r))|r„=rn+Arn « (VrnU(r))

N

U (rk + Ark) « U (rk) + J] (VrnU(r))

n=1

• Ark.

n

n — * n

n

2

0

0

1

2

rn=r

rn=r

Then, the second and third equations in (5) give us the cubic equation in one unknown parameter A:

2 x _ r "I 2

Eo - U (rk) .

1 \ 2 A N

1 K (pk) - A V [(V„U(r))

1 + A

n=1

We are looking for approximate solution of this equation assuming that A is close to zero because corrections should be relatively small. Then, we arrive at the following solution:

U (rk) + K (pk) - Eo A — -

2K (pk) + kS EN=1 [(Vr„U(r))|pn=p

Taking into account first and second equations in (5), we arrive at the following algorithm (ECI2 - energy conserving integrator 2):

„fc+l/2

rn + Tfc-pn

Pn

'2mn' pn - Tk (Vr„ U (r))

r n = p*

k+1/2 , _ pn

k+1/2 '

r

+ Tfc

2mn

„fc+i

rn - Ak (Vrn U (r)) |r =r

U (r*) + K (p*) - Eo

(6)

pfc+1

n

1 + kA

n = 1,...,N; k = 0,1,

2ksK (p*) + EN=i (VrnU(r)) 1 *

2.1. Harmonic oscillator

Consider the motion of a harmonic oscillator. The Hamiltonian equations in that case can be written as:

dp 2 dq (7)

dt— —wq, dt— p' (7) q(0), p(0) are given,

where p is the momentum, q — mx is the canonical coordinate, w — (k/m)1/2, m is the mass and k is the stiffness. The total energy of the harmonic oscillator is E — K(p) + U(q) — 0.5(p2 + q2)/m. The exact solution of the problem (7) is:

cos (wt) w-1 sin (wt) q (0) —w sin (wt) cos (wt) p (0)

p (t)_

Application of symplectic method (2) to the harmonic oscillator problem gives the following scheme:

qk+i/2 — qfc + 0.5Tpk, pk+1 — pk — Tw2qk+1/2,

qk+1 — qk+1/2 + 0.5rpk+1, (8)

k — 0,1,....

This scheme is stable if tw < 2.

Figures 1-3 show the phase trajectory and time dependence of the position of the harmonic oscillator, as calculated with the use of methods (4), (6) and (8). Initial conditions

2

n

rn=r

are q (0) = 0 and p (0) = 1. The oscillator parameters are m = 1 and k = 1. A relatively long time step was chosen with rw = 0.25n (eight time steps per oscillation period). To estimate the energy conservation, we introduced the relative energy error as:

1 M

(AE ) = eTgtm D Ie (tk) - E (0)l, (9)

V ' k=1

where E(0) is the initial total energy and M is the number of time steps. Fig. 1(a) indicates that the total energy is conserved, the relative energy error (AE) = 3.7 ■ 10_17, while Fig. 1(b) indicates that the position of harmonic oscillator slightly deviates from the exact solution.

Fig. 1. Phase trajectory and time dependence of the position of harmonic oscillator, method (4); a) Phase trajectory in the phase space q—p, solid line is the exact trajectory, dots indicate approximate solution; b) position of harmonic oscillator versus time, solid line is the exact solution, dashed line indicates approximate solution.

Fig. 2. Phase trajectory and time dependence of the position of harmonic oscillator, method (6) with = 5; a) Phase trajectory in the phase space q — p, solid line is the exact trajectory, dots indicate approximate solution; b) position of harmonic oscillator versus time, solid line is the exact solution, dashed line indicates approximate solution.

a) i b) í

Fig. 3. Phase trajectory and time dependence of the position of harmonic oscillator, symplectic method (8); a) Phase trajectory in the phase space q — p, solid line is the exact trajectory, dots indicate approximate solution; b) position of harmonic oscillator versus time, solid line is the exact solution, dashed line indicates approximate solution.

As is seen from Fig. 2(a), the approximate phase trajectory slightly deviates from the exact one. Method (6) conserves the total energy much worse than method (4), the relative energy error being (AE) = 7.9 ■ 10_3, but provides similar results to method (4) for the time dependence of the harmonic oscillator's position, as shown in Fig. 2(b). Symlectic method (8) does not conserve the total energy, but rather the value of the modified Hamiltonian, as can be seen from Fig. 3(a). In this example, the relative energy error was (AE) = 7.7 ■ 10_2. The discrepancy between the calculated and exact position of the harmonic oscillator grows more quickly with time in comparison to methods (4) and (6).

The presented results demonstrated the accuracy and efficiency of the proposed method for a simple, analytically solvable, test case. However, actual MD simulations involve considerably more complicated dynamics, and it is therefore worthwhile to examine a more complex example.

2.2. MD simulation of the Lennard-Jones fluid

Let us analyze the time dependence of the energy for a fluid of N particles interacting via the well-known Lennard-Jones (LJ) pair potential:

, i d\12 /dx 6 (r) = 4e

where r is the distance between particles, e and d are potential parameters. Particles are contained in a box with rigid walls and V is the volume of the box. The potential energy of the system is given by:

1 N N

U (r) 2 ^ y ^ ^ ^ (rnm) , rnm |rn rm| .

n=1 m=1 m=n

In order to choose the time step of the simulation, it is necessary to determine the highest oscillating frequency wmax in the physical system and satisfy the condition for stability, Twmax < 2. The second condition to be satisfied is that the integration of the trajectories for

most particles should be sufficiently accurate. Most particles will have an oscillation frequency wa associated with the average separation ra = (V/N)1/3. The oscillation frequency estimate can be obtained with a simple one-dimensional argument considering the case when two neighboring particles move towards each other to a position of closest approach [8]. Then, we propose the following estimates that are sufficient for the choice of a reasonable time step:

^ = (-^ irk -k = (!_ dV (r )\

-max Vmp dr2 ' Vmp dr2 (r"V '

where mp is the minimal reduced mass and rm 1n is the closest approach at some time tk. They are:

. f TOiffl^ k . ( k X

mp = min - , rm1n = mm (rj .

p Vmi + mj J min v j

Another constraint on the time step arises ^om the requirement that a particle should not move more than the minimal distance between particles. Hence, the condition 2rk^max < rm.n should be satisfied, where ^max is maximal particle velocity at some time tk defined as:

Vrnax = max( — (pn ■ Pn)1/2) .

n \m/ J

The condition on the time step may be summarized as follows:

/ k \ 1 M

. / ai a2 a3rm1n \ 1 V^

Tk = min -j—— > t =J7 VTk.

V—max Vraaw M ^

Parameters a1, a2 and a3 define the accuracy of the integration of particle motion, t is the average time step and M is the number of time steps. Because the probability of small separations occurring between particles is small, it is not necessary to accurately integrate the motion of such particles. The trajectories for most particles should be integrated more accurately. Then, we have 2 > a1 > a2 and parameter a3 < 0.1.

We considered a fictitious fluid with potential parameters e = 0.015 eV and d =3 A. All particles have equal mass m =20 u. Then, an atomic time unit = (u ■ A2/eV)1/2 corresponds to 1.019■ 10-14 s. We simulated N = 100 LJ particles at reduced density p* = 0.728 and reduced temperature T* = 1.376, a moderate-pressure liquid state. In MD simulations, interactions beyond a certain distance rc are usually ignored. In what follows, rc is taken to be equal to 2.5d. To obtain initial conditions, we used simplectic method (2), with t « 0.001, and integrated the equations of motion for 105 time steps, saving the resulting coordinates and momenta. These coordinates and momenta were then used as initial conditions for all simulations. The resulting value of the total energy of the system was taken as the initial value of energy E(0) that must be conserved. The dynamics were simulated in double-precision arithmetic.

The parameters of three simulations are presented in Table 1.

Table 1. Parameters of simulations

Method ai a3 T, atomic units M

1 Symplectic (2) 0.025n 6.25-10-3n 0.1 0.0081 4105

2 Symplectic (2) 0.5n 0.125n 0.1 0.161 2104

3 ECI1 (4) 0.5n 0.125n 0.1 0.162 2104

In all simulations, we considered the relative energy error (9) and some basic thermo-dynamic quantities, such as, temperature T, pressure P and constant-volume heat capacity CV. These thermodynamic quantities can be calculated as averages in a microcanonical ensemble as follows [9]:

N N

T = (K) P = NkBT + ± yyr d*(r )

T 3NkB (K) , P V + 3V \ ^ rnm dr (rnm)

\n=1 m>n

3 ( 3 ((SK)2nX-1

cv = 3 NkB ^ - ^ ^

where angle brackets (■) mean the time average, SK (tk) = K (tk) — (K) is the fluctuation of kinetic energy at time tk and kB is the Boltzmann constant. Calculated thermodynamic quantities are presented in Table 2.

Table 2. Calculated properties for a LJ fluid

(AE > T, K P, eV/A3 CV, eV/K

1 9.610-3 240.3 9.310-4 0.0174

2 0.12 249.2 10-3 0.0173

3 1.110-15 239.2 9.510-4 0.0175

The dependence of the energy of the system of LJ particles on time for three simulations presented in Table 1 is shown in Figs. 4-6. Symplectic methods do not conserve energy exactly along trajectories. However, although they do not conserve energy, symplectic methods have been observed to maintain system energy in a narrow band near the true energy for long periods of time (see Fig. 4). That is, the energy of the system fluctuates about some value with very little apparent long-term drift. However, over longer time periods there is a slow drift in the energy away from this range. If relatively long time step is used then a drift in the energy becomes even more apparent (see Fig. 5). We considered the system of particles contained in a box with rigid walls. In this case, the integration of the motion of the system with the use of symplectic method (2) resulted in an increase of the system's kinetic energy. The proposed method (4) conserved the total energy of the system over the interval of simulation (see Fig. 6); the relative energy error is close to round-off errors.

The average deviation (|1 — ak|) is equal to 4.8 ■ 10_5 for simulation 3 (ak is momentum correction parameter from (4)), so very small correction of momenta is needed to conserve the total energy.

3. Constant temperature simulation

The definition of thermodynamic temperature can be obtained through the fundamental equation of state [10], from which we get:

1 - dS

T = dE,

where S is the entropy and T is the thermodynamic temperature recorded by a thermometer in thermal equilibrium with the system. The kinetic temperature of the system is related to system kinetic energy through the standard relationship:

Fig. 4. The evolution of the energy for a LJ fluid, simulation 1.

Fig. 5. The evolution of the energy for a LJ fluid, simulation 2. Dashed line indicates the initial total energy.

3

2 NkBTk = <K>. (10)

For a system in equilibrium, it is well known that the kinetic temperature is identical to the thermodynamic temperature, i.e. T = Tk. With the use of the maximum entropy formalism, it was found that even though a system is far enough from equilibrium for nonlinear effects, Tk/T - 1 < 0.01 [11]. In what follows, we assume that T ~ Tk.

Simulations at constant temperature are important for studying the behavior of systems at different temperatures. Several different methods of prescribing the temperature in a molecular

Fig. 6. The evolution of the energy for a LJ fluid, simulation 3.

dynamics simulation exist [9]. We propose a new approach to adapt MD so as to sample a constant-temperature ensemble.

Let us consider a system of N particles and T0 is the prescribed temperature that corresponds to the kinetic energy K0 through relation (10). Given coordinates and momenta at some time ik-1, (rfc-1, pfc-1), the one step given by equation (2) generates (rfc, pfc). Then, the kinetic energy of the system K (pfc) and particle kinetic energy Kn (p^) at time tk are defined as:

N k k

K (pk ) = £ pb^, K (p£

n=1

2mn

pra pra

2m„

The kinetic energy (temperature) is not conserved: K (pk) = K0. Let us introduce kinetic energy difference per particle as:

ak (pk) = N (K - K (pk))

To achieve the conservation of kinetic energy, we introduce corrections Apk to momenta pk with the proviso that K (pk + Apk) = K0. We assume that:

pn+Apn=pn+pn

(1 + Yn) pn.

Then, we have:

Kn (pn + Apn) - Kn (pn) = AK (pk) = (1 + 7n)2 Kn (pn) - Kn (pn). Solving the last equation and assuming that y is a real and positive parameter, we obtain:

1 + Yn (p^)

(1

+ AK (pk) /Kra (pn))1/2,

if AK (pk) /Kra (p^) > -1; otherwise.

In summary, the above results lead to the following algorithm (KECI - kinetic energy conserving integrator):

„fc+1/2

rn + rk

*

P n

Pn - Tk (Vr

„fc+1

r

fc+1/2

pn

2 mn' U(r)) '

+ Tfc-

1 n

(11)

pfc+1

n

1,.

2mn'

(1 + Yn (pn)) pn,

n = 1,...,N; k = 0,1,....

Let us consider a system of N particles interacting via the Lennard-Jones pair potential defined in section 2.2. The target temperature T0 is set to 300 K. The system was integrated with the use of method (11) for 2-104 time steps with the average time step t = 0.152. Fig. 7 shows the evolution of the instantaneous temperature deviation (T0 — T) /T0 over the interval of the simulation. During the simulation, the speed distribution histogram was constructed. As can be seen from Fig. 8, the obtained speed distribution of the particles is close to the Maxwell distribution. Pearson's chi-squared test confirms this observation.

Fig. 7. The evolution of relative temperature deviation for a LJ fluid.

n

n

4. Conclusions

Approximate solutions for the equations of motion in MD simulations are typically computed by some symplectic method, however any known symplectic numerical method does not conserve energy exactly along trajectories. This is because the various types of errors have an effect on the accuracy of the computation. Our investigation showed that the energy can be accurately conserved by introducing relatively small corrections to the particle momenta. These corrections bring the system back to a statistically-admissible state on the constant energy surface in a phase space.

Fig. 8. Speed distribution. Histogram is obtained from simulation; dashed line corresponds to the Maxwell distribution.

The proposed kinetic energy conserving integrator is simple to implement and is very efficient for maintaining the desired temperature. It also does not depend upon any additional parameters. This method produces a statistical ensemble that is close to canonical ensemble.

References

[1] C.J. Budd, M.D. Piggott. Geometric integration and its applications. Handbook of Numerical Analysis, XI, North-Holland, Amsterdam, P. 35-139 (2003).

[2] E. Hairer, C. Lubich, G. Wanner. Geometric Numerical Integration. Springer, Berlin, 2002, 644 p.

[3] J.M. Sanz-Serna, M.P. Calvo. Numerical Hamiltonian Problems. Chapman and Hall, London, 1994, 207 p.

[4] L. Verlet. Computer 'experiments' on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules. Phys. Rev., 159, P. 98-103 (1967).

[5] G. Rowlands. A numerical algorithm for Hamiltonian systems. J. Comput. Phys., 97, P. 235-239 (1991).

[6] J. Gans, D. Shalloway. Shadow mass and the relationship between velocity and momentum in symplectic numerical integration. Phys. Rev. E, 61, P. 4587-4592 (2000).

[7] S. Toxvaerd, O.J. Heilmann, J.C. Dyre. Energy conservation in molecular dynamics simulations of classical systems. J. Chem. Phys., 136, P. 224106(8) (2012).

[8] R.W. Hockney, J.W. Eastwood. Computer Simulation Using Particles. Adam Hilger, Bristol and New York, 1988, 523 p.

[9] M.P. Allen, D.J. Tildesley. Computer Simulation of Liquids. Clarendon Press, Oxford, 1989, 385 p.

[10] L.D. Landau, E.M. Lifshitz. Statistical Physics, Part 1 (Course of Theoretical Physics, Volume 5), Third Edition, Elsevier Butterworth Heinemann, Oxford, 1980, 544 p.

[11] R.E. Nettleton. On the relation between thermodynamic temperature and kinetic energy per particle. Canadian Journal of Physics, 72(3-4), P. 106-112 (1994).