Advances in nonequilibrium molecular dynamics simulations of lubricants and additives

Nonequilibrium molecular dynamics (NEMD) simulations have provided unique insights into the nanoscale behaviour of lubricants under shear. This review discusses the early history of NEMD and its progression from a tool to corroborate theories of the liquid state, to an instrument that can directly evaluate important fluid properties, towards a potential design tool in tribology. The key methodological advances which have allowed this evolution are also highlighted. This is followed by a summary of bulk and confined NEMD simulations of liquid lubricants and lubricant additives, as they have progressed from simple atomic fluids to ever more complex, realistic molecules. The future outlook of NEMD in tribology, including the inclusion of chemical reactivity for additives, and coupling to continuum methods for large systems, is also briefly discussed.


Introduction
In order to reduce energy consumption in engineering systems and thus CO 2 emissions, better lubricants need to be developed [1]. As in many areas of science and engineering, computer simulation has become a powerful tool which could help to accelerate future lubricant design. In tribology, processes at the smallest scales often drive the macroscopic friction [2] and wear [3] behaviour that we observe. As a result, there is growing interest in the application of molecular simulation to tribological problems. This review provides a historical and critical overview of the most important developments in this field of study, highlights some of the outstanding issues currently faced by researchers, and suggests future areas of development.
Molecular simulation, a field of research that started in the 1950s, uses computers to generate a series of representative snapshots of a group of interacting molecules from which experimental observables can be calculated. A system of interest is simulated by a collection of individual molecules, represented at the atomic level, and the locations of the individual atoms are governed by the interaction force laws between them. There are two main methods for performing molecular simulation, one stochastic and one deterministic [4]. The first to be developed was (Metropolis [5]) Monte Carlo (MC), which evolves the molecules by moving (usually) one molecule at a time by a trial random displacement. This new position is accepted or rejected on the basis of the change of potential energy of the system (more specifically its exponential or the "Boltzmann factor"). The second, and nowadays much more widely used method is molecular dynamics (MD), in which the atomic trajectories are generated by solving Newton's equations of motion using a finite difference scheme over a series of short time steps. Both MC and MD can be used to compute thermodynamic and other static properties such as infinite or zero frequency elastic moduli. Time-dependent processes and transport coefficients such as the viscosity arise more naturally from MD simulations, making it the more appropriate method to investigate tribological phenomena.
In tribology and specifically for liquid lubricants, there is considerable interest in predicting bulk transport properties of fluids (particularly the viscosity) under a wide range of conditions (temperature, pressure, shear rate) relevant to real components. As early as 1969, Goldstein [6] envisaged the use of MD to study the viscosity of liquids in the following prescient remarks: "The modern molecular approach to viscosity proceeds via the formulation of the shear viscosity in terms of appropriate integrals over the stress-time correlation function in a fluid. There is no difficulty, in principle, in imagining computer experiments generating this correlation function in a finite sample of molecules set into motion according to Newton's laws. It is not yet clear how large a sample of molecules is needed, nor how long a time of study is required, to achieve the equivalent of an infinite number of molecules studied for an infinite time, but one might guess that such a calculation will be within our reach in the immediate future, for liquids having the simplest, spherically symmetric, pair potentials. " Goldstein was percipient and the first MD study of the viscosity of a real fluid (argon) followed only four years later [7]. Very significant progress has been made since then and the number of studies applying MD to tribological systems has increased dramatically in recent years, for a variety of reasons. Firstly, reliable simulation methodologies have now been developed which allow researchers to study the shear of bulk fluid systems, as well as those confined by solid surfaces. Secondly, most lubricating oils are composed of relatively large organic molecules (C 20 −C 40 ) and force-fields that govern the interactions between such molecules in MD simulations have now been refined to a stage where they can reproduce realistic viscous behaviour [8][9][10]. Finally, recent increases in computational performance coupled with the development of highly parallelised simulation algorithms to exploit multiprocessor, high performance computing (HPC) architectures, have enabled larger molecules, system sizes, and timescales to be modelled. Nonequilibrium molecular dynamics (NEMD), where a shear force is applied to the system, has significantly improved our understanding of both dry friction, where the sliding surfaces are in intimate contact, and wet friction, where the surfaces are separated by a liquid lubricant. The focus of this article is NEMD simulations of wet friction, with a particular emphasis on (non-aqueous) lubricants and additives. In 2004, Urbakh et al. summarised experimental and theoretical contributions to the fundamental understanding of both wet and dry friction [2]. Szlufarska et al. reviewed experimental and theoretical studies of wet and dry friction specifically in the area of nanotribology in 2008 [11]. In 2013, Vanossi et al. [12] summarised advances in nanoscale to mesoscale modelling of wet and dry friction while and in the same year Dong et al. [13] discussed NEMD studies of dry friction in the context of atomic force microscopy (AFM) experiments. In 2014, Sawyer et al. [14] provided a more general overview of recent advances in the understanding of dry friction, including the contribution of NEMD. The reader should refer to these contributions for specific details on these topics.
How did we get to this point? We will now briefly summarise the stages that the subject has gone through to arrive at the present state where molecular simulation has become an important and powerful tool in the field of tribology (Section 1.1). We then provide some general guidelines for NEMD simulations of lubricants and additives (Section 1.2). We then describe a number of key step-changing methodological advances which are often taken for granted, but have revolutionised the subject (Section 1.3).

History of NEMD simulations
MD was invented in the 1950s, using the hard sphere (HS) model in which the "molecules" were represented as perfectly elastic and infinitely repulsive spheres ("billiard balls"). The HS MD approach is very computationally efficient (particularly at low density), and would today be referred to as an "event driven" method, as it follows the collisions between the spheres generated by classical mechanics in chronological order. Probably the most significant paper of this pioneering work was by Alder and Wainwright [15], who in 1957 showed that an HS liquid crystallised as the density was increased above a certain value. This discovery had wide-ranging consequences, as it was widely thought at the time that it was necessary for molecules to have an attractive component to their interaction potential to show a freezing transition and Alder and Wainwright's work showed that this was not the case. Today this would be called an "entropydriven" freezing transition, of which there are abundant examples in colloid and polymeric liquid systems, and its existence plays a major role in our understanding of the phase behaviour and physical properties of a wide range of soft condensed matter chemical systems [16].
The 1960s was the period in which model molecules represented by a continuous interaction potential were first simulated by MD. The first paper was by Rahman in 1964 [17], who calculated the diffusion coefficient in a system of argon atoms interacting through a Lennard−Jones (L−J) potential. This was also a groundbreaking paper in that it showed the diffusion and structural evolution of small molecules takes place by a series of small, highly coordinated motions of neighbouring molecules. Until this MD paper, it had been generally held that diffusion in liquids took place through a series of rare, large distance jumps of molecules from "lattice sites" into vacancies or "holes" in a crystal-like lattice, the jump distance being roughly a molecular diameter. Historically, liquids were viewed as a crystal with many unoccupied sites [18], but Rahman's paper convincingly disproved this supposition. However, erratic, localised displacements have more recently been observed in the MD simulations of supercooled liquids near the glass transition, a phenomenon known as "heterogeneous dynamics" [19], so the defect-containing crystal view of a liquid is not entirely without merit. The diffusion coefficient is closely related to the viscosity, indeed the Stokes-Einstein equation [20] suggests that the product of the two is almost constant in liquids [21]. In 1970, Alder et al. [22] calculated the viscosity of the HS fluid at various densities up to crystallisation. They used a mean square displacement of the shear stress route, which is now called the Einstein-Helfand method [23].
Verlet and co-workers in 1973 [7] were the first to calculate the shear viscosity of a real liquid (argon) represented by the Lennard−Jones (L−J) potential. They used the Green-Kubo time correlation method [24,25], where the Newtonian viscosity is defined as the integral of the shear stress relaxation function, C s (t). In the small shear rate limit, C s (t) can be shown by statistical mechanics to be the product of the shear stress at one time multiplied by its value at a later time. This function decays monotonically to zero, as liquids subjected to a step in strain can only sustain the initial jump in shear stress for a finite time. Therefore, perhaps paradoxically at first sight, the Newtonian viscosity can be calculated without actually imposing a shear rate across the system. Near the melting line, the viscosity of the HS fluid increases rapidly because of a slowing down in the decay of the shear stress time correlation function (caused by crowding effects of the molecules packed closely together), an effect that was called the "molasses tail" [26].
Periodic boundary conditions (PBCs) are a very important tool which allow MD simulations to study essentially bulk systems in a computationally efficient manner. PBCs create an infinite lattice of periodically replicated molecular simulation cells which fills all space ( Fig. 1) [4]. If a molecule leaves the cell through one boundary, it is reintroduced with the same velocity at the same point of departure through the opposite face of the cell. This construction effectively eliminates surface effects that would otherwise dominate the behaviour of the system, due to the large surface-tovolume ratio of a typically-sized simulation cell [4].
Until the 1970s, the liquid simulations were carried out mainly using equilibrium MD which by default follows the microcanonical ensemble (an ensemble is simply a probability distribution for the state of the system). The microcanonical ensemble is also sometimes called the NVE ensemble because the total number of particles in the system, N, the volume of the system, V, as well as the total energy in the system, E, remain constant. Then, in the early 1970s, the first nonequilibrium molecular dynamics (NEMD) simulations of imposed shear flows at large shear rates were carried out. A variety of techniques were developed, and these continue to be invented and refined today (Section 1.3). The version perhaps most relevant to tribology is the "sliding wall" method invented by Ashurst and Hoover in 1975 [27]. This mimics the shearing of a fluid of infinite extent in the x-and y-directions using PBCs, confined by bounding walls whose plane normal is in the z-direction (see Fig. 1(a)). The two walls (which were initially simply the outermost regions of the fluid) are constrained to translate in opposite directions in the x-direction to shear the fluid. This was followed by simulations by Bitsanis et al. in 1987 in which solid walls were used to confine and shear the fluid [28]. The study of confined fluids under shear by NEMD has undergone a rapid expansion in recent years, and is now a common geometry for simulations in many different fields of research.
Up to the end of the 1970s, NEMD simulations were quasi-2D in that they used PBCs in only two directions ( Fig. 1(a)), but in the 1980s equations of motion for imposing shear flow in systems with PBCs in all three Cartesian directions were invented ( Fig. 1(b)). These are known as the DOLLS and SLLOD methods [30]. By good fortune, PBCs which can accommodate shear flow of infinite extent within MD had been invented by Lees and Edwards (LE) several years previously [31], and these new non-Hamiltonian equations of motion are compatible with these "sliding brick" PBCs (see Fig. 1(b)). Recent books by Evans and Morriss [32], and Todd and Daivis [33], provide a comprehensive description of the field from that time to the present.
An early NEMD simulation using the LE PBCs, performed by Heyes et al. in 1980 [34], modelled the start-up of shear flow in a supercooled Lennard−Jones (L−J) fluid. In these simulations, "stress overshoot" behaviour was observed shortly after the initiation of shear, a phenomenon widely seen in experimental studies of polymeric and glassy systems. The importance of the 1980s cannot be overestimated as it was then that the statistical mechanics or governing physics of liquid systems arbitrarily far from equilibrium started to be derived, largely by Evans and Morriss [32]. Since the 1980s, NEMD has progressed to a point where it has become a predictive tool in tribology, as described in Section 2.

Computational guidelines
In order for NEMD simulations to be useful in terms of predictive capabilities, a number of aspects must be considered, from the inherent limitation of the time and length scales that can be modelled, to the accuracy of the representation of the intramolecular and intermolecular interactions [35]. The computational time required to conduct MD simulations depends on a number of factors including the computer hardware and software available, the system size, as well as the complexity of the force-fields used. The number of atoms required in NEMD simulations of tribological systems varies depending on the particular problem being studied, but is generally in the thousands to hundreds of thousands range. Classical MD forcefields require a time step of approximately 1 fs to ensure energy conservation [4]. This means that, even when performed on multiprocessor high performance computing (HPC) systems, using highly parallelised MD software, only ns or in extreme cases μs timescales are accessible in MD simulations [36], since these correspond to millions/billions of time steps.
In NEMD simulations, where shear is applied, the main consequence of the short accessible timescales is that relatively high shear rates are required (generally > 10 7 s −1 ) to ensure that the properties of interest (most commonly the viscosity) reach a steady state [37,38]. Moreover, at very low shear rates, the systematic nonequilibrium response becomes comparable to the equilibrium fluctuations in the phase variables of interest, resulting in a low signal-to-noise ratio [39,40]. The ability to simulate lower shear rates has been a long-term goal of NEMD simulations to facilitate direct overlap with experiments and real components [32]. For comparison, the high pressure viscosity of lubricants can generally only be measured up to shear rates of approximately 10 4 s −1 [37,38], tribology experiments can reach up to around 10 6 s −1 [41], and high-performance engine components can extend up to 10 8 s −1 [42].

Methodological developments
Following the initial development of bulk and confined NEMD simulation methodologies (Section 1.1), many aspects have been altered and refined to allow different conditions to be studied and more accurate results obtained. Significant examples include improved temperature and pressure control methods, classical force-fields, as well as new simulation methods; these are briefly discussed below.

Temperature and pressure control
The high shear rates that are usually imposed in NEMD simulations create heat, which must be dissipated to achieve a nonequilibrium steady state [32]. In parallel with the theoretical advances described in Section 1.1, a number of different thermostats have been invented which can be used in conjunction with the NEMD equations of motion. In bulk MD simulations, the thermostat is usually applied globally (to all of the atoms in the system). The temperature in MD simulations can be controlled through a number of methods including scaling the atomic velocities (e.g., Andersen [43]), coupling atoms to an external bath (e.g., Berendsen et al. [44]), or applying a random "friction" force to the atoms (e.g., Langevin [45] and dissipative particle dynamics (DPD) [46]). The Nosé-Hoover (NH) thermostat [47,48], in which a Hamiltonian with an extra degree of freedom for a heat reservoir is introduced to control the temperature, is the most widely used thermostat for bulk NEMD simulations. Under equilibrium conditions, the NH thermostat can be proven to generate a canonical or NVT ensemble. The NH thermostat can also be easily implemented in conjunction with the SLLOD equations of motion [49]. However, some researchers have raised concerns with thermostats, such as NH, which require a predefined flow profile (e.g., linear/Couette) to be assumed [50]. So-called profile-biased thermostats can result in artefacts in some systems, for example the "string phase" which was often observed in NEMD simulations of L−J fluids at high shear rates [51]. Moreover, they are of course incorrect for systems where the flow profile is not Couette-like. To overcome these problems, configurational thermostats, which make no assumptions regarding the flow profile, were developed [52,53]. These have not yet been applied in bulk NEMD simulations of lubricant-sized molecules (C 20-40 ), probably because they are more difficult to implement for large molecular systems, so the standard NH thermostat is typically still employed [37]. However, configurational thermostats may be useful for bulk NEMD simulations under high pressure and high shear rate conditions, where deviations from linear velocity profiles have been observed experimentally (see Section 2.2.4).
For confined NEMD simulations, alternative thermostatting strategies are required; Yong and Zhang [54] evaluated some of these in 2013. In confined NEMD simulations, the thermostat is usually only applied to the wall atoms [54], allowing a thermal gradient to develop through the thickness of the lubricant film, as occurs in real experiments [55]. The | https://mc03.manuscriptcentral.com/friction localised thermostat is often only applied in the Cartesian directions perpendicular to the shearing direction in order to avoid significantly affecting the dynamics. Applying the thermostat directly to the confined fluid molecules has been shown to artificially influence their structure, flow and friction behaviour under shear [56,57]. For example, the slip behaviour of thermostatted molecules confined between rigid walls disagrees with that observed experimentally and in simulations where only flexible wall atoms are thermostatted [57] (see Section 2.2.3). In confined NEMD simulations of tribological systems, stochastic thermostats (e.g., Langevin [45]) are generally more popular than the NH thermostat since they can remove heat more efficiently [58]. For the relatively thin films usually studied in confined NEMD simulations, the temperature rise, which is generally largest in the centre of the film, only becomes significant at high shear rates [58]. More specialised thermostatting strategies have also been developed, for example, thermostats can be applied to "virtual" wall atoms rather than the wall atoms themselves [59]. This technique is particularly useful for large, complex systems where the wall atoms can be "frozen" to preserve structural stability and reduce the computational cost [59]. Stochastic boundary conditions (SBC) [60], in which the border atoms experience random friction forces (Langevin) to control the flow of energy through the system border, have also been shown to be effective in temperature control of confined NEMD simulations of tribological systems [61].
Pressure control is also important for NEMD simulations of tribological systems. In order to mimic experimental conditions, in bulk NEMD simulations, pressure is most commonly controlled using the Nosé-Hoover (NH) barostat [47,48], which continuously alters the unit cell vectors, and thus the cell volume, to maintain a target pressure. This method can be used to sample the isothermal-isobaric or NPT ensemble [49]. For confined NEMD simulations, pressure can be controlled by setting a fixed separation distance between the confining walls, or by applying a constant normal force to one or both of the walls (e.g., Fig. 1(a)) [62]. The latter technique is usually preferable since the former has been shown to lead to unphysical behaviour in certain systems [62]. For incompressible systems, such as L−J fluids, it can also be beneficial to add a damping term to the barostat in order to prevent excessive oscillations after the normal force is applied [62].
An alternative approach for confined NEMD simulations is grand-canonical molecular dynamics (GCMD) [63,64], in which the chemical potential, temperature, and pressure are held constant. This technique was developed from GCEMC (grandcanonical ensemble Monte Carlo) [65] to investigate interesting results from surface force apparatus (SFA) experiments of confined liquid films [66]. In GCMD simulations, a reservoir of fluid molecules surrounds the confined system [64], thus allowing molecules to be "squeezed out" from between the solid surfaces in one or more directions at high pressure [67]. GCMD is still less widely used than conventional confined NEMD, primarily because of the additional computational expense to simulate the reservoir molecules; however, it has recently been applied to study the behaviour of ionic liquid lubricants under shear [68]. Liquid−vapour molecular dynamics (LVMD) [69], in which the confined fluid is surrounded by a vapour rather than a liquid, has also been developed in order to trivialise the calculation of the pressure tensor from GCMD simulations. LVMD has been shown to yield more realistic squeeze-out behaviour compared to GCMD [70], and has been successfully applied to study the stick-slip behaviour of atomic fluids under shear [71] (Section 2.2.2).

Classical force-fields
A common concern in classical MD simulations is that their accuracy is heavily dependent on the force-fields used to describe the intermolecular and intramolecular forces between a system of interacting molecules. The functional forms of most classical molecular forcefields are quite similar, capturing both bonded (bond stretching, bending and torsion) and non-bonded (van der Waals, electrostatic) interactions [10]. These terms can be empirically parametrised to reproduce important experimental properties for a set of training compounds, derived from first principles calculations, or computed using a combination of the two. Force-field parameterisation is a complex and time-consuming process, and most of the force-fields used to study liquid lubricants were originally developed for biological applications [10]. The computational expense of classical MD simulations is usually dominated by modelling the van der Waals and particularly the long range electrostatic interactions between atomic partial charges. The former is usually modelled by a Lennard− Jones (L−J) potential and is generally "cut off" at certain distance where the energy tends to zero (≈ 10 Å) for computational efficiency. Conversely, electrostatic interactions should not be cut off, but rather treated with a long range solver such as the Ewald summation [72] or particle−particle, particle−mesh (PPPM) [73] methods, since truncation can lead to unphysical results [74]. Note that confined NEMD simulations require a slab implementation of these algorithms [75].
A key target of NEMD simulations of tribological systems is to yield realistic viscous behaviour of lubricant-sized molecules, which requires the use of highly accurate force-fields. Quantitative comparisons between several popular force-fields for NEMD simulations of lubricants can be found elsewhere [10,76,77], but a brief overview is also provided here for completeness. Most historic simulations of tribological systems which include alkane molecules have employed united-atom (UA) force-fields where the nonpolar hydrogens are grouped with the carbon atoms to generate CH, CH 2 and CH 3 pseudo-atoms [10]. The popularity of UA force-fields has stemmed from the fact that they decrease the number of interaction sites by around 2/3 and computational expense by up to an order of magnitude compared to all-atom (AA) forcefields [78]. However, UA force-fields have been shown to lead to dramatic viscosity under-prediction for linear, long-chain alkanes (≈ 50% for n-hexadecane [10]) compared to experiment. For alkanes with multiple short branches, UA force-fields can give reasonably accurate viscosity results (≈ 15% for squalane [79]), but for those with fewer, longer branches they are less accurate (≈ 50% for 9-n-octyldocosane).
AA force-fields are much more computationally expensive and, as originally developed, many of them resulted in anomalous solid−liquid phase behaviour for the long-chain alkanes relevant to tribology [80]. Fortunately, the latter issue has been overcome by re-parameterising AA force-fields specifically for long-chain alkanes [8,9], which has led to far more accurate density and viscosity prediction (≈ 10% for n-hexadecane [10]), even under high temperature and high pressure (HTHP) conditions [10]. More modest improvements in the density prediction of lubricantsized alkanes using AA force-fields can be made simply by reducing the intramolecular 1-4 scaling parameters for the L−J and electrostatic interactions using otherwise unmodified force-fields [80]. An increasing number of tribological NEMD studies are now employing accurate AA force-fields in order to yield more realistic viscous behaviour; however, for some types of molecules (multiple small branches, e.g., squalane) the additional computational expense compared to UA force-fields may not be justified.

Modified NEMD algorithms
Many modifications of the original NEMD methodologies have been made over the last three decades to improve their signal-to-noise ratio, simplify their implementation, and extend them to study different experimental conditions. Noteworthy in this set of theoretical tools is the so-called "transient time correlation function" (TTCF) method [39,40], derived by Morriss and Evans in the late 1980s. TTCF can be used to generalise the Green-Kubo method [24,25], which strictly is only applicable to equilibrium systems, to systems forced out of equilibrium by the imposition of an external field (e.g., shear rate). TTCF involves averaging correlation products from many independent NEMD trajectories, making it rather computationally expensive. However, the main advantage of TTCF is that it allows very low shear rates to be simulated compared to NEMD, allowing direct overlap with those used in experiments and real components [42]. Indeed, it has been successfully applied to study the viscosity of short alkanes (n-decane) at shear rates as low as 10 5 s −1 [81]. TTCF was initially used to study the rheology of bulk systems but it has also been used to study friction in confinement [82,83]. The extent to which the TTCF can be used to calculate the friction coefficient of confined fluids has been much discussed, and alternative methods have also been proposed [84][85][86].
In 1999, Müller-Plathe [87] introduced an alternative nonequilibrium method for calculating the shear viscosity, so-called "reverse" NEMD or R-NEMD. The R-NEMD method reverses the cause-and-effect picture customarily used in NEMD in which the velocity gradient or shear rate is imposed and the momentum flux or stress determined, to a stress being imposed and the consequent velocities observed. The method involves a simple exchange of particle momenta between boundary layer molecules, which is straightforward to implement. Moreover, it can be made to conserve the total energy as well as the total linear momentum, so no coupling to an external temperature bath is needed. The method was initially tested for an L−J fluid near its triple point and yielded a viscosity in agreement with literature results [87]. The method has since been used to predict the viscosity behaviour of simple alkanes [88] as well as more complex multicomponent systems such as lipid bilayers [89], ionic liquids [90], and polymer viscosity modifier additives [91,92]. R-NEMD simulations are simple to implement, versatile and can be used to study transport properties in bulk or confined systems. However, unlike conventional NEMD, the external forces are not applied to the system in a manner consistent with experiments in R-NEMD simulations, and it has also been shown to under-predict the viscosity compared to conventional NEMD at high shear rates [93].
Another notable advance was the extension of NEMD methodology to simulate planar elongational flow (PEF) by Todd and Daivis in 1998 [94]. This is important in tribology since PEF occurs in the entrance region of high pressure, elastohydrodynamic contacts [95]. PEF simulations have been widely used to study polymers flow [96,97] but more rarely for lubricant-sized molecules. Baig et al. [98] studied the rheological and structural properties of linear liquid n-alkanes (C 10 , C 12 , C 20 ) using NEMD simulations under PEF. They showed that these molecules behave very differently to the way they do under shear flow, with strong alignment of fully stretched chains at high elongation rates suggesting the formation of a liquidcrystal-like, nematic structure [98]. However, it should be noted that concerns have raised with the validity of the p-SLLOD algorithm used in these simulations [99]. Hunt et al. [100] developed an algorithm capable of imposing mixed Couette flow and PEF. Furthermore, Hartkamp et al. [101] applied TTCF to these mixed flows in order to capture the behaviour of fluids at lower shear rates.
Another method is driving spring-block NEMD, which was developed by Thompson and Robbins [102] to study the stick-slip behaviour observed in SFA experiments [103] (Section 2.2.2). In this technique, the sliding velocity is not applied directly to the wall atoms, but rather through a stage (or virtual atom) linked to the wall atoms via a spring. As the stage is moved, the spring is stretched and, if the confined film is liquid, the top plate accelerates until the steadystate spring force balances the viscous dissipation. If the film has crystallized, stick-slip behavior is observed; initially, the film responds elastically and the wall remains stuck, but eventually the spring force exceeds the yield stress of the film, and the top plate slides until the plate re-sticks and the process repeats. This technique has since been used to mimic SFA experiments of a range of strongly confined systems, including atomic fluids [71] and water [104].

Summary
Since their conception, NEMD simulation methodologies have been continually varied and refined in order to more faithfully recreate experimental conditions and provide more accurate results. Significant progressions have been made in temperature and pressure control, force-field accuracy, as well as the simulation algorithms themselves to allow a more diverse range of tribology problems to be investigated more accurately. Despite considerable advances, classical MD simulations are still limited in both the accessible time and length scales and their inability to study chemical reactivity. Thus, in the future, the most interesting methodological advances are envisaged to be in coupling MD simulations to smaller (quantum) and larger (continuum) scales; this will be discussed in Section 4.

NEMD simulations of liquid lubricants
Many tribology problems are well-suited to study by NEMD simulations. Shear can be applied to systems of lubricants and additives to give insights into their flow and friction behaviour. As discussed in Section 1, shear forces can be added either directly to the fluid molecules or by momentum transfer from sliding surfaces which confine the fluid (Fig. 1). While there are clearly significant differences between NEMD simulations of strongly confined, inhomogeneous systems and the bulk [105][106][107] (as justified below), comparisons between the viscosities from these methods for homogenous flows have shown good agreement [55,108,109]. Many lubricant properties, both in the bulk and under confinement, have been investigated using NEMD simulations, providing information at a scale that is difficult to access experimentally [4]. Indeed, NEMD simulations have predicted some important lubricant behaviour many years before experimental techniques capable of providing comparisons on the same scales were developed. The following sections will discuss the use of NEMD to investigate the behaviour of lubricants in the bulk (Section 2.1), as well as when confined by solid surfaces (Section 2.2).

Bulk NEMD simulations of liquid lubricants
Bulk NEMD simulations are frequently used to investigate the change in viscosity with important variables such as the shear rate, temperature and pressure. As the shear rate is increased, lubricants often undergo a transition from a Newtonian regime (where the viscosity, η, is independent of shear rate,  ) to shear thinning behaviour (where η decreases with increasing  ). The critical shear rate transition,  c , typically occurs at      1 , where  is the longest relaxation time at equilibrium (in the absence of shear). This is usually the rotational relaxation time,  rot [37,38], which increases with the chain length of the molecule [110]. Thus, longer simulations are required to obtain the Newtonian viscosity from NEMD simulations of larger molecules. As a result, early bulk NEMD studies were limited to studying the rheology of atomic Lennard− Jones (L−J) fluids and later short chain alkanes (≤ C 10 ) using UA force-fields [111]. Only more recently have they been applied to study lubricant-sized molecules (C 20-40 ). Some significant examples of bulk NEMD simulations of lubricant-sized molecules are discussed below. A list of bulk NEMD simulations applied to a wider range of fluids was compiled by Todd and Daivis [112] in 2007.

Viscosity prediction from bulk NEMD
Morriss et al. [113] were the first to study lubricantsized molecules with bulk NEMD simulations in 1991. They used a UA model of n-eicosane (C 20 ) with a purely repulsive, Weeks, Chandler, Andersen (WCA) potential [114] and, as is common in NEMD simulations, they progressively increased the shear rate to study its effect on viscosity. At the lowest shear rates tested they observed a shear thinning regime, followed by shear thickening at higher shear rates. However, this shear thickening regime was later shown to be an unphysical consequence of the simplified force-field employed [49]. Moreover, the lowest shear rates they used were not low enough to capture the Newtonian plateau.
In 1998, Khare et al. [115] became the first to use NEMD to investigate the rheological behaviour of lubricant-sized alkanes using a UA force-field which included repulsive and attractive interactions through a Lennard−Jones (L−J) potential [116]. They were able to simulate lower shear rates and thus capture the Newtonian plateau for several linear (C 16 , C 22 , C 28 ) and branched (5,12-dipropylhexadecane) alkanes. The Newtonian viscosity of n-decane from the NEMD simulations was in good agreement with experiment; however, the viscosity results for longer alkanes were somewhat below the corresponding experimental values (e.g., 44% for C 28 ). However, the UA force-field employed was able to capture qualitatively the change in viscosity with lubricant molecular structure.

Viscosity temperature and pressure dependence
More recently, several bulk NEMD simulations of lubricant-sized alkanes focussed on calculating the Newtonian viscosity and also probing its change with temperature (viscosity index, VI or viscosity number, VN) and pressure (α). The VI, VN, and α are important indicators of lubricant performance, so there is considerable interest in their prediction by molecular simulation.
In 2000, Moore et al. [79] performed bulk NEMD simulations on three C 30 isomers using a UA forcefield [117]. They found that, while the viscosity of the linear (48%) and single-branched (36%) isomers were significantly underpredicted, viscosity prediction for squalane (15%), which contains multiple short methyl branches along its backbone, was more accurate [79]. They also analysed the viscosity at different temperatures to calculate the VI, which was in excellent agreement with experiment for all of the isomers (≈ 10%) [79]. McCabe et al. [118] performed similar simulations to calculate the VN of 9-octylheptadecane, 9-octyldocosane, and squalane using a UA force-field [117]. The viscosity was again significantly underpredicted but the VN values were within the statistical error of the experimental results for all of the molecules studied.
McCabe et al. [119] and Liu et al. [120] performed NEMD simulations of 9-octylheptadecane and a model poly-α-olefin (PAO) base oil (1-decene trimer) respectively from ambient up to GPa pressure to predict their α values. They both used a UA force-field [117] and although the Newtonian viscosity was significantly underpredicted (30%−70%), the calculated α values were in good agreement with experimental results (≈ 10%). It should also be noted that Ramasamy et al. [121] accurately predicted α values for lubricant-sized molecules from their pressure-volume behavior in simple EMD simulations.

Shear thinning
Another focus of bulk NEMD simulations of lubricantsized moleules has been investigating their shear thinning behaviour. There is still much debate regarding the most suitable shear thinning model to describe lubricant behaviour under high pressure, high shear rate, and elastohydrodynamic lubrication (EHL) conditions [41,122,123]. NEMD simulations have been used to investigate this behaviour and can give unique insights into behaviour under conditions which are challenging to probe experimentally [124][125][126]. Some investigations of shear thinning behaviour in lubricantsized alkanes are discussed below.
In 1998, Cui et al. [110] used NEMD and a UA force-field [117] to study the rheological properties of a variety of linear (C 10 , C 16 , C 24 ) and branched (10-nhexylnonadecane, squalane) alkanes. At low shear rates, the viscosity behaviour transitioned to a Newtonian plateau and the zero-shear viscosity agreed with that from independent EMD simulations using the Green-Kubo method [24,25]. At high shear rates, the viscosity showed a power-law shear thinning behaviour over several orders of magnitude in shear rate. This powerlaw shear thinning was shown to be closely related to the ordering of the molecules and the molecular architecture had a significant influence on the powerlaw exponent [110].
In 2002, Bair et al. [37,38] compared bulk NEMD simulations of squalane using a UA force-field [117] with experimental measurements under EHL conditions in both the Newtonian and shear thinning regimes. This was the first comparison of the nonlinear rheology predicted by NEMD with experiment, and was thus the first experimental test of NEMD simulations in the shear thinning regime. The experimental and simulation data were separated by several orders of magnitude in shear rate; however, they collapsed onto the same time-temperature superposition master curve [127]. This master curve was successfully fit using the Carreau shear thinning model [128].
Liu et al. [129] proposed an approach to correlate shear thinning with the change in the molecular conformation, characterised by the radius of gyration of the molecule from NEMD simulations using a UA force-field [117]. This approach was tested by analysing the critical shear rate for squalane and 1-decene trimer, and then extended to study the behaviours of different molecular weight PAO structures. The changes in viscosity and radius of gyration with shear rate were adequately fit with the Carreau equation [128]. Their analysis led to a relationship between molecular weight and critical shear rate for PAO structures, and the calculated critical shear rates agreed well with the Einstein-Debye equation [130] for low molecular weight PAOs.
Recently, Jadhao and Robbins [131] used NEMD simulations to determine the nonequilibrium shear response of squalane over a wide range of pressures and temperatures using a UA force-field [117]. The viscosity-shear rate behaviour from the NEMD simulations agreed well with experimental data conducted at lower shear rates (10 2 −10 4 s −1 ). The Eyring shear thinning model [132] could be used to extract an accurate prediction for the Newtonian viscosity from NEMD simulations at high shear rates (> 10 5 s −1 ) as long as there was a substantial linear region in the plots of stress vs. log(shear rate). This condition was met whenever the Newtonian viscosity was ≥ 1 Pa·s, where the fluid was in a "glassy" regime. At lower viscosity, the response deviated from the simple Eyring form, indicating that the motion could not be described by a single activation barrier or rearrangement mechanism. The response in this "low viscosity" |www.Springer.com/journal/40544 | Friction http://friction.tsinghuajournals.com regime was better described by models with a wide range of characteristic times, such as the modified power-law in the Carreau shear thinning model [128].

Confined NEMD simulations of liquid lubricants
In addition to bulk lubricants, NEMD simulations have provided significant insights into the behaviour of thin lubricant films confined between solid surfaces. Confined lubricants behave very differently to their bulk counterparts under shear, and density inhomogeneity can lead to the breakdown of continuum approximations such as Navier-Stokes hydrodynamics [133]. Density oscillations in confined fluids were first identified in Monte−Carlo (MC) simulations of L−J liquids [65,134]. Shortly afterwards, surface forces apparatus (SFA) experiments by Horn and Israelachvili on confined octamethylcyclotetrasiloxane (OMCTS) [135] (a quasi-spherical molecule) showed that the force between two surfaces does not evolve smoothly with the separation distance as predicted by the van der Waals theory. Instead, in films below six to ten molecular diameter thickness, the wall-wall interaction oscillated between repulsion and attraction states, with a periodicity roughly equal to the molecular diameter. This behaviour, known as the "solvation force", was the first experimental evidence of density inhomogeneity in confined fluids [135]. Subsequently, the first MD simulations of confined fluids by Magda et al. [136] agreed with earlier MC and experimental results, showing that strongly confined films have inhomogeneous density profiles across the channel. NEMD simulations have also given many additional insights into the nonequilibrium behaviour of confined films, as discussed below.

Density and viscosity inhomogeneity
The first NEMD simulations of fluid films confined by walls were performed in 1987 by Bitsanis et al. [28]. They used the NEMD procedure developed by Ashurst and Hoover [27], but also included a continuous wall potential to confine an L−J fluid. They hypothesised that the inhomogeneous films could be viewed as a set of homogeneous layers and developed the local average density model (LADM), where each layer is given its own density and thus its own local viscosity. The LADM gave satisfactory results in predicting the flow of moderately-confined fluids [28]. However, the LADM model failed to predict many of the properties of strongly-confined films (≤ 4 molecular layers), most notably the enhanced viscosity, where the density variation alone cannot explain the local dynamic properties of the film [137].
In conventional rheology models, the resistance to flow of a lubricant is characterized by its shear viscosity,  . For bulk systems (Section 2.1),  depends purely on the fluid and its response to variations in pressure, temperature, and shear rate. Conversely, the dynamics of confined fluids are also influenced by factors that are not intrinsic to the fluid, so the definition of a bulk viscosity is not strictly appropriate [138]. However, the ratio of the shear stress to the shear rate still gives insight into the resistance to flow of confined films and can be compared directly to experiments [137]. Most experiments and NEMD simulations of confined films report this value as the "effective" shear viscosity,  eff .
When a fluid film is thicker than approximately ten molecular diameters, at which point its centre becomes a structureless medium, the experimentally measured  eff is usually within around 10% of the bulk viscosity [103,139]. Experiment and NEMD simulations have shown that when the wall separation is decreased below this level,  eff of confined L−J fluids [137,140], linear alkanes [141,142] as well as many other organic liquids [143] begins to increase significantly. At this point, continuum approximations such as Navier-Stokes hydrodynamics begin to break down [133]. When the film thickness is reduced to around six molecular layers,  eff has been shown to increase by several orders of magnitude for a range of fluids [141][142][143][144]. For L−J fluids, this viscosity enhancement has been attributed to the crystallisation of the entire film [145,146]. Compared to L−J fluids, linear alkanes are less susceptible to crystallisation and so either a second order glass transition [145] or the formation of "crystal bridges" [141] have been attributed to the viscosity enhancement. However, more recent MD simulations by Docherty and Cummings [147] using an all-atom force-field suggested that both cyclic (cyclohexane) and long-chain (n-dodecane) alkanes undergo a rapid, first order phase transition to a layered and ordered (i.e., crystalline) system (see at surface separations equivalent to fewer than nine molecular layers. These simulations support experiments which showed an abrupt transition from liquid-like to solid-like behaviour upon increasing confinement [143,144], rather than a more gradual transition to a glassy state (i.e., vitrification) [142,148]. It is noteworthy that most evidence for vitrification originated from SFA experiments which were later shown to be contaminated with platinum nanoparticles [149], which could go some way to explaining differences between experiments. A more in-depth discussion of the discrepancies between different experimental techniques and MD simulations of various levels of complexity was given in 2010 by Cummings et al. [150]. It has also been shown through experiments [148] and NEMD simulations [145] that elevated pressure can induce such transitions, as well as the associated increase in  eff , in thicker films.
The transition from liquid-like to solid-like behaviour in confined fluids also depends strongly on molecular shape and surface structure. MD simulations have shown that branched or flexible molecules are more resistant to solidification and remain liquid-like under more severe confinement and higher pressures compared to linear or inflexible molecules [151,152]. Thus branched molecules can give lower  eff under strong confinement than linear molecules, despite their higher bulk viscosity [103]. It should be noted most NEMD simulations of confined film behaviour employ atomically-smooth surfaces, but it has been found that surface roughness frustrates the ordering of the neighbouring fluid molecules so that films are more liquid-like on surfaces with nanoscale roughness com-pared to those which are atomically-smooth [153,154].
When the film is thinner or when the pressure is increased, the Newtonian plateau is displaced to lower shear rates. Thus, at the shear rates accessible to NEMD simulations (generally > 10 7 s −1 ), confined fluids are often highly shear thinning [155]. A viscosityshear rate power law was identified for L−J fluids at different pressures and thicknesses [145]. Both the Carreau [90] and Eyring [27] shear thinning models have been successfully used to describe the viscosityshear rate behaviour of confined fluids from NEMD simulations.
In stark contrast to the elevated viscosity results discussed above, some NEMD simulations [146] and SFA experiments [156] have reported an almost vanishing friction coefficient (superlubricity) in strongly confined lubricant films. This low friction was initially attributed to boundary slip (Section 2.2.3) at incommensurate solid−fluid interfaces, but subsequent NEMD simulations run for longer simulation times [157][158][159] have suggested that slip does not occur at the interface, but rather shear-induced alignment causes slip between well-defined molecular layers within the film itself. Surprisingly, there have been relatively few experimental or MD studies which have investigated this phenomena in the last decade, and reductions in the effective viscosity due to confinement certainly remain the exception rather than the rule.
NEMD results for the change in viscosity with confinement have also been used to improve the predictive capabilities of continuum models. For example, Martini et al. [160] used NEMD simulations to extract the film thickness and velocity dependence of the viscosity of n-decane confined and sheared between atomically-smooth gold surfaces. This information was then fed into a full numerical solution for EHL problems [161] to assess the impact of these viscosity changes on the predicted film thickness.
As well as the enhanced effective viscosity, confined NEMD has been used to study several tribological phenomena where the flow becomes inhomogeneous or discontinuous. Examples include stick-slip (Section 2.2.2), boundary slip (Section 2.2.3), and shear localisation (Section 2.2.4). These behaviours can arise directly as a result of strong confinement or from extreme applied pressure and shear rate conditions in thicker films [162]. NEMD simulations have and continue to play a central role in improving our understanding of all of these phenomena, which is critical to help design better lubricants.

Stick-slip
Stick-slip occurs when, instead of sliding over each other smoothly, solid surfaces alternately stick together and then slip past each other [102]. In 1990, SFA experiments of several non-polar fluids strongly confined and sheared between atomically-smooth mica surfaces showed a saw-tooth pattern in the friction coefficient over time [103]. This stick-slip friction response was also seen in spring-block NEMD simulations (Section 1.3.3) of a confined L-J fluid performed by Thompson and Robbins [102] the same year. The authors suggested that stick-slip motion involved periodic shear-melting transitions (slip) and recrystallization (stick) of the thin fluid film [102]. However, more recent spring-block, LVMD simulations of a confined L−J fluid by Lei and Leng [71] were in contrast with this shear-melting hypothesis, as they suggested that the ordered, solidified film persisted during both the stick and slip phases. They proposed that periodic slip events, either within the fluid film itself or at the wall-fluid interface (Section 2.2.3) were actually the cause of the stick-slip friction response [71]. Recent surface force balance (SFB) experiments of OMCTS from Rosenhek-Goldian et al. [163] supported the NEMD results of Lei and Leng [71], as they observed no dilation of the intersurface gap, which the authors suggest discounts the shear melting hypothesis. However, the resolution of the experiments of Rosenhek-Goldian et al. [163] was questioned by Jee et al. [164] who suggested that the "missing" dilatancy behaviour was actually observed experimentally on the same system several years earlier [165]. Thus there is clearly still much to learn regarding the stick-slip friction behaviour observed for strongly confined fluids, and NEMD simulations are likely to continue to play a central role in improving our understanding of this phenomenon.

Boundary slip
The no-slip boundary condition is the common assumption in continuum hydrodynamics which states that fluid adjacent to a sliding solid surface moves at the same velocity as the surface [166]. However, several experiments have shown that this is not always the case for fluids confined by smooth surfaces with oleophobic coatings (usually self-assembled monolayers or SAMs), as reviewed in 2005 by Neto et al. [167].
Boundary slip occurs when the fluid cohesion is relatively stronger than its adhesion to the sliding surfaces [105], meaning that the surfaces do not transfer sufficient momentum to fully shear the fluid [168]. This can lead to large reductions in friction compared to the Couette case [146]. Slip has been commonly observed in NEMD simulations of thin fluid films confined between atomically-smooth surfaces [57,[168][169][170][171]. Many of these studies have quantified the slip length, which is defined as an extrapolated distance into the wall where the tangential velocity component vanishes (see Fig. 3) [172]. There is significant disagreement between the slip length results from different experimental procedures, and several mechanisms for slip have been proposed, from enhanced liquid mobility at the surface [173,174] to bubble formation at the fluid-solid interface [175]. NEMD simulations offer a useful opportunity to test these mechanisms, and they have led to more accurate models for boundary slip behaviour in a range of systems [57,176].
NEMD simulations have also been used to study the change in slip length under a wide range of conditions. These simulations have shown that the slip length generally increases with increasing sliding velocity and pressure [168][169][170][171]. At very high sliding velocity, the slip length asymptotes towards a constant value [57]. An important observation is that realistic velocity-slip behaviour can only be captured in NEMD simulations when flexible, rather than rigid walls are employed [57]. NEMD simulations have shown that the slip length increases with increasing molecular chain length and degree of branching, in line with the viscosity [151]. Inflexible molecules, which can conform to the surface only with a greater expense of conformational energy generally exhibit more slip [177]. Conversely, the slip length decreases with increasing surface-fluid attraction strength and film thickness [169][170][171], and is virtually eliminated in the presence of atomic-scale surface roughness, which frustrates molecular layering [153,154]. This final observation is critical since most experimental surfaces will contain fractal roughness, even down to the atomic scale [178], which is likely to discourage boundary slip on conventional surfaces. Moreover, large-scale NEMD simulations of > 100 nm n-hexane films between atomically-smooth iron surfaces under realistic EHL conditions (F N = 1 GPa,  = 10 5 s −1 ) confirmed the suitability of a macroscopic no-slip boundary condition for this particular system [170]. NEMD simulations of n-alkanes (C 3 −C 10 ) confined between surfaces with patterned stick and slip regions by Savio et al. [179] showed that such heterogeneous surfaces can lead to cavitation under sufficient confinement.
Continuum models for fluid flow, such as the Reynolds equation [180], generally assume a no-slip boundary condition [181]. Slip lengths extracted from NEMD simulations have been fed into a modified Reynolds equation to study the effect of slip at one [182] or both [168] of the surfaces on flow behaviour, film thickness, and EHL friction. These results have shown that slip can lead to a significant reduction in friction and film thickness compared to the no-slip case, particularly at high shear rates.
Boundary slip has been directly observed in real tribological contacts under EHL conditions for a viscous polybutadiene lubricant between smooth glass surfaces with an oleophobic coating using [183]. Here the slip length was shown to increase with pressure and slip was also shown to substantially reduce friction under EHL conditions [183]. Slip has also been inferred from optical interferometry experiments on glass and even steel surfaces under EHL conditions for viscous polybutadiene and polyphenyl ether (5P4E) lubricants [184]. However, slip has not yet been observed experimentally for realistic lubricants confined between conventional surfaces and thus its significance in tribology remains uncertain [185].

Shear localisation
In the last few years, NEMD simulations of confined L−J liquids sheared at GPa-level pressures have revealed a range of nonequilibrium phase behaviour [186][187][188]. These simulations were designed to mimic the conditions found in the EHL regime, which is particularly important for machine components that roll and slide together, for example; rolling bearings, gears, constant velocity joints and cam/follower systems [41]. Although shear in these simulations was applied through the walls, the thickness of the films was above the ten molecular-layer limit in which fluid films experience enhanced viscosity as a direct result of confinement [103,139]. These systems and conditions could also be investigated through bulk NEMD simulations; however, in this case most thermostats enforce a linear (Couette) flow profile [50] (Section 1.3.1), which may not be the case in reality, as discussed below.
Above a certain pressure, the NEMD simulation [186][187][188] flow profiles transitioned from Couette-like to those indicative of shear localisation (or banding). Nonequilibrium phase diagrams were constructed for L−J fluids to map the changes in behaviour with pressure and shear rate (or sliding velocity). Figure 4 shows the different types of flow behaviour which have been observed in NEMD simulations of fluids under EHL conditions. Unlike the simulations described in Section 2.2.3, in these high pressure simulations, slip does not occur at the solid−liquid interface (PS in Fig. 4), but rather within the fluid itself. At low shear rates, shear was localised close to the walls, while the central part of the fluid formed a solid-like, unsheared region. This behaviour is usually referred to as "plug-slip" (PS) (PS in Fig. 4). At higher shear rates, the localisation of shear velocity gradient moved to the centre of the film, while the fluid near the walls became solid-like. This behaviour is typically known as "central localisation" (CL) (CL in Fig. 4). As the shear rate is increased in the CL region, the central liquid region widens until eventually a Couette profile is recovered [186][187][188]. Under the lowest shear rates studied, stick-slip friction behaviour was observed due to the commensurability of the walls and the solid-like fluid [188]. This mirrors the stick-slip behaviour [102] observed for more strongly confined films at ambient pressure (Section 2.2.2), supporting the postulate due to Robbins and Smith [162] that confinement-induced phase transitions observed in SFA experiments and highly confined NEMD simulations may be closely related to the pressure-induced transitions observed in macroscopic tribology experiments.
Slip within the fluid itself, rather than at the solid-fluid boundary, was first observed in NEMD simulations of an L−J fluid in 1990 by Thompson and Robbins [189]. In these simulations, only the first one or two layers of fluid were "locked" to the surface (due to strong wall interaction) and the rest of the fluid was sheared. This is distinct from the shear localisation behaviour observed in the high pressure NEMD simulations described above [186][187][188] where the "unsheared" region extends well beyond the direct influence of the wall atoms. Subsequent confined NEMD simulations of L−J fluids between one crystalline and one amorphous wall by Butler and Harrowell [190][191][192], L−J glasses (mixture of different sized L−J atoms to frustrate crystallization) by Varnik et al. [193], as well as L−J glasses and amorphous polymers by Rottler and Robbins [194] have shown various forms of shear localisation. Thus the onset of shear localisation at high pressure is probably related to localised phase changes within the lubricant (either vitrification or strong crystallisation), similar to those occurring due to confinement (Section 2.2.1) [162].
The discovery of such nonequilibrium phase behaviour by NEMD simulations of L−J fluids demonstrates the advantage of carrying out simulations using simple model systems. By doing so, it helps to reveal whether the phenomenon is "universal" to the liquid state and arises for all (or at least the vast majority of) liquids under certain conditions (which vary between molecules), or is specific to certain types of molecule [195].
In the field of tribology, shear localisation has only been directly observed experimentally in viscous model lubricants, such as polybutadiene and polyphenyl ethers (e.g., 5P4E) [196][197][198][199][200][201]. However, it has also been inferred from film thickness measurements in  | https://mc03.manuscriptcentral.com/friction realistic mineral oils [202,203] and there is currently a growing consensus that this phenomenon may also be important for more realistic lubricants under extreme EHL conditions [185]. NEMD simulations provide a useful complement to these experiments since they can more readily study the challenging conditions of interest.
Recently, extensive NEMD simulations of weakly confined L−J fluids under GPa pressure have revealed that the nonequilibrium phases correlate well with corresponding friction maps [188]. The characteristics of friction on the atomic scale can deviate quite significantly from the laws of classical friction. For example, when the applied conditions yield shear localisation, the friction coefficient of L−J fluids can decrease with increasing pressure and shear rate. Such behaviour is not commonly observed experimentally for lubricants, but is more similar to how traction fluids behave [204] above their critical shear stress [185]. The nonequilibrium phase diagrams and friction maps for L−J fluids were also found to be sensitive to the degree of the wall roughness on the atomic scale [188].
In similar high-pressure NEMD simulations of L−J glasses [187], crystallisation was suppressed and only liquid or weak CL phases were observed. In these systems, the friction coefficient increased linearly with log(shear rate), as has commonly been observed for realistic lubricants, such as PAO, under EHL conditions [41]. In this case, the friction-shear rate behaviour can be adequately described using stress-augmented thermal activation theory [205], according to the rheological models of Eyring [132]. An underlying assumption in current models of EHL friction is that, in the absence of thermal effects, or until a critical shear rate is reached, the rheological properties of the lubricant do not vary through its thickness. This is implicit in applying shear thinning models to predict friction in the EHL regime since they generally assume that the film is subject to planar Couette flow. Thus significant modifications of these models may be required where there are deviations from Couette flow, such as shear localisation [41].
Recently, a combination of tribological experiments and confined NEMD simulations has been used to investigate the effect of base fluid molecular structure on nonequilibrium phase behaviour and friction by Ewen et al. [29]. An extensive parameter study, including several lubricant and traction fluid molecules subjected to pressures (0.5−2.0 GPa) and shear rates (10 4 −10 10 s −1 ) typical of the EHL regime, revealed clear relationships between the friction and flow behaviour. There was good agreement between the friction-shear rate behaviour from tribology experiments and the NEMD simulations conducted at higher shear rates. Lubricants, which are flexible, broadly linear molecules, gave low friction coefficients that increased with shear rate and pressure in both the experiments and the simulations. Conversely, traction fluids, which are based on inflexible cycloaliphatic groups, gave high friction coefficients that were less sensitive to shear rate and pressure. The observed differences in friction behaviour was rationalised through the stronger shear localisation which was observed for the traction fluids (CL and PS) compared to the lubricants (Couette and CL) in the NEMD simulations [29]. Washizu et al. [206] also reported shear localisation (CL) in NEMD simulations on n-hexane confined between gold slabs at high pressure (0.1−8.0 GPa). At the transition from Couette flow to CL, they observed limiting shear stress behaviour, where the friction coefficient became insensitive to increasing pressure.

Summary
Many bulk NEMD simulation studies have accurately the Newtonian viscosity of realistic lubricant molecules as well as its temperature and pressure dependence. The accuracy of these simulations has reached a point where they can be used as a predictive tool to compare the rheology of different lubricant molecular structures. Such structure-property relationships can be extremely useful in lubricant design, and NEMD simulations provides the opportunity to probe a wide parameter space cheaply and efficiently. Bulk NEMD simulations have also given atomic-scale insights into lubricant shear thinning behaviour under extreme conditions and they are expected to play a central role in helping to improve macroscopic shear thinning models.
Confined NEMD simulations have also proved extremely useful in understanding the behaviour of thin lubricant films between solid surfaces. This includes building fundamental understanding of the density and viscosity inhomogeneity in confined films, as well as knowledge of important tribological phenomena under shear such as stick-slip, boundary slip, and shear localisation. Confined NEMD simulations are particularly useful since they mirror the layout and physics of real tribology experiments and require fewer assumptions with respect the flow profile and thermostatting procedure compared to bulk NEMD simulations. Thus, confined NEMD simulations are being increasingly used to study relatively thick lubricant films in tribology. Interestingly, there seem to be inherent similarities between the tribological phenomena observed in strongly confined lubricant films and thicker films subjected to high pressures.
Critical to any simulation technique is validation through direct comparison with experiment. As discussed above, this is becoming possible for NEMD simulations thanks to significant improvements in simulation algorithms and computer hardware, which have allowed larger systems and lower shear rates, of direct relevance to tribology experiments, to be investigated.

NEMD simulations of lubricant additives
A more recent area of research has been the application of MD to study the behaviour of lubricant additives. The function of many lubricant additives depends on chemical reactions, which are generally not considered in classical MD simulations. Despite this, classical MD simulations have given unique insights into the behaviour of a range of lubricant additives including detergents [207,208], dispersants [209,210], viscosity modifiers [91,92,211], anti-wear additives [212][213][214][215], and corrosion inhibitors [216][217][218]. The most widely studied class of lubricant additives with MD are friction modifiers [219], which are a wellsuited application of confined NEMD simulations.
The need for greater energy efficiency to reduce energy consumption has led to a shift toward lowerviscosity lubricants, which means that an increasing number of engineering components operate under boundary lubrication conditions [219]. As a result, lubricant additives that reduce friction and wear under boundary conditions (i.e., friction modifiers) are of increasing importance [219]. Several classes of friction modifier additives exist, the main ones being organic friction modifiers (OFMs), functionalised polymers, organo-molybdenum additives and dispersed nanoparticles [219]. An overview of the experimental advances on friction modifier additives is given by in 2015 by Spikes [219]. Sections 3.1 and 3.2 will focus on classical MD simulations of OFMs and nanoparticle friction modifiers respectively, and a brief discussion of "reactive" MD simulations of other additives is given in Section 4.1. In general, the simulation methods described in Section 1.3 (e.g., thermostats, barostats, force-fields) are also applicable to study these more complex additive systems.

Organic friction modifiers
OFMs are amphiphilic surfactant molecules that contain nonpolar hydrocarbon tail groups attached to polar head groups. Their main friction reduction mechanism involves the adsorption of the polar head group to metal, ceramic, or carbon-based surfaces, with strong, cumulative van der Waals forces between proximal nonpolar tails leading to the formation of incompressible monolayers that prevent contact between solid surfaces to reduce adhesion and friction [220]. Many different head and tail groups have been tested in the literature, as discussed elsewhere [219].
Classical NEMD simulations can be used to simultaneously probe the nanoscale structure and friction of OFM films, making them a valuable complement to experiments. Adsorption of amphiphiles from solution to solid surfaces using realistic additive concentrations is a relatively slow process, and close-packed monolayers can take hours to form in the absence of shear [221][222][223]. As a result, MD simulation of the adsorption process using realistic concentrations and all-atom force-fields are currently unattainable [224]. However, extensive experimental analysis has shown that OFMs do eventually form close-packed monolayers on a range of contact surfaces, and that the formation of these films is accelerated by the shearing process [221][222][223]. Therefore, MD simulations are generally constructed with pre-formed OFM monolayers by placing additive molecule head groups close to solid surfaces at the beginning of the simulation. This is similar to the Langmuir-Blodgett experimental approach to quickly form surfactant films with a controllable packing density (or coverage) on solid surfaces [225]. MD simulations generally consist of OFM monolayers adsorbed on two solid surfaces, sometimes also with 366 Friction 6(4): 349-386 (2018) | https://mc03.manuscriptcentral.com/friction a layer of confined lubricant in between. EMD simulations have been performed to study the structure of OFM films on solid surfaces. Confined NEMD simulations have also been performed ( Fig. 1(a)) to study the structure, flow, and friction behaviour of the films. The results of these OFM simulations are discussed below.

Film structure
Moller et al. [226] were the first to use EMD simulations to study the structure of OFM films (stearic acid) on solid (graphite) surfaces in 1991. They varied the size of the head group to study its effect on the tilt angle from the surface normal. Their simulations predicted that, under ambient conditions, OFM molecules were normal to the surface at high surface coverage (head group areas below 21 Å 2 ), but became tilted from the normal at lower surface coverage. They suggested that titling maximises the van der Waals attraction between the chains by increasing the packing efficiency [227]. The tilting transition was attributed to the packing of the hydrogen atoms belonging to methylene groups on neighbouring molecules, meaning that the correct behaviour could only be accurately reproduced using all-atom (AA) and not united-atom (UA) force-fields [226]. Recent NEMD simulations by Ewen et al. [10] showed that the use of AA force-fields is also critical to give accurate flow and particularly friction behaviour for OFM films under shear.
Davidson et al. [228] performed EMD simulations to help elucidate structure-activity relationships of mono-, di-, and tri-oleyl glyceride OFMs. Mono-glyceride molecules with two hydroxyl groups and one alkenyl chain packed more efficiently than di-glycerides which contain one hydroxyl group and two alkenyl chains. This was because the mono-glycerides occupied approximately half the area (22 Å 2 ) needed by di-glycerides (46 Å 2 ) and formed significantly more intermolecular hydrogen bonds. The density profiles showed significant mixing of the hydrophobic tail groups with the non-polar solvent and the distribution of tail group torsion angles suggested the gylcerides formed liquidlike rather than solid-like films at room temperature. Experimentally measured friction coefficients of equimolar glyceride solutions showed that the efficacy as friction modifiers varied in the order mono-, diand then tri-oleyl glyceride, which was consistent with the EMD film formation results.
Greenfield and Ohtani [229] used EMD simulations to investigate the orientation of monolayers of model OFMs (head group area 17 Å 2 ) confined between atomically-smooth surfaces. The normal load, fluid layer thickness, and additive surface concentration dependencies agreed with those measured experimentally using the SFA [230]. Greenfield and Ohtani also performed NEMD simulations of the same systems and found that the OFM film structure remained virtually unchanged under a wide range of normal pressures (3.2−22 MPa) and sliding velocities (0.5− 7.5 m·s −1 ) [231].
Doig et al. [232] also used NEMD simulations to gain insights into sum-frequency spectroscopy (SFS) and polarised neutron reflectometry (PNR) experiments [233] of the structure of hexadecylamine films adsorbed on iron oxide surfaces in dodecane and hexadecane. At the highest surface coverage investigated in the NEMD simulations (1.8 nm −2 ), the film thickness was about 15-20 Å, and the molecules had an average tilt angle from the surface normal of 40°-50°, in agreement with experiment. The in-layer ordering of the hexadecylamine head group atoms was found to be dictated by the crystalline structure of the iron oxide surface but, at the coverages studied, this influence quickly diminished further along the tail group.
Polar OFM molecules generally form dimers or reverse micelles [219] in non-polar lubricants since they are more thermodynamically stable. For some OFMs, such as glycerides, it remains unclear whether close-packed monolayers form on the surface, as experimentally observed for fatty acids [221][222][223], or whether they have a different friction reduction mechanism. Bradley-Shaw et al. [234] used NEMD simulations to study the stability of reverse micelles, formed by glycerol mono-oleate (GMO) OFMs in organic solvents [235], confined and sheared between mica surfaces. The GMO reverse micelles were found to be metastable in the absence of shear for all of the systems studied. When shear was applied, the reverse micelles remained intact and adsorbed onto one of the surfaces in n-heptane, thus becoming surface micelles. In dry toluene, the reverse micelles broke apart under shear, forming disordered monolayers on the surface. However, in toluene with water present, they survived and became surface micelles.

Friction behaviour
In 1993, Glosli et al. [236] became the first to use NEMD to investigate the friction of OFM-like systems. They studied the friction between two close-packed monolayers (head group area 17 Å 2 ) bound at their ends to two rigid substrates (similar to OFM molecules absorbed on solid surfaces). Depending on the interfacial interaction strength between the two layers, energy dissipation either occurred by a discontinuous "plucking" mechanism (strong interaction) or a continuous "viscous" mechanism (weak interaction). Friction during the plucking mechanism followed the simple stress-augmented thermal activation model proposed by Briscoe and Evans [225] from experimental friction data for Langmuir-Blodgett films. Friction during the viscous mechanism was generally much lower but increased at the rotator transition temperature of the monolayers [236] (T at which molecules become rotationally ordered in a herringbone-like pattern).
Kong et al. [237] performed NEMD simulations of adsorbed monolayers of dialkylammonium chloride surfactants between charged solid surfaces. They found that the friction coefficient decreased with increasing normal force (30 to 450 MPa), decreased with increasing surfactant coverage (77-50 Å 2 ), and increased with increasing sliding velocity (1-100 m·s −1 ). Moreover, the friction between the OFM films correlated well with the amount of layer overlap or interpenetration, as defined by the common area under the chain density profiles. Each dialkylammonium chloride surfactant molecule contains two alkyl chains and, in a subsequent NEMD study [238], the effect of monolayers with asymmetrical (C 10 C 18 ) and symmetrical (C 18 C 18 and C 10 C 10 ) tail groups was investigated. The friction between the layers of asymmetrical surfactants was greater than that between layers of symmetrical surfactants, suggesting that films with uneven heights were less effective in reducing friction.
Eder et al. [239] performed NEMD simulations of several OFMs on iron surfaces, varying the OFM type (stearic acid, oleic acid, methyl stearate), OFM surface coverage (2.54−4.57 nm −2 ), the pressure (0.07−1.60 GPa), and the surface roughness. The resulting load versus friction behaviour was then analysed to study these parameters affected the friction coefficient and the extrapolated friction force at zero load (Derjaguin offset). A smooth-particle-based evaluation method [240] was used to visualize the sliding interface between the two OFM layers. This allowed the definition and calculation of a dimensionless normalized sliding resistance area, which was then related to the Derjaguin offset. The friction coefficient and Derjaguin offset were found to generally increase with increasing surface roughness, and decrease with increasing OFM coverage. The Derjaguin offset was eliminated by the addition of a thin n-hexadecane film between the OFM layers.
Doig et al. [241] used NEMD simulations to examine the structure and friction of stearic and oleic acid OFM films adsorbed on iron oxide surfaces and lubricated by squalane under hydrodynamic conditions (0.1 GPa). At high surface coverage (2.59 nm −2 ), the measured properties of stearic acid and oleic acid films were very similar. At low (0.58 nm −2 ) and intermediate (1.30 nm −2 ) surface coverages, the presence of Z-unsaturation (in oleic acid) resulted in less penetration of lubricant in the surfactant film and less layering of the lubricant near the film. The friction coefficient increased linearly with log(shear rate) within the hydrodynamic lubrication regime, consistent with a stress-augmented thermal activation model [225] and macroscopic friction experiments [242]. Lubricant penetration and layering were found to affect the hydrodynamic friction coefficient. Oleic acid friction depended more weakly on surface coverage, while stearic acid allowed more lubricant penetration, and thus friction increased significantly with decreasing surface coverage.
Recently, Ewen et al. [243] used NEMD simulations to examine the atomistic structure and friction properties of commercially relevant organic friction modifier (OFM) monolayers adsorbed on iron oxide surfaces and lubricated by a thin, separating layer of n-hexadecane (Fig. 5). The hexadecane thickness was determined at the studied pressure (0.5 GPa) in preliminary LVMD squeeze-out simulations. Carboxylic acid, amide, and glyceride OFMs, with saturated and Z-unsaturated hydrocarbon tail groups, were simulated at various surface coverages (1.52−4.56 nm −2 ) and sliding velocities (0.5−20 m·s −1 ). At low and medium coverage, the OFMs formed liquid-like and amorphous monolayers, respectively, which were | https://mc03.manuscriptcentral.com/friction significantly interdigitated with the hexadecane lubricant, resulting in relatively high friction coefficients. At high coverage, solid-like monolayers were formed which, during sliding, resulted in slip planes between well-defined OFM and hexadecane layers, giving a significant reduction in the friction coefficient [243]. In agreement with experiment [242], OFMs with glyceride head groups gave significantly lower friction coefficients than amide and particularly carboxylic acid head groups. For all of the OFMs and coverages simulated, the friction coefficient was found to increase linearly with log(sliding velocity); however, the gradient of this increase depended on the coverage [243]. The friction-sliding velocity behaviour from the NEMD simulations was consistent with a stress-augmented thermal activation model [225] and macroscopic boundary friction experiments [242].
Ewen et al. [244] also used NEMD simulations to examine the structure and friction of OFM (stearic acid) films adsorbed on iron surfaces with random nanoscale roughness (0.2−0.8 nm root-mean-squared). The OFMs were strongly adsorbed, were not squeezed out, and prevented direct contact of asperities under the high pressures simulated (0.5−2.0 GPa). Increased coverage generally resulted in lower lateral (friction) forces due to reductions in both the friction coefficient and Derjaguin offset. Rougher surfaces led to more liquid-like, disordered films, but the friction coefficient and Derjaguin offset were only slightly increased. These results suggested that OFM films are almost as effective on contact surfaces with nanoscale roughness as those which are atomicallysmooth. Further NEMD simulations of OFM films using rough surfaces, or perhaps GCMD simulations [63,64], would be useful to improve understanding of their resilience to high local pressures experienced in the boundary lubrication regime.
Although the focus of this review is on lubricants for macroscale engineering components such as engines and transmissions, NEMD simulations of confined films have also been performed to study friction modifier additives for microelectromechanical systems (MEMS) and nanoelectromechanical systems (NEMS). In this context, NEMD simulations have been used to investigate the friction behaviour of monolayers of carboxylic/perfluorocarboxylic acids [245][246][247] and alkylsilanes [248][249][250][251] on silica surfaces. The tribological behaviours of model self-assembled monolayers (SAMs) such as thiols on gold surfaces have also been widely studied with NEMD [252][253][254]. These simulations have revealed the effect of head group and tail group of the nanoscale structure and friction behaviour of the monolayers.

Nanoparticle friction modifiers
In recent years, there has been growing interest in the possibility of using colloidal solid particles in the size range 1−500 nm as friction modifier and anti-wear additives [219]. Many types of nanoparticle additives have been tested, from various carbon allotropes to metals and inorganic salts [219]. The development of fullerenes [255] and inorganic fullerenes has generated particular interest since their graphitic-like structures bear some resemblance to lamellar materials which are known to effectively reduce friction [219].
Despite this interest, there is still considerable uncertainty as to the friction reduction mechanisms of nanoparticle additives [219]. When dispersed in a lubricant, nanoparticles can modify the lubricant viscosity and potentially also the flow behaviour [256]. Under boundary lubrication conditions, when surfaces come into direct contact, several mechanisms have been proposed for friction and wear by nanoparticles including [257] sliding (a), rolling (b), exfoliation (c), polishing (d), and mending (e), as shown in Fig. 6.
If nanoparticles are able to maintain the separation of surface asperities under boundary conditions they could give large reductions in friction and wear [257]. This requires nanoparticles to be sufficiently hard to not be plastically deformed under high contact pressures. However, if the nanoparticles are excessively hard, they may indent into the surfaces to such an extent that they no longer separate the asperities. Nanoparticles that are able to prevent asperity contact may roll ( Fig. 6(a)) or slide ( Fig. 6(b)) between the surfaces in order to reduce friction. Experimental studies have suggested that unlike macroscopic systems, there is little (if any) benefit of the rolling action for nanoparticles compared to sliding [258,259]. Hard nanoparticles may also reduce friction through polishing ( Fig. 6(d)), whereby asperities are smoothed by abrasion by the nanoparticles to reduce roughness and shift the system from boundary, towards mixed lubrication conditions [260]. While this may initially provide a large benefit, it will probably diminish in the presence of other running-in and smoothing processes [261]. It has also been suggested that hard nanoparticles may also "mend" the surface by filling the gaps between asperities on rough surfaces (Fig. 6(e)) [260]. For soft or layered nanoparticles, there is also the possibility that friction may be reduced through exfoliation (Fig. 6(c)), whereby the nanoparticles break down into a solid lubricating film, which prevents contact between the surfaces and reduces friction by forming a low shear strength layer [262].
The prevalence of each of these mechanisms is expected to be highly situation-specific [219], and depend on variables including the material, size, shape, hardness, and coverage of the nanoparticle as well as the temperature, pressure, and sliding velocity [257]. For example, in-situ high resolution scanning/ transmission electron microscope (HRSEM/HRTEM) [262,263] images of inorganic fullerene (MoS 2 ) nanoparticles during sliding have suggested that their friction reduction mechanism consists of a mixture of sliding/rolling at low pressure and exfoliation to lamellar sheets of MoS 2 at high pressure. Exfoliation of hollow core MoS 2 fullerenes was also observed in MD simulations by Lahouij et al. [263] under high uniaxial loads.
The simulation of nanoparticles dispersed in a lubricant and confined between solid surfaces is a complex problem for study by MD simulations. The nanoparticles tested in experiments (generally > 10 nm diameter) are relatively large for atomistic MD, and thus only a few small nanoparticles are usually  | https://mc03.manuscriptcentral.com/friction included in the system. In addition, although accurate force-field parameters are available for each component in isolation, they are generally not transferable between the different components. Moreover, the addition of a lubricant between atomically-smooth surfaces with PBCs in the lateral directions means that the lubricant supports the load. This means that only the viscosity/flow modification of the nanoparticles, rather than their boundary friction reduction mechanisms (Fig. 6) can be probed [257]. Therefore, most NEMD studies have studied nanoparticles confined between solid surfaces in the absence of a lubricant, even though they are designed as additives (Section 3.2.1), although some studies have also been performed with nanoparticles dispersed in a lubricant (Section 3.2.2).

Nanoparticles as dry lubricants
There have been a considerable number of NEMD simulations to investigate nanoparticle behaviour in confined systems in the absence of a liquid lubricant. For example, Joly-Pottuz et al. [264] conducted friction and wear experiments and NEMD simulations of carbon nano-onions (2 nm diameter) between sliding diamond-like carbon (DLC) surfaces at high contact pressure (1−5 GPa). Both the experiments and simulations suggest that the carbon nano-onions remained intact under compression and sliding and did not exfoliate to graphitic planes [264], unlike MoS 2 fullerenes [263]. Further simulations by Bucholz et al. [265] indicated that the ability of the nano-onions to roll is inhibited both by increased contact pressure and the presence of a diamond core within the nanoparticles. The transition from rolling to sliding behaviour was accompanied by a significant increase in the friction coefficient due to the formation of interfacial bonds. Bucholz et al. [266] also investigated the behaviour of amorphous carbon nanoparticles (2−5 nm) between hydrogen-terminated diamond surfaces with NEMD. The simulations showed that during sliding, hydrogen passivated the unsaturated carbon atoms near the surface of the nanoparticles, which prevented interfacial bond formation and allowed them to roll within the interface. It was demonstrated that amorphous carbon nanoparticles will only provide good tribological performance when sufficient chemical passivation of the surfaces and nanoparticles is maintained.
Hu et al. [267] investigated the effect of soft metallic copper nanoparticles between relatively hard iron surfaces using NEMD simulations. The results suggested that the nanoparticles deformed to generate a low shear strength film between the contact surfaces, which accommodated the velocity gradient through plastic deformation. The nanoparticles were most effective in reducing friction at low loads and sliding velocities (0.01 GPa, 10 m·s −1 ), but still provided a small benefit under more extreme conditions (0.5 GPa, 500 m·s −1 ). Hu et al. [268] also used NEMD simulations to investigate the tribology behaviour of hard diamond and silicon dioxide (SiO 2 ) nanoparticles. At low velocity and low load (10 m·s −1 , 0.5 GPa), the nanoparticles separated the two contact surfaces from each other and significantly reduced friction. At high velocity and high load (500 m·s −1 , 1.0 GPa), they either plastically deformed (SiO 2 ) or indented into the contact surfaces (diamond), leading to smaller reductions in friction.
Eder et al. [269] simulated the abrasion process of an atomically rough iron surface with multiple hard, abrasive particles. By monitoring the nanoscopic wear depth over time, the authors showed that the Barwell macroscopic wear law [270] holds even at the atomic scale. They found that in this multiasperity contact system, the Bowden-Tabor term [271], which describes the friction force as a function of the real nanoscopic contact area, can predict the friction even when wear is involved.
Ewen et al. [257] used NEMD simulations to examine the friction and wear reducing mechanisms of carbon nanoparticles. Specifically, the friction and wear behaviour of carbon nanodiamonds (CNDs) and carbon nano-onions (CNOs) with a diameter of approximately 3 nm, confined between iron slabs was probed at a range of coverages, pressures, and sliding velocities. At high coverage and low pressure, the nanoparticles did not indent into the iron slabs during sliding, leading to zero wear and very low friction. At low coverage and high pressure, the nanoparticles indented into, and ploughed through the slabs during sliding, leading to atomic-scale wear and a much higher friction coefficient. This contribution to the friction coefficient was well predicted by an expression developed for macroscopic indentation by Bowden and Tabor [271]. Even at the highest pressure and lowest coverage simulated, both types of nanoparticle were able to maintain separation of the opposing slabs and reduced friction by approximately 75% compared to when no nanoparticles were present, in agreement with experimental observations [272]. CNO nanoparticles showed a lower indentation (wear) depth and lower friction coefficients than CND nanoparticles at the same coverage and pressure, making them more attractive friction modifier additives.

Nanoparticles as lubricant additives
Although NEMD simulations of nanoparticles can provide interesting insights into their friction and wear reduction under boundary lubrication conditions, it is clearly also important to study them in the presence of a liquid lubricant. NEMD simulations have been performed on atomically-smooth surfaces to study their effect on friction [273] and flow [256] under hydrodynamic lubrication conditions. Compressive LVMD simulations have also been used to study the load carrying capacity of copper nanoparticles in an n-octane lubricant [274].
Recently, NEMD simulations have been used to complement macroscopic wear experiments to study the effect of titanium dioxide nanoparticle lubricant additives in a PAO lubricant [275]. Results from both the experiments and simulations suggested that there was an optimal particle size that will minimize friction and wear for a given surface roughness [275]. The results support the "mending" mechanism for the titanium dioxide nanoparticle additives (Fig. 6(e)). In this context, particles that are smaller than the characteristic roughness of the surfaces are most likely to perform this function. Given the considerable interest in nanoparticle additives [219], we expect to see further NEMD studies of nanoparticles in the presence of a lubricant on rough surfaces.

Summary
NEMD simulations of a range of friction modifier additives have given atomic-scale insights into the films they form on solid surfaces and the way in which they alter the flow and friction behaviour under a wide range of conditions. Due to the long timescales involved, MD simulation of film formation from solution remains a challenge. However, by assuming film formation on the surface, they have reproduced important experimental observations, and helped to explain these from atomic-scale behaviour. They have highlighted important differences between various additive molecular structures and coverages. To more faithfully reproduce boundary lubrication conditions, NEMD simulations of additives between realistic rough surfaces or using GCMD, which both allow squeeze out from high pressure regions, are expected to provide interesting results. An important aspect of lubricant additive behaviour which cannot be studied using conventional MD simulations is chemical reactivity. However, newly developed methods have begun to allow such investigations using molecular simulations and these are discussed in Section 4.1.

Beyond classical NEMD: Current trends and future perspectives
NEMD simulations have made remarkable advances in the past few decades both in terms of methodology development and application-driven research. Despite this, some complex tribology problems (e.g., tribochemistry, triboemmission, adhesion, plasticity) are influenced by physical/chemical processes that occur at time and length scales which are not accessible through classical NEMD simulations. Our understanding of such problems would significantly benefit from models that span multiple spatial and/or temporal scales with different levels of complexity. Such multiscale/ multiphysics models use information from several levels (that investigate a phenomenon over a specific temporal and spatial window) to calculate material properties or model system behaviour. The following main levels are usually distinguished (Fig. 7): quantum mechanical models (electronic scale), MD models (atomic scale), and continuum models (macroscopic scale). Figure 7 provides a map to illustrate the different modelling techniques available to describe tribological phenomena at different scales. The focus of this review is on molecular simulations, and classical MD is the centrepiece of the map. Techniques which include electronic effects, e.g., density functional theory (DFT),  can be used to study processes at smaller scales, while continuum methods, such as finite element methods (FEM) and computational fluid dynamics (CFD), can describe processes at larger scales. A full review of all these methodologies is certainly outside the focus of this review; however, here we will highlight successful multiscale/multiphysics approaches that already exist, and suggest how the furthering of such methodologies could play a fundamental role in addressing unsolved issues in the field of tribology. One area where multiscale/multiphysics modelling has been successfully applied is integrated computational materials engineering (ICME), which allows the prediction of material properties or system behaviour based on knowledge of the process-structure-property relationships [276]. These same attributes make these models interesting in tribology, particularly in the design of synergetic surfaces and lubricants [277]. Breakthroughs are thus expected from the future development of versatile and efficient multiscale/ multiphysics tools designed specifically for tribology problems.
Critical to these advances will be to link MD with simulations performed at smaller and larger scales (see Fig. 7). In doing so, modelling strategies can be designed to minimise the computational resources required while maximising the insight provided by simulations performed either simultaneously or sequentially at different scales. In the next two subsections, an outlook is provided which describes the most promising multiscale/multiphysics techniques currently available to link MD simulations to quantum mechanics (Section 4.1) and continuum models (4.2) for tribology problems.

Linking MD to smaller scales
The function of many lubricant additives depends on chemical reactivity both in solution and at solid surfaces. Due to their large computational expense, conventional quantum mechanical calculations are static and include only a few molecules. Modelling chemical reactivity is beyond the capabilities of classical MD simulations; however, techniques which allow this to be studied have been developed and applied to study lubricant additives. One promising technique is the use of "reactive" force-fields, such as ReaxFF [278] which can approximate chemical reactivity in classical MD simulations through so-called bond order terms [279]. Liang et al. [280] gave a general overview of the different types of reactive force-field available in 2013. A major drawback of all reactive force-fields is the limited availability of reliable parameter sets to study the materials, conditions, and properties of interest. The generation of high-quality parameter sets is a very difficult and time consuming task due to the huge number of parameters which need to be fitted [281]. Despite this, ReaxFF parameters have been developed which can be used to investigate some important lubricant and additive molecules, such as alkanes and alcohols [282,283], phosphoric acids [284], organosulphur compounds [285], and molybdenum disulphide [286]. Using such force-fields, relatively large molecular systems (nm) can be compressed and sheared [284] in a similar manner to conventional confined NEMD simulations ( Fig. 1(a)), and the consequent chemical reactions monitored over reasonable simulations time (ns).
Alternatively, the reactivity of lubricant additives can be explicitly modelled using first principles or ab initio MD techniques. Examples of this include Car-Parrinello MD (CP-MD) [287][288][289][290][291] and tight-binding |www.Springer.com/journal/40544 | Friction http://friction.tsinghuajournals.com quantum chemical MD (TB-QCMD) [214,[292][293][294][295], which have given important insights into the reactivity of common friction modifier (molybdenum disulphide, molybdenum dithiocarbamate) and anti-wear (zinc dialkyldithiophosphate) additives. Ab initio NEMD simulations have also been performed on additive systems subjected to confinement and shear. For example, TB-QCMD has been used to investigate OFM adsorption on iron oxide surfaces under sliding conditions [296][297][298]. The main advantage of ab initio MD simulations over MD with reactive force-fields is that they do not require extensive force-field parameterisation in order to ensure accurate results. However, they increase the computational expense by more than an order of magnitude compared to classical MD simulations with reactive force-fields, meaning that the system sizes and timescales accessible are more limited [299]. One potential solution for confined NEMD simulations, where the "reactive" region is well-defined, is spatial coupling of classical MD for the "non-reactive" region to ab initio MD for the "reactive" region [300]. This technique is commonly known as quantum mechanics/molecular mechanics (QM/MM), and is widely used in materials science [301,302] and the natural sciences [303].

Linking MD to larger scales
Linking MD to larger scales is also of considerable interest to tribology researchers. As shown in Fig. 7, the next scale up from atomistic MD is coarse-grained MD (CG-MD), which involves grouping of neighbouring atoms together to form "beads" in order to reduce the number of interaction sites. It is commonly used to simulate large polymers, where monomers or even a combination of monomers can be represented as a single CG bead. CG-MD has not been extensively applied in tribology; however, it has been used to study friction in polymer brushes [304,305], model friction modifier additives [306] and ionic liquids [68]. The use of CG-MD is likely to remain rare in tribology, since the importance of atomistic detail has been highlighted in accurately reproducing the viscous behaviour of lubricants and additives [10]. Possible applications of CG-MD in tribology are to study the behaviour of large polymer additives [307] or relatively slow processes such as surfactant self-assembly at solid surfaces, whose timescales are inaccessible to study using all-atom force-fields [224].
Most tribological processes are governed by physical phenomena localised in an interfacial layer which controls the macroscopic behaviour of the system. With this is mind, an attractive possibility is to capture the interfacial laws of friction, heat transfer, and other relevant phenomena at small, accurate scales and pass this information to continuum models to study macroscopic tribological systems. Choosing how to pass information across the scales is not always straightforward; some issues concern what assumptions can be made and the accuracy required in linking the two levels. Two options are usually available, depending whether it is possible to establish scale-separation between the atomic level and the continuum. If such separation exists, then a hierarchical (or sequential) multiscale model is acceptable, otherwise a hybrid (or concurrent) scheme is required. Both of these MD-continuum coupling methods have already been applied to a variety of problems in tribology. Hierarchical coupling involves using MD simulations to gain information about the boundary conditions and constitutive laws to input into separate continuum models [308]. This method is relatively straightforward to implement, but requires careful planning to establish which regions of phase space can be reliably probed with MD to provide useful inputs for the continuum models. For example, Martini et al. [160] and Savio et al. [168,182] respectively used local viscosity and slip information obtained using MD simulations to improve continuum predictions of film thickness and friction for strongly confined systems. These are examples of phenomena that can be studied at a single relevant scale and integrated at a larger scale in a hierarchical manner.
Alternatively, the MD and continuum models can be run concurrently on different spatial regions; this is the preferred method when atomic interactions dictate the observed macroscopic tribological behaviour but are themselves affected by the whole macroscopic domain. Modelling surfaces with realistic levels of roughness lubricated by a fluid is a case in point. Such a problem is still beyond the scales accessible to full MD simulation [269]. Alternatively, atomistic details of the surface topography and fluid close to 374 Friction 6(4): 349-386 (2018) | https://mc03.manuscriptcentral.com/friction the walls can be explicitly modelled and intimately linked to the macroscopic response of the fluid further from the walls. In this case MD-continuum coupling strategies involve the transfer of information between the MD and continuum regions, and particular care must be taken when the two descriptions merge [309,310]. A number of schemes exist for this [311][312][313][314], with various procedures enabling hybrid coupling [315].
Models dealing with concurrent coupling are generally far more complex than their hierarchical alternatives, but are of high importance in understanding the physics of certain phenomena in contact interfaces, such as the effect of surface roughness, surface coatings and additive layers on fluid flow. One major limitation of concurrent coupling techniques is that the simulation time step is limited to the shortest timescale which needs to be modelled at the smallest scale. Despite this, considerable performance enhancements can be achieved compared to modelling the entire system at the atomic level [313]. The future success of concurrent coupling schemes will depend on the development of efficient coupling libraries and algorithms that facilitate hybrid coupling simulations and support their implementation on different HPC platforms (see, e.g., a recent software made available by the Imperial College Tribology Group [316]).

Summary
For the study of reactive lubricant additives, conventional MD simulations are deficient in their inability to model chemical reactivity. Moreover, for computational reasons, conventional MD simulations are limited to the study of relatively small systems and timescales which continues to limit their ability to adequately model situations of specific engineering significance. Moving forwards, we expect classical MD simulations to be supplemented by coupling techniques to smaller scales, to study chemical reactivity, as well as larger scales, to extend the accessible time and length scales. Tribology is an ideal subject for such multiscale/multiphysics models since many macroscopic tribological phenomena have been shown to driven by atomic-scale processes, often localised at an interface. Such tools will shed further light on these processes and could even transform lubricant formulation by identifying process-structure-property relationships as it has in other fields of research.

Conclusions
Molecular dynamics is, at its heart, a relatively simple simulation method, yet it has become one of the most powerful tools in science and engineering to understand the atomic-scale behaviour of fluids. Historically, it has played a central role in corroborating theories of the liquid state, but it has now become capable of directly evaluating physical properties in industrially-important systems. Moving forwards, we expect NEMD simulations to play an increasingly important role in building a fundamental or "bottomup" understanding of tribological phenomena from the atomic scale. Moreover, it is reaching the point at which it can be a truly useful tool for lubricant design. Key to this process will be direct comparison of NEMD to experiment, something which is beginning to become possible thanks to improvements in simulation algorithms, classical force-fields, and computer hardware, as well as experimental advances.
Many bulk NEMD studies have predicted the Newtonian viscosity of realistic lubricant molecules, as well as its temperature and pressure dependence. The accuracy of these simulations has progressed to a point where they can be used as a predictive tool to compare different lubricant molecular structures. Such structure-property relationships can be extremely useful in lubricant design. Bulk NEMD simulations have also given unique insights into lubricant shear thinning behaviour under extreme conditions and they are expected to play a key role in helping to improve macroscopic shear thinning models.
Confined NEMD simulations have revolutionised our fundamental understanding of the behaviour of very thin lubricant films between solid surfaces. This includes the density and viscosity inhomogeneities in confined films, as well as important tribological phenomena such as stick-slip and boundary slip. It is also being increasingly employed to study shear localisation behaviour in thicker films subjected to high pressures. There seem to be parallels between the behaviour of lubricants under extreme confinement and at thicker films at high pressure which certainly seem worthy of further study.
The atomic-scale behaviour of lubricant additives has also been successfully investigated with NEMD. There has been a particular focus on organic friction modifier and nanoparticle additives, in part because simulation of these additives does not necessarily need to account for chemical reactivity. These simulations have highlighted important differences between various additive molecular structures and coverages, which have direct relevance to lubricant formulation. The tribochemical behaviour of friction modifier and anti-wear additives under pressure and shear has also been given by NEMD simulations using reactive force-fields as well as ab initio MD techniques. Even though accessible system sizes and timescales of these simulations are more limited than for conventional MD, fundamental behaviour of reactive lubricant additives can be explored, potentially enabling more effective additives to be designed.
Finally, it should be noted that, while molecular dynamics can now simulate systems of millions or even billions of individual atoms, MD is unlikely ever to be suitable for explicitly modelling macro-scale engineering components (≈ 10 24 atoms). However, by coupling MD to continuum models, MD can be used just to capture key processes which drive the overall behaviour of macroscopic tribological systems.