Short history of liquid simulations

Milestones

- Metropolis Monte Carlo (Los Alamos, 1953) hard sphere potential

- Mote Carlo (MC) MC applied to Lennard-Jones (L-J) potential (Wood and Parker, 1957)

- Molecular Dynamics (Alder and Wainwright, 1957) hard sphere

- MD of L-J potential (Rahman, 1964, Verlet 1967)

- Diatomic molecular liquid (Harp and Berne, 1968)

Motivation

Only a few problems in statistical mechanics are exactly soluble. Computers simulations provide essentially exact results for problems that are not soluble analytically.

- Baxter, R. J. Exactly solved problems in statistical mechanics Academic Press

The force field

Molecular dynamics solves Newton’s equations of motion on a potential energy surface. We define the equations of motion and a potential energy surface for this classical calculation. In the calculation only the positions and momenta (velocities) of the nuclei will appear since it is assumed that we have applied the Born-Oppenheimer approximation. The significance of this is that the rapid motion of the electrons is averaged out. The nuclear positions are given by q = (q₁, q₂, q₃, …q_N) and p = (p₁, p₂, p₃, …p_N). The hamiltonian H is the sum of the kinetic and potential energies for the particles in the system. H(q,p) = K(p) + V(q). The q are usually the Cartesian coordinates of the individual nuclei.

where m_i is the molecular mass and a refers to the x, y, or z coordinate.

The potential may be divided into terms that depend on the coordinates of individual atoms, pairs, triples, etc.

Where the index j>i is used to ensure that we do not count any pairwise interactions twice. Pair potentials can be determined experimentally. However, the potential between an isolated pair of atoms or molecules is not the potential that one uses in a simulation. In fact the potential that is used in a simulation is an effective potential that accounts in part for the higher order terms in the expansion of the potential energy (three-body interactions etc.). These are likely less than 10% the magnitude of the two-body interactions, and are therefore not negligibly small, but are most often ignored in practical potential energy calculations.

The forcefield is an empirical fit to the potential. The forcefield employs internal coordinates (stretches, bends, and torsions) to describe the bonding part of the potential. The non-bonding interactions are described using Lennard-Jones or other non-bonding potential energy functions and Coulomb’s law for charge interactions and to mimic dipolar interactions.

Non-bonding or van der Waal’s terms

The Lennard-Jones potential is a commonly used effective potential in computer simulations.

The potential has a long-range attractive tail –1/r⁶, and negative well depth e, and a steeply rising repulsive wall at r = s. There are other potential functions such as the hard sphere, square well, etc.

For ions the Coulomb potential must by included.

This does not include polarization. A charge on a molecule can induce a dipole an on neighboring atom or molecule. The charges used in parameter sets however are not required to be only ionic charges. For example, the dipole or quadrupole of a molecule or a residue such as an amino acid or a nucleic acid can be represented by charges placed at the nuclear position (or elsewhere).

Bonding parameters

The energy expression for bonding parameters requires a defined equilibrium bond length, angle etc. and a force constant that describes the amount of energy required to distort from that equilibrium value. Covalent bonds are described using a harmonic oscillator approximation. In this approximation the potential energy surface is a quadratic function of the displacement.

Stretching	K( r – r0)^2
Bending	H( q – q0 )^2
Torsion	T[1 + cos(nf)]

In addition, there are cross terms that represent the interaction between any two of these quadratic functions. Even though internal coordinates are used for the forcefield terms, the molecular motions are described in Cartesian coordinates.

Atom types

Different bonding geometries may require different types of atoms to be defined. For example, sp³ and sp² hybrized carbons have different equilibrium bond angles and bond lengths. In addition, there are a number of different environments for each type of carbon, sp³ bonded to hydrogens, sp³ bonded to heavy atoms, sp³ in a five membered ring etc. The philosophy behind the force field is that it should be universal for a given atom type. The weakness of forcefields is generally that they have a large number of atom types to try to account for the large number of different environments an atom can be found in.

Molecular systems

The rigid molecule approach has been used for diatomics and triatomics. The internuclear distance is held fixed for atoms within the molecule. Interactions are calculated using hard sphere or L-J between each pair of atoms. However, even for alkanes the low frequency torsional motions of the molecule are important and therefore a potential energy function should be provided. In modern simulations using AMBER or DISCOVER the SHAKE routine can be used to fix the bond lengths of certain atoms. For example, hydrogen bond lengths with other atoms are often fixed to reduce the computational requirements. This is particularly useful because the frequency of C-H, N-H, and O-H stretches are significantly higher than other frequencies of vibration.

Instead of placing charges at the nuclei to represent the Coulombic interactions, several schemes have been developed to account for the higher order terms in the multipole expansion of a molecule. The charge in MD simulations is simply a parameter and there is no requirement that the charge be located at the position of a nucleus. The charges can have any position with respect to the axes of a molecular entity defined in the MD program.

Constructing and intermolecular potential

First we select reasonable parameters. For example, for the L-J potential, we use values of e (well depth) and s (van der Waal’s radius), available from experiment or previous computer simulations. These values represent a first guess that is refined during the course of the simulation. The values for an atom such as carbon can be quite different in different molecules. The cross terms are calculated using the Lorentz-Bertholet mixing rules. For example, in CS₂, the cross terms are given by

s_CS = ½[ s_CC + s_SS] and e_CS = [e_CCe_SS]^1/2

The first guess can be used to calculate a number of properties in the gas, liquid, and solid phase. A&T discuss calculation of the second virial coefficient on page 22. This is a gas phase approach to optimizing the parameters. Of course, there are many methods for further optimizing these parameters using a training set of crystal structures and data from known liquids and solids and testing the results of the force field against these known data.

For example, a series of computer simulations of the liquid state could be carried out. The densities and temperatures can be chosen to be near the orthobaric curve (i.e. the vapor-liquid coexistence line). The internal energy and pressure can be compared to experiment. The pressures are known experimentally (they are also given by the Clausius-Clapeyron equation). The internal energy can be determined from the latent heat of vaporization.

The solid state offers a sensitive test of any potential model. The energy minimized structure using the given force field can be compared with the experimental structure. In addition, if lattice dynamics have been determined experimentally these can also be compared (phonons, librational modes, and dispersion curves of the solid). In a solid or extended macromolecule the energy minimized structure can be thought of as being at zero Kelvin. Although this structure is static, it should be born in mind that there will be zero point motion in the solid.

Molecular systems

Isolated macromolecules

Frequently, studies of proteins, nucleic acids, or polymers use a gas phase approximation. If the structures are not solvated the isolated macromolecule dynamics can be studied. In such a simulation the relative motions of the macromolecule are of interest. For example, in a classic study of the conformational substates of myoglobin Elber and Karplus (Science, 1987, Vol. 235, page 318) used dynamics to demonstrate the motion of a-helices with respect to one another in myoglobin on the picosecond time scale. It is of interest to address the question of how solvation would change the structures and dynamics. To place a macromolecule in solution a solvent "box" must be created.

Periodic boundary conditions

Periodic boundary conditions refer to simulation of a molecular system in a periodic lattice of identical subunits. Imagine a cubic box replicated throughout space to form an infinite lattice. As a molecule moves in the original box, its periodic image moves in all of the boxes in the lattice in exactly the same way. As a molecule moves through one of the faces of the box, one of its images will enter through the opposite face.

We must compare the properties of the small infinitely periodic system and the macroscopic system. For a fluid of L-J atoms the box should be at least 6s without a particle being able to sense the symmetry of the periodic lattice. If a potential is long ranged (i.e. V(r) > 1/r³) there will be substantial interaction of a particle with its image in neighboring boxes unless there is a cutoff or the boxes are large enough. In practical terms the Coulombic terms in the force field are of the greatest concern since the potential decreases as V(r) µ 1/r. To see the magnitude of the Coulombic term we use the relation V(r ) » -332 (kcal/mole)Å/r. Thus, at a distance of 10 Å the interaction energy between two charges is 33 kcal/mole and even at 100 Å the interaction energy is 3.3 kcal/mole. This is a substantial interaction at such a large distance. Of course, dielectric screening will tend to reduce these interactions by a factor of 1/e where e is the low frequency dielectric constant. For this reason in such calculations it is necessary to truncate the potential (i.e. to apply a cutoff distance to the calculated potential) when periodic boundary conditions are used.

Potential truncation

Of course, if we set out to compute the pairwise interactions of a particle with all of the other particles in the system, this becomes an infinite number under periodic boundary conditions. The means to avoid this problem is to define a box around any given atom in a simulation and to set the potential function to zero at a certain radius around that atom. These two steps are known as the minmum image convention and potential truncation. The truncation of the potential can lead to artifacts due to the abrupt switching to zero. These artifacts can be minimized by use of a switching function.

A. Minimum image convention

Use a box surrounding any given atom that is identical in size and shape to the periodic box and includes the N – 1 other atoms in the simulation.

B. Potential cutoff

A second way to circumvent the fundamental problem of counting pairwise interactions is to introduce a cutoff. One arbitrarily fixes a distance at which the potential is set equal to zero.

V = v(r ) r < r_cutoff

V = 0 r > r_cutoff

The cutoff distance may be no greater than ½ L (L is the length of the box) for consistency with the minimum image convention.

C. Switching functions

A simple cut-off of the potential energy function leads to discontinuities in the energy and its derivatives. Molecular dynamics programs typically use a spline function to smoothly turn off the interactions over a range of distance to minimize this effect. The switching function is defined by two variables: the point where the function reaches zero and the range over which the function decreases from one to zero.