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.

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 = PK + PC, 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 kB 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, PK, is well defined, as noted above, the configurational part, PC is not because the intermolecular forces themselves are nonlocal. Thus, while the force between molecules i and j, located at positions ri and rj, 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 dS 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 σr,t is usually used in place of the pressure tensor.1 These two tensors differ only in sign,1,3

σr,t=Pr,t.
(1)

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–liquid11–13 or liquid–solid14,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 energy14 and for understanding the distinct structure of the nucleus15 [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 cylindrical21 interfaces). For confined systems, the knowledge of the microscopic pressure tensor paves the way for understanding phase transitions in nanopores22–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 bulk28,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.

FIG. 1.

Number of papers that include the keywords “molecular dynamics” or “Monte Carlo” and either “stress tensor” or “pressure tensor” in the last 50 years, obtained from Google Scholar.58 

FIG. 1.

Number of papers that include the keywords “molecular dynamics” or “Monte Carlo” and either “stress tensor” or “pressure tensor” in the last 50 years, obtained from Google Scholar.58 

Close modal
FIG. 2.

Illustration of applications for the microscopic pressure/stress tensor. (a) Local lateral (PL) and normal stress (PN) field in lipid membrane showing strong lipid unsaturation in the tail group regions. Adapted with permission from Vanegas et al., J. Chem. Theory Comput. 10(2), 691–702 (2014). Copyright 2014 American Chemical Society. (b) Tangential (PT) and normal (PN) pressure profile in a solid–liquid nucleation system, showing strikingly lower pressure in the solid nucleus than the liquid environment. Data taken from Ref. 15. (c) Snapshots and tangential and normal pressure profiles of Lennard-Jones argon (pink) confined in an atomistic carbon slit pore (cyan) at 87.3 K and 1 bar bulk pressure with pore width of 51 Å. High tangential pressure near the surface indicates the strong compression effect inside the physisorbed layers, shedding light on the understanding of high-pressure phenomena in more complex systems. Adapted from Long et al., J. Chem. Phys. 139(14), 144701 (2013) with the permission of AIP Publishing. (d) The liquid–vapor–solid moving contact line: molecules from a molecular dynamics simulation shown with a Chebyshev function fitted to the liquid–vapor interface to be used in the stress calculation. The setup is a liquid bridge between two sliding molecular walls,59 where the interface is split into bilinear patches and the stress is obtained as the force acting over each patch of the area.60 

FIG. 2.

Illustration of applications for the microscopic pressure/stress tensor. (a) Local lateral (PL) and normal stress (PN) field in lipid membrane showing strong lipid unsaturation in the tail group regions. Adapted with permission from Vanegas et al., J. Chem. Theory Comput. 10(2), 691–702 (2014). Copyright 2014 American Chemical Society. (b) Tangential (PT) and normal (PN) pressure profile in a solid–liquid nucleation system, showing strikingly lower pressure in the solid nucleus than the liquid environment. Data taken from Ref. 15. (c) Snapshots and tangential and normal pressure profiles of Lennard-Jones argon (pink) confined in an atomistic carbon slit pore (cyan) at 87.3 K and 1 bar bulk pressure with pore width of 51 Å. High tangential pressure near the surface indicates the strong compression effect inside the physisorbed layers, shedding light on the understanding of high-pressure phenomena in more complex systems. Adapted from Long et al., J. Chem. Phys. 139(14), 144701 (2013) with the permission of AIP Publishing. (d) The liquid–vapor–solid moving contact line: molecules from a molecular dynamics simulation shown with a Chebyshev function fitted to the liquid–vapor interface to be used in the stress calculation. The setup is a liquid bridge between two sliding molecular walls,59 where the interface is split into bilinear patches and the stress is obtained as the force acting over each patch of the area.60 

Close modal

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 equilibrium43,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 simulation52,53 where stress coupling is useful in both fluid simulation54 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 boundaries56 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.

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

Pr,t=PKr,t+PCr,t,
(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,

PKr,t=i=1Npipimiδrri,
(3)

where the angular bracket ⟨⋯⟩ denotes ensemble average, N is the total number of particles in the system, pi is the momentum of the particle i and pipi is an outer product, mi is the mass of particle i, and δr=δxδyδz is the delta function for a vector position in a Cartesian coordinate system. Here, the momentum may include a streaming component, pi/mi=ṙiu, due to the streaming velocity u at point r being nonzero, where ṙi 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 PC can be obtained from the tensor product of pair forces between particles, Fij, and the line of interactions rij (rij = rjri),

PCr,t=12i,jNFijrijOijδrri,
(4)

where the pre-factor 1/2 accounts for double counting and the notation i,jNFij=i=1NjiNFij has been introduced as shorthand for the double summation over all pairs of interactions. The Oij term is known as the Irving–Kirkwood (IK) operator,1,61

Oij=112!rijr++1n!rijrn1+,
(5)

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

δrriδrrj=rijrOijδrri.
(6)

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

PIK1Cr,t=12i,jNFijrijδrri.
(7)

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 (Oij = 1) is exact.61 For inhomogeneous fluids, such as those confined in nanopores or near interfaces, the IK1 approximation violates the mechanical equilibrium condition63,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

rijOijδrri=Cijδrd,
(8)

where Cij 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 

PCr,t=12i,jNFijCijδrd.
(9)

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 dS, Harasima67 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 dS if the joining straight line (contour) between two molecules passes through dS. 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 Buff70 and later referred to as the Harasima definition. It says a pair force contributes to the pressure tensor at the surface dS if one of the molecules lies in the cylinder whose base is dS (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 dS. 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 = ri + λrij with 0 ≤ λ ≤ 171 and Eq. (9) becomes

PIKCr,t=12i,jNFijrij01δrriλrijdλ.
(10)
FIG. 3.

Four possible definitions of contour connecting the particles i and j for the local pressure tensor in a system of planar geometry. (a) Irving–Kirkwood (IK) definition. (b) Harasima (H) definition. (c) A variation of the Irving–Kirkwood (IK-VR) definition. (d) A variation of the Harasima (H-VR) definition. Due to the indistinguishability of particles, contours Cij (in black, from particles i to j) and Cji (in orange, from particles j to i) are equivalent and symmetric. Since the Cij and Cji overlap for the IK and the H-VR definitions, only the contour Cij is plotted for clarity. All contours are projected onto the xz-plane, and the z-direction is perpendicular to the planar surface. These possible contour definitions illustrate the nonunique nature of the local pressure tensor. Reprinted from Shi et al., J. Chem. Phys. 154(8), 084502 (2021) with the permission of AIP Publishing.

FIG. 3.

Four possible definitions of contour connecting the particles i and j for the local pressure tensor in a system of planar geometry. (a) Irving–Kirkwood (IK) definition. (b) Harasima (H) definition. (c) A variation of the Irving–Kirkwood (IK-VR) definition. (d) A variation of the Harasima (H-VR) definition. Due to the indistinguishability of particles, contours Cij (in black, from particles i to j) and Cji (in orange, from particles j to i) are equivalent and symmetric. Since the Cij and Cji overlap for the IK and the H-VR definitions, only the contour Cij is plotted for clarity. All contours are projected onto the xz-plane, and the z-direction is perpendicular to the planar surface. These possible contour definitions illustrate the nonunique nature of the local pressure tensor. Reprinted from Shi et al., J. Chem. Phys. 154(8), 084502 (2021) with the permission of AIP Publishing.

Close modal

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 field72 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 (spherical73 and cylindrical74 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 

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.

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 

P=0
(11)

which has two implications:76 (1) the tangential pressure PT has no local gradient in directions parallel to the xy-plane that induce flow anywhere; (2) the normal pressure PN should be a constant throughout the system,

PN=Pzz=const.
(12)

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).

FIG. 4.

Pressure tensor definition in systems having different geometries. (a) Cartesian coordinates for a planar geometry. (b) Spherical coordinates for a spherical geometry. (c) Cylindrical coordinates for a cylindrical geometry. (d) An example of a more general coordinate system in terms of vector surface normal nα (α = x, y, z) where one (or more) surface is a general function ñz=nz(x,y), so that normal and tangential components of the pressure tensor will depend on surface position.

FIG. 4.

Pressure tensor definition in systems having different geometries. (a) Cartesian coordinates for a planar geometry. (b) Spherical coordinates for a spherical geometry. (c) Cylindrical coordinates for a cylindrical geometry. (d) An example of a more general coordinate system in terms of vector surface normal nα (α = x, y, z) where one (or more) surface is a general function ñz=nz(x,y), so that normal and tangential components of the pressure tensor will depend on surface position.

Close modal
FIG. 5.

Pressure tensor profile across a planar gas–liquid interface for the argon–krypton mixture at T = 115.77 K, NAr/N = 0.5. These are smoothed molecular dynamics results for the truncated, shifted LJ model, taken from Ref. 80. The normal pressure PN is a constant (within statistical uncertainty) across the interface, and the tangential pressure PT by the IK and Harasima definitions are similar.

FIG. 5.

Pressure tensor profile across a planar gas–liquid interface for the argon–krypton mixture at T = 115.77 K, NAr/N = 0.5. These are smoothed molecular dynamics results for the truncated, shifted LJ model, taken from Ref. 80. The normal pressure PN is a constant (within statistical uncertainty) across the interface, and the tangential pressure PT by the IK and Harasima definitions are similar.

Close modal

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

PNz=nzkBT+12Szi,jNzij2rij1zijFijHzzizijHzjzzij,
(13)

where nz is the local number density at z-position. The surface area is Sz = LxLy. If the sampling of the pressure tensor is carried out over the entire xy-plane, Lx and Ly will be the simulation box size in the x- and y-directions, respectively; Lx and Ly 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=durij/drij, where urij is the pair potential between particles i and j separated by a scalar distance rij. zij = zjzi is the z-component of vector rij, 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 

pα22m=kBT2,
(14)

where pα is the molecule’s momentum in the α-direction with α = x, y, z78 and we use nz=i=1Nδzzi. 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 PT can be obtained by averaging over Pxx and Pyy. 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 by16,79

PIKTz=nzkBT+14Szi,jNxij2+yij2rij1zijFijHzzizij×Hzjzzij.
(15)

Other equally valid contour definitions are possible. For example, assuming the Harasima (H) definition leads to16,67

PHTz=nzkBT+14Szi,jNxij2+yij2rijFijδzzi,
(16)

where the Dirac delta function can be approximated as

δzzi=limΔz01ΔzΛzzi=limΔz01ΔzHziz+Δz2HzizΔz2,
(17)

where Λzzi is the so-called boxcar or top-hat function of zi for the interval from z − Δz/2 to z + Δz/2, which checks if zi 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,θ,ϕ [Fig. 4(b)]. The local spherical pressure tensor can be written as

PR=PRRRR̂R̂+PθθRθ̂θ̂+PϕϕRϕ̂ϕ̂,
(18)

where PRR = PN is the normal pressure in the radial direction; Pθθ and Pϕϕ are the equivalent tangential components (PT) due to symmetry in the polar and azimuthal directions, respectively; and R̂, θ̂, and ϕ̂ are unit vectors that are orthogonal to each other. Mechanical equilibrium [Eq. (11)] dictates the following:

PTR=PNR+R2dPNRdR.
(19)

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),

PIKRRR=nRkBT+18πR2i,jNk=12rijR̂λkrijFijHλkH1λk,
(20)

where

R̂λk=xi+λkxij/Ryi+λkyij/Rzi+λkzij/R
(21)

and λk are the roots of a quadratic equation, rij2λ2+2λririj+ri2R2=0. The polar pressure is given by15 

PIKθθR=nRkBT+18πR2i,jNk=12rijθ̂λk2rijR̂λkFijrijHλkH1λk,
(22)

where

θ̂λk=xi+λkxijzi+λkzij/Rdxyyi+λkyijzi+λkzij/Rdxydxy/R
(23)

and dxy=xi+λkxij2+yi+λkyij2. The azimuthal component is given by15 

PIKϕϕR=nRkBT+18πR2i,jNk=12rijϕ̂λk2rijR̂λkFijrijHλkH1λk,
(24)

where

ϕ̂λk=yi+λkyij/dxyxi+λkxij/dxy0.
(25)

In practice, to enhance the statistics, the tangential pressure is calculated as the average of the polar and azimuthal components, PT=Pθθ+Pϕϕ/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)]

PR=PRRRR̂R̂+PϕϕRϕ̂ϕ̂+PzzRẑẑ,
(26)

where PRR = PN is the normal pressure in the radial direction (R̂); Pϕϕ is the tangential pressure in the azimuthal direction (ϕ̂); and Pzz is the tangential pressure in the axial direction (ẑ) with PzzPϕϕ. Mechanical equilibrium [Eq. (11)] dictates that

PϕϕR=PNR+RdPNRdR.
(27)

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 by21,74

PIKRRR=nRkBT+14πRi,jNk=12rijR̂λkrijLFijHλkH1λk,
(28)

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

R̂λk=xi+λkxij/Ryi+λkyij/R0
(29)

and λk are the roots of the equation R2=xi+λxij2+yi+λyij2. The azimuthal component is written as74 

PIKϕϕR=nRkBT+14πRi,jNk=12rijϕ̂λk2rijR̂λkFijrijLHλkH1λk,
(30)

where

ϕ̂λk=yi+λkyij/Rxi+λkxij/R0.
(31)

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

PIKzzR=nRkBT+14πRi,jNk=12zij2rijR̂λkFijrijLHλkH1λk.
(32)

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 volume66 [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.

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 dS.” 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),

P=AVN,T,
(33)

where A = −kBT 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 as82 

Q=1Λ3NN!expβUrNdrN,
(34)

where rNr1, r2, …, rN 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/kBT. Following the pioneering work of Eppenga and Frenkel92 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) becomes93 

PVPAV+ΔVAVΔV=kBTΔVln1+ΔVVNexpβΔUV,
(35)

where ΔU=UV+ΔVUV 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 

Pαα=AVN,T,Lβα.
(36)

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ββα 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

PVPααz=nzkBTUzVzN,T,LβαnzkBT+kBTΔVzlnexpβΔUzV,
(37)

where Uz is the total configurational energy in an infinitesimally thin slab at a z-position and Vz=LxLyΔ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 ΔVz=SαΔLαΔz/Lz=SαζLαΔ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 ζ. Δ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.1Figure 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

Uz=12i,jNurijHziz+Δz2HzizΔz2.
(38)

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 μrij into slabs between two interacting particles i and j,34 

Uz=12i,jNΔzzijμrijHzzizijHzjzzij.
(39)

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 coordinates73,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

FIG. 6.

A schematic diagram showing the nonuniqueness involved in the distribution of potential energy in space for a slit-shaped pore system. A: Half of the pair potential contributes to the thin slab containing particle i and the other half contributes to the slab containing particle j. B: Pair potential is equally distributed in the space between two interacting particles. This illustration reflects the nonuniqueness of the local pressure tensor from a thermodynamic perspective.

FIG. 6.

A schematic diagram showing the nonuniqueness involved in the distribution of potential energy in space for a slit-shaped pore system. A: Half of the pair potential contributes to the thin slab containing particle i and the other half contributes to the slab containing particle j. B: Pair potential is equally distributed in the space between two interacting particles. This illustration reflects the nonuniqueness of the local pressure tensor from a thermodynamic perspective.

Close modal

Extension of the thermodynamic route to curved interfaces, such as spherical76 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 R0, R0=1+ζR0, where ζ is a positive, infinitesimal number. The volume change of the entire system is

ΔV=VV=3ζV,
(40)

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,

ΔA=PRR+Pθθ+Pϕϕ0R04πr2ζdr =−PRR+Pθθ+PϕϕζV.
(41)

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

P=14πR2ARN,T =PN+2PT3,
(42)

where PRR = PN is the normal pressure and Pθθ = Pϕϕ = PT 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.

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,

ρutTimeEvolution+ρuuConvection=P,
(43)

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.

Let BrN,pN be some function of the 6N phase space variables. The ensemble average ⟨B⟩ is then the 6N-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=i=1Nmiṙiδrrit, the time evolution of a phase space averaged momentum is

ti=1Nmiṙiδrrit=PK+PC,
(44)

where ṙi 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 6N dimensional integral of the ensemble average). The work of Noll109 integrates the Dirac delta function over phase space and uses Noll’s lemma to address the Oij 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 

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=i=1Nmiṙiδrrit is as before. Applying the time derivative of momentum results, after some manipulation,41,42 in the following:

ti=1Nmiṙiδrrit=i=1Nmiṙiṙiδrrit+i=1Nmir̈iδrrit=PK+PC.
(45)

Here, Newton’s law, Fi=mir̈i, 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 function78,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” form115 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 distributions117 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.

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,

u(r,t)=ρ(r,t)u(r,t)ρ(r,t)=i=1Nmiṙiδrrii=1Nmiδrri.
(46)

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 pi/mi, the particle motion not contributing to the net velocity field,

pimi=ṙiur,t.
(47)

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,

i=1Nmiṙiṙiδrri=i=1Npipimiδrri+ρ(r,t)u(r,t)u(r,t).
(48)

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 flow119 or parabolic in pressure-driven Poiseuille flow u(y) ∼ y2.120 When the expected velocity profile is known, such as in Couette or Poiseuille flow, it makes the definition of peculiar velocity, pi/mi, 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 instability131,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.

FIG. 7.

Turbulent Couette flow in a molecular simulation. (a) The evolution of iso-surfaces of turbulent kinetic energy (TKE) within a regeneration cycle with positive TKE in blue and negative in orange. (b) The average pressure moving across the channel from channel center to the top wall (y = [−1, 1] in the wall-normal direction), showing the contribution of kinetic pressure, PxyK, and configurational pressure, PxyC, which add to give total pressure, Pxy, and with the turbulent shear stress ρuv̄ give a constant (black line) at all points in the channel (satisfying mechanical equilibrium). Note that all quantities are normalized by τ0, the shear stress between the wall and first layer of the fluid.

FIG. 7.

Turbulent Couette flow in a molecular simulation. (a) The evolution of iso-surfaces of turbulent kinetic energy (TKE) within a regeneration cycle with positive TKE in blue and negative in orange. (b) The average pressure moving across the channel from channel center to the top wall (y = [−1, 1] in the wall-normal direction), showing the contribution of kinetic pressure, PxyK, and configurational pressure, PxyC, which add to give total pressure, Pxy, and with the turbulent shear stress ρuv̄ give a constant (black line) at all points in the channel (satisfying mechanical equilibrium). Note that all quantities are normalized by τ0, the shear stress between the wall and first layer of the fluid.

Close modal

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 flow136,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 × 106 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, to define a long-time-average streaming velocity ū. This average velocity allows us to define a turbulent fluctuating part, denoted here by the prime u′. This u′ 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 pi. The typical cycle of fluctuations is shown in Fig. 7(a) as the iso-surfaces of squared velocity |u′|2 = u2 + v2 + w2 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,

ρuū=ρuuρuū.
(49)

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., [ρuv̄+PxyK+PxyC]/τ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 P+ρuū=0, and as a result dū/dt=0. This turbulent equilibrium is a property of the relatively simple channel flow and would not be true in general.

FIG. 8.

The volume average (VA) form of pressure tensor, showing (a) the molecules 2 and 4 are inside the local volume and contribute to the kinetic pressure; (b) the length of line lij in the local volume used for the configurational pressure (contributions are shown in red); and (c) a more general implementation using Voronoi volumes, adapted from H. W. Hatch and P. G. Debenedetti, J. Chem. Phys. 137(3), 035103 (2012) with the permission of AIP Publishing. In both (b) and (c), the IK contour is used.

FIG. 8.

The volume average (VA) form of pressure tensor, showing (a) the molecules 2 and 4 are inside the local volume and contribute to the kinetic pressure; (b) the length of line lij in the local volume used for the configurational pressure (contributions are shown in red); and (c) a more general implementation using Voronoi volumes, adapted from H. W. Hatch and P. G. Debenedetti, J. Chem. Phys. 137(3), 035103 (2012) with the permission of AIP Publishing. In both (b) and (c), the IK contour is used.

Close modal

The off-diagonal components of the Reynolds stress tensor in Eq. (49), e.g., ρuv̄ 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,

i=1Nmiṙiṙiδrri=i=1NpipimiδrritPK+ρuū+ρuū.
(50)

Here, we have included a time averaging, ⟨⋯⟩t, over a shorter time scale tMD. The overbar denotes an average over the length of the simulation, in this case tTotal. 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 PK and the Reynolds stress will decrease, with the limit tMDtTotal 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

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,

V[PKr,t+PCr,t]dV=i=1NpipimiVδrridV+12i,jNFijrijV01δrriλrijdλdV=i=1Npipimiϑi+12i,jNFijrij01ϑλdλ,
(51)

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 ϑi=ΛxxiΛyyiΛzzi, where Λαα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, VPdV=PVAΔV, Eq. (51) results in the so-called volume average (VA) pressure tensor,

PVA=1ΔVi=1Npipimiϑi+12i,jNFijrijlij,
(52)

where ΔV is the local volume, and the shorthand lij=01ϑλdλ 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)=i=1Nmiṙiϑi/i=1Nmiϑi. 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 lij can be obtained exactly in a cube from plane–line intersections shown in Fig. 8(b) or in a sphere87,144 or cylinder21 from surface–line intersections. For more complicated local volume shapes, lij 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 Debenedetti8 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 Vg 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 

FIG. 9.

The local surface form of pressure, showing (a) the molecules 1 is entering the volume over a surface with x-normal and molecules 4 is leaving the y-normal surface, with the momentum carried contributing to kinetic pressure. (b) the molecules interacting over the x-normal surface as red squares and y-normal surface as yellow squares contribute to the configurational pressure Px and Py, respectively. (c) A more general volume with an arbitrary surface requiring a ray-tracing style solution to identify crossings with additional terms for curvature and surface movement.60 In both (b) and (c), the IK contour is used.

FIG. 9.

The local surface form of pressure, showing (a) the molecules 1 is entering the volume over a surface with x-normal and molecules 4 is leaving the y-normal surface, with the momentum carried contributing to kinetic pressure. (b) the molecules interacting over the x-normal surface as red squares and y-normal surface as yellow squares contribute to the configurational pressure Px and Py, respectively. (c) A more general volume with an arbitrary surface requiring a ray-tracing style solution to identify crossings with additional terms for curvature and surface movement.60 In both (b) and (c), the IK contour is used.

Close modal

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 at145 

limΔx0PxVAr,t=limΔx01ΔxΔyΔzi=1Npipixmiϑi+12i,jNFijxijlij=1ΔSxi=1NpipixmiδxxiΛyyiΛzzi+12i,jNFijxij01δxxiλxijΛyyλΛzzλdλ,
(53)

where ΔSx = Δ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ΔyLyΛy=limΔzLzΛz=1, ΔSx = Sx, and Eq. (53) simplifies to the method of planes (MoP) pressure,61 

PMopxx=1Sxi=1Npipixmiδxxit+14i,jNFijsgnxxisgnxxj.
(54)

The integral along λ has been evaluated to give signum functions (for a contour, use xixjδxxf(x)dx=f(x)HxxiHxxj), where sgnx=1 for x > 0 and sgnx=1 for x < 0, and sgnx=0 otherwise. The interpretation of the differences in signum functions is that xi and xj 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 Tsai65 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 yz 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 t1 to t2 and using a change of variables dt=dxi/ẋi to write

tt+ΔtPMoPxKxdt=1Sxi=1Nt1t2pipixmiδxxitdt=1Sxi=1Nxit1xit2pipixmiẋiδxxitdxi=1Sxi=1NpiHxxit1Hxxit2,PMoPxKx=12ΔtSxi=1Npisgnxxi(t1)sgnxxi(t2).
(55)

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 tt+ΔtPMoPxKdtΔtPMoPxK. Finally, we note that, to be consistent, the momentum flux in these surface definitions pi is expressed relative to a streaming velocity measured over the same surface, so if dSix=sgnxxi(t1)sgnxxi(t2) then u(x,t)=i=1NmiṙidSix/i=1NmidSix.

FIG. 10.

Conservation of momentum for a control volume in a molecular dynamics simulation, where the sum of the molecules crossing the surface (advection) and forces acting over the surface (forcing) is equal to momentum change inside the volume (accumulation).

FIG. 10.

Conservation of momentum for a control volume in a molecular dynamics simulation, where the sum of the molecules crossing the surface (advection) and forces acting over the surface (forcing) is equal to momentum change inside the volume (accumulation).

Close modal

The original MoP formulation61 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 Lee147 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., Pyx, Pyy, and Pyz as

PLPyK=12ΔtΔSyi=1Npisgnyyit1sgnyyit2Λx(xi)Λz(zi)dSiy,PLPyC=14ΔSyi,jNFijsgnyyisgnyyjΛx(xk)Λz(zk)dSijy,
(56)

where Δt = t2t1 and the function denoted as dSiy checks if molecule i is crossing the plane [Fig. 9(a)]. Function dSijy 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ΔSySyPdSy. 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(xk) in Eq. (56) identifies if a crossing xk = xi + λkxij 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,

V[PKr,t+PCr,t]dV=PLPx+PLPxΔSx+PLPy+PLPyΔSy+PLPz+PLPzΔSz=S[PKr,t+PCr,t]dS=ddtVρudV,
(57)

where PLPα+ and PLPα (α = 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 literature144 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 quantities148 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.

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 Λx=Hxix+Δx/2+ξ(y,z,t)HxixΔx/2+ξ(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 

ρuux+PSFxK+=1ΔtΔSxi=1Nmiṙit1t2×ẋi+ẏiξi+yi+żiξi+ziKineticCurvature+ξi+tSurfaceEvolutiondSix+dt,PSFxC+=12ΔSxi,jNFij01xij+yijξλ+yλ+zijξλ+zλConfigurationalCurvaturedSλx+dλ.
(58)

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 ξi+ (or ξλ+). The dSix+anddSλ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 ñx=αξxα/αξxα 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

ρuux+PSFxK+=1ΔtΔSxi=1Nmiṙiri12ñxri12ñxdS++1ΔSxi=1Nmiṙiϑt,PSFxC+=12ΔSxi,jNFijrijñxrijñxdS+.
(59)

The surface evolution term is contained in the function ϑt=Hξt2xit2Hξt1xit2Λy(yit2)Λ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 t1 and t2. For the kinetic term, ri12 = ri2ri1 is the line of time evolution of a molecule i between t1 and t2 mirroring the configurational term’s intermolecular (IK) contour rij = rjri. Assuming the IK contour, the equation for a line is rλ = ri1 + λri12 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:

dS+=k=1NrootsH1λkHλkΛyykΛzzk,
(60)

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, yk and zk, 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 ñx in Eq. (59) behaves like the radial (R̂) or azimuthal (ϕ̂) 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., ñx=ñxy,z. The function H1λkHλk in Eq. (60) collects only the interactions crossing the surface as the function HλkH1λ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.

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.

FIG. 11.

An example of the use of mapping in obtaining a pressure relevant to the geometry of a problem. The example is a single snapshot in time of a NEMD boiling simulation, where a bubble starts from a square notch on a wall and expands into a circular bubble shown here. The bubble is clear in the density field of (a) with blue vapor density and yellow liquid density. The identified interface points are shown as black pixels and the green line shows a fitted circle arc used to get radius and angle. (b) The pressure field of Pxx component with blue indicating large negative values. (c) The pressure rotated using the angle from the fitted circle and Eq. (61) to give tangential pressure Pθθ, with green lines shown at θb = −90°, −45°, 0°, 45°, and θt = 90° to guide the eye. The arrow shows the whole 180° arc. (d) The mapped field using a Cartesian to polar mapping, where the green lines correspond to the ones from (c) between θb and θt. (e) The average normal and tangential pressure over the 180 arc.

FIG. 11.

An example of the use of mapping in obtaining a pressure relevant to the geometry of a problem. The example is a single snapshot in time of a NEMD boiling simulation, where a bubble starts from a square notch on a wall and expands into a circular bubble shown here. The bubble is clear in the density field of (a) with blue vapor density and yellow liquid density. The identified interface points are shown as black pixels and the green line shows a fitted circle arc used to get radius and angle. (b) The pressure field of Pxx component with blue indicating large negative values. (c) The pressure rotated using the angle from the fitted circle and Eq. (61) to give tangential pressure Pθθ, with green lines shown at θb = −90°, −45°, 0°, 45°, and θt = 90° to guide the eye. The arrow shows the whole 180° arc. (d) The mapped field using a Cartesian to polar mapping, where the green lines correspond to the ones from (c) between θb and θt. (e) The average normal and tangential pressure over the 180 arc.

Close modal

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 Pxx in Fig. 11(b),153 

PRR=Pxxcos2θ+Pyysin2θ+Pxysin2θ,Pθθ=Pxxsin2θ+Pyycos2θPxysin2θ.
(61)

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 PRR, 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.

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 formula43,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.

FIG. 12.

Comparison of the probability distribution function of Pxx using the VA form (ij is length of line in a volume shown schematically in red on the top right) and LP form (dS̲ij is surface crossing shown schematically as red crosses with surface normal on the bottom right). For a reduced cell side length of 13.7 σ, the blue solid line on the main plot is the VA pressure and the dotted yellow line is the LP pressure. For a reduced cell side length of 6.8 σ, the green solid line is the corresponding VA pressure and the dotted red line is the LP pressure. Data are taken from Ref. 157.

FIG. 12.

Comparison of the probability distribution function of Pxx using the VA form (ij is length of line in a volume shown schematically in red on the top right) and LP form (dS̲ij is surface crossing shown schematically as red crosses with surface normal on the bottom right). For a reduced cell side length of 13.7 σ, the blue solid line on the main plot is the VA pressure and the dotted yellow line is the LP pressure. For a reduced cell side length of 6.8 σ, the green solid line is the corresponding VA pressure and the dotted red line is the LP pressure. Data are taken from Ref. 157.

Close modal

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 ⟨Pxx⟩ ≈ 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.

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 Clausius81 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 articles159,160 with Gray160 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 Buff70 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 constraints78,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 Tadmor78 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 Sun171 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 kth bin (k = 1, 2, 3, …) along the z-axis is given by

PCGk=1ΔrkΔrkPrdr=1ΔzkΔzkPzdz,
(62)

where the local (averaging) volume of the kth bin is Δrk = LxLyΔzk, Lx and Ly are the constant lateral dimension of the pore surface in the x- and y-directions, respectively, and Δzk is the characteristic length (width) of the kth bin that leads to a unique CG pressure tensor. This CG pressure tensor has the same appearance as the conventional VA definitions110,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λij, where λij is the linearly scaled z-distance from particle i to particle j (assuming zi < zj); thus λij = 0 amounts to z = zi and λij = 1 is z = zj. 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λ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 Δz1, Δz2 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.

FIG. 13.

Simulation evidence for the uniqueness of the CG pressure tensor defined in Eq. (62).35 (a) Graphs showing ten functional forms for fCλij with λij=zminzi,zj/zij, corresponding to ten unique contour definitions. (b) Reduced number density profile (upper panel) for LJ argon adsorbed in a structureless carbon slit pore and the integral of the (configurational) local tangential pressure PxxCz over z-direction using the ten contour definitions (bottom panel). Only the adsorbate–adsorbate interactions contribute to the tangential pressure calculations. In this case, the characteristic lengths Δzk that can lead to a unique CG pressure tensor were chosen to be those that are comparable to the thickness of the adsorbed layer, and they are marked in the plot as Δz1, Δz2, etc. Adapted from Shi et al., J. Chem. Phys. 154(8), 084502 (2021) with the permission of AIP Publishing.

FIG. 13.

Simulation evidence for the uniqueness of the CG pressure tensor defined in Eq. (62).35 (a) Graphs showing ten functional forms for fCλij with λij=zminzi,zj/zij, corresponding to ten unique contour definitions. (b) Reduced number density profile (upper panel) for LJ argon adsorbed in a structureless carbon slit pore and the integral of the (configurational) local tangential pressure PxxCz over z-direction using the ten contour definitions (bottom panel). Only the adsorbate–adsorbate interactions contribute to the tangential pressure calculations. In this case, the characteristic lengths Δzk that can lead to a unique CG pressure tensor were chosen to be those that are comparable to the thickness of the adsorbed layer, and they are marked in the plot as Δz1, Δz2, etc. Adapted from Shi et al., J. Chem. Phys. 154(8), 084502 (2021) with the permission of AIP Publishing.

Close modal

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.

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 algorithm180). The molecular pressure tensor takes the molecule to be a rigid body (i.e., rigid body approximation181), 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 constraints74,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.

FIG. 14.

Central force (CF) decomposition for a three-body potential. In a CF decomposition,78 the total force on an atom Fi is decomposed into pairwise central terms Fij; the (unknown) magnitude of Fij can be obtained by solving a system of linear equations in Eq. (63). In a non-CF decomposition, such as the scheme proposed by Goetz and Lipowsky,4 the pairwise term is simply FjFi/m with m = 3 for a three-body potential; this pairwise term is noncentral, i.e., pair force is not parallel to the vector rij.

FIG. 14.

Central force (CF) decomposition for a three-body potential. In a CF decomposition,78 the total force on an atom Fi is decomposed into pairwise central terms Fij; the (unknown) magnitude of Fij can be obtained by solving a system of linear equations in Eq. (63). In a non-CF decomposition, such as the scheme proposed by Goetz and Lipowsky,4 the pairwise term is simply FjFi/m with m = 3 for a three-body potential; this pairwise term is noncentral, i.e., pair force is not parallel to the vector rij.

Close modal

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 rij),

Fi=riUmrm =jiFijrijrij,
(63)

where rmr1, r2, …, rm 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 linear175,189,191 or angular4 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 Fij are obtained by solving a system of linear equations given in Eq. (63). Here, the number of independent force equations is 3m − 6 (as forces satisfy the conservation of linear and angular momentum), and the number of unknown pairwise central terms is mm1/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 method193 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 interface6,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 Debenedetti8 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-potential197 or a shifted-force198 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.

TABLE I.

Pressure or stress tensor functionalities available in the LAMMPS software. Note that the stress tensor is defined as the negative of the pressure tensor. More detailed explanations and restrictions of these commands are available in the LAMMPS manual (https://docs.lammps.org/Manual.html, accessed on May 16, 2022, LAMMPS version 4 May 2022).

LAMMPS commandExplanation
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-range178 and many-body interactions.206  
Compute pressure Computation of a scalar pressure and a global pressure tensor of the entire system206  
Compute pressure/uef Computation of the pressure tensor in the reference frame of the applied flow field 
Compute stress/mopa 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/Cartesiana,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/sphericala,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/cylindera,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 commandExplanation
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-range178 and many-body interactions.206  
Compute pressure Computation of a scalar pressure and a global pressure tensor of the entire system206  
Compute pressure/uef Computation of the pressure tensor in the reference frame of the applied flow field 
Compute stress/mopa 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/Cartesiana,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/sphericala,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/cylindera,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.

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 

FIG. 15.

Molecular simulation and experimental results for the normal pressure of (a) CCl4 and (b) H2O adsorbed in slit carbon pores of different pore width H, at 300 K and 1 bar bulk pressure. The experimental data were estimated using Young’s equation based on the change of interplanar distance of the activated carbon fibers. Adapted with permission from Long et al., Colloids Surf., A 437, 33–41 (2013). Copyright 2013 Elsevier.

FIG. 15.

Molecular simulation and experimental results for the normal pressure of (a) CCl4 and (b) H2O adsorbed in slit carbon pores of different pore width H, at 300 K and 1 bar bulk pressure. The experimental data were estimated using Young’s equation based on the change of interplanar distance of the activated carbon fibers. Adapted with permission from Long et al., Colloids Surf., A 437, 33–41 (2013). Copyright 2013 Elsevier.

Close modal

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.

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, PN, 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 CCl4 and H2O adsorbed in carbon slit pores.

For the tangential pressure, PT, 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 P2D (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 leff in the direction normal to the surface. The effective tangential pressure estimated by the 2D route is210 

P2DTP2DT,ρ2Dleff,
(64)

where T is temperature and ρ2D is the 2D density of the adsorbed film. The 2D pressure P2D is a function of T and ρ2D, and the relation can be established by a 2D equation of state.209 Although P2D is well defined, P2DT is not unique due to the arbitrary choice of leff. To decide on a sensible choice of leff, it is instructive to rewrite Eqs. (64) as a spatial (volume) average,210 i.e., P2DTleffPTzdz/leff. Based on the results in Fig. 13, it is motivating to choose leff to be the characteristic length Δzk 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 leff. 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 leff 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.

FIG. 16.

A schematic knowledge flowchart relating different configurational pressure forms discussed in this Perspective. See Fig. 17 for a similar flowchart for the kinetic pressure. Any form without angular brackets can be applied away from equilibrium.

FIG. 16.

A schematic knowledge flowchart relating different configurational pressure forms discussed in this Perspective. See Fig. 17 for a similar flowchart for the kinetic pressure. Any form without angular brackets can be applied away from equilibrium.

Close modal
FIG. 17.

A schematic knowledge flowchart relating different kinetic pressure forms discussed in this Perspective. See Fig. 16 for a similar flowchart for the configurational pressure.

FIG. 17.

A schematic knowledge flowchart relating different kinetic pressure forms discussed in this Perspective. See Fig. 16 for a similar flowchart for the configurational pressure.

Close modal

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=βT1, where the isothermal compressibility βT is defined as
    βT=1VVPT.
    (65)
    β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.

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.

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.

The authors have no conflicts to disclose.

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).

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

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

PIKNCz=12Szi,jNzij2rij1zijFijHzzizijHzjzzij,
(A1)
PIKTxCz=12Szi,jNxij2rij1zijFijHzzizijHzjzzij,
(A2)
PIKTyCz=12Szi,jNyij2rij1zijFijHzzizijHzjzzij.
(A3)

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

PMoPzCz=14Szi,jNFijsgnzzisgnzzj.
(A4)

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

PIKzCz=12Szi,jNFijrij1zijHzzizijHzjzzij,
(A5)

where Fijrij denotes the element-wise product of force vector Fij and distance vector rij. The Heaviside function expressed as a product gives identical behavior to the difference of two Heaviside functions,

HzzizijHzjzzij=HzzizijHzzjzij.
(A6)

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

HzzizijHzzjzij=12sgn1zijsgnzzi+112sgn1zijsgnzzj+1=12sgn1zijsgn(zzi)sgn(zzj).
(A7)

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

PIKNC=PIKzzC=12Szi,jNzijzijFijz12sgn1zijsgnzzisgnzzj=14Szi,jNFijzsgn(zzi)sgn(zzj)=PMoPzzC.
(A8)

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

PIKTyC=PIKyyC=14Szi,jNyijzijFijysgn(zzi)sgn(zzj),
(A9)

where we have an extra factor of yij/zij 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 yij 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 Fijyyij to get Pyy, but in the VA form, this is weighted by the fraction of interaction lij 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 zij 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.

FIG. 18.

Schematic showing Pyy and Pzz on a single z-normal plane from the IK contour method, while the MoP pressure defines a y-normal plane in order to get Pyy.

FIG. 18.

Schematic showing Pyy and Pzz on a single z-normal plane from the IK contour method, while the MoP pressure defines a y-normal plane in order to get Pyy.

Close modal

2. The Noll form of pressure

The Noll reformulation replaces the Dirac delta function δrrif|ri=r. This notation denotes an integral of probability density function f over all phase space except ri, 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

PNOLLCr,t=12i,jNzz01Fijf|ri=r+λz,rj=r(1λ)zdλdz.
(A10)

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 ri to rj 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 = rij = rjri, we can rearrange both equalities in the brackets to show they are the same condition rj − (1 − λ)rij = ri + λrij,

PNOLLCr,t=12i,jNrijrij01Fijf|ri+λrij=rdλdrij,
(A11)

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 average41,110,112 as discussed in Sec. IV.

1.
J. H.
Irving
and
J. G.
Kirkwood
, “
The statistical mechanical theory of transport processes. IV. The equations of hydrodynamics
,”
J. Chem. Phys.
18
(
6
),
817
829
(
1950
).
2.
P.
Schofield
and
J. R.
Henderson
, “
Statistical mechanics of inhomogeneous fluids
,”
Proc. R. Soc. A
379
(
1776
),
231
246
(
1982
).
3.
J. S.
Rowlinson
and
B.
Widom
,
Molecular Theory of Capillarity
(
Clarendon Press
,
Oxford
,
1982
).
4.
R.
Goetz
and
R.
Lipowsky
, “
Computer simulations of bilayer membranes: Self-assembly and interfacial tension
,”
J. Chem. Phys.
108
(
17
),
7397
7409
(
1998
).
5.
R. H.
Templer
,
S. J.
Castle
,
A.
Rachael Curran
,
G.
Rumbles
, and
D. R.
Klug
, “
Sensing isothermal changes in the lateral pressure in model membranes using di-pyrenyl phosphatidylcholine
,”
Faraday Discuss.
111
,
41
53
(
1999
).
6.
J.
Sonne
,
F. Y.
Hansen
, and
G. H.
Peters
, “
Methodological problems in pressure profile calculations for lipid bilayers
,”
J. Chem. Phys.
122
(
12
),
124903
(
2005
).
7.
J. M.
Vanegas
,
A.
Torres-Sánchez
, and
M.
Arroyo
, “
Importance of force decomposition for local stress calculations in biomembrane molecular simulations
,”
J. Chem. Theory Comput.
10
(
2
),
691
702
(
2014
).
8.
H. W.
Hatch
and
P. G.
Debenedetti
, “
Molecular modeling of mechanical stresses on proteins in glassy matrices: Formalism
,”
J. Chem. Phys.
137
(
3
),
035103
(
2012
).
9.
A. L.
Kolesnikov
,
H.
Uhlig
,
J.
Möllmer
,
J.
Adolphs
,
Y. A.
Budkov
,
N.
Georgi
,
D.
Enke
, and
R.
Gläser
, “
Pore size distribution of MCM-41-type silica materials from pseudomorphic transformation—A minimal input data approach based on excess surface work
,”
Microporous Mesoporous Mater.
240
,
169
177
(
2017
).
10.
C.
Balzer
,
A. M.
Waag
,
S.
Gehret
,
G.
Reichenauer
,
F.
Putz
,
N.
Hüsing
,
O.
Paris
,
N.
Bernstein
,
G. Y.
Gor
, and
A. V.
Neimark
, “
Adsorption-induced deformation of hierarchically structured mesoporous silica—Effect of pore-level anisotropy
,”
Langmuir
33
(
22
),
5592
5602
(
2017
).
11.
S. M.
Thompson
,
K. E.
Gubbins
,
J. P. R. B.
Walton
,
R. A. R.
Chantry
, and
J. S.
Rowlinson
, “
A molecular dynamics study of liquid drops
,”
J. Chem. Phys.
81
(
1
),
530
542
(
1984
).
12.
P. R.
ten Wolde
and
D.
Frenkel
, “
Computer simulation study of gas–liquid nucleation in a Lennard-Jones system
,”
J. Chem. Phys.
109
(
22
),
9901
9918
(
1998
).
13.
A.
Malijevský
and
G.
Jackson
, “
A perspective on the interfacial properties of nanoscopic liquid drops
,”
J. Phys.: Condens. Matter
24
(
46
),
464121
(
2012
).
14.
K. G. S. H.
Gunawardana
and
X.
Song
, “
Theoretical prediction of crystallization kinetics of a supercooled Lennard-Jones fluid
,”
J. Chem. Phys.
148
(
20
),
204506
(
2018
).
15.
P.
Montero de Hijes
,
K.
Shi
,
E. G.
Noya
,
E. E.
Santiso
,
K. E.
Gubbins
,
E.
Sanz
, and
C.
Vega
, “
The Young–Laplace equation for a solid–liquid interface
,”
J. Chem. Phys.
153
(
19
),
191102
(
2020
).
16.
J. P. R. B.
Walton
,
D. J.
Tildesley
,
J. S.
Rowlinson
, and
J. R.
Henderson
, “
The pressure tensor at the planar surface of a liquid
,”
Mol. Phys.
48
(
6
),
1357
1368
(
1983
).
17.
M.
Sega
,
B.
Fábián
, and
P.
Jedlovszky
, “
Layer-by-layer and intrinsic analysis of molecular and thermodynamic properties across soft interfaces
,”
J. Chem. Phys.
143
(
11
),
114709
(
2015
).
18.
C.
Braga
,
E. R.
Smith
,
A.
Nold
,
D. N.
Sibley
, and
S.
Kalliadasis
, “
The pressure tensor across a liquid-vapour interface
,”
J. Chem. Phys.
149
(
4
),
044705
(
2018
).
19.
M.
Orselly
,
J.
Devémy
,
A.
Bouvet-Marchand
,
A.
Dequidt
,
C.
Loubat
, and
P.
Malfreyt
, “
Molecular interactions at the metal–liquid interfaces
,”
J. Chem. Phys.
156
(
23
),
234705
(
2022
).
20.
J.
Vrabec
,
G. K.
Kedia
,
G.
Fuchs
, and
H.
Hasse
, “
Comprehensive study of the vapour–liquid coexistence of the truncated and shifted Lennard–Jones fluid including planar and spherical interface properties
,”
Mol. Phys.
104
(
9
),
1509
1527
(
2006
).
21.
C. K.
Addington
,
Y.
Long
, and
K. E.
Gubbins
, “
The pressure in interfaces having cylindrical geometry
,”
J. Chem. Phys.
149
(
8
),
084109
(
2018
).
22.
K.
Koga
,
G. T.
Gao
,
H.
Tanaka
, and
X. C.
Zeng
, “
Formation of ordered ice nanotubes inside carbon nanotubes
,”
Nature
412
(
6849
),
802
805
(
2001
).
23.
V. V.
Chaban
,
V. V.
Prezhdo
, and
O. V.
Prezhdo
, “
Confinement by carbon nanotubes drastically alters the boiling and critical behavior of water droplets
,”
ACS Nano
6
(
3
),
2766
2773
(
2012
).
24.
K.
Mochizuki
and
K.
Koga
, “
Solid–liquid critical behavior of water in nanopores
,”
Proc. Natl. Acad. Sci. U. S. A.
112
(
27
),
8221
8226
(
2015
).
25.
M.
Sadeghi
and
G. A.
Parsafar
, “
Toward an equation of state for water inside carbon nanotubes
,”
J. Phys. Chem. B
116
(
16
),
4943
4951
(
2012
).
26.
V.
Bråten
,
D. T.
Zhang
,
M.
Hammer
,
A.
Aasen
,
S. K.
Schnell
, and
Ø.
Wilhelmsen
, “
Equation of state for confined fluids
,”
J. Chem. Phys.
156
(
24
),
244504
(
2022
).
27.
K. E.
Gubbins
,
K.
Gu
,
L.
Huang
,
Y.
Long
,
J. M.
Mansell
,
E. E.
Santiso
,
K.
Shi
,
M.
Śliwińska-Bartkowiak
, and
D.
Srivastava
, “
Surface-driven high-pressure processing
,”
Engineering
4
(
3
),
311
320
(
2018
).
28.
K. S.
Vasu
,
E.
Prestat
,
J.
Abraham
,
J.
Dix
,
R. J.
Kashtiban
,
J.
Beheshtian
,
J.
Sloan
,
P.
Carbone
,
M.
Neek-Amal
,
S. J.
Haigh
,
A. K.
Geim
, and
R. R.
Nair
, “
Van der Waals pressure and its effect on trapped interlayer molecules
,”
Nat. Commun.
7
,
12168
(
2016
).
29.
D.
Srivastava
,
C. H.
Turner
,
E. E.
Santiso
, and
K. E.
Gubbins
, “
The nitric oxide dimer reaction in carbon nanopores
,”
J. Phys. Chem. B
122
(
13
),
3604
3614
(
2018
).
30.
K.
Urita
,
Y.
Shiga
,
T.
Fujimori
,
T.
Iiyama
,
Y.
Hattori
,
H.
Kanoh
,
T.
Ohba
,
H.
Tanaka
,
M.
Yudasaka
,
S.
Iijima
,
I.
Moriguchi
,
F.
Okino
,
M.
Endo
, and
K.
Kaneko
, “
Confinement in carbon nanospace-induced production of KI nanocrystals of high-pressure phase
,”
J. Am. Chem. Soc.
133
(
27
),
10344
10347
(
2011
).
31.
T.
Fujimori
,
A.
Morelos-Gómez
,
Z.
Zhu
,
H.
Muramatsu
,
R.
Futamura
,
K.
Urita
,
M.
Terrones
,
T.
Hayashi
,
M.
Endo
,
S.
Young Hong
,
Y.
Chul Choi
,
D.
Tománek
, and
K.
Kaneko
, “
Conducting linear chains of sulphur inside carbon nanotubes
,”
Nat. Commun.
4
(
1
),
2162
(
2013
).
32.
C. K.
Addington
,
J. M.
Mansell
, and
K. E.
Gubbins
, “
Computer simulation of conductive linear sulfur chains confined in carbon nanotubes
,”
Mol. Simul.
43
(
7
),
519
525
(
2017
).
33.
Y.
Long
,
J. C.
Palmer
,
B.
Coasne
,
M.
Śliwinska-Bartkowiak
, and
K. E.
Gubbins
, “
Pressure enhancement in carbon nanopores: A major confinement effect
,”
Phys. Chem. Chem. Phys.
13
(
38
),
17163
(
2011
).
34.
Y.
Long
,
J. C.
Palmer
,
B.
Coasne
,
M.
Śliwinska-Bartkowiak
,
G.
Jackson
,
E. A.
Müller
, and
K. E.
Gubbins
, “
On the molecular origin of high-pressure effects in nanoconfinement: The role of surface chemistry and roughness
,”
J. Chem. Phys.
139
(
14
),
144701
(
2013
).
35.
K.
Shi
,
E. E.
Santiso
, and
K. E.
Gubbins
, “
Can we define a unique microscopic pressure in inhomogeneous fluids?
,”
J. Chem. Phys.
154
(
8
),
084502
(
2021
).
36.
S.
Brunauer
and
P. H.
Emmett
, “
The use of low temperature van der Waals adsorption isotherms in determining the surface areas of various adsorbents
,”
J. Am. Chem. Soc.
59
(
12
),
2682
2689
(
1937
).
37.
G. L.
Aranovich
and
M. D.
Donohue
, “
Surface compression in adsorption systems
,”
Colloids Surf., A
187-188
(
188
),
95
108
(
2001
).
38.
S.
Abaza
,
G. L.
Aranovich
, and
M. D.
Donohue
, “
Adsorption compression in surface layers
,”
Mol. Phys.
110
(
11–12
),
1289
1298
(
2012
).
39.
Y.
Long
,
J. C.
Palmer
,
B.
Coasne
,
K.
Shi
,
M.
Śliwińska-Bartkowiak
, and
K. E.
Gubbins
, “
Reply to the ‘Comment on ‘Pressure enhancement in carbon nanopores: A major confinement effect’’ by D. van Dijk, Phys. Chem. Chem. Phys., 2020, 22, DOI: 10.1039/C9CP02890K
,”
Phys. Chem. Chem. Phys.
22
(
17
),
9826
9830
(
2020
).
40.
L. D.
Landau
and
E. M.
Lifshitz
,
Fluid Mechanics
, 1st ed. (
Pergamon Press
,
1969
).
41.
D. J.
Evans
and
G. P.
Morriss
,
Statistical Mechanics of Nonequilibrium Liquids
, 2nd ed. (
ANU Press
,
2007
).
42.
B. D.
Todd
and
P. J.
Daivis
,
Nonequilibrium Molecular Dynamics
(
Cambridge University Press
,
2017
).
43.
M. S.
Green
, “
Markoff random processes and the statistical mechanics of time‐dependent phenomena. II. Irreversible processes in fluids
,”
J. Chem. Phys.
22
(
3
),
398
413
(
1954
).
44.
R.
Kubo
, “
Statistical-mechanical theory of irreversible processes. I. General theory and simple applications to magnetic and conduction problems
,”
J. Phys. Soc. Jpn.
12
(
6
),
570
586
(
1957
).
45.
W. T.
Ashurst
and
W. G.
Hoover
, “
Argon shear viscosity via a Lennard-Jones potential with equilibrium and nonequilibrium molecular dynamics
,”
Phys. Rev. Lett.
31
(
4
),
206
208
(
1973
).
46.
J. P.
Ewen
,
C.
Gattinoni
,
J.
Zhang
,
D. M.
Heyes
,
H. A.
Spikes
, and
D.
Dini
, “
On the effect of confined fluid molecular structure on nonequilibrium phase behaviour and friction
,”
Phys. Chem. Chem. Phys.
19
(
27
),
17883
17894
(
2017
).
47.
S.
Kumar Kannam
,
B. D.
Todd
,
J. S.
Hansen
, and
P. J.
Daivis
, “
Slip length of water on graphene: Limitations of non-equilibrium molecular dynamics simulations
,”
J. Chem. Phys.
136
(
2
),
024705
(
2012
).
48.
C. D.
Williams
,
Z.
Wei
,
M. R. b.
Shaharudin
, and
P.
Carbone
, “
A molecular simulation study into the stability of hydrated graphene nanochannels used in nanofluidics devices
,”
Nanoscale
14
(
9
),
3467
3479
(
2022
).
49.
P. A.
Thompson
and
M. O.
Robbins
, “
Simulations of contact-line motion: Slip and the dynamic contact angle
,”
Phys. Rev. Lett.
63
(
7
),
766
769
(
1989
).
50.
T.
Qian
,
X.-P.
Wang
, and
P.
Sheng
, “
Molecular hydrodynamics of the moving contact line in two-phase immiscible flows
,”
Commun. Comput. Phys.
1
,
1
52
(
2005
); arXiv:cond-mat/0510403.
51.
R.
Adkins
,
I.
Kolvin
,
Z.
You
,
S.
Witthaus
,
M. C.
Marchetti
, and
Z.
Dogic
, “
Dynamics of active liquid interfaces
,”
Science
377
(
6607
),
768
772
(
2022
).
52.
S. T.
O’Connell
and
P. A.
Thompson
, “
Molecular dynamics–continuum hybrid computations: A tool for studying complex fluid flows
,”
Phys. Rev. E
52
(
6
),
R5792
(
1995
).
53.
K. M.
Mohamed
and
A. A.
Mohamad
, “
A review of the development of hybrid atomistic–continuum methods for dense fluids
,”
Microfluid. Nanofluid.
8
(
3
),
283
302
(
2010
).
54.
E. G.
Flekkøy
,
G.
Wagner
, and
J.
Feder
, “
Hybrid model for combined particle and continuum dynamics
,”
Europhys. Lett.
52
(
3
),
271
276
(
2000
).
55.
D.
Davydov
,
J.-P.
Pelteret
, and
P.
Steinmann
, “
Comparison of several staggered atomistic-to-continuum concurrent coupling strategies
,”
Comput. Methods Appl. Mech. Eng.
277
,
260
280
(
2014
).
56.
Y.
Cui
and
H. B.
Chew
, “
Machine-learning prediction of atomistic stress along grain boundaries
,”
Acta Mater.
222
,
117387
(
2022
).
57.
S.
Xiao
,
R.
Hu
,
Z.
Li
,
S.
Attarian
,
K.-M.
Björk
, and
A.
Lendasse
, “
A machine-learning-enhanced hierarchical multiscale method for bridging from molecular dynamics to continua
,”
Neural Comput. Appl.
32
(
18
),
14359
14373
(
2020
).
58.
V.
Strobel
(
2022
). “
Academic-keyword-occurrence
,” Zenodo. , https://zenodo.org/record/1218409#.YsYeX3bMKbg; accessed July 6, 2022.
59.
E. R.
Smith
,
E. A.
Müller
,
R. V.
Craster
, and
O. K.
Matar
, “
A Langevin model for fluctuating contact angle behaviour parametrised using molecular dynamics
,”
Soft Matter
12
(
48
),
9604
9615
(
2016
).
60.
E. R.
Smith
and
C.
Braga
, “
Hydrodynamics across a fluctuating interface
,”
J. Chem. Phys.
153
(
13
),
134705
(
2020
).
61.
B. D.
Todd
,
D. J.
Evans
, and
P. J.
Daivis
, “
Pressure tensor for inhomogeneous fluids
,”
Phys. Rev. E
52
(
2
),
1627
1638
(
1995
).
62.
G. C.
Rossi
and
M.
Testa
, “
The stress tensor in thermodynamics and statistical mechanics
,”
J. Chem. Phys.
132
(
7
),
074902
(
2010
).
63.
F.
Varnik
,
J.
Baschnagel
, and
K.
Binder
, “
Molecular dynamics results on the pressure tensor of polymer films
,”
J. Chem. Phys.
113
(
10
),
4444
4453
(
2000
).
64.
J.-G.
Weng
,
S.
Park
,
J. R.
Lukes
, and
C.-L.
Tien
, “
Molecular dynamics investigation of thickness effect on liquid films
,”
J. Chem. Phys.
113
(
14
),
5917
5923
(
2000
).
65.
D. H.
Tsai
, “
The virial theorem and stress calculation in molecular dynamics
,”
J. Chem. Phys.
70
(
3
),
1375
1382
(
1979
).
66.
E. R.
Smith
, “
The importance of reference frame for pressure at the liquid–vapour interface
,”
Mol. Simul.
48
(
1
),
57
72
(
2022
).
67.
A.
Harasima
, “
Molecular theory of surface tension
,” in
Advances in Chemical Physics
(
John Wiley & Sons, Inc.
,
1958
), Vol. 64, pp.
203
237
.
68.
I.
Newton
,
Philosophiae Naturalis Principia Mathematica
, 1st ed. (
Jussu Societatis Regiae ac Typis Josephi Streater. Prostat Apud Plures Bibliopolas
,
1687
).
69.
A.
Motte
,
The Mathematical Principles of Natural Philosophy
, 1st ed. (
1729
).
70.
J. G.
Kirkwood
and
F. P.
Buff
, “
The statistical mechanical theory of surface tension
,”
J. Chem. Phys.
17
(
3
),
338
343
(
1949
).
71.
C. G.
Gray
,
K. E.
Gubbins
, and
C. G.
Joslin
, “
Pressure tensor
,” in
Theory of Molecular Fluids. Vol. 2: Applications
(
Oxford University Press
,
2011
), pp.
928
942
, Sec. 8.3.
72.
T.
Nakamura
,
S.
Kawamoto
, and
W.
Shinoda
, “
Precise calculation of the local pressure tensor in Cartesian and spherical coordinates in LAMMPS
,”
Comput. Phys. Commun.
190
,
120
128
(
2015
).
73.
B.
Hafskjold
and
T.
Ikeshoji
, “
Microscopic pressure tensor for hard-sphere fluids
,”
Phys. Rev. E
66
(
1
),
011203
(
2002
).
74.
K.
Shi
,
Y.
Shen
,
E. E.
Santiso
, and
K. E.
Gubbins
, “
Microscopic pressure tensor in cylindrical geometry: Pressure of water in a carbon nanotube
,”
J. Chem. Theory Comput.
16
(
9
),
5548
5561
(
2020
).
75.
E. N.
Parker
, “
Tensor virial equations
,”
Phys. Rev.
96
(
6
),
1686
1689
(
1954
).
76.
Y.
Long
, “
Pressure tensor of adsorbate in nanoporous materials: Molecular simulation studies
,” Ph.D. dissertation (
North Carolina State University
,
2012
).
77.
C. G.
Gray
and
K. E.
Gubbins
, “
Basic statistical mechanics
,” in
Theory of Molecular Fluids. Vol. 1: Fundamentals
(
Oxford University Press
,
1984
), pp.
143
240
, Sec. 3.
78.
N. C.
Admal
and
E. B.
Tadmor
, “
A unified interpretation of stress in molecular systems
,”
J. Elasticity
100
(
1–2
),
63
143
(
2010
).
79.
M.
Rao
and
B. J.
Berne
, “
On the location of surface of tension in the planar interface between liquid and vapour
,”
Mol. Phys.
37
(
2
),
455
461
(
1979
).
80.
D. J.
Lee
,
M. M.
Telo da Gama
, and
K. E.
Gubbins
, “
The vapour-liquid interface for a Lennard-Jones model of argon-krypton mixtures
,”
Mol. Phys.
53
(
5
),
1113
1130
(
1984
).
81.
R.
Clausius
, “
XVI. On a mechanical theorem applicable to heat
,”
London, Edinburgh, Dublin Philos. Mag. J. Sci.
40
(
265
),
122
127
(
1870
).
82.
J.-P.
Hansen
,
I. R.
McDonald
,
J.-P.
Hansen
, and
I. R.
McDonald
, “
Chapter 2—Statistical mechanics
,” in
Theory of Simple Liquids
(
Elsevier
,
2013
), pp.
13
59
.
83.
S. M. A.
Malek
,
F.
Sciortino
,
P. H.
Poole
, and
I.
Saika-Voivod
, “
Evaluating the Laplace pressure of water nanodroplets from simulations
,”
J. Phys.: Condens. Matter
30
(
14
),
144005
(
2018
).
84.
A.
Emelianova
,
M. A.
Maximov
, and
G. Y.
Gor
, “
Solvation pressure in spherical mesopores: Macroscopic theory and molecular simulations
,”
AIChE J.
67
(
3
),
e16542
(
2021
).
85.
J. I.
Cutler
,
E.
Auyeung
, and
C. A.
Mirkin
, “
Spherical nucleic acids
,”
J. Am. Chem. Soc.
134
(
3
),
1376
1391
(
2012
).
86.
E. M.
Blokhuis
and
D.
Bedeaux
, “
Pressure tensor of a spherical interface
,”
J. Chem. Phys.
97
(
5
),
3576
3586
(
1992
).
87.
G. J.
Tjatjopoulos
and
J. A.
Mann
, “
The pressure tensor of an inhomogeneous fluid
,”
Mol. Phys.
60
(
6
),
1425
1432
(
1987
).
88.
X.
Wang
,
G.
Guerin
,
H.
Wang
,
Y.
Wang
,
I.
Manners
, and
M. A.
Winnik
, “
Cylindrical block copolymer micelles and co-micelles of controlled length and architecture
,”
Science
317
(
5838
),
644
647
(
2007
).
89.
H.
Qiu
,
Z. M.
Hudson
,
M. A.
Winnik
, and
I.
Manners
, “
Multidimensional hierarchical self-assembly of amphiphilic cylindrical block comicelles
,”
Science
347
(
6228
),
1329
1332
(
2015
).
90.
A.
Statt
,
P.
Virnau
, and
K.
Binder
, “
Finite-size effects on liquid-solid phase coexistence and the estimation of crystal nucleation barriers
,”
Phys. Rev. Lett.
114
(
2
),
026101
(
2015
).
91.
P.
Montero de Hijes
,
J. R.
Espinosa
,
V.
Bianco
,
E.
Sanz
, and
C.
Vega
, “
Interfacial free energy and Tolman length of curved liquid–solid interfaces from equilibrium studies
,”
J. Phys. Chem. C
124
(
16
),
8795
8805
(
2020
).
92.
R.
Eppenga
and
D.
Frenkel
, “
Monte Carlo study of the isotropic and nematic phases of infinitely thin hard platelets
,”
Mol. Phys.
52
(
6
),
1303
1334
(
1984
).
93.
V. I.
Harismiadis
,
J.
Vorholz
, and
A. Z.
Panagiotopoulos
, “
Efficient pressure estimation in molecular simulations without evaluating the virial
,”
J. Chem. Phys.
105
(
18
),
8469
8470
(
1996
).
94.
E.
de Miguel
and
G.
Jackson
, “
The nature of the calculation of the pressure in molecular simulations of continuous models from volume perturbations
,”
J. Chem. Phys.
125
(
16
),
164109
(
2006
).
95.
G. J.
Gloor
,
G.
Jackson
,
F. J.
Blas
, and
E.
De Miguel
, “
Test-area simulation method for the direct determination of the interfacial tension of systems with continuous or discontinuous potentials
,”
J. Chem. Phys.
123
(
13
),
134703
(
2005
).
96.
E.
De Miguel
and
G.
Jackson
, “
Detailed examination of the calculation of the pressure in simulations of systems with discontinuous interactions from the mechanical and thermodynamic perspectives
,”
Mol. Phys.
104
(
22–24
),
3717
3734
(
2006
).
97.
G. V.
Lau
,
I. J.
Ford
,
P. A.
Hunt
,
E. A.
Müller
, and
G.
Jackson
, “
Surface thermodynamics of planar, cylindrical, and spherical vapour-liquid interfaces of water
,”
J. Chem. Phys.
142
(
11
),
114701
(
2015
).
98.
J. G.
Sampayo
,
A.
Malijevský
,
E. A.
Müller
,
E.
de Miguel
, and
G.
Jackson
, “
Communications: Evidence for the role of fluctuations in the thermodynamics of nanoscale drops and the implications in computations of the surface tension
,”
J. Chem. Phys.
132
(
14
),
141101
(
2010
).
99.
C.
Ibergay
,
A.
Ghoufi
,
F.
Goujon
,
P.
Ungerer
,
A.
Boutin
,
B.
Rousseau
, and
P.
Malfreyt
, “
Molecular simulations of the n-alkane liquid-vapor interface: Interfacial properties and their long range corrections
,”
Phys. Rev. E
75
(
5
),
051602
(
2007
).
100.
A.
Ghoufi
and
P.
Malfreyt
, “
Local description of surface tension through thermodynamic and mechanical definitions
,”
Mol. Simul.
39
(
8
),
603
611
(
2013
).
101.
A. J. C.
Ladd
and
L. V.
Woodcock
, “
Interfacial and co-existence properties of the Lennard-Jones system at the triple point
,”
Mol. Phys.
36
(
2
),
611
619
(
1978
).
102.
B. D.
Todd
,
P. J.
Daivis
, and
D. J.
Evans
, “
Heat flux vector in highly inhomogeneous nonequilibrium fluids
,”
Phys. Rev. E
51
(
5
),
4362
4368
(
1995
).
103.
E. R.
Smith
,
P. J.
Daivis
, and
B. D.
Todd
, “
Measuring heat flux beyond Fourier’s law
,”
J. Chem. Phys.
150
(
6
),
064103
(
2019
).
104.
P. J.
Daivis
and
B. D.
Todd
, “
A simple, direct derivation and proof of the validity of the SLLOD equations of motion for generalized homogeneous flows
,”
J. Chem. Phys.
124
(
19
),
194103
(
2006
).
105.
A. J.
Chorin
and
O. H.
Hald
,
Computational Statistical Mechanics
(
Elsevier Science
,
Oxford
,
2013
), pp.
157
170
.
106.
A. W.
Lees
and
S. F.
Edwards
, “
The computer study of transport processes under extreme conditions
,”
J. Phys. C: Solid State Phys.
5
(
15
),
1921
1928
(
1972
).
107.
S. Y.
Liem
,
D.
Brown
, and
J. H. R.
Clarke
, “
Investigation of the homogeneous-shear nonequilibrium-molecular-dynamics method
,”
Phys. Rev. A
45
(
6
),
3706
3713
(
1992
).
108.
J.
Petravic
and
P.
Harrowell
, “
The boundary fluctuation theory of transport coefficients in the linear-response limit
,”
J. Chem. Phys.
124
(
1
),
014103
(
2006
).
109.
W.
Noll
, “
Die herleitung der grundgleichungen der thermomechanik der kontinua aus der statistischen mechanik
,”
J. Ration. Mech. Anal.
4
,
627
(
1955
).
110.
R. J.
Hardy
, “
Formulas for determining local properties in molecular‐dynamics simulations: Shock waves
,”
J. Chem. Phys.
76
(
1
),
622
628
(
1982
).
111.
A. I.
Murdoch
, “
A critique of atomistic definitions of the stress tensor
,”
J. Elasticity
88
(
2
),
113
140
(
2007
).
112.
A. I.
Murdoch
, “
On molecular modelling and continuum concepts
,”
J. Elasticity
100
(
1–2
),
33
61
(
2010
).
113.
L. B.
Lucy
, “
Numerical approach to the testing of the fission hypothesis
,”
Astron. J.
82
(
12
),
1013
(
1977
).
114.
T.
Weinhart
,
R.
Hartkamp
,
A. R.
Thornton
, and
S.
Luding
, “
Coarse-grained local and objective continuum description of three-dimensional granular flows down an inclined surface
,”
Phys. Fluids
25
(
7
),
070605
(
2013
).
115.
C.
Hirsch
, “
Introduction
,” in
Numerical Computation of Internal and External Flows
(
Elsevier
,
2007
), pp.
1
5
.
116.
R. B.
Lehoucq
and
A.
Von Lilienfeld-Toal
, “
Translation of Walter Noll’s ‘derivation of the fundamental equations of continuum thermodynamics from statistical mechanics
,’”
J. Elasticity
100
(
1–2
),
5
24
(
2010
).
117.
O. C.
Zienkiewicz
,
R. L.
Taylor
, and
J. Z.
Zhu
,
The Finite Element Method: Its Basis and Fundamentals
, 7th ed. (
Elsevier
,
2013
).
118.
A. I.
Murdoch
and
D.
Bedeaux
, “
On the physical interpretation of fields in continuum mechanics
,”
Int. J. Eng. Sci.
31
(
10
),
1345
1373
(
1993
).
119.
K. P.
Travis
,
B. D.
Todd
, and
D. J.
Evans
, “
Departure from Navier-Stokes hydrodynamics in confined liquids
,”
Phys. Rev. E
55
(
4
),
4288
4295
(
1997
).
120.
K. P.
Travis
and
K. E.
Gubbins
, “
Poiseuille flow of Lennard-Jones fluids in narrow slit pores
,”
J. Chem. Phys.
112
(
4
),
1984
1994
(
2000
).
121.
D. C.
Rapaport
,
The Art of Molecular Dynamics Simulation
, 2nd ed. (
Cambridge University Press
,
Cambridge
,
2004
).
122.
E. M.
Gosling
,
I. R.
Mcdonald
, and
K.
Singer
, “
On the calculation by molecular dynamics of the shear viscosity of a simple fluid
,”
Mol. Phys.
26
,
1475
1484
(
1973
).
123.
D. C.
Rapaport
, “
Microscale hydrodynamics: Discrete-particle simulation of evolving flow patterns
,”
Phys. Rev. A
36
(
7
),
3288
3299
(
1987
).
124.
S. T.
Cui
and
D. J.
Evans
, “
Molecular dynamics simulation of two dimensional flow past a plate
,”
Mol. Simul.
9
(
3
),
179
192
(
1992
).
125.
J.
Sun
,
Y. L.
He
,
Y. S.
Li
, and
W. Q.
Tao
, “
Study on nanoscale obstructed flow with molecular dynamics simulation method
,”
Prog. Comput. Fluid Dyn.
10
(
1
),
51
(
2010
).
126.
D. C.
Rapaport
, “
Molecular-dynamics study of Rayleigh-Bénard convection
,”
Phys. Rev. Lett.
60
(
24
),
2480
2483
(
1988
).
127.
D. C.
Rapaport
, “
Hexagonal convection patterns in atomistically simulated fluids
,”
Phys. Rev. E
73
(
2
),
025301
(
2006
).
128.
K. H.
Han
,
C.
Kim
,
P.
Talkner
,
G. E.
Karniadakis
, and
E. K.
Lee
, “
Molecular hydrodynamics: Vortex formation and sound wave propagation
,”
J. Chem. Phys.
148
(
2
),
024506
(
2018
).
129.
D.
Hirshfeld
and
D. C.
Rapaport
, “
Molecular dynamics simulation of Taylor-Couette vortex formation
,”
Phys. Rev. Lett.
80
(
24
),
5337
5340
(
1998
).
130.
D. J.
Trevelyan
and
T. A.
Zaki
, “
Wavy Taylor vortices in molecular dynamics simulation of cylindrical Couette flow
,”
Phys. Rev. E
93
(
4
),
043107
(
2016
).
131.
K.
Kadau
,
T. C.
Germann
,
N. G.
Hadjiconstantinou
,
P. S.
Lomdahl
,
G.
Dimonte
,
B. L.
Holian
, and
B. J.
Alder
, “
Nanohydrodynamics simulations: An atomistic view of the Rayleigh–Taylor instability
,”
Proc. Natl. Acad. Sci. U. S. A.
101
(
16
),
5851
5855
(
2004
).
132.
K.
Kadau
,
J. L.
Barber
,
T. C.
Germann
,
B. L.
Holian
, and
B. J.
Alder
, “
Atomistic methods in fluid simulation
,”
Philos. Trans. R. Soc., A
368
(
1916
),
1547
1560
(
2010
).
133.
S.
Root
,
R. J.
Hardy
, and
D. R.
Swanson
, “
Continuum predictions from molecular dynamics simulations: Shock waves
,”
J. Chem. Phys.
118
(
7
),
3161
3165
(
2003
).
134.
W. G.
Hoover
and
C. G.
Hoover
, “
Tensor temperature and shock-wave stability in a strong two-dimensional shock wave
,”
Phys. Rev. E
80
(
1
),
011128
(
2009
).
135.
M.
Mareschal
and
B. L.
Holian
,
Microscopic Simulations of Complex Hydrodynamic Phenomena
, NATO ASI Series Vol. 292, edited by
M.
Mareschal
and
B. L.
Holian
(
Springer US
,
Boston, MA
,
1992
).
136.
J.
Jiménez
and
P.
Moin
, “
The minimal flow unit in near-wall turbulence
,”
J. Fluid Mech.
225
,
213
240
(
1991
).
137.
J. M.
Hamilton
,
J.
Kim
, and
F.
Waleffe
, “
Regeneration mechanisms of near-wall turbulence structures
,”
J. Fluid Mech.
287
,
317
348
(
1995
).
138.
E. R.
Smith
, “
A molecular dynamics simulation of the turbulent Couette minimal flow unit
,”
Phys. Fluids
27
(
11
),
115105
(
2015
).
139.
O.
Reynolds
, “
IV. On the dynamical theory of incompressible viscous fluids and the determination of the criterion
,”
Philos. Trans. R. Soc. London, Ser. A
186
,
123
164
(
1895
).
140.
O.
Reynolds
,
Papers on Mechanical and Physical Subjects: Volume 3
, 1st ed. (
Cambridge University Press
,
Cambridge
,
1903
).
141.
A. A.
Townsend
,
The Structure of Turbulent Shear Flow
, 2nd ed., Cambridge Monographs on Mechanics (
Cambridge University Press
,
1980
).
142.
S. B.
Pope
,
Turbulent Flows
, 1st ed. (
Cambridge University Press
,
2000
).
143.
J.
Cormier
,
J. M.
Rickman
, and
T. J.
Delph
, “
Stress calculation in atomistic simulations of perfect and imperfect solids
,”
J. Appl. Phys.
89
(
1
),
99
104
(
2001
).
144.
D. M.
Heyes
,
E. R.
Smith
,
D.
Dini
, and
T. A.
Zaki
, “
The method of planes pressure tensor for a spherical subvolume
,”
J. Chem. Phys.
140
(
5
),
054506
(
2014
).
145.
D. M.
Heyes
,
E. R.
Smith
,
D.
Dini
, and
T. A.
Zaki
, “
The equivalence between volume averaging and method of planes definitions of the pressure tensor at a plane
,”
J. Chem. Phys.
135
(
2
),
024512
(
2011
).
146.
P. J.
Daivis
,
K. P.
Travis
, and
B. D.
Todd
, “
A technique for the calculation of mass, energy, and momentum densities at planes in molecular dynamics simulations
,”
J. Chem. Phys.
104
(
23
),
9651
9653
(
1996
).
147.
M.
Han
and
J. S.
Lee
, “
Method for calculating the heat and momentum fluxes of inhomogeneous fluids
,”
Phys. Rev. E
70
(
6
),
061205
(
2004
).
148.
E. R.
Smith
,
D. M.
Heyes
,
D.
Dini
, and
T. A.
Zaki
, “
Control-volume representation of molecular dynamics
,”
Phys. Rev. E
85
(
5
),
056705
(
2012
).
149.
S. D.
Ramsey
,
K.
Potter
, and
C.
Hansen
, “
Ray bilinear patch intersections
,”
J. Graphics Tools
9
(
3
),
41
47
(
2004
).
150.
H.
Inaoka
and
N.
Ito
, “
Numerical simulation of pool boiling of a Lennard-Jones liquid
,”
Physica A
392
(
18
),
3863
3868
(
2013
).
151.
Y.
Chen
,
J.
Li
,
B.
Yu
,
D.
Sun
,
Y.
Zou
, and
D.
Han
, “
Nanoscale study of bubble nucleation on a cavity substrate using molecular dynamics simulation
,”
Langmuir
34
(
47
),
14234
14248
(
2018
).
152.
A. D.
Lavino
,
E.
Smith
,
M.
Magnini
, and
O. K.
Matar
, “
Surface topography effects on pool boiling via non-equilibrium molecular dynamics simulations
,”
Langmuir
37
(
18
),
5731
5744
(
2021
).
153.
J.
Wen
,
D.
Dini
,
H.
Hu
, and
E. R.
Smith
, “
Molecular droplets vs bubbles: Effect of curvature on surface tension and Tolman length
,”
Phys. Fluids
33
(
7
),
072012
(
2021
).
154.
J. S.
Rowlinson
and
F. L.
Swinton
, “
The thermodynamic properties
,” in
Liquids and Liquid Mixtures
(
Elsevier
,
London
,
1982
), pp.
11
58
.
155.
N. G.
Hadjiconstantinou
,
A. L.
Garcia
,
M. Z.
Bazant
, and
G.
He
, “
Statistical error in particle simulations of hydrodynamic phenomena
,”
J. Comput. Phys.
187
(
1
),
274
297
(
2003
).
156.
Y.
Li
and
G. J.
Wang
, “
How to produce confidence intervals instead of confidence tricks: Representative sampling for molecular simulations of fluid self-diffusion under nanoscale confinement
,”
J. Chem. Phys.
156
(
11
),
114113
(
2022
).