The importance of reference frame for pressure at the liquid–vapour interface

ABSTRACT

The local pressure tensor is non-unique, a fact which has generated confusion and debate in the 70 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 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 confusion in the literature. Using an explicit reference frame with a more general derivation of pressure clarifies some of the issues of uniqueness, providing a pressure tensor which is defined at any instant in time and valid away from thermodynamic equilibrium.

KEYWORDS:

• Pressure tensor
• stress tensor
• liquid–vapour interface
• intrinsic interface

View correction statement:

Correction

1. 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-orderterm 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–7]. These solid mechanics stress tensors, where stress is simply the negative of the pressure tensor, are largely in the form of 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 occur in an infinite number of ways. This is explored for various contours in recent work [19]. 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 [20,21]. In addition, the use of Harasima and Irving Kirkwood surface forms led to the conclusion that Marangoni flow cannot be measured at the molecular scale [22]. Malijevsý and Jackson [23] 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 determine 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, $\mathrm{d}S$ , 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 measurement 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 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 that we can track the liquid–vapour interface instantaneously down to the intermolecular spacing [28–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 $\xi =\xi (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 time step. 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 [23]. The long-range contributions to the pressure tensor [35] as well as the inclusion of three body [10] and greater interactions are not considered in this work, although the exact balance on a control volume might prove useful in the extension to these cases.

This article 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 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.

2. Theory

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

2.1. 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 $\mathit{r}$ in the fluid is of the form (1) $\begin{array}{rl}& \mathit{P}(\mathit{r},t)=\\ & \sum _{i=1}^N〈m_i\left({\dot{\mathit{r}}}_i-\mathit{u}\right)\left({\dot{\mathit{r}}}_i-\mathit{u}\right)\delta ({\mathit{r}}_i-\mathit{r})+\frac{1}{2}\sum _{j\ne i}^N{\mathit{f}}_{ij}{\mathit{r}}_{ij}{O}_{ij}\delta ({\mathit{r}}_i-\mathit{r});f〉,\end{array}$(1) where angular brackets $〈a;f〉$ denote the phase space average of a, the velocity of molecule i is given by ${\dot{\mathit{r}}}_i$ and the streaming velocity is $\mathit{u}$. The distance between the positions of molecules i and j is defined by ${\mathit{r}}_{ij}={\mathit{r}}_i-{\mathit{r}}_j$ and ${\mathit{f}}_{ij}=-\mathrm{\partial }{\varphi }_{ij}/\mathrm{\partial }{r}_{ij}$ 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 $O_{ij}$ term is known as the IK operator [11], obtained from the difference between the Dirac delta functions for molecule i and j, (2) $\delta ({\mathit{r}}_i-\mathbf{r})-\delta ({\mathit{r}}_j-\mathit{r})=-{\mathit{r}}_{ij}\cdot \frac{\mathrm{\partial }}{\mathrm{\partial }\mathit{r}}{O}_{ij}\delta ({\mathit{r}}_i-\mathit{r}),$(2) where $O_{ij}$ denotes a Taylor expansion operator, (3) $O_{ij}=\left(1-\frac{1}{2}{\mathit{r}}_{ij}\cdot \frac{\mathrm{\partial }}{\mathrm{\partial }{\mathit{r}}_i}+\cdots +\frac{1}{n!}{\left(-{\mathit{r}}_{ij}\cdot \frac{\mathrm{\partial }}{\mathrm{\partial }{\mathit{r}}_i}\right)}^{n-1}+\cdots \right).$(3) 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 ${\mathit{p}}_i/m_i={\dot{\mathit{r}}}_i-\mathit{u}$, Equation (1(1) $\begin{array}{rl}& \mathit{P}(\mathit{r},t)=\\ & \sum _{i=1}^N〈m_i\left({\dot{\mathit{r}}}_i-\mathit{u}\right)\left({\dot{\mathit{r}}}_i-\mathit{u}\right)\delta ({\mathit{r}}_i-\mathit{r})+\frac{1}{2}\sum _{j\ne i}^N{\mathit{f}}_{ij}{\mathit{r}}_{ij}{O}_{ij}\delta ({\mathit{r}}_i-\mathit{r});f〉,\end{array}$(1) ), therefore, becomes (4) $\mathit{P}(\mathit{r},t)=\sum _{i=1}^N\frac{{\mathit{p}}_i{\mathit{p}}_i}{m_i}\delta ({\mathit{r}}_i-\mathit{r})+\frac{1}{2}\sum _{i,j}^N{\mathit{f}}_{ij}{\mathit{r}}_{ij}{O}_{ij}\delta ({\mathit{r}}_i-\mathit{r}),$(4) where the sum notation $\sum _{i,j}^N=\sum _{i=1}^N\sum _{j\ne j}^N$ has been introduced for conciseness.

2.2. Irving–Kirkwood one (IK1) tensor

First, we consider Equation (4(4) $\mathit{P}(\mathit{r},t)=\sum _{i=1}^N\frac{{\mathit{p}}_i{\mathit{p}}_i}{m_i}\delta ({\mathit{r}}_i-\mathit{r})+\frac{1}{2}\sum _{i,j}^N{\mathit{f}}_{ij}{\mathit{r}}_{ij}{O}_{ij}\delta ({\mathit{r}}_i-\mathit{r}),$(4) ) with only the first term in the IK operator of Equation (3(3) $O_{ij}=\left(1-\frac{1}{2}{\mathit{r}}_{ij}\cdot \frac{\mathrm{\partial }}{\mathrm{\partial }{\mathit{r}}_i}+\cdots +\frac{1}{n!}{\left(-{\mathit{r}}_{ij}\cdot \frac{\mathrm{\partial }}{\mathrm{\partial }{\mathit{r}}_i}\right)}^{n-1}+\cdots \right).$(3) ), known as the IK1 approximation. The result is the tensorial form of the virial pressure [2] but applied locally (5) $\stackrel{\mathrm{I}\mathrm{K}1}{\mathrm{P}}(\mathbf{r},t)=\sum _{i=1}^N\frac{{\mathbf{p}}_i{\mathbf{p}}_i}{m_i}\delta \left(\mathbf{r}-{\mathbf{r}}_i\right)+\frac{1}{2}\sum _{i,j}^N{\mathbf{f}}_{ij}{\mathbf{r}}_{ij}\delta \left(\mathbf{r}-{\mathbf{r}}_i\right).$(5) 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 Equation (3(3) $O_{ij}=\left(1-\frac{1}{2}{\mathit{r}}_{ij}\cdot \frac{\mathrm{\partial }}{\mathrm{\partial }{\mathit{r}}_i}+\cdots +\frac{1}{n!}{\left(-{\mathit{r}}_{ij}\cdot \frac{\mathrm{\partial }}{\mathrm{\partial }{\mathit{r}}_i}\right)}^{n-1}+\cdots \right).$(3) ) so effects of local inhomogeneity in the fluid are lost. Away from equilibrium, a localised description is required and full equation (3(3) $O_{ij}=\left(1-\frac{1}{2}{\mathit{r}}_{ij}\cdot \frac{\mathrm{\partial }}{\mathrm{\partial }{\mathit{r}}_i}+\cdots +\frac{1}{n!}{\left(-{\mathit{r}}_{ij}\cdot \frac{\mathrm{\partial }}{\mathrm{\partial }{\mathit{r}}_i}\right)}^{n-1}+\cdots \right).$(3) ) expression must be retained.

2.3. The contour forms – Irving Kirkwood and Harasima

In order to retain the higher order terms of Equation (3(3) $O_{ij}=\left(1-\frac{1}{2}{\mathit{r}}_{ij}\cdot \frac{\mathrm{\partial }}{\mathrm{\partial }{\mathit{r}}_i}+\cdots +\frac{1}{n!}{\left(-{\mathit{r}}_{ij}\cdot \frac{\mathrm{\partial }}{\mathrm{\partial }{\mathit{r}}_i}\right)}^{n-1}+\cdots \right).$(3) ), the delta functions in Equation (2(2) $\delta ({\mathit{r}}_i-\mathbf{r})-\delta ({\mathit{r}}_j-\mathit{r})=-{\mathit{r}}_{ij}\cdot \frac{\mathrm{\partial }}{\mathrm{\partial }\mathit{r}}{O}_{ij}\delta ({\mathit{r}}_i-\mathit{r}),$(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, (6) $\begin{array}{rl}\delta (\mathbf{r}-{\mathbf{r}}_i)-\delta (\mathbf{r}-{\mathbf{r}}_j)& =\oint \frac{\mathrm{\partial }}{\mathrm{\partial }\mathit{\ell }}\cdot \delta \left(\mathbf{r}-{\mathit{r}}_i-\mathit{\ell }\right)\mathrm{d}\mathit{\ell }\\ & =\oint \frac{\mathrm{\partial }{\mathbf{r}}_i+\mathit{\ell }}{\mathrm{\partial }\mathit{\ell }}\cdot \frac{\mathrm{\partial }}{\mathrm{\partial }\left[{\mathbf{r}}_i+\mathit{\ell }\right]}\delta \left(\mathbf{r}-{\mathit{r}}_i-\mathit{\ell }\right)\mathrm{d}\mathit{\ell }\\ & =\frac{\mathrm{\partial }}{\mathrm{\partial }\mathbf{r}}\cdot \oint \delta \left(\mathbf{r}-{\mathit{r}}_i-\mathit{\ell }\right)\mathrm{d}\mathit{\ell },\end{array}$(6) where the integral along $\mathit{\ell }$ 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 $\mathit{\ell }=\lambda {\mathit{r}}_{ij}$ with $0<\lambda <1$, (7) $\frac{\mathrm{\partial }}{\mathrm{\partial }\mathbf{r}}\cdot \oint \delta \left(\mathbf{r}-{\mathit{r}}_i-\mathit{\ell }\right)\mathrm{d}\mathit{\ell }=\frac{\mathrm{\partial }}{\mathrm{\partial }\mathbf{r}}\cdot {\mathit{r}}_{ij}{\int }_0^1\delta \left(\mathbf{r}-{\mathit{r}}_i-\lambda {\mathit{r}}_{ij}\right)\mathrm{d}\lambda ,$(7) where $\mathrm{d}\mathit{\ell }={\mathit{r}}_{ij}\mathrm{d}\lambda$. So the configurational pressure is (8) ${\mathit{P}}^C(\mathit{r},t)=\frac{1}{2}\sum _{i,j}^N{\mathit{f}}_{ij}{\mathit{r}}_{ij}{\int }_0^1\delta \left(\mathbf{r}-{\mathit{r}}_i-\lambda {\mathit{r}}_{ij}\right)\mathrm{d}\lambda .$(8) Now as $\delta (\mathbf{r}-{\mathit{r}}_i-\lambda {\mathit{r}}_{ij})$ = $\delta (x-x_i-\lambda x_{ij})\delta (y-y_i$ - $\lambda y_{ij}\left)\delta \right(z-z_i-\lambda z_{ij})$, 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, (9) ${\int }_0^1\delta \left(z-z_i-\lambda z_{ij}\right)\mathrm{d}\lambda =\frac{1}{\left|z_{ij}\right|}H\left(\frac{z-z_i}{z_{ij}}\right)H\left(\frac{z_j-z}{z_{ij}}\right)$(9) 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 (10) $\stackrel{1K}{{\mathrm{P}}^C}(\mathbf{r},t)=\frac{1}{2}\sum _{i,j}^N{\mathbf{f}}_{ij}\frac{{\mathbf{r}}_{ij}}{|z_{ij}|}H\left(\frac{z-z_i}{z_{ij}}\right)H\left(\frac{z_j-z}{z_{ij}}\right),$(10) with components (11) $\stackrel{\mathrm{I}\mathrm{K}}{\mathrm{P}}\underset{N}{\overset{C}{}}=\stackrel{\mathrm{I}\mathrm{K}}{\mathrm{P}}\underset{zz}{\overset{C}{}}=\frac{1}{2}\sum _{i,j}^N{f}_{zij}\frac{z_{ij}}{|z_{ij}|}H\left(\frac{z-z_i}{z_{ij}}\right)H\left(\frac{z_j-z}{z_{ij}}\right),$(11) (12) $\stackrel{\mathrm{I}\mathrm{K}}{\mathrm{P}}\underset{T}{\overset{C}{}}=\frac{1}{2}\left[\stackrel{\mathrm{I}\mathrm{K}}{\mathrm{P}}\underset{xx}{\overset{C}{}}+\stackrel{\mathrm{I}\mathrm{K}}{\mathrm{P}}\underset{yy}{\overset{C}{}}\right]=\frac{1}{4}\sum _{i,j}^N\left[\frac{f_{xij}{x}_{ij}+f_{yij}{y}_{ij}}{|z_{ij}|}\right]H\left(\frac{z-z_i}{z_{ij}}\right)H\left(\frac{z_j-z}{z_{ij}}\right).$(12) This form of pressure is known in the literature as the Irving Kirkwood contour [15]. It has been shown that normal pressure equation (11(11) $\stackrel{\mathrm{I}\mathrm{K}}{\mathrm{P}}\underset{N}{\overset{C}{}}=\stackrel{\mathrm{I}\mathrm{K}}{\mathrm{P}}\underset{zz}{\overset{C}{}}=\frac{1}{2}\sum _{i,j}^N{f}_{zij}\frac{z_{ij}}{|z_{ij}|}H\left(\frac{z-z_i}{z_{ij}}\right)H\left(\frac{z_j-z}{z_{ij}}\right),$(11) ) is equivalent to the normal part of the MOP expression [14]; however, the tangential parts of Equation (12(12) $\stackrel{\mathrm{I}\mathrm{K}}{\mathrm{P}}\underset{T}{\overset{C}{}}=\frac{1}{2}\left[\stackrel{\mathrm{I}\mathrm{K}}{\mathrm{P}}\underset{xx}{\overset{C}{}}+\stackrel{\mathrm{I}\mathrm{K}}{\mathrm{P}}\underset{yy}{\overset{C}{}}\right]=\frac{1}{4}\sum _{i,j}^N\left[\frac{f_{xij}{x}_{ij}+f_{yij}{y}_{ij}}{|z_{ij}|}\right]H\left(\frac{z-z_i}{z_{ij}}\right)H\left(\frac{z_j-z}{z_{ij}}\right).$(12) ) appears to be a hybrid of the VA and MOP forms. The pressure expression $f_{xij}{x}_{ij}$ of Equation (12(12) $\stackrel{\mathrm{I}\mathrm{K}}{\mathrm{P}}\underset{T}{\overset{C}{}}=\frac{1}{2}\left[\stackrel{\mathrm{I}\mathrm{K}}{\mathrm{P}}\underset{xx}{\overset{C}{}}+\stackrel{\mathrm{I}\mathrm{K}}{\mathrm{P}}\underset{yy}{\overset{C}{}}\right]=\frac{1}{4}\sum _{i,j}^N\left[\frac{f_{xij}{x}_{ij}+f_{yij}{y}_{ij}}{|z_{ij}|}\right]H\left(\frac{z-z_i}{z_{ij}}\right)H\left(\frac{z_j-z}{z_{ij}}\right).$(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 (13) $\stackrel{H}{P}{}_N=\stackrel{IK}{P}{}_N,$(13) (14) $\stackrel{H}{P}{}_T=\frac{1}{4}\sum _{i,j}^N\left[\frac{f_{xij}{x}_{ij}+f_{yij}{y}_{ij}}{|z_{ij}|}\right]\delta (z-z_i).$(14) For more complicated surfaces, with a position in the z-direction described by a function $\xi =\xi (x,y,t)$, the normal $\mathit{n}=\mathbf{\nabla }\xi /\left|\mathbf{\nabla }\xi \right|$ and tangents ${\mathit{t}}_1$ and ${\mathit{t}}_2$ would be used to give an interaction path moving along the tangential and then the normal direction. The general contour of Equation (6(6) $\begin{array}{rl}\delta (\mathbf{r}-{\mathbf{r}}_i)-\delta (\mathbf{r}-{\mathbf{r}}_j)& =\oint \frac{\mathrm{\partial }}{\mathrm{\partial }\mathit{\ell }}\cdot \delta \left(\mathbf{r}-{\mathit{r}}_i-\mathit{\ell }\right)\mathrm{d}\mathit{\ell }\\ & =\oint \frac{\mathrm{\partial }{\mathbf{r}}_i+\mathit{\ell }}{\mathrm{\partial }\mathit{\ell }}\cdot \frac{\mathrm{\partial }}{\mathrm{\partial }\left[{\mathbf{r}}_i+\mathit{\ell }\right]}\delta \left(\mathbf{r}-{\mathit{r}}_i-\mathit{\ell }\right)\mathrm{d}\mathit{\ell }\\ & =\frac{\mathrm{\partial }}{\mathrm{\partial }\mathbf{r}}\cdot \oint \delta \left(\mathbf{r}-{\mathit{r}}_i-\mathit{\ell }\right)\mathrm{d}\mathit{\ell },\end{array}$(6) ) can be considered to be $\mathit{\ell }={\mathit{\ell }}_t+{\mathit{\ell }}_n$ with ${\mathit{\ell }}_t$ and ${\mathit{\ell }}_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 the 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 [8, 9], which leads to the so-called Volume Average (VA) pressure formulation.

2.4. Volume average pressure

The volume average (VA) pressure is obtained from a volume integral of Equation (4(4) $\mathit{P}(\mathit{r},t)=\sum _{i=1}^N\frac{{\mathit{p}}_i{\mathit{p}}_i}{m_i}\delta ({\mathit{r}}_i-\mathit{r})+\frac{1}{2}\sum _{i,j}^N{\mathit{f}}_{ij}{\mathit{r}}_{ij}{O}_{ij}\delta ({\mathit{r}}_i-\mathit{r}),$(4) ), assuming linear inter-molecular interactions using Equations (6(6) $\begin{array}{rl}\delta (\mathbf{r}-{\mathbf{r}}_i)-\delta (\mathbf{r}-{\mathbf{r}}_j)& =\oint \frac{\mathrm{\partial }}{\mathrm{\partial }\mathit{\ell }}\cdot \delta \left(\mathbf{r}-{\mathit{r}}_i-\mathit{\ell }\right)\mathrm{d}\mathit{\ell }\\ & =\oint \frac{\mathrm{\partial }{\mathbf{r}}_i+\mathit{\ell }}{\mathrm{\partial }\mathit{\ell }}\cdot \frac{\mathrm{\partial }}{\mathrm{\partial }\left[{\mathbf{r}}_i+\mathit{\ell }\right]}\delta \left(\mathbf{r}-{\mathit{r}}_i-\mathit{\ell }\right)\mathrm{d}\mathit{\ell }\\ & =\frac{\mathrm{\partial }}{\mathrm{\partial }\mathbf{r}}\cdot \oint \delta \left(\mathbf{r}-{\mathit{r}}_i-\mathit{\ell }\right)\mathrm{d}\mathit{\ell },\end{array}$(6) ) and (7(7) $\frac{\mathrm{\partial }}{\mathrm{\partial }\mathbf{r}}\cdot \oint \delta \left(\mathbf{r}-{\mathit{r}}_i-\mathit{\ell }\right)\mathrm{d}\mathit{\ell }=\frac{\mathrm{\partial }}{\mathrm{\partial }\mathbf{r}}\cdot {\mathit{r}}_{ij}{\int }_0^1\delta \left(\mathbf{r}-{\mathit{r}}_i-\lambda {\mathit{r}}_{ij}\right)\mathrm{d}\lambda ,$(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, (15) $\stackrel{VA}{P}=\frac{1}{\mathrm{\Delta }V}\left[\sum _{i=1}^N\frac{{\mathbf{p}}_i{\mathbf{p}}_i}{m_i}{\vartheta }_i+\frac{1}{2}\sum _{i,j}^N{\mathbf{f}}_{ij}{\mathbf{r}}_{ij}\underset{0}{\overset{1}{\int }}{\vartheta }_{\lambda }\mathrm{d}\lambda \right],$(15) where $\mathrm{\Delta }V$ is the local volume; the ${\vartheta }_i$ function is only non-zero for a molecule inside the volume and the ${\vartheta }_{\lambda }$ is non-zero when a fraction of the interaction is inside the averaging volume. The $\vartheta$ function is simply the integral of the Dirac delta function over a volume [33, 38, 39], ${\vartheta }_{\alpha }={\mathrm{\Lambda }}_{x\alpha }{\mathrm{\Lambda }}_{y\alpha }{\mathrm{\Lambda }}_{z\alpha }$, where ${\mathrm{\Lambda }}_{x\alpha }=H(x^{+}-x_{\alpha })-H(x^{-}-x_{\alpha })$ with H the Heaviside function, $x^{+}$ the top of a volume in x, $x^{-}$ the bottom and $x_{\alpha }$ the position of either molecules when $\alpha =i$ or a point on the line of interaction when $\alpha =\lambda$ so $x_{\lambda }=x_i+\lambda x_{ij}$. The difference between the two Heaviside functions is known as a boxcar function, shown in Figure 1, a function which is one between two points and zero outside.

Figure 1. (Colour online) 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, ${\lambda }_k$, is checked to be between the two molecules using the difference $H(1-{\lambda }_k)-H(-{\lambda }_k),$ while the various Λ expressions 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 an arbitrary surface) or along the z-axis.

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 ${\mathrm{\Lambda }}_{z\alpha }=H(z^{+}+\xi (x_{\alpha },y_{\alpha },t)$ - $z_{\alpha })-H(z^{-}+\xi (x_{\alpha },$ $y_{\alpha },t)-z_{\alpha })$. Mathematically the value of $z_{\alpha }$, with molecule $\alpha =i$ or line $\alpha =\lambda$, can be interpreted as shifted based on the surface ξ at the same $x_{\alpha }$ and $y_{\alpha }$ 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 ${\mathrm{\Lambda }}_{\lambda x}$ and ${\mathrm{\Lambda }}_{\lambda y}$ functions give the line segment in the volume in x, which is $\mathrm{\Delta }x$, and in y, which is $\mathrm{\Delta }y$, and taking $P_T^C$ as the sum of the xx and yy components, (16) $P_T^C=\frac{1}{4}\sum _{i,j}^N\left[\frac{f_{xij}{x}_{ij}+f_{yij}{y}_{ij}}{\mathrm{\Delta }z}\right]{\mathrm{\Lambda }}_{iz}.$(16) 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 (17) $\stackrel{VA}{P}=\frac{1}{\mathrm{\Delta }V}\left[\sum _{i=1}^N\frac{{\mathbf{p}}_i{\mathbf{p}}_i}{m_i}{\vartheta }_i+\frac{1}{2}\sum _{i,j}^N{\mathbf{f}}_{ij}{\mathbf{r}}_{ij}\oint {\vartheta }_{\ell }\mathrm{d}\ell \right],$(17) where the length of a contour in a volume is obtained from the integral over ${\vartheta }_{\mathit{\ell }}={\mathrm{\Lambda }}_{x\ell }{\mathrm{\Lambda }}_{y\ell }{\mathrm{\Lambda }}_{z\ell }$ function with ${\mathrm{\Lambda }}_{x\ell }=H(x^{+}-x_i-{\ell }_x)-H(x^{-}-x_i-{\ell }_x),$ where ${\ell }_x$ is the x component of the contour, similarly for y with the moving z surface, ${\mathrm{\Lambda }}_{z\ell }=H(z^{+}+\xi ({\ell }_x,{\ell }_y,t)-z_i-{\ell }_z)$ - $H(z^{-}+\xi ({\ell }_x,{\ell }_y,t)-z_i-{\ell }_z)$. In practice, this would be obtained in a piecewise manner in a computer code, stepping along a contour and evaluating ${\vartheta }_{\mathit{\ell }}$ 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.

2.5. Method of Planes and Surface Flux Pressure

The surface pressure approach avoids the expansion in delta functions of Equation (3(3) $O_{ij}=\left(1-\frac{1}{2}{\mathit{r}}_{ij}\cdot \frac{\mathrm{\partial }}{\mathrm{\partial }{\mathit{r}}_i}+\cdots +\frac{1}{n!}{\left(-{\mathit{r}}_{ij}\cdot \frac{\mathrm{\partial }}{\mathrm{\partial }{\mathit{r}}_i}\right)}^{n-1}+\cdots \right).$(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, Equation (1(1) $\begin{array}{rl}& \mathit{P}(\mathit{r},t)=\\ & \sum _{i=1}^N〈m_i\left({\dot{\mathit{r}}}_i-\mathit{u}\right)\left({\dot{\mathit{r}}}_i-\mathit{u}\right)\delta ({\mathit{r}}_i-\mathit{r})+\frac{1}{2}\sum _{j\ne i}^N{\mathit{f}}_{ij}{\mathit{r}}_{ij}{O}_{ij}\delta ({\mathit{r}}_i-\mathit{r});f〉,\end{array}$(1) ) in [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 (18) $\begin{array}{rl}& \stackrel{MOP}{P}{}_x=\\ & \frac{1}{\mathrm{\Delta }{A}_x}\left[\sum _{i=1}^N\frac{{\mathbf{p}}_i{p}_{ix}}{m_i}\delta (x_i-x)+\frac{1}{4}\sum _{i,j}^N{\mathbf{f}}_{ij}\left[sgn(x-x_j)-sgn(x-x_i)\right]\right],\end{array}$(18) where $sgn\left(x\right)$ is the signum function. These two signum functions have the interpretation of including inter-molecular crossings over the surface of the plane. The normal component $P_{xx}$ uses the normal components $p_{ix}$ and $f_{xij}$, and is equivalent to the IK form of pressure in the normal direction given in Equation (11(11) $\stackrel{\mathrm{I}\mathrm{K}}{\mathrm{P}}\underset{N}{\overset{C}{}}=\stackrel{\mathrm{I}\mathrm{K}}{\mathrm{P}}\underset{zz}{\overset{C}{}}=\frac{1}{2}\sum _{i,j}^N{f}_{zij}\frac{z_{ij}}{|z_{ij}|}H\left(\frac{z-z_i}{z_{ij}}\right)H\left(\frac{z_j-z}{z_{ij}}\right),$(11) ) (as $sgn\left(z_{ij}\right)=z_{ij}/\left|z_{ij}\right|$ is cancelled by the $z_{ij}$ denominator in Equation (11(11) $\stackrel{\mathrm{I}\mathrm{K}}{\mathrm{P}}\underset{N}{\overset{C}{}}=\stackrel{\mathrm{I}\mathrm{K}}{\mathrm{P}}\underset{zz}{\overset{C}{}}=\frac{1}{2}\sum _{i,j}^N{f}_{zij}\frac{z_{ij}}{|z_{ij}|}H\left(\frac{z-z_i}{z_{ij}}\right)H\left(\frac{z_j-z}{z_{ij}}\right),$(11) ), see [38]). The other two components of momentum and force provide the two other pressure components, $P_{xy}\mathrm{a}\mathrm{n}\mathrm{d}{P}_{xz}$, and return 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 $P_{yy},P_{yx}\mathrm{a}\mathrm{n}\mathrm{d}{P}_{yz}$ as (19) $\begin{array}{rl}\stackrel{MOP}{P}{}_y& =\frac{1}{\mathrm{\Delta }{A}_x}\left[\sum _{i=1}^N\frac{{\mathbf{p}}_i{p}_{ix}}{m_i}\delta (y_i-y){\mathrm{\Lambda }}_{xi}{\mathrm{\Lambda }}_{zi}\right.\\ & \left.+\frac{1}{4}\sum _{i,j}^N{\mathbf{f}}_{ij}\left[sgn(y-y_j)-sgn(y-y_i)\right]{\mathrm{\Lambda }}_x\left({\lambda }_k\right){\mathrm{\Lambda }}_z\left({\lambda }_k\right)\right],\end{array}$(19) where ${\mathrm{\Lambda }}_{xi}$ and ${\mathrm{\Lambda }}_x\left({\lambda }_k\right)$ are boxcar functions which localise the plane to a patch in space. These check if the molecular i and point of crossing $\left({\lambda }_k\right)$ are between the limits of the surface in x (similar for y) with ${\mathrm{\Lambda }}_x\left({\lambda }_k\right)=H(x^{+}-x_i$ + ${\lambda }_k{x}_{ij})-H(x^{-}-x_i+{\lambda }_k{x}_{ij})$. It is instructive to compare this pressure, where $P_{yy}$ is the y component of force on the y surface, with the Irving Kirkwood contour of Equation (12(12) $\stackrel{\mathrm{I}\mathrm{K}}{\mathrm{P}}\underset{T}{\overset{C}{}}=\frac{1}{2}\left[\stackrel{\mathrm{I}\mathrm{K}}{\mathrm{P}}\underset{xx}{\overset{C}{}}+\stackrel{\mathrm{I}\mathrm{K}}{\mathrm{P}}\underset{yy}{\overset{C}{}}\right]=\frac{1}{4}\sum _{i,j}^N\left[\frac{f_{xij}{x}_{ij}+f_{yij}{y}_{ij}}{|z_{ij}|}\right]H\left(\frac{z-z_i}{z_{ij}}\right)H\left(\frac{z_j-z}{z_{ij}}\right).$(12) ) which measures the y force times $y_{ij}$ but on the z surface. This difference will be small for thin volumes but could become significant as $\mathrm{\Delta }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 $\mathbf{\nabla }\cdot \mathit{P}=0$ is satisfied, shear components, e.g. $\mathrm{\partial }{P}_{xy}/\mathrm{\partial }x,$ maybe 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 a 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 time step, so becomes $\mathbf{\nabla }\cdot \mathit{P}=\mathrm{d}\rho \mathit{u}/\mathrm{d}t+\mathbf{\nabla }\cdot \rho \mathit{u}\mathit{u}$ 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 $\mathrm{d}{S}_{xi}^{+}=\delta (x_i-x^{+}){\mathrm{\Lambda }}_{iy}{\mathrm{\Lambda }}_{iz}$ and $\mathrm{d}{S}_{x\lambda }^{+}=\delta (x_{\lambda }-x^{+}),{\mathrm{\Lambda }}_{\lambda y}{\mathrm{\Lambda }}_{\lambda z}$ for the top surface (denoted by the + superscript) in the x -direction, (20) $\rho \mathbf{u}{u}_x+\stackrel{SF}{P_x^{+}}=\frac{1}{\mathrm{\Delta }{S}_x}\sum _{i=1}^N{m}_i{\dot{r}}_i{\dot{x}}_i\mathrm{d}{S}_{xi}^{+}+\frac{1}{2\mathrm{\Delta }{S}_x}\sum _{i,j}^N{\mathbf{f}}_{ij}\underset{0}{\overset{1}{\int }}{x}_{ij}\mathrm{d}{S}_{x\lambda }^{+}\mathrm{d}\lambda ,$(20) and it can be shown that (21) ${\int }_0^1{x}_{ij}\mathrm{d}{S}_{x\lambda }^{+}\mathrm{d}\lambda =\left[sgn(x^{+}-x_j)-sgn(x^{+}-x_i)\right]{\mathrm{\Lambda }}_y\left({\lambda }_k\right){\mathrm{\Lambda }}_z\left({\lambda }_k\right)$(21) for a flat interface. The use of peculiar momentum ${\mathit{p}}_i$ has been dropped and the convective terms moved to the left-hand side to avoid making the assumption that streaming velocity is constant. Equation (20(20) $\rho \mathbf{u}{u}_x+\stackrel{SF}{P_x^{+}}=\frac{1}{\mathrm{\Delta }{S}_x}\sum _{i=1}^N{m}_i{\dot{r}}_i{\dot{x}}_i\mathrm{d}{S}_{xi}^{+}+\frac{1}{2\mathrm{\Delta }{S}_x}\sum _{i,j}^N{\mathbf{f}}_{ij}\underset{0}{\overset{1}{\int }}{x}_{ij}\mathrm{d}{S}_{x\lambda }^{+}\mathrm{d}\lambda ,$(20) ) can, therefore, be seen to be the localised MOP as shown in Equation (19(19) $\begin{array}{rl}\stackrel{MOP}{P}{}_y& =\frac{1}{\mathrm{\Delta }{A}_x}\left[\sum _{i=1}^N\frac{{\mathbf{p}}_i{p}_{ix}}{m_i}\delta (y_i-y){\mathrm{\Lambda }}_{xi}{\mathrm{\Lambda }}_{zi}\right.\\ & \left.+\frac{1}{4}\sum _{i,j}^N{\mathbf{f}}_{ij}\left[sgn(y-y_j)-sgn(y-y_i)\right]{\mathrm{\Lambda }}_x\left({\lambda }_k\right){\mathrm{\Lambda }}_z\left({\lambda }_k\right)\right],\end{array}$(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 ${\mathrm{\Lambda }}_{\alpha 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 Figure 1).

For interactions crossing the arbitrary surface itself, ξ, the SF pressure is expressed as follows: (22) $\begin{array}{rl}& \rho \mathbf{u}{u}_z+\stackrel{SF}{P_x^{+}}=\\ & \frac{1}{\mathrm{\Delta }{S}_z}\sum _{i=1}^N{m}_i{\dot{r}}_i\left[\stackrel{\text{Kinetic Curvature}}{\overbrace{{\dot{x}}_i\frac{\mathrm{\partial }{\xi }_i^{+}}{\mathrm{\partial }{x}_i}+{\dot{y}}_i\frac{\mathrm{\partial }{\xi }_i^{+}}{\mathrm{\partial }{y}_i}}}+{\dot{z}}_i+\stackrel{\text{Surface Evolution}}{\overbrace{\frac{\mathrm{\partial }{\xi }_i^{+}}{\mathrm{\partial }t}}}\right]\mathrm{d}{S}_{zi}^{+}\\ & +\frac{1}{2\mathrm{\Delta }{S}_z}\sum _{i,j}^N{\mathbf{f}}_{ij}{\int }_0^1\left[\underset{\text{Configurational Curvature}}{\underbrace{x_{ij}\frac{\mathrm{\partial }{\xi }_{\lambda }^{+}}{\mathrm{\partial }{x}_{\lambda }}+y_{ij}\frac{\mathrm{\partial }{\xi }_{\lambda }^{+}}{\mathrm{\partial }{y}_{\lambda }}}}+z_{ij}\right]\mathrm{d}{S}_{z\lambda }^{+}\mathrm{d}\lambda ,\end{array}$(22) where $\mathrm{d}{S}_{z\alpha }^{+}=\delta ({\xi }_{\alpha }^{+}-z_{\alpha }){\mathrm{\Lambda }}_x{\mathrm{\Lambda }}_y$ with ${\xi }_{\alpha }^{+}=z^{+}+\xi (x_{\alpha },y_{\alpha },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 the 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 Equation (22(22) $\begin{array}{rl}& \rho \mathbf{u}{u}_z+\stackrel{SF}{P_x^{+}}=\\ & \frac{1}{\mathrm{\Delta }{S}_z}\sum _{i=1}^N{m}_i{\dot{r}}_i\left[\stackrel{\text{Kinetic Curvature}}{\overbrace{{\dot{x}}_i\frac{\mathrm{\partial }{\xi }_i^{+}}{\mathrm{\partial }{x}_i}+{\dot{y}}_i\frac{\mathrm{\partial }{\xi }_i^{+}}{\mathrm{\partial }{y}_i}}}+{\dot{z}}_i+\stackrel{\text{Surface Evolution}}{\overbrace{\frac{\mathrm{\partial }{\xi }_i^{+}}{\mathrm{\partial }t}}}\right]\mathrm{d}{S}_{zi}^{+}\\ & +\frac{1}{2\mathrm{\Delta }{S}_z}\sum _{i,j}^N{\mathbf{f}}_{ij}{\int }_0^1\left[\underset{\text{Configurational Curvature}}{\underbrace{x_{ij}\frac{\mathrm{\partial }{\xi }_{\lambda }^{+}}{\mathrm{\partial }{x}_{\lambda }}+y_{ij}\frac{\mathrm{\partial }{\xi }_{\lambda }^{+}}{\mathrm{\partial }{y}_{\lambda }}}}+z_{ij}\right]\mathrm{d}{S}_{z\lambda }^{+}\mathrm{d}\lambda ,\end{array}$(22) ) will simplify to the localised form of the MOP pressure given in Equation (20(20) $\rho \mathbf{u}{u}_x+\stackrel{SF}{P_x^{+}}=\frac{1}{\mathrm{\Delta }{S}_x}\sum _{i=1}^N{m}_i{\dot{r}}_i{\dot{x}}_i\mathrm{d}{S}_{xi}^{+}+\frac{1}{2\mathrm{\Delta }{S}_x}\sum _{i,j}^N{\mathbf{f}}_{ij}\underset{0}{\overset{1}{\int }}{x}_{ij}\mathrm{d}{S}_{x\lambda }^{+}\mathrm{d}\lambda ,$(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, (23) $\begin{array}{rl}& {\int }_{t_1}^{t_2}\rho \mathbf{u}{u}_z+\stackrel{SF}{P_z^{c+}}\mathrm{d}t=\frac{1}{\mathrm{\Delta }{S}_z}\sum _{i=1}^N{m}_i{\dot{r}}_i\frac{{\mathbf{r}}_{12}\cdot {\mathbf{n}}_z}{|{\mathbf{r}}_{12}\cdot {\mathbf{n}}_z|}\mathrm{d}{S}^{+}\\ & +\frac{1}{\mathrm{\Delta }{S}_z}\sum _{i=1}^N{m}_i{\dot{r}}_i{\vartheta }_t\stackrel{SF}{P_z^{c+}}=\frac{1}{2\mathrm{\Delta }{S}_z}\sum _{i,j}^N{\mathbf{f}}_{ij}\frac{{\mathbf{r}}_{ij}\cdot {\mathbf{n}}_z}{|{\mathbf{r}}_{ij}\cdot {\mathbf{n}}_z|}\mathrm{d}{S}^{+},\end{array}$(23) where ${\mathit{r}}_{12}={\mathit{r}}_2-{\mathit{r}}_1$ is the line of time evolution of a molecule between t and $t+\mathrm{\Delta }t$, ${\vartheta }_t=\left[H\right(\xi \left(t\right)-z_i\left(t\right))-H(\xi (t-\mathrm{\Delta }t)-z_i\left(t\right)\left)\right]{\mathrm{\Lambda }}_{ix}{\mathrm{\Lambda }}_{iy}$ counting how many molecules have left or entered the volume due to surface movement between time $t-\mathrm{\Delta }t$ and t, and the surface normal includes all of the surface curvature terms of Equation (22(22) $\begin{array}{rl}& \rho \mathbf{u}{u}_z+\stackrel{SF}{P_x^{+}}=\\ & \frac{1}{\mathrm{\Delta }{S}_z}\sum _{i=1}^N{m}_i{\dot{r}}_i\left[\stackrel{\text{Kinetic Curvature}}{\overbrace{{\dot{x}}_i\frac{\mathrm{\partial }{\xi }_i^{+}}{\mathrm{\partial }{x}_i}+{\dot{y}}_i\frac{\mathrm{\partial }{\xi }_i^{+}}{\mathrm{\partial }{y}_i}}}+{\dot{z}}_i+\stackrel{\text{Surface Evolution}}{\overbrace{\frac{\mathrm{\partial }{\xi }_i^{+}}{\mathrm{\partial }t}}}\right]\mathrm{d}{S}_{zi}^{+}\\ & +\frac{1}{2\mathrm{\Delta }{S}_z}\sum _{i,j}^N{\mathbf{f}}_{ij}{\int }_0^1\left[\underset{\text{Configurational Curvature}}{\underbrace{x_{ij}\frac{\mathrm{\partial }{\xi }_{\lambda }^{+}}{\mathrm{\partial }{x}_{\lambda }}+y_{ij}\frac{\mathrm{\partial }{\xi }_{\lambda }^{+}}{\mathrm{\partial }{y}_{\lambda }}}}+z_{ij}\right]\mathrm{d}{S}_{z\lambda }^{+}\mathrm{d}\lambda ,\end{array}$(22) ) with ${\mathit{n}}_z=\frac{{\mathbf{\nabla }}_{\alpha }(\xi -z_{\alpha })}{\parallel {\mathbf{\nabla }}_{\alpha }(\xi -z_{\alpha })\parallel }$. This normal is the one used in the Harasima contour for a general surface, which could be combined with the tangential contour of Equation (16(16) $P_T^C=\frac{1}{4}\sum _{i,j}^N\left[\frac{f_{xij}{x}_{ij}+f_{yij}{y}_{ij}}{\mathrm{\Delta }z}\right]{\mathrm{\Lambda }}_{iz}.$(16) ). We have generalised the treatment of both kinetic and configurational terms in Equation (23(23) $\begin{array}{rl}& {\int }_{t_1}^{t_2}\rho \mathbf{u}{u}_z+\stackrel{SF}{P_z^{c+}}\mathrm{d}t=\frac{1}{\mathrm{\Delta }{S}_z}\sum _{i=1}^N{m}_i{\dot{r}}_i\frac{{\mathbf{r}}_{12}\cdot {\mathbf{n}}_z}{|{\mathbf{r}}_{12}\cdot {\mathbf{n}}_z|}\mathrm{d}{S}^{+}\\ & +\frac{1}{\mathrm{\Delta }{S}_z}\sum _{i=1}^N{m}_i{\dot{r}}_i{\vartheta }_t\stackrel{SF}{P_z^{c+}}=\frac{1}{2\mathrm{\Delta }{S}_z}\sum _{i,j}^N{\mathbf{f}}_{ij}\frac{{\mathbf{r}}_{ij}\cdot {\mathbf{n}}_z}{|{\mathbf{r}}_{ij}\cdot {\mathbf{n}}_z|}\mathrm{d}{S}^{+},\end{array}$(23) ) by recognising that molecules must move from a starting point ${\mathit{r}}_1={\mathit{r}}_i\left(t\right)$ to an endpoint ${\mathit{r}}_2={\mathit{r}}_i(t+\mathrm{\Delta }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, ${\lambda }_k$, is obtained, it can be inserted in the following expression (24) $\mathrm{d}{S}^{+}=\sum _{k=1}^{N_{roots}}\left[H\left(1-{\lambda }_k\right)-H\left(-{\lambda }_k\right)\right]{\mathrm{\Lambda }}_x\left({\lambda }_k\right){\mathrm{\Lambda }}_y\left({\lambda }_k\right),$(24) with the first two Heaviside functions selecting if the root ${\lambda }_k$ is between the endpoints ${\mathit{r}}_1$ and ${\mathit{r}}_2$. 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, this approach should be applicable for any interaction contour and surface, simply requiring the point on the contour that the surface crossing ${\lambda }_k$ occurs to evaluate $\mathrm{d}{S}^{+}$. It seems reasonable that this approach could therefore be applied to get the intersection of Equation (6(6) $\begin{array}{rl}\delta (\mathbf{r}-{\mathbf{r}}_i)-\delta (\mathbf{r}-{\mathbf{r}}_j)& =\oint \frac{\mathrm{\partial }}{\mathrm{\partial }\mathit{\ell }}\cdot \delta \left(\mathbf{r}-{\mathit{r}}_i-\mathit{\ell }\right)\mathrm{d}\mathit{\ell }\\ & =\oint \frac{\mathrm{\partial }{\mathbf{r}}_i+\mathit{\ell }}{\mathrm{\partial }\mathit{\ell }}\cdot \frac{\mathrm{\partial }}{\mathrm{\partial }\left[{\mathbf{r}}_i+\mathit{\ell }\right]}\delta \left(\mathbf{r}-{\mathit{r}}_i-\mathit{\ell }\right)\mathrm{d}\mathit{\ell }\\ & =\frac{\mathrm{\partial }}{\mathrm{\partial }\mathbf{r}}\cdot \oint \delta \left(\mathbf{r}-{\mathit{r}}_i-\mathit{\ell }\right)\mathrm{d}\mathit{\ell },\end{array}$(6) ), with any surface. The pressures tensors defined from a closed volume of bounding surfaces would have the property already discussed which ensures any missed contributions from differing contour paths would be included as shear components on another surface.

2.6. 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. VA pressure equation (15(15) $\stackrel{VA}{P}=\frac{1}{\mathrm{\Delta }V}\left[\sum _{i=1}^N\frac{{\mathbf{p}}_i{\mathbf{p}}_i}{m_i}{\vartheta }_i+\frac{1}{2}\sum _{i,j}^N{\mathbf{f}}_{ij}{\mathbf{r}}_{ij}\underset{0}{\overset{1}{\int }}{\vartheta }_{\lambda }\mathrm{d}\lambda \right],$(15) ) is commonly obtained by integrating Equation (4(4) $\mathit{P}(\mathit{r},t)=\sum _{i=1}^N\frac{{\mathit{p}}_i{\mathit{p}}_i}{m_i}\delta ({\mathit{r}}_i-\mathit{r})+\frac{1}{2}\sum _{i,j}^N{\mathit{f}}_{ij}{\mathit{r}}_{ij}{O}_{ij}\delta ({\mathit{r}}_i-\mathit{r}),$(4) ) over a volume [8, 33]. However, we derive the SF equations from the time evolution of momentum inside a control volume, (25) $\frac{\mathrm{d}}{\mathrm{d}t}\sum _{i=1}^N{m}_i{\dot{r}}_i{\vartheta }_i=\sum _{i=1}^N{m}_i{\dot{r}}_i\frac{\mathrm{d}{\vartheta }_i}{\mathrm{d}t}+\sum _{i,j}^N{\mathit{f}}_{ij}{\mathit{r}}_{ij}\cdot {\int }_0^1\frac{\mathrm{\partial }{\vartheta }_{\lambda }}{\mathrm{\partial }{\mathit{r}}_{\lambda }}\mathrm{d}\lambda ,$(25) which results in extra terms for both the time evolution of the moving volume and the surface curvature. For the first term, (26) $\begin{array}{rl}& \sum _{i=1}^N{m}_i{\dot{r}}_i\frac{\mathrm{d}{\vartheta }_i}{\mathrm{d}t}=\\ & \frac{\mathrm{\partial }}{\mathrm{\partial }\mathbf{r}}\cdot \stackrel{\stackrel{VA}{P}k\mathrm{\Delta }V}{\overbrace{\sum _{i=1}^N{m}_i{\dot{r}}_i{\dot{r}}_i{\vartheta }_i}}+\stackrel{\text{Extra Kinetic Curvature }\left(\mathrm{K}\mathrm{C}\right)\text{ Terms}}{\overbrace{\sum _{i=1}^N{m}_i{\dot{r}}_i{\dot{r}}_i\cdot \left[\frac{\mathrm{\partial }{\xi }_i^{+}}{\mathrm{\partial }{\mathbf{r}}_i}\mathrm{d}{S}_{zi}^{+}-\frac{\mathrm{\partial }{\xi }_i^{-}}{\mathrm{\partial }{\mathbf{r}}_i}\mathrm{d}{S}_{zi}^{-}\right]}}\\ & +\underset{\text{Extra Time Evolving }\left(\mathrm{T}\mathrm{E}\right)\text{ Terms}}{\underbrace{\sum _{i=1}^N{m}_i{\dot{r}}_i\left[\frac{\mathrm{\partial }{\xi }_i^{+}}{\mathrm{\partial }t}\mathrm{d}{S}_{zi}^{+}-\frac{\mathrm{\partial }{\xi }_i^{-}}{\mathrm{\partial }t}\mathrm{d}{S}_{zi}^{-}\right]}}\end{array}$(26) where we have used $\mathrm{\partial }{\xi }_i^{\pm }/\mathrm{\partial }{z}_i=0$ to write the expression concisely in a vector form. We identify the kinetic part with reference to Equation (15(15) $\stackrel{VA}{P}=\frac{1}{\mathrm{\Delta }V}\left[\sum _{i=1}^N\frac{{\mathbf{p}}_i{\mathbf{p}}_i}{m_i}{\vartheta }_i+\frac{1}{2}\sum _{i,j}^N{\mathbf{f}}_{ij}{\mathbf{r}}_{ij}\underset{0}{\overset{1}{\int }}{\vartheta }_{\lambda }\mathrm{d}\lambda \right],$(15) ) using an overbrace, under the assumption that everything inside the $\mathrm{\partial }/\mathrm{\partial }\mathit{r}$ operator is VA pressure, and identify the two additional terms, KC and TE, missing from the VA treatment.

The second term in Equation (25(25) $\frac{\mathrm{d}}{\mathrm{d}t}\sum _{i=1}^N{m}_i{\dot{r}}_i{\vartheta }_i=\sum _{i=1}^N{m}_i{\dot{r}}_i\frac{\mathrm{d}{\vartheta }_i}{\mathrm{d}t}+\sum _{i,j}^N{\mathit{f}}_{ij}{\mathit{r}}_{ij}\cdot {\int }_0^1\frac{\mathrm{\partial }{\vartheta }_{\lambda }}{\mathrm{\partial }{\mathit{r}}_{\lambda }}\mathrm{d}\lambda ,$(25) ) is the configurational part, which for the VA in a moving interface, swapping $\mathrm{\partial }{\vartheta }_{\lambda }/\mathrm{\partial }{r}_{\lambda }$ to $\mathrm{\partial }{\vartheta }_{\lambda }/\mathrm{\partial }r$ results in additional terms [34]; so the full expression with correction term CC is (27) $\begin{array}{rl}\sum _{i,j}^N{\mathbf{f}}_{ij}{\mathbf{r}}_{ij}\cdot \underset{0}{\overset{1}{\int }}\frac{\mathrm{\partial }{\vartheta }_{\lambda }}{\mathrm{\partial }{\mathbf{r}}_{\lambda }}\mathrm{d}\lambda & =-\frac{\mathrm{\partial }}{\mathrm{\partial }\mathbf{r}}\cdot \stackrel{\stackrel{VA}{P^C}\mathrm{\Delta }V}{\overbrace{\sum _{i,j}^N{\mathbf{f}}_{ij}{\mathbf{r}}_{ij}\underset{0}{\overset{1}{\int }}{\vartheta }_{\lambda }\mathrm{d}\lambda }}\\ & +\underset{\text{Extra Configurational Curvature }\left(\mathrm{C}\mathrm{C}\right)\text{ Term}}{\underbrace{\sum _{i,j}^N{\mathbf{f}}_{ij}{\mathbf{r}}_{ij}\cdot \underset{0}{\overset{1}{\int }}\left[\frac{\mathrm{\partial }{\xi }_{\lambda }^{+}}{\mathrm{\partial }{\mathbf{r}}_{\lambda }}\mathrm{d}{S}_{z\lambda }^{+}-\frac{\mathrm{\partial }{\xi }_{\lambda }^{-}}{\mathrm{\partial }{\mathbf{r}}_{\lambda }}\mathrm{d}{S}_{z\lambda }^{-}\right]\mathrm{d}\lambda }},\end{array}$(27) again using $\mathrm{\partial }{\xi }_{\lambda }^{\pm }/\mathrm{\partial }{z}_{\lambda }=0$ to write in a vector form. So, Equation (25(25) $\frac{\mathrm{d}}{\mathrm{d}t}\sum _{i=1}^N{m}_i{\dot{r}}_i{\vartheta }_i=\sum _{i=1}^N{m}_i{\dot{r}}_i\frac{\mathrm{d}{\vartheta }_i}{\mathrm{d}t}+\sum _{i,j}^N{\mathit{f}}_{ij}{\mathit{r}}_{ij}\cdot {\int }_0^1\frac{\mathrm{\partial }{\vartheta }_{\lambda }}{\mathrm{\partial }{\mathit{r}}_{\lambda }}\mathrm{d}\lambda ,$(25) ) shows the link between SF and VA form with the extra terms naturally included in the SF derivation (28) $\frac{\mathrm{d}}{\mathrm{d}t}\sum _{i=1}^N{m}_i{\dot{r}}_i{\vartheta }_i=-\sum _{\alpha \in \{x,y,z\}}[\stackrel{SF}{P_{\alpha }^{+}}-\stackrel{SF}{P_{\alpha }^{-}}]\mathrm{\Delta }{S}_{\alpha }\approx -\frac{\mathrm{\partial }}{\mathrm{\partial }\mathbf{r}}\cdot \stackrel{VA}{P}\mathrm{\Delta }V+KC+TE+CC,$(28) 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 .¹

Pressure in the VA form is the average inside a divergence operator $\mathrm{\partial }/\mathrm{\partial }\mathit{r}$, so it is natural to write the extra terms in Equations (26(26) $\begin{array}{rl}& \sum _{i=1}^N{m}_i{\dot{r}}_i\frac{\mathrm{d}{\vartheta }_i}{\mathrm{d}t}=\\ & \frac{\mathrm{\partial }}{\mathrm{\partial }\mathbf{r}}\cdot \stackrel{\stackrel{VA}{P}k\mathrm{\Delta }V}{\overbrace{\sum _{i=1}^N{m}_i{\dot{r}}_i{\dot{r}}_i{\vartheta }_i}}+\stackrel{\text{Extra Kinetic Curvature }\left(\mathrm{K}\mathrm{C}\right)\text{ Terms}}{\overbrace{\sum _{i=1}^N{m}_i{\dot{r}}_i{\dot{r}}_i\cdot \left[\frac{\mathrm{\partial }{\xi }_i^{+}}{\mathrm{\partial }{\mathbf{r}}_i}\mathrm{d}{S}_{zi}^{+}-\frac{\mathrm{\partial }{\xi }_i^{-}}{\mathrm{\partial }{\mathbf{r}}_i}\mathrm{d}{S}_{zi}^{-}\right]}}\\ & +\underset{\text{Extra Time Evolving }\left(\mathrm{T}\mathrm{E}\right)\text{ Terms}}{\underbrace{\sum _{i=1}^N{m}_i{\dot{r}}_i\left[\frac{\mathrm{\partial }{\xi }_i^{+}}{\mathrm{\partial }t}\mathrm{d}{S}_{zi}^{+}-\frac{\mathrm{\partial }{\xi }_i^{-}}{\mathrm{\partial }t}\mathrm{d}{S}_{zi}^{-}\right]}}\end{array}$(26) ) and (27(27) $\begin{array}{rl}\sum _{i,j}^N{\mathbf{f}}_{ij}{\mathbf{r}}_{ij}\cdot \underset{0}{\overset{1}{\int }}\frac{\mathrm{\partial }{\vartheta }_{\lambda }}{\mathrm{\partial }{\mathbf{r}}_{\lambda }}\mathrm{d}\lambda & =-\frac{\mathrm{\partial }}{\mathrm{\partial }\mathbf{r}}\cdot \stackrel{\stackrel{VA}{P^C}\mathrm{\Delta }V}{\overbrace{\sum _{i,j}^N{\mathbf{f}}_{ij}{\mathbf{r}}_{ij}\underset{0}{\overset{1}{\int }}{\vartheta }_{\lambda }\mathrm{d}\lambda }}\\ & +\underset{\text{Extra Configurational Curvature }\left(\mathrm{C}\mathrm{C}\right)\text{ Term}}{\underbrace{\sum _{i,j}^N{\mathbf{f}}_{ij}{\mathbf{r}}_{ij}\cdot \underset{0}{\overset{1}{\int }}\left[\frac{\mathrm{\partial }{\xi }_{\lambda }^{+}}{\mathrm{\partial }{\mathbf{r}}_{\lambda }}\mathrm{d}{S}_{z\lambda }^{+}-\frac{\mathrm{\partial }{\xi }_{\lambda }^{-}}{\mathrm{\partial }{\mathbf{r}}_{\lambda }}\mathrm{d}{S}_{z\lambda }^{-}\right]\mathrm{d}\lambda }},\end{array}$(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 (29) $\begin{array}{rl}& \frac{\mathrm{\partial }}{\mathrm{\partial }z}\int \left[KC+TE+CC\right]\mathrm{d}z=\\ & \frac{\mathrm{\partial }}{\mathrm{\partial }z}\left[\sum _{i=1}^N\left(m_i{\dot{z}}_i{\dot{r}}_i\cdot \frac{\mathrm{\partial }{\xi }_i}{\mathrm{\partial }{\mathit{r}}_i}+m_i{\dot{z}}_i\frac{\mathrm{\partial }{\xi }_i}{\mathrm{\partial }t}\right){\vartheta }_i+\sum _{i,j}^N{f}_{zij}{\mathit{r}}_{ij}\cdot {\int }_0^1\frac{\mathrm{\partial }{\xi }_{\lambda }}{\mathrm{\partial }{\mathit{r}}_{\lambda }}{\vartheta }_{\lambda }\mathrm{d}\lambda \right]\end{array}$(29) so everything inside the derivative is in the form of a pressure contribution inside a volume. As with other VA terms, the functions ${\vartheta }_{\alpha }$ assigns pressure contributions only when molecules, $\alpha =i$, or points on interactions path, $\alpha =\lambda$, are inside the volume. The interface movement $\mathrm{\partial }{\xi }_i/\mathrm{\partial }t$ and surface curvature $\mathrm{\partial }{\xi }_i/\mathrm{\partial }{\mathit{r}}_i$ at the location of the molecules is used in the kinetic part, while the curvature integrated along the contour between molecules $\mathrm{\partial }{\xi }_{\lambda }^{}/\mathrm{\partial }{\mathit{r}}_{\lambda }$ 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 Equation (22(22) $\begin{array}{rl}& \rho \mathbf{u}{u}_z+\stackrel{SF}{P_x^{+}}=\\ & \frac{1}{\mathrm{\Delta }{S}_z}\sum _{i=1}^N{m}_i{\dot{r}}_i\left[\stackrel{\text{Kinetic Curvature}}{\overbrace{{\dot{x}}_i\frac{\mathrm{\partial }{\xi }_i^{+}}{\mathrm{\partial }{x}_i}+{\dot{y}}_i\frac{\mathrm{\partial }{\xi }_i^{+}}{\mathrm{\partial }{y}_i}}}+{\dot{z}}_i+\stackrel{\text{Surface Evolution}}{\overbrace{\frac{\mathrm{\partial }{\xi }_i^{+}}{\mathrm{\partial }t}}}\right]\mathrm{d}{S}_{zi}^{+}\\ & +\frac{1}{2\mathrm{\Delta }{S}_z}\sum _{i,j}^N{\mathbf{f}}_{ij}{\int }_0^1\left[\underset{\text{Configurational Curvature}}{\underbrace{x_{ij}\frac{\mathrm{\partial }{\xi }_{\lambda }^{+}}{\mathrm{\partial }{x}_{\lambda }}+y_{ij}\frac{\mathrm{\partial }{\xi }_{\lambda }^{+}}{\mathrm{\partial }{y}_{\lambda }}}}+z_{ij}\right]\mathrm{d}{S}_{z\lambda }^{+}\mathrm{d}\lambda ,\end{array}$(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.

3. 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, (30) $\varphi \left(r_{ij}\right)=ϵ\left[{\left(\frac{\sigma }{r_{ij}}\right)}^{12}-{\left(\frac{\sigma }{r_{ij}}\right)}^6\right]-\varphi \left(r_c\right);r_{ij}<r_c,$(30) with a cut-off $r_c=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 a good experimental agreement [40]. Time integration uses the velocity Verlet scheme (31) $\begin{array}{rl}{\mathit{r}}_i(t+\mathrm{\Delta }t)& ={\mathit{r}}_i\left(t\right)+\mathrm{\Delta }t{\mathit{v}}_i(t+\mathrm{\Delta }t/2),\\ {\mathit{v}}_i(t+\mathrm{\Delta }t/2)& ={\mathit{v}}_i(t-\mathrm{\Delta }t/2)+\mathrm{\Delta }t{\mathit{F}}_i\left(t\right)\end{array}$(31) with a time step $\mathrm{\Delta }t=0.005$ and force on i calculated from ${\mathit{F}}_i=\sum _{j\ne i}^N{\mathit{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].

Figure 2. (Colour online) A schematic of the three cases considered in this work: (a) a solid-liquid interface with a fixed grid, (b) a liquid–vapour interface with a fixed grid and (c) a liquid–vapour interface with a grid moving with the intrinsic interface. The profiles of normal pressure are shown below with kinetic (

), configurational (

) and total (

) pressure contributions.

For the solid–liquid simulation, the domain is set up 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 $k_4=5\times {10}^3$ and $k_6=5\times {10}^6$. The system is initialised using an FCC lattice of density $\rho =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 $T_{ave}=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\text{\%}$ designated as liquid with a density of ${\rho }_l=0.79$ and remaining gas at density ${\rho }_g=0.0002$ giving N=2635 molecules. The system is equilibrated for 100,000 time steps in a Nosè Hoover NVT ensemble with set point $T_s=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 cluster analysis with Stillinger cut-off length $r_d=1.5$ and requiring more than three neighbours per molecule. The functional form of surface $\xi (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 (32) $\xi (x,y,t)=\sum _{\mu =-k_u}^{k_u}\sum _{\nu =-k_u}^{k_u}{a}_{\mu \nu }\left(t\right)f_{\mu }\left(x\right)f_{\nu }\left(y\right),$(32) where $f_{\mu }\left(x\right)=\cos \left(k_{x}x\right)$, $f_{-\mu }\left(x\right)=\sin \left(k_{x}x\right)$ and $a_{\mu \nu }$ is the matrix of surface wavenumbers. The number of wavelengths is calculated from the system size $k_u=nint\left(\sqrt{LxLy}/{\lambda }_u\right)$ with ${\lambda }_u=1$ the minimum wavelength. This function is fitted by minimising the least-square difference between surface molecules $z_p$ and the intrinsic surface function at these positions $\xi =\xi (x_p,y_p,t)$, (33) $W=\frac{1}{2}\sum _{p=1}^{N_p}{\left[z_p\left(t\right)-\xi (x_p,y_p,t)\right]}^2+\psi \tilde{A},$(33) with extra constraint $\psi =1\times {10}^{-8}$ to prevent overfitting by ensuring intrinsic area $\tilde{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 $n_s=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 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 $\mathrm{\Delta }x=\mathrm{\Delta }y=0.198$ tangential to the surface and $\mathrm{\Delta }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 $\mathrm{\Delta }x$ by $\mathrm{\Delta }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 tested every time step, by summing up Equation (23(23) $\begin{array}{rl}& {\int }_{t_1}^{t_2}\rho \mathbf{u}{u}_z+\stackrel{SF}{P_z^{c+}}\mathrm{d}t=\frac{1}{\mathrm{\Delta }{S}_z}\sum _{i=1}^N{m}_i{\dot{r}}_i\frac{{\mathbf{r}}_{12}\cdot {\mathbf{n}}_z}{|{\mathbf{r}}_{12}\cdot {\mathbf{n}}_z|}\mathrm{d}{S}^{+}\\ & +\frac{1}{\mathrm{\Delta }{S}_z}\sum _{i=1}^N{m}_i{\dot{r}}_i{\vartheta }_t\stackrel{SF}{P_z^{c+}}=\frac{1}{2\mathrm{\Delta }{S}_z}\sum _{i,j}^N{\mathbf{f}}_{ij}\frac{{\mathbf{r}}_{ij}\cdot {\mathbf{n}}_z}{|{\mathbf{r}}_{ij}\cdot {\mathbf{n}}_z|}\mathrm{d}{S}^{+},\end{array}$(23) ) for all six surfaces and checking it has a difference of less than $1\times {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]

4. Results and discussion

In this section, we compare three types of pressure measurements, IK1 equation (5(5) $\stackrel{\mathrm{I}\mathrm{K}1}{\mathrm{P}}(\mathbf{r},t)=\sum _{i=1}^N\frac{{\mathbf{p}}_i{\mathbf{p}}_i}{m_i}\delta \left(\mathbf{r}-{\mathbf{r}}_i\right)+\frac{1}{2}\sum _{i,j}^N{\mathbf{f}}_{ij}{\mathbf{r}}_{ij}\delta \left(\mathbf{r}-{\mathbf{r}}_i\right).$(5) ), VA equation (15(15) $\stackrel{VA}{P}=\frac{1}{\mathrm{\Delta }V}\left[\sum _{i=1}^N\frac{{\mathbf{p}}_i{\mathbf{p}}_i}{m_i}{\vartheta }_i+\frac{1}{2}\sum _{i,j}^N{\mathbf{f}}_{ij}{\mathbf{r}}_{ij}\underset{0}{\overset{1}{\int }}{\vartheta }_{\lambda }\mathrm{d}\lambda \right],$(15) ) and SF equation (22(22) $\begin{array}{rl}& \rho \mathbf{u}{u}_z+\stackrel{SF}{P_x^{+}}=\\ & \frac{1}{\mathrm{\Delta }{S}_z}\sum _{i=1}^N{m}_i{\dot{r}}_i\left[\stackrel{\text{Kinetic Curvature}}{\overbrace{{\dot{x}}_i\frac{\mathrm{\partial }{\xi }_i^{+}}{\mathrm{\partial }{x}_i}+{\dot{y}}_i\frac{\mathrm{\partial }{\xi }_i^{+}}{\mathrm{\partial }{y}_i}}}+{\dot{z}}_i+\stackrel{\text{Surface Evolution}}{\overbrace{\frac{\mathrm{\partial }{\xi }_i^{+}}{\mathrm{\partial }t}}}\right]\mathrm{d}{S}_{zi}^{+}\\ & +\frac{1}{2\mathrm{\Delta }{S}_z}\sum _{i,j}^N{\mathbf{f}}_{ij}{\int }_0^1\left[\underset{\text{Configurational Curvature}}{\underbrace{x_{ij}\frac{\mathrm{\partial }{\xi }_{\lambda }^{+}}{\mathrm{\partial }{x}_{\lambda }}+y_{ij}\frac{\mathrm{\partial }{\xi }_{\lambda }^{+}}{\mathrm{\partial }{y}_{\lambda }}}}+z_{ij}\right]\mathrm{d}{S}_{z\lambda }^{+}\mathrm{d}\lambda ,\end{array}$(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 (colour online, see caption of Figure 3 for black and white equivalents in print). 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 2(a), 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 2(b) 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.

Figure 3. (Colour online) Comparing wall-normal, $P_N$, pressure measurements near a solid–liquid interface showing (a) half channel and (b) near-wall region. The kinetic components are shown for IK1 (

), VA (

) and SF (

), the configurational part for IK1 (

), VA (

) and SF (

) while the total pressure is IK1 (

), VA (

) and SF (

). The density (

) with the zero axis shown by a horizontal black line (

) and the shaded region on (a) is the section shown in (b).

Figure 4. (Colour online) (a) Normal $P_N$ and (b) tangential $P_T$ pressure near a liquid–vapour interface using a fixed reference frame (a uniform grid). The kinetic components for IK1 (

), VA (

) and SF (

), the configurational part for IK1 (

), VA (

) and SF (

) and the total pressure for IK1 (

), VA (

) and SF (

) with the zero axis shown by a horizontal black line (

). The Irving Kirkwood contour would give the same normal and tangential pressure as the VA and SF curves.

Finally, Figure 2(c) 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 Equation (22(22) $\begin{array}{rl}& \rho \mathbf{u}{u}_z+\stackrel{SF}{P_x^{+}}=\\ & \frac{1}{\mathrm{\Delta }{S}_z}\sum _{i=1}^N{m}_i{\dot{r}}_i\left[\stackrel{\text{Kinetic Curvature}}{\overbrace{{\dot{x}}_i\frac{\mathrm{\partial }{\xi }_i^{+}}{\mathrm{\partial }{x}_i}+{\dot{y}}_i\frac{\mathrm{\partial }{\xi }_i^{+}}{\mathrm{\partial }{y}_i}}}+{\dot{z}}_i+\stackrel{\text{Surface Evolution}}{\overbrace{\frac{\mathrm{\partial }{\xi }_i^{+}}{\mathrm{\partial }t}}}\right]\mathrm{d}{S}_{zi}^{+}\\ & +\frac{1}{2\mathrm{\Delta }{S}_z}\sum _{i,j}^N{\mathbf{f}}_{ij}{\int }_0^1\left[\underset{\text{Configurational Curvature}}{\underbrace{x_{ij}\frac{\mathrm{\partial }{\xi }_{\lambda }^{+}}{\mathrm{\partial }{x}_{\lambda }}+y_{ij}\frac{\mathrm{\partial }{\xi }_{\lambda }^{+}}{\mathrm{\partial }{y}_{\lambda }}}}+z_{ij}\right]\mathrm{d}{S}_{z\lambda }^{+}\mathrm{d}\lambda ,\end{array}$(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 is 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.

Figure 5. (Colour online) (a) Normal $P_N$ and (b) tangential $P_T$ pressure for a reference moving with the liquid–vapour interface, showing the kinetic components for IK1 (

), VA (

) and SF (

) (note the dashed line exactly follows the solid line) and the configurational part for IK1 (

), VA (

) and SF (

) with density (

) with the zero axis shown by a horizontal black line (

). The Harasima contour would give the same normal pressure as the SF curves in (a) and the same tangential pressure as the IK1 in (b).

Figure 6. (Colour online) The total pressure components for a reference moving with the interface: (a) normal $P_N$ and (b) tangential $P_T$, with IK1 (

), VA (

) and SF (

) with the kinetic and configurational constituent parts from Figure 5 shown faintly for reference and zero axis shown by a horizontal black line (

). The Harasima contour would give the same normal pressure as the SF curves in (a) and same tangential pressure as the IK1 in (b).

Figure 7. (Colour online) The cumulative integral of the total pressures in a reference moving with the interface shown in Figure 6 with normal $P_N$ on (a) and the negative of the tangential $-P_T$ pressure shown on (b), starting from z=−5 up to the current z on the axis, for the IK1 (

), VA (

) and SF (

) measurements, and the horizontal black line (

) is the zero axis. The inserts zoom in on the values that the integrals have converged to by the upper limit displayed on the plot, z=2, where the sum of normal and tangential components would give the total surface tension.

Figure 8. (Colour online) The terms required to balance normal pressure in a reference moving with the interface, where (a) shows the normal component of the Volume Average (VA) pressure, both kinetic (

) and configurational (

) with additional correction terms from Equations (26(26) $\begin{array}{rl}& \sum _{i=1}^N{m}_i{\dot{r}}_i\frac{\mathrm{d}{\vartheta }_i}{\mathrm{d}t}=\\ & \frac{\mathrm{\partial }}{\mathrm{\partial }\mathbf{r}}\cdot \stackrel{\stackrel{VA}{P}k\mathrm{\Delta }V}{\overbrace{\sum _{i=1}^N{m}_i{\dot{r}}_i{\dot{r}}_i{\vartheta }_i}}+\stackrel{\text{Extra Kinetic Curvature }\left(\mathrm{K}\mathrm{C}\right)\text{ Terms}}{\overbrace{\sum _{i=1}^N{m}_i{\dot{r}}_i{\dot{r}}_i\cdot \left[\frac{\mathrm{\partial }{\xi }_i^{+}}{\mathrm{\partial }{\mathbf{r}}_i}\mathrm{d}{S}_{zi}^{+}-\frac{\mathrm{\partial }{\xi }_i^{-}}{\mathrm{\partial }{\mathbf{r}}_i}\mathrm{d}{S}_{zi}^{-}\right]}}\\ & +\underset{\text{Extra Time Evolving }\left(\mathrm{T}\mathrm{E}\right)\text{ Terms}}{\underbrace{\sum _{i=1}^N{m}_i{\dot{r}}_i\left[\frac{\mathrm{\partial }{\xi }_i^{+}}{\mathrm{\partial }t}\mathrm{d}{S}_{zi}^{+}-\frac{\mathrm{\partial }{\xi }_i^{-}}{\mathrm{\partial }t}\mathrm{d}{S}_{zi}^{-}\right]}}\end{array}$(26) ) and (27(27) $\begin{array}{rl}\sum _{i,j}^N{\mathbf{f}}_{ij}{\mathbf{r}}_{ij}\cdot \underset{0}{\overset{1}{\int }}\frac{\mathrm{\partial }{\vartheta }_{\lambda }}{\mathrm{\partial }{\mathbf{r}}_{\lambda }}\mathrm{d}\lambda & =-\frac{\mathrm{\partial }}{\mathrm{\partial }\mathbf{r}}\cdot \stackrel{\stackrel{VA}{P^C}\mathrm{\Delta }V}{\overbrace{\sum _{i,j}^N{\mathbf{f}}_{ij}{\mathbf{r}}_{ij}\underset{0}{\overset{1}{\int }}{\vartheta }_{\lambda }\mathrm{d}\lambda }}\\ & +\underset{\text{Extra Configurational Curvature }\left(\mathrm{C}\mathrm{C}\right)\text{ Term}}{\underbrace{\sum _{i,j}^N{\mathbf{f}}_{ij}{\mathbf{r}}_{ij}\cdot \underset{0}{\overset{1}{\int }}\left[\frac{\mathrm{\partial }{\xi }_{\lambda }^{+}}{\mathrm{\partial }{\mathbf{r}}_{\lambda }}\mathrm{d}{S}_{z\lambda }^{+}-\frac{\mathrm{\partial }{\xi }_{\lambda }^{-}}{\mathrm{\partial }{\mathbf{r}}_{\lambda }}\mathrm{d}{S}_{z\lambda }^{-}\right]\mathrm{d}\lambda }},\end{array}$(27) ) including kinetic curvature $\mathrm{\partial }{\xi }_i/\mathrm{\partial }{\mathit{r}}_i$, KC, (

), time-evolving $\mathrm{\partial }{\xi }_i/\mathrm{\partial }t$, TE, (

) as well as the negative of the configurational curvature $\mathrm{\partial }{\xi }_{\lambda }/\mathrm{\partial }{\mathit{r}}_{\lambda }$, CC, (

) (shown as negative to allow comparison with the VA configurational pressure). In figure (b), the VA pressure is displayed with the correction terms added compared to the SF forms, which naturally include all of these terms. These include the kinetic VA with KC and TE added (

) compared to kinetic SF (

), configurational VA with CC added (

) against configurational SF (

) and total corrected VA (

) against total SF (

). The zero axis is shown by a horizontal black line (

).

Figure 9. (Colour online) Total tangential pressure divided by density for quantities in a reference frame moving with the interface, VA (

) and total SF (

). The kinetic, SF (

), VA (

), and configurational, SF (

), VA (

) components are shown faintly in the background for reference. A fitting is shown (

) with functional form and fitting coefficients annotated on the figure with the zero axis the horizontal black line (

).

We start with the case shown in the schematic of Figure 2(a). The wall normal ($P_N=P_{zz}$) 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 3(a) shows the area focused on in Figure 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 $\mathbf{\nabla }\cdot \mathit{P}=0$. The IK1 pressure has identical kinetic components 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 $\mathbf{\nabla }\cdot \mathit{P}\ne 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 give similar results, while the IK1 shows a difference in the normal configurational pressure in Figure 4 (a). The result again shows a flat line for the VA and SF forms, indicating the required $\mathbf{\nabla }\cdot \mathit{P}=0$ condition is satisfied. The normal component of SF pressure, shown by the solid lines, is equivalent to the Irving Kirkwood contour. The IK1 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 (Figure 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 Equation (12(12) $\stackrel{\mathrm{I}\mathrm{K}}{\mathrm{P}}\underset{T}{\overset{C}{}}=\frac{1}{2}\left[\stackrel{\mathrm{I}\mathrm{K}}{\mathrm{P}}\underset{xx}{\overset{C}{}}+\stackrel{\mathrm{I}\mathrm{K}}{\mathrm{P}}\underset{yy}{\overset{C}{}}\right]=\frac{1}{4}\sum _{i,j}^N\left[\frac{f_{xij}{x}_{ij}+f_{yij}{y}_{ij}}{|z_{ij}|}\right]H\left(\frac{z-z_i}{z_{ij}}\right)H\left(\frac{z_j-z}{z_{ij}}\right).$(12) ), but it is a combination of the VA form measured on an 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 $P_N$ contribution above and below the axis is 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 $P_N-P_T$ 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 $\mathbf{\nabla }\cdot \mathit{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 5(a), the IK1 and VA agree, showing peaks which mirror the density shown in grey, representing the location of the molecules relative to the surface. The SF term is shown and naturally includes the surface evolution term, $\mathrm{\partial }\xi /\mathrm{\partial }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 (Figure 5(a)), 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 (Figure 5(a)), so the sum of kinetic and configurational terms gives a flat pressure profile in Figure 6(a). We, therefore, see from Figure 6(a), only the SF form satisfies the equilibrium force-balance condition. The extra term in Equation (27(27) $\begin{array}{rl}\sum _{i,j}^N{\mathbf{f}}_{ij}{\mathbf{r}}_{ij}\cdot \underset{0}{\overset{1}{\int }}\frac{\mathrm{\partial }{\vartheta }_{\lambda }}{\mathrm{\partial }{\mathbf{r}}_{\lambda }}\mathrm{d}\lambda & =-\frac{\mathrm{\partial }}{\mathrm{\partial }\mathbf{r}}\cdot \stackrel{\stackrel{VA}{P^C}\mathrm{\Delta }V}{\overbrace{\sum _{i,j}^N{\mathbf{f}}_{ij}{\mathbf{r}}_{ij}\underset{0}{\overset{1}{\int }}{\vartheta }_{\lambda }\mathrm{d}\lambda }}\\ & +\underset{\text{Extra Configurational Curvature }\left(\mathrm{C}\mathrm{C}\right)\text{ Term}}{\underbrace{\sum _{i,j}^N{\mathbf{f}}_{ij}{\mathbf{r}}_{ij}\cdot \underset{0}{\overset{1}{\int }}\left[\frac{\mathrm{\partial }{\xi }_{\lambda }^{+}}{\mathrm{\partial }{\mathbf{r}}_{\lambda }}\mathrm{d}{S}_{z\lambda }^{+}-\frac{\mathrm{\partial }{\xi }_{\lambda }^{-}}{\mathrm{\partial }{\mathbf{r}}_{\lambda }}\mathrm{d}{S}_{z\lambda }^{-}\right]\mathrm{d}\lambda }},\end{array}$(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 Figure 6(a) 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 Figure 5(a) and, therefore, the total in Figure 6(a). 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 the 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 surface crossing measurements are offset either side by $\mathrm{\Delta }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 Figure 6(b), 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 [23]. 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 $\int (P_N-P_T)\mathrm{d}z$ , it is instructive to split them into these two additive contributions, as is done in Figure 7 with (a) the integral of the normal $\int P_N\mathrm{d}z$ and (b) the integral of the negative of the tangential $-\int P_T\mathrm{d}z$. 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 (Figure 7(a)) and minus tangential (Figure 7(b)), converges to roughly the same value of $\gamma =\int (P_N-P_T)dz\approx 0.6$ for all methods. Differences are attributed to statistical noise, which integrated quantities like surface tension are more susceptible to.

Figure 8(a) shows the VA pressure and the three additional terms, KC, 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 8(b), as required for mechanical equilibrium in a moving interface, as well as a good agreement with the SF pressure which naturally includes these terms. For the kinetic pressure, the normal VA contribution is shown in Figure 8(a) 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 in Figure 8(b) includes the surface movement correction in 8(a), and the kinetic curvature (both extra terms from Equation (26(26) $\begin{array}{rl}& \sum _{i=1}^N{m}_i{\dot{r}}_i\frac{\mathrm{d}{\vartheta }_i}{\mathrm{d}t}=\\ & \frac{\mathrm{\partial }}{\mathrm{\partial }\mathbf{r}}\cdot \stackrel{\stackrel{VA}{P}k\mathrm{\Delta }V}{\overbrace{\sum _{i=1}^N{m}_i{\dot{r}}_i{\dot{r}}_i{\vartheta }_i}}+\stackrel{\text{Extra Kinetic Curvature }\left(\mathrm{K}\mathrm{C}\right)\text{ Terms}}{\overbrace{\sum _{i=1}^N{m}_i{\dot{r}}_i{\dot{r}}_i\cdot \left[\frac{\mathrm{\partial }{\xi }_i^{+}}{\mathrm{\partial }{\mathbf{r}}_i}\mathrm{d}{S}_{zi}^{+}-\frac{\mathrm{\partial }{\xi }_i^{-}}{\mathrm{\partial }{\mathbf{r}}_i}\mathrm{d}{S}_{zi}^{-}\right]}}\\ & +\underset{\text{Extra Time Evolving }\left(\mathrm{T}\mathrm{E}\right)\text{ Terms}}{\underbrace{\sum _{i=1}^N{m}_i{\dot{r}}_i\left[\frac{\mathrm{\partial }{\xi }_i^{+}}{\mathrm{\partial }t}\mathrm{d}{S}_{zi}^{+}-\frac{\mathrm{\partial }{\xi }_i^{-}}{\mathrm{\partial }t}\mathrm{d}{S}_{zi}^{-}\right]}}\end{array}$(26) )). With these additional corrections, the kinetic part of the VA has an identical shape to the SF pressure which naturally includes these terms. For the configurational pressure on 8(a) shows the magnitude of the configurational curvature contributions, which are the extra terms missing from the VA measurements in Equation (27(27) $\begin{array}{rl}\sum _{i,j}^N{\mathbf{f}}_{ij}{\mathbf{r}}_{ij}\cdot \underset{0}{\overset{1}{\int }}\frac{\mathrm{\partial }{\vartheta }_{\lambda }}{\mathrm{\partial }{\mathbf{r}}_{\lambda }}\mathrm{d}\lambda & =-\frac{\mathrm{\partial }}{\mathrm{\partial }\mathbf{r}}\cdot \stackrel{\stackrel{VA}{P^C}\mathrm{\Delta }V}{\overbrace{\sum _{i,j}^N{\mathbf{f}}_{ij}{\mathbf{r}}_{ij}\underset{0}{\overset{1}{\int }}{\vartheta }_{\lambda }\mathrm{d}\lambda }}\\ & +\underset{\text{Extra Configurational Curvature }\left(\mathrm{C}\mathrm{C}\right)\text{ Term}}{\underbrace{\sum _{i,j}^N{\mathbf{f}}_{ij}{\mathbf{r}}_{ij}\cdot \underset{0}{\overset{1}{\int }}\left[\frac{\mathrm{\partial }{\xi }_{\lambda }^{+}}{\mathrm{\partial }{\mathbf{r}}_{\lambda }}\mathrm{d}{S}_{z\lambda }^{+}-\frac{\mathrm{\partial }{\xi }_{\lambda }^{-}}{\mathrm{\partial }{\mathbf{r}}_{\lambda }}\mathrm{d}{S}_{z\lambda }^{-}\right]\mathrm{d}\lambda }},\end{array}$(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 give 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 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 $P_T\left(z\right)$ and density $\rho \left(z\right)$ correspond exactly as a function of z, so the ratio is constantly 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 on either side of zero, as the intrinsic fitting process results in a large density peak at the interface and no molecules on either side. This means $P_T/\rho \to \mathrm{\infty }$ as $\rho \to 0$ but a finite value is observed at zero due to the molecules which sit on the intrinsic surface. The shape of the profile represents 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\pi z+d)$, where $a,b,c\mathrm{a}\mathrm{n}\mathrm{d}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 (34) $P_T\left(z\right)\approx \rho \left(z\right)e^{-z}\left[1-2.5\sin \left(2\pi z+0.5\right)\right]\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}-\mathrm{\infty }<z<-1.$(34) 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 coarse-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 Figures 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 Equation (16(16) $P_T^C=\frac{1}{4}\sum _{i,j}^N\left[\frac{f_{xij}{x}_{ij}+f_{yij}{y}_{ij}}{\mathrm{\Delta }z}\right]{\mathrm{\Lambda }}_{iz}.$(16) ) is, therefore, a similar quantity to the tangential IK1 components shown in Figures 5(b) and 6(b). This shows much smaller oscillations in Figure 6(b) 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 Figure 6(b) 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 (Equation (23(23) $\begin{array}{rl}& {\int }_{t_1}^{t_2}\rho \mathbf{u}{u}_z+\stackrel{SF}{P_z^{c+}}\mathrm{d}t=\frac{1}{\mathrm{\Delta }{S}_z}\sum _{i=1}^N{m}_i{\dot{r}}_i\frac{{\mathbf{r}}_{12}\cdot {\mathbf{n}}_z}{|{\mathbf{r}}_{12}\cdot {\mathbf{n}}_z|}\mathrm{d}{S}^{+}\\ & +\frac{1}{\mathrm{\Delta }{S}_z}\sum _{i=1}^N{m}_i{\dot{r}}_i{\vartheta }_t\stackrel{SF}{P_z^{c+}}=\frac{1}{2\mathrm{\Delta }{S}_z}\sum _{i,j}^N{\mathbf{f}}_{ij}\frac{{\mathbf{r}}_{ij}\cdot {\mathbf{n}}_z}{|{\mathbf{r}}_{ij}\cdot {\mathbf{n}}_z|}\mathrm{d}{S}^{+},\end{array}$(23) )) and then the IK1 for the tangential components using Equation (16(16) $P_T^C=\frac{1}{4}\sum _{i,j}^N\left[\frac{f_{xij}{x}_{ij}+f_{yij}{y}_{ij}}{\mathrm{\Delta }z}\right]{\mathrm{\Lambda }}_{iz}.$(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 the calculation of the SF interaction and surface normal is 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 would lose this property of exact control volume balance.

5. 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 $\xi =\xi (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 of 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 constant normal pressure over the interface. 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 a 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 equation linking SF pressure and instantaneous momentum change in the volume. As the SF pressures can be defined on all surfaces which 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.

Data Availability

The data that support the findings of this study are freely available from https://doi.org/10.24435/materialscloud:5b-0g saved in the open Python pickle format with all required Python scripts to load the data and reproduce the figures in this work (as well as allowing further analysis of the data). The input files to simulate this exact case are also included in the above link with a README.txt explaining the process of building and running Flowmol. The latest version of the Flowmol code is freely available from https://github.com/edwardsmith999/flowmol, however, to ensure reproducibility, the version of Flowmol used to generate the data in this work is available at the following persistent link https://doi.org/10.5281/zenodo.4639546.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

The author is grateful to the UK Materials and Molecular Modelling Hub for computational resources, which is partially funded by Engineering and Physical Sciences Research Council [EP/T022213/1].

Notes

1 The exact equality is lost because the VA form assumes pressure is constant in a volume, while the exact balance of the control volume requires fluxes averaged over the bounding surfaces of the volume.