🔬 physics

Estimating Irving-Kirkwood pressure profiles from per-atom virial and density profiles: a closure for the contour-dependent term in planar liquid films

This paper introduces a closure relation that estimates the contour-dependent Irving-Kirkwood pressure profile from computationally inexpensive per-atom virial data and density profiles, significantly improving accuracy and enabling the resolution of surface tensions in asymmetric liquid films without requiring costly pair-wise calculations.

Original author: Jose L. Rivera

Published 2026-11-02✓ Author reviewed ⓘ

📖 5 min read🧠 Deep dive

CC BY 4.0

Original author: Jose L. Rivera

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written by the author. For technical accuracy, refer to the original paper. Read full disclaimer

In the world of fluids, pressure is not just a single number pushing against a wall; it is a complex, shifting landscape that changes from point to point, especially where a liquid meets a gas. When scientists study a thin sheet of liquid floating in a vacuum, they need to understand exactly how the forces are distributed across that sheet to determine if the film will hold together or break apart. This local pressure is crucial for calculating surface tension, the force that holds the liquid together, but there is a long-standing puzzle in how to measure it. Different ways of calculating this pressure from the motion of individual atoms lead to different maps of the landscape. One method, which is the standard reference, requires tracking every single pair of atoms and the invisible line connecting them as they move, a process that is incredibly slow and computationally expensive. Another method, which modern computers can do almost instantly, simply assigns the force of each atom to its own location. While both methods agree on the total surface tension of the film, they disagree on the local details, creating a discrepancy that matters when scientists try to understand the stability of very thin films or the behavior of different liquids mixed together.

A researcher at the Universidad Michoacana de San Nicolás de Hidalgo has found a way to bridge this gap without paying the heavy computational price. The scientist developed a method to take the fast, cheap measurements and correct them so they match the slow, expensive standard with high precision. The core of the discovery is that the difference between the two methods is not random noise but a predictable pattern driven by the geometry of the liquid. By looking at how dense the liquid is at different points and using a simple model of how atoms usually arrange themselves in a bulk liquid, the researcher created a mathematical correction. This correction acts like a filter that transforms the rough, fast measurement into the precise, standard one. The result is a way to get the high-quality data that was previously only available through expensive simulations, but now it can be done with a fraction of the computing time.

The researcher tested this idea on twenty different simulations of liquid films made of simple atoms, varying the thickness of the films, the temperature, and the range of the forces between atoms. In every case, the corrected fast measurements matched the expensive standard measurements almost perfectly. The remaining difference between the two was so small that it was comparable to the natural statistical noise of the simulation itself. In the liquid interior, the correction was indistinguishable from random fluctuation, and at the surfaces where the liquid meets the gas, the method captured over eighty-six percent of the difference that usually separates the two approaches. This level of accuracy means that for the first time, scientists can study the local pressure profiles of thin films using the fast method and trust that the results are physically real, not just an artifact of how the calculation was done.

One of the most significant applications of this work involves films that have two different surfaces, such as a liquid sheet with a layer of a second substance clinging to only one side. In these cases, the total surface tension of the film is the same regardless of how it is calculated, but the tension of each individual side is not. The choice of calculation method can shift the calculated tension from one side to the other by a significant amount, up to ten percent of the total force. This shift is critical because it determines whether the film is stable or if it will rupture. The new method successfully predicts exactly how much tension is transferred between the two sides, allowing researchers to see the true tension of each face. In simulations where a clean side of a film was compared to a side covered in a different substance, the method correctly identified that the clean side had the same tension as a pure liquid film, while the covered side had a lower tension, a result that matches physical expectations.

The method also works for the forces acting perpendicular to the film, which are often messy and full of artificial patterns when calculated with the fast method. By applying the same correction logic, the researcher was able to remove nearly ninety percent of these spurious patterns, leaving a flat, stable profile that matches the physical requirement that pressure should be constant across the film. This provides a built-in test: if the corrected profile is not flat, the method has failed, but in all the tests performed, the profile became flat, confirming the validity of the approach. The researcher notes that while the method relies on specific assumptions about how atoms interact, it works remarkably well across a wide range of conditions, from very cold, dense liquids to warmer, more diffuse ones.

This breakthrough changes how scientists can approach problems involving thin films and interfaces. Instead of being forced to choose between a fast calculation that gives the wrong local details or a slow calculation that is too expensive to run for many different scenarios, they can now use the fast method and apply this correction to get the best of both worlds. The approach requires no new simulations of the film itself, only the density profile from the fast run and a separate, quick calculation of how atoms are arranged in a standard liquid. This makes it possible to study complex systems, such as mixtures of different liquids or films with adsorbed layers, with a level of detail that was previously out of reach. The work confirms that the local pressure is a real, measurable quantity that can be recovered from standard data, provided the right geometric correction is applied, opening the door to more accurate models of fluid stability and interfacial phenomena.

Technical Summary: Estimating Irving–Kirkwood Pressure Profiles from Per-Atom Virial and Density Profiles

Problem Statement
The local pressure tensor of an inhomogeneous fluid is not unique; it depends on the choice of the "contour" along which the interatomic force is distributed in space. In planar molecular dynamics (MD) simulations, two tangential pressure profiles are commonly used:

Irving–Kirkwood (IK) Profile: The reference choice, where the pair virial is distributed uniformly along the straight line connecting two atoms. Calculating this requires a dedicated, computationally expensive pass over all atom pairs for every slab, often dominating simulation wall time (e.g., 55% of the run time in the author's tests).
Per-Atom Virial (PAV) Profile: Computed at negligible cost by standard MD codes (e.g., LAMMPS) by binning the per-atom virial. For tangential components, this coincides with the Harasima contour.

While integrated quantities like total surface tension ( $\gamma$ ) are invariant to the contour choice, local profiles differ by the derivative of a pair quantity, $Q(z)$ , such that $P_{IK}(z) - P_{PAV}(z) = dQ/dz$ . This difference becomes critical for thin films with non-equivalent interfaces (e.g., asymmetric adsorption), where the contour choice redistributes tension between the two faces. The paper addresses the need to recover the accurate IK profile from the cheap PAV profile without incurring the computational cost of a full IK calculation.

Methodology
The author proposes a closure-based estimator to predict the gauge term $Q(z)$ and its derivative $dQ/dz$ using only the density profile $\rho(z)$ and the bulk radial distribution function $g(r)$ of the coexisting liquid.

Theoretical Formulation:
- The difference between IK and PAV weights is expressed as a total derivative of a ramp function $q(z)$ .
- The quantity $Q(z)$ is defined as the virial-weighted sum of these ramps over all pairs.
- Under the closure approximation $\rho^{(2)}(\mathbf{r}_i, \mathbf{r}_j) \approx \rho(z_i)\rho(z_j)g(r_{ij})$ , $Q(z)$ is reduced to a double integral of the density profile against a kernel $K(\zeta)$ constructed from the pair potential derivative $u'(r)$ and $g(r)$ .
Refinements:
- Pair-Local Interpolation: To handle the interface where the environment varies, the kernel $K$ is interpolated between a "liquid" kernel (using bulk $g(r)$ ) and a "dilute" kernel (using the low-density limit $g(r) = \exp[-u(r)/k_BT]$ ). The interpolation weight depends on the mean density of the pair.
- Mixture Extension: For multi-component systems, the closure is applied to each species pair. For cross-pairs lacking a bulk reference (e.g., dilute adsorbates), a conformal approximation is used, scaling the majority species' $g(r)$ by the pair diameter.
- Normal Component Test: The same closure is applied to the normal component ($zz$). Since the IK normal pressure must be constant (equal to vapor pressure) in mechanical equilibrium, the correction $dQ_1/dz$ serves as a self-consistency test: if the corrected PAV profile is not flat, the closure error is quantified without needing an IK run.
Estimator:
The estimated IK tangential profile is calculated as:
$P_{IK, est}(z) = P_{PAV}(z) + \frac{dQ}{dz}$
where $P_{PAV}$ is the measured profile and $dQ/dz$ is derived from the density profile and bulk $g(r)$ .

Key Results
The method was tested on 20 Lennard-Jones free-standing films with thicknesses ranging from 7 to 25 $\sigma$ , three cutoff radii (3.5, 5.0, 7.5 $\sigma$ ), and two temperatures (0.72 and 0.90).

Accuracy of Tangential Profiles:
- The root-mean-square (rms) residual between the estimated and measured IK profiles is 0.001–0.002 $\epsilon/\sigma^3$ , a tenfold reduction compared to the uncorrected PAV difference (0.010–0.021).
- In the liquid interior, the residual is within a factor of two of the statistical noise at $T=0.72$ and matches the noise level at $T=0.90$ .
- At the interface, the estimator captures 86–93% of the contour difference. The remaining error is a systematic doublet on the liquid flank of the interface, attributed to short pairs straddling the steep density gradient where the closure is less accurate.
- The central-plane anisotropy ( $\Delta P = P_{IK} - P_{PAV}$ ), crucial for mechanical stability criteria, is reproduced to within 0.0012 over a range of 0.03. Uncorrected PAV values underestimate this anisotropy by 22–28% in thin films.
Normal Component:
- Applying the closure to the normal component removes 83–92% of the spurious structure in the PAV $zz$ profile, restoring the flatness required by mechanical equilibrium. This provides a validation metric that requires no IK simulation.
Surface Tension and Partitioning:
- Total surface tension remains invariant across all profiles, as expected.
- For films with non-equivalent interfaces (simulated with a surface-active adsorbate), the method successfully resolves the tension of each face. A clean face recovers the tension of the pure film within statistical error.
- The gauge term $Q(z)$ , representing the tension transferred between interfaces by the contour choice, is recovered to within 0.6% of the total surface tension.

Significance and Claims
The paper claims that the proposed estimator allows researchers to obtain the physically motivated Irving–Kirkwood tangential pressure profile from standard, low-cost Per-Atom Virial (PAV) simulations.

Computational Efficiency: The method eliminates the need for the expensive, CPU-bound IK pair-pass. For the systems studied, it reduces the sampling stage wall time by a factor of roughly 2.2 (from 3.84h to 1.73h), as the IK calculation previously consumed 55% of the runtime.
Physical Meaning: The author argues that for thin films without a bulk interior, the local partition of surface tension is contour-dependent. The IK contour is identified as the physically meaningful choice because it locates stress between the atoms carrying it. The estimator makes this local partition accessible for asymmetric films (e.g., with preferential adsorption) without additional computational cost.
Validation: The method is validated not only against direct IK measurements but also via the "flatness test" of the normal component, which requires no IK data.
Limitations: The author notes that the method relies on the closure approximation, which performs best when the interface is wide (higher temperatures) and the liquid is less structured. The residual error at the interface (approx. 10% of the difference) is attributed to short-range correlations not captured by the simple closure, but this error is deemed small enough not to affect conclusions regarding mechanical stability or surface pressure.

In summary, the paper provides a practical, low-cost computational tool to correct the ubiquitous PAV pressure profiles, enabling accurate local stress analysis and mechanical stability criteria for thin films and complex interfaces.

More like this