The Importance of Reference Frame for Pressure at the Liquid-Vapour Interface

The local pressure tensor is non-unique, a fact which has generated confusion and debate in the seventy years since the seminal work by Irving Kirkwood. This non-uniqueness is normally attributed to the interaction path between molecules, especially in the interfacial-science community. In this work we reframe this discussion of non-uniqueness in terms of the location, or reference frame, used to measure the pressure. By using a general mathematical description of the liquid-vapour interface, we obtain a reference frame that moves with the interface through time, providing a new insight into the pressure. We compare this instantaneous moving reference frame with the fixed Eulerian one. Through this process, we show the requirement that normal pressure balance at the moving surface is satisfied by surface fluxes, however an additional corrective term based on surface curvature is required for the average pressure in a volume. We make the case that a focus on the path of integration is the cause of much of the confusion in the literature. Using an explicit reference frame with a more general derivation of pressure clarifies some of the issues of uniqueness in the pressure tensor and provides a pressure tensor which is defined at any instant in time and valid away from thermodynamic equilibrium.


Introduction
Since the pioneering work of Irving and Kirkwood [1] we have had a firm theoretical foundation for the pressure tensor in statistical mechanics. In practice, however, there remains an ongoing debate about what form of pressure tensor is the right one to use in a molecular dynamics (MD) simulation. Irving and Kirkwood [1] present the mathematics of a point in space, a direct consequence of using continuum quantities in a discrete system. Mathematically this gives us a Dirac delta function, which is an infinitesimal point and so must be adapted for use in a molecular dynamics simulation, expanding it over a finite area in order to collect measurements. Different communities of researchers have done this in different ways. The simplest approach is to set the delta functions to one and use the tensor version of the Virial [2] at a local point in space to get the pressure. This is the IK1 approximation, so called because it is a first order term in the full Irving Kirkwood expressions [3]. In the solid mechanics community, the delta function is sometimes replaced by defining a non-infinite function with a finite width and compact support as in [4][5][6][7]. These solid mechanics stress tensors, where stress is simply the negative of the pressure tensor, are largely in the form of edward.smith@brunel.ac.uk averages inside volumes, so have become known as the volume average (VA) expressions [8]. For the non-equilibrium molecular dynamics (NEMD) community, an approach treating the functions in Fourier space [3,9,10] has been followed to yield mathematically tractable forms. The most well-known form is the method of planes (MOP) [11] and its extensions [10,12]. Both the MOP and VA forms of pressure can be shown to be equivalent in the limit of small volumes [13], both demonstrating a pressure gradient of zero normal to a solid-liquid interface, a requirement for mechanical equilibrium. This is worth noting as the IK1 form of pressure fails to satisfy this requirement for mechanical equilibrium [14]. In addition to the VA, IK1 and MOP pressures from the solid mechanics and NEMD communities, researchers interested in the interface between two fluids, in particular between a liquid and a vapour, tend to work with pressure obtained relative to that interface. There are two common forms obtained by integrating relative to the interface area itself, called the Irving Kirkwood [15] and Harasima [16] contours. The difference between the Harasima and Irving Kirkwood forms is attributed to the different interaction paths, or interaction contours, assumed between the molecules [17]. This is generalised in the work of Schofield and Henderson [18] who show that the path of interaction between two particles, which defines their interaction force, can occurs in an infinite number of ways. This is explored for various contours in recent work [23]. Far from a philosophical point, this non-uniqueness has profound implications for measurements in MD systems. Published results were recently called into question due to the non-uniqueness of pressure [19,20]. In addition, use of Harasima and Irving Kirkwood surface forms led to the conclusion that Marangoni flow cannot be measured at the molecular scale [21]. Malijevský and Jackson [22] concluded the mechanical route to pressure is unreliable and should be avoided in systems with thermal fluctuations. In an attempt to address the confusion in the literature, Varnik et al. [14] compared different forms of pressure, considering the Harasima and Irving Kirkwood contour along with the method of planes on a flat surface. In this work we aim to extend this comparison to an arbitrary liquid-vapour interface, working with the VA pressure and a generalisation of the MOP pressure we call the surface flux (SF) form, we then discuss the relationship between these forms and the Harasima and Irving Kirkwood contours.
This consideration of pressure on a general interface shifts the focus away from the interaction contour between two molecules. A focus on the contour of interaction ignores the range of other factors in defining the pressure tensor. For example, the split of the kinetic term into streaming and fluctuating part is not free from ambiguity [24], a point which takes on new importance for a reference frame moving with the interface. In addition, the measured pressure is only correct to a gauge, meaning any divergence-free term could be added and the equations of motion would still hold [11]. Most importantly for interfaces, both the interaction contour and the location of measurement are required to determining the pressure, a point highlighted in the original work of Irving and Kirkwood [1] who note the force acting and the location it acts across, dS, are quite arbitrary. In comparing Harasima and Irving Kirkwood forms, Varnik et al. [14] note the interactions will be counted upon crossing an area, showing the contour path definitions are intertwined with the location of measurement. A contour could therefore be chosen which avoids being counted because it misses the area used to measure pressure. Here the SF formulation has a clear advantage, as it measures pressure on six connected surfaces enclosing a volume, for a molecule in the volume it is impossible for any choice of contour to not cross at least one surface. The location of SF pressure measurment is therefore not a plane but a volume defined by the set of connected patches which encloses it. The measurement location is typically defined to be a located at the centre of the liquid-vapour interface [15], the surface of a sphere [25,26], cylinder [27], and so on. A number of recent publications have shown we can track the liquid-vapour interface instantaneously down to the intermolecular spacing, [28][29][30][31][32], moving with the molecules as they evolve in time. Instead of considering the average of a plane (or sphere, cylinder, etc), in this work we use a pressure measurement which follows the instantaneous surface described by an arbitrary function ξ = ξ(x, y, t). This is mathematically tedious, but can be done for both VA [33] and SF pressure [34]. Using a measurement relative to an arbitrary surface gives the most general possible form of interface pressure, where ξ could then be chosen to be a flat surface, spherical segment, etc. As a result, these are valid for any surface away from equilibrium at every timestep. We provide a comparison of the resulting profiles for these SF pressures to the expression obtained inside a volume (IK1, VA) using both a flat plane and one moving with the instantaneous surface. Results for both the solid-liquid interface and liquid-vapour interface are compared and it is shown that Harasima and Irving Kirkwood contours are a special case of the more general VA and SF pressures. Through this comparison, the importance of explicit consideration of the reference frame is highlighted. The use of a reference frame fluctuating with the instantaneous interface, changes the pressure definition from thermodynamic to purely mechanical by accounting for every single force and flux. With every crossing counted, measured SF pressures on all surfaces of a volume can therefore be shown to exactly balance the momentum change in the volume, to machine precision [34]. This addresses at least one concern of thermal fluctuations invalidating the pressure measurement [22]. The long-range contributions [35] as well as three body [10] and greater interactions are not considered in this work, although the exactly balance on a control volume might prove useful in deriving these cases.
This manuscript is structured as follows, we start with an overview of the various mathematical forms of the pressure tensor in section 2, with a novel correction term for the VA pressure derived at the end. Next, the molecular setup is discussed, for the two cases simulated in this work, a solid-liquid interface for reference and the liquid-vapour interface in section 3. The results and discussion from these cases are presented next in section 4, with particular focus on the results taken in a reference frame moving with the interface, followed by the conclusions in section 5.

Theory
We briefly present the theoretical developments that lead to the different forms of pressure tensor.

Irving Kirkwood Pressure Tensor
The derivation of the pressure tensor is given in the work of Irving and Kirkwood [1], taking the time evolution of momentum and comparing forms with the expected continuum equations. The resulting pressure tensor at point r in the fluid is of the form, where angular brackets a; f denote the phase space average of a, velocity of molecule i is given byṙi and the streaming velocity is u. The distance between the positions of molecules i and j is defined by rij = ri − rj and f ij = −∂φij/∂rij denotes the inter-molecular force (the force exerted by i on j). The first term is the kinetic part of the pressure tensor, which includes the contributions due to molecular velocity. The second term is the configurational pressure, due to the interactions between molecules. The Oij term is known as the IK operator [11], obtained from the difference between the Dirac delta functions for molecule i and j, where Oij denotes a Taylor expansion operator, In what follows, we drop the phase space average so we can obtain expressions that are valid instantaneously in a molecular dynamics simulation. Introducing the definition p i /mi =ṙi − u, Eq. (1) therefore becomes, Where the sum notation N i,j = N i=1 N j =j has been introduced for conciseness.

Irving Kirkwood One (IK1) Tensor
First we consider Eq. (4) with only the first term in the IK operator of Eq. (3), know as the IK1 approximation. The result is the tensorial form of the virial pressure [2] but applied locally, The IK1 form of pressure is the most widely used in the molecular dynamics literature due to its simplicity, and at the time of writing is the default in LAMMPS, calculated from a binning of the per-atom stress [36]. Using the virial pressure locally in an inhomogeneous system is incorrect, as interactions with the surrounding fluid are not included [37]. This can also be interpreted as a consequence of neglecting terms of higher order in Eq. (3) so effects of local inhomogeneity in the fluid are lost. Away from equilibrium a localised description is required and the full Eq. (3) expression must be retained.

The Contour Forms -Irving Kirkwood and Harasima
In order to retain the higher order terms of Eq. (3), the delta functions in Eq. (2) can be reformulated as a contour integral between the two molecules. In order to do this, we rewrite the IK operator using the fundamental theorem of contour integration, where the integral along represents an infinite number of paths of integration between the atoms i and j [18]. In the classical paradigm of molecular dynamics, especially for the pairwise potentials considered in this work, Occam's Razor dictates that the interaction is a straight line between molecules. This is also in line with Newton's assumption of impressed force between two points. This assumptions results in = λrij with 0 < λ < 1, where d = rijdλ. so the configurational pressure is, Now as δ (r − ri − λrij) = δ (x − xi − λxij) δ (y − yi − λyij) δ (z − zi − λzij) , the Irving Kirkwood contour can be seen to be the special case when delta functions in the x and y directions are set to one and we get the resulting form of pressure over a flat interface in z by integrating along the line of interaction, with H the Heaviside function. This choice of a straight line contour in z, crossing a fixed area, appears in the appendix of the original paper by Irving and Kirkwood [1] to derive, with components, This form of the pressure is known in the literature as the Irving Kirkwood contour [15]. It has been shown that the normal pressure Eq. (11) is equivalent to the normal part of the MOP expression [14], however, the tangential parts of Eq.
(12) appears to be a hybrid of the VA and MOP forms. The pressure expression fxijxij of Eq. (12) is the same as the tangential VA component but counted like the normal MOP contribution, when a crossing occurs on the z plane. The other common formulation of pressure in the interface literature is based on the work of Harasima [16], which adjusts the interaction contour between the molecules based on the location of the interface of interest. The contour splits the path into a tangential and a normal component to a surface, so for a flat interface, the resulting form of stress is, For more complicated surfaces, with position in the z direction described by a function ξ = ξ(x, y, t), the normal n = ∇ξ/|∇ξ| and tangents t1 and t2 would be used to give an interaction path moving along the tangential and then normal direction. The general contour of Eq. (6) can be considered to be = t + n with t and n defined as tangent and normal to the interface function ξ. This choice of contour has a side effect. The inter-molecular interaction between the molecules is based on the current shape of the liquid-vapour interface and the reference frame we fit to it. It is argued here that it is more natural to consider the inter-molecular interaction to be the same regardless of system and instead explicitly consider the reference frame changing based on the interface. We will then obtain expressions which are similar to the Harasima contour but derived in terms of the reference frame. To do this, we integrate over a volume as in Lutsko [9] and Cormier et al. [8], which leads to the so-called Volume Average (VA) pressure formulation.

Volume Average (VA) Pressure
The volume average (VA) pressure is obtained from a volume integral of Eq. (4), assuming linear inter-molecular interactions using Eqs. (6) and (7). The location, size and shape of the volume is arbitrary and can be chosen wherever we want to measure the pressure. The kinetic pressure includes the kinetic contributions of molecules inside the volume, while the configurational term is based on the fraction of an inter-molecular interaction that passes through the given volume in space. This has the advantage that it is equivalent to the virial when all volumes in space are added up, as all fractions are included, where ∆V is the local volume; the ϑi function is only non-zero for a molecule inside the volume and the ϑ λ is non zero when a fraction of the interaction is inside the averaging volume. The ϑ function is simply the integral of the Dirac delta function over a volume [33,38,39], ϑα = ΛxαΛyαΛzα where with H the Heaviside function, x + the top of a volume in x, x − the bottom and xα the position of either a molecules when α = i or a point on the line of interaction when α = λ so x λ = xi + λxij. The difference between two Heaviside functions is known as a boxcar function, shown in figure 1, a function which is one between two points and zero outside. An advantage of the VA approach is we can choose any volume, including one which is moving with the interface ξ in the z direction, so this z boxcar function becomes Mathematically the value of zα, with molecule α = i or line α = λ, can be interpreted as shifted based on the surface ξ at the same xα and yα location , where this form follows directly from the mathematical derivation [33]. The shape of a typical volume with the z surface following the interface can be seen in figure  1. As a result, it has some similarities to the Harasima approach, expressing terms parallel to the current interface. We can simplify the integral along λ in the Harasima case by noting that if molecule i is in the volume, the integration of a contour over all components tangential to the surface must also be inside. Assuming the integral of the Λ λx and Λ λy functions give the line segment in the volume in x, which is ∆x, and in y, which is ∆y, and taking P C T as the sum of the xx and yy components, This is similar to the tangential pressure part of the Harasima contour, but derived assuming the reference frame is moving with the surface. Finally, it is worth mentioning that nothing in the integration process over a volume depends on the contour, so we could write the general form of contour integral, where the length of a contour in a volume is obtained from the integral over is the x component of the contour, similarly for y with the moving z surface, . In practice, this would be obtained in a piecewise manner in a computer code, stepping along a contour and evaluating ϑ for each segment. Different contours would give different results, for example a few contours are given in Figure 1, which each show different lengths inside the volume. The next section expresses pressure as surface crossings, which avoids this ambiguity, counting contributions entirely based on which surface is crossed.

Method of Planes (MOP) and Surface Flux (SF) Pressure
The surface pressure approach avoids the expansion in delta functions of Eq. (3) by considering the interaction over a plane; introduced empirically by Tsai [37] and derived formally by a Fourier transform of the Irving Kirkwood expression, Eq. (1) in Todd et al. [11]. This MOP form has the ability to deal with systems arbitrarily far from equilibrium [11]. It also gives three components of the stress tensor acting across a plane, when the plane normal is aligned along the x axis this is, where sgn(x) is the signum function. This two signum functions have the interpretation of including inter-molecular crossings over the surface of the plane. The normal component Pxx uses the normal components pix and fxij, and is equivalent to the IK form of pressure in the normal direction given in Eq. (11) (as sgn(zij) = zij/|zij| is cancelled by the zij denominator in Eq. (11), see [38]). The other two components of momentum and force provide the two other pressure components, Pxy and Pxz, and returns a single value for a plane covering the whole domain. Han and Lee [12] used three mutually perpendicular planes converging at a point to obtain all nine components of stress, with planes limited to a local region of interest. This allows us to write the tangential components of pressure Pyy, Pyx and Pyz as, where Λxi and Λx(λ k ) are boxcar functions which localise the plane to a patch in space. These check if the molecular i and point of crossing (λ k ) are between the limits of the surface in x (similar for y) with It is instructive to compare this pressure, where Pyy is the y component of force on the y surface, with the Irving Kirkwood contour of Eq. (12) which measures the y force times yij but on the z surface. This difference will be small for thin volumes but could become significant as ∆z increases. The local MOP pressure proposed by Han and Lee [12] occurs naturally from the derivation in the original paper by Irving and Kirkwood [1], by working in integrated, or control volume, form [38]. This proceeds by integrating the Irving and Kirkwood [1] momentum over a volume and evaluating the time evolution of this volume [38]. The resulting equations have a number of advantages, the pressure on the six volume surfaces can be used to get all nine components as in the Han and Lee [12] form. The derivation provides a natural link between VA and MOP pressures [38] while the use of a closed volume means any choice of contour must be included on one of the surfaces (see Figure 1). This directly addresses the non-uniqueness of interaction contour, as all interactions must be included provided one molecule is inside the volume and the other outside, so contour choice simply changes which surface an interaction is counted on. This has the effect of shifting it from a shear pressure contribution to a direct contribution (or vice versa), and becomes similar to choosing a different slice through a material. A notable corollary of this is that while ∇ · P = 0 is satisfied, shear components e.g. ∂Pxy/∂x may be non-zero for certain choices of contour. We can therefore verify direct components are constant and shear components are zero (required in a Newtonian liquid by definition), as a test that the choice of contour is meaningful. The other important advantage of the control volume approach is by defining the surfaces enclosing a volume, the sum over all surfaces is exactly equal to momentum change inside the volume. The pressure equilibrium condition is generalised to a form valid every timestep, so becomes ∇ · P = dρu/dt + ∇ · ρuu in control volume form. This allows us to derive a general surface pressure expression we are sure includes all terms for surface movement and curvature, by checking the sum of forces and momentum on a volume balance exactly. It is this property that has been used to validate the result for pressure presented in this work.
In order to generalise the MOP form to a surface, we introduce the notation dS + xi = δ(xi − x + )ΛiyΛiz and dS + xλ = δ(x λ − x + ), Λ λy Λ λz for the top surface (denoted by the + superscript) in the x direction, and it can be shown that, for a flat interface. The use of peculiar momentum p i has been dropped and the convective terms moved to left-hand side to avoid making the assumption that streaming velocity is constant. Equation (20) can therefore be seen to be the localised MOP shown in Eq. (19), which we call the Surface Flux (SF) pressure, because it can be generalised to get fluxes on any surface, including ones which are not flat. The Λαz term is a function of the moving surface ξ, so the size of the rectangular part of the plane it selects is constant, but moves as the surface changes. As a result, it is always the same distance from the surface and the same width (see Fig 1). For interactions crossing the arbitrary surface itself, ξ, the SF pressure is expressed as follows, where here the dS + zα = δ(ξ + α − zα)ΛxΛy with ξ + α = z + + ξ(xα, yα, t) catching the point of crossing if it is on the interface surface. We have a number of additional terms due to local surface curvature at the location of crossing. These include the crossing of a molecular trajectory and the crossing of an intermolecular interaction for kinetic and configurational parts respectively, as well as a surface time-evolution term.
For the case when ξ is constant or zero, the expression of Eq. (22) will simplify to the localised form of the MOP pressure given in Eq. (20), so this expression can be used to get pressure on every surface of a volume. The different surfaces of a volume moving with an interface are shown in Figure 1. The expression to actually calculate the pressure on a surface in an MD simulation, here the z + surface, is, where r12 = r2 − r1 is the line of time evolution of a molecule between t and t + ∆t, ϑt = [H(ξ(t) − zi(t)) − H(ξ(t − ∆t) − zi(t))] ΛixΛiy counting how many molecules have left or entered the volume due to surface movement between time t − ∆t and t and the surface normal includes all of the surface curvature terms of Eq. (22) with nz = ∇α(ξ−zα) ||∇α(ξ−zα)|| . This normal is the one used in the Harasima contour for a general surface, which could be combined with tangential contour of Eq. (16). We have generalised the treatment of both kinetic and configurational terms in Eq. (23) by recognising that molecules must move from a starting point r1 = ri(t) to an end point r2 = ri(t + ∆t) for surface crossings to occur. These movements can be written as an integral between two points in space (the molecule at different times), which makes it identical in form to the configurational term, an integral between two points (different molecules) in space. Both kinetic and configuration terms are treated identically, the crossing of a moving molecule or the crossing of an interaction are mathematically and algorithmically equivalent. In the case of a straight line and general surface ξ here, this becomes a ray-tracing problem which is common in computer graphics. Once the location of crossing, λ k , is obtained, it can be inserted in the following expression, with the first two Heaviside functions selecting if the root λ k is between the endpoints r1 and r2 . The Λ functions then check if that crossing is in a rectangle area of the surface, see Figure 1 for a graphical illustration of this. Multiple crossings are possible as the surface function ξ is general. Indeed, the form of crossing need make no assumptions about interaction contour or shape of surface, simply requiring a surface crossing λ k to be determined to evaluate dS + . It seems reasonable that this approach could also be applied to get the intersection of any contours, of the form given in Eq. (6), with any surface. The pressures tensors defined from a closed volume of bounding surfaces could then take advantage of the property already discussed which ensures any missed contributions from differing contour paths would be included as shear components on another surface.

A Correction to the Volume Average
In this section we present three extra terms which must be included to get the correct form of VA pressure on a moving interface. The VA pressure Eq. (15) is commonly obtained by integrating Eq. (4) over a volume [8,33] . However, we derive the SF equations from the time evolution of momentum inside a control volume, which results in extra terms for both the time evolution of the moving volume and the surface curvature. For the first term, where we have used ∂ξi ± /∂zi = 0 to write the expression concisely in vector form. We identify the kinetic part with reference to Eq. (15) using an overbrace, under the assumption that everything inside the ∂/∂r operator is VA pressure, and identify the two additional terms, KC and T E, missing from the VA treatment.
The second term in Eq. (25) is the configurational part, which for the VA in a moving interface, swapping ∂ϑ λ /∂r λ to ∂ϑ λ /∂r results in additional terms [34], so Figure 1.: A schematic showing a single control volume and the action of the different mathematical terms in the SF expression. The location of a surface crossing, λ k , is checked to be between the two molecules using the difference H(1−λ k )−H(−λ k ) while the various Λ expression act to limit the crossing to within a rectangular region of the surface. The normal at the point of crossing is used in the calculation of the surface term. The top right shows different example contours between molecules including the Irving Kirkwood as a solid line and two interpretations of the Harasima as different dotted lines assumed to move tangentially to the surface until it gets to a line either normal (which is consistent with the definition but can have multiple solutions for a arbitrary surface) or along the z axis the full expression with correction term CC is, again using ∂ξ ± λ /∂z λ = 0 to write in vector form. So, Eq. (25) shows the link between SF and VA form with the extra terms naturally included in the SF derivation, where the time evolution is exactly equal to the sum of SF over all six faces (convective terms are set to zero for simplicity) but only approximately equal to the divergence of the volume average with correction terms 1 . Pressure in the VA form is the average inside a divergence operator ∂/∂r, so it is natural to write the extra terms in Eqs. (26) and (27) in the same way. Terms in the x and y directions are expected to be zero so we take the derivative of the integral in just the z direction, so everything inside the derivative is in the form of a pressure contribution inside a volume. As with other VA terms, the functions ϑα assigns pressure contributions only when molecules, α = i, or points on interactions path, α = λ, are inside the volume. The interface movement ∂ξi/∂t and surface curvature ∂ξi/∂ri at the location of the molecules is used in the kinetic part, while the curvature integrated along the contour between molecules ∂ξ λ /∂r λ appears in the configurational part. Obtaining curvature at varying locations while moving along an interaction contour will likely be complex and computationally expensive. For simplicity, these terms will be estimated using the surface flux terms from Eq. (22) in this work as these are already obtained when collecting SF expressions and are expected to be similar for the thin averaging volumes employed.

Methodology
In this section, we outline the molecular dynamics (MD) setup for the three cases shown schematically in Figure 2. All cases use the simple, pairwise shifted Lennard Jones potential, with a cutoff rc = 2.5, where all numbers are presented in reduced LJ units. This is chosen to allow efficient simulations but is shorter than required to give surface tension with good experimental agreement [40]. Time integration uses the velocity Verlet scheme, ri(t + ∆t) = ri(t) + ∆tvi(t + ∆t/2) vi(t + ∆t/2) = vi(t − ∆t/2) + ∆tF i(t) (31) with a timestep ∆t = 0.005 and force on i calculated from F i = N j =i f ij . The simulations are run using Flowmol which has been verified in a range of previous publications [34,41,42] and has recently been made open-source [43].
For the solid-liquid simulation, the domain is setup with height Lz = 34.2 and walls of thickness 3 from the top and bottom of the domain. The system is periodic in the other directions, with Lx = Ly = 13.6 giving a total of N = 5120 molecules. The walls are fixed using tethered molecules with anharmonic spring constants of Petravic and Harrowell [44] with strength k4 = 5 × 10 3 and k6 = 5 × 10 6 . The system is initialised using an FCC lattice of density ρ = 0.8 throughout and initial temperature of T = 1, with the untethered region allowed to melt while the system is equilibrated. After melting, the average temperature of the system is around Tave = 0.6. Statistics are then collected over 100,000 steps.
For the liquid-vapour coexistence, the setup is identical to previous work [34], with periodic boundaries in all directions of size Lx = Ly = 12.7 and Lz = 47.62. The setup also uses an FCC lattice and deletes molecules, with the middle 40% designated as liquid with a density of ρ l = 0.79 and remaining gas at density ρg = 0.0002 giving N = 2635 molecules. The system is equilibrated for 100,000 timesteps in a Nosè Hoover NVT ensemble with setpoint Ts = 0.7. The main runs are restarted in an NVE ensemble run until well-resolved statistics are obtained.
The intrinsic interface is fitted to the outer molecules of a cluster of liquid molecules, identified by a cluster analysis with Stillinger cutoff length r d = 1.5 and requiring more than three neighbours per molecule. The functional form of surface ξ(x, y, t) is fitted using the intrinsic surface method (ISM) [28] where the interface is described by an arbitrary number of sines and cosines, where fµ(x) = cos(kxx), f−µ(x) = sin(kxx) and aµν is the matrix of surface wavenumbers. The number of wavelengths is calculated from the system size ku = nint( √ LxLy/λu) with λu = 1 the minimum wavelength. This function is fitted by minimising the least-square difference between surface molecules zp and the intrinsic surface function at these positions ξ = ξ(xp, yp, t), with extra constraint ψ = 1 × 10 −8 to prevent overfitting by ensure intrinsic areã A does not become too large. This process proceeds in stages, starting with a 3 by 3 grid of the most extreme molecules, the fitting is repeated, using new pivots added in batches based on proximity to the current surface, until the density of the surface reaches a value of ns = 0.7. Once fitted, this surface is then used to determine the grid of cells, or bins, which we use to collect averaged statistics. This is the Lagrangian reference frame used in to obtain pressure measurements which are always the same distance from the current position of the moving surface. The domain is split into bins with dimensions ∆x = ∆y = 0.198 tangential to the surface and ∆z = 0.0875 in the normal direction. For the SF term, the intrinsic interface is converted to a piecewise set of bilinear squares each of size ∆x by ∆y which are used to calculate the interaction with the surface. This employs an efficient algorithm for ray-tracing to obtain crossings [45] and assigns these to the appropriate cells. The conservation of every control volume in the domain is  has a difference of less than 1 × 10 −9 compared to the momentum change inside that volume. For the VA pressure, a mapping based on molecular positions, or points on the interaction contour between molecules, is applied using the corresponding surface position. The mapped coordinate can then be efficiently binned using integer division, as if on a uniform grid. The full details of these algorithms are described extensively in previous work [34]

Results and Discussion
In this section, we compare three types of pressure measurements, the IK1 Eq. (5), VA Eq. (15) and SF Eq. (22), split into kinetic and configurational components. These three forms also provide the Irving Kirkwood and Harasima contours, and the link to these is discussed. The convention is the same in all figures, kinetic terms are shown in orange, configurational in blue and total pressures in green. The IK1 pressure is shown as a dotted line, the VA as filled circles and the SF as a solid line. The three cases are shown in Figure 2, first 2a) the solid-liquid case is shown in Figure 3, which is used to demonstrate the flat normal pressure profile obtained by both the VA and SF methods, as required for mechanical equilibrium, together with the varying IK1. The same set of pressure measurements of case 2b) the liquid-vapour interface using a fixed averaging grid are shown in Figure 4, demonstrating the same flat normal pressure profile from both VA and SF terms and the same variation in the IK1 expression. Finally, 2c) is the focus of the remaining work, using interface tracking to obtain pressure in a reference moving with an intrinsic liquid-vapour interface. All of the IK1, VA and SF forms of pressure give different behaviour, and the kinetic and configurational terms are compared in Figure 5, before discussing the total pressure in Figure 6. In order to understand the differences, the normal and tangential parts are presented individually and the resulting contributions to surface tension from their respective integrals are displayed in Figure 7. A full breakdown of all SF contributions to normal pressure from Eq. (22), together with the surface curvature and movement which are needed as corrections to the VA form, then given in Figure 8. This is shown to satisfy the momentum equilibrium in a moving control volume, proving the SF form must contain all possible contributions and that the curvature and surface movement are essential to ensure VA expressions satisfied mechanical equilibrium in the reference frame moving with the interface. The shape of tangential pressure is shown normalised by density in Figure 9 with a fitting proposed which could be used to approximate pair interactions in terms of density. The link to the Irving Kirkwood and Harasima contours is discussed. We start with the case shown in the schematic of Figure 2 a). The wall normal (PN = Pzz) pressure near a solid-liquid interface is presented first in Figure 3 with half of the channel shown in a) and a focus on just the near-wall region in b). The shaded region in Figure 3a) shows the area focused on in 3 b), as it contains the majority of the variation, and in all subsequent figures only a small region near the interface is displayed. The kinetic and configurational parts of the pressure are shown to balance for both the VA and SF forms of pressure, with the average giving a constant pressure value. This satisfies the condition for equilibrium ∇ · P = 0. The IK1 pressure has identical kinetic component to both the VA and SF, but the configurational part is seen to have larger shifted peaks which results in a non-zero average pressure so ∇ · P = 0. These peaks in the pressure cannot be correct in an equilibrium system, as they would induce a flow, so this suggests the IK1 pressure is not correct, a well documented result in the literature [11,14].
The same quantities discussed for the solid-liquid case are shown for the liquid-vapour interface using a fixed grid in Figure 4, the schematic of Figure 2 b). As with the solid-liquid cases, the VA and SF gives similar results while the IK1 shows a difference in the normal configurational pressure, the blue line in fails to give the correct normal component, although the difference can be seen to have roughly equal positive and negative areas. This means the Kirkwood and Buff [46] integral used to get surface tension will be the same as the other methods. The tangential pressure 4 b) is almost identical for IK1, VA and SF measurements. The IK1 shows a slight shift toward the interface when compared to the other forms of pressure. We do not explicitly calculate the tangential Irving Kirkwood contour of Eq. (12), but it is a combination of the VA form measured on a SF style z plane, so we would expect it to give the same results as both the VA and SF expressions, which are very similar in Figure 4 b).
So, the conclusion from Figures 3 and 4 are that the VA and SF satisfy force balance in the normal direction while the IK1 does not. However, the IK1 PN contribution above and below the axis are roughly equal so this will cancel in an integral. As the tangential components of the IK1 is similar to both the VA and SF, the integral of PN − PT used in surface tension will be identical. For the normal pressure, the Irving Kirkwood contour is represented by the SF results and satisfy the equilibrium ∇ · P = 0. For the tangential pressure, the Irving Kirkwood contour is a combination of the VA and SF, which give very similar results. There is very little difference in these various measures for a fixed reference frame, so next we consider the same pressure measurements in a reference frame moving with the interface, the case shown in Figure 2 c).
Both kinetic and configurational components of pressure are shown in Figure 5 with the total pressure shown in Figure 6, again split into normal a) and tangential components b). As before, the IK1, VA and SF terms are compared, but here we see notable differences in all three. For the kinetic terms in Figure 5a), the IK1 and VA agree, showing peaks which mirror the density shown in gray, representing the location of the molecules relative to the surface. The SF term is shown and naturally includes the surface evolution term, ∂ξ/∂t and kinetic curvature which are not included in the VA expression. The kinetic curvature terms are zero on average, but the inclusion of the surface evolution is important, without it the kinetic SF pressure profile would be identical to the VA and IK1 [34]. As a result of including the surface movement, the SF profile has a flat region on the liquid side (orange line Figure 5a)), between the interface and first layer, as well as a smoother transition to zero on the gas side. The profile is exactly mirrored by the SF configurational pressure, the solid blue line on 5a), so the sum of kinetic and configurational terms gives a flat pressure profile on Figure 6a), the solid green line. We therefore see from Figure 6a), only the SF form satisfies the equilibrium force-balance condition. The extra term in Eq. (27) for the configurational curvature is also missing from the configurational term of the VA pressure, a point we return to at the end of this section. As a result, the VA pressure does not satisfy the equilibrium force-balance condition in 6a) in a moving reference frame. The IK1 form is very different to the VA and SF, showing large oscillations in normal configurational pressure which results in significant peaks in the configurational pressure of Figures 5a) and therefore the total in 6a). As the IK1 pressure fails to satisfy the balance condition in the simpler case of a solid-liquid interface, as well as the fixed reference liquid-vapour case, there is no reason to suspect it will perform any better for a moving-reference case. In fact, the oscillatory nature of the IK1 pressure appears to show more prominent peaks which is a feature of the IK1 observed in the solid-liquid case of Figure 3. This is interesting as the IK1 departs further as the system becomes more inhomogeneous and structured [14], so the IK1 performs worse here as a moving reference frame exposes the inhomogeneity of the molecular structure near the intrinsic interface. These molecular peaks are apparently smoothed by the continuous operation of taking a fraction of line inside a volume for the VA case, while the SF expression depends on crossings and so is less dependent on the absolute position of the molecules. The measuring volumes are also centred on the interface itself so surfaces crossing measurements are offset either side by ∆z/2. The SF therefore avoids the peak at the interface in the normal pressure measurement seen in Figure 5 a). The tangential components are identical for VA and SF in 6b), as they are both calculated using a uniform grid in the surface tangent direction, so we would expect them to agree in the same way any fixed volume cases do. The IK1 actually shows smaller oscillation than the VA/SF terms in the tangential direction with a peak at the interface and a smaller, more smeared contribution for the first and second fluid layers.
Despite the very different profiles for the three pressure measurements, it is well-documented that the Kirkwood and Buff [46] integral will give the same surface tension for all three [22]. This integral of the difference between normal and tangential terms will also be the same as using a fixed reference frame [34]. One way to understand this: all definitions of pressure considered here use the same forces but vary how they distribute these measurements to different bins or cells, so the integral must be the same. As the Kirkwood and Buff [46] integral is using the difference between normal and tangential components (PN − PT )dz, it is instructive to split them into these two additive contributions, as is done in Figure 7 with a) the integral of the normal PN dz and b) the integral of the negative of the tangential − PT dz. In Figure 7 a), the flat SF normal pressure profile gives a linearly increasing integral over the interface which agrees with the oscillating VA and IK1 after the integral has moved far enough away from the interface into the gas phase (shown approaching z = 2 in the insert). The integral of the normal contribution is small, contributing less than 10% (0.05 by z = 2) to the total surface tension integral which is approximately 0.6 in this case). The main contribution is from the tangential pressure which shows an identical trend for both the VA and SF in Figure 7 b). The tangential part of the IK1 pressure can be seen to have fewer fluctuations and therefore a simpler integrated contribution. Again, all three methods converge to approximately the same integrated result for the tangential contributions, with a value of about 0.55 by z = 2. Therefore, the measured surface tension, obtained from the sum of the normal, 7a) and minus tangential 7b), converges to roughly the same value of γ = PN − PT dz ≈ 0.6 for all methods. Differences are attributed to statistical noise, which integrated quantities like surface tension are more susceptible to. Figure 8a) shows the VA pressure and the three additional terms, KE, TE and CC, which were derived as corrections to the VA pressure in section 2.6. These terms are due to the moving interface and would be zero for a fixed reference frame, which is why the VA and SF normal pressures have identical profiles in Figures 3  and 4 but not in Figures 5 and 6. The corrected VA forms are shown to give a constant normal pressure in Figure 8b), as required for mechanical equilibrium in a moving interface, as well as good agreement with the SF pressure which naturally include these terms. For the kinetic pressure, the normal VA contribution is shown by orange circles on Figure 8a) and the additional correction terms due to kinetic and configurational interface curvature as well as the movement of the surface itself, are shown with lines because they are obtained from SF terms. The VA shown by orange circles on Figure 8b) includes the surface movement correction, the purple line on 8a), and the kinetic curvature which is the red line (both extra terms from Eq. (26)). With these additional corrections, the kinetic part of the VA has identical shape to the SF pressure which naturally includes these terms. For the configurational pressure, the brown line on 8a) shows the magnitude of the configurational curvature contributions, which are the extra terms missing from the VA measurements in Eq. (27). The total for both SF and corrected VA is shown in green as a flat normal pressure which is indicative of mechanical equilibrium. However, the VA shows a few values which are not flat near the interface, attributed to using the SF form of pressure to get the correction terms instead of the VA form. This was expected to gives a small error for the thin bins used here, but near the rapid change at the interface the limits of this assumption appear to show. This indicates how important getting an exactly a consistent form of pressure is, with only the SF able to give a mechanical balance at all points.
Finally the tangential pressure is considered in Figure 9, where the value is plotted normalised by density. The peaks in the kinetic pressure PT (z) and density ρ(z) correspond exactly, as a function of z, so the ratio is constant giving a flat profile The shape of the total pressure is therefore due to the configurational part, which depends on interactions passing through a location in space, and will be non-zero even where no molecules are located. As a result, a diverging profile is observed either side of zero, as the intrinsic fitting process results in a large density peak at the interface and no molecules either side. This mean PT /ρ → ∞ as ρ → 0 but a finite value is observed at zero due to the molecules which sit on the intrinsic surface. The shape of the profile represent the pressure per molecule, so it is useful to understand the functional form it takes. Fitting is obtained using a decaying exponential e −az times an oscillating function 1 − b sin(2cπz + d), where a, b, c and d are coefficients to be fitted, where a reasonable fit is observed for values a = 1, b = 2.5, c = 1 and d = 0.5, so tangential pressure in the liquid can be approximated using, This form seems to be a reasonable approximation for the first two peaks and captures the decay to zero moving into the bulk. This fitted expression could be used in a course-grained approximation or mean-field approach to model the pressure near the interface. Finally we consider the link between results obtained here in a reference frame tracking the interface and the contour forms. The Irving Kirkwood contour assumes a fixed reference, so is not trivially related to the forms of pressure shown in Figure 5 and 6 for a moving reference frame. However, the Harasima contour is obtained by moving first tangentially to the surface, then in the normal direction. The tangential part from Eq. (16) is therefore a similar quantity to the tangential IK1 components shown in Figures 5b) and 6b). This shows much smaller oscillations in 6b) than the VA and SF methods. The normal component of the Harasima contour is mathematically similar to the SF form, which takes the contribution dotted with the surface normal. It can be seen in Figure 7 that the integral of the normal and tangential terms all converge to the same value, so mixing the tangential components of the IK1 pressure with the normal components of the SF to obtain the Harasima contour would result in the same surface tension. The Harasima contour would therefore satisfy the equilibrium condition in the normal direction, give the simple tangential profile of 6b) and return the correct surface tension. This work therefore gives the mathematical process to measure   this contour for a general surface, by taking the SF in the normal direction Eq. (23) and then the IK1 for the tangential components using Eq. (16). This process was first reported in the work of Sega et al. [47] with differences between the Harasima and IK tangential pressure explored in later work [35] However, for a general surface, the split of the Harasima contour into a normal and tangential part is non-unique, as highlighted in Figure 1 where different combinations of tangent and normal to the surface would be possible. As calculation of the SF interaction and surface normal are required to get the Harasima contour for a general interface, simply using the SF is preferable in this case, although for more complex potentials and long-range interaction the Harasima contour may offer additional advantages [27,35]. The SF provides three components per surface, so all nine pressure terms are available, also giving the curvature contributions and surface movement shown in Figure 8. It also provides an explicit mathematical link between any contour crossing the closed surface and momentum evolution inside that closed control volume. By mixing the tangential IK1 and SF, we lose this property of exact control volume balance and the guarantee all pressure components are captured. However, using the SF form moving with the interface provides a pressure which naturally includes all terms and satisfies the required condition of mechanical equilibrium.

Conclusions
The appropriate definition of pressure in molecular dynamics (MD) simulation has been the subject of some debate. In this work, it is argued that a consideration of the reference frame is more illuminating in defining pressure than the interaction contour. This is especially true near an interface. By considering a reference frame described by an arbitrary function ξ = ξ(x, y, t) fitted to the surface every time, we get a general form of interface pressure. The derivation provides an expression in terms of the surface fluxes (SF), a generalisation of the method of planes (MOP) pressure to an arbitrary surface. This is compared to a range of other pressure definitions, with the differences discussed. These include a contour between two molecules averaged in a volume called the volume average (VA) pressure, the virial applied locally with interaction split into two called the IK1 pressure as well as the Irving Kirkwood and Harasima contours. All forms of pressure give the same surface tension but result in different pressure distributions over an interface. Results show the IK1 form does not satisfy mechanical equilibrium in the interface normal direction, even in the simplest cases, while the VA fails to satisfy mechanical equilibrium in the case moving with a reference frame. Correction terms are derived for the VA form and shown to account for moving interface and curvature, giving the expected results. Combining the SF and IK1 results is shown to provide a generalisation of the Harasima contour for an arbitrary surface. As a result, we show choice of reference frame provides a clear link between the various forms of pressure. The SF is argued to be the preferred form as it naturally includes all curvature and surface movement terms, provides an exact link between SF pressure and instantaneous momentum change in the volume and as six SF pressures bound a closed volume, it ensures any of the infinite possible contour choices must cross one of the surfaces and be included in the pressure tensor.