Simulation of the Liquid State

BY D.M. HEYES

1 Introduction

The molecular simulation of liquids is a now vast field of research, and as with many others in recent years it is becoming increasingly difficult to keep abreast of all of the significant developments that are taking place. The ready availability of fast computers has meant that there are many more researchers working in this ever expanding field of applications, producing ever larger amounts of work to assimilate! This poses something of a problem when it comes to writing a review, especially one with the rather ambitious title of 'Computer Simulation of Liquids'. I am not going to attempt to cover all the branches of this field. Rather, in my review of the developments between 1999 and 2001 I have restricted my discussion to a few areas that have interested me. The choice is inevitably somewhat subjective, but hopefully by adopting this approach I will have a better chance of producing a useful document, rather than a gallop through many topics with only the briefest of discussion about each, which I am sure would be of little use to the scientific community. There are a number of molecular simulation books that describe the standard techniques, and these are recommended as background material for the present article, e.g. refs. 1-7.

I am therefore not going to discuss the 'nuts and bolts' of molecular simulation, except to mention an often overlooked fact, which is the reason for much of the success of these approaches. Most simulations are carried out still typically for less than a thousand molecules, and if it was not for the use of periodic boundary conditions (PBC) it would not be possible to simulate bulk systems with this number of molecules. These systems would have such a high surface to volume ratio that the results would be dominated by surface effects. The PBC procedure is illustrated in Figure 1, which shows a two-dimensional square cell in which the molecules are surrounded by image cells. A molecule near the cell boundary interacts with the 'real' molecules in the central cell and with image cell molecules. Molecules leaving the cell re-enter through the opposite face with the same velocity.

2 Simple Liquids

The study of simple liquids can be said to be the beginning of Molecular Dynamics and Monte Carlo in the 1950s and 60s. Although the scope of molecular simulation, as a field or discipline, has widened dramatically since then, there is nevertheless a continual interest in simple liquids. In fact, this is partly due to the fact that the so-called 'simple' liquids are far from simple! One of the motivations for the continual interest in the simple liquids is that, because of the basic nature of the interparticle interactions, an improved understanding of these systems should lead to better theoretical models, which can be extended to more complex molecular liquids. Also, the rapid growth of interest in colloids and polymers (so-called 'complex' liquids) in recent years has provided new areas where the theories of simple liquids can be applied, especially those associated with local structure and thermodynamics. In the latter case, phase equilibria and the location of phase boundaries feature prominently. In this section, some of the recent advances in our understanding of simple liquids are covered.

2.1 Dynamics. - Of course, within the category of 'simple liquids' studied by statistical mechanics and molecular simulation, there are model liquids that are, strictly speaking, not found in nature. For example, the ubiquitous hard sphere fluid, where the pair potential has the form

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2.1)

is a case in point. The energy is infinite on contact of the spheres (at σ) and zero for larger separations. This is the energy of interaction which would approximate that of two macroscopically sized elastic spheres with high elastic modulus, say two snooker or billiard balls. In these cases the length-scale of the particle interactions is many orders of magnitude smaller than the particle diameter. One of the main features of hard spheres is the co-ordination number. It has been shown recently that for spheres at random close packing, the mean number of particle contacts is 4.8, which is somewhat lower than has often been assumed before (6, and even 12, have been used). Interestingly these authors also performed a computer 'simulation' in which they took a random test sphere, and placed immobile point contacts on its surface. They determined the mean number of points required on the surface of the sphere to eradicate the possibility of translation of the test particle, which was found to be 2D + 1, where D is the space dimension. Therefore, in 3D this 'co-ordination number' is 7, which is lower than the value of 4.8, indicating that in states where the contacts are 'correlated' (i.e. in a dense liquid or glassy state) translational diffusion can be removed by fewer contacts. The procedure for carrying out Metropolis Monte Carlo of hard spheres is particularly simple, as the Boltzmann factor does not require specific evaluation – just overlap detection. Jater proposed an improved Metropolis Monte Carlo algorithm to simulate hard core systems, in which they replaced the usual sequence of single particle trial displacements by a collective trial 'move' of a chain of particles.

The hard-sphere system is widely used in statistical mechanics as a reference state in theories of liquids and solids. Its uses have traditionally been quite broad, extending from equations of state, the structure of molecular liquids and dynamical properties. As mentioned already, it has also found a new lease of life as a model reference system for some colloids and granular materials. Simple molecules (e.g. water) interact with an appreciable attractive tail extending beyond the hard core, and which usually has more or less the same range as the core. It is not possible to find a simple molecule that does not have an appreciable attractive or van der Waals region as well as a hard repulsive core. In contrast, on the micron and larger scale, for these systems, the hard-sphere can be an even more realistic representation of the effective pair potential, which can be steeply repulsive and have a negligible attractive component. It must be borne in mind, however, that the hard-sphere particle potential is fundamentally unrealistic in that its pair potential is discontinuous and non-differentiable unlike those for all real systems. Considerable care is therefore required in extrapolating from any steeply repulsive potential, made progessively steeper, to the hard-sphere potential. This is because many quantities diverge either to zero or infinity according to the order in which the limit is made, of potential steepness and number of particles etc., any of which factors may be significant. Non-physical results such as purely exponential and delta function time correlation functions may be generated. These issues have been considered recently by Powles, Rickayzen and Heyes in several publications. These authors have considered, amongst other forms, the generalised soft-sphere potential

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2.2)

where σ and ε represent the particle diameter and energy scale, respectively. With n = 12, this potential is no other than the repulsive part of the Lennard-Jones potential. If n is taken much larger than 12, above 100 say, then the fluid behaves in many respects more like the hard-sphere fluid. In the large n limit they referred to this repulsive potential as the 'steeply repulsive potential', or SRP for short. As n [right arrow] ∞ the properties approach those of the hard sphere. In the publications, the dynamical relaxation and transport coefficients as a function of n were explored, and in particular the behaviour of these systems as they become more like hard spheres. In addition to being of intrinsic interest, these fluids furnish new perspectives on the behaviour of the hard-sphere system itself. A useful link between the SRP results and those of a hard-sphere system can be established by ascribing an effective hard-sphere diameter, σ HS, which can then be used to specify equivalent hard-sphere packing fractions in the hard-sphere analytic formulae (e.g. Carnahan-Starling equation of state). The Barker and Henderson formula for σ HS , which is based on a free energy optimisation, was the prescription used in these studies,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2.3)

where ITLβITL= 1/kBT as usual. They showed that at short times the time correlation functions relevant to the collective-property transport coefficients (i.e. shear and bulk viscosities and thermal conductivity) all obey the same simple scaling behaviour. This is in terms of a renormalised 'time' constructed by multiplying time by the parameter characterising the steepness of the potential, n. The time t is replaced by nt.

An example of this scaling is shown in Figure 2. The reason for this can be shown to be rigorous in the hard-sphere limit by formal expansion of the time correlation function about t= 0. For a time correlation function, C(t) (normalised so that C(0) = 1) we find

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2.4

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and m is the mass of the SRP particle, which is a non-dimensional time that includes the temperature and the potential stiffness parameter, n, as already mentioned. Note that the correlation function is independent of density, which is not surprising as in the large n limit the short time decay of C(t) is dominated by binary collisions. It is tempting to extrapolate this series in the following closure,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2.5)

which fits the C(t) data for the heat flux, shear stress and pressure correlation functions down to at least C(t) = 0.5 (especially for the pressure case) quite well. Using the Green-Kubo formulae (relating these time correlation functions to the transport coefficients) one can write the transport coefficient as the product of a modulus, M∞, and a relaxation time, τ which is the integral of C(t) from t= 0 to ∞. M∞ can be calculated as a static average, and in the hard-sphere limit has a simple analytic form involving the equation of state of the equivalent hard-sphere fluid, and n. M ∞ diverges as n and τ goes as n-1 as n-1 [right arrow] ∞. Although the hard-sphere fluid has been very successful as a reference fluid, for example, in developing analytical equations of state it has serious deficiencies in accounting for the dynamical relaxation processes in real systems. The hard-sphere fluid, in fact, is not a good reference fluid for the short time ('β') viscoelastic relaxation aspects of rheology.

There are potential technical problems in simulating steeply repulsive potential fluids by MD, and in response a new molecular dynamics algorithm capable of integrating the equations of motion of impulsive-continuous potential systems (e.g. such as a hard sphere core combined with an r-n tail) has been derived using operator splitting techniques by Houndonougbo et al.

Anento et al. used MD to investigate the dependence of the dynamical properties of simple liquids on the softness of the potential core. They found that the longitudinal modes associated with density fluctuations propagate to higher wave numbers in liquids with softer potential cores. In contrast, the propagating transverse modes are weakly influenced by the softness of the potential core. Verdaguer and Padro' also examined velocity cross-correlations and momentum transfer between concentric shells around a typical particle. They concluded that velocity cross-correlations are not only the consequence of direct interactions, but also indirect interactions caused by the neighbours countering local density fluctuations. Wu et al. investigated by MD the instantaneous normal modes of simple Lennard-Jones-like and repulsive LJ liquids, focussing on the effects of the potential on localised resonant modes. Itagaki et al. used MD on a tapered LJ potential systems to calculate the self-part of van Hove's spatio-temporal distribution function,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2.6)

where Δri(t) is the displacement of particle i over a time t. In the intermediate time regime between 1 and 500 ps they found a 'pre-Gaussian' distribution of the general form,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2.7)

The function γ(t) decreases to zero as time elapses and passes over into the so-called macroscopic time regime.

Ohta et al. performed MD simulations of the wave dispersion relationships of Yukawa systems. The Yukawa potential has a screened Coulomb form,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2.8)

where κD is the inverse screening length, q is the charge and ε0 is the dielectric constant of free space. The Yukawa potential could represent, for example, a dust particle particle in a plasma or a colloid particle in an electrolyte, so this potential form is quite ubiquitous. The Yukawa system is characterised by two dimension-less parameters, [Kappa] = κDa, where a is the mean interparticle distance, a = (3/4πρ)1/3 where ρ is the average number density. There is also the inverse temperature, Γ = q2 /4[πε]0aT, where T is the temperature. The system is said to be strongly coupled if [GAMMA* = Γ exp(-kappa) is greater than unity. MD simulations of the wave dispersion relationships showed that there is a truncation of the transverse wave dispersion at long wavelength.

The search for an analytic expression for the velocity autocorrelation function, Z(t) = (1/3)(v(0).v(t)), where v(t) is the velocity of an arbitrary molecule at time t, has been of interest since the start of molecular simulation. The self-diffusion coefficient can be obtained from Z(t) using the Green-Kubo formula,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2.9)

A quasi-solid model of the liquid has proved popular. This idea is illustrated well by the treatment of Zwanzig, who derived the following expression for Z(t) in a simple liquid based on the 'hopping' of molecules between potential energy minima in the energy landscape,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2.10)

where p(ω) is the density of normal modes and the exponential pre-factor represents the hopping of the molecules between neighbouring sites with an associated characteristic relaxation timescale, τzw. The normal mode distribution is most appropriately determined by quenching the system to the nearest local energy minimum, and expanding the potential energy to second order to determine the frequency distribution for that configuration. Chisolm et al. used this approach more directly, from the set of normal modes and frequencies, ω[lamba]. Their final formula for Z(t) was

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2.11)

where the ω[lamba] are arbitrarily chosen vectors. This formula gave excellent agreement with the Z(t) determined by MD for liquid sodium. Molecular Dynamics has been used to calculate the self-diffusion coefficient and shear viscosity of liquid transition and noble metals using effective pair potentials obtained from the embedded atom model, with reasonable success.