Computer Simulation Algorithms

The ensemble averages employed in Table 6.1, are usually taken to be canonical while the time dependence of the correlation functions is generated by field free Newtonian equations of motion. In §5.4, we proved that, except for bulk viscosity, thermostatted equations of motion can also be used to generate the equilibrium time correlation functions. For bulk viscosity the correlation function involves functions of the kinetic energy of the system. We cannot therefore use Gaussian isokinetic equations of motion (see equation (5.86) and (5.87)). This is because, for these equations, the kinetic energy is a constant of the motion.

To calculate thermal transport coefficients using computer simulation we have the same two options that were available to us in the mechanical case. We could use equilibrium molecular dynamics to calculate the appropriate equilibrium time correlation functions, or we could mimic experiment as closely as possible and calculate the transport coefficients from their defining constitutive relations. Perhaps surprisingly the first technique to be used was equilibrium molecular dynamics (Alder and Wainwright, 1956). Much later the more efficient nonequilibrium approach was pioneered by Hoover and Ashurst (1975). Although the realistic nonequilibrium approach proved more efficient than equilibrium simulations it was still far from ideal. This was because for thermal transport processes appropriate boundary conditions are needed to drive the system away from equilibrium - moving walls or walls maintained at different temperatures. These boundary conditions necessarily make the system inhomogeneous. In dense fluids particles pack against these walls, giving gives rise to significant number dependence and interpretative difficulties.

The most effective way to calculate thermal transport coefficients exploits the formal similarities between susceptibilities for thermal and mechanical transport coefficients. We invent a fictitious external field which interacts with the system in such a way as to precisely mimic the linear thermal transport process. The general procedure is outlined in Table 6.2. These methods are called 'synthetic' because the invented mechanical perturbation does not exist in nature. It is our invention and its purpose is to produce a precise mechanical analogue of a thermal transport process.

With regard to step 3 in Table 6.2, it is not absolutely necessary to invent equations of motion which satisfy AIΓ (see §5.3). One can generalise response theory so that AIΓ is not required. However it is simpler and more convenient to require AIΓ and thus far it has always proved possible to generate algorithms which satisfy AIΓ. Although AIΓ is satisfied, most sets of equations of motion used in synthetic NEMD are not derivable from a Hamiltonian. The preferred algorithms for thermal conductivity and shear viscosity are not derivable from Hamiltonians. In the case of thermal conductivity the Hamiltonian approach must be abandoned because of conflicts with the periodic boundary condition convention used in simulations. For shear viscosity the breakdown of the Hamiltonian approach occurs for deeper reasons.

Equations of motion generated by this procedure are not unique, and it is usually not possible a priori to predict which particular algorithm will be most efficient. It is important to realise that the algorithms generated by this procedure are only guaranteed to lead to the correct linear (limit F_e→0) transport coefficients. We have said nothing so far about generating the correct nonlinear response.

Many discussions of the relative advantages of NEMD and equilibrium molecular dynamics revolve around questions of efficiency. For large fields, NEMD is orders of magnitude more efficient than equilibrium molecular dynamics. On the other hand one can always make NEMD arbitrarily inefficient by choosing a sufficiently small field. At fields which are small enough for the response to be linear, there is no simple answer to the question of whether NEMD is more efficient than equilibrium MD. The number dependence of errors for the two methods are very different - compared to equilibrium MD, the relative accuracy of NEMD can be made arbitrarily great by increasing the system size.

These discussions of efficiency ignore two major advantages of NEMD over equilibrium molecular dynamics. Firstly, by simulating a nonequilibrium system one can visualise and study the microscopic physical mechanisms that are important to the transport processes (this is true both for synthetic and realistic NEMD). One can readily study the distortions of the local molecular structure of nonequilibrium systems. For molecular systems under shear, flow one can watch the shear induced processes of molecular alignment, rotation and conformational change (Edberg, Morriss and Evans, 1987). Obtaining this sort of information from equilibrium time correlation functions is possible but it is so difficult that no one has yet attempted the task. It is likely that no one ever will. Secondly, NEMD opens the door to studying the nonlinear response of systems far from equilibrium.

We will now give an extremely brief description of how one performs molecular dynamics simulations. We refer the reader to far more detailed treatments which can be found in the excellent monograph by Allen and Tildesley (1987) and in the review of NEMD by the present authors (Evans and Morriss, 1984a).

Consider the potential energy, Φ, of a system of N interacting particles. The potential energy can always be expanded into a sum of pair , triplet, etc., interactions:

(6.3)

For the inert gas fluids it is known that the total potential energy can be reasonably accurately written as a sum of effective pair interactions with an effective pair interaction potential denoted Φ(r_i,r_j). The Lennard-Jones potential, Φ^LJ, is frequently used as an effective pair potential,

(6.4)

The potential energy of the two particles i,j is solely a function of their separation distance r_ij and is independent of the relative orientation of their separation vector r_ij. The Lennard-Jones potential is characterised by a well depth ε, which controls the energy of the interaction, and a distance σ, which is the distance at which the potential energy of the pair changes sign due to the cancellation of the Van der Waals attractive forces by the short ranged quantum repulsive forces. If ε∕k_B=119.8K and , the Lennard-Jones potential forms a surprisingly accurate representation of liquid argon (Hansen and Verlet, 1969). For proper scaling during simulations, all calculations are performed in reduced units where ε∕k_B=σ=m=1. This amounts to measuring all distances in units of σ, all temperatures in units of ε∕k_B and all masses in units of m. The Lennard-Jones potential is often truncated at a distance, r_c=2.5σ. Other potentials that are commonly used include the Weeks-Chandler-Andersen potential, usually written as WCA, which is the Lennard-Jones potential truncated at the position of minimum potential energy (2^1/6σ) and then shifted up so that the potential is zero at the cutoff.

(6.5)

The main advantage of this potential is its extremely short range of interaction. This permits simulations to be carried out much more quickly than is possible with the longer ranged Lennard-Jones potential. Another short ranged potential than is often used is the soft sphere potential which omits the r^-6 term from the Lennard-Jones potential. The soft sphere potential is often truncated at 1.5σ.

In molecular dynamics one simply solves the equations of motion for a system of (N≅100-100000) interacting particles. The force on particle i, due to particle j, F_ij, is evaluated from the equation,

(6.6)

The N interacting particles are placed in a cubic cell which is surrounded by an infinite array of identical cells - so-called periodic boundary conditions. To compute the force on a given particle in the primitive cell one locates the closest (or minimum) image positions of the other N-1 particles. The minimum image of particle i, may be within the primitive cell, or in one of the surrounding image cells (see Figure 6.1). One then finds all the minimum images particles for i, that lie within the potential cutoff distance r_c and uses (6.6) to compute the contributions to the force on i, F_i=∑F_ij.

Finally one solves Newton’s or Hamilton’s equations of motion for the system

(6.7)

If, during the course of the motion, particle i leaves the primitive cell it will be replaced under the periodic boundary condition convention by an image of itself, travelling with exactly the same momentum, one lattice vector distant. We prefer to use Hamilton’s form for the equations of motion because this form is much more convenient than the Newtonian form both for NEMD and for equilibrium molecular dynamics with velocity dependent forces (such as thermostats). We often solve these equations of motion using a 5^th order Gear predictor-corrector method. In studies of the transient response of systems to external fields we use the less efficient Runge-Kutta methods. Unlike the Gear algorithms, Runge-Kutta methods are self-starting, achieving full accuracy in the first timestep.

We will now give a summary of some of the synthetic NEMD algorithms that have been used to calculate Navier-Stokes transport coefficients.