Over the past several decades, molecular simulation has become increasingly important for chemists, physicists, bio-scientists, and engineers, and plays a role in applications such as rational drug design and the development of new types of materials. While many levels of detail can be incorporated if desired, in most cases work is performed with a simplified atomic model, consisting of a large number of mass points interacting in various types of forces, i.e. an N-body problem.

There are essentially two principal types of simulation methodology in common use. In Monte-Carlo (MC) methods, random steps are taken in order to achieve a rapid sampling of the most likely states of the molecule. In molecular dynamics (MD), the idea is to construct approximate trajectories for the N-body problem and to use these to gain an understanding of how the molecule evolves in time, for example in response to a stimulus, during a transition between states, or as a means for calculating averages. It should be stressed that only MD and not MC methods allow the theoretical possibility of obtaining time-dependent quantities from simulation, while both schemes can in principle be used for the same statistical-mechanical calculations. Increasingly, one finds that MD and MC schemes are combined in various ways to seek improved efficiency. In this chapter we will focus only on (pure) MD methods, and in particular on the geometric integration issues associated to computing MD trajectories. For a more complete perspective on molecular simulation, the reader is referred to a text on the subject such as that of Schlick, Allen and Tildesley, Rappaport, or Frenkel and Smit.

Hamiltonian systems often exhibit dynamical phenomena covering a vast range of different time scales. In this chapter, we will discuss systems with two well separated time scales. More specifically, we consider systems for which the fast motion is essentially oscillatory. Such systems can arise from very different applications such as celestial or molecular dynamics and they might manifest themselves in very different types of Hamiltonian equations. Hence, the discussion in this chapter is necessarily limited to special cases. However, the basic principles and ideas have a much wider range of applicability.

A standard integrator, whether symplectic or not, will, in general, have to use a stepsize that resolves the oscillations in the fast system and, hence, one might be forced to use very small timesteps in comparison to the slow dynamics which is of primary interest. However, in special cases, one might be able to individually exactly solve the fast oscillatory and the slow system. Following the idea of splitting methods, this suggests to compose these two exact propagators and to apply a stepsize that is large with respect to the period of the fast oscillations. Such a method is then called a large timestep (LTS) method. Often the fast oscillations cannot be integrated analytically. A natural idea for the construction of an LTS method is then to assign different timesteps to different parts of the system. This approach is called multiple timestepping (MTS) and can often even be implemented such that the overall timestepping procedure still generates a symplectic map. We will explain the basic idea of symplectic LTS/MTS methods in Section 10.1.

In Chapter 2, we introduced the concept of a numerical integrator as a mapping which approximates the flow-map of a given system of differential equations. We have also seen a few instances of how such integrators behave, demonstrating concepts such as convergence and order of accuracy. We observed that the typical picture is a locally accurate approximation that gradually drifts further from the true trajectory (see Fig. 2.3, Fig. 2.5 and the left panel of Fig. 2.7); the rate of drift can be reduced by reducing the stepsize (and thereby also increasing the amount of computational work), but the qualitative picture does not change in any significant way.

What stands out as remarkable, therefore, is the behavior, illustrated in the right panel of Fig. 2.7, of the Euler-B method, which retains bounded trajectories when applied to the harmonic oscillator. In Chapter 2, we provided an explanation for this in the form of a linear stability analysis showing that certain methods, including Störmer–Verlet and Euler-B, have eigenvalues on the unit circle when applied to the harmonic oscillator (or any other oscillatory linear system), if the stepsize is below some threshold value. The Euler-B and Störmer-Verlet methods (among others) possess a strong asymptotic stability property for linear systems.

It is interesting to note that a related long-term stability property extends to nonlinear models. If we apply, for example, the Störmer–Verlet methods to the Lennard-Jones oscillator, we obtain the results illustrated in Fig. 4.1 (compare with Fig. 2.3 and Fig. 2.5).

In Chapter 4, we introduced several first- and second-order symplectic integration methods for Hamiltonian systems. In this chapter, we will discuss the construction of “higher-order” symplectic methods (with order greater than two), focusing in particular on those types of schemes that have been found to be most useful for practical computations. In traditional practice, higher-order integrators are employed for solving problems with relatively smooth solutions, such as gravitational simulations (solar system simulations, satellite trajectories). They are also traditionally used for many types of computations when very high accuracy (for example near the rounding error of the computer) is desired.

As we have seen in Chapter 2, the appropriateness of a given numerical method for a given computational task is a complicated issue. In some cases, the principles of geometric integration are in contradiction with the demand for high accuracy. If the purpose of simulation is to reconstruct, as exactly as possible, a particular trajectory segment, it may not matter what sort of qualitative features the integrator possesses: the efficiency of the integration method in terms of solution error per unit work is instead of paramount importance. Since the development of symplectic integrators adds a number of additional constraints on the design of the method, such schemes typically sacrifice something in efficiency compared with their nonsymplectic counterparts at similar accuracy, for example requiring an extra force evaluation or two at each timestep. Thus the problem of correctly determining the precise entry point and time instant that a space probe arrives at the Martian atmosphere is a task best handled by a standard integration method, for example, a high-order multistep integrator (for example, Diva) or explicit Runge–Kutta method (for example, RKSUITE).

In this chapter, we discuss formulation issues and symplectic integration methods for simulating the motion of a rigid body. Rigid bodies arise frequently in engineering, chemistry, and physics. For example, they occur in molecular simulation when the flexibility of small polyatomic units such as the water molecule, or CH4 is ignored. Cogwheels, space vehicles, and the planets are some other objects that are commonly modeled by rigid bodies.

Even in the absence of external applied forces, any rigid body more complicated than a uniform sphere will exhibit complicated motion, as defined by the moments of inertia of the body. A hint of the potential complexity of the motion is provided by the classic illustration using a hardbound book, which typically has three unequal moments of inertia I1 < I2 < I3 with I1 corresponding to an axis drawn along the binding, I2 to an axis across the cover, and I3 to an axis through the pages of the book (see Fig. 8.1). As the book is tossed up and spinning around each of the axes, the following dynamics are observed: around the first and third axes, the motion combines a stable periodic rotation with the rising and falling motion due to gravity, whereas the rotation with respect to the middle axis is much more complicated. (See, for example, for more explanation.) (It helps to place a rubber band around the book's cover to keep it closed while conducting experiments.)

Developing a method to simulate general rigid body motions, especially for long-term integration, proves an interesting and challenging task. The first issue we must confront is the selection of a set of coordinates that describe body orientation and spatial position.

About geometric integration

This book is about simulating dynamical systems, especially conservative systems such as arise in celestial mechanics and molecular models. We think of the integrator as the beating heart of any dynamical simulation, the scheme which replaces a differential equation in continuous time by a difference equation defining approximate snapshots of the solution at discrete timesteps. As computers grow in power, approximate solutions are computed over ever-longer time intervals, and the integrator may be iterated many millions or even billions of times; in such cases, the qualitative properties of the integrator itself can become critical to the success of a simulation. Geometric integrators are methods that exactly (i.e. up to rounding errors) conserve qualitative properties associated to the solutions of the dynamical system under study.

The increase in the use of simulation in applications has mirrored rising interest in the theory of dynamical systems. Many of the recent developments in mathematics have followed from the appreciation of the fundamentally chaotic nature of physical systems, a consequence of nonlinearities present in even the simplest useful models. In a chaotic system the individual trajectories are by definition inherently unpredictable in the exact sense: solutions depend sensitively on the initial data. In some ways, this observation has limited the scope and usefulness of results obtainable from mathematical theory. Most of the common techniques rely on local approximation and perturbation expansions, methods best suited for understanding problems which are “almost linear,” while the new mathematics that would be needed to answer even the most basic questions regarding chaotic systems is still in its infancy.

We have seen in the previous chapter that integrators preserving symplectic structure and/or first integrals can often be constructed in a straightforward way. In this chapter, we consider the properties of those methods and the implications for long-term simulations.

The traditional approach of numerical analysis generally assumes that the purpose of simulation is the faithful reproduction of a particular solution or trajectory, but individual trajectories typically are not of primary interest in most modern, scientific research; rather, the scientist typically treats the trajectory as a particular realization of a fundamentally stochastic evolution modelling in some way the myriad undetermined perturbations present in a “real-world” environment. It was the important discovery of Lorenz that differential equations can exhibit a chaotic solution behavior that includes an essentially stochastic or “random” component. The scientist views the model being analyzed as representative of a class of nearby models based on parameters which are typically only empirically (and approximately) determined. Furthermore, exact initial conditions are also typically not available. Some classical examples of such a scenario are molecular dynamics and numerical weather prediction.

It is now apparent that most modern large-scale simulations are conducted with timesteps and time intervals such that the numerical solution cannot be thought of as close to any particular model solution. The purpose of wedding the development of integrators to the standard axiomatic principle of timestepping – that one is attempting to approximate a particular trajectory – is thus called into question. Although high accuracy often is not needed in nonlinear dynamics computations, we must recognize certain important constraints imposed by the laws of nature.

In this chapter, we examine the feasibility of the algorithms in Part Two of this book using the implementations described in the last chapter. To make our experiments meaningful, we use real-world data from a variety of different sources, in a variety of different sizes. We time each algorithm to examine its running time behavior in practice. We also gather statistics on significant structural information contained in the data, such as number of conflicts or collisions in the persistence algorithm.

We begin by introducing the three-dimensional data for α-complex filtrations in Section 12.1. This is the data we use for timings and experiments on the persistent algorithm for ℤ2 coefficients in Section 12.2, topological simplification in Section 12.4, and the linking number algorithm in Section 12.6. We introduce alternate data for the persistence algorithm for fields in Section 12.3 as well as the Morse-Smale complex algorithm in Section 12.5. When appropriate, we also discuss additional implementation details not included in the last chapter.

Three-Dimensional Data

In Chapter 1, we motivated the study of topological spaces through a few diverse examples. It is appropriate, therefore, that the experimental data be from disparate sources. We use data that range in scale from nanometers to centimeters. The data will include proteins and inorganic molecules, resolved molecular structures, designed synthetic molecules, acquired samples from real world objects, and sampled mathematical functions.