Computer simulation of simple and complex atomistic fluids by nonequilibrium molecular dynamics techniques

Computer Physics Communications 142 (2001) 14–21 www.elsevier.com/locate/cpc

Computer simulation of simple and complex atomistic fluids by nonequilibrium molecular dynamics techniques B.D. Todd Centre for Molecular Simulation and School of Information Technology, Swinburne University of Technology, P.O. Box 218, Hawthorn, Victoria 3122, Australia

Abstract Nonequilibrium molecular dynamics (NEMD) is introduced. Examples are presented to show how thermodynamic and transport properties of simple and complex fluids under a variety of flow fields are calculated. Its usefulness as a predictive tool is demonstrated by applications that throw into question standard constitutive equations used in computational fluid dynamics calculations. Its promise as a tool for the rational design of new polymer materials is briefly discussed.  2001 Elsevier Science B.V. All rights reserved. PACS: 44.10.+i; 44.30.+v; 47.50.+d; 61.20.Ja; 61.25.Hq; 66.20.+d Keywords: Nonequilibrium molecular dynamics; Computational fluid dynamics; Transport properties of fluids; Constitutive equations

1. Introduction High performance computing has impacted significantly on almost every field of scientific endeavor over the past decade. Of particular interest to the community of fluid dynamists (though yet relatively unknown to this same community) are the rapid developments in the field of nonequilibrium molecular dynamics (NEMD), which may be viewed as computational statistical mechanics for nonequilibrium fluids. In short, NEMD involves the application of classical mechanics and nonequilibrium statistical mechanics to determine the rheological properties of fluids. Due to rapid improvement of computational power over the last few years, the fluids now accessible are no longer simple atomic fluids (such as liquid argon). These days the simulation of liquid alkanes to determine their transport properties is progressing rapidly, as is the appliE-mail address: [email protected] (B.D. Todd).

cation of NEMD to polymer melts and solutions under shear and extensional flows. These methodologies have increasing significance for the lubrication, polymer processing and nanotechnology industries. This paper serves as a review of the basic concepts of NEMD, as well as some specific applications to demonstrate its power and usefulness to the fluid dynamics community. As a first step, it is shown how a simple fluid may be homogeneously sheared (in essence, a microscopic shearing cell of a bulk fluid flow) via the so-called SLLOD equations of motion. From these simulations one may extract all the relevant thermodynamic properties of the fluid, as well as its transport coefficients, such as the diffusion, viscosity and thermal conductivity coefficients. NEMD simulations are not only confined to idealized homogeneous fluids. An example of a fluid confined between microchannel atomistic walls undergoing planar Poiseuille flow will be presented, and ev-

0010-4655/01/$ – see front matter  2001 Elsevier Science B.V. All rights reserved. PII: S 0 0 1 0 - 4 6 5 5 ( 0 1 ) 0 0 3 0 4 - 6

B.D. Todd / Computer Physics Communications 142 (2001) 14–21

idence given that Fourier’s law of heat conduction is an inadequate description of heat transport. Its consequences for the formulation of an appropriate constitutive equation for heat transfer, and the Navier–Stokes heat equation will be discussed. Finally, some simulation results for polymer melts under shear and extensional flows will be presented, opening up the exciting possibility of rational materials design at the molecular level by computer simulation.

2. Nonequilibrium molecular dynamics and homogeneous flow NEMD involves the simulation of a classical N body system of atoms or molecules interacting via interatomic potential forces. However, in addition to interatomic forces, an external field drives the system away from thermodynamic equilibrium. Generally NEMD is used to study fluid systems, though in principle there is no restriction for using it to study either gases or solids under the influence of external forces. For this talk, only fluid systems are considered. The equations of motion for a system of N single component atoms may be formulated as follows: r˙ i =

pi , mi




 pi p˙ i = Fi + FE − α − u(r, t) , mi

(1)

where ri represent the laboratory positions of atom i and pi the laboratory momenta. Fi is the sum of all interatomic forces acting on atom i due to all other N − 1 atoms, FE is an external field that drives the system away from thermodynamic equilibrium (e.g., a gravitational field or a driving pressure head), α is a thermostat multiplier that may be used to keep the temperature of the system constant, and u(r, t) is the streaming velocity of the fluid. These equations of motion are numerically integrated at each timestep to generate the particle phase-space trajectories of the system. From these trajectories one may use the principles of nonequilibrium statistical mechanics to calculate all the thermodynamic properties of interest of the fluid (e.g., elements of the pressure tensor or heat flux vector), as well as its transport properties. Of particular usefulness is the determination of transport

15

coefficients, such as diffusion, viscosity and thermal conductivity coefficients. In the case of molecular fluids, the equations of motion need to be formulated such that they include additional bond constraint forces. Furthermore, the thermostatting mechanism needs to be correctly formulated. Thermostats act only on the peculiar, or thermal, component of an atom’s velocity, so the streaming component of the fluid’s flow must be subtracted out of the total velocity. This is seen in Eq. (1) for an atomic fluid. However, molecules have rotational as well as translational degrees of freedom, and a correct formulation of a thermostat requires an accurate knowledge of both the local instantaneous translational and angular streaming velocities. The latter is a computationally daunting task, even with current high performance computing resources [1,2]. If rotational degrees of freedom are ignored or incorrectly determined the result is an artificially induced non-zero antisymmetric stress. This, as well as incorrect assumptions of the translational streaming velocity, can also lead to enhanced alignment of the molecules with respect to the flow field and the onset of so-called string phases. The simplest and most feasible solution for the time being has been to re-cast Eq. (1) such that the thermostat acts only on the center of mass of each molecule [3–5]. In this case the angular velocity of the center of mass is zero, and one only needs to calculate peculiar velocities with respect to the translational streaming velocity. Such a formulation ensures that the total antisymmetric stress tensor is zero, and there are no enhanced ordering effects due to incorrect assumptions of the streaming angular velocity. A recent alternative formulation is to thermostat the so-called configurational temperature [6,7]. As this temperature is determined only from atomic coordinates, the momenta do not feature in its calculation. This obviates the need to explicitly calculate the local streaming velocity. A thorough comparison of both formalisms remains to be done. Of relevance to computational fluid dynamicists is how NEMD complements computational fluid dynamics (CFD). It does this primarily in two ways: first, NEMD allows a direct calculation of the transport coefficients of a fluid, assuming only knowledge of the interatomic forces acting between atoms. This means that the output of a typical NEMD simulation (i.e. the transport coefficients) may be used as input for

16

B.D. Todd / Computer Physics Communications 142 (2001) 14–21

a typical CFD calculation. The transport coefficient is thus calculated ab-initio, without recourse to semiempirical approximations or experimental data. Secondly, the ab-initio nature of NEMD simulations allows one to extract thermodynamic fluxes directly and to relate them to their conjugate thermodynamic forces. This in turn allows one to examine in greater detail the validity of a host of constitutive equations used for closing hydrodynamic equations. Such equations may be as fundamental as the standard Navier– Stokes equations for momentum and heat transport, or may invoke more sophisticated constitutive equations for modeling the viscoelasticity of complex polymer solutions or melts [8]. We will examine the former situation in Section 3, by looking at the example of planar Poiseuille flow, and how this throws into question the validity of the standard Navier–Stokes heat transfer equation. As an example of the direct calculation of a transport coefficient, consider a homogeneous fluid of atoms under planar shear flow in the x-direction, where the velocity gradient is constant in the ydirection. If we are interested in the shear viscosity, this may be extracted directly as η=

σxy + σyx  , 2γ˙

(2)

where η is the shear viscosity, σxy and σyx are the xy and yx elements of the stress tensor, and γ˙ is the strain rate imposed by the flow field, given as γ˙ = ∂ux /∂y. The angle brackets indicate that the stresses are time averaged, and we note that at steady-state, σxy  = σyx . The steady-state stresses are calculated directly from the appropriate statistical mechanical expressions [9]. Eq. (2) is of course none other than an expression of Newton’s law of viscosity. Thus, for a variety of thermodynamic state points, one may calculate the shear viscosity as functions of temperature, density and strain rate. This information could conceivably be fed back into a CFD calculation, which needs η(ρ, T , γ˙ ) as input. One may thus generate a three-dimensional grid of input viscosities for a CFD calculation. The output of such a CFD calculation may be, for example, the stress and velocity profiles of a shearing fluid where the density, temperature and strain rate may vary as a function of position. Eq. (1) as currently formulated poses some problems. In order to simulate boundary driven shear flow,

for example, one would need suitable boundaries, such as a solid wall. On microscopic length scales this leads to very strong inhomogeneities in the fluid density, which makes analysis of simulation data complex. An alternative may be to apply a spatially varying external field with suitable periodic boundary conditions to generate shear flow, but this also induces strong inhomogeneities making data analysis difficult [10]. Ideally we would like to simulate a small volume of bulk fluid, where boundary effects may be ignored and the fluid is genuinely homogeneous. Indeed, Eq. (1) may be reformulated such that these conditions are precisely met. The resulting equations are known as the SLLOD equations of motion [11], and for an atomic fluid are given as pi + ri .∇u, mi p˙ i = Fi − pi .∇u − αpi .

r˙ i =

(3)

Now, pi is already defined by Eq. (3a) to be peculiar with respect to the streaming velocity u. If the flow is planar shear flow, then (3) reduces to pi + iγ˙ yi , mi p˙ i = Fi − iγ˙ pyi − αpi . r˙ i =

(4)

These equations now have the salutary property of transforming a boundary driven flow (such as planar Couette flow) into one that is driven by a smooth mechanical force. When implemented with suitable periodic boundary conditions, such as those devised by Lees and Edwards [12], they generate homogeneous isothermal shear flow, making extraction of the shear viscosity a relatively simple task. For planar elongational flow, where the fluid expands at a rate of ε˙ in the x-direction, and contracts at the same rate in the y-direction (the fluid remains static in the z-direction), Eq. (3) reduces to pi + i˙εxi − j˙εyi , mi p˙ i = Fi − i˙εpxi + j˙εpyi − αpi . r˙ i =

(5)

For a system of molecules, Eq. (3) becomes piδ + ri .∇u, miδ miδ miδ C p˙ iδ = Finter pi .∇u − αM pi , iδ + Fiδ − Mi Mi

r˙ iδ =

(6)

B.D. Todd / Computer Physics Communications 142 (2001) 14–21

where riδ is the position of atomic site δ of molecule i, ri is the center of mass of molecule i, piδ is the peculiar momentum of atom δ of molecule i, Finter is iδ the total intermolecular force on atom δ of molecule i, FC iδ is the total bond constraint force acting on atom δ of molecule i, miδ is the mass of atom δ of molecule i, Mi is the center of mass of molecule i, pi is the center of mass momentum of molecule i, and αM is a molecular thermostat multiplier that acts on the center of mass momentum of molecule i. Details of these equations of motion, and the proper formulation of the molecular thermostat, are given elsewhere [1–5]. The use of the SLLOD equations of motion to investigate the rheology of molecular fluids is increasing rapidly. It has been particularly useful in studying the rheology of lubricants (e.g., see Ref. [13]), but is now on the verge of being routinely used to study the rheology of polymer melts and solutions (e.g., see Ref. [14]).

3. Inhomogeneous flow. Case study: planar Poiseuille flow Consider a 3-dimensional fluid of atoms confined by atomistically detailed walls, such as depicted in Fig. 1. In this geometry, an external field (e.g., gravity) acts on fluid atoms in the x-direction, generating a streaming velocity u(r, t) = iux (r, t). The walls are separated in the y-direction, and all thermodynamic quantities we are concerned with are functions of this coordinate. The z-axis is normal to the page. The interatomic potential used is the shifted and truncated

17

Lennard-Jones potential [15]. This system is described in detail in Ref. [16]. As the fluid is now inhomogeneous, the SLLOD equations of motion cannot be used. Instead, we recast Eq. (1) for both wall and fluid atoms. The wall equations are pi r˙ i = , mi (7) p˙ i = Fi − κ(ri − qi ) − αpi − jζLj . Here κ is a spring constant, and the second term in the momentum equation represents an harmonic tethering potential used to confine wall atoms to vibrate about their equilibrium lattice sites qi . α is the thermostat multiplier, and ζ is a constraint multiplier used to fix the center of mass position of each wall layer. This is necessary so that the walls do not expand as the fluid heats up due to viscous heat generation. As viscous heat is removed only through the walls via the wall thermostat, the equations of motion for the fluid atoms are purely Newtonian, i.e. pi , r˙ i = mi (8) p˙ i = Fi + iFE , where FE is the external driving field. When solved, the classical Navier–Stokes heat equation predicts a quartic temperature profile. However, this assumes that the thermal conductivity and shear viscosity are constant. Strictly speaking, this is not true, as both these transport coefficients are state-point dependent. If we ignore the first few monolayers of fluid atoms near the wall surfaces, then the characteristic wavelengths of the transport properties of the fluid

Fig. 1. Schematic diagram of the geometry of the Poiseuille flow simulation. The z-axis is into the page.

18

B.D. Todd / Computer Physics Communications 142 (2001) 14–21

are large compared to the range of the interatomic potential. Thus we may assume localization of momentum and energy transport, and write the governing constitutive equations as: σxy (y) = η(y)γ˙ (y), JQ (r) = jJQy (y) = −jλ(y)

∂T (y) . ∂y

(9)

The first constitutive equation is just Newton’s law of viscosity, while the second is Fourier’s law of heat conduction. Both equations now include the explicit position (hence state point) dependence of the transport coefficients. T (y) is the position dependent temperature of the fluid. The steady-state Navier– Stokes heat equation is thus ∂JQy (y) − σxy (y)γ˙ (y) ∂y   ∂T (y) ∂ λ(y) − η(y)γ˙ 2 (y) = 0 =− ∂y ∂y

(10)

which, when solved explicitly for T (y), yields T (y) = T0 + T4 y 4 + T8 y 8 + · · · .

(11)

The coefficients T0 , T4 , T8 , etc., are functions of the density, viscosity, thermal conductivity and driving field. Details of this derivation, and the assumptions made about the viscosity and thermal conductivity profiles, may be found in [16]. Even with the explicit y dependence of the transport coefficients included, the Navier–Stokes (N–S) heat equation still predicts a temperature profile that is quartic in its leading term. In Fig. 2 the results of our NEMD generated temperature profiles are plotted and compared with the theoretical prediction of Eq. (11) [16]. Clearly the prediction is poor in the center of the channel and does not correctly capture the profile on either side. Furthermore, the value of the thermal conductivity at the center of the channel (i.e. y = 0) is predicted to be λ(0) = 3.89 (reduced units). Independent NEMD simulations using the Evans heat flux algorithm to compute the thermal conductivity [17] give a known value of 6.89 ± 0.05 at the state point of (ρ(0), T (0)) = (0.836, 0.955). The classical N–S equation thus predicts the wrong temperature profile as well as the wrong value of λ(0)! One is tempted to dispute the simulation results, given the long and successful history of classical N–S

Fig. 2. Temperature profile for Poiseuille flow. The data points are obtained from simulation, whereas the solid curve is the predicted Navier–Stokes profile.

fluid dynamics. However, it should be noted that microscopic flows are not the common experience of classical macroscopic flows, where the traditional N–S solution appears undisputed. Microscopic flows involve length and timescales that often invalidate the central assumptions of macroscopic fluid dynamics, and so must be treated appropriately at the microscopic level. It is this level which holds the key to new nanotechnologies. In 1992, Baranyai, Evans and Daivis performed simulations of a fluid under the influence of a sinusoidal transverse force field [10]. By use of an appropriate thermostat, they nulled out all the temperature harmonics but still found a heat flux transverse to the flow velocity, in stark contradiction to Fourier’s law, which predicts zero heat flux for zero temperature gradient. They then postulated that the heat flux vector couples to an additional term proportional to the gradient of the strain rate tensor twice contracted into its transpose, i.e.   JQ = −λ∇T − ξ ∇ ∇u : (∇u)T , (12) where ξ is a phenomenological strain rate coupling coefficient. While they were able to explain their results with this assumption, the work was criticized for the synthetic thermostatting mechanism used to regulate the fluid temperature. In our Poiseuille flow simulations, we suffer from no such artefact. In fact, the fluid is not thermostatted at all. Viscous heat is entirely removed by the walls, and the fluid is described

B.D. Todd / Computer Physics Communications 142 (2001) 14–21

19

application to complex molecular fluids unattempted. The latter in particular remains a significant challenge and may hold the key to its significance in real world applications. 4. Polymer melts under planar shear and elongational flows

Fig. 3. As with Fig. 3, but now the curve represents the theoretical prediction with the inclusion of strain-rate coupling.

by Newton’s equations without any constraint forces, thus more closely representing nature. Substituting Eq. (12) into the energy continuity equation one obtains a modified N–S heat equation:     ∂ ∂T (y) ∂ ∂ γ˙ (y)2 λ(y) + ξ(y) ∂y ∂y ∂y ∂y + η(y)γ˙ 2(y) = 0

(13)

Rapid increases in computational power implies that the study of molecular rheology by simulation methods on the verge of practical realization. NEMD simulation of polymer melts and solutions, while limited at this time, will most likely become a fundamental tool in the study of polymer viscoelasticity. As was demonstrated by the previous example, NEMD can be used to test the validity of assumed constitutive equations. This is potentially a rich source of study in polymer rheology, where the number and type of constitutive equations are proliferating (see, for example, Ref. [21]). One example would be to use NEMD to check the predictions of the tube model of viscoelasticity, otherwise known as the Doi–Edwards model [8]. In this model, the constitutive equation for the stress tensor is given as: t

which when solved for the temperature yields T (y) = T0 + T2 y 2 + T4 y 4 + T6 y 6 + · · · .

(14)

The leading term is now quadratic in y. In Fig. 3 the simulation data is now compared with the profile predicted by Eq. (14). The agreement is excellent. Furthermore, the thermal conductivity at the center of the channel is predicted to be λ(0) = 6.93, which is within the range of errors of the exact known value of 6.89. Our NEMD simulations have thrown into question the validity of Fourier’s law of heat conduction. This has now been confirmed and validated by recent independent work [18–20]. These results particularly sound a warning that the use of Fourier’s law in the standard N–S heat equation is likely to lead to significant errors in the prediction of thermal properties of fluids. This may be extremely important in the prediction of microporous/microchannel flows, where the effect is far more pronounced than in macroscopic flows [16,20]. While the existence of this new transport coefficient appears confirmed, its state-point dependence is not yet characterized, and its

σαβ (t) = Ge −∞

dt






 ∂ψ(t − t ) (I A)  Qαβ E(t, t ) , ∂t (15)

where Ge is a constant, ψ is related to the time correlation function of the end-to-end vector of a prim(I A) itive chain, and Qαβ is the orientational tensor, itself a function of the deformation gradient E. The socalled independent alignment approximation is made in this formulation of the stress tensor. The parameters that appear in these terms, such as the diffusion coefficients, step length of the primitive chain and characteristic relaxation times, may be extracted directly from NEMD simulation, and then fed back into Eq. (15). The resulting predicted stresses (hence viscosities) can then be directly compared with those from the simulation. A number of variants of the original Doi– Edwards model exist, and NEMD may be used to critically assess their merits and deficiencies, as well as lead to more accurate theoretical formulations. This work is currently being undertaken at Swinburne. While simulations of planar shear flow are now routine, simulations of planar elongational flow have only

20

B.D. Todd / Computer Physics Communications 142 (2001) 14–21

Fig. 4. Shear and elongational viscosities of a system of 256 20-site molecules. Viscosities are plotted against the second scalar invariant of the strain rate tensor. Circles represent the shear viscosity for molecules undergoing planar Couette (shear) flow, squares represent elongational viscosities for planar elongational flow at constant volume, and diamonds represent the same for constant pressure simulations.

recently been possible by implementation of periodic boundary conditions (pbcs) devised by Kraynik and Reinelt [22]. These pbcs were first used in NEMD simulations in 1998, and the details may be found in Refs. [23–26]. We model our linear polymers by Lennard-Jones tangent chains with excluded volume, and the details of the model, as well as other simulation details, may be found in Ref. [27]. In Fig. 4 [28] the shear and elongational viscosities of 256 20-site chains are plotted against the second scalar invariant of the strain rate tensor, I2 . I2 is effectively a measure of the rate of viscous heat dissipation and is a natural quantity to compare shear and elongational viscosities, at least in the linear regime [29], though it has been shown to be the case for simple atomic fluids even in the nonlinear range if there is no first normal stress difference for shear flow [30]. The elongational viscosities are calculated for simulations run at constant volume and constant pressure. In an actual physical experiment, one free surface of the polymer is exposed to the atmosphere. For mechanical equilibrium to be maintained, the pressure normal to the surface (e.g., Pzz ) must be constant and equal to atmospheric pressure. It may thus be more useful to perform simulations at constant pressure, rather than constant volume, for comparison with experimental data. There are three particularly interesting features in these results:

(i) The Newtonian regime, where viscosity is independent of I2 , is not yet reached even for the weakest strain and elongational rates used. As the onset of nonlinear effects occur at weaker rates for longer chain systems, simulations at weaker fields need to be performed to capture the full Newtonian behavior of the fluid. Unfortunately statistical fluctuations in the shear and normal stresses also increase with decreasing I2 , making direct calculation of the viscosities notoriously unreliable. A likely solution to this problem is to use nonlinear response theory to calculate the shear stress. Such calculations have been used in the past to predict the weak field dependence of the shear and elongational viscosities of simple atomic fluids with impressive statistical accuracy [11,31]. (ii) The shear viscosity exhibits expected shear thinning with increased strain rate, whereas the elongational viscosity, defined as ηE = (Pyy − Pxx )/ 4˙ε, first thickens and then thins for both constant V and constant P simulations. However, for constant volume simulations the elongational viscosity again thickens at higher elongation rates. (iii) Constant pressure simulations do not exhibit a second region of thickening. The viscosity profile for constant pressure simulations is in qualitative agreement with experimental data in the existing literature [32], suggesting that indeed this is a more physically realistic simulation ensemble than constant volume. A fuller discussion of these results and others not yet published will appear in the literature shortly [33]. We are currently performing detailed NEMD simulations to study the Newtonian and non-Newtonian rheology of polymer melts as functions of molecular architecture and flow field, with the ultimate goal of using such simulations as a basis for the rational design of real polymer materials.

5. Conclusion Nonequilibrium molecular dynamics has been described with the specific intention of demonstrating its usefulness to, and complementarities with, computational fluid dynamics. It was shown how simulations of bulk homogeneous isothermal shear and elongational

B.D. Todd / Computer Physics Communications 142 (2001) 14–21

flows lead to a simple extraction of the shear and elongational viscosities, which in turn could be fed as input in a CFD calculation. The complementary nature of NEMD to CFD was demonstrated through an application to planar Poiseuille flow. In this example it was shown that for a microscopically confined fluid the predicted Navier–Stokes temperature profile does not agree with the simulation results. The theory also fails to give an accurate estimate of the thermal conductivity at the center of the channel. It is suggested that Fourier’s law of heat conduction is an incomplete constitutive equation to use in the energy conservation equation. An additional term coupling the heat flux vector to the gradient of the square of the strain rate leads to a correctly predicted temperature profile and thermal conductivity. This new transport phenomenon may be very important for correctly determining the heat transfer of microscopically confined fluids under non-uniform shear, and much work remains to be done for a greater understanding of its relevance. Finally, the use of NEMD in the study of the molecular rheology of polymer melts was briefly described. This work is actively being undertaken by several groups in the world today, and holds the promise of a powerful tool for the rational, costeffective design of new polymer based materials.

Acknowledgements The author wishes to thank Dr P.J. Daivis and Mr M. Matin for the use of unpublished data in Fig. 4.

References [1] K.P. Travis, P.J. Daivis, D.J. Evans, J. Chem. Phys. 103 (1995) 1109. [2] K.P. Travis, P.J. Daivis, D.J. Evans, J. Chem. Phys. 103 (1995) 10 638; Erratum, J. Chem. Phys. 105 (1996) 3893.

21

[3] R. Edberg, D.J. Evans, G.P. Morriss, J. Chem. Phys. 84 (1986) 6933. [4] R. Edberg, G.P. Morriss, D.J. Evans, J. Chem. Phys. 86 (1987) 4555. [5] G.P. Morriss, D.J. Evans, Comput. Phys. Commun. 62 (1991) 267. [6] B.D. Butler, G. Ayton, O.G. Jepps, D.J. Evans, J. Chem. Phys. 109 (1998) 6519. [7] L. Lue, D.J. Evans, Phys. Rev. E (2000). [8] M. Doi, S.F. Edwards, The Theory of Polymer Dynamics, Oxford University Press, Oxford, 1986. [9] J.H. Irving, J.G. Kirkwood, J. Chem. Phys. 18 (1950) 817. [10] A. Baranyai, D.J. Evans, P.J. Daivis, Phys. Rev. A 46 (1992) 7593. [11] D.J. Evans, G.P. Morriss, Statistical Mechanics of Nonequilibrium Liquids, Academic, London, 1990. [12] A.W. Lees, S.F. Edwards, J. Phys. C 5 (1972) 1921. [13] J.D. Moore, S.T. Cui, P.T. Cummings, H.D. Cochran, AIChE J. 43 (1997) 3260. [14] M. Kroger, S. Hess, Phys. Rev. Lett. 85 (2000) 1128. [15] J.D. Weeks, D. Chandler, H.C. Andersen, J. Chem. Phys. 54 (1971) 5237. [16] B.D. Todd, D.J. Evans, Phys. Rev. E 55 (1997) 2800. [17] D.J. Evans, Phys. Lett. A 91 (1982) 457. [18] P. Cordero, D. Risso, Physica A 257 (1998) 36. [19] G. Ayton, O.G. Jepps, D.J. Evans, Mol. Phys. 96 (1999) 915. [20] P.J. Daivis, J.L.K. Coelho, Phys. Rev. E 61 (2000) 6003. [21] Proceedings of 13th International Congress on Rheology, British Society of Rheology, 2000. [22] A.M. Kraynik, D.A. Reinelt, Int. J. Multiphase Flow 18 (1992) 1045. [23] B.D. Todd, P.J. Daivis, Phys. Rev. Lett. 81 (1998) 1118. [24] B.D. Todd, P.J. Daivis, Comput. Phys. Commun. 117 (1999) 191. [25] A. Baranyai, P.T. Cummings, J. Chem. Phys. 110 (1999) 42. [26] B.D. Todd, P.J. Daivis, J. Chem. Phys. 112 (2000) 40. [27] M. Matin, P.J. Daivis, B.D. Todd, J. Chem. Phys. 113 (2000) 9122. [28] M. Matin, P.J. Daivis, B.D. Todd, unpublished data. [29] M.N. Hounkonnou, C. Pierleoni, J.-P. Ryckaert, J. Chem. Phys. 97 (1992) 9335. [30] A. Baranyai, P.T. Cummings, J. Chem. Phys. 103 (1995) 10 217. [31] B.D. Todd, Phys. Rev. E 56 (1997) 6723. [32] H.A. Barnes, J.H. Hutton, K. Walters, Introduction to Rheology, Elsevier, Amsterdam, 1989. [33] M. Matin, P.J. Daivis, B.D. Todd (to be published).