Hydrodynamics Across a Fluctuating Interface
aa r X i v : . [ phy s i c s . c o m p - ph ] A ug Hydrodynamics Across a Fluctuating Interface
Edward R. Smith a) and Carlos Braga b)1) Department of Mechanical and Aerospace Engineering,Brunel University London School of Engineering and Materials Science, Queen Mary University of London
ABSTRACT
Understanding what happens inside the rippling and dancing surface of a liquidremains one of the great challenges of fluid dynamics. Using molecular dynamics(MD) we can pick apart the interface structure and understand surface tension. In thiswork we derive an exact mechanical formulation of hydrodynamics for a liquid-vapourinterface using a control volume which moves with the surface. This mathematicalframework provides the local definition of hydrodynamic fluxes at any point on thesurface. These are represented not only by the flux of molecules and intermolecularinteractions acting across the surface, but also as a result of the instantaneous localcurvature and movement of the surface itself. By explicitly including the surfacedynamics in the equations of motion, we demonstrate an exact balance between kineticand configurational pressure normal to the surface. The hydrodynamic analysis makesno assumptions regarding the probability distribution function, so is valid for anysystem arbitrarily far from thermodynamic equilibrium. The presented equationsprovide a theoretical basis for the study of time-evolving interface phenomena suchas bubble nucleation, droplet dynamics and liquid-vapour instabilities. a) [email protected] b) [email protected] . INTRODUCTION The liquid-vapour interface stands out as one of the great modelling challenges in engi-neering and the physical sciences. Central to this challenge is the observation that, whilebulk fluids and gases are well described by continuum models, this approach breaks downat the interface, where large changes in the physical properties, localised to a very smallregion, invalidate any continuum assumptions. Molecular dynamics (MD) is an ideal modelto study a problem of this type, modelling the liquid-vapour coexistence with no more thanthe solution of Newton’s law for a system of molecules.However, using a discrete MD model introduces a new challenge, translation of the in-formation at the molecular level to relevant continuum hydrodynamic properties. Thermalfluctuations, characteristic of the dynamics at the atomistic scale, are present at the interface,blurring any property of interest. Capillary wave theory (CWT) provides a framework todescribe these fluctuations, by representing the instantaneous shape of the fluid surface asan intrinsic surface in parametric form. Thermodynamic profiles across the interface rep-resent the convolution of an assumed sharp intrinsic profile with the Gaussian distributioncharacterising the height of the intrinsic surface.A mechanical route is also possible, based on instantaneous time evolving equations withno need for thermodynamics averaging, equivalent to the Newtonian mechanics which un-derpin molecular dynamics. Using the sharp intrinsic surface directly as a moving referenceframe, mechanical equations can be obtained from the molecular data with no thermody-namic blurring. The most important quantity describing the mechanical properties of theinterfacial region is the pressure tensor. For a homogeneous system of particles with periodicboundaries, a single virial pressure tensor can be defined for the whole system . In inho-mogeneous systems such as the liquid-vapour interface we require a local definition. Theseminal work of Irving and Kirkwood provides a localisation to give the pressure at a pointin space using the Dirac delta function. There are two considerations with the Irving andKirkwood stress, the practical and mathematical consequences of the Dirac delta functionin an MD simulation and the non-uniqueness in the definition of the pressure tensor .We consider the problematic Dirac delta functions first, a direct consequence of the math-2matical idealisation of a continuum. Several approaches have been applied to solve this, in-cluding mollifying the Delta function , integrating to get ”Volume Average” (VA) stress or reformulating in terms of the pressure over a surface , known as the method of planes(MOP) pressure. Of these approaches, the pressure over the surface has the advantage ofbeing conceptually the simplest form, namely the force acting over the surface divided bythe surface area Surf Π = F /A . It is also the only form to give exactly conservative equa-tions at every MD timestep, shown through the concept of a control volume . The controlvolume equations are important in fluid dynamics, providing the basis for the so-called con-servative finite volume approach in computational fluid dynamics (CFD) . By expressinga weakened statement of the equations of motion, conserved in an average sense over anarbitrary volume, they no longer demand a continuous field. This makes the control volumeapproach ideally suited to molecular dynamics, expressing everything as an average over avolume. This solves the problem of the Dirac delta functions by expressing the Irving andKirkwood equations in terms of integrated quantities inside a volume whose evolution isexactly equal to the sum of pressure, Surf Π , over the bounding surface.Next we consider the non-uniqueness, argued to make mechanical definitions unreliable .The ambiguity in pressure tensor consists of at least three possible components: i ) the in-clusion of kinetic terms , ii ) the intermolecular interaction path and iii ) the choice ofmeasuring reference frame . By providing an exact mechanical link between surface pressureand time evolution of momentum, the conservative nature of the control volume formulationhelps to sidestep some issues associated with the non-uniqueness of the pressure tensor. Theinclusion of kinetic terms i ) in the pressure has been shown to be essential so surface pressureis equal to the momentum change inside the control volume , a result extended in this workto include movement of the surface itself. The second source of non-uniqueness ii ), stemsfrom the infinite number of possible interaction paths between atoms. In this work, a linearpath is used, sometimes called the Irving Kirkwood contour, which is consistent with thedefinition of ”impressed force” used by Newton . Another common interaction path is theHarasima contour , known to give unphysical results for spherical coordinates and requireadjustment for cylindrical coordinates . The Harasima contour also combines path withreference frame, linking issues ii ) and iii ), as the tangential contour follows the interface.3he Harasima contour would follow the intrinsic interface in this work with varying normalsat each point, and may be expected to experience similar problems to the spherical andcylindrical cases. It would also be difficult and computationally expensive to implement asan integral following the intrinsic surface.The measuring reference frame, iii ), used in this work is a closed bounding surface ofa control volume. A direct consequence of using MOP style pressure defined over thisbounding surface, H Surf Π · d S = P SurfacesI F I /A I , is it simply adds forces between particlesover a surface, F I = P N Isurf i,j f ij , so we retain the exact balance of forces to momentumchange from Newtons laws. As noted by Schofield and Henderson , individual pressureis non-unique but the total pressure over the closed surfaces of a volume is invariant tochoice of contour. The result here extends this idea away from equilibrium, mathematicallylinking the invariant pressure over all surfaces to momentum change inside the volume, d/dt P N V ol i m i v i , and including deformation of the volume with movement of the intrinsicinterface. It is hoped that this exact link between the invariant closed surface pressure andthe resulting particle motion on an interface shown in this work could provide a way forwardin addressing the arbitrary nature of the pressure tensor and concerns that it renders resultsmeaningless . This treatment only considers two body interactions in this work, forthree body interaction or long-range contributions ensuring conservation would requiremore care.In addition to conservation, we take advantage of an additional benefit of the controlvolume approach, the ability to define them with arbitrary shape. We use the function formof the liquid-vapour interface, ξ , and define the face of our control volumes so the edge of ouraveraging grid follows the 2D surface. To define this surface, we apply the intrinsic interfaceapproach first proposed by Chac´on and Tarazona , fitting a set of Fourier components to theoutermost molecules in a liquid cluster. A range of other approaches to define an interfacehave been proposed , and the surface pressure presented in this work can, in principle,be applied to any other surface, not just the Chac´on and Tarazona functional form.An assumption of CWT is that the density and pressures profiles are uncorrelated withthe intrinsic surface from which the profiles are obtained. However, in order to balance thecontrol volume equations of motion, it will be shown that it is necessary to explicitly include4he pressure change due to the surface time evolution. The surface itself becomes part ofthe equations of motion, an addition which is demonstrated here to be essential in obtainingthe correct form of the equations of motion. Constant normal pressure over a liquid-vapourinterface is expected and has been shown in the literature with a fixed reference frame .Near solid-liquid interfaces a constant normal pressure must exist for momentum balance,which the IK1 form of pressure fails to show , an observation which motivates use of the VAand MOP formulations . Applying the same approach to the liquid-vapour interface, thenormal component of the intrinsic pressure does not show momentum balance even usingthe VA form of pressure .We show that the surface pressure equations, derived in this workgives, give the expected constant normal pressure. This flat profile requires the pressure toinclude both the instantaneous normal to the intrinsic surface at every time and a term for thesurface movement itself. The resulting equations instantaneously include all time evolvingterms balancing forces and fluxes at each timestep, meaning that the commonly assumedthermodynamic or average equilibrium condition ∇ · Π = 0 is neither required norvalid. As a result, the presented equations will allow an exploration of the instantaneoushydrodynamics of the interface, known to persist even with a very small number of moleculesmoving over very short timescales .The structure of this manuscript is as follows, in the theory section II the mathematicalform of the control volume is derived to give the expression for surface fluxes on a volumemoving with the intrinsic interface. The next section outlines the details of implementationIII, giving an overview of how to actually use these expressions in molecular simulationsand the details of the setup. The result and discussion are included in section V beforeconcluding remarks in section VI. II. THEORY
We start with a high level description of the process in this section. The control volumeformula for surface pressure are derived by starting from the definitions of Irving and Kirk-wood . These link the continuum expressions for density, momentum and energy to theirmolecular equivalents at any point in space. The process is then a formal integral of these5efinitions over a volume in space, followed by the evaluation of their time evolution to getthe mass and energy conservation as well as the momentum balance equations.The volume integral follows the interface, with functional form obtained from the in-trinsic surface method by fitting to the outer molecules at the interface. For the readeronly interested in applying these equations, the key formula to get pressure are given insection III B as Eq. (45) and Eq. (46). These have the simple interpretation of a gridwhich deforms and moves following the liquid-vapour interface. We obtain all intersectionsof molecular interactions or particle trajectories with surfaces of the cells in that grid. Theinteractions and forces over a surface are in the form of a force over area, F /A , the surfacepressure. When summed over all surfaces of any enclosed volume they exactly define thechange inside, a requirement of the validity of the conservation laws expressed in control vol-ume form. A piecewise bilinear approximation of the moving interface is used to get surfacecrossings more efficiently, and mapping applied to simplify calculations. The contributiondue to the movement of the interface itself is obtained by considering particles before andafter movement. A. The Irving Kirkwood Equations
The density at a given point in space can be obtained using Irving-Kirkwood’s procedure ,defined here without the ensemble average , to give the instantaneous quantity obtained atany time in a molecular dynamics simulation, ρ ( r , t ) = N X i =1 m i δ ( r − r i ) . (1)Here ρ is the continuum density at point r in three dimensional space and time t . The sumon the right adds the mass m i over all N molecules in the system, with the Dirac deltafunction δ only non-zero when r is equal to r i , that is molecule i is located at point r . Wecan define momentum, ρ ( r , t ) u ( r , t ) = N X i =1 m i ˙ r i δ ( r − r i ) , (2)6nd energy, ρ ( r , t ) E ( r , t ) = N X i =1 m i e i δ ( r − r i ) , (3)in an analogous manner, with ρ u and E continuum momentum and energy respectively. En-ergy here is e i = v i + 1 / P j = i φ ij based on half the energy of the intermolecular interaction. B. Control Volume integral
The continuum approximation represent reality as a continuous field, obtained by takingthe zero limit of an infinitesimal volume at each point. This same concept applied to adiscrete system results in, an infinitely thin and infinite large function at a point, the Diracdelta (formally a generalised function). The Dirac delta can be thought of as a usefulplaceholder representing the continuum assumption in a molecular system, but has limiteduse in practice, particularly in software implementation. A more tractable form of the Diracdelta function, for use in a discrete system, is obtained by the integration over an arbitraryvolume to get the control volume form . This has the advantage that a conservative setof equations can be obtained in a molecular system. These are directly relatable to theequivalent control volume expressions in the continuum. More importantly for this work,the control volume shape can be chosen based on the geometry of interest. In this work, thevolume is chosen to follow the intrinsic surface as it varies in time. The density at a point,Eq. (1), can be integrated over a volume as follows, Z V ρ ( r , t ) dV = N X i =1 m i Z V δ ( r − r i ) dV, (4)similar for momentum Eq. (2), Z V ρ ( r , t ) u ( r , t ) dV = N X i =1 m i ˙ r i Z V δ ( r − r i ) dV, (5)and energy Eq. (3), Z V ρ ( r , t ) E ( r , t ) dV = N X i =1 m i e i Z V δ ( r − r i ) dV. (6)7o evaluate the integral of Eqs. (1 - 3), only the Dirac delta function must be integrated. Thevolume described by the triple integral is between four flat surfaces and two faces followingthe shape of the intrinsic surface ξ , Z V δ ( r − r i ) dV = Z x + x − Z y + y − Z ξ + ξ − δ ( x − x i ) δ ( y − y i ) δ ( z − z i ) dzdydx, (7)here the arbitrary volume is cuboidal in x and y directions denoted by plus and minussuperscripts for top and bottom surfaces. The location of these surfaces are x + ≡ x + ∆ x/ x − ≡ x − ∆ x/ x the x width of the CV with centre at point x , witha similar definition using ∆ y as width in y . The surface in the z directions is described byposition z ± ≡ z ± ∆ z/ ξ ( x, y, t ), a continuous function of x , y and time t , so surfaceposition in z denoted by ξ ± ( x, y, t ) ≡ z ( t ) ± + ξ ( x, y, t ). In general, ξ can be any function andwe will make no assumption about its form until section IV A where a sum of trigonometricfunction will be used to describe the intrinsic surface between a liquid and vapour phase.As the ξ ± limits are a function of x and y , the z integral must be evaluated first, Z V δ ( r − r i ) dV = Z x + x − Z y + y − m i δ ( x − x i ) δ ( y − y i ) × (cid:2) H ( ξ + ( x, y, t ) − z i ) − H (cid:0) ξ − ( x, y, t ) − z i (cid:1)(cid:3) dydx. (8)The Dirac delta is the Heaviside function upon integration, with the finite limits ξ − and ξ + inserted. The next two integrals over x and y use the sifting property of the Dirac delta,namely R δ ( x − a ) f ( x ) = f ( a ) so the function ξ ( x, y, t ) which describes the surface roughnessbecomes expressed in terms of molecular position, ϑ i ≡ Z V δ ( r − r i ) dV = (cid:2) H (cid:0) x + − x i (cid:1) − H (cid:0) x − − x i (cid:1)(cid:3) × (cid:2) H (cid:0) y + − y i (cid:1) − H (cid:0) y − − y i (cid:1)(cid:3) × (cid:2) H ( ξ + ( x i , y i , t ) − z i ) − H (cid:0) ξ − ( x i , y i , t ) − z i (cid:1)(cid:3) = Λ xi Λ yi ˜Λ zi , (9)where ϑ i a function which selects molecules inside the Heavisides, called the control volumefunction . The Λ notation is also introduced, a box car, or Bracewell , function for eachdirection, so e.g. Λ xi = H ( x + − x i ) − H ( x − − x i ) and the tilde on Λ zi indicates the functionmoves with the intrinsic surface ˜Λ zi = H (cid:0) ξ i + − z i (cid:1) − H (cid:0) ξ i − − z i (cid:1) where the subscript i on ξ ξ i ± = ξ ± ( x i , y i , t ). As each Λ isone if the particle is between the limits in that coordinate direction, the product Λ xi Λ yi ˜Λ zi is therefore one if the particle is inside the volume; located between two intrinsic surfacefunctions in the z direction and bounded by two planes in the x and y directions. Anymolecules outside of this volume will return a value of zero. This control volume describedby ϑ i has a constant width and height of ∆ x , ∆ y respectively with the depth ∆ z alwaysconstant at any x, y location, as the same intrinsic interface function ξ is used for top andbottom surfaces, so ξ i + − ξ i − = ∆ z . C. Mass
The expression linking the control volume form of density in a continuum and molecularsystem of Eq. (4) can therefore be written concisely using Eq. (9) as, Z V ρ ( r , t ) dV = N X i =1 m i ϑ i . (10)That is, the mass of any molecule in the enclosed region between two intrinsic surfaces in z and four planes in x and y is contributing to the density in that control volume.We now use the control volume function to derive expressions for the fluxes over thesurface of a volume, the so-called flux forms of the equations of motion. In the continuum,the control volume analysis involves taking the time derivative of the mass in a volume toobtain the flux over the surfaces of that volume , ddt Z V ρ ( r , t ) dV = − I S ρ u · d S = Z S + ρ u · d S + − Z S − ρ u · d S − , (11)where the d S = n dS expresses an infinitesimal surface element dS with surface normal n and H indicates the integral is over all surfaces of the volume, which here represents six piecewisesurface integrals. The molecular control volume has been shown previously derived for auniform cuboid in space , a sphere , and is extended here for a general volume betweenintrinsic surfaces. Taking the time derivative of Eq. (10), ddt Z V ρdV = ddt N X i =1 m i ϑ i = N X i =1 m i dϑ i dt = d Λ xi dt Λ yi ˜Λ zi + Λ xi d Λ yi dt ˜Λ zi + Λ xi Λ yi d ˜Λ zi dt . (12)9he derivative of each of the three Λ functions will generate two terms for the top and bottomsurfaces corresponding to the six surfaces of the volume. For example in x the derivative d Λ xi /dt = dH ( x + − x i ) /dt − dH ( x − − x i ) /dt can be seen to give two terms for the top andbottom surface in the x direction. Consider just the top, or +, surface in x , ddt H (cid:0) x + − x i (cid:1) = − ˙ x i δ (cid:0) x + − x i (cid:1) , (13)can be seen to be the particle’s x velocity ˙ x i localised by a delta function to the x + surface,the mass flux of a particle over the surface. The same process can be applied for x − and y ± surfaces. We obtain the top surface in z from d ˜Λ zi /dt where both the molecular positionsand the intrinsic surface ξ i ± depend on time, ddt H (cid:0) ξ i + − z i (cid:1) = (cid:20) dξ i + dt − dz i dt (cid:21) δ (cid:0) ξ i + − z i (cid:1) = (cid:20) ˙ x i ∂ξ i + ∂x i + ˙ y i ∂ξ i + ∂y i + ∂ξ i + ∂t − ˙ z i (cid:21) δ (cid:0) ξ i + − z i (cid:1) . (14)The time derivative of ξ i + ( t ) = z + ( t ) + ξ ( x i ( t ) , y i ( t ) , t ) depends on particle positions x i and y i which are themselves a function of time, as well as the explicit surface time dependence.Each of the terms have a physical interpretation, the Dirac delta function is only non zerowhen molecules are crossing the surface, counting at the point of crossing 1) the z velocitycomponents of the molecule ˙ z i , 2) the surface curvature ˙ x i ∂ξ + /∂x i and ˙ y i ∂ζ + /∂y i times the x and y particle’s velocity at the location of a surface crossing and 3) the crossings due tosurface time evolution itself ∂ξ + /∂t .We evaluate Eq. (12) using the derivatives as shown in Eq. (13) and Eq. (14) to get thetime evolution of density in a molecular control volume between intrinsic surfaces, ddt N X i =1 m i ϑ i = − N X i =1 m i ˙ x i (cid:2) δ (cid:0) x + − x i (cid:1) − δ (cid:0) x − − x i (cid:1)(cid:3) Λ yi ˜Λ zi + ˙ y i (cid:2) δ (cid:0) y + − y i (cid:1) − δ (cid:0) y − − y i (cid:1)(cid:3) Λ xi ˜Λ zi + " (cid:18) ˙ z i − dξ i + dt (cid:19) δ (cid:0) ξ i + − z i (cid:1) − (cid:18) ˙ z i − dξ i + dt (cid:19) δ (cid:0) ξ i − − z i (cid:1) Λ xi ˜Λ zi ! . (15)10o simplify this expression, we introduce notation for the Dirac delta terms analogous to thecontinuum surface element dS used in the integral, dS ± xi ≡ δ (cid:0) x ± − x i (cid:1) S xi ; dS ± yi ≡ δ (cid:0) y ± − y i (cid:1) S yi ; dS ± zi ≡ δ (cid:0) ξ i ± − z i (cid:1) S zi . (16)with S αi is the product of boxcar functions in the other two directions, S xi ≡ Λ yi ˜Λ zi ; S yi ≡ Λ xi ˜Λ zi ; S zi ≡ Λ xi Λ yi ; (17)which can be seen to define an area in space between the surfaces, e.g. S xi is only one if aparticle is between y + and y − in y and ξ i + and ξ i + . Using the definitions of Eq. (16) in thetime evolution Eq. (15), ddt N X i =1 m i ϑ i = N X i =1 m i Surface Evolution z }| { ∂ξ i + ∂t dS + zi − ∂ξ i − ∂t dS − zi − X β ∈{ x,y,z } ˙ β i " dS + βi − dS − βi − ∂ξ i + ∂β i dS + zi + ∂ξ i − ∂β i dS − zi | {z } Curvature = − I S ρ u · d S , (18)where ∂ξ i ± /∂z i = 0 is used to write concisely and the last equality, linking continuum andmolecuale expressions, follows from Eq. (11). The sum over all surfaces in Eq. (18), allowsthe continuum and molecular expressions to be compared surface by surface, for β = x andtaking just the top + surface this is, Z S + x ρ u · d S + x = N X i =1 m i ˙ x i dS + xi , (19)and for β = z , against considering just the top surface, Z S + z ρ u · d S + z = N X i =1 m i h − ˙ x i ∂ξ i + ∂x i − ˙ y i ∂ξ i + ∂y i | {z } Curvature + ˙ z i − ∂ξ i + ∂t |{z} Surface Evolution i dS + zi . (20)So, the z + surface flux of mass is made up of direct fluxes, surface curvature and surfaceevolution components which must all be evaluated. The shorthand notation for flux overeach surface is introduced, ddt N X i =1 m i ϑ i = − N X i =1 m i h ˙ β + i · d S + i − ˙ β − i · d S − i i = − N X i =1 m i ˙ β i · d S i , (21)11here ˙ β ± i = ˙ x i ˙ y i ˙ z i − ∂ξ i ± ∂x i − ∂ξ i ± ∂y i − ∂ξ i ± ∂t and d S ± i = dS ± xi dS ± yi dS ± zi (22)and the ± superscript is omitted in Eq. (21) when expressing the difference d S i = d S + i − d S − i or more generally could be a shorthand to denote flux over an arbitary surface. D. Momentum
The same integration process used for density can be applied to get the momentum in acontrol volume, written in integrated form as, Z V ρ ( r , t ) u ( r , t ) dV = N X i =1 m i ˙ r i ϑ i . (23)The equivalent time evolution of momentum for a continuum control volume gives the surfaceflux of momentum (convection) ρ uu and pressure tensor denoted by Π , ddt Z V ρ u dV = − I s [ ρ uu + Π ] · d S = − I s (cid:2) ρ uu + Π k + Π c (cid:3) · d S , (24)where the total pressure tensor can be split into kinetic pressure, Π k , and configurationalpressure, Π c , contributions. It is perhaps more natural to talk about the compressed stateof molecules in the configurational contributions as stress, but as this is simply the negativeof pressure, the term pressure is used here for both kinetic and configurational contributions.We could also explicitly include the term for surface tension H γd ℓ in these continuum equa-tions, but it is more convenient to include all terms in the pressure tensor and calculate thesurface tension from the stress tensors using the Kirkwood and Buff approach.As with the mass equation, we evaluate the time evolution of Eq. (23) to obtain expres-sions for the molecular surface flux and stresses. obtained from manipulation of the Diracdelta functions, ddt Z V ρ u dV = ddt N X i =1 m i ˙ r i ϑ i = N X i =1 m i ˙ r i dϑ i dt | {z } Kinetic + N X i =1 m i ¨ r i ϑ i | {z } Configurational . (25)12he Kinetic term proceeds along the same lines as density, above, which gives, N X i =1 m i ˙ r i dϑ i dt = − N X i =1 m i ˙ r i ˙ β i · d S i = − I S (cid:2) ρ uu + Π k (cid:3) · d S . (26)We consider the Configurational term next, N X i =1 m i ¨ r i ϑ i = N X i =1 F i ϑ i = 12 N X i,j f ij [ ϑ i − ϑ j ] , (27)where the i, j notation is shorthand for a double sum over all i and all j indices. Equa-tion (27) is the integral of the differences of two Dirac delta functionals, the infamous IKoperator . Using the fundamental theorem of the calculus, this can be expressed in a muchmore convenient form, ϑ i − ϑ j = Z V [ δ ( r − r i ) − δ ( r − r j )] dV = Z V Z ∂∂λ δ ( r − r λ ) dλdV = Z ∂∂λ Z V δ ( r − r λ ) dV dλ = Z ∂ϑ λ ∂λ dλ = Z ∂ r λ ∂λ · ∂ϑ λ ∂ r λ dλ = Z r ij · ∂ϑ λ ∂ r λ dλ, (28)where r λ = r i + λ r ij , representing an integration along the line of interaction between r i and r j . The intrinsic surface ξ ( x, y, t ) is independent of the dummy integral along the pathbetween molecules, λ , so we can change the order of integration allowing a volume integralof the Delta functions which follows the same process as Eq. (7) to Eq. (9) with molecularposition replaced by point on line of interaction r λ → r i , giving, ϑ λ ≡ (cid:2) H (cid:0) x + − x λ (cid:1) − H (cid:0) x − − x λ (cid:1)(cid:3) (29) × (cid:2) H (cid:0) y + − y λ (cid:1) − H (cid:0) y − − y λ (cid:1)(cid:3) × (cid:2) H ( ξ + ( x λ , y λ , t ) − z λ ) − H (cid:0) ξ − ( x λ , y λ , t ) − z λ (cid:1)(cid:3) = Λ x λ Λ y λ ˜Λ z λ . Here the control volume function selects molecular interactions through Heavisides to obtainthe length of line in a control volume, by integrating along that line where any point is onewhen r λ is inside and zero otherwise. We denote the intrinsic interface function in terms of r λ as ξ ± λ ≡ z ± + ξ ( x λ , y λ , t ). The functions Λ x λ , Λ y λ and ˜Λ z λ are analogous to the previousdefinitions with r i replaced by r λ . 13he expression r ij · ∂ϑ λ /∂ r λ in Eq. (28) is a sum over all three directions, consideringthe z component first, z ij ∂ϑ λ ∂z λ = − z ij (cid:2) dS + zλ − dS − zλ (cid:3) , (30)where the dS ± αλ term is defined following the convention used for Eq. (16) to be, dS ± xλ ≡ δ (cid:0) x ± − x λ (cid:1) S xλ ; dS ± yλ ≡ δ (cid:0) y ± − y λ (cid:1) S yλ ; dS ± zλ ≡ δ (cid:0) ξ ± λ − z λ (cid:1) S zλ , (31)with S xλ = Λ y λ ˜Λ z λ , S yλ = Λ x λ ˜Λ z λ and S zλ = Λ x λ Λ y λ . The derivatives of the x and y components are slightly more complicated due to the x λ and y λ dependency in ζ , so for x we have, x ij ∂ϑ λ ∂x λ = − x ij (cid:20) dS + xλ − dS − xλ − ∂ξ + λ ∂x λ dS + zλ + ∂ξ - λ ∂x λ dS − zλ (cid:21) . (32)A similar result can also be obtained for y . The derivative in x gives the same surface cur-vature terms ∂ξ/∂x seen previously for molecular flux but in this case due to intermolecularinteractions. Combining the six surfaces of the control volume, the configurational term ofEq. (27) can be written as,12 N X i,j f ij [ ϑ i − ϑ j ] = 12 N X i,j f ij Z β ij · d S λ dλ = I S Π c · d S , (33)which is the total configurational pressure over all surfaces. Defining an analogous β ij and d S λ to the kinetic flux case, β ± ij = x ij y ij z ij − ∂ξ ± λ ∂x λ − ∂ξ ± λ ∂y λ and d S ± i = dS ± xλ dS ± yλ dS ± zλ (34)We can therefore take any of the surfaces, for example the configuration pressure on the x surface is, Z S + x Π c · dS + x = 12 N X i,j f ij Z x ij dS + xλ dλ, (35)14nd the z surface is, Z S + z Π c · dS + z = 12 N X i,j f ij Z h − x ij ∂ξ + λ ∂x λ − y ij ∂ξ + λ ∂y λ | {z } Curvature + z ij i dS + zλ dλ. (36)We are now in a position to write Eq. (25), the time derivative of the control volume interms of the kinetic Eq. (26) and configurational Eq. (33) parts, ddt N X i =1 m i ˙ r i ϑ i = − N X i =1 m i ˙ r i ˙ β i · d S i − N X i,j f ij Z β ij · d S λ dλ = − I S [ ρ uu + Π ] · d S , (37)which is the total stress over every surface of the control volume. Each of the faces of thecontrol volume, top + or bottom − can be seen to define three of the components of stresstensor. Considering the top x surface, with dS + xi and S + xλ and assuming an average pressureon the surface area R S + x [ ρ uu + Π ] · d S + x ≈ ∆ S x [ ρ u u x + Π x ], with Π x = [Π xx , Π xy , Π xz ] T sothe pressure on the x surface can be written as, ρ u u x + Π x = 1∆ S x N X i =1 m i ˙ r i ˙ x i dS + xi + 12∆ S x N X i,j f ij Z x ij dS + xλ dλ, (38)which is identical to the flat surface obtained in previous work and is consistent with alocalised method of planes stress (see section III). The pressure on the z surface is, ρ u u z + Π z = 1∆ S z N X i =1 m i ˙ r i h Kinetic Curvature z }| { ˙ x i ∂ξ i + ∂x i + ˙ y i ∂ξ i + ∂y i + ˙ z i + Surface Evolution z}|{ ∂ξ i + ∂t i dS + zi + 12∆ S z N X i,j f ij Z h x ij ∂ξ + λ ∂x λ + y ij ∂ξ + λ ∂y λ | {z } Configurational Curvature + z ij i dS + zλ dλ, (39)where the equation describes the components of pressure on a control volume surface whichis following the intrinsic interface, including terms due to curvature and surface evolution.Three connected faces form a tetrahedron which is consistent with Cauchy’s original defini-15ion of the stress tensor , giving,[ Π x , Π y , Π z ] = Π xx Π yx Π zx Π xy Π yy Π zy Π xz Π yz Π zz (40) |{z} x surf |{z} y surf |{z} z surf Section III will discuss how to implement this equation.
E. Energy
For completeness, the expressions for energy are stated here as they require no additionalmathematics. The time derivative of Eq. (6) for energy results in a similar kinetic andconfigurational term to the momentum equation, ddt N X i =1 e i ϑ i = N X i =1 e i dϑ i dt + N X i =1 de i dt ϑ i = N X i =1 e i ˙ r i · dϑ i d r + 12 N X i,j ˙ r i · f ij [ ϑ i − ϑ j ]= N X i =1 e i ˙ β i · d S i + 12 N X i,j ˙ r i · f ij Z β ij · d S λ dλ = I S [ ρ u E + Π · u + q ] · d S . (41)Where any surface of the control volume gives the energy equation in that coordinate di-rection. Isolating and obtaining contributions for heat flux q could proceed as outlined inprevious work , although careful consideration of the evolving surface and its contributionto heat flux may be required. This is left for future work. III. IMPLEMENTATION
In this section, the mathematical equations are manipulated to obtain expressions whichcan be coded in an MD simulation. The similarity between the kinetic and configurationalterms will be highlighted, as well as the similarity in operation required for both intrinsicand flat surfaces of the volume. It will be shown that the problem of obtaining the surfacestress is reduced to obtaining the intersection of a line and a surface, a common problem16n computer graphics and ray tracing. From equation (39), it is apparent the form of both
Kinetic Curvature and
Conf igurational Curvature terms are similar, with the differencebetween them the surface time evolution ∂ξ + /∂t , which we consider first. A. Time Evolving Interface
This term, ∂ξ + /∂t , describes the change in mass, momentum or energy in a controlvolume as new molecules are absorbed or left behind when the intrinsic surface moves. Thiscan be thought of as a 2D function sweeping through space and crossing the position ofthe particles. To evaluate this term, we first consider the time integration process appliedin an MD simulation, i) particle positions at time t are used for the force calculation, ii)the surface ξ ( t ) ≡ ξ ( x i ( t ) , y i ( t ) , t ) is fixed at time t and the evolution of the particles from r i ( t ) to r i ( t ) is used to get surface flux and iii) the particles are fixed at t while the surfaceis evolved from ξ ( t ) to ξ ( t ). Considering the top surface in z , this proceeds as follows, Z t t ∂ξ i + ∂t dS + zi dτ = Z t t ∂ξ i + ∂t δ (cid:0) ξ i + ( t ) − z i (cid:1) S zi dτ = Z t t N roots X k =1 ∂ξ i + ∂t δ ( t − t k ) | ∂ξ i + ( t k ) /∂t | S zi dτ = N roots X k =1 sgn (cid:18) ∂ξ i + ( t k ) ∂t (cid:19) [ H ( t − t k ) − H ( t − t k )] S zi , (42)using the roots of the Dirac delta function given in the Appendix Eq. (A7). The derivativeof the surface can be expressed, using the definition of the partial derivative in time with∆ t = t − t , ∂ξ i + ∂t = lim ∆ t → ξ i + ( x i ( t ) , y i ( t ) , t + ∆ t ) − ξ i + ( x i ( t ) , y i ( t ) , t )∆ t ≈ ξ i + ( t ) − ξ i + ( t ) t − t , so, the expression sgn ( ∂ξ i + /∂t ) = sgn (cid:0) ξ i + ( t ) − ξ i + ( t ) (cid:1) as t > t and this term can beseen to simply determine the crossing direction. Obtaining the multiple potential roots, t k ,of a time evolving polynomial crossing position in space is a non-trivial exercise, but theexpression sgn (cid:0) ξ i + ( t ) − ξ i + ( t ) (cid:1) [ H ( t − t k ) − H ( t − t k )] can be seen to be equivalent to asimple check if points z i ( t ) is crossed as the surface moves from ξ i + ( t ) to ξ i + ( t ), which can17e achieved by the following functional, sgn (cid:0) ξ i + ( t ) − ξ i + ( t ) (cid:1) [ H ( t − t k ) − H ( t − t k )] S zi = (cid:2) H (cid:0) ξ i + ( t ) − z i ( t ) (cid:1) − H (cid:0) ξ i + ( t ) − z i ( t ) (cid:1)(cid:3) S zi ≡ ϑ t , (43)Note that unlike the previous control volume functionals ϑ i and ϑ λ , it is possible for ϑ t tobe negative. B. Summary of Equations for Pressure
In this section, the new surface pressure equations derived in the previous section arepresented in a form that can be implemented in an MD simulation. These equations will becompared to the volume average (VA) form of pressure derived in previous work , obtainedby integrating the Irving and Kirkwood expressions over a volume following an intrinsicsurface, to give, Z V (cid:2) ρ uu + VA Π (cid:3) dV = N X i =1 m i ˙ r i ˙ r i ϑ i + 12 N X i,j f ij r ij Z ϑ λ dλ. (44)A more general derivation of these volume average expression follows from the time evolutionof Eq. (25), as shown in appendix C.The equations for pressure on a flat surface Eq. (38) and an intrinsic surface Eq. (39)must be integrated so they can be used in a molecular simulation. The process of takingthis integral is given in appendix A, with just the final forms stated here. The equationfor pressure on the flat control volume faces, here the y + surface is chosen, is shown in theappendix A to be, Z t t (cid:20) ρ u u y + Surf Π ky (cid:21) dt = 1∆ S y N X i =1 m i ˙ r i r i · n y | r i · n y | (cid:20) H (cid:18) y + − y i y i (cid:19) − H (cid:18) y i − y + y i (cid:19)(cid:21) Λ x ( t k ) ˜Λ z ( t k ) Surf Π cy = 12 1∆ S y N X i,j f ij r ij · n y | r ij · n y | (cid:20) H (cid:18) y + − y j y ij (cid:19) − H (cid:18) y i − y + y ij (cid:19)(cid:21) Λ x ( λ k ) ˜Λ z ( λ k ) . (45)Equation (45) is written in this form to emphasise the expression is the molecular pressuretensor m i ˙ r i r i and f ij r ij dotted with the surface normal n y . The r i = r i − r i term is18he vector from the position of molecule i at time t to its position at t . This Heavisidefunctions H in the square brackets check if the position y i before moving and y i after areon opposite sides of the surface at y + . These are equivalent to the method of planes (MOP)form of stress, as demonstrated in the Appendix A, but localised to a control volume surfaceby Λ x ˜Λ z . Recall that the Lambda functions are boxcar or Bracewell functions checking ifa point is between two points using Heaviside functions Λ x ( a ) = H ( x + − a ) − H ( x − − a )and ˜Λ z ( a ) = H ( ξ i + − a ) − H ( ξ i − − a ). This control volume surfaces and associated normalsare shown in Figure 1 a ) for an intermolecular interaction crossing the y + surface by greencrosses. An example of the actual interactions from an MD simulation at a single timestepare also shown in Figure 1 b ) with blue and green crosses denoting the x and y surfaces,respectively. The expressions for the roots t k and λ k can be obtained analytically for thisflat surface case, with expression given in the Appendix A. The expression for stress on theintrinsic surface z + is, Z t t ρ u u z + Surf Π kz dt = 1∆ S z N X i =1 m i ˙ r i r · ˜ n z | r · ˜ n z | N roots X k =1 [ H (1 − t k ) − H ( − t k )] Λ x ( t k )Λ y ( t k )+ 1∆ S z N X i =1 m i ˙ r i ϑ t Surf Π cz = 12∆ S z N X i,j f ij r ij · ˜ n z | r ij · ˜ n z | N roots X k =1 [ H (1 − λ k ) − H ( − λ k )] Λ x ( λ k )Λ y ( λ k ) , (46)where, as in the flat surface, the expression for m i ˙ r i r i and f ij r ij are dotted with the surfacenormal. In this case, the normal to the intrinsic interface ˜ n z ( x, y ) = ∇ α ( ξ − z α ) || ∇ α ( ξ − z α ) || . The normal˜ n z is shown on the schematic of Figure 1 a ) and the MD simulation of Figure 1 b ) with themany surface intersections shown as red crosses. The intersection of the interaction line andthe intrinsic surface is required to evaluate the expression of Eq. (46), an identical roots-finding process for crossing time t k or point on intermolecular interaction λ k . Given themultiple possible crossings, a closed form expression is not possible. An approach based onapproximating the surface as a set of bilinear patches is used here which allows applicationof an efficient root calculation from the ray-tracing literature . The computations are thenfast enough to be part of a molecular simulation. Once the roots t k and λ k have been found,19hey can be put into the Λ functions to check if they are within the limits of a given controlvolume surface. In practice, integer division can be used to speed this process up, assigningthe crossing to a cell, as discussed in the next section. Finally, the surface evolution termis obtain by checking if a particles positions r i ( t ) is in a given volume before and after thesurface has evolved.The result is a concise expression for the momentum flux over all surfaces, written as, Accumulation z }| { N X i =1 [ m i ˙ r i ( t ) ϑ i ( t ) − ˙ r i ( t ) ϑ i ( t )] = − X α = { x,y,z } S α " Advection z }| { N X i =1 m i ˙ r i (cid:0) r α (cid:2) dS + αt k − dS − αt k (cid:3) + ϑ t (cid:1) + ∆ t N X i,j f ij r αij (cid:2) dS + αλ k − dS − αλ k (cid:3)| {z } Forcing , (47)where both sides are integrated over time with the force term shown integrated with themidpoint rule for simplicity R dt ≈ ∆ t , consistent with the Verlet integration scheme usedto propagate molecular positions and velocities. The function to get crossings on a surfaceis denoted for time s = t and space s = λ with, dS ± αs k = ˜ n α | r · ˜ n α | N roots X k =1 [ H (1 − s k ) − H ( − s k )] Λ β ( s k )Λ γ ( s k ) , (48)where the β and γ are the orthogonal directions to α , so if α = x , β = y and γ = z , andthe normal ˜ n α can be for a flat or intrinsic surface with roots s k obtained for that surface.This same expression is valid for both kinetic and configurational terms on both curved andflat surfaces, with the simple interpretation of checking if a crossing is on the surface of acontrol volume.The six surface pressures summed for any arbitrary control volume consist of terms dueto intermolecular forces, labelled F orcing in Eq. (47) as well as kinetic molecular crossingsdue to both molecular motion and surface movement labelled
Advection . These are exactlyequal, to machine precision, to the change in momentum in the volume, called
Accumulation ,demonstrated in section V. 20
IG. 1. a ) A schematic showing surface normals n y and ˜ n z on an intermolecular interaction line r ij contributing to pressure on both the flat surfaces via Eq. (45) and the intrinsic surface withEq. (46) and b ) a random snapshot of an arbitrary control volume inside the liquid phase of anMD simulation, showing all the contributions used in the calculation of configurational pressure;with blue crosses on the x surface, green crosses on the y and red crosses for a bilinear patch of theintrinsic surface with the top ξ + surface coloured in light blue and the bottom ξ − coloured lightred. We now discuss the process of how to obtain these surface crossing terms efficiently aspart of an MD simulation.
IV. METHOD
We use molecular dynamics (MD) simulation to model the liquid vapour interface, witha shifted Lennard Jones potential, φ ( r ij ) = ǫ "(cid:18) σr ij (cid:19) − (cid:18) σr ij (cid:19) − φ ( r c ); r ij < r c , (49)with cutoff r c = 2 .
5, which is shorter than required to give good agreement with the experi-mental measurements of the surface tension in argon but chosen to allow more efficientsimulations. The force on particle i is obtained from the sum of the gradient of this potentialdue to interaction with all N other particles F i = P Nj = i f ij = − P Nj = i ∇ φ ( r ij ) The simula-tions are run using the Flowmol MD code which has been validated extensively in previous21ork The time integration is achieved by velocity Verlet with a timestep ∆ t = 0 . r i ( t + ∆ t ) = r i ( t ) + ∆ t v i ( t + ∆ t/ v i ( t + ∆ t/
2) = v i ( t − ∆ t/
2) + ∆ t F i ( t ) (50)The system is initialised as a liquid-vapour coexistence by creating an FCC lattice andremoving molecules until the desired density is obtained. The system is periodic in alldirections, with Lx = Ly = 12 . Lz = 47 .
62. The middle40% of the domain is designated to be initialised as liquid with a density of ρ l = 0 .
79 andthe remaining domain is set to a gas density of ρ g = 0 . N = 2635 molecules. The system is then run in the NVT ensemble, controlled bya Nos`e Hoover thermostat at a temperature setpoint of T s = 0 . a ). In this simulation,the liquid region tends to be located in the centre of the domain and the cluster analysisidentifies all connected molecules which are within the Stillinger cutoff length r d = 1 . N ℓ liquid particles. The cluster isthen used to fit the intrinsic surface, as detailed in section IV A, with the fitting performedon the surface on the right-hand side of the cluster. A. Intrinsic Surface
Up until this section, no assumption has been made about the functional form of surface ξ ( x, y, t ). To give a general form, the intrinsic surface method (ISM) approximates theliquid-vapour interface using a Fourier series representation, ξ ( x, y, t ) = X k The control volume functions ϑ i and ϑ λ mathematically determine if a given point r α islocated in any given cell as it moves with the intrinsic interface. In practice, a contiguousgrid of cells is used, so instead of testing every point with every cell, the point locations areassigned to an appropriate cell using integer division. In pseudocode for this is just, function g e t _ c e l l ( r )cell (1) = c e i l i n g ( r (1)/ c e l l s i d e (1)) Λ x cell (2) = c e i l i n g ( r (2)/ c e l l s i d e (2)) Λ y cell (3) = c e i l i n g ( r (3)/ c e l l s i d e (3)) Λ z r e t u r n cellend function g e t _ c e l l where r ( α ) is r α , the point in three dimensions and rounded division by cellside( α ) returnsan integer which is the index of the cell where the point is located. Quantities can then beadded to a cell using this index, including mass, velocity, energy when the point is a molecular24 IG. 2. A schematic outlining the process of defining the surface fluxes starting by a ) getting theliquid cluster, b ) fitting the initial surface to a number of outermost particles, with 3 show hereof the 3 × c ) refining by including all closest molecules until surface density is a targetvalue ρ s = n /A , d ) defining a set of layers with uniform offset ∆ z from the surface e ) Defininga grid by dividing the domain into a uniform grid of ∆ x and ∆ y cells and f ) approximating theintrinsic interface as bilinear patches for each cell, with the calculation of crossings shown for flatsurface as green crosses, the shaded region showing the cells between these flat surface crossingsthat are checked for intrinsic surface crossings with the two possible crossings shown by red crosses.A snapshot of the actual intrinsic surface ontained in the molecular dynamics simulation is shownin g ) with red pivot molecules and the intrinsic surface shown as a black grid split into 32 × x , Λ y and Λ z defined in Eq. (9), with no tilde in the z direction as this is for a flat surface. The intrinsicsurface form ˜Λ z can then be implemented by a simple mapping of the point followed byrounding to obtain the cell index. function g e t _ i n t r i n s i c _ c e l l ( r )cell (1) = c e i l i n g ( r (1)/ c e l l s i d e (1)) Λ x cell (2) = c e i l i n g ( r (2)/ c e l l s i d e (2)) Λ y ˜Λ z cell (3) = c e i l i n g (( r (3) - xi ( r (1) , r ( 2 ) ) ) / c e l l s i d e (3))r e t u r n cellend function g e t _ i n t r i n s i c _ c e l l where xi is ξ at each timestep, a function which returns the z position of the intrinsic surfacegiven inputs of position x and y . These can be x i , y i , x λ , y λ , or x k , y k shown in Figure 3where this mapping is applied to molecular positions, intermolecular interaction modelled ina piecewise manner or surface crossings, respectively. This has the convenient property thatthe intrinsic surfaces mapping can be switched on in a code by simply swapping a functionpointer from the get cell function to the get intrinsic cell function. For large domains in thesurface tangential directions, the time taken by function xi to evaluate surface position canbe prohibitive as it requires a sum over the two-dimensional Fourier surface, a calculationof order O (4 k u ) with k u scaling as L (assuming L x = L y = L ) so each surface evaluationis O ( L ). The position of the surface at the location of every molecule is pre-calculatedfor efficiency to get density, momentum, temperature, potential energy and the kinetic partof the pressure tensor. Interactions between molecules must be obtained on the fly bydiscretising the line of interaction between two molecules and using the surface value at eachsegment of the line, as shown in Fig 3 by the dotted interaction lines. Some optimisationscan be introduced, for example a quick estimate of the number of bins crossed betweentwo points is used to minimise the number of segments in the line and the intermolecularstress and work calculation are performed at the same time. However, this still requires a26arger number of evaluations of the intrinsic surface as each interaction crosses O (20) bins,for an average of O (100) interactions per particle (with r c = 2 . 5) required for each of the N particles. This is even worse for the control volume flux calculation, which requires theintersection of the intrinsic interface and the line of interaction, computational prohibitive asthe intersection point is required to machine precision if we want to verify exact conservationof momentum. Even an efficient root finding algorithm with good initial estimate requires thesurface to be evaluated many times to converge for every one of the O (20) crossed surface byevery one of the O (100) interactions. These rough numbers clearly depend on bin resolution,cut off length as well as the density of liquid and vapour states, but even in low resolutioncases appear to represent a limiting step computationally. A fast evaluation is therefore ofparamount importance for intrinsic interface tracking calculation to be used as part of anMD simulation, achieved here using a bilinear approximation. C. Bilinear Approximation for Surface Crossings The surface is replaced by a piecewise continuous bilinear surface, sampling the intrinsicsurface at an arbitrary accuracy. This uses a bilinear patch of the form, ξ BL ( x, y ) = a + a x x + a y y + a xy xy (54)where the coefficients are obtained by solving, a a x a y a xy = x y x y x y x y x y x y x y x y − ξ ( x , y ) ξ ( x , y ) ξ ( x , y ) ξ ( x , y ) (55)for each patch of the intrinsic interface. Once we have the expression for the bilinear surface,this can be used in the mapping for both particle position and lines between molecules. Inorder to get surface crossings, we use the fast bilinear-patch line intersection algorithm which equates the solution of a line of interaction r s = r + sr and the bilinear patch ξ BL to obtain a quadratic equation. Solving this provides up to two crossings for a given surface,27ith some care required for zero or complex roots. As the surface is split into piecewisebilinear patches, this also has the advantage that each control volume has its own bilinearsurface and the conservation check of Eq. (47) can be performed on each cell. The derivativescan also be obtained analytically ∂ξ BL /∂x = a x + a xy y and ∂ξ BL /∂y = a y + a xy x , althoughthese are discontinuous at the interface of contiguous bilinear patches. In practice this doesnot represent a problem for the surface pressures as the Λ α functions ensure all quantitiesare per cell. The advantage of supplementing the full intrinsic surface with the bilinearapproximation, instead of directly fitting a bilinear surface, is the calculation of quantitiesfrom the full intrinsic surface ξ are available analytically, if needed, with the bilinear solutionused as a starting point. Bicubic or higher order surface approximations could also be usedbut would not be able to take advantage of the fast interaction calculation for line andbilinear patch.Determining the number of crossings of an intermolecular interaction and the 3D gridfollowing the intrinsic surfaces (see Fig 2 f )) turns out to be a non-trivial exercise. Thisis because it is not possible to predict a priori how many times the intermolecular lineswill cross an arbitrary intrinsic surface split into interconnected ξ BL patches from just theknowledge of the starting and ending points. To solve this, for each interaction we first obtainall the crossings locations on all flat surfaces in x and y , made possible as we have the closedform expression from Eq. (45). By ordering all intermediate flat surface ( x and y ) crossingsalong the line of interaction, including the two particle positions at the end, we can takesuccessive pairs and are guaranteed any intrinsic surface crossings between a pair must be onthe same ξ BL patch. This is shown schematically in Fig 2 f ) with green crosses denoting apair of flat plane crossings, the shaded region shows the row of cells in between these two flatcrossings to be checked for intrinsic surface intersections and the red crosses denoting twofound intrinsic crossings. As each grid volume must have at least two crossings or contain anend point (molecule position), by stepping in pairs we ensure all surface crossings have beenobtained by this process. This has the added advantage that the coefficients of the bilinearsurface, ξ BL , can be loaded once for each row, the shaded region in Fig 2 f ), and used tocheck all line-plane crossings with the efficient Ramsey et al. algorithm. The straight lineto check is then between the current pair of flat surface crossings and patch ξ BL is shifted28 IG. 3. A schematic view of the three operations required to obtain density, momentum andpressures including i) binning molecular position z α → z i , ii) binning lines between moleculesapproximated as a series of points z α → z λ and iii) obtaining the location of surface crossings z α → z k . The two Heaviside functions are essentially a binning operation, checking if a point isbetween ξ − and ξ + with ξ ± = z ± + ξ , shown in red, which can be mapped to a simpler function onthe right by subtracting the intrinsic surface value at the particle locations ξ ( x i , y i , t ) or point online ξ ( x λ , y λ , t ). This mapping changes the calculation to a simple check if molecules are betweenflat surfaces z − and z + by integer division. Intersections on flat surfaces, denoted by a green cross,can be mapped in the same way to find the surface of the CV they apply, obtained by Eq. (45).Mapping is not possible for crossings of the intrinsic surface, denoted by a red cross before mapping,but a question mark after to emphasise a root finding must be employed. in multiples of ∆ z to get the the successive bin surfaces between the z coordinates of thesurface crossings . The pseudocode for this process is as follows, function s u r f a c e _ f l u x e s ( r1 , r2 , quantity , f l u x e s ) C h e c k if same cell in both x and y f ) if ( c e l l A (1) == c e l l B (1) & c e l l A (2) == c e l l B (2) V. RESULTS AND DISCUSSION In this section, the results for pressure on the moving surface are presented using the newsurface flux equations derived in Sections II, summarised for implementation in III with setupand methodology in IV. We start with a validation of the intrinsic density calculation andparameterise the bilinear approximation for varying resolutions. Next, the novel surface fluxforms are shown to be exactly conservative, before being compared to the volume averagedexpressions. The full balance of the surface flux contributions is shown next, highlighting30he importance of surface evolution in the equations of motions and showing the normalpressure should be constant over the interface. Finally, the calculation of the surface tensionusing the surface pressure is discussed.A study of bilinear resolutions is shown in Figure 4, with the density profile used toassess the accuracy of the varying surface approximations. This is calculated as described insection IV B by mapping the molecular positions based on the intrinsic interface and thenbinning molecules to the appropriate volumes with width ∆ z = 0 . λ r = λ u = σ and compared to theresults of Chac´on and Tarazona included as crosses in Fig 4 a ). This surface with λ r = σ is then sampled to define a piecewise bilinear approximation of the surface. The cases with5 and 10 bilinear bins per σ are indistinguishable from the intrinsic surface and are omittedin Fig 4 a ). The case where each σ unit has 2 . r = 0 . 4, is shown as a red line,where despite some decrease in the peaks at zero, the density profile is largely identical. Asincreased sampling requires large memory requirements and slows calculation, the optimalcase should provide good agreement for minimal resolution. At 1 . 25 bins per σ shown ingreen on Figure 4 a ), the sampling clearly gives a smeared density peak and is deemed notsufficiently accurate to provide a good representation of the interface. The case of 0 . 625 binsis also shown in light blue with the intrinsic interface for λ r = 4 σ included for comparisonin dark blue. This highlights that poor binsize Resolution has the same effect as a lowerwavelength fitting of the intrinsic surface. The L Error are shown in Figure 4 b ) defined asthe absolute sum of the density binning obtained for a given resolution, minus the λ r = σ case at every bin between − r , whilethe red points are different intrinsic minimum wavelength values, λ r , including λ r = 2 σ and λ r = 4 σ giving a blurred density from underfitting as well as the case with λ r = 0 . σ which,perhaps counter intuitively, gives a less sharp density profile due to overfitting. The pointsfor varying intrinsic and binsize resolution in Figure 4 b ) are reasonably well fitted by thelines which are of the form A /Resolution and B /Resolution with A = 6 . 42 and B = 2 . r < . − − z . . . . . . ρ Resolution0246810 L E rr o r a ) b )FIG. 4. The effect on intrinsic density of the bilinear approximation is shown for different resolutionsdefined as Resolution = σ/ ∆ r which is the number of times the square bilinear patches of sidelength ∆ r fit in an intermolecular spacing σ and for intrinsic surfaces the number of wavelengths λ r per σ with Resolution = σ/λ r . In a ) the full intrinsic solution is shown as a yellow line forthe case where λ r = λ u = σ compared to yellow crosses showing the solution from Chac´on andTarazona . Bin Resolution of 10 and 5 are omitted as indistinguishable from the yellow line, redis a Resolution of 2 . 5, green is 1 . 25 and dark blue is 0 . 63. The intrinsic interface with λ r = 4 isshown for comparison in light blue. In b ) the L Error = R − ρ ∆ r = σ ( z ) − ρ ∆ r ( z ) dz obtained from theintegral over z from − r = { . , . , . , . , . } and intrinsic (red) cases λr = { . , . , . } are shown, where the λ u = 1 . A /Resolution with A = 6 . 42 for binsizes and B /Resolution with B = 2 . 46 for theintrinsic wavelength. r < . 4. The effect of changing bilinear resolution can be seen to be similar to usinga larger minimum wavelength when fitting the intrinsic surface, with roughly a two folddifference in trends ( A / B ≈ . 6) so we need twice the increase in ∆ r resolution to matcha change in λ r . In this work, the presented plots use ∆ r = 0 . σ or Resolution = 5 forthe bilinear approximation to ensure measure quantities such as pressure are free from anyartefacts.Having parameterised the effect of bilinear resolution, we move on to checking exactconservation for an arbitrary volume consisting of 10 × 10 bilinear segments of size ∆ r = 0 . ∼ × x and y , centred on the intrinsic surface molecules,so it is the flux either side of the interface at ± ∆ z/ ± . Advection ), both from molecules movingover the surface in red and the surface moving over molecules in yellow, plus the interactionforce ( F orcing ) in green between the molecules all add up to the blue line for the changeof momentum in the control volume ( Accumulation ). The Accumulation is shown in Fig 5as a filled area under a curve to emphasise the fact that it is the integral of this area thatdetermines the changing momentum inside the control volume. This average momentum isconstantly changing due to forces acting on the volume as well as molecules entering andleaving. Exact conservation allows us to be sure that the implementation is correct, mollifiesthe issues associated with the non-uniqueness of the pressure tensor, providing an exact linkbetween surface pressure on the volume and momentum change inside, as well as guaranteeingall possible terms which could contribute to surface tension have been included. This exactconservation is valid for any arbitrary volume in space and is checked during all calculationsin this work, providing a thorough validation of the mathematics and the implementation ofthe ray tracing on a complex moving surface. The contributions due to forcing can be seento appear as a continuous line in Fig 5, as each molecule in the interface interacts with aneighbourhood of surrounding molecules within distance r c . The sum of all interactions onall molecules in the control volume, as well as between molecules on either side of the ∼ × IG. 5. Components of momentum labelled in Eq. (47) for a control volume of width ∆ z = 0 . x and y extents of 3 . 97, following the intrinsic surface with Advection including fluxesof molecules over the surface (red) and molecule crossings due to surface movement (yellow), F orcing the contributions due to forces between molecules crossing the volume surface (green)and Accumulation the resulting change in momentum inside the control volume which is exactlyequal to F orcing + Advection ( checked to machine precision). molecules entering or leaving the volume represents an evaporation or condensation event.These occur both due to molecular motion and when the fitting process of the interface findsa closer molecule than the current set, which can manifest in either the interface moving pastthe molecule or the molecule moving over the interface which is no longer following. Moleculescan also leave the volume in Fig 5 by diffusing along the intrinsic surface and leaving the 4 by34 control volume. This detailed balance has potential applications in designing improvementsto the intrinsic fitting process, for example to minimise evaporation events or track surfacetransport. The intrinsic fitting to a target density ρ s is equivalent to the mass control volumeEq. (21) having zero Accumulation and therefore ensuring zero net Advection for a controlvolume the size of the whole interface.The new surface flux forms of pressure are compared to the volume average (VA) formsEq. (44) in Figure 6. The volume average pressures are shown as points and surface pres-sures are shown by lines with kinetic pressure in red, tangential configurational pressure, Surf Π cT ≡ (cid:20) Surf Π cxx + Surf Π cyy (cid:21) , shown in yellow and surface normal pressure Surf Π cN ≡ Surf Π czz in blue.The kinetic pressure term of Eq. (46) is split into a surface evolution component ∂ξ/∂t andthe remaining kinetic term denoted with a prime, Surf Π k ′ N , so that, Surf Π kN = Surf Π k ′ N + ∂ξ∂t (56)and convection is assumed to be zero ρ u u z = 0. The kinetic pressure calculated usingthe surface pressure definition is visually identical to the volume average one in Figure 6.Note the normal component of surface kinetic pressure is shifted by ∆ z/ 2, as cell surfacespressure is obtained on surfaces while volume average pressure is at the cell centres. Thedensity is shown in light grey for comparison, to highlight the kinetic pressure contributionsare due to kinetic energy of molecules so directly correlated with their locations and identicalin both surface and volume average measures. The similarity in tangential components ofconfigurational pressure, shown by the yellow points and lines, is consistent with past work ,which shows VA and MOP give identical results in the limit of small bins. As the tangentialcontributions are calculated on a flat surface, using the form of Eq. (45) they are expectedto be identical to the VA pressure expressions. The normal components shown in blue isthe only quantity that is seen to be different between surface and volume average pressure.This is a direct result of the form of Eq. (46) which calculates the pressure dotted with thesurface normal ˜ n z at the location of every interaction for the interface at every time. Incontrast, the volume average tensor f zij r zij is independent of the surface normal so remainsaligned with the z Cartesian axis. This also makes clear why the tangential componentsare the same, the surface pressure is shown to be dotted with the n y in the flat surface35 − − z − . − . − . . . . P r e ss u r e a nd D e n s i t y FIG. 6. A comparison of pressure as a function of z position, including surface tangential pressureEq. (45) and surface normal pressure Eq. (46) (lines) compared to the volume average (points)pressure Eq. (44) from previous work . The kinetic pressure Surf Π k ′ N and Surf Π kT are shown as red lines,where the prime on the normal kinetic component denotes it does not include the surface movement, ∂ξ/∂t , term (where kinetic normal Surf Π k ′ N has an identical profile to the tangential component). Redpoints are VA Π kT and density of particles is shown as a grey line for reference. The tangentialconfigurational pressure Surf Π cT is shown as a yellow line with VA Π cT as yellow points and the normalcomponent of configurational pressure Surf Π cN is shown as a blue line with blue points for the VAterm VA Π cN . The position of the interface is plotted as a black line with a semi-transparent maskregion to hide the peak at the intrinsic interface. − − z − . − . . . . . C V P r e ss u r e FIG. 7. All terms which contribute to the normal component of the total CV Pressure includingconfigurational pressure Surf Π cN shown in blue, kinetic pressure Surf Π k ′ N shown in red and surface move-ment ∂ξ/∂t shown in green, with the total kinetic contribution Surf Π kN = Surf Π k ′ N + ∂ξ/∂t shown inyellow. The total pressure Surf Π N = Surf Π kN + Surf Π cN , black line, is a small constant value over the surfaceas required to ensure the surface is not moving. case which aligns with the volume average f yij r yij components so giving identical results. Asimilar observation is true for the other tangential component dotted with n x .Using surface pressure normal to the instantaneous surface can be shown to be essentialto the exact balance of momentum over the surface. To see this, in Figure 7 all contributionsto the surface in Eq. (46) are plotted on the same graph. The configurational (blue, Surf Π cN )and kinetic (red, Surf Π k ′ N ) pressure terms are identical to the ones plotted on figure 6, but when37he surface fluctuations (green, ∂ξ/∂t ) are added to give total kinetic pressure (yellow, Surf Π kN ),it can be seen that the resulting profile perfectly mirrors the configurational pressure. Theresult is the sum of the extended set of kinetic terms and configurational pressure givesa perfectly flat normal pressure over the surface, a required result for the interface to bestationary. Demonstration of a constant pressure profile near a wall in a molecular systemwas shown to be an important reason for using the VA or MOP form of pressure instead ofthe virial or IK1 expressions .The shape of the configurational pressure profile Surf Π cN has a flat region from z = 0 . . z ≈ 1. The total kinetic part, Surf Π kN , including ∂ξ/∂t , inyellow exactly mirrors the configurational with this same flat region. It can be seen thisflat region, that is, Surf Π k ′ N + ∂ξ/∂t = Constant from z = 0 . . z direction. This is a consequence of the intrinsic surface fitting process, bychoosing our reference frame to follow the interface molecules the measured kinetic pressureexposes the liquid structure in a similar way to the radial distribution function, in particularthe liquid’s tendancy to separate between molecular layers results in a drop to zero in thegap between the interface and first fluid layer. The interface movement term ∂ξ/∂t capturesthe movement of the intrinsic interface, a reference frame which tracks the surface molecules,and that allows the plot to identify the liquid structure peaks observed Put another way, theΠ k ′ term captures the structure inside the moving liquid cluster interface, the ∂ξ/∂t captureshow that structure moves.It is worth noting that Figure 7 shows the common equilibrium assumption ∇ · Π = 0applied in statistical mechanical derivation of surface tension is only valid if the contri-bution from surface movement and local curvature are correctly included. The VA terms inFig 6 will not satisfy this condition due to missing curvature term shown in Appendix C.This may suggest the equilibrium assumption is a source of error for spherical and cylindricalvolumes. A full mechanical approach may provide insight by including kinetic fluxes, surfacemovement and changing momentum of the volume as shown in Figure 5.Finally, we consider the surface tension in Figure 8. In order to explore the various38 − − z − . − . . . . . . . . P r e ss u r e a ndSu r f a ce T e n s i o n FIG. 8. Calculation of the surface tension with normal minus tangential pressure Π N − Π T as redlines for the full VA pressure, yellow lines for the surface configurational part only and blue linesthe entire surface pressure including the surface evolution, where the minimum at z = 0 of theyellow and blue lines are not shown as they go to − . − . ( ρ l + ρ g ). The corresponding cumulative integral of each curve to agiven z value, γ ( z ) = R z − [Π N ( z ′ ) − Π T ( z ′ )] dz ′ , is shown in the same colour, with red circles VApressure, yellow circles configurational surface, blue circles the full surface pressure and green thefixed grid MOP pressure. style formula are plotted, VA γ = Z z − (cid:20) VA Π N − VA Π T (cid:21) dz ′ (57) Surf γ c ( z ) = Z z − (cid:20) Surf Π cN − Surf Π cT (cid:21) dz ′ (58) Surf γ ( z ) = Z z − (cid:20) Surf Π N − Surf Π T (cid:21) dz ′ ≈ Surf γ c ( z ) + Z z − ∂ξ ( z ′ ) ∂t dz ′ (59)which includes the surface tension from the volume average pressure Eq. (57), the config-uration part of the surface pressure Eq. (58) and from the full surface pressure Eq. (59).In previous work , the contributions to surface tension from the volume average terms werediscussed, and the red lines and circles in Figure 8 are presented as a basis for comparison,together with the method of planes pressure measurements obtain using a fixed grid shownby green lines and circles. The configurational part is important as the normal component ofconfigurational surface pressure was the only term which showed a difference when comparedto the volume average pressure in Figure 6. As the kinetic normal and tangential componentsare identical for both volume average and surface pressure these cancel in the Kirkwood andBuff formula and the configurational terms contain all contributions to surface tension.As a result, any difference between the surface tension calculated from the VA pressure andsurface pressure would be expected to come from the Surf Π CN term. The red and yellow linesin Figure 8 show the difference between these two pressure measurements is mostly locatedbetween the intrinsic interface and first layer in the liquid region from z = 0 to 1. Theresulting integral to give surface tension from Eq. (58) (yellow circles) can be seen to stillconverge to the same value as the volume average Eq. (57) (red circles) in Figure 8 (and thesurface tension obtained from the fixed references frame pressure shown by green circles).However, the surface pressure sees a greater contribution to surface tension in the liquidregion, from z = − . z = 0, than the VA. This is offset by a much larger negative con-tribution from the intrinsic surface molecules themselves at z = 0. This is more pronouncedfor the full surface pressure, the blue line, and resulting surface tension using Eq. (59),shown by blue circles in Figure 8. There is an almost linear contribution from z = − z = 0 and an even larger negative contribution from the interface molecules at z = 0 when40ompared to the configurational part of Eq. (58). The interface control volume sits on theinterfacial molecules themselves, so this peak represents the tangential forces between them,as well as a contribution due to the movement of the surface. The key difference betweenthe configurational Eq. (58) and the full surface pressure Eq. (59) is the inclusion of surfaceevolution in time ∂ξ/∂t . This can be seen to give a net zero contribution to surface tension inthis equilibrium case, but redistributes giving a large negative contribution from the surfaceitself at z = 0 and equal positive contribution between z = ± ξ ( x, y, t ), withthe control volume formulation providing exact conservation every single timestep. Thisexact momentum balance on a complex time evolving interface could provide useful insightsin bubble nucleation, film rupture or contact line dewetting. VI. CONCLUSIONS In this work, we derive a new formulation of surface pressure in a reference frame movingwith the interface between a liquid and a vapour. This derivation starts from statisticalmechanical definitions of density, momentum and energy equations but without ensembleaverage. The resulting equations are integrated over a volume fitted to the interface as itevolves in time, giving instantaneous surface fluxes on the moving surface. These surfaceflux equations are shown to include two extra contributions, one due to the instantaneoussurface curvature at the point of surface crossings and one due to the movement of the surfacein time. By including all contributions for curvature and surface movement, the equationscan be shown to be exactly conservative in a molecular dynamics (MD) simulation. Theseextra terms are also shown to be essential to provide an exact force balance over the movingliquid-vapour interface. The derived equations are presented in a form which can be easily41mplemented in MD, applying a ray tracing approach commonly used for computer graphics.Several insights into the pressure and surface tension are presented for the simplest case ofa flat interface between the liquid and vapour. Despite being tested in a simple system, thederived equations make no assumption other than mass conservation and Newton’s laws.As a result, they provide exact expressions for the equations of hydrodynamics which cantake any local region on the surface and follow the instantaneous surface shape as it evolves.Therefore, the new equations are expected to have great potential applications to understanda range of hydrodynamics phenomena including Marangoni effects, growing bubbles, movingcontact lines and deforming interfaces. VII. DATA AVAILABILITY STATEMENT The data that support the findings of this study are available from the correspondingauthor upon reasonable request. Appendix A: Integrating the Equations This appendix details the full mathematics of the process used to provide the equationsfor use in a molecular dynamics simulation. First to highlight the similarity between theconfigurational and kinetic part we integrate the kinetic expression over time in Eq. (39).This can be interpreted in two ways, either as part of time averaging or as part of theevolution of the system, Z t t (cid:2) ρ u u z + Π kz (cid:3) dt = 1∆ S z N X i =1 m i ˙ r i Z t t h ˙ x i ∂ξ i + ∂x i + ˙ y i ∂ξ i + ∂y i + ˙ z i + ∂ξ i + ∂t i dS + zi dt Π cz = 12 1∆ S z N X i,j f ij Z h x ij ∂ξ + λ ∂x λ + y ij ∂ξ + λ ∂y λ + z ij i dS + zλ dλ. (A1)42he expression for pressure in Eq. (A1) can then be written as follows, Z t t (cid:2) ρ u u z + Π kz (cid:3) dt = 1∆ S z N X i =1 m i ˙ r i Z h dx τ dτ ∂ξ i + ∂x i + dy τ dτ ∂ξ i + ∂y i + dz τ dτ i dS + zi dτ + ϑ t ! Π cz = 12 1∆ S z N X i,j f ij Z h ∂x λ ∂λ ∂ξ + λ ∂x λ + ∂y λ ∂λ ∂ξ + λ ∂y λ + ∂z λ ∂λ i dS + zλ dλ, (A2)where we use the substitution τ = [ t − t ] / [ t − t ] in the integral over time to get the same formas the interaction over path λ , denoting r ≡ r i ( t ) and r ≡ r i ( t ) so r i ( t ) = r τ = r + τ r with r = r − r and τ from 0 to 1. Now both expressions are in the form r α = r + s r with α = { τ, λ } and s = { τ, λ } allowing us to write the generalised control volume function, ϑ s = (cid:2) H (cid:0) x + − x α ( s ) (cid:1) − H (cid:0) x − − x α ( s ) (cid:1)(cid:3) × (cid:2) H (cid:0) y + − y α ( s ) (cid:1) − H (cid:0) y − − y α ( s ) (cid:1)(cid:3) × (cid:2) H (cid:0) ξ i + − z α ( s ) (cid:1) − H (cid:0) ξ i − − z α ( s ) (cid:1)(cid:3) = Λ x ( s )Λ y ( s ) ˜Λ z ( s ) , (A3)where each directional difference between two Heaviside functions is written using shorthandΛ β with β = { x, y, z } and just the functional dependence on s shown. The tilde on the Λ z term denotes the intrinsic interface component is included on this surface. Both integralexpressions in Eq. (A2) are in the form of an integral of the derivative of Eq. (A3) withrespect to s . More generally, the total expression for all surfaces can be written concisely as, Z ∂ϑ s ∂s ds = Z ∂∂s h Λ x ( s )Λ y ( s ) ˜Λ z ( s ) i ds = Z " ∂ Λ x ∂s Λ y ˜Λ z + Λ x ∂ Λ y ∂s ˜Λ z + Λ x Λ y ∂ ˜Λ z ∂s ds, (A4)for example, with s = λ the last term describes the difference between top and bottomconfigurational term in z , Z Λ x Λ y ∂ ˜Λ z ∂s ds = Z h ∂x λ ∂λ ∂ξ + λ ∂x λ + ∂y λ ∂λ ∂ξ + λ ∂y λ + ∂z λ ∂λ i dS + zλ − h ∂x λ ∂λ ∂ξ - λ ∂x λ + ∂y λ ∂λ ∂ξ - λ ∂y λ + ∂z λ ∂λ i dS − zλ dλ. (A5)For the x and y surfaces, the intersection of a plane and line can be obtained directly,considering the first term on the right-hand side of Eq. (A4), Z ∂ Λ x ∂s Λ y ˜Λ z ds = Z dx α ds (cid:2) δ (cid:0) x + − x α ( s ) (cid:1) − δ (cid:0) x − − x α ( s ) (cid:1)(cid:3) Λ y ( s ) ˜Λ z ( s ) ds. (A6)43e apply the property of the Delta function to express it as the sum of its roots, δ (cid:0) x + − x α (cid:1) = N roots X k =1 δ ( s − s k ) | dx α ( s k ) /ds | , (A7)so positive surface in x can then be expressed, Z ∂x α ( s ) ∂s δ (cid:0) x + − x α ( s ) (cid:1) Λ y ( s ) ˜Λ z ( s ) ds = Z ∂x α ( s ) ∂s N roots X k =1 δ ( s − s k ) | dx α ( s k ) /ds | Λ y ( s ) ˜Λ z ( s ) ds = ∂x α ( s k ) /∂s | ∂x α ( s k ) /∂s | | {z } i ) Direction of crossing N roots X k =1 [ H (1 − s k ) − H ( − s k )] | {z } ii ) Crossing between limits Λ y ( s k ) ˜Λ z ( s k ) | {z } iii ) Crossing on yz CV Surface . (A8)The three annotated terms include i ) a signum function which determines the direction ofcrossing with ∂x α ( s k ) /∂s | ∂x α ( s k ) /∂s | = x | x | = n x · r | n x · r | , (A9)using dx α /ds = x so i ) is expressed in terms of the normal to the x surface n x = [1 , , r and r . The expression for this root s k on the flatsurfaces is the intersection of a plane and a line , obtained by equating the surface x + andline x + sx to solve for the value of s at crossing s k = ( x + − x ) /x ≡ x + k allowing thecrossing between limits term of Eq. (A8) to be expressed as,[ H (1 − s k ) − H ( − s k )] = H (cid:18) x + − x x (cid:19) − H (cid:18) x − x + x (cid:19) = 12 sgn (cid:18) x (cid:19) (cid:2) sgn ( x + − x ) − sgn ( x − x + ) (cid:3) , (A10)which is the expression found in the method of planes form of stress and determines if r and r are on opposite sides of the plane . Finally, iii) the location of the crossing s k on the flat x + surface gives a value of zero for Λ( s k ) when not within the limits of the control volumesurface in the y and z directions, where the z CV surface moves as the intrinsic surfacesmoves, shown in Fig 3. The location of crossing in each direction is, x α ( s k ) = x + ; y α ( s k ) = y + s k y ; z α ( s k ) = z + s k z , (A11)44nd the ˜Λ z ( s k ) function is,˜Λ z ( s k ) = H (cid:0) z + + ξ (cid:0) x + , y + y s k (cid:1) − z − z s k (cid:1) − H (cid:0) z − + ξ (cid:0) x + , y + y s k (cid:1) − z − z s k (cid:1) . (A12)The stress on the x + surface is therefore, Z t t (cid:2) ρ u u x + Π kx (cid:3) dt = 1∆ S x N X i =1 m i ˙ r i r i · n x | r i · n x | (cid:20) H (cid:18) x + − x i x i (cid:19) − H (cid:18) x i − x + x i (cid:19)(cid:21) Λ y ( t k ) ˜Λ z ( t k ) Π cx = 12 1∆ S z N X i,j f ij r ij · n x | r ij · n x | (cid:20) H (cid:18) x + − x j x ij (cid:19) − H (cid:18) x i − x + x ij (cid:19)(cid:21) Λ y ( λ k ) ˜Λ z ( λ k ) , note the inclusion of the index i so r → r i to emphasise it is molecule i which is evolvingfrom time t to t The link to the common method of planes (MOP) expression in theliterature can be seen using x/ | x | = sgn ( x ) and Eq. (A10) to give, r · n x | r · n x | (cid:20) H (cid:18) x + − x x (cid:19) − H (cid:18) x − x + x (cid:19)(cid:21) = 12 sgn ( x ) sgn (cid:18) x (cid:19) (cid:2) sgn ( x + − x ) − sgn ( x − x + ) (cid:3) = 12 (cid:2) sgn ( x + − x ) − sgn ( x − x + ) (cid:3) (A13)so the expressions for stress on the x surface are the well know MOP form, Z t t (cid:2) ρ u u x + Π kx (cid:3) dt = 12 1∆ S x N X i =1 m i ˙ r i (cid:2) sgn ( x + − x i ) − sgn ( x i − x + ) (cid:3) Λ y ( t k ) ˜Λ z ( t k ) Π cx = 14 1∆ S z N X i,j f ij (cid:2) sgn ( x + − x j ) − sgn ( x i − x + ) (cid:3) ]Λ y ( λ k ) ˜Λ z ( λ k ) . (A14)Obtaining the expression for crossing on the other flat surfaces x − and y ± follow theprocess just outlined. It is complicated in the z direction by the intrinsic surface, which wecan evaluate as follows. The z surface is the last term on the right of Eq. (A4), Z Λ x Λ y ∂ ˜Λ z ∂s ds = Z Λ x ( s )Λ y ( s ) ∂∂s (cid:20) H (cid:0) z + − ξ ( x α ( s ) , y α ( s )) − z α ( s ) (cid:1) − H (cid:0) z − − ξ ( x α ( s ) , y α ( s )) − z α ( s ) (cid:1) (cid:21) ds = Z ∂ ( ξ − z α ) ∂s (cid:2) δ (cid:0) z + − ξ − z α (cid:1) − δ (cid:0) z − − ξ − z α (cid:1)(cid:3) Λ x Λ y ds. (A15)45gain, using the property of the Delta function to express the sum of the roots of inter-section, the crossings of the intrinsic interface by the line of interaction, δ (cid:0) z + − ξ − z α (cid:1) = N roots X k =1 δ ( s − s k ) | dds ( ξ ( s k ) − z α ( s k )) | (A16)allows Eq. (A15) to be written as, Z ∂ ( ξ − z α ) ∂s δ (cid:0) z + − ξ − z α (cid:1) Λ x Λ y ds = Z ∂∂s ( ξ ( s ) − z α ( s )) N roots X k =1 δ ( s − s k ) | ∂∂s ( ξ ( s k ) − z α ( s k )) | Λ x ( s )Λ y ( s ) ds = ∂∂s ( ξ ( s k ) − z α ( s k )) | ∂∂s ( ξ ( s k ) − z α ( s k )) | | {z } i ) Direction of crossing N roots X k =1 [ H (1 − s k ) − H ( − s k )] | {z } ii ) Crossing between limits Λ x ( s k )Λ y ( s k ) . | {z } iii ) Crossing on xy CV Surface . (A17)The expression is broadly similar to the flat surface case of Eq. (A8), with a more complexexpression for the direction of crossing in terms of the intrinsic surface i). This can beunderstood using a change of variable, ∂∂s ( ξ − z α ) = ∂ r α ∂s · ∂∂ r α ( ξ − z α ) = r · ∇ α ( ξ − z α ) , (A18)which allows term i) to be written as, ∂∂s ( ξ − z α ) | ∂∂s ( ξ − z α ) | = r · ∇ α ( ξ − z α ) | r · ∇ α ( ξ − z α ) | || ∇ α ( ξ − z α ) |||| ∇ α ( ξ − z α ) || = r · ˜ n z | r · ˜ n z | , (A19)with || a || denoting vector magnitude, which must be positive so can be moved inside theabsolute value on the denominator allowing the expression to be written in terms of thenormal to the intrinsic surface, ˜ n z ≡ ∇ α ( ξ − z α ) || ∇ α ( ξ − z α ) || , (A20)this can be seen to be in the same form as the flat surface, obtaining the direction of thevector between r and r dotted with the surface normal. The crossing between limits ofii) and surface area term iii) are identical in form to the Eq. (A8) case, although s k can no46onger be obtained as a closed form solution as in the flat surface case of Eq. (A10). Thesolution of line and surface z i + s k z ij = ξ ( x, y ) to obtain s k requires the application of a rootfinding process such as Newton Raphson.So, the implementation of Eq. (A1) is therefore, Z t t ρ u u z + Π kz dt = 1∆ S z N X i =1 m i ˙ r i r i · ˜ n z | r i i · ˜ n z | N roots X k =1 [ H (1 − t k ) − H ( − t k )] Λ x ( t k )Λ y ( t k )+ 1∆ S z N X i =1 m i ˙ r i ϑ t Π cz = 12∆ S z N X i,j f ij r ij · ˜ n z | r ij · ˜ n z | N roots X k =1 [ H (1 − λ k ) − H ( − λ k )] Λ x ( λ k )Λ y ( λ k ) , (A21)where again the use of r i to emphasise this is per molecule. As the surface is no longerflat, the roots t k and λ k in this expression must be obtained using a form of algorithmic raytracing, discussed in more detail in the body of the text. Appendix B: Linear Algebra The least square minimisation of the weighting function to fit the intrinsic surface hasbeen discussed in the literature, most extensively in the supplementary material of Longfordet al. . However, a slightly different notation is given here to clarify a few steps, startingfrom, W ( a µν ) = 12 N p X p =1 " z p − k u X µ = − k u k u X ν = − k u a µν ( t ) f µ ( x p ) f ν ( y p ) + ψ ˜ A. (B1)47his can be seen to be a linear algebra problem by defining the vector f ( x, y ) where eachrow is a wavevector, f ( x, y ) = f − k u ( x ) f − k u ( y ) f − k u ( x ) f − k u +1 ( y ) . . .f − k u ( x ) f k u ( y ) f − k u +1 ( x ) f − k u ( y ) . . .f k u ( x ) f k u ( y ) (B2)so the total matrix F can be defined by stacking matrices f T for all pivot locations, F = f T ( x , y ) f T ( x , y ) . . . f T ( x N p , y N p ) (B3)and defining z = [ z , z , ..., z N p ] T , a = a νµ and ˜ A = a T a B with B = 4 π diag( ν + µ ) amatrix with values only on the diagonal and zero elsewhere. So, W in Eq. (53) becomes, W ( a ) = 12 || z − F a || + ψ ˜ A , (B4)and the minimum is ontained by setting the derivative with respect to a to zero, ∇ a W = − F T ( z − F a ) + ψ ∇ a ˜ A = 0 , (B5)and ∇ a ˜ A = B a so the optimal value of a is obtained by solving this equation, a = ( F T F − ψ B ) − F T z . (B6)To maximise efficiency, LAPACK is used to solve this equation. In practice this means thematrix F of size M × N p with M = 4 k u wavenumbers and N p pivot positions, is multipliedby its transpose (using e.g. Lapack DGEMM). The constraint is then applied by subtractinga matrix of size M × M with just diagonal elements B = 4 π ψ [ µ + ν ] that are non-zero.This can then be used in a linear algebra solver (e.g. Lapack DGESV) with right-hand side F T z to get the values of a . 48 ppendix C: A Derivation of the Volume Average Form In this section, we consider how to obtain the volume average (VA) expressions from Eq.(25), ddt Z V ρ u dV = ddt N X i =1 m i ˙ r i ϑ i = N X i =1 m i ˙ r i dϑ i dt + N X i =1 m i ¨ r i ϑ i , (C1)where the continuum left-hand side is expressed using the divergence theorem, ddt Z V ρ u dV = − I S [ ρ uu + Π ] · d S = − Z V ∂∂ r · [ ρ uu + Π ] dV, (C2)while the right-hand side is expanded using Eq. (14), Eq. (27) and Eq. (28), ddt N X i =1 m i ˙ r i ϑ i = N X i =1 m i ˙ r i (cid:18) ˙ r i · ∂ϑ i ∂ r i + ∂ξ∂t ∂ϑ i ∂ξ (cid:19) + N X i,j f ij Z r ij · ∂ϑ λ ∂ r λ dλ, (C3)by assuming ∂ϑ i /∂ r i = − ∂ϑ i /∂ r , ∂ϑ λ /∂ r λ = − ∂ϑ λ /∂ r and ∂ξ/∂t = 0, we can expresseverything as a derivative in terms of r , Z V ∂∂ r · [ ρ uu + Π ] dV = ∂∂ r · " N X i =1 m i ˙ r i ˙ r i ϑ i + N X i,j f ij r ij · Z ϑ λ dλ , (C4)and so, comparing the expressions inside the derivative and assuming an average value forthe volume yields the VA pressure given in previous work, ρ uu + VA Π = 1∆ V " N X i =1 m i ˙ r i ˙ r i ϑ i + 12 N X i,j f ij r ij Z ϑ λ dλ . (C5)However, the assumptions ∂ϑ α /∂ r α = − ∂ϑ α /∂ r and ∂ξ/∂t = 0 are not valid, as will beshown here for the configurational term, resulting in an extra term in the VA expressionfor a volume between curved interfaces. Pressure is defined in the Irving and Kirkwood process by collecting terms inside the gradient with respect to r and comparing forms tothe continuum expression ∂/∂ r · Π . In the pointwise Irving and Kirkwood , this uses theproperty of the Dirac delta ∂/∂ r λ δ ( r − r λ ) = − ∂/∂ r δ ( r − r λ ). For an integrated controlvolume between intrinsic surfaces, to see if the same process can be applied, we compare thederivatives, ∂ϑ λ ∂ r λ and ∂ϑ λ ∂ r . Starting with the derivative of ϑ λ with respect to r , ∂ϑ λ ∂ r = i ∂ Λ x ∂x Λ y ˜Λ z + j Λ x ∂ Λ y ∂y ˜Λ z + k Λ x Λ y ∂ ˜Λ z ∂z , (C6)49here e.g. ∂ Λ x ∂x = δ ( x + − x λ ) − δ ( x − − x λ ) and ∂ ˜Λ z ∂z = δ ( ξ + λ − z λ ) − δ ( ξ − λ − z λ ) as ξ ± λ = z ± ∆ z + ξ so dξ ± λ /dz = 1.Next, consider the derivative of ϑ λ with respect to r λ ∂ϑ λ ∂ r λ = i ∂ Λ x ∂x λ Λ y ˜Λ z + j Λ x ∂ Λ y ∂y λ ˜Λ z + Λ x Λ y " i ∂ ˜Λ z ∂x λ + j ∂ ˜Λ z ∂y λ + k ∂ ˜Λ z ∂z λ . (C7)Noting that ∂ Λ x ∂x λ = − [ δ ( x + − x λ ) − δ ( x − − x λ )] = − ∂ Λ x ∂x and similar for y , while the deriva-tive of ˜Λ z with respet to z is ∂ ˜Λ z ∂z λ = δ ( ξ + λ − z λ ) − δ ( ξ − λ − z λ ), so ∂ ˜Λ z ∂z = − ∂ ˜Λ z ∂z λ . Using theseequivalences, ∂ϑ λ ∂ r λ = − i ∂ Λ x ∂x Λ y ˜Λ z − j ∂ Λ y ∂y Λ x ˜Λ z − k ∂ ˜Λ z ∂z Λ x Λ y + Λ x Λ y " i ∂ ˜Λ z ∂x λ + j ∂ ˜Λ z ∂y λ = − ∂ϑ λ ∂ r + (cid:20) i (cid:18) ∂ξ + λ ∂x λ dS + zλ − ∂ξ − λ ∂x λ dS − zλ (cid:19) + j (cid:18) ∂ξ + λ ∂y λ dS + zλ − ∂ξ − λ ∂y λ dS − zλ (cid:19)(cid:21) , (C8)where ∂ ˜Λ z ∂x λ = (cid:20) ∂ξ + λ ∂x λ δ ( ξ + λ − z λ ) − ∂ξ − λ ∂x λ δ ( ξ − λ − z λ ) (cid:21) , (C9)and similar for the y derivative, with the surface notation used, dS ± zλ = δ ( ξ ± λ − z λ )Λ x Λ y . Thefull expression is therefore, N X i,j f ij [ ϑ i − ϑ j ] = N X i,j f ij r ij · Z ∂ϑ λ ∂ r λ dλ = − ∂∂ r · N X i,j f ij r ij Z ϑ λ dλ | {z } VA Π c + N X i,j f ij r ij · Z (cid:20) ∂ξ + λ ∂ r λ dS + zλ − ∂ξ − λ ∂ r λ dS − zλ (cid:21) dλ | {z } Extra Term , (C10)using ∂ξ ± λ /∂z λ = 0 to write the extra term in vector form. REFERENCES R. Evans. The nature of the liquid-vapour interface and other topics in the statisticalmechanics of non-uniform, classical fluids. Advances in Physics , 28(2):143–200, 1979. ISSN0001-8732. 50 J S Rowlinson and B Widom. Molecular Theory of Capillarity . Dover books on chemistry.Dover Publications, 2002. ISBN 9780486425443. E. N. Parker. Tensor Virial Equations. Phys. Rev. , 96:1686, 1954. J. H. Irving and J. G. Kirkwood. The Statistical Mechanics Theory of Transport Processes.IV. The Equations of Hydrodynamics. J. Chem. Phys. , 18:817, 1950. P. Schofield and J. R. Henderson. Statistical Mechanics of Inhomogeneous Fluids. Proc.R. Soc. Lond. A , 379:231, 1982. W. Noll. Die Herleitung der Grundgleichungen der Thermomechanik der Kontinua aus derStatistischen Mechanik. J. Ration. Mech. Anal. , 4:627, 1955. N. C. Admal and E. B. Tadmor. A Unified Interpretation of Stress in Molecular Systems. J. Elast. , 100:63, 2010. A. I. Murdoch. A Critique of Atomistic Definitions of the Stress Tensor. J. Elast. , 88:113,2007. J. F. Lutsko. Stress and elastic constants in anisotropic solids: Molecular dynamics tech-niques. J. Appl. Phys. , 64:1152, 1988. J. Cormier, J. M. Rickman, and T. J. Delph. Stress calculation in atomistic simulations ofperfect and imperfect solids. J. Appl. Phys. , 89:99, 2001. D. H. Tsai. The virial theorem and stress calculation in molecular dynamics. J. Chem.Phys. , 70:1375, 1978. B. D. Todd, D. J. Evans, and P. J. Daivis. Pressure Tensor for inhomogeneous fluids. Phys.Rev. E , 52:1627, 1995. M. Han and J. S. Lee. Method for calculating the heat and momentum fluxes of inhomo-geneous fluids. Phys. Rev. E , 70:061205, 2004. E. R. Smith, D. M. Heyes, D. Dini, and T. A. Zaki. Control-volume representation ofmolecular dynamics. Phys. Rev. E. , 85:056705, 2012. C. Hirsch. Numerical Computation of Internal and External Flows . Elsevier, Oxford, 2ndedition, 2007. Alexandr Malijevsk´y and George Jackson. A perspective on the interfacial properties ofnanoscopic liquid drops. Journal of Physics: Condensed Matter , 24(46):464121, oct 2012. David van Dijk. Comment on pressure enhancement in carbon nanopores: a major con-51nement effect by y. long, j. c. palmer, b. coasne, m. liwinska-bartkowiak and k. e. gubbins,phys. chem. chem. phys., 2011, 13, 17163. Phys. Chem. Chem. Phys. , 22:9824–9825, 2020. M. Zhou. A new look at the atomic level virial stress: on continuum-molecular systemequivalence. Proc. R. Soc. Lond. , 459:2347, 2003. I. Newton. Philosophiæ Naturalis Principia Mathematica . 3rd edition, 1726. A. Motte. The Mathematical Principles of Natural Philosophy . 1st edition, 1729. A. Harasima. Molecular theory of surface tension. Advances in Chemical Physics , 1:203–237, 1958. Bjørn Hafskjold and Tamio Ikeshoji. Microscopic pressure tensor for hard-sphere fluids. Phys. Rev. E , 66:011203, Jul 2002. Kaihang Shi, Yifan Shen, Erik E. Santiso, and Keith E. Gubbins. Microscopic pressuretensor in cylindrical geometry: Pressure of water in a carbon nanotube. Journal of ChemicalTheory and Computation , 0(0), 0. doi:https://dx.doi.org/10.1021/acs.jctc.0c00607. Yun Long, Jeremy C. Palmer, Benoit Coasne, Kaihang Shi, Magorzata liwiska Bartkowiak,and Keith E. Gubbins. Reply to the comment on pressure enhancement in carbonnanopores: a major confinement effect by d. van dijk, phys. chem. chem. phys., 2020,22, doi: 10.1039/c9cp02890k. Phys. Chem. Chem. Phys. , 22:9826–9830, 2020. Billy D. Todd and Peter J. Daivis. Nonequilibrium Molecular Dynamics: Theory, Algo-rithms and Applications . Cambridge University Press, 2017. E. Chac´on and P. Tarazona. Intrinsic Profiles beyond the Capillary Wave Theory: A MonteCarlo Study. Physical Review Letters , 91(16):166103, oct 2003. ISSN 0031-9007. L´ıvia B. P´artay, Gy¨orgy Hantal, P´al Jedlovszky, ´Arp´ad Vincze, and George Horvai. A newmethod for determining the interfacial molecules and characterizing the surface roughnessin computer simulations. application to the liquid–vapor interface of water. Journal ofComputational Chemistry , 29(6):945–956, 2008. Adam P. Willard and David Chandler. Instantaneous liquid interfaces. The Journal ofPhysical Chemistry B , 114(5):1954–1958, 2010. PMID: 20055377. Miguel Jorge, P´al Jedlovszky, and M. Nat´alia D. S. Cordeiro. A critical assessment ofmethods for the intrinsic analysis of liquid interfaces. 1. surface site distributions. TheJournal of Physical Chemistry C , 114(25):11169–11179, 2010.52 Marcello Sega, Sofia S. Kantorovich, P´al Jedlovszky, and Miguel Jorge. The generalizedidentification of truly interfacial molecules (itim) algorithm for nonplanar interfaces. TheJournal of Chemical Physics , 138(4):044110, 2013. M V Berry. The molecular mechanism of surface tension. Physics Education , 6(2):79, 1971. Marcello Sega, Bal´azs F´abi´an, George Horvai, and P´al Jedlovszky. How is the surfacetension of various liquids distributed along the interface normal? Journal of PhysicalChemistry C , 120(48):27468–27477, 2016. ISSN 19327455. J.P.R.B. Walton, D.J. Tildesley, J.S. Rowlinson, and J.R. Henderson. The pressure tensorat the planar surface of a liquid. Molecular Physics , 48(6):1357–1368, 1983. F. Varnik, J. Baschnagel, and K. Binder. Molecular dynamics results on the pressure tensorof polymer films. The Journal of Chemical Physics , 113(10):4444–4453, 2000. D. M. Heyes, E. R. Smith, D. Dini, and T. A. Zaki. The equivalence between VolumeAveraging and Method of Planes definitions of the pressure tensor at a plane. J. Chem.Phys. , 135:024512, 2011. Carlos Braga, Edward R. Smith, Andreas Nold, David N. Sibley, and Serafim Kalliadasis.The pressure tensor across a liquid-vapour interface. The Journal of Chemical Physics ,149(4):044705, 2018. J. S. Rowlinson. Thermodynamics of inhomogeneous systems. Pure and Applied Chemistry ,65(5), 1993. R. Delgado-Buscalioni, E. Chacon, and P. Tarazona. Hydrodynamics ofNanoscopic Capillary Waves. Physical Review Letters , 101(10):106102,sep 2008. ISSN 0031-9007. doi:10.1103/PhysRevLett.101.106102. URL http://link.aps.org/doi/10.1103/PhysRevLett.101.106102 . D J Evans and G Morriss. Statistical Mechanics of Nonequilibrium Liquids . Theoreticalchemistry. Cambridge University Press, 2008. ISBN 9781139471930. R. N. Bracewell. Fourier Transform and Its Applications . McGraw-Hill Companies, 3rdedition, 1999. ISBN 0073039381. M. C. Potter and D. C. Wiggert. Mechanics of Fluids . Brooks/Cole, California, 3rd edition,2002. D. M. Heyes, E. R. Smith, D. Dini, and T. A. Zaki. The method of planes pressure tensor53or a spherical subvolume. J. Chem. Phys. , 140(5):054506, 2014. John G. Kirkwood and Frank P. Buff. The statistical mechanical theory of surface tension. The Journal of Chemical Physics , 17(3):338–343, 1949. D. J. Evans and G. P. Morriss. Statistical Mechanics of Non-Equilibrium Liquids . Aus-tralian National University Press, Canberra, 2nd edition, 2007. Note1. Half of of the control volume is a tetrahedron, so the top and bottom surfaces givea different pressure tensor. In this work, one surface is no longer flat, instead following theintrinsic surface, which means we have departed from the Cauchy definition. E. R. Smith, P. J. Daivis, and B. D. Todd. Measuring heat flux beyond fouriers law. TheJournal of Chemical Physics , 150(6):064103, 2019. Shaun D. Ramsey, Kristin Potter, and Charles Hansen. Ray bilinear patch intersections. Journal of Graphics Tools , 9(3):41–47, 2004. W.C. Reynolds. Thermodynamic Properties in SI: Graphs, Tables, andComputational Equations for Forty Substances . Department of Mechani-cal Engineering, Stanford University, 1979. ISBN 9780917606052. URL https://books.google.co.uk/books?id=6nogAQAAIAAJ . Bo Shi. Molecular dynamics simulation of the surface tension and contact angle of argonand water . PhD thesis, University of California Los Angeles. USA, 2006. E. R. Smith, E. A. Mller, R. V. Craster, and O. K. Matar. A langevin model for fluctuatingcontact angle behaviour parametrised using molecular dynamics. Soft Matter , 12:9604–9615, 2016. E. R. Smith. On the Coupling of Molecular Dynamics to Continuum Computational FluidDynamics . PhD thesis, Imperial College London, 2014. Fernando Bresme, E. Chac´on, and P Tarazona. Molecular dynamics investigation of theintrinsic structure of water–fluid interfaces via the intrinsic sampling method. PhysicalChemistry Chemical Physics , 10(32):4676–7, aug 2008. Francis G. J. Longford, Jonathan W. Essex, Chris-Kriton Skylaris, and Jeremy G. Frey.Surface reconstruction amendment to the intrinsic sampling method. The Journal of Chem-ical Physics , 149(23):234705, 2018. E. Chac´on, E. M. Fern´andez, D. Duque, R. Delgado-Buscalioni, and P. Tarazona. Com-54arative study of the surface layer density of liquid surfaces. Phys. Rev. B , 80:195403, Nov2009. E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz,A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen.