By Hans H. Diebner, 10-März-2000
In 1967 Orban and Bellemans pointed in a brief paper  to numerical instabilities with respect to molecular dynamics simulations (MDS). Such instabilities question the reliability and significance of MDS that are carried out for entropic studies. Needless to say, if the so-called time-reversibility paradoxon is under investigation purely numerical dissipation must be avoided. In principle, this problem can be tackled nowadays by using an algorithm out of the class of so-called "symplectic algorithms".
For an understanding of the problem please recall the fact that the fundamental equation of motion of a conservative systems, the Newtonian equation of motion, is invariant under time reversal. For a typical N-particle problem this means that no primary direction of time exists. Thus, a change of the sign of all particles' momenta allows to exactly retrace the trajectory in the reverse time direction. This implicates that also an arbitrary function of the microstates must be retraceable under the exchange of momenta. For example, Orban and Bellemans computed the Boltzmann H-function in their simulation of a gas yielding the following results:
Fig. taken from .
The H-Funktion decreases in the time course and asymptotically reaches a minimum - its equilibrium value. The latter is reached after a total relaxation of the gas which is, in general, a very fast process. After 10 time units in a first and after 20 time units in a second run, respectively, the authors exchanged the sign of all momenta. The shape of the H-function then should show a symmetry with respect to the vertical line throught t=10 or t=20, respectively. In other words: the initial state of the gas should be exactly reproduced. Obviously, the later the reversal takes place (t=10 and t=20, respectively) the more the theoretically expected symmetry is violated. Additionally, a superposition of noise to each interation (relative error:(a)10-8, (b)10-5 and (c)10-2) leads to a even more destroyed symmetry. The noise mimics a surrounding which leads in effect to an open system.
The fact, that the evolution of the system back to its initial non-equilibrium state is supressed after a superposition of noise can be well understood in the light of a paper by James Hurley . A discussion of that phenomenon is far beyond the scope of that brief treatise. Instead, we prefer to focus on a phenomenon which is of fundamental importance for reliable MDS. If one uses one of the established standard algorithms for solving a differential equation implemented on a floating-point representation on a digital computer the superposition of random noise is unavoidable. Indeed, this was the real shock released by Orban and Belleman's paper. Rounding processes are of non-linear nature a manifest themselves as noise. A conservative Hamiltonian system becomes an open system, i.e., a superpositon of a numerical dissipation takes place. The analytical results correspond to that open system and not to the conservative one under investigation.
The problem at hand is closely related to problems that are encountered in chaos research which resulted in controversial debates on the reliability of results within that field of research that is extensively supported by numerical computations. There are still some critics who speak of "computer generated artifacts". However, it could be shown analytically as well as by means of "real" experiments that chaotic phenomena indeed exist. A rehabilitation of "computer experiments", however, could not be given until the publications by the numerics expert Peter Kloeden and his colleagues . They were able to show that there exists a class of algorithms that are capable of reproducing invariant structures of the phase space like chaotic attractors, for example, within an arbitrary small vicinity. An analogous derivation of so-called symplectic algorithms  can be performed by using the same recipe as above. Such algorithms are constructed so that certain invariant structures are conserved up to an arbitrary small deviation.
The decisive invariant entities in the MDS context are energy and momentum. An important contribution towards algorithmic consistency is the derivation of a symplectic algorithm from the principle of least action . The application of a variational method to the discrete space-time of a digital computer yields a "discrete Newtonian equation of motion". The result is an algorithm that exactly conserves momentum, prevents from a secular drift in energy and is exactly reversible. The seminal idea stems from Gillilan and Wilson , who applied, as a slight shortcoming, the action principle only to discretized time. Levesque and Verlet  and independently Diebner  recogized that furthermore the algorithm has to be implemented using integer arithmetik.
To conclude, Hamiltonian N-particle systems can now be consistently investigated by numerical means. The problem of non-retracability of Boltzmann's H-function after a reversal of momenta can be tackled which has been shown in [7,8]. The MDS and entropic investigation in part I of this brief treatise on this web-site has been done by using the symplectic exactly reversible algorithm and, thus, finds an additional justification. A detailed treatise of this topic and related themes can be found in .
Literature: J. Orban and A. Bellemans:
 Diese Web-Site: Molekular-Dynamik-Simulationen I
 Hans H. Diebner: