The Quantum Method of Planes -- Local Pressure Definitions for Machine Learning Potentials

TL;DR

The paper addresses the challenge of defining local stress in inhomogeneous systems modeled with quantum‑inspired machine learning potentials. By deriving a local MoP stress from the Irving–Kirkwood framework and implementing it for the MACE ACE‑based potential, it connects Planes‑on‑a‑plane stress to exact momentum and energy conservation in non‑equilibrium MD. The key contributions are a consistent MoP stress and MoP energy flux formulation for MACE, demonstrations of momentum and energy conservation in a water–ZrO2 interface, and open‑source ASE code to reproduce the validations. This work lays the groundwork for incorporating quantum‑level physics into NEMD and molecular fluid dynamics with ML potentials, enabling accurate, locally defined stress and heat flux measurements at interfaces and in transport problems.

Abstract

Stress, or pressure, is a central quantity in engineering and remains vital in molecular modelling. However, the commonly used virial stress tensor is invalid for an inhomogeneous fluid, which is essential in fluid dynamics and non-equilibrium molecular dynamics (NEMD) simulation. This is solved by using the method of planes (MoP), a mechanical form of pressure, simply interpreted as the force divided by area, yet is derived from the firm foundations of statistical mechanics. We present an extension of MoP stress1 to the MACE potential, a particular form of machine learning (ML) potentials allowing the incorporation of quantum mechanical (QM) physics into classical simulation. We present the derivation of this local stress for the MACE potential using the theoretical framework set out by Irving and Kirkwood 2 . For the test case of an interface between water and Zirconium Oxide, we show that the MoP measures the correct force balance while the virial form fails. Further, we demonstrate that this planar definition of stress is valid arbitrarily far from equilibrium, showing exact conservation every timestep in a control volume bounded by MoP planes. This links the stress directly to the conservation equations and demonstrates the validity in non equilibrium molecular dynamics (NEMD) systems. All code to reproduce these validations for any MACE system, together with ASE accelerated code to calculate the MoP, is provided as open source. This work helps build the foundation to extend the ML revolution in materials to NEMD and molecular fluid dynamics modelling.

Figures (4)

• Figure 1: The domain with Zirconium dioxide(ZrO$_2$) walls and water (H$_2$O) in the channel showing all dimensions. The top is shown on the left highlighting the triclinic nature of the system, with the view angle on the right along the tricilinc domain angle of about $\theta = 9^{\circ}$ to highlight the crystal structure. The domain is directly taken from the DFT work of Yang_et_al2021 with the fluid region doubled in size by copying the molecules and re-equilibrating the larger system.
• Figure 2: Pressure in a pure water NPT simulation controlled to 10,000 bars comparing virial pressure ($\mathop{\boldsymbol{\mathop{\mathrm{\boldsymbol{P}}}\limits}} \nolimits^{ IK1}$ from Eq. (\ref{['virial']})) to Method of Planes over 30 bins and 31 planes respectively. Plots are an average over 30,000 samples in time taken every 10 timesteps steps from MD simulations.
• Figure 3: Plot of spatially localised virial pressure ($\mathop{\boldsymbol{\mathop{\mathrm{\boldsymbol{P}}}\limits}} \nolimits^{ IK1}$ from Eq. (\ref{['virial']})) compared to Method of Planes over 400 bins and 401 planes respectively. The full channel is shown in $a)$ with MoP plotted every 5th bin below the molecular locations from a trajectory snapshot (Zr black, oxygen red and hydrogen white) demonstrating how pressure changes from the solid ZrO$_2$ to the liquid H$_2$O regions. The zoomed in near-wall region is shown in $b)$, where the colours/legend is consistent in both plots. Plots are an average over 30,000 samples in time taken every 10 timesteps steps from MD simulations, with symmetry assumed to improve statistics in $b)$, to 60,000 samples. The kinetic pressure and configurational pressure must sum to a constant for $\partial P_{zz} / \partial z = 0$ to be valid, seen as a constant value for the red point with MoP but not for the IK1/virial sum (dashed red line) as the measured pressure is biased by being assigned to the location of the molecules.
• Figure 4: Control volume conservation plots, $a)$ is momentum conservation for a simulation with $\Delta t = 0.5$, where the measured configurational forces and momentum fluxes are equal to momentum change, shown by the sum which is zero to machine precision. $b)$ Energy, at smaller timestep $\Delta t = 0.125$ so crossings are shown as arrows and both pairwise force are compared: with $\widetilde{\mathop{\mathrm{\boldsymbol{J}}}\limits{}^C}_{\!\!\!z}$ based on $\left[\frac{dU}{d\boldsymbol{r}_{ij}} -\frac{dU}{d\boldsymbol{r}_{ji}} \right] \cdot \frac{\boldsymbol{p}_j}{m_j}$ and $\mathop{\mathrm{\boldsymbol{J}}}\limits{}^C_{z}$ based on $\frac{dU_i}{d\boldsymbol{r}_j} \cdot \frac{\boldsymbol{p}_j}{m_j}$. Sum is based on the more accurate $\frac{dU_i}{d\boldsymbol{r}_j} \cdot \frac{\boldsymbol{p}_j}{m_j}$ form. Insert shows full scale, highlighting molecular surface crossings with magnitudes$\to \infty$ as $\Delta t \to 0$. The plot is for volume number 211, roughly in the middle of the channel, although similar plots can be obtained for any volume.