The pressure tensor (equivalent to the negative stress tensor) at both microscopic and macroscopic levels is fundamental to many aspects of engineering and science, including fluid dynamics, solid mechanics, biophysics, and thermodynamics. In this Perspective, we review methods to calculate the microscopic pressure tensor. Connections between different pressure forms for equilibrium and nonequilibrium systems are established. We also point out several challenges in the field, including the historical controversies over the definition of the microscopic pressure tensor; the difficulties with many-body and long-range potentials; the insufficiency of software and computational tools; and the lack of experimental routes to probe the pressure tensor at the nanoscale. Possible future directions are suggested.

## I. INTRODUCTION

Pressure, *P*, defined as the average force *F* per unit area acting on a surface of area *S*, is one of the primary state variables, together with the temperature, *T*, and composition, that determine the thermodynamic properties of a homogeneous system of molecules at equilibrium. In such a system, the average force and the pressure are the same in all directions, *P* is a scalar and is well defined even at the microscale. Statistical mechanics informs us that *P* = *P*^{K} + *P*^{C}, where superscript *K* indicates the contribution due to the kinetic energy of the molecules and *C* indicates the configurational contribution, i.e., that from the intermolecular forces and any external field. For condensed phases, the configurational contribution to *P* is usually dominant, especially at low temperatures. For a perfect gas at equilibrium in the absence of an external field, or a real equilibrium gas at low enough density that the influence of intermolecular forces can be neglected, the configurational contribution is negligible, and $P=PK=nrkBT$, where $nr$ is the number density at position **r** and *k*_{B} is the Boltzmann constant. For such a gas, this relation holds even when the gas is nonuniform in density or composition as long as the concept of a local average density, $nr$, is valid.

For more general situations, for example an inhomogeneous dense fluid, a nanoscale fluid or solid, or a nonequilibrium system, the pressure **P** is a second-order tensor, depending on the direction of both the force and of the surface it acts on. In general, **P** has nine components *P*_{αβ}, where *P*_{αβ} is the force per unit area in the *β*-direction acting on a surface element normal to the *α*-direction. These components depend on position, **r**, and for nonequilibrium systems they will also depend on time, *t*. Off-diagonal components are the shear pressures (shear stress) and the diagonal components are the direct pressures. In some specific types of systems, the number of nonvanishing components of **P** may be less than 9. For an inhomogeneous fluid that is at equilibrium and not under strain, for example, the off-diagonal components vanish and there are only three nonvanishing components. Moreover, the condition of mechanical (hydrostatic) equilibrium (the average rate of change of linear momentum vanishes) often provides relations between the remaining nonvanishing components.

A difficulty in many applications is that the local (microscopic) pressure tensor at some point **r** is not uniquely defined in nonequilibrium systems or in equilibrium ones that are inhomogeneous. Although the kinetic contribution, **P**^{K}, is well defined, as noted above, the configurational part, **P**^{C} is not because the intermolecular forces themselves are nonlocal. Thus, while the force between molecules *i* and *j*, located at positions **r**_{i} and **r**_{j}, is well defined in general, there is no well-defined way to assign a contribution to the force acting on a surface element at some position **r**. This arbitrariness in the force acting across the surface element d*S* seems to have first been stated explicitly by Irving and Kirkwood (IK) in 1950.^{1} There, they stated the matter succinctly as follows:

“

…all definitions (of the configurational pressure tensor) must have this in common – that the stress between a pair of molecules be concentrated near the line of centers. When averaging over a domain large compared with the range of intermolecular force, these differences are washed out, and the ambiguity remaining in the macroscopic stress tensor is of negligible order.”

This point was discussed in a footnote to an appendix to Irving and Kirkwood’s paper and so was not noticed widely at the time. It was of little consequence to these authors, who were primarily interested in transport processes at the macroscale. However, it is important for nanoscale systems, such as small nanoparticles, drops and fluids confined within nanoporous materials or living cells, and we expand on this later in this Perspective. The nonuniqueness of the local pressure tensor, and its consequences for inhomogeneous fluids and the properties of gas–liquid interfaces, was investigated in the 1982 paper of Schofield and Henderson.^{2}

While in this Perspective we shall mainly focus our discussion on the local pressure tensor, in fields where the primary interest is in nonequilibrium phenomena and solid mechanics, the stress tensor $\sigma r,t$ is usually used in place of the pressure tensor.^{1} These two tensors differ only in sign,^{1,3}

The negative sign ensures that the stress tensor definition is consistent with Newton’s law of viscous flow and that the viscosity is positive. The terms “pressure tensor” and “stress tensor” are used interchangeably in this Perspective and in much of the literature, the sign change being understood.

The molecular level local pressure/stress tensor has been the key to the depiction of the mechanical and thermodynamic picture of many important phenomena. Its important role is evidenced by a rapid growth of publications mentioning it (Fig. 1). In biophysics, the local stress tensor has been applied to understand the mechanical properties of lipid bilayer membranes [see Fig. 2(a)].^{4–7} The structure and mechanics of the lipid membrane play a critical role in the function of proteins involved in processes of transport, signaling, and mechano-transduction. The local stress tensor also enables the quantification of the mechanical state of proteins in glassy matrices.^{8} Such information is pivotal to a sophisticated design and control of the lyophilization (freeze-drying) process for long-term storage and stabilization of labile biomolecules in the food and pharmaceutical industries. In material science, the stress tensor has been related to the structural deformation of materials upon adsorption, and such a connection is useful for materials characterization.^{9,10} For gas–liquid^{11–13} or liquid–solid^{14,15} nucleation, the pressure tensor profile provides a mechanical picture of the nucleus interfacing with the surrounding environment; such a profile is useful for calculating the Tolman length for interfacial free energy^{14} and for understanding the distinct structure of the nucleus^{15} [see Fig. 2(b)]. The pressure tensor profile across the interfacial region is also essential in a virial (or mechanical) route to the surface tension (see examples for planar,^{3,16–19} spherical,^{3,11,12,20} and cylindrical^{21} interfaces). For confined systems, the knowledge of the microscopic pressure tensor paves the way for understanding phase transitions in nanopores^{22–24} and for developing sophisticated equations of state for confined fluids.^{25,26} Recently, the microscopic pressure tensor has provided a mechanistic understanding of high-pressure phenomena in confinement or near strongly wetting surfaces for advanced materials synthesis and enhanced chemical processing.^{27} The high-pressure phenomena include enhanced chemical reactions in pores that normally require a high pressure in the bulk^{28,29} and the formation and stabilization of high-pressure phases in nanopores.^{30–32} For simple non-reacting adsorption systems, high (tangential) pressures that are about three to four orders of magnitude larger than the bulk pressure were found in the adsorbed layers on carbon surfaces [see Fig. 2(c)].^{33–35} This compression effect is caused by strong attractive forces exerted on the adsorbate molecules by the surface, which leads to a tightly packed adsorbed layer near the surface, and strong repulsions between adsorbate molecules.^{36–39} For more complex systems, the mechanism behind the induced high pressure in confinement is under active investigation.

The stress tensor also underpins fluid dynamics. It is the driving force in continuum equations,^{40} is necessary for the understanding of nonequilibrium molecular dynamics (NEMD),^{41,42} and is essential for calculating the viscosity in Newtonian fluids using both equilibrium^{43,44} and nonequilibrium methods.^{45} Pressure/stress can provide insight into tribology,^{46} reveal the slip length,^{47,48} and capture many aspects of multiphase flows, including the moving contact line for the wetting behavior [see Fig. 2(d)]^{49,50} and dynamics of active liquid interfaces.^{51} These molecular details can be included in continuum-based engineering simulation, such as computational fluid dynamics (CFD), through coupled simulation^{52,53} where stress coupling is useful in both fluid simulation^{54} and for solid mechanics.^{55} Recently, machine learning models have focused on stresses as the central property that should be predicted for predicting atomic stress at grain boundaries^{56} or in coupling to continuum models.^{57}

In the applications cited above, the range of length and time of interest can be very different and we comment on the effect of these differences in later sections of this Perspective. We treat equilibrium and nonequilibrium systems separately in what follows. The Perspective is organized as follows: In Sec. II, we introduce the fundamental pointwise forms of the pressure tensor. In Sec. III, we describe the local pressure tensor for inhomogeneous systems that are at thermodynamic equilibrium. We introduce a thermodynamic route to the local pressure tensor, alternative to the common mechanical route. In Sec. IV, we consider the definition of the pressure tensor for nonequilibrium systems, where **P** depends on time as well as position in space. In Sec. V, we point out several challenges associated with the microscopic pressure tensor along with our perspectives on future developments. In Sec. V A, the historical controversies over the microscopic pressure tensor are discussed in detail. These include the nonuniqueness of the local pressure tensor due to the arbitrary contour integral, the possibility of defining a unique coarse-grained (CG) pressure through integration of the local pressure over some spatial domain, and the questions over the existence of the kinetic term in the stress tensor. In this review, we focus mainly on the formalisms of the local pressure/stress tensor for discrete particles that interact with short-range, isotropic, pairwise potentials. Extensions to complex molecular and material systems interacting with many-body and long-range potentials are discussed in Sec. V B. Other practical aspects that are also examined include current availability of software and computational tools (Sec. V C) and the possible experimental routes to the validation of the microscopic pressure tensor (Sec. V D). Finally, concluding remarks are provided.

## II. FUNDAMENTAL EQUATIONS FOR THE POINTWISE PRESSURE TENSOR

The local pressure tensor for a particle (point-mass) system with both spatial and temporal dependence can be written as^{1,2}

where $PKr,t$ and $PCr,t$ are kinetic and configurational contributions, respectively, a tensor field defined at any given point in space **r** and time *t*. The kinetic contribution describes the flux of momentum,

where the angular bracket ⟨⋯⟩ denotes ensemble average, *N* is the total number of particles in the system, **p**_{i} is the momentum of the particle *i* and **p**_{i}**p**_{i} is an outer product, *m*_{i} is the mass of particle *i*, and $\delta r=\delta x\delta y\delta z$ is the delta function for a vector position in a Cartesian coordinate system. Here, the momentum may include a streaming component, $pi/mi=r\u0307i\u2212u$, due to the streaming velocity ** u** at point

**r**being nonzero, where $r\u0307i$ is the derivative with respect to time of the particle position. The mechanical definition of the kinetic pressure tensor in Eq. (3) corresponds to the ideal gas contribution in an equilibrium system (see Sec. III A) but must be defined in terms of streaming velocity in a nonequilibrium system (see Sec. IV C).

For notational simplicity, we assume pairwise interactions here (additional terms due to many-body interactions are considered later in Sec. V B). The configurational pressure **P**^{C} can be obtained from the tensor product of pair forces between particles, **F**_{ij}, and the line of interactions **r**_{ij} (**r**_{ij} = **r**_{j} − **r**_{i}),

where the pre-factor 1/2 accounts for double counting and the notation $\u2211i,jNFij=\u2211i=1N\u2211j\u2260iNFij$ has been introduced as shorthand for the double summation over all pairs of interactions. The *O*_{ij} term is known as the Irving–Kirkwood (IK) operator,^{1,61}

which is obtained as the Taylor expansion in space of the difference between two Dirac delta functions for molecule *i* and *j*,

If the expansion in Eq. (5) is simply truncated at the zeroth-order term, *O*_{ij} = 1, we reach the so-called IK1 expression for the configurational part of the pressure tensor,

This is the *virial* pressure applied locally at a point in space.^{62} For bulk homogeneous fluids where density is uniform throughout the system, the IK1 approximation (*O*_{ij} = 1) is exact.^{61} For inhomogeneous fluids, such as those confined in nanopores or near interfaces, the IK1 approximation violates the mechanical equilibrium condition^{63,64} and will lead to erroneous results as interactions with the surrounding fluids are not included.^{65} An exact form for the configurational pressure tensor can be reached by rewriting the IK operator using the fundamental theorem of contour integration,^{2,66}

where *C*_{ij} denotes an arbitrary contour from particle *i* to particle *j* and ** ℓ** is the corresponding contour vector. Substituting Eq. (8) into Eq. (4), we arrive at the contour form of the configurational pressure tensor,

^{2}

This equation is exact and allows the interaction between molecules to be included at an arbitrary location unrelated to the molecules’ positions.

While the kinetic pressure tensor [Eq. (3)] is well defined, the configurational part is nonunique due to the arbitrary contour involved in the calculation [Eq. (9)].^{2} In practice, a particular contour definition needs to be chosen to arrive at an operational form of the local pressure tensor. By introducing a concept of surface element d*S*, Harasima^{67} elegantly depicted two ways to assign a force contribution across such a surface of atomic dimension, which corresponds to two distinct definitions of the contour path. The first contour definition is a straight line of interactions between the molecules, known in the literature as the IK definition. It says that a pairwise force contributes to the pressure tensor at a surface element d*S* if the joining straight line (contour) between two molecules passes through d*S*. This definition is consistent with Newton’s assumption of impressed force between two points.^{68,69} The second definition was first implicitly adopted by Kirkwood and Buff^{70} and later referred to as the Harasima definition. It says a pair force contributes to the pressure tensor at the surface d*S* if one of the molecules lies in the cylinder whose base is d*S* (the axial direction of the cylinder is either parallel or perpendicular to the planar surface) and the other molecule is located on the other side of the plane of d*S*. Figures 3(a) and 3(b) illustrate the IK and Harasima definitions of the contour for a planar interface. Assuming the IK definition, the contour vector in Eq. (9) is simply ** ℓ** =

**r**

_{i}+

*λ*

**r**

_{ij}with 0 ≤

*λ*≤ 1

^{71}and Eq. (9) becomes

Alternative but equally valid contour definitions are possible. Figures 3(c) and 3(d) show possible variations of the IK and Harasima contour definitions. In general, the IK contour definition is arguably the most convenient and natural choice. It can be readily implemented for a general three-dimensional (3D) pressure field^{72} and in systems having an arbitrary geometry.^{8} Compared to the other widely adopted contour choice of Harasima, the IK definition has been shown to be physically consistent in different coordinate systems (spherical^{73} and cylindrical^{74} coordinates). If long-range Coulombic interactions are present in the system, however, the Harasima contour is preferred due to its compatibility and optimal computational efficiency with the Ewald summation method.^{6,74} No contour choice is more correct than the other in general. We will discuss this nonuniqueness problem further in detail in Sec. V A. For a bulk homogeneous system, the local pressure tensor is invariant to the choice of the contour definition. A spatial average of Eq. (9) over the entire system simplifies to the virial pressure tensor.^{75}

## III. MICROSCOPIC PRESSURE TENSOR IN EQUILIBRIUM SYSTEMS

For systems that are at thermodynamic equilibrium, the temporal dependence of the microscopic pressure tensor can be averaged out in the corresponding ensemble, and the spatial dependence of the pressure tensor is the main interest. Here we assume there is no shearing, as is the case for the equilibrium system, so that the off-diagonal elements in the pressure tensor simply vanish. In this section, we first discuss how the local pressure tensor is formulated in different geometries (coordinate systems) that are of practical importance for inhomogeneous systems at equilibrium. We then introduce a thermodynamic route to the microscopic pressure tensor. The equivalence between the thermodynamic route and the conventional mechanical (or virial) route is clarified.

### A. Local pressure tensor in different geometries

#### 1. Planar geometry

Systems that have planar interfaces are common and have been extensively investigated. Examples include fluids confined in a slit-shaped pore,^{35} two phases (e.g., gas and liquid) separated by a planar interface,^{16} and planar self-assembled layers (e.g., planar lipid bilayers).^{6,7} The Cartesian coordinate system is a convenient choice for planar geometries [Fig. 4(a)]. Here we take the *z*-axis to be the direction normal to the planar surface and assume homogeneity in the *xy*-plane, so that the local pressure tensor is a function of the *z*-position only. In the absence of external fields, the condition of mechanical (hydrostatic) equilibrium must be satisfied:^{2}

which has two implications:^{76} (1) the tangential pressure *P*_{T} has no local gradient in directions parallel to the *xy*-plane that induce flow anywhere; (2) the normal pressure *P*_{N} should be a constant throughout the system,

That is for a two-phase system in equilibrium (phases *α* and *β*) separated by a planar interface, the normal pressure in the bulk phase *α* is equal throughout the interfacial region and into the bulk phase *β* (Fig. 5). For the slit-pore system, if we consider the material as a part of the system rather than an external field (no gravity is considered here), the normal pressure is constant across the entire pore.^{35} It is worth noting that, due to the incommensurate packing of adsorbed layers in a slit pore, the normal pressure oscillates as the pore size increases.^{33} The normal pressure converges to the pressure of the bulk phase that is in equilibrium with the adsorbed phase, in the limit of an open surface (i.e., pore size is infinitely large).

The normal pressure is independent of the contour definition in Eq. (9).^{35} By taking the normal (*zz*) component out of the pressure tensor in Eq. (9) and integrating (averaging) over the *x*- and *y*-directions (i.e., *xy*-surface), the normal pressure is given by

where $nz$ is the local number density at *z*-position. The surface area is *S*_{z} = *L*_{x}*L*_{y}. If the sampling of the pressure tensor is carried out over the entire *xy*-plane, *L*_{x} and *L*_{y} will be the simulation box size in the *x*- and *y*-directions, respectively; *L*_{x} and *L*_{y} can also represent the size in the *x*- and *y*-directions for a specified region in the simulation box, respectively, if the pressure tensor is only sampled over that space.^{35} The scalar *ij*-pair force is $Fij=\u2212durij/drij$, where $urij$ is the pair potential between particles *i* and *j* separated by a scalar distance *r*_{ij}. *z*_{ij} = *z*_{j} − *z*_{i} is the *z*-component of vector **r**_{ij}, and $Hx$ is the Heaviside step function [$x>0,Hx=1;x<0,Hx=0;H0=1/2$]. We note that although the normal pressure is often written as a function of *z* as in Eq. (13) for use in molecular simulations, it is essentially a constant according to Eq. (12). The first term on the right-hand side of Eq. (13) is the kinetic (ideal gas) contribution to the pressure in equilibrium systems. It can be related to the mechanical definition in Eq. (3) at the limit of thermal equilibrium, by using the equipartition theorem,^{77}

where *p*_{α} is the molecule’s momentum in the *α*-direction with *α* = *x*, *y*, *z*^{78} and we use $nz=\u2211i=1N\delta z\u2212zi$. We note that Eq. (14) is valid in any classical system at thermal equilibrium, where the momenta and the coordinates are uncorrelated. For systems where quantum effects are significant, the equipartition theorem does not hold and a quantum version of the theorem should be applied.^{77} The second term in Eq. (13) is the configurational contribution arising from intermolecular interactions.

While the pressure normal to the *xy*-plane is well defined, the tangential pressure parallel to the plane is not uniquely defined and depends on the intermolecular interaction contour definition. Because the system is homogeneous in the *xy*-plane, the local tangential pressure *P*_{T} can be obtained by averaging over *P*_{xx} and *P*_{yy}. Unlike the normal pressure, which is independent of *z*, the tangential pressure does depend on *z*. The local tangential pressure based on the IK contour definition [Eq. (10)] is given by^{16,79}

Other equally valid contour definitions are possible. For example, assuming the Harasima (H) definition leads to^{16,67}

where the Dirac delta function can be approximated as

where $\Lambda zzi$ is the so-called boxcar or top-hat function of *z*_{i} for the interval from *z* − Δ*z*/2 to *z* + Δ*z*/2, which checks if *z*_{i} is less than *z* + Δ*z*/2 and greater than *z* − Δ*z*/2. In practice, we choose Δ*z* ∼ 0.001*σ* to numerically approximate the delta function where *σ* is the Lennard-Jones (LJ) diameter of the particle. We note that using Λ_{z} from Eq. (17) without taking the limit, Eq. (16) is equivalent to the 1D volume average (VA) form of the local pressure tensor. More discussions on the VA formalism will be presented in Sec. IV D.

For gas–liquid interfaces, molecular simulation results show that the IK and Harasima definitions yield tangential pressures that differ by less than 10% (Fig. 5).^{16,80} For liquid–solid interfaces, the difference between these two contour definitions can be considerably larger.^{34,35,63} The difference vanishes when the local pressure tensor is integrated over the entire system (or at least over the inhomogeneous interfacial region), which is the case for the surface tension calculations.^{16} The hydrostatic (averaged) pressure of the system can be calculated as the average of the trace of the pressure tensor, $Pz=Pxxz+Pyyz+Pzzz/3$. Spatially averaging $Pz$ over the entire system gives the bulk pressure, the configurational part of which is consistent with the virial theorem of Clausius.^{81,82} This bulk pressure is unique and independent of the contour definition.

#### 2. Spherical geometry

Examples of systems having a spherical geometry include gas/liquid/solid nuclei in the nucleation and crystallization process,^{11–13,15,83} fluids confined in spherical pores,^{84} and many spherical or quasi-spherical nanoparticles such as core-filled spherical nucleic acid.^{85} Due to the symmetry of the system, it is convenient to calculate the local pressure tensor in spherical coordinates, $R,\theta ,\varphi $ [Fig. 4(b)]. The local spherical pressure tensor can be written as

where *P*_{RR} = *P*_{N} is the normal pressure in the radial direction; *P*_{θθ} and *P*_{ϕϕ} are the equivalent tangential components (*P*_{T}) due to symmetry in the polar and azimuthal directions, respectively; and $R\u0302$, $\theta \u0302$, and $\varphi \u0302$ are unit vectors that are orthogonal to each other. Mechanical equilibrium [Eq. (11)] dictates the following:

In spherical coordinates, the Harasima-like definition of the pressure tensor leads to unphysical results in the bulk system due to the singularity near the origin.^{73,86} Therefore, the IK contour definition is commonly adopted in the literature. The normal (radial) pressure based on the IK contour definition is given by Ref. 87 (detailed derivations for all components are available in the supplementary material of Ref. 15),

where

and *λ*_{k} are the roots of a quadratic equation, $rij2\lambda 2+2\lambda ri\u22c5rij+ri2\u2212R2=0$. The polar pressure is given by^{15}

where

and $dxy=xi+\lambda kxij2+yi+\lambda kyij2$. The azimuthal component is given by^{15}

where

In practice, to enhance the statistics, the tangential pressure is calculated as the average of the polar and azimuthal components, $PT=P\theta \theta +P\varphi \varphi /2$.

#### 3. Cylindrical geometry

Compared to the local pressure tensor for planar and spherical geometries, the theoretical development in a cylindrical geometry is generally overlooked. A complete derivation for the cylindrical pressure tensor was made available very recently.^{21,74} The cylindrical pressure tensor is useful for understanding the behavior of systems having cylindrical interfaces. Examples include self-assembled micelles of a cylindrical shape,^{88,89} a solid nucleus having a cylindrical shape,^{90,91} and molecules confined in cylindrical pores.^{74} Most of the synthetized porous materials with a well-defined geometry have cylindrical or quasi-cylindrical pores, such as carbon nanotubes and porous silica materials. The local cylindrical pressure tensor can be written as [Fig. 4(c)]

where *P*_{RR} = *P*_{N} is the normal pressure in the radial direction ($R\u0302$); *P*_{ϕϕ} is the tangential pressure in the azimuthal direction ($\varphi \u0302$); and *P*_{zz} is the tangential pressure in the axial direction ($z\u0302$) with *P*_{zz} ≠ *P*_{ϕϕ}. Mechanical equilibrium [Eq. (11)] dictates that

Similar to the situation in spherical coordinates, the Harasima-like contour definition leads to a cylindrical pressure tensor with an unrealistic radial dependence near the origin in a bulk system.^{74} In general, we suggest that any construction of the contour should not depend on polar coordinates (i.e., radius and polar angles).^{74} While a valid alternative to the Harasima-like definition in cylindrical coordinates is possible,^{74} the IK contour definition naturally satisfies this polar-coordinate-independence condition, being a widely adopted choice of contour. The normal component of the cylindrical pressure tensor based on the IK contour definition is given by^{21,74}

where *L* is the height of the cylinder. The unit radial vector is

and *λ*_{k} are the roots of the equation $R2=xi+\lambda xij2+yi+\lambda yij2$. The azimuthal component is written as^{74}

where

Moreover, the axial component of the cylindrical pressure tensor is given by

We note that Eqs. (30) and (32) are equivalent to the corresponding ones in Ref. 21 despite the difference in their appearance. Unlike equations in Ref. 21 that should be averaged over a number of *ϕ* and *z* values in the simulation for the azimuthal and axial pressure, respectively, here we have already averaged over all possible *ϕ* and *z* analytically through integration.^{74}

#### 4. Arbitrary geometry

While many systems of practical interest can be approximated with the three aforementioned simple geometries (planar, spherical, and cylindrical), there are still cases where systems have complex shapes without well-defined symmetries, and special treatments are needed. Typical methods for handling an arbitrary geometry involve the discretization of the system into local micro-volumes of molecular dimensions, and the local pressure tensor is evaluated either on the surface of the local volume^{66} [Fig. 4(d)] or as a spatial average over such local space.^{8,72} A more general form of this kind is discussed in Sec. IV D for nonequilibrium systems.

### B. Thermodynamic route to the pressure tensor and its equivalence to the mechanical route in thermodynamic equilibrium

So far, we have only talked about the mechanical route to the local pressure tensor, which is derived from the concept of “the force acting across a surface element d*S*.” Equivalently, the pressure tensor can also be derived from a thermodynamic definition. It is instructive to start with the definition of the bulk (scalar) pressure in a canonical (NVT) ensemble (generalization to other ensembles is straightforward),

where *A* = −*k*_{B}*T* ln *Q* is the Helmholtz free energy of the system with *N* particles and a volume *V* at temperature *T*. The canonical partition function *Q* is defined as^{82}

where **r**^{N} ≡ **r**_{1}, **r**_{2}, …, **r**_{N} represent the positions of all particles in the system, Λ is the de Broglie wavelength, $U$ is the total configurational energy of the system, and *β* = 1/*k*_{B}*T*. Following the pioneering work of Eppenga and Frenkel^{92} on hard-core particles and Harismiadis *et al.*^{93} on systems with continuous potentials, the bulk pressure can be computed by considering a volume perturbation (VP) from *V* to *V*′ = *V* + Δ*V* (particle coordinates are scaled accordingly), with Δ*V* > 0 being an infinitesimal, *isotropic* change of the volume, and Eq. (33) becomes^{93}

where $\Delta U=UV+\Delta V\u2212UV$ is the energy associated with the increase in volume and the angular bracket with subscript *V* denote a configurational average in the canonical ensemble over the unperturbed system of volume *V*. We note that Eq. (35) includes both kinetic and configurational contributions to the pressure. Equation (35) also assumes the volume change is positive; in practice, a central finite-difference approximation with both positive and negative volume changes is recommended for better statistics.^{94} Equation (35) is the VP method or thermodynamic route to the pressure. A similar perturbation scheme for the surface (test-area method) has also been developed to compute the surface tension.^{95–97} This test-area method is advantageous over the conventional pressure tensor route to the surface tension, especially for small drops, because it can capture the entropic contributions due to the fluctuations in the energy of deformation.^{98}

Similarly, Eq. (33) can be rewritten to represent the diagonal Cartesian components of the pressure tensor *P*_{αα} (*α* = *x*, *y*, *z*),^{94}

In this case, instead of performing an *isotropic* change of the volume as for the bulk pressure, only the simulation box dimension *L*_{α} in the *α*-direction is perturbed to *L*_{α} + Δ*L*_{α} while all other dimensions $L\beta \beta \u2260\alpha $ are kept fixed. Now, we apply a common planar symmetry, i.e., the *z*-axis is normal to the planar surface and the surface lies in the *xy*-plane, and we assume Eq. (36) is locally valid (i.e., *Q*, *A*, and *V* can be localized to a thin slab at *z*). By taking the partial derivative in Eq. (36) exactly, we can derive a local form of the pressure tensor,^{34,94,99}

where $Uz$ is the total configurational energy in an infinitesimally thin slab at a *z*-position and $Vz=LxLy\Delta z$ is the corresponding volume of the thin slab of width Δ*z*. The first term on the right of Eq. (37) is the kinetic contribution and the second term is the configurational contribution due to intermolecular interactions. The second line in Eq. (37) is a VP approximation similar to Eq. (35), but here the infinitesimal virtual expansion is performed only in the *α*-direction, and $\Delta Vz=S\alpha \Delta L\alpha \Delta z/Lz=S\alpha \zeta L\alpha \Delta z/Lz$, where *S*_{α} is the surface area normal to the *α*-direction and *ζ* is a positive, infinitesimal number. In practice, the coordinates of all particles in the system (including those not in the thin slab at *z*) should be scaled by *ζ*. $\Delta Uz$ is the change of the total potential energy in the thin slab at *z* due to the corresponding virtual volume expansion. For confined fluids, for example, the total potential energy change should include contributions from both fluid–fluid and fluid–solid interactions.

Different criteria of assigning potential energy in slabs will lead to different local pressure tensor values,^{34,100} again reflecting the nonuniqueness of the configurational contribution to the local pressure tensor. The ambiguity associated with the localization of the potential energy was also realized by Irving and Kirkwood.^{1} Figure 6 illustrates two ways of assigning a pair interaction potential in space, labeled methods A and B. Method A assigns half of the pair potential energy, $urij/2$, into the thin slab that contains particle *i*, and the other half into the slab that contains particle *j*,^{34,101}

Using Eqs. (17) and (38), it is straightforward to show that the first line of Eq. (37) (which is exact) for pairwise interactions will lead to the Harasima definition of the local tangential pressure in Eq. (16) by taking Δ*z* → 0 (see Ref. 74 for similar derivations in cylindrical coordinates). Method B equally distributes the pair potential $\mu rij$ into slabs between two interacting particles *i* and *j*,^{34}

This energy attribution rule corresponds to the IK contour definition for the local pressure tensor [Eq. (15)]. Similar to the integral contour whose choice is restricted in polar coordinates (for example, Harasima contour leads to invalid results in polar coordinates^{73,74}), the ways to assign the local energy should also be regulated under certain conditions. Further studies are required to elucidate such restrictions. Finally, we note that the nonuniqueness in assigning the local energy may also explain the nonunique nature of the local heat flux.^{102,103}

Extension of the thermodynamic route to curved interfaces, such as spherical^{76} and cylindrical shapes,^{21,76} is possible. However, caution should be exercised due to the conjugate nature of polar components in the pressure tensor. For example, in spherical coordinates, perturbing the *R*-coordinates will *not* result in radial pressure because such rescaling also leads to perturbations in *ϕ*- and *θ*-directions. We will illustrate this point as follows: Let us assume we perform a small perturbation in the radial direction of a sphere having a radius of *R*_{0}, $R0\u2032=1+\zeta R0$, where *ζ* is a positive, infinitesimal number. The volume change of the entire system is

where in the last step we omit higher-order terms involving *ζ*^{2} and *ζ*^{3}. By perturbing the radius, the reversible work (free energy change) is done in all three polar directions,

The integral in Eq. (41) is carried out over the entire sphere because the reversible work is done on the entire system (i.e., all *R*-positions are rescaled with *ζ*). Using Eqs. (40) and (41), and the thermodynamic definition of the pressure in Eq. (33), we get

where *P*_{RR} = *P*_{N} is the normal pressure and *P*_{θθ} = *P*_{ϕϕ} = *P*_{T} is the tangential pressure due to symmetry. Equation (42) clearly shows that perturbing the radial direction leads to a hydrostatic (bulk) pressure of the system defined as an average of the trace of the spherical pressure tensor. To decouple the normal and tangential pressure, we take advantage of the mechanical equilibrium condition in Eq. (19). Interested readers should refer to Ref. 76 for derivations of the local version of Eq. (42).

In short, the thermodynamic route (VP method) gives an equivalent form of pressure tensor to the one obtained from the mechanical route for equilibrium fluid systems (with off-diagonal elements being zero). Molecular simulation results indeed confirm this equivalence.^{34,76,94,100} While both of them are useful for the calculation of the pressure tensor profile in inhomogeneous systems, they each have strengths and limitations. While the mechanical route is useful in studying the solid phase, the thermodynamic route is invalid in such a case. For a solid, the off-diagonal elements may not be zero due to internal strain. That means if the solid is perturbated in one direction, as the thermodynamic route suggests, the results will be a coupling of the shear modes and the direct pressure. Nevertheless, the thermodynamic route is arguably more convenient for systems interacting with complex intermolecular potentials (e.g., many-body interactions) where an explicit evaluation of the forces might be computationally challenging. Moving away from thermodynamic equilibrium, the mechanical route is the preferred choice, as discussed in Sec. IV.

## IV. MICROSCOPIC PRESSURE TENSOR IN NONEQUILIBRIUM SYSTEMS

In this section, we extend the discussions of Sec. III to consider the definition of pressure with temporal evolution, convection, and flow in a moving reference frame. Moving away from equilibrium, we expect variation of the pressure tensor to drive flows and for the pressure/stress tensor to depend on time as well as space (inhomogeneous). A simple example of this is a dynamic, or hydrodynamic, equilibrium in steady-state Couette flow, where the boundary conditions (or forcing/molecular walls) drive a flow giving a linear velocity profile and a constant shear stress. Many different computational techniques have been developed to study this system in nonequilibrium molecular dynamics (NEMD), including SLLOD,^{104,105} Lees–Edwards boundary conditions,^{106} and the application of shearing through tethered or fixed walls.^{107,108} More generally, the onset of turbulence introduces convective terms that are nonzero even in an average sense. Finally, the most general case of nonequilibrium is an unsteady flow, where the gradient of convective transport and pressure together are equal to the time evolution of momentum in the system,

where *ρ* is the mass density and ** u** is the streaming velocity. This time-evolving flow can occur in bubble growth, moving contact lines, onset of instabilities, and many other areas of fluid dynamics.

We start by discussing the importance of temporal evolution on the definition of the stress tensor for a statistical mechanical ensemble in Sec. IV A before simplifying to the case of a single trajectory evolving in time in Sec. IV B. We discuss the problem of defining kinetic pressure as we move away from equilibrium in Sec. IV C. For the time evolution and convection, we need a clear mathematical framework for spatial localization of pressure valid away from equilibrium, and these are discussed in Sec. IV D. In Sec. IV E, we present recent studies and consider the most general case where the framework itself can move in time, which is useful in multiphase fluid flow with deforming interfaces. We then outline a more pragmatic way to get this time evolution using mapping in Sec. IV F. Finally, we discuss the statistical uncertainty involved in the pressure tensor calculations in Sec. IV G.

### A. Ensemble average and the time-evolving phase space

Let $BrN,pN$ be some function of the 6*N* phase space variables. The ensemble average ⟨*B*⟩ is then the 6*N*-dimensional integration over all position vectors and over all momenta vectors for an ensemble having a probability density function $frN,pN;t$. Using the assumption that phase space is bounded,^{1} and noting that for momentum $B=\u2211i=1Nmir\u0307i\delta r\u2212rit$, the time evolution of a phase space averaged momentum is

where $r\u0307i$ is the first-order derivative of the particle position with respect to time *t*. The right-hand side is the pressure introduced in Eqs. (3) and (4).

There are at least two purposes served by the ensemble average: the first is the practical reduction of noise and the second is to ensure the validity of the Dirac delta function, which is not practically meaningful outside an integral (here the 6*N* dimensional integral of the ensemble average). The work of Noll^{109} integrates the Dirac delta function over phase space and uses Noll’s lemma to address the *O*_{ij} operator of Eq. (5), giving a form with similarities to the line integral of Eq. (8), but in Noll’s phase space integrated notation. The form is included in Subsection 2 of the Appendix and interested readers should refer to the review paper by Admal and Tadmor.^{78}

### B. A single trajectory in time

If the ensemble average is dropped, equivalent forms of the pressure can still be obtained.^{41} This step is essential in order to derive general equations free from the requirement of ensemble averaging, *∂*⟨*B*⟩/*∂t* → *∂B*/*∂t*. Here, the form of $B=\u2211i=1Nmir\u0307i\delta r\u2212rit$ is as before. Applying the time derivative of momentum results, after some manipulation,^{41,42} in the following:

Here, Newton’s law, $Fi=mir\u0308i$, is applied before expressing the force in a pairwise manner and the delta functions are expanded as in the original work of Irving and Kirkwood.^{1} The final right-hand side is the pressure at any instant, as in Eq. (44), but without the ensemble average. This has the advantage that it is a purely mechanical form of the pressure. More importantly, it is the starting point for a more general treatment.

However, we have introduced the formal mathematical problem to this pointwise pressure form that a Dirac delta function now exists outside an integral, which makes it poorly defined. To solve this problem, approaches in the literature approximate the delta function as a mollified weighting function^{78,110–112} or approximate kernel.^{4,108,109} A weighting function can be chosen to have any desired mathematical property, such as compact support or normalization to unity. In this work, we do not approximate the delta function, instead we evaluate the formal integral of the delta function over a local volume in space.^{113,114,116} This control volume or “finite volume” form^{115} can be written in terms of surface tractions and ensures the exact conservation of momentum during the single-trajectory time evolution as shown in Sec. IV D. Our treatment is equivalent to a choice of a uniform weighting function (a 3D boxcar function)^{77} without any mollification. A weighting function may not satisfy momentum conservation as it smears the momentum average over its functional form. This is analogous to the finite element method where different shape functions give different momentum distributions^{117} and only the zeroth-order element (a finite volume) is conservative.^{115} In Sec. IV C, the importance and difficulty of identifying and subtracting the hydrodynamic or streaming velocity are discussed, particularly when we move away from the ensemble picture.

### C. Streaming velocity and the kinetic term

The dynamics of a fluid manifests itself through a velocity field ** u**(

**r**,

*t*) coupled with a scalar pressure

*P*(

**r**,

*t*) at every point in space. This fluid velocity can be thought of as the average coherent motion of a stream of molecules, which we can define in terms of an instantaneous form of the Irving–Kirkwood momentum and density,

Here, *ρ*(**r**, *t*) is the mass density of the fluid and ** u**(

**r**,

*t*) is the average velocity of the molecules, often known as the streaming velocity. The kinetic pressure can then be defined by introducing the peculiar velocity

**p**

_{i}/

*m*

_{i}, the particle motion not contributing to the net velocity field,

A clear problem with the instant trajectory is apparent here: For a single time step in a molecular simulation, there is no way to split fluctuation and streaming parts. Of a molecular kinetic motion, the contribution that becomes velocity and the contribution that is kinetic pressure can only be determined with both spatial or temporal averaging. We discuss the spatial averaging in Sec. IV D. The length of temporal averaging will adjust measurements of velocity and requires some pragmatism in choice. However, this does not mean that the use of a single trajectory is wrong, as will be discussed below. In fact, chaotic trajectory divergence means an ensemble approach is not always possible in situations away from an attractor state,^{78,118} so piecewise temporal averaging may be the only approach to get velocity in highly nonlinear systems.

Taking the kinetic term in Eq. (45), substituting in the peculiar velocity of Eq. (47), and applying the definition of momentum and density allow the total kinetic tensor to be written in terms of a kinetic pressure and a convection term,

In an equilibrium system, the second term on the right (convection) is equal to zero. In the case of hydrodynamic equilibrium, for example, in channel flow, convection is also negligible, and the velocity varies in a well-defined way. This is true even in molecular systems down to 5–10 atomic diameters,^{42} e.g., a linear function *u*(*y*) ∼ *y* in Couette flow^{119} or parabolic in pressure-driven Poiseuille flow *u*(*y*) ∼ *y*^{2}.^{120} When the expected velocity profile is known, such as in Couette or Poiseuille flow, it makes the definition of peculiar velocity, **p**_{i}/*m*_{i}, straightforward.^{121} This idea of constructing quantities with a known velocity inspired some of the earliest nonequilibrium molecular dynamics (NEMD) simulation methods, by applying a known forcing function, so that the resulting velocity is as expected.^{122} However, as molecular simulation increasingly pushes to more complex nonequilibrium cases, a well-defined velocity profile is no longer possible. Examples of nonequilibrium systems include complex flow patterns around obstacles,^{123–125} rolling or cells formed by thermal gradients,^{126,127} vortex formation,^{128} Taylor Couette flow,^{129,130} the Rayleigh–Taylor instability^{131,132} or shock wave instability.^{133,134} There is also a recent explosion in papers using molecular simulations in multiphase flows for films, droplets, bubbles, and contact lines. We defer a consideration of these to Sec. IV E, where a moving interface reference frame is presented specifically for these types of problems. In all nonlinear single-phase examples, the evolution of velocity is chaotic, which makes the concept of an ensemble averaging problematic, as different simulations would diverge given enough time.

Perhaps the most interesting and general example of splitting the streaming velocity from the kinetic contribution is for turbulent flow.^{135} This is shown in Fig. 7 for the smallest known cases of time-steady turbulent flow^{136,137} but simulated in a molecular dynamics simulation.^{138} This minimal-channel Couette flow exists at a Reynolds number of 400, which requires *N* ≈ 300 × 10^{6} molecules taking an average density of 0.3 in LJ units and temperature chosen to minimize viscosity. Convection is nonzero and velocity varies in time, so we use Reynolds decomposition, $u=u\u0304+u\u2032$, to define a long-time-average streaming velocity $u\u0304$. This average velocity allows us to define a turbulent fluctuating part, denoted here by the prime ** u**′. This

**′ is like the peculiar velocity in NEMD, a velocity that is not contributing to the mean flow, but the fluctuations are large-scale eddying motions instead of molecular fluctuations of**

*u***p**

_{i}. The typical cycle of fluctuations is shown in Fig. 7(a) as the iso-surfaces of squared velocity |

**′|**

*u*^{2}=

*u*′

^{2}+

*v*′

^{2}+

*w*′

^{2}where

*u*′,

*v*′, and

*w*′ are the three velocity components. As with the kinetic pressure tensor, knowing the average velocity allows the average fluctuations to be identified, which is called the Reynolds stress tensor,

This is a dimensionally and physically the same form as the kinetic pressure, the outer product of fluctuating velocity components, but on a larger scale. Instead of small molecular fluctuations, it is the average momentum carried by turbulent eddies that make up the Reynolds stress. Indeed, Osborne Reynolds apparently defined this stress, inspired by the subtraction of streaming velocity seen in his earlier work on kinetic theory.^{139,140} The mean flow is time stationary on a longer time scale, with the flow going through a regeneration cycle, where the streaks break down the energizing flow vortices before a regeneration occurs [Fig. 7(a)].This cycle repeats indefinitely and a long-time average can be collected. The average pressure and Reynolds shear stress in the top half of the symmetrical channel are shown in Fig. 7(b) where the total shear contribution is constant, i.e., $[\rho u\u2032v\u2032\u0304+PxyK+PxyC]/\tau 0=1$. Interestingly, the kinetic and configurational contributions to shear stress are similar in magnitude in Fig. 7(b). It is similar to the case of hydrodynamic stability in laminar Couette flow, where adding the average Reynolds shear stresses together with the shear pressure gives a constant, i.e., it satisfies $\u2207\u22c5P+\rho u\u2032u\u2032\u0304=0$, and as a result $du\u0304/dt=0$. This turbulent equilibrium is a property of the relatively simple channel flow and would not be true in general.

The off-diagonal components of the Reynolds stress tensor in Eq. (49), e.g., $\rho u\u2032v\u2032\u0304$ are like the shear stress $PxyK$ in that they include *x*-momentum carried in the *y*-direction. As a result, the kinetic pressure can be split into velocity existing on three length and time scales,

Here, we have included a time averaging, ⟨⋯⟩_{t}, over a shorter time scale *t*_{MD}. The overbar denotes an average over the length of the simulation, in this case *t*_{Total}. Equation (50) highlights a very interesting insight into the time averaging process; as the averaging period of ⟨⋯⟩_{t} increases relative to the overbar, the kinetic contribution will increasingly be assigned to molecular kinetic shear stress **P**^{K} and the Reynolds stress will decrease, with the limit *t*_{MD} → *t*_{Total} seeing all turbulent fluctuations counted as kinetic pressure. In Eq. (50), we see kinetic pressure as simply the first-order term in a series expansion of fluctuations at increasing length and time scales. This is a source of nonuniqueness in the kinetic pressure away from equilibrium, intimately linked to the turbulent cascade, where fluctuations ** u**′ form a continuous spectrum of scales. In the turbulence literature, understanding and modeling the range of turbulent scales are main focuses of research.

^{141,142}This section shows that the uncertainty in splitting molecular motion into kinetic pressure and streaming velocity is part of a larger picture in fluid dynamics. Here, flow is multi-scale over many orders of magnitude and care is needed to decide which scales must be modeled. In experimental aerodynamics, the pressure measured by a pitot tube, which physically slows the flow velocity to zero, or stagnation, is known as the total pressure. This is to distinguish it from static pressure due only to kinetic collisions measured in a static flow. A complete picture of the kinetic pressure and fluid flow will require a collaboration between the fluid dynamics and NEMD communities.

^{141,142}

### D. Localization of the pressure tensor in space

The field of NEMD closely mirrors fluid dynamics, where fluid properties are defined on a field, which is expressed as a discrete grid of values in CFD. For this reason, expressing the pressure on a tessellating grid of cells is a natural choice to explore fluid phenomena. This process starts by integrating the Irving–Kirkwood form of pressure over a volume *V* in space, including the kinetic pressure of Eq. (3) and the configurational pressure of Eq. (10) using the IK contour,

where *ϑ*_{i} and *ϑ*_{λ} are the integral of the Dirac delta functions over a finite volume, in the case of a cuboid centered at $r=x,y,z$ and having a dimensional Δ**r**. The function *ϑ*_{i} can be written as $\u03d1i=\Lambda xxi\Lambda yyi\Lambda zzi$, where $\Lambda \alpha \alpha i$ with *α* = *x*, *y*, *z* is the boxcar function introduced in Eq. (17) and *ϑ*_{i} = 1 when the particle *i* is inside a cuboid and *ϑ*_{i} = 0 otherwise. Assuming a single average value of pressure in a volume, $\u222bVPdV=PVA\Delta V$, Eq. (51) results in the so-called volume average (VA) pressure tensor,

where *ΔV* is the local volume, and the shorthand $lij=\u222b01\u03d1\lambda d\lambda $ gives the fraction of the interaction contour ** ℓ** inside the averaging volume. We note that the streaming velocity $ur,t$, defined in Eq. (46), is averaged over the same volume in the VA form [Fig. 8(a)], i.e., $u(r,t)=\u2211i=1Nmir\u0307i\u03d1i/\u2211i=1Nmi\u03d1i$. This VA pressure was originally proposed in a 1D form for shockwaves by Hardy in 1982.

^{110}It was then extended by Cormier

*et al.*

^{143}to a spherical volume and made more formal by Murdoch.

^{111,112}For the case of the IK contour, the fraction of line

*l*

_{ij}can be obtained exactly in a cube from plane–line intersections shown in Fig. 8(b) or in a sphere

^{87,144}or cylinder

^{21}from surface–line intersections. For more complicated local volume shapes,

*l*

_{ij}can be obtained by splitting the line into segments and binning each segment numerically if it is inside the volume. As an example, the Voronoi decomposition of Hatch and Debenedetti

^{8}is shown in Fig. 8(c). They developed a formalism that enables the calculation of the local stress tensor on an atom or an arbitrary group of atoms by averaging over the local volume of this group. The local volume for a certain group was obtained using the Voronoi decomposition method. In their formalism, the contour segment that is within the volume

*V*

_{g}of the targeting group of atoms contributes to the local stress of this group. For instance, in Fig. 8(c), the red segment of the line connecting particles 3 and 4 contributes to the local stress of the group composed of particles 2 and 5. In this way, it is identical to the VA pressure of Eq. (52) using length of interaction inside the volume, but, with the volumes chosen based on molecular structure. The presence of a molecule at the center of a volume changes the measured pressure to include insight into the material structure, with some similarity to the molecular centric radial distribution function.

^{144}

On the other hand, the pressure tensor can be expressed in a plane form. For example, by taking the three components of the VA pressure tensor in Eq. (52) acting on the surface that is normal to the *x*-direction, i.e., $Px=Pxx,Pxy,Pxz$, and evaluating the limit that the volume tends to zero in the *x*-direction, we arrive at^{145}

where Δ*S*_{x} = *Δy**Δz* is the surface element normal to the *x*-direction. In the limiting case, Eq. (17) is used to convert a boxcar function to a delta function. If we take limits of Δ*y* and Δ*z* to be the edges of the simulation box with periodic boundaries, we have $lim\Delta y\u2192Ly\Lambda y=lim\Delta z\u2192Lz\Lambda z=1$, Δ*S*_{x} = *S*_{x}, and Eq. (53) simplifies to the method of planes (MoP) pressure,^{61}

The integral along *λ* has been evaluated to give signum functions (for a contour, use $\u222exixj\delta x\u2212\u2113xf(\u2113x)d\u2113x=f(x)Hx\u2212xi\u2212Hx\u2212xj$), where $sgnx=1$ for *x* > 0 and $sgnx=\u22121$ for *x* < 0, and $sgnx=0$ otherwise. The interpretation of the differences in signum functions is that *x*_{i} and *x*_{j} have to be on opposite sides of a *yz*-plane at *x* for the expression to be nonzero. This is the condition that the interaction is crossing the surface, and the force contribution is included. In this way, the MoP is most clearly related to the force over area definition of mechanical stress. The early work of Tsai^{65} postulated this form of stress in a molecular system, but the work of Todd *et al.*^{61} derived it from statistical mechanics for the first time and provided a convenient form to use in molecular simulations. This was originally obtained through a Fourier transform of the original Irving–Kirkwood equations.^{1,61} Working in Fourier space has the effect of averaging in the lateral directions and so the pressure is for an infinite plane. Regardless of the contour between two molecules, the contour must cross the plane if these two molecules are located on either side of the plane. The Fourier transform assumes an infinite periodic domain in *y* − *z* and so avoids the need to choose an integral contour. However, we can prove that, due to mechanical equilibrium, the ensemble average of the MoP form in Eq. (54) is equivalent to Eq. (13), which is derived directly from the contour form of the pressure tensor (see Subsection 1 of the Appendix for derivation). It is worth noting that the kinetic term in Eq. (54) still has the Dirac delta function, which requires some treatment before it can be used in a molecular simulation. Previous work writes the delta function as the sum of its roots,^{146} which physically correspond to any crossing of the plane. However, a kinetic pressure form that is more consistent with the configurational term can be obtained by taking the integral of the kinetic pressure over a time interval from *t*_{1} to *t*_{2} and using a change of variables $dt=dxi/x\u0307i$ to write

The final line uses $H(x)=1/2(sgnx+1)$ to show the similarity of form to the common MoP configurational part and takes the average of the time integral on the left-hand side $\u222btt+\Delta tPMoPxKdt\u2248\Delta tPMoPxK$. Finally, we note that, to be consistent, the momentum flux in these surface definitions **p**_{i} is expressed relative to a streaming velocity measured over the same surface, so if $dSix=sgnx\u2212xi(t1)\u2212sgnx\u2212xi(t2)$ then $u(x,t)=\u2211i=1Nmir\u0307idSix/\u2211i=1NmidSix$.

The original MoP formulation^{61} only provides three pressure components as in Eq. (54), on an infinite plane. As the plane has a single normal component and is infinite in *y* and *z*, this returns only a single pressure vector per plane. Han and Lee^{147} used three mutually perpendicular planes converging at a point to obtain all nine components of the pressure tensor and also limited the planes to a local region of interest. For example, using the boxcar function Λ of Eq. (17), we can write this Local MoP (LP) for the kinetic and configurational pressure components on a local plane that is normal to the *y*-direction, i.e., *P*_{yx}, *P*_{yy}, and *P*_{yz} as

where Δ*t* = *t*_{2} − *t*_{1} and the function denoted as *dS*_{iy} checks if molecule *i* is crossing the plane [Fig. 9(a)]. Function *dS*_{ijy} checks the crossing of the plane for intermolecular interaction path between *i* and *j* [Fig. 9(b)]. Here, on the left-hand side of Eq. (56), we are using the shorthand notation for the average pressure integrated over a plane, $PLPy\Delta Sy\u2248\u222bSyP\u22c5dSy$. By localizing the pressure to a region of a plane, we are forced to choose a contour between molecules as different contours may or may not cross the corresponding plane subsection. Using the IK contour, the boxcar function Λ_{x}(*x*_{k}) in Eq. (56) identifies if a crossing *x*_{k} = *x*_{i} + *λ*_{k}*x*_{ij} is on the surface, with *λ*_{k} being the value of *λ* along the line of contour integration at the point of crossing of the plane/surface. It is interesting to note that the kinetic pressure in this form is also now dependent on the trajectory “contour” of the molecules as they evolve in time, so that different integration methods or time steps could result in different measured crossings and therefore a nonuniqueness in the kinetic pressure. In practice, molecules move by very small amounts in a time step, so any difference is unlikely to be apparent.

For now, we have derived the surface pressure forms of Eq. (54) (MoP) and Eq. (56) (LP) by taking the zero volume limit of the VA form of Eq. (52). However, the localized surface pressure is more rigorously derived by taking the derivative of the control volume pressure in Eq. (51),^{148} which gives the six local surfaces pressures bounding an enclosed region,

where $PLP\alpha +$ and $PLP\alpha \u2212$ (*α* = *x*, *y*, *z*) denote the LP pressure tensor on upper (+) and lower (−) *α*-surfaces, corresponding to front/back, right/left, and top/bottom pairs in the *x*-, *y*- and *z*-directions, respectively [see Fig. 4(a) for the coordinate system]. The pressure on all six surfaces of an enclosed volume is exactly equal to the momentum change in time (to machine precision) due to the conservative property of the finite volume form.^{115} This is denoted by the final line in Eq. (57), which states that the integral around the bounding surface is equal to the change inside that volume. It is also possible to use this technique for any volume, for example, the surface pressure for spherical volumes have been derived in the literature^{144} and we consider more general surfaces in Sec. IV E. To show this conservation, we introduce the shorthand *P*^{±} = *P*^{+} − *P*^{−} and define Advection to include the convective term and kinetic pressure, a forcing term including configurational pressure on a surface, and plot both against the change inside the volume, called "accumulation” in Fig. 10.

This process of surface flux (SF) derivation can also allow more general volumes, providing a form of pressure on any arbitrary surface (e.g., a rippling surface), as shown in Sec. IV E. Mechanically, each of the terms can be viewed as the tractions in a Cauchy tetrahedron, so using the three orthogonal planes on the top surfaces with a traction force vector on each, cf. Fig. 4(a), the three surface traction vector can be assembled to give $Pr,t=PLPx+,PLPy+,PLPz+T$, which is a nine-component pressure tensor. A similar tensor could be defined for the bottom set of surfaces defining the other tetrahedron that together makes up the cube volume. It is worth noting that these surface pressures can be derived directly as surface localization of phase space quantities^{148} that satisfy the requirements of statistical mechanics.^{78} However, the real strength in these instantaneous stresses on a bounding surface is that we can get the stress at any instant that is directly responsible for a corresponding momentum change inside the volume.

### E. A moving reference frame

For interactions crossing an arbitrary surface that is itself changing in time, we can define a volume that follows a feature of the simulation. For example, we consider a two-dimensional example of a liquid–vapor interface (projected onto the *yz*-plane) defined by a generalized function *ξ*(*y*, *z*, *t*). Using this moving interface function, the volume average equation of Eq. (51) can be repurposed with the boxcar function redefined as $\Lambda x=Hxi\u2212x+\Delta x/2+\xi (y,z,t)\u2212Hxi\u2212x\u2212\Delta x/2+\xi (y,z,t)$. Taking the derivative of Eq. (51) gives the surface pressure, as before, but it is now on a curved and moving surface and we refer to it as the surface flux (SF) pressure $PSF$. Working through the mathematics, the kinetic and configurational SF pressure tensor can be expressed as follows:^{66}

Equation (58) is identical to the LP pressure of Eq. (56) but with all the extra terms accounting for interface curvature and movement. The surface location that the molecular trajectory (or intermolecular interaction) crosses is denoted by $\xi i+$ (or $\xi \lambda +$). The $dSix+anddS\lambda x+$ are functions which are nonzero only if these crossing locations are on the particular patch of surface [a generalization of the signum functions of Eq. (56)]. The underbraces highlight the physical meaning of the various terms. The surface evolution term accounts for the movement of the interface in time, and the kinetic and configurational curvature terms ensure that the *y* and *z* components are included on the *x*-surface. The convective term is included on the left-hand side as surface movement makes it more difficult to determine what is convection and what is kinetic pressure. By introducing the definition for the surface normal vector $n\u0303x=\u2207\alpha \xi \u2212x\alpha /\u2207\alpha \xi \u2212x\alpha $ and evaluating the integrals in Eq. (58), we obtain a form of pressure purely in terms of the surface normal (for a full derivation please see Ref. 60). We note that the kinetic term in Eq. (58) has been written as a time integral, which provides the symmetry with the configurational term discussed in relation to Eq. (56) and makes the (numerical) implementation identical to the configurational part. By taking the integral in Eq. (58), we arrive at

The surface evolution term is contained in the function $\u03d1t=H\xi t2\u2212xit2\u2212H\xi t1\u2212xit2\Lambda y(yit2)\Lambda z(zit2)$, where molecule positions are fixed and we count how many molecules have left or entered the volume due to the movement of the surface in the time interval between *t*_{1} and *t*_{2}. For the kinetic term, **r**_{i12} = **r**_{i2} − **r**_{i1} is the line of time evolution of a molecule *i* between *t*_{1} and *t*_{2} mirroring the configurational term’s intermolecular (IK) contour **r**_{ij} = **r**_{j} − **r**_{i}. Assuming the IK contour, the equation for a line is **r**_{λ} = **r**_{i1} + *λ***r**_{i12} where the value of *λ* at the point of crossing of the surface is *λ*_{k}. There is no closed form equation to get *λ*_{k} in general, but we can triangulate or split the surface into patches and use a ray-tracing process to obtain the point of intersection of a line and the surface.^{149} Once we obtain this crossing, it can be inserted in the following expression for the use of Eq. (59) in molecular simulations:

which is nonzero only if the line crossing the surface is between the start *λ*_{k} = 0 and finish *λ*_{k} = 1 of the line, and the point of crossing in the *y*- and *z*-direction, *y*_{k} and *z*_{k}, respectively, falls between the surface patch limits. Figure 9(c) illustrates the use of Eq. (59) on a general surface.

If we neglect the surface evolution term in Eq. (58), the general form of surface pressure is similar to the spherical and cylindrical pressure tensors presented in Eqs. (20), (22), (28), and (30). The surface normal $n\u0303x$ in Eq. (59) behaves like the radial ($R\u0302$) or azimuthal ($\varphi \u0302$) unit vectors, providing the pressure tensor aligned to the normal or tangential vector for a general surface varying in both *y* and *z*, i.e., $n\u0303x=n\u0303xy,z$. The function $H1\u2212\lambda k\u2212H\u2212\lambda k$ in Eq. (60) collects only the interactions crossing the surface as the function $H\lambda kH1\u2212\lambda k$ does in Eqs. (20), (22), (28), and (30). The extra localization of the Λ functions could have been included in the cylindrical and spherical forms of pressure if inhomogeneity along the surface were of interest. The general form presented in Eq. (59) can provide a detailed picture of molecular stacking and its effect on pressure near a much more complex surface, for example, the intrinsic interface, obtained by refitting each time to a set of molecules where a liquid meets a vapor.^{66} This uses sine and cosine functions with wavelengths chosen to allow fitting down to the molecular spacing. To capture instantaneous fluctuations about the spherical shape, spherical harmonics or similar functions could be used to capture the details of the molecular stress structure. These approaches quickly become cumbersome in general, so in Sec. IV F, we present a more pragmatic approach.

### F. Coordinate transforms

Evaluating the pressure form in Eq. (59) on a time-evolving surface requires an interface definition, together with an interaction calculation for every pair of molecules crossing that interface. As MD becomes a more common simulation tool, researchers are increasingly tackling more complex interface geometries, which require a more general fitted surface. This makes interface tracking increasingly more complex, and the resulting averaging grid can become impossibly deformed. Instead, in this section we discuss a process of collecting pressure values on a uniform grid and performing a transform afterward. Transforming the pressure requires two steps. The first is a rotation of the pressure tensor so the pressure is aligned with the normal to the surface. The second is a mapping so the pressure we collect can be obtained either at a distance from the surface or as a function moving along that interface. Through this process, we also highlight an important subtlety in pressure tensor studies: that the correct pressure tensor for a problem is dependent on measuring location and alignment.

As an example, we consider the NEMD boiling simulation shown in Fig. 11, a case of great practical interest.^{150–152} Here, a solid wall of tethered molecules is heated by a thermostat from the bottom and a phase change occurs in the liquid, starting at the bottom of a nanoscale square pore in the wall. The liquid is set up as a finite film with a large gas region above, and a bubble is allowed to nucleate and grow. The simulation is pseudo-two-dimensional with a nominal thickness in *z* and with periodic boundaries. The bubble grows initially inside the square pore, before extending beyond and forming a roughly circular shape. Figure 11(a) shows the density field where the yellow region is liquid with density of 0.7 in LJ units and the dark blue region is vapor with density of 0.05 in LJ units. There is no unique way to identify the liquid–vapor interface. In this work, a simple thresholding operation *ρ* > 0.3 (LJ units) is performed to identify the liquid, followed by a gradient operation to get the interface location where the density gradient is nonzero. These interface pixels are fitted using a least-square algorithm, with a circular arc to the part of the bubble above the wall. The fitted arc has a radius and angle used to rotate the pressure field *P*_{xx} in Fig. 11(b),^{153}

It is worth noting that the pressure tensor is a measure of the alignment of force with a given coordinate axis, so any rotation of the tensor is equally valid. Upon transform, this exposes a clear polar (tangential) pressure *P*_{θθ} around the interface in Fig. 11(c), similar to a hoop stress in a solid pressure vessel, which is what holds the bubble’s shape. It is this tangential contribution that will be significant in a Kirkwood–Buff surface tension calculation. In order to get *P*_{θθ} at a given radius, we need to average all values at the same radius over *θ*. Here, we apply this on the top half of the bubble, to avoid the near-wall region. The integration limits are shown in Fig. 11(c) by the double ended arrow between the angle at the bottom denoted by *θ*_{b} and the angle at the top denoted by *θ*_{t}. We use a mapping or projection from the Cartesian grid to a polar grid. Interestingly, instantaneous local fluctuations about the perfect circle can be seen in this projection in the blue line of tangential pressure in Fig. 11(d). We average along the interface, in this case between *θ*_{b} = −90° and *θ*_{t} = 90° for both radial and tangential pressure to obtain the plot in Fig. 11(e). This shows both the imbalance of radial pressure *P*_{RR}, which is driving the interface to grow, and the contribution to surface tension due to the tangential pressure *P*_{θθ}. The resulting pressure is similar to the one that would be obtained using Eqs. (20) and (22), but we have obtained it using data collected from the simulation on a uniform grid and by fitting the interface afterward. This technique is ideal for existing molecular dynamics software packages that do not have a rich selection of pressure calculation methods inbuilt, as discussed in Sec. V C, allowing uniform fields to be repurposed.

### G. Statistical uncertainty of different pressure methods

One of the main difficulties of pressure tensor calculations is the high level of noise relative to quantities such as *ρ*, ** u** and

*T*. This is clearly observed in Fig. 11(b) where a noisy pressure field gives a grainy appearance, which is absent from the density field in Fig. 11(a). It is worth observing that what we term “noise” here in NEMD are more concretely fluctuations about the time-evolving averaged quantities of interest, which tends to obscure them. In equilibrium systems, this noise can actually be the quantity of interest, such as diffusion or viscosity obtained from the Green–Kubo formula

^{43,44}or certain thermodynamic functions (e.g., heat capacity) that are a measure of fluctuations.

^{154}As discussed in Sec. IV C, these fluctuations might contribute to the turbulent eddies, and so they are important for an overall understanding of the flow. Hadjiconstantinou

*et al.*

^{155}estimated that collecting pressure is orders of magnitude worse in terms of statistics, which, they argue, makes coupling, passing averaged MD pressure values to be used as boundary conditions in a continuum solver, untenable. For NEMD simulations, getting good statistics becomes increasingly problematic as time-evolving events can depend on a chaotic trajectory, making it difficult to create a consistent ensemble. The example of boiling in Sec. IV F makes this clear, where nucleation occurs at different times in the members of an ensemble and the bubble growth proceeds in a varied and apparently stochastic way. A general discussion on the topic of pressure noise is difficult compared to other microscopic properties,

^{156}as the statistical requirements are case specific, depending on rate of time evolution in the system and choice of averaging volume size. Often, we are forced to use a coarse spatial resolution in order to provide sufficiently well-behaved pressure measurement.

To get a sense of the magnitude of this noise, Fig. 12 shows the distribution of measured pressure in a sub-volume of a large periodic simulation box of reduced density 0.8.^{157} The distributions are Gaussians for boxes of this size. Two cubic volume sizes are considered: The large volume has a reduced cell side length of *L* = 13.7*σ* and will contain about *N* = 2000 molecules, while the smaller volume has a reduced side length of 6.8*σ* with around *N* = 250. The VA pressure using Eq. (52) and the LP pressure using Eq. (56) are shown by solid and dotted lines, respectively. The average pressure is roughly the same for all cases with value ⟨*P*_{xx}⟩ ≈ 0.9 reduced units shown by the similar peak locations of the Gaussians in Fig. 12. Interestingly, the VA measures show a much lower spread in the distribution, e.g., standard deviations for the large volume $stdPVAxx=0.16$ vs $stdPLPxx=0.68$, both in LJ units. This means standard deviation is more than four times higher in the surface pressure measurement than the volume average for the same number of samples. For the smaller volume case, where the surface to volume ratio is improved by a factor of two in favor of the surface measurement, the standard deviation of the LP pressure is still about three times that of the volume average. One reason for this is apparent from the way pressures are obtained, as shown schematically in Fig. 12 (top right and bottom right). The VA scheme uses a continuous fraction of all interactions based on the length of line in the volume, which provides some smoothing as well as averaging from all interactions. Meanwhile, the LP pressure includes a contribution only if the surface is crossed, so fewer interactions are counted and these change abruptly as molecules move until interactions suddenly no longer cross a surface. Moreover, for $PLPxx$ only a single surface is used compared to the entire volume for the VA, further reducing the samples obtained in practice.

In conclusion, the VA pressure performs better at reducing noise, giving between three and four times lower standard deviation than a surface definition in this example. However, the conservative properties of the surface definition, as well as the ability to track complex geometries, make surface pressure preferable in some cases.

## V. CHALLENGES AND FUTURE DIRECTIONS

### A. Controversies over the microscopic pressure tensor

The microscopic pressure or stress tensor has been a controversial topic since the 1950s, largely due to the arbitrary contour involved in the formalism [Eq. (9)]. In fact, even the concept of the virial pressure introduced in 1870 by Clausius^{81} was not free of controversy. The form of the virial pressure was questioned in 1895, by a Basevi,^{158} who argued that by performing the integration in time that the kinetic and configurational parts should cancel. This argument was refuted by at least two articles^{159,160} with Gray^{160} noting that “*Basevi has, it seems to me, overlooked the fact that in the theorem it is the forces acting on each particle relatively to the assumed axes, and the corresponding motions that must be taken into account*.” Such confusions about the assumed axis and corresponding motion are apparently still a subject of confusion today together with the nonuniqueness of the pressure tensor itself.

The nonuniqueness problem of the microscopic pressure tensor is well known in the statistical mechanical community. It was implicit in the paper of Kirkwood and Buff^{70} in 1949, and the arbitrariness of the force acting across a surface element was then discussed in an appendix of the seminal paper by Irving and Kirkwood.^{1} Irving and Kirkwood’s warning did not attract much attention, although it was apparently noticed by Harasima *et al.*^{67,161} when they studied the surface tension of liquids. The nonunique nature of the microscopic pressure tensor became a focus in the 1980s, a time when molecular simulation was emerging as a technique to study complex systems that are critical in engineering, biology, and physics. Schofield and Henderson,^{2} for the first time, crystalized the ambiguity in the microscopic pressure tensor as an emergent property of the arbitrary contour shown in Eq. (9). Since then, many attempts have been made to find reasonable arguments and additional constraints^{78,162–167} that limit the choice of the contour. Until now, no consensus has been reached in the field,^{39,168} and there is no convincing justification for choosing one contour definition over the other for general cases. The nonuniqueness of the microscopic pressure tensor reflects the fact that there is no unequivocal way to assign a force [mechanical route, Eq. (9)] or potential energy [thermodynamic route, Eq. (37)] to a point **r** in space. In practice, the pressure is not measured at a point but over a volume or surface, so the functional form of this kernel also matters together with its location and shape. Any rotation or coordinate transform of a stress tensor changes the relative magnitude of the tensor components according to the orientation of the measurement. This has led to the concept of principal components in stress analysis,^{169} i.e., a rotation of the tensor, so that shear components are zero, providing an invariant or unique stress. From a fluid dynamical perspective, the local pressure tensor is subject to the so-called “gauge transform” where one can add a constant, or even the curl of any vector field, to the momentum density without affecting the system dynamics.^{2,61} This is because it is the gradient of the stress tensor that is well defined and not the stress tensor itself.

The other noteworthy controversy over the definition of the stress tensor was raised by Zhou in 2003,^{170} where the inclusion of the kinetic term in the stress definition was questioned. This controversy is in part due to the different definitions of pressure/stress tensor in the solid mechanics, thermodynamics, and fluid mechanics literature. In the solid mechanics literature, the Cauchy stress tensor is defined in terms of forces at zero temperature (i.e., no kinetic part). Often, the temperature dependence is included using an extra term in the continuum. Zhou’s argument, however, has been refuted by multiple studies. Admal and Tadmor^{78} showed that Zhou’s conclusions result from not considering the difference between absolute and relative velocities. They also demonstrated that the kinetic contribution to the stress is significant even for solid systems at a finite temperature. Subramaniyan and Sun^{171} showed the importance of temperature on stress in a thermoelastic study using molecular dynamics. Hoover *et al.*^{172} achieved an excellent agreement between atomistic mechanics and continuum mechanics provided that both kinetic and configurational contributions to the stress tensor are considered. It has been shown that the kinetic contribution to the stress is a direct consequence of the canonical transformation.^{8,173} We also note that the kinetic term is essential in the pressure and stress tensor definitions for thermodynamic consistency in the ideal gas limit (Sec. III). Away from thermodynamic equilibrium, kinetic pressure is defined in terms of peculiar velocity, the molecular velocity left after subtracting the streaming velocity of the flow. This contribution is essential in the pressure, gradients of which can drive flow in fluid dynamics, as well as in the shear stress between the molecules, which underpins fluid viscosity. The importance of kinetic pressure is most apparent when molecular simulation includes turbulent flow (Sec. IV C), where the kinetic contribution is shown in Fig. 7(b) to be as large as the configurational shear stress in turbulent flow and essential as a direct continuation of Reynolds stress below the scale of the measuring grid.

As molecular simulations find wider use both in fundamental science and, increasingly, in industrial applications, the necessity of finding an agreed definition of the microscopic pressure/stress tensor has become more critical than ever. We now consider a few promising contributions in this direction.

Motivated by the thermodynamic concept of pressure that is conjugate to a finite volume instead of to a point, Shi *et al.*^{35} showed that by spatially averaging the local (nonunique) pressure tensor over a small region of space of molecular dimensions, it is possible to define a coarse-grained (CG) microscopic pressure tensor that is unique and free from ambiguities in the definition of the local pressure tensor. In the case of fluids confined in a slit-shaped pore with *z*-axis perpendicular to the flat surface, such unique CG pressure tensor in the *k*th bin (*k* = 1, 2, 3, …) along the *z*-axis is given by

where the local (averaging) volume of the *k*th bin is Δ**r**_{k} = *L*_{x}*L*_{y}Δ*z*_{k}, *L*_{x} and *L*_{y} are the constant lateral dimension of the pore surface in the *x*- and *y*-directions, respectively, and Δ*z*_{k} is the characteristic length (width) of the *k*th bin that leads to a unique CG pressure tensor. This CG pressure tensor has the same appearance as the conventional VA definitions^{110,143,174,175} in Eqs. (51) and (52). However, the CG pressure is essentially different from the VA pressure in terms of the choice of the averaging region. For the CG pressure tensor, the averaging volume is constrained to give a unique CG pressure, while for the VA pressure tensor, the averaging volume is unrestricted and thus the resulting VA pressure can still be subject to the arbitrary choice of the integral contour.

To find the proper averaging region that will lead to a unique CG pressure, Shi *et al.*^{35} carried out the integration of the local tangential pressure over the *z*-direction analytically. They found that the contour path connecting particles *i* and *j*, upon integration, is fully dictated by a function $fC\lambda ij$, where *λ*_{ij} is the linearly scaled *z*-distance from particle *i* to particle *j* (assuming *z*_{i} < *z*_{j}); thus *λ*_{ij} = 0 amounts to *z* = *z*_{i} and *λ*_{ij} = 1 is *z* = *z*_{j}. Taking advantage of the symmetry of the contour path due to the indistinguishability of particles, ten contour definitions were designed, equivalent to ten unique functional forms for $fC\lambda ij$, as shown in Fig. 13(a). These ten contours include IK, H, IK-VR, and H-VR definitions that are introduced in Fig. 3. Using these ten types of contour definitions, they found that integrating the nonunique local tangential pressure over a certain *z*-distance leads to convergent integral results (within numerical uncertainty) that are independent of the arbitrary contour definitions [Fig. 13(b)].^{35} The CG pressure tensor that is defined between these convergence points [see characteristic length Δ*z*_{1}, Δ*z*_{2} marked in Fig. 13(b) for example] appears to be unique. Because these characteristic lengths are comparable to the thickness of the adsorbed layer [see density profile in Fig. 13(b)], the CG pressure tensor has direct physical significance, representing the microscopic pressure in an adsorbed layer.

The proposed CG scheme may serve as a unified solution toward a unique microscopic pressure tensor that is free from ambiguities in contour definitions, averaging volume and shape, and measurement locations. Future studies should focus on providing further simulation evidence and, ideally, rigorous mathematical proof to support the existence of this unique CG pressure tensor. Systems that are of particular interest for testing include those having curved interfaces and arbitrary geometries and those with moving reference frames.

Another possible aid to the problem of nonuniqueness are the conservative properties of the control volume form of pressure outlined in Sec. IV D and demonstrated in Fig. 10. This simply states that the pressure measured on surfaces that form an enclosed volume must exactly equal the momentum change inside.^{148,176} This is valid arbitrarily far from equilibrium and can be checked for any form of interaction contour, and with deforming or moving volumes,^{60} ensuring the pressure measurements satisfy Newton’s law. A different choice of volume or contour will simply redistribute the contribution to different terms in the tensor (as contributions are counted on different faces, for example). This exact equality between pressure and momentum change can help restrict the definition of the microscopic pressure tensor to the one that ensures that we do not violate Newton’s law.

### B. Complex systems interacting with many-body and long-range potentials

So far, we have only focused on the pressure tensor for systems of discrete particles that interact with short-range pairwise potentials. It is of practical interest to extend these formalisms to molecular and material systems that are controlled by more realistic interaction potentials. For example, in biological systems, the internal structure of the molecule is mainly dictated by many-body intramolecular interactions, such as angular, torsion, and improper potentials. As for intermolecular potentials, in addition to the short-range dispersion interactions, long-range Coulombic interactions are commonly present in the system due to the uneven distribution of charges in the molecule.

In general, the pressure tensor can be expressed in two equivalent forms: an *atomic* form and a *molecular* form. The *atomic* pressure tensor is defined in such a way that all forces on each atom should be evaluated explicitly;^{177–179} these include contributions from intermolecular interactions and intramolecular interactions (e.g., bonded and bending forces, and constraint forces imposed by the SHAKE algorithm^{180}). The *molecular* pressure tensor takes the molecule to be a rigid body (i.e., rigid body approximation^{181}), and only (nonbonded) intermolecular interactions contribute to the pressure tensor.^{74,83,177–179} Compared to the atomic pressure tensor, this molecular formalism neglects the intramolecular interactions and thus loses the mechanical details of the internal structure of the molecule. In addition, implementation of the molecular pressure tensor requires defining the center of mass (COM) at which to localize the momentum of a molecule.^{179} This definition of COM is straightforward for small molecules, but it is by no means intuitive for solid materials or polymer chains as they are usually modeled as infinitely extended structures under periodic boundary conditions. Therefore, the atomic pressure tensor is generally considered to be more useful, and it has been widely applied to study the mechanical properties of crystalline polymers,^{179} lipid bilayers,^{4,6,7} liposomes,^{182} proteins,^{8} metals, and alloys.^{183} In cases where molecules are small and behave like rigid bodies (such as water, methane, or ethylene), however, the molecular pressure tensor is more convenient. This is because the molecular formalism circumvents the evaluation of rigid constraints^{74,83,184} and thus potentially avoids the complexity of the anisotropic kinetic term in the atomic pressure tensor near interfaces.^{185} Sega *et al.*^{185} showed that, without considering such an anisotropic kinetic term for rigid molecules, the surface tension, calculated with only the configurational part of the atomic pressure tensor, tends to deviate from the true value by about 9%–15% for water systems.

The technical development for the atomic pressure/stress tensor has been centered on the decomposition of a many-body potential into pairwise components, so that the typical pairwise formalism based on Eq. (9) can be implemented. This force decomposition is not unique,^{186} and a number of decomposition methods have been proposed.^{4,78,175,182,187–190} A noticeable development is the so-called central force (CF) decomposition by Admal and Tadmor.^{78} In a CF decomposition, the total force on atom *i* due to a *m*-body potential $Umrm$ is written as a summation of pairwise central forces (i.e., forces that are parallel to the vector **r**_{ij}),

where **r**^{m} ≡ **r**_{1}, **r**_{2}, …, **r**_{m} is a collection of particle positions in a *m*-body cluster. Figure 14 illustrates the CF decomposition for a three-body potential. While some popular force decomposition methods violate the balance of linear^{175,189,191} or angular^{4} momentum, the CF decomposition yields a symmetric stress tensor by construction and satisfies the balance of both linear and angular momentum.^{191} In practice, unknown parameters *F*_{ij} are obtained by solving a system of linear equations given in Eq. (63). Here, the number of independent force equations is 3*m* − 6 (as forces satisfy the conservation of linear and angular momentum), and the number of unknown pairwise central terms is $mm\u22121/2$. For *m* = 3, 4, solving this system of linear equations is a well-posed math problem. For potentials beyond four-body interactions, however, the number of unknown parameters is larger than the number of independent equations, and the CF decomposition becomes nonunique.^{7} Torres-Sánchez *et al.*^{190} developed a covariant central force (cCF) decomposition based on the Doyle–Ericksen relation of continuum mechanics, rather than on the statement of balance of linear momentum, as in the classical Irving–Kirkwood–Noll approach. The cCF decomposition is consistent with the CF approach for three- and four-body potentials but allows for many-body interactions of arbitrarily high order.^{190,191} For pairwise potentials, all force decomposition schemes result in the same pressure or stress tensor. The nonunique scheme of the force decomposition is related to the nonuniqueness of the local pressure tensor due to arbitrary contour. One can argue that, using the thermodynamic route (Sec. III B), no force decomposition is needed, but the problem now becomes the ambiguity in assigning many-body potential energy into a local space.

The other technical challenge for both the atomic and molecular representations of the pressure tensor involves the consideration of long-range Coulombic interactions. Unlike the short-range LJ potential, the Coulombic potential decays very slowly in space, and it is impossible to use a simple tail correction to account for the missing long-range part.^{192} Compared to the direct Coulomb sum, which scales as $ON2$ (where *N* is number of atoms or particles in the system), the Ewald summation method^{193} is the standard method used in molecular simulations to efficiently handle the Coulombic interactions, scaling as $ON3/2$. Better efficiency can be achieved by modern, mesh-based Ewald methods,^{194} which scale as $ONlogN$. While the algorithm for computing the bulk (macroscopic) pressure tensor in the presence of the Coulombic interactions is well established in the field,^{178,195} it is still a challenge to handle the Coulombic or any long-range interaction in an efficient manner for the local pressure tensor. Previous studies have focused on finding a suitable scheme for assigning the local force that is compatible with the Ewald summation method. Since the *k*-space (Fourier space) part of the Coulombic energy in the Ewald method is the most computationally efficient in a non-pairwise form, the Harasima contour definition turns out to be a better choice than the IK contour. This can be understood by the delta function in the Harasima formulation [Eq. (16)]; the delta function indicates that the configurational part only contributes to the tangential pressure at planes where molecules are present. This feature allows the per-atom form of the *k*-space energy term to be naturally incorporated into the Harasima formulation. Compatibility of the Ewald summation method with the Harasima contour has been developed for the local tangential pressure across a planar interface^{6,184,196} and for the local axial pressure in a cylindrical geometry.^{74} Nevertheless, it is still possible to use the IK contour with the Ewald method. Hatch and Debenedetti^{8} successfully captured the full Coulombic energy using the IK contour definition by writing the Ewald sum in an explicit pairwise form. However, as expected, a pairwise form of the Ewald sum is computationally expensive. A computationally amenable alternative using the IK contour is to consider the Coulombic potential up to a certain cutoff radius,^{7,182} but such a treatment cannot guarantee a consistent pressure profile because the Coulombic potential was treated differently (with Ewald-based method) in the molecular simulations and (with simple cutoff) in the pressure calculation. It is possible to replace the bare Coulombic potential with a damped, shifted-potential^{197} or a shifted-force^{198} one, which is generally called the Wolf potential. The Wolf potential can reproduce the energetics and dynamics of various systems to an acceptable accuracy compared to the (exact) Ewald method. However, the error that might be introduced to the pressure tensor by this (approximate) Wolf potential is still unclear. We note that the preference of the Harasima contour in the treatment of the Coulombic interactions does not imply that the Harasima contour is more correct than the IK contour, but the Harasima/Ewald method clearly has advantages in its computational efficiency.

LAMMPS command . | Explanation . |
---|---|

Compute stress/atom | Computation of per-atom stress tensor. A virial contribution produced by a m-body potential is equally assigned to each atom in the set, e.g., 1/4 of the dihedral virial to each of the four atoms. We note that this function is commonly used for visualization purpose; caution should be exercised when interpreting this per-atom value in a local fashion, which has some similarities to the IK1 approximation [Eq. (7)]. This command works for long-range^{178} and many-body interactions.^{206} |

Compute pressure | Computation of a scalar pressure and a global pressure tensor of the entire system^{206} |

Compute pressure/uef | Computation of the pressure tensor in the reference frame of the applied flow field |

Compute stress/mop^{a} | Computation of local stress tensor using the method of planes [Eq. (54)]; specifically, it computes three components in directions αx, αy, and αz, where α is the direction normal to the plane.^{61} The profile of the stress can be computed with “compute stress/mop/profile” command. |

Compute stress/Cartesian^{a}^{,}^{b} | Computation of coarse-grained profiles of the diagonal components of the local stress tensor in Cartesian coordinates; the output stress tensor is averaged over a small local volume [Eq. (52) with a slab-like local volume].^{174} |

Compute stress/spherical^{a}^{,}^{b} | Computation of profiles of the diagonal components of the local stress tensor in spherical coordinates [Eq. (52) with a spherical-shell local volume]^{174} |

Compute stress/cylinder^{a}^{,}^{b} | Computation of profiles of the diagonal components of the local pressure tensor in cylindrical coordinates [Eqs. (28), (30), and (32)];^{21} this command does not consider periodic boundary conditions, so the system should be large enough to ensure the boundary effect is negligible |

LAMMPS command . | Explanation . |
---|---|

Compute stress/atom | Computation of per-atom stress tensor. A virial contribution produced by a m-body potential is equally assigned to each atom in the set, e.g., 1/4 of the dihedral virial to each of the four atoms. We note that this function is commonly used for visualization purpose; caution should be exercised when interpreting this per-atom value in a local fashion, which has some similarities to the IK1 approximation [Eq. (7)]. This command works for long-range^{178} and many-body interactions.^{206} |

Compute pressure | Computation of a scalar pressure and a global pressure tensor of the entire system^{206} |

Compute pressure/uef | Computation of the pressure tensor in the reference frame of the applied flow field |

Compute stress/mop^{a} | Computation of local stress tensor using the method of planes [Eq. (54)]; specifically, it computes three components in directions αx, αy, and αz, where α is the direction normal to the plane.^{61} The profile of the stress can be computed with “compute stress/mop/profile” command. |

Compute stress/Cartesian^{a}^{,}^{b} | Computation of coarse-grained profiles of the diagonal components of the local stress tensor in Cartesian coordinates; the output stress tensor is averaged over a small local volume [Eq. (52) with a slab-like local volume].^{174} |

Compute stress/spherical^{a}^{,}^{b} | Computation of profiles of the diagonal components of the local stress tensor in spherical coordinates [Eq. (52) with a spherical-shell local volume]^{174} |

Compute stress/cylinder^{a}^{,}^{b} | Computation of profiles of the diagonal components of the local pressure tensor in cylindrical coordinates [Eqs. (28), (30), and (32)];^{21} this command does not consider periodic boundary conditions, so the system should be large enough to ensure the boundary effect is negligible |

^{a}

The command works only for short-range pair interactions; i.e., if any bond, angle, dihedral, etc., contributions and *k*-space contributions (in Ewald summation method) for the long-range Coulombic interactions are present in the system, the results will be incorrect.

^{b}

The calculation is based on the IK contour definition.

Future breakthrough in this regard would rely on a better understanding of the physical nature of the microscopic pressure tensor in the presence of multi-body and Coulombic interactions. The possibility to define a unique CG pressure tensor (Sec. V A) will benefit the calculation of Coulombic contributions to the pressure tensor as we are free to choose the most convenient contour (e.g., Harasima contour) without worrying about the physical and numerical arbitrariness of the final results. Due to the convenience in calculating the short-range interaction, it is also worthwhile to systematically evaluate the uncertainty that a spherically truncated or damped Coulombic potential will introduce to the pressure in general, in comparison to the exact treatment using direct Coulomb sum (without potential cutoff) or more efficient Ewald-based methods.

### C. Software and computational tools

Currently, the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) molecular dynamics software provides native functionalities for calculating the local and global stress/pressure tensor on the fly, and these *compute* commands are summarized in Table I. In contrast, the Nanoscale Molecular Dynamics (NAMD) software only provides a basic function to compute the pressure tensor profile along the *z*-axis for systems having a planar geometry.^{199} This specific implementation adopts the Harasima contour definition, which allows an efficient computation of the Coulombic contribution to the local pressure tensor based on the Ewald method (see Sec. V B for details).^{6} The limited pressure or stress tensor functionality in major molecular simulation software has motivated the development of dedicated analysis tools. Nakamura *et al.*^{72,200} prepared a patch file for the LAMMPS software that enables the calculation of the local pressure tensor in Cartesian and spherical coordinates. Vanegas *et al.*^{190,201} developed a computational tool “GROMACS-LS” for the GROMACS software to calculate the local stress in molecular systems. Admal *et al.*^{78,202} developed a post-analysis program “MDStressLab” that takes the input data (particle coordinates, velocities, species, etc.) in a general format, so it is compatible with different simulation software provided that the scripts for converting software-specific file format are available. All of these dedicated analysis tools are capable of calculating the 3D pressure or stress tensor field (as a volume average value, Hardy definition) in systems having arbitrary geometries. They differ in the availability of pressure/stress tensor definitions, force decomposition schemes for many-body potentials, and supported interaction potentials. Interested readers are referred to the program’s user guide for details. Other pressure tensor codes that were designed for specific coordinate systems and intermolecular potentials are also available.^{203–205}

At the moment, all available computational tools are more or less limited in their definitions of the local pressure tensor, in applicability to certain system geometries, and in the availability of intermolecular potentials. Therefore, there is a strong motivation to develop a general-purpose software in the future with a well-documented user guide. As users may have personal preference for molecular simulation packages, it is imperative to develop compatibility of this pressure analysis software with different simulation packages. This could be realized through support for a range of input formats (or at least, scripts for converting file formats) or a generalized cross-language functional interface that can be called from any of the codes during the force calculation to tally the pressure. Calculating the local pressure tensor is computationally expensive, therefore parallelization of the computation using high-performance central/graphics processing unit (CPU/GPU) would be essential in the future.

The development of a standard computational package for the microscopic pressure tensor will also standardize the measurement of pressure and enable the generation of a database of high-fidelity pressure/stress data for systems of both engineering and scientific interest. The availability of a large amount of data organized in a curated database will help advance the process optimization and materials design using machine learning and data science.

### D. Experimental measurements of microscopic pressure tensor

At the moment, it is still a challenge to directly measure or estimate the pressure or stress tensor at the nanoscale from experiments, and nearly all microscopic pressure or stress tensor results reported in the literature so far are theoretical values. The difficulties come from the limitation in experimental techniques to approach the molecular level stress and from the nonuniqueness in the definition of the microscopic pressure tensor. Advances in either side will significantly benefit the other.

Gubbins *et al.* have demonstrated that the pressure tensor in a molecular scale system could be estimated from experimental input in simple equilibrium systems, such as for fluids adsorbed in a carbon slit-shaped pore. The normal pressure component, *P*_{N}, in such a system is a constant, and Śliwińska-Bartkowiak *et al.*^{207,208} demonstrated that it is possible to estimate this pressure by measuring the resulting changes in the interplanar distance of the activated carbon fibers, using X-ray diffraction. The in-pore normal pressure can then be estimated using Young’s equation, provided that the transverse compressive modulus is known. Figure 15 shows that an agreement has been achieved between the simulated normal pressure and the experimental estimations, within the (rather large) uncertainties of the latter, for CCl_{4} and H_{2}O adsorbed in carbon slit pores.

For the tangential pressure, *P*_{T}, the experimental estimation is more challenging because the local tangential pressure is nonunique, and the force does not act directly on the adsorbent material but in a direction parallel to the wall. Thus, a noninvasive method is needed. Molecular simulation and experimental results show that for adsorbates that wet the pore walls, the adsorbed layers of molecules very near to the wall are quasi-two-dimensional; this being particularly pronounced for the contact layers next to the pore walls. This observation enabled the development of a “2D route” to the effective tangential pressure inside a single adsorbed layer.^{209,210} In this 2D route, the behavior of a single adsorbed layer near the surface is related to that of a strictly 2D reference film by projecting the center of mass positions of the molecules in the layer onto the surface plane. The 2D pressure *P*_{2D} (in units of force per unit length) in the reference film is then mapped back to the 3D pressure (in units of force per unit area) by being divided by an effective length scale *l*_{eff} in the direction normal to the surface. The effective tangential pressure estimated by the 2D route is^{210}

where *T* is temperature and *ρ*_{2D} is the 2D density of the adsorbed film. The 2D pressure *P*_{2D} is a function of *T* and *ρ*_{2D}, and the relation can be established by a 2D equation of state.^{209} Although *P*_{2D} is well defined, $P2DT$ is not unique due to the arbitrary choice of *l*_{eff}. To decide on a sensible choice of *l*_{eff}, it is instructive to rewrite Eqs. (64) as a spatial (volume) average,^{210} i.e., $P2DT\u2248\u222bleffPTzdz/leff$. Based on the results in Fig. 13, it is motivating to choose *l*_{eff} to be the characteristic length Δ*z*_{k} that is comparable to the thickness of an adsorbed layer, so that the spatial average appears to be unique and has a clear physical meaning. In this 2D route, experimental input parameters are *T*, *ρ*_{2D}, and *l*_{eff}. The 2D density *ρ*_{2D} can be estimated from molecular simulations or adsorption theories or obtained from particular experiments such as small-angle neutron diffraction. The effective thickness of the layer *l*_{eff} can be estimated from molecular simulations and theories or from optical sensing experiments.^{211} Limitations of the 2D route include the following: (1) It neglects the interactions between the layer of interest and its neighboring layers; (2) It may fail when the adsorbed layer deviates from a quasi-2D structure, for example, in weakly wetting systems or for materials with a rough surface.

Advanced experimental techniques have been developed to estimate the scalar pressure at very small scales. These techniques could serve as a foundation for future routes to approach the microscopic pressure in its correct tensor form:

**Direct methods:**The pressure can be directly estimated by measuring the mechanical response of the system to a physical probe, such as an atomic force microscopy (AFM) tip.^{212–214}This method has been commonly adopted to estimate the scalar pressure inside a nanobubble.^{212,214}For example, the interaction force between the nanobubble and the AFM tip can be dynamically captured based on the AFM cantilever’s instantaneous deflection. The internal pressure of the nanobubble can be obtained by fitting the recorded force–displacement curve to a theory that separates the internal pressure of the confined fluids from the elastic deformation of the solid materials.^{212}Another possible way is to relate the pressure to the elastic properties of the system, such as the elastic modulus, $KT=\beta T\u22121$, where the isothermal compressibility*β*_{T}is defined as(65)$\beta T=\u22121V\u2202V\u2202PT.$*β*_{T}of the confined fluids can be measured by ultrasonic experiments.^{215,216}However, the measured elastic property is commonly interpreted as a macroscopic (scalar) quantity averaged over all directions and the system volume, instead of as a tensor.**Indirect methods:**The scalar microscopic pressure can be sensed by molecular probes that are (mechanically, optically, thermally, etc.) susceptible to the pressure change of the environment. Vasu*et al.*^{28}studied the effect of confinement between two graphene layers on molecules that are susceptible to deformation under pressure. By comparing Raman spectra for the confined molecules with those for the same molecules in the bulk phase, they estimated the “van der Waals” (scalar) pressure in the confined phase to be 1–1.5 GPa. Another relevant example is using molecular rotors as sensors to measure the local viscosity of a fluid under extreme confinement conditions.^{217}The probe molecular rotor changes its fluorescent properties due to the pressure and shear of the surrounding fluids. This information could allow an estimate of the pressure tensor in a fluid away from equilibrium at very small scales.

An immediate challenge is the lack of rigorous relations that connect the experimental measurements to the microscopic pressure tensor in statistical mechanics. Future studies should focus on establishing these fundamental relations to bridge the gap between experiments and theories.

## VI. CONCLUDING REMARKS

In this Perspective, we have reviewed several routes to calculate the microscopic pressure tensor (equivalent to the negative stress tensor) in both equilibrium and nonequilibrium systems. These formulas can be generally divided into two types, depending on whether they were derived via mechanical or thermodynamic routes. The mechanical (or “virial”) route follows the mechanical concept of “the force acting across a surface,” and it can be used in both equilibrium and nonequilibrium systems. By contrast, the thermodynamic route uses the thermodynamic definition of pressure, which is the negative of the change in the Helmholtz free energy with respect to volume. The thermodynamic route can only be used in equilibrium fluid systems where no shearing is present, but it is arguably preferable to the mechanical route for systems interacting with complex (multi-body) potentials.

We have categorized available pressure equations into different forms, based on where and how the microscopic pressure tensor is measured. These include macroscopic (bulk), pointwise, volume, and surface forms. We then attempt to show the underlying connection between the different forms, highlighting the inherent assumptions with the limitations of each of these choices. The equations and connections between different forms are summarized in Fig. 16 for the configurational part of the pressure tensor and in Fig. 17 for the corresponding kinetic part.

We have also pointed out four aspects that currently face challenges and need further investigations. In brief, they are as follows:

**Historical controversies over the definition of the microscopic pressure tensor:**Controversies are centered on the nonuniqueness of the microscopic pressure tensor at a point in space, resulting from the fact that the forces between molecules do not act at a unique point in space; this difficulty manifests itself in the equations as the arbitrary contour involved in the calculations. However, coarse-graining over a relatively small spatial region of space may result in a well-defined pressure tensor, as has been shown recently for a simple system.^{35}A breakthrough in this regard may open the door to a thermodynamically and mechanically consistent picture of nanoscale systems.**Difficulties with many-body and long-range potentials:**This technical difficulty is outstanding for complex systems, such as those in biology. It is suggested that future research focus on better understanding the physical nature of the pressure tensor in the presence of many-body interactions and on the development of convenient, accurate, and efficient algorithms to account for long-range interactions in the pressure tensor.**Inadequate software and computational tools for the calculation of the local pressure/stress tensor:**To our knowledge, no general-purpose software is available. Such software and computational tools are essential to avoid confusion and to overcome the knowledge barriers for nonexpert users and are needed to accelerate the process optimization and materials design using machine learning and data science. The ease of calculating the virial pressure results in its use when running exiting software packages, often in cases where it is invalid (e.g., confined flows).**Lack of experimental methods to measure the pressure tensor at the nanoscale:**Advances in this regard require a combined effort from both the experimental and theoretical/computational communities. A breakthrough in measuring the microscopic pressure tensor will enable the determination of a wide range of thermodynamic and transport properties at the nanoscale in the correct tensor form.

The microscopic pressure tensor is a pivotal property for a wide range of disciplines and technologies, including fluid dynamics, solid mechanics, biophysics, thermodynamics, nucleation and crystallization, chemical manufacturing, and chemical separations. We invite both experimentalists and theorists to contribute to this field to enable a thorough and consistent understanding of tensorial features of inhomogeneous systems.

## ACKNOWLEDGMENTS

K.S. thanks Dr. Yun Long for helpful discussions about the thermodynamic route for curved interfaces. K.S. also thanks Dr. Harold (“Wick”) Hatch for pointing out important references on the canonical transformation. E.S. acknowledges funding from the EPSRC Project No. EP/S019545/1 and is grateful to the UK Materials and Molecular Modeling Hub for computational resources, which is partially funded by EPSRC (Grant Nos. EP/P020194/1 and EP/T022213/1). K.E.G. thanks the U.S. National Science Foundation for partial support under Grant No. CBET-1603851. E.E.S. thanks the U.S. National Science Foundation for support under Grant No. CBET-1855465.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Kaihang Shi**: Conceptualization (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Edward R. Smith**: Conceptualization (equal); Funding acquisition (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Erik E. Santiso**: Conceptualization (equal); Funding acquisition (equal); Writing – review & editing (equal). **Keith E. Gubbins**: Conceptualization (equal); Funding acquisition (equal); Resources (equal); Supervision (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding authors upon reasonable request.

### APPENDIX: ADDITIONAL DERIVATIONS

#### 1. Linking the IK contour pressure to the MoP form

In this appendix, we derive the mathematical relationship between the configurational part of the IK contour pressure [Eqs. (13) and (15)] and the MoP pressure [Eq. (54)]. This is significant as the MoP pressure is derived in Fourier space and valid away from equilibrium, while the IK contour pressure is obtained by integrating (averaging) over the *x*- and *y*-directions and generally used only in an equilibrium system. Considering only the configurational part of Eq. (13) and splitting Eq. (15) into two the different tangents, with averaging notation ⟨⋯⟩ omitted for simplicity, we have

For comparison, we rewrite the MoP pressure in Eq. (54) for three components on the *z*-surface (the surface that is normal to the *z*-direction); the configurational term gives

Noting that the scalar force in Eqs. (A1)–(A3) can be related to its vector form by **F**_{ij} = *F*_{ij}**r**_{ij}/*r*_{ij} and the projection of the pair force in the *α*-direction (*α* = *x*, *y*, *z*) is *F*_{ijα} = *F*_{ij}*α*_{ij}/*r*_{ij}, we write Eqs. (A1)–(A3) in a unified vector form as

where **F**_{ij}◦**r**_{ij} denotes the element-wise product of force vector **F**_{ij} and distance vector **r**_{ij}. The Heaviside function expressed as a product gives identical behavior to the difference of two Heaviside functions,

Using the property $H(ax)=1/2(sgnasgnx+1)$ to transform Eq. (A6) to signum functions, we get

Using the definition of the signum function $sgn(zij)=zij/zij$, the product is $sgn1/zijsgnzij=1$, so Eq. (A1) can be rewritten as

Here, we obtain the equivalence between the normal MoP pressure on a *z*-surface and the corresponding *zz*-component in the IK contour pressure tensor (and note *zz*-component is independent of the contour definition).

The tangential components, however, highlight a fundamental difference between the IK contour pressure and the MoP pressure. For example, considering the *yy*-component of Eq. (A3) and rewriting it in the same way as Eq. (A8), it is trivial to show

where we have an extra factor of *y*_{ij}/*z*_{ij} when compared to the normal MoP pressure on a *y*-surface [similar to Eq. (56)]. The IK contour pressure therefore uses extra molecular information by taking the *y*_{ij} distance to calculate the direct *y*-pressure on a *z*-plane. This is, however, a departure from the Cauchy definition of pressure as shown in Fig. 4(a), which requires the three different normal pressures to be defined on orthogonal planes. Other form of the pressure, for example, the VA pressure, uses this same quantity *F*_{ijy}*y*_{ij} to get *P*_{yy}, but in the VA form, this is weighted by the fraction of interaction *l*_{ij} inside the volume. In contrast, the IK contour form accumulates interactions crossing a *z*-plane and weights these by the length of interaction in the surface normal direction *z*_{ij} as shown in Fig. 18. We also note that averaging the IK contour pressure over a volume should lead to the corresponding VA pressure form.

#### 2. The Noll form of pressure

The Noll reformulation replaces the Dirac delta function $\delta r\u2212ri\u2192f|ri=r$. This notation denotes an integral of probability density function *f* over all phase space except **r**_{i}, which is then replaced by **r**. Conditions are then placed on *f*, including a similar condition to the phase space bounded assumption of Irving and Kirkwood. Equation (4) can be reformulated using Noll's lemma to replace the IK operator as

This states that the stress tensor at point **r** is a superposition of expectation forces from all possible bonds that might pass point **r**, while the integral over *λ* “slides” from **r**_{i} to **r**_{j} along a vector ** z** between the molecules (see Ref. 78 for a sketch). This form is potentially more general than the pair potential that is assumed by Eq. (8) where

**=**

*z**ℓ*can be a general contour.

^{78}For the case of pairwise interactions,

**z**=

**r**

_{ij}=

**r**

_{j}−

**r**

_{i}, we can rearrange both equalities in the brackets to show they are the same condition

**r**

_{j}− (1 −

*λ*)

**r**

_{ij}=

**r**

_{i}+

*λ*

**r**

_{ij},

which is similar to the form of line integral used throughout this work. Note in general that we depart from the approach of Noll as the equations we consider are derived without the use of an ensemble average^{41,110,112} as discussed in Sec. IV.