Numerical irreversibility in time-reversible molecular dynamics simulation

Physica D 195 (2004) 391–397

Numerical irreversibility in time-reversible molecular dynamics simulation Nobuyoshi Komatsua,∗ , Takashi Abea,b a

Department of Aeronautics and Astronautics, University of Tokyo, Hongo 7-3-1, Bunkyo-ku, Tokyo, Japan b The Institute of Space and Astronautical Science, Yoshinodai 3-1-1, Sagamihara, Kanagawa, Japan Received 6 January 2004; received in revised form 11 May 2004; accepted 11 May 2004 Communicated by V. Frisch

Abstract It is demonstrated that in the molecular dynamics simulation for the N-body system, even a negligibly small deliberate noise can cause irreversibility, or the loss of reversibility. This demonstration is unique and clear in that it is based on the “bit-reversible algorithm” which is, unlike the standard molecular dynamics simulation being unable to avoid round-off errors, completely time-reversible and is free from any round-off error. It is also confirmed that the irreversibility is closely related to the extent of the instability of the system, in that the more unstable the system is, the more quickly the irreversibility appears and the more rapidly the irreversibility prevails. © 2004 Elsevier B.V. All rights reserved. PACS: 02.70.Ns; 05.70.Ln; 05.40.Ca; 05.45.Pq Keywords: Molecular dynamics; Round-off errors; Irreversibility

1. Introduction Why does a time-asymmetric (or irreversible) behavior of macroscopic systems arise from a time-symmetric (or reversible) microscopic law? To illustrate this question more concretely, let us consider, as a simple example, a system consisting of a box ∗ Corresponding author. Tel.: +81 42 759 8255; fax: +81 42 759 8460. E-mail addresses: [email protected] (N. Komatsu); [email protected] (T. Abe).

0167-2789/$ – see front matter © 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.physd.2004.05.004

divided into two regions by a partition. Suppose that, initially, the left half region is filled with gas composed of molecules and the right half one is empty. Based on our experience, as soon as the partition is removed, the gas expands into the whole container and the molecules never go back to the left half region of the box, although such a motion would be possible because of the time-symmetric microscopic law. In the middle of 19th century, Boltzmann [1] tried to verify this irreversible process of macroscopic systems using his elegant technique based on the kinetic equation and H-theorem. But his approach drew two major

392

N. Komatsu, T. Abe / Physica D 195 (2004) 391–397

criticisms. The first one was Loschmidt’s reversibility paradox [2] based on time-reversible dynamics and the second one was Zermelo’s recurrence paradox [3] based on Poincar´e’s recurrence theorem. Therefore, the original question, on the time-asymmetric behavior of macroscopic systems arising from the time-symmetric microscopic law, remains unanswered [4,5]. To confirm the Boltzmann’s hypothesis, the N-body problem must be solved. On the other hand, as Poincar´e proved, the N-body problem, in general, cannot be solved analytically. To resolve this difficulty, in the 1950s, molecular dynamics (MD) methods using classical mechanics were devised, which were quite useful to simulate a non-equilibrium process because, in the MD methods, a time evolution of microscopic process was calculated based on the deterministic Newton’s law. For example, Loschmidt’s reversibility paradox was investigated by Orban and Bellemans [6] in 1967 by using the hard disk MD method. They considered a system of interactive molecules, stating from a nonequilibrium state. In accordance with the Boltzmann’s hypothesis, they observed the evolution of Boltzmann’s H-function H(t) in the system during a certain time interval trev , and confirmed that the value of H(t) decreased with time t and approached to the equilibrium state realizing the Maxwell distribution. Then they suddenly reversed all velocities of the molecules at time t = trev , forcing the system to return to the initial non-equilibrium state and to exhibit an ‘anti-kinetic’ behavior [6]: the value of H(t) increases with time. They found that the value of H(t) tended to return to the initial one but not exactly, and this tendency was amplified as trev increased and by adding some noise at the reversal. Based on these findings, they claimed that (1) the unavoidable round-off error is the cause of the failure of the complete regain of H(t) and (2) anti-kinetic situations present some kind of instability and, due to this instability, the anti-kinetic behavior of the system tends to be reduced or even suppressed. Their first remark was claimed repeatedly by many researchers [5,7,8], while the second one has been extensively studied in relation to the chaotic behavior of the system [5,9–15]. However, since the MD simulation was not able to avoid the round-off error by nature, they could not succeed to verify quantitatively how round-off errors had such an effect on the MD simulation. Furthermore, they even could not succeed to verify the complete regain of H(t) when any round-off error was avoided.

In 1993, Levesque and Verlet succeeded in verifying the complete regain of H(t) based on their newly developed ‘bit-reversible algorithm’ (Bit MD) [8,16]. Just as the standard MD method, in the ‘bit-reversible algorithm’, the time evolution of the microscopic process is calculated based on Newton’s law. The main exception is that the locations of the particles are discretized in the ‘bit-reversible algorithm’. However, since the interval of the discretization is small enough to be neglected in a physical point of view, both of them are sufficiently equivalent. The ‘bit-reversible algorithm’ is an ideal time-reversible algorithm, in that it can avoid round-off errors completely. This success means that, if a quantitatively-controlled noise is added to the MD method by the bit-reversible algorithm, it may be possible to quantitatively investigate numerical irreversibility which may arise. Accordingly, in the present paper, both irreversibility and stability in many-particle systems are investigated, using the Bit MD added with the quantitativelycontrolled noise. First, in Section 2, we give a brief review for the bit-reversible algorithm in order to make the present paper self-consistent. In Section 3, we describe the time-reversible simulation method in order to study irreversibility and stability. In Section 4, based on Bit MD with the quantitatively-controlled noise, we study behavior of both the loss of reversibility which is caused by the controlled noise and the stability of the systems.

2. Bit-reversible algorithm In the classical molecular dynamics technique, a time-reversible algorithm that has been widely used is the Verlet algorithm, which employs a central difference scheme for time derivatives in Newton’s equation. Although the scheme itself is timereversible, exact time-reversibility of its algorithm is violated. The cause of this is believed to be numerical round-off errors which are unavoidable in the standard MD based on a floating-point arithmetic (Float MD). On the other hand, the bit-reversible algorithm, which was developed by Levesque and Verlet, was demonstrated to be completely time-reversible [8]. The bit-reversible algorithm employs the Verlet algorithm for time derivatives but, unlike the standard Float MD

N. Komatsu, T. Abe / Physica D 195 (2004) 391–397

technique, employs a discrete coordinate space, instead of a continuous coordinate space. Since the bit-reversible algorithm may not be wellknown very much, we will give a brief review for this algorithm in order to make this paper self-consistent. Now, we will consider a system of N particles of mass m enclosed in a cubical box of size L. First, the minimum lattice distance L (= L/M), by which the coordinate space is discretized, is defined based on the side length of the cubical box L and some integer value M. Although we can choose any integer value M, the maximum value is 2n , if n-bit integers are employed. Because of the discretization, the position of a particle in the discrete coordinate system is represented by integers. In this context, the equation of motion of the particles in the discrete coordinate system can be written as ∂ 2 Xi (t)2 = Xi (t + t) − 2Xi (t) + Xi (t − t) ∂t 2  
 fij (t)  (t)2 = (1) m L j

Integer

where fij (t) is a partial force from the j-th particle on the i-th molecule, located at position Xi at a time t. It should be noted that the location of the i-th particle in the discrete coordinate space Xi is an integer although, in the physical space, its position is a real number and designated as Xi L. The term in the summation on the right hand side of Eq. (1) is calculated based on the particle locations, and the value in the bracket [·]Integer is converted to an integer since, in general, the value calculated in this way is not in integer but in real number. Because of both the space discretization and this conversion process to integer, the bit-reversible algorithm keeps the time-reversibility [8,16,17]. In the Bit MD simulation, the total momentum is exactly conserved but, due to the conversion process to integer, the energy is not exactly conserved. Nevertheless, the exact time-reversibility of the integration algorithm precludes any systematic drift in the total energy [8]. Through simulations in the present study, we have confirmed that the fluctuation of the total energy by Bit MD simulations is equivalent to that by Float MD simulations.

393

3. Time-reversible simulation In the present simulation, we consider a system of N particles of mass m enclosed in a closed region, and assume the two-dimensional system for simplicity. Hence, for the closed region, we consider a square cell whose side length is L, assuming periodic boundary conditions. The particles interact through the Lennard– Jones potential:  12 6 σ  4ε σ , (r ≤ 21/6 σ) − r r P(r) = (2)  1/6 0, (r > 2 σ) where only the repulsive part of the potential with a cut-off radius rc of 21/6 σ is taken into account. Accordingly, in general, the system considered in the present study is unstable as pointed out by Dellago and Posch [15]. Hereafter, we employ σ, ε, m, ε/κB , σ(m/ε)1/2 , as units of length, energy, mass, temperature, and time, respectively. Here κB is Boltzmann’s constant. Initially, in the square cell (L = 80), all the particles are distributed on all the cross points of mesh grids of size 2, and are assumed to have a velocity equal in magnitude V0 (= 2.0) but with a random direction, keeping the total momentum 0. Hence, all together, 1600 particles are considered in the square cell, and the number density ρ is equal to 0.25. The time increment t for integral of the equation of motion is set to be 0.01. In the present Bit MD simulation, is employed 30-bit integers to represent the discrete coordinate of the particles. That is, the side of the square cell is divided into 230 discrete segments and the particles move around in the discrete coordinate space discretized by this segment length. It should be noted that the segment length L is small enough in that it equals to 7.5 × 10−8 in contrast to a particle diameter d of 1. The system of particles in the square cell is, initially, in a highly non-equilibrium state and temporally evolves towards an equilibrium state. At a certain time trev in the middle of the equilibrium state, is executed the time-reversal operation in which all the particles reverse their velocity instantaneously. Therefore, the system evolves reversibly from then and, if the system is reversible, the initial state appears again after a passage of time which is equal to a time interval from 0 to trev . On the other hand, if the system contains any irreversibility, it is expected that the initial state ap-

394

N. Komatsu, T. Abe / Physica D 195 (2004) 391–397

pears only incompletely or does not appear anymore, depending on the extent of the irreversibility. In order to make such a non-equilibrium behavior visible, we employ Boltzmann’s H-function

∞ H(t) ≡ f (t) ln f (t) dv (3) −∞

where f (t) is a probability distribution for the particle velocity v at a time t. According to Boltzmann’s Htheorem, the H-function, in the isolated systems, varies monotonically from the initial value towards its thermal equilibrium value just as entropy increases. For comparison with the bit-reversible MD method, the standard MD method (Float MD) is employed. In Float MD, a single precision floating point real number is employed, which has a precision equivalent to that of the 30-bit integers Bit MD in that the minimum value in the floating point real number is equivalent to L or 7.5 × 10−8 . Fig. 1 shows the time evolution of the H-function for the system when the time-reversal operation is applied at various timings designated by a to d which corresponds to trev = 3, 4, 5, and 6, respectively. As expected, in Bit MD, the H-functions return to the initial value exactly after a passage of time, trev , from the time-reversible operations. That is, at t = 2trev , the initial state appears again. On the other hand, in Float MD, the time evolution of the H-function is strongly influenced by the timings of the time-reversal operation and, in any case, the H-functions do not go back to the initial state completely. In fact, the degree of the recovery to the initial state, which can be evaluated by the recovered value of the H-function or the recovery rate, RR , defined by, RR ≡

dH H

(4)

becomes worse for the larger value of the timing, trev , of the time-reversal operation. Here dH and H are defined in Fig. 1(b). This behavior is the one pointed out to be due to round-off errors by Orban and Bellemans [6]. If this is the case, we can expect that adding some numerical noise to the bit-reversible simulation may induce a similar irreversibility and that such an investigation may clarify the root cause of the irreversibility observed in Float MD. The noise added in the following simulation is a deliberate displacement of the particles. That is, the parti-

Fig. 1. Time evolution of the H-function obtained by time-reversible simulations of 1600-soft-disc system. Time-reversal operations are executed at the timings designated by a, b, c and d. dH and H at trev = b are also defined for Eq. (4). (a) Bit MD (30-bit integers). (b) Float MD (single precision: 32-bit).

cles change their positions instantaneously at a certain time by a certain amount. Since the noise should be as small as possible for the present purpose to demonstrate the irreversibility caused by the numerical noise, an arbitrary single particle out of the whole particles is chosen to be displaced into an arbitrary direction by a minimum lattice distance, L, in the discrete coordinate-space, only once at the timing of the timereversal operation. In the present study, we call it the ‘minimum controlled noise’. Since the irreversible process and the stability of the system are expected to be closely related, we will examine not only the recovery rate RR but also the extent of the instability included in the present system as well. For this purpose, the maximum Lyapunov exponent based on the chaos theory is useful. The maximum Lyapunov exponent, which can evaluate the extent of the instability inherent to the dynamical systems, indi-

N. Komatsu, T. Abe / Physica D 195 (2004) 391–397

395

zero [11,12], since the present systems are the Hamilton system at a constant energy in which the equation of motion is time-reversible. In this context, the maximum Lyapunov exponent obtained from this method is equal to the maximum Lyapunov exponent of the ‘original system’, i.e. before the time-reversal operation. To compute the maximum Lyapunov exponent based on those two trajectories, we employ the following equation proposed in [18] for simplicity, although several methods have been proposed [11], λ≡ Fig. 2. Sketch of trajectory distance. The origin of ti is the timing trev of the time-reversal operation, and the minimum controlled noise is added to an arbitrary single particle only once at t = trev .

cates a separation rate of a distance between two nearby trajectories of the system in the phase space. To realize such a nearby trajectory, the small displacement of a particle, which is intended to be a numerical noise, can be made use of. That is, the temporal evolution of the system with this small displacement of a particle gives the nearby trajectory adjacent to the original trajectory. In the present study, the maximum Lyapunov exponent is computed as follows. As shown in Fig. 2, the original trajectory is taken as the trajectory before the time-reversal operation, i.e. from t = 0 to t = trev . Then, at t = trev , the minimum controlled noise is added to the system. Also, the nearby trajectory is taken as the trajectory after the time-reversal operation, i.e. from t = trev to t = 2trev . Mathematically, the systems before and after the time-reversal operation can be discriminated as the ‘original system’ and the ‘reverse system’, respectively. It should be noted that, as shown in Fig. 2, the origin of the time for the computation is taken as the timing of t = trev , and the temporal integration is calculated from t = trev to t = 2trev . Accordingly, the temporal integration of the original trajectory is in the reversed direction comparing to the standard definition. However, since the trajectory is time-reversible because of the time-reversibility of the equation of motion, this definition is consistent with the standard definition. Therefore, from the present method, we obtain the maximum Lyapunov exponent of the ‘reverse system’, i.e. the system after the time-reversal operation. On the other hands, the spectrum of Lyapunov exponents of the present systems is symmetrical around

n−1 1
S(ti+1 ) ln nt S(ti )

(5)

i=0

where S(ti ) is the distance between the original trajectory and the nearby trajectory at a time ti , and subscript i represents each time step of the time increment t. It should be noted that, in Eq. 5, the origin of ti is the timing trev of the time-reversal operation. The validity of such a simplified method has been confirmed through a comparison with the result presented in [15].

4. Irreversibility of the system and its stability In the following simulations, the time-reversal operations are executed at several timings trev ranging from 2 to 15, and the recovery rate RR is determined for each time-reversal operation. It should be noted that, if the system is reversible, the system goes back completely to the initial state after a passage of time, trev , from the time-reversal operation, or the recovery rate RR becomes 100%. As shown in Fig. 3(a) where the recovery late RR is depicted versus trev , the recovery rate of the H-function RR is strongly influenced by even the weakest noise assumed in the present computation when the time trev becomes larger. That is, the recovery rate deteriorates when the timing for the time-reversal operation becomes late. Since the larger the time trev is, the more time is require for the system to recover to the initial state, it can be concluded that the more time is required for the system to recover to the initial state, the more significantly the recovery rate deteriorates. This deterioration of the recovery rate RR means that the added noise as small as the one assumed in the present paper is strong enough to bring some irreversibility to the system when the time trev becomes larger. This is an extraordinary result since the weak-

396

N. Komatsu, T. Abe / Physica D 195 (2004) 391–397

Fig. 4. Effects of collision frequency on the recovery rate of Hfunction RR . The elapsed time is normalized by the mean collision time τ using Eq. (6).

Fig. 3. Effects of temperature. (a) Recovery rate of H-function RR . (b) Maximum Lyapunov exponent λ.

est noise induced in the simulation is considered to be almost negligible since the displacement is as small as 7.5 × 10−8 times of the particle diameter and is assumed for a single particle only once at the timereversal operation. The behavior of the recovery rate of the H-function, however, is influenced also by the system parameter such as temperature. To see the effect of the temperature, the simulations with various initial particle speed V0 are carried out since the temperature T can be defined as V 2 /2 in terms of the averaged particle speed V. According to Fig. 3(a), the higher the temperature is, the more quickly the recovery rate starts to deteriorate and the more rapidly the recovery rate deteriorates. That is, the higher the temperature is, the more quickly

appears the irreversibility and the more rapidly prevails the irreversibility. The maximum Lyapunov exponents are positive in all the conditions and increases with increasing temperature as shown in Fig. 3(b), where Trev is the temperature based on the averaged particle speed V at t = trev . In fact, λ is nearly proportional to the averaged particle speed. That is, the system is unstable and the extent of the instability increases with the temperature. Therefore, it can be concluded that, in the unstable system such as the present one, even a negligibly small noise can cause a significant irreversibility and that the more unstable the system is, the more quickly appears the irreversibility and the more rapidly prevails the irreversibility. The maximum Lyapunov exponent has been intensively investigated by using the standard MD method for various dynamical systems [15]. The behavior of the maximum Lyapunov exponent obtained in the present study is quite consistent with their results. The behavior of the irreversibility and the extent of the instability appearing in the system should be root-caused by the interactions between the particles or the particle thermal collisions. Therefore, it can be expected that their temporal behavior may be correlated with the mean collision time. For a low density system, the following elementary definition for a hard sphere

N. Komatsu, T. Abe / Physica D 195 (2004) 391–397

model can be applied, τ=

1 √ 2σ0 Trev ρ

(6)

where τ is a collision time, and σ0 is a collision crosssection area which is equal to 2d in the two-dimensional analysis. Even though, in the present system, the number density ρ (= 0.25) is rather large for application of Eq. (6), the equation may still describe qualitatively the effect of mean collision time. As expected, when we replot the recovery rate against the time normalized by the mean collision time, all the results are quite consistent with each other as shown in Fig. 4. This suggests behavior of the irreversibility, or the loss of reversibility, in the present system, is root-caused by the thermal collisions between particles, and that the irreversibility appears as a propagation of the effect of noise through the thermal collisions between particles.

397

which combines the classical Newton’s law with the quantum random force. That is, his approach to the quantum mechanics suggests that there is a possibility for the classical dynamics to be applicable not only to the macroscopic system but also to the microscopic system when an appropriate quantum noise is assumed. According to his theory, the uncertainty of the particle trajectory can be estimated as, 
δx ∼ δt (7) 2m where δt is a time increment for Nelson’s equation and δx is a quantum noise corresponding to δt. As a matter of fact, Eq. (7) gives 4.1 × 10−12 m when we assume Argon and δt = t, while the magnitude of the present numerical noise, L, is as small as 2.5 × 10−17 m for Argon. That is, the present numerical noise is sufficiently smaller than the quantum noise. This suggests that the irreversibility in the present study is not caused by unphysically large numerical noise.

5. Discussion As revealed above, the numerical irreversibility, which is similar to the effect of round-off errors in the standard MD, does arise from the deliberate numerical noise even if the magnitude of the noise is extremely small, and the appearance of the irreversibility is related to the fact that the system is unstable in terms of the maximum Lyapunov exponent. However, there is still a doubt that this finding is based on the unphysically large noise. Therefore, it is useful to confirm that the irreversibility in the present study is not caused by unphysically large numerical noise. First of all, we must recall that the present approach, which is based on the application of classical mechanics to the microscopic process, has an essential difficulty because, in general, microscopic laws are based on quantum mechanics. In quantum mechanics, however, an observable in the Schr¨odinger equation is timereversible and, therefore, it seems that the difficulty, for the macroscopic irreversibility to be explained based on the reversible microscopic equation, still exists. However, according to Nelson [19], the basic equation for the quantum mechanics can be formulated by a new stochastic differential equation, i.e. Nelson’s equation,

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19]

L. Boltzmann, Wiener Berichte 63 (1872) 275. J. Loschmidt, Wiener Berichte 73 (1876) 128. E. Zermelo, Wiedemann Annalen 57 (1896) 485. J.L. Lebowitz, Phys. Today (September) (1993) 32. W.G. Hoover, Time Reversibility, Computer Simulation, and Chaos, World Scientific Publishing Co., 1999. J. Orban, A. Bellemans, Phys. Lett. 24A (1967) 620. J.M. Haile, Molecular Dynamics Simulation: Elementary Methods, John Wiley & Sons, Inc., 1997. D. Levesque, L. Verlet, J. Stat. Phys. 72 (1993) 519. W.G. Hoover, J. Chem. Phys. 109 (1998) 4164. N.S. Krylov, Works on the Foundations of Statistical Physics, Princeton University Press, 1979. H.A. Posch, W.G. Hoover, Phys. Rev. A 38 (1988) 473. W.G. Hoover, Computational Statistical Mechanics, Elsevier Science, 1991. W.G. Hoover, H.A. Posch, Chaos 8 (1998) 366. W.G. Hoover, Phys. D 112 (1998) 225. Ch. Dellago, H.A. Posch, Phys. A 230 (1996) 364. O. Kum, W.G. Hoover, J. Stat. Phys. 76 (1994) 1075. W. Nadler, H.H. Diebner, O.E. R¨ossler, Z. Naturforsch. 52a (1997) 585. G. Benettin, L. Galgani, J.M. Strelcyn, Phys. Rev. A 14 (1976) 2338. E. Nelson, Phys. Rev. 150 (1966) 1079.