Introduction

Practically all the simulations described so far have involved systems that are either in equilibrium or in some time-independent stationary state; while individual results are subject to fluctuation, it is the well-defined averages over sufficiently long time intervals that are of interest. In this chapter we extend the MD approach to a class of problem in which the behavior is not only time dependent, but the properties themselves are also spatially dependent in ways that are not always predictable. The analysis of the behavior of such systems cannot be carried out following the methods described earlier, which were generally based on the evaluation of system-wide averages or correlations, and one is therefore compelled to resort to graphical methods. Here we focus on MD applications in fluid dynamics, a subject in which atomic matter is conventionally replaced by a continuous medium for practical purposes; recovering the atomic basis is part of trying to understand more complex fluid behavior of the type studied in rheology. For more on the microscopic approach to hydrodynamics see [mar92].

Open systems

Most current MD applications deal with closed systems; this implies either total isolation from the outside world, or coupling to the environment in a way described by one of the ensembles of statistical mechanics. The coupling can occur, for example, with the aid of a thermostat (§6.3), in which case the equilibrium properties are those of the canonical ensemble.

Introduction

Some internal degrees of freedom are important to molecular motion, while others can be regarded as frozen. Classical mechanics allows geometrical relations between coordinates to be included as holonomic constraints. We have already encountered constraints in connection with non-Newtonian modifications of the dynamical equations (Chapter 6); here the constraints occur in a Newtonian context, so that there is little doubt as to the physical nature of the trajectories.

In this chapter we focus on a class of model where constraints play an important role, namely, the polymer models used for studying alkane chains and more complex molecules, in which a combination of geometrical constraints and internal motion is required. The treatment of constraints is not the only new feature of such models; the interactions responsible for bond bending and torsion are essentially three- and four-body potentials, and some rather intricate vector algebra is required to determine the forces. The particular alkane model described here incorporates one further simplification, namely, the use of the often encountered ‘united atom’ approximation – the hydrogen atoms attached to each carbon atom in the backbone are absorbed into the backbone atoms and are thereby eliminated from the problem.

Geometric constraints

Role of constraints

The notion of a constraint acting at the molecular level is merely an attempt at simplification; the justification for assuming that certain bond lengths and angles are constant is that, at the prevailing temperature, there is insufficient energy to excite the associated vibrational degrees of freedom (or modes) out of their quantum ground states.

Historical background

The origins of molecular dynamics – MD – are rooted in the atomism of antiquity. The ingredients, while of more recent vintage, are not exactly new. The theoretical underpinnings amount to little more than Newton's laws of motion. The significance of the solution to the many-body problem was appreciated by Laplace [del51]: ‘Given for one instant an intelligence which could comprehend all the forces by which nature is animated and the respective situation of the beings who compose it – an intelligence sufficiently vast to submit these data to analysis – it would embrace in the same formula the movements of the greatest bodies of the universe and those of the lightest atom; for it, nothing would be uncertain and the future, as the past, would be present to its eyes’. And the concept of the computer, without which there would be no MD, dates back at least as far as Babbage, even though the more spectacular hardware developments continue to this day. Thus MD is a methodology whose appearance was a foregone conclusion, and indeed not many years passed after digital computers first appeared before the first cautious MD steps were taken [ald57, gib60, rah64].

The N-body problem originated in the dynamics of the solar system, and the general problem turns out to be insoluble for three or more bodies. Once the atomic nature of matter became firmly established, quantum mechanics took charge of the microscopic world, and the situation became even more complicated because even the constituent particles seemed endowed with a rather ill-defined existence.

Introduction

In the study of equilibrium behavior, MD is used to probe systems that, at least in principle, are amenable to treatment by statistical mechanics. The fact that statistical mechanics is generally unable to make much headway without resorting to simplification and approximation is merely a practical matter; the concepts and general relationships are extremely important even in the absence of closed-form solutions. When one departs from equilibrium, very little theoretical guidance is available and it is here that MD really begins to fill the role of an experimental tool.

There are many nonequilibrium phenomena worthy of study, but MD applications have so far tended to concentrate on relatively simple systems, and the case studies in this chapter will focus on the simplest of problems. To be more specific, we will demonstrate two very different approaches to questions related to fluid transport. The first approach uses genuine Newtonian dynamics applied to spatially inhomogeneous systems, in which the boundaries play an essential role: simulations of fluids partly constrained by hard walls will be used to determine both shear viscosity and thermal conductivity. The second approach is based on a combination of modified equations of motion and fully homogeneous systems: the same transport coefficients will be measured, but since there are no explicit boundaries the dynamics must be altered in very specific ways to compensate for their absence.

Homogeneous and inhomogeneous systems

As computational tools, both homogeneous and inhomogeneous nonequilibrium methods have their strengths and weaknesses.

Introduction

In earlier chapters, polymer chains were represented as series of atoms coupled by customized springs (Chapter 9), or atoms coupled by rigid links whose length and angle constraints are handled by computations that supplement the timestep integration (Chapter 10). It is also possible to formulate the problem so that the only internal coordinates of the molecule are those actually corresponding to the physical degrees of freedom. Though the formalism involved, which is based on techniques used in robot dynamics, is more complex than the previous methods, the elegance of the approach and the fact that it provides an effective solution to the problem cannot be denied.

Chain coordinates

Consider a linear polymer chain of monomers. While in principle, each monomer (assumed to be a rigid object) contributes six mechanical degrees of freedom – abbreviated DOFs – to the chain, we use the argument of §10.2 to justify freezing the DOFs associated with variations in bond length and bond angle. Thus, apart from the first monomer which has six DOFs, each additional monomer contributes just a single DOF to the chain. Each such DOF corresponds to torsional motion, or twist, around the appropriate bond axis and is represented by a dihedral angle.

If each torsional DOF is regarded as a mechanical joint with a single rotational DOF that is associated with the site at one end of the link, then the system corresponds to a standard problem in the field of robotic manipulators for which techniques are available that express the dynamical equations of motion in a particularly effective manner [jai91, rod92].

Introduction

The importance of understanding the processes governing the transport of granular materials [jae96] has long been recognized, particularly because of its industrial relevance. Methods analogous to MD modeling turn out to be appropriate for the study of granular matter, although the constituent particles are, of course, no longer the atoms and molecules of MD.

Mere inspection reveals the complexity of granular matter. The grains themselves are irregularly shaped, often covered with asperities, and are normally polydisperse. Grain collisions are highly inelastic and friction is important for forming heaps. The wear and tear of collisions can alter the shape of the grains to some extent; electrostatic forces, moisture, adhesion and the presence of air can all affect the behavior. Which of these, and other, characteristics must be incorporated into the model to reproduce the key features of the behavior can only be established empirically.

The goal of this brief departure from simulation at the molecular scale is to demonstrate the wider applicability of the approach, but not to provide a survey of either granular dynamics simulation techniques or applications; reviews of the subject include [bar94, her95]. The discussion of this chapter deals with methods based on soft-particle MD and, while there are many fascinating granular systems to choose from, the examples here deal with vibrating layers [biz98, rap98], mainly because of the visual impact of the results. The methods are readily extended to other kinds of problem.

Introduction

The rigid molecule approach described in Chapter 8 is limited in its applicability, because it is really only appropriate for small, compact molecules. Here we consider the opposite extreme, namely, completely flexible molecules of a type used in certain kinds of polymer studies. No new principles are involved, since the intramolecular forces that maintain structural integrity by holding the molecule together, as well as providing any other necessary internal interactions, are treated in the same way as intermolecular forces. Later, in Chapters 10 and 11, we will consider more complex models, in which molecules exhibit a certain amount of flexibility but are also subject to various structural constraints that restrict the internal motion. The first case study in this chapter deals with the configurational properties of a single chain molecule in solution. The second deals with a model of a surfactant in solution, in which very short, interacting chain molecules are just one of the components of a three-component fluid; this very simple system is capable of producing coherent structures on length scales greatly exceeding the molecular size, as the results will demonstrate.

Description of molecule

Polymer chains

Owing to the central role played by polymers in a variety of fields, biochemistry and materials science are just two examples, model polymer systems have been the subject of extensive study, both by MD and by other methods such as Monte Carlo [bin95].

Introduction

The equations of motion used in MD are based on Newtonian mechanics; in this way MD mimics nature. If one adopts the purely mechanical point of view there is little more to be said, but if a broader perspective is permitted and MD is regarded as a tool for generating equilibrium states satisfying certain specified requirements, then it is possible to modify the dynamics and address a broader range of problems. But at the outset it must be emphasized that no physical meaning is attributed to the actual dynamics, and the approach is merely one of computational convenience for generating particular equilibrium thermodynamic states, although – and this is not an attempt to extract any such meaning – the deviations of the motion from the truly Newtonian may in fact be extremely small.

Conventional MD differs from most experimental studies in that it is the energy and volume that are fixed, rather than temperature and pressure. In statistical mechanical terms, MD produces microcanonical (NVE) ensemble averages, whereas constant-temperature experiments correspond to the canonical (NVT) ensemble; if constant pressure is imposed as well, as is generally the case in the laboratory, it is the isothermal–isobaric (NPT) ensemble that is the relevant one. While the choice of ensemble is usually one of convenience at the macroscopic level since (away from the critical point) thermal fluctuations are small, for the microscopic systems studied by MD the fluctuations of nonregulated quantities can be sufficiently large to make precise measurement difficult.

Molecular dynamics simulation provides the methodology for detailed microscopic modeling on the molecular scale. After all, the nature of matter is to be found in the structure and motion of its constituent building blocks, and the dynamics is contained in the solution to the N-body problem. Given that the classical N-body problem lacks a general analytical solution, the only path open is the numerical one. Scientists engaged in studying matter at this level require computational tools to allow them to follow the movement of individual molecules and it is this need that the molecular dynamics approach aims to fulfill.

The all-important question that arises repeatedly in numerous contexts is the relation between the bulk properties of matter – be it in the liquid, solid, or gaseous state – and the underlying interactions among the constituent atoms or molecules. Rather than attempting to deduce microscopic behavior directly from experiment, the molecular dynamics method – MD for short – follows the constructive approach in that it tries to reproduce the behavior using model systems. The continually increasing power of computers makes it possible to pose questions of greater complexity, with a realistic expectation of obtaining meaningful answers; the inescapable conclusion is that MD will – if it hasn't already – become an indispensable part of the theorist's toolbox. Applications of MD are to be found in physics, chemistry, biochemistry, materials science, and in branches of engineering.

This is a recipe book. More precisely, it is a combination of an introduction to MD for the beginner, and a cookbook and reference manual for the more experienced practitioner.

Introduction

The previous chapters concentrated on translating physical problems into practical simulations. Computational efficiency, beyond the use of cells and neighbor lists (as well as hierarchical subdivision when appropriate), received little attention. For ‘conventional’ computers, there is not a great deal more that the average user can do in this respect, assuming that a reasonably effective programming style has been adopted. This attitude is no longer adequate when modern, high performance, multiprocessor machines are to serve as the platforms for large-scale simulation.

In this chapter we focus on ways of adapting the basic MD approach to take advantage of advanced computer architectures; since enhanced performance comes not only from a faster processor clock cycle, but also from a number of fundamental changes in the way computers process data, this is a subject that cannot be ignored. The subject is also a relatively complex one and, at best, only peripheral to the goals of the practicing simulator. We will therefore not delve too deeply into the issues involved, but will merely focus on three examples, all of which can be valuable for large-scale MD simulation; the first employs message-passing parallelism, the second involves parallelism achieved by the use of computational threads and shared memory, and the third demonstrates how to rearrange data to achieve effective vector processing.

The quest for performance

It comes as no surprise to learn that spreading a computational task over several processors is a way to complete the job sooner.

Introduction

In this chapter we focus on a number of techniques used in MD simulation, primarily the methods for computing the interactions and integrating the equations of motion. The goal is to generate the atomic trajectories; subsequent chapters will deal with the all-important question of analyzing this raw ‘experimental’ data. We continue to work with the simplest atomic systems, in other words, with monatomic fluids based on the LJ potential, not only because we want to introduce the methodology gradually, but also because a lot of the actual qualitative (and even quantitative) behavior of many-body systems is already present in this simplest of models. Models of this kind are widely used in MD studies of basic many-body behavior, examples of which will be encountered in later chapters.

Equations of motion

While Newton's second law suffices for the dynamics of the simple atomic fluid discussed in this chapter, later chapters will require more complex forms of the equations of motion. The Lagrangian formulation of classical mechanics provides a general basis for dealing with these more advanced problems, and we begin with a brief summary of the relevant results. There are, of course, other ways of approaching the subject, and we will also make passing reference to Hamilton's equations. A full treatment of the subject can be found in textbooks on classical mechanics, for example [go180].

Introduction

Practically all the case studies in this book involve systems whose interactions are expressed in terms of continuous potentials. As a consequence, the dynamical equations can be solved numerically with constant-timestep integration methods. If one is prepared to dispense with this continuity another route is available that offers several advantages, although it has its weak points as well. The alternative method is based on step potentials; hard spheres are the simplest example, but the method can be extended to include potentials that have the shape of square wells or barriers, and even flexible ‘molecules’ can be built. Quantitative comparisons with specific real substances are obviously not the goal here, although comparisons with simple analytical models are possible. In fact, the earliest MD simulations [ald58, ald59, ald62] were of this kind, motivated by a desire to test basic theory.

A limitation of the methods used for continuous potentials, all of which involve a constant timestep Δt, is that they require the changes in interactions over each timestep to be small, otherwise uncontrolled numerical errors can suddenly appear. While this does not usually affect equilibrium studies, because Δt can be made sufficiently (but not too) small that for a particular simulation (namely, a given potential function, temperature and density) the results are predictably stable, systems that are inhomogeneous because of, for example, a large imposed temperature gradient, may prove problematic unless Δt is made unacceptably small.