Systematic incorporation of nuclear quantum effects into atomistic simulations by smoothed trajectory analysis

TL;DR

This work presents PIGSTA, a post-processing framework that inserts nuclear quantum effects into atomistic simulations by applying bead-number dependent convolution kernels $g_P(t)$ to (P)IMD trajectories without modifying the dynamics. By deriving kernels from the harmonic reference via $g_P(t) = rac{1}{ pi} ext{int}_0^ty ext{sqrt}igl{w_P( u)igr{)}} ext{cos}( u t) ext{d} u$, PIGSTA achieves exact quantum limits for harmonic systems at any bead number $P$ and markedly improves convergence of thermodynamic and structural observables at finite $P$. It also introduces internally consistent convergence and harmonicity diagnostics based on energy and force estimators, enabling reference-free assessment of bead-number convergence. Applied to ultralow-temperature Zundel cation and ambient liquid water, PIGSTA reproduces converged PIMD results with reduced bead counts and identifies regimes where reduced beads are physically meaningful, all while incurring negligible extra computational cost. Overall, PIGSTA offers a robust, general, and cost-effective tool for incorporating NQEs into large-scale atomistic simulations with built-in convergence diagnostics.

Abstract

Nuclear quantum effects (NQEs) play an essential role in many atomistic systems, yet their explicit inclusion in molecular simulations remains challenging. Path-integral molecular dynamics (PIMD) provides a rigorous framework for incorporating NQEs, but its practical applicability is often limited by the slow and strongly system-dependent convergence with respect to the number of beads. Here we introduce path-integral generalized smoothed trajectory analysis (PIGSTA), a post-processing framework for the systematic incorporation of NQEs into atomistic simulations, using either classical or path-integral molecular dynamics trajectories. By applying analytically defined convolution kernels to simulation trajectories, PIGSTA corrects the frequency-dependent discretization error associated with a finite number of beads, without modifying the underlying dynamics. For harmonic systems, PIGSTA recovers the exact quantum-mechanical limit at any bead number, whereas standard PIMD becomes exact only in the infinite-bead limit. More generally, the method significantly improves the convergence of thermodynamic and structural observables at finite bead numbers and provides an internal, reference-free diagnostic of bead-number convergence based on the consistency of energy and force estimators. We assess PIGSTA for ambient liquid water and for the Zundel cation at ultralow temperature, representing a particularly demanding case for bead-number convergence. In both systems, PIGSTA reproduces the converged PIMD limit and enables physically consistent results at reduced bead numbers when convergence criteria are satisfied. Owing to its post-processing nature and negligible additional computational cost, PIGSTA offers a practical and broadly applicable approach for incorporating NQEs into atomistic simulations.

Figures (10)

• Figure 1: Schematic overview of the trajectory convolution framework underlying GSTA and its path-integral extension. The same trajectory-based convolution framework applies to both classical molecular dynamics (MD/GSTA) and path-integral molecular dynamics (PIMD/PIGSTA). Left: representative nuclear configurations from a (PI)MD simulation of the Zundel cation (H$_5$O$_2^+$) at finite bead number. Middle: schematic representation of the frequency-dependent trajectory convolution operator. Right: the resulting PIGSTA-filtered configurations, yielding a statistical representation closer to the converged quantum distribution without modifying the underlying dynamics. All configurations are aligned by fitting the oxygen atoms prior to visualization.
• Figure 2: PIGSTA convolution kernels. Convolution kernels $g_P(t)$ evaluated on a discrete time grid for $P=1,2,4,6,8,16,$ and $32$ at a temperature of $300$ K using a time step of $0.25$ fs. Colors increase in brightness with increasing $P$.
• Figure 3: Kinetic energy convergence of the Zundel cation. Convergence of the kinetic energy as a function of the number of beads $P$ at $T = 1.67$ K. Results from PIGLET and PIQTB simulations are taken from Schran et al.Schran2018-py. Error bars indicate 95% confidence intervals.
• Figure 4: Potential energy convergence of the Zundel cation. Convergence of the potential energy as a function of the number of beads $P$. For PIGSTA, both the energy-balance estimator $U^{\mathrm{EB}}$ and the recomputed potential energy on the filtered coordinates $U(\tilde{x})$ are shown. PIQTB and PIGLET data are taken from Schran et al.Schran2018-py. Error bars indicate 95% confidence intervals.
• Figure 5: Convergence of the proton-sharing coordinate fluctuations. Convergence of the standard deviation of the proton-sharing coordinate, $\sigma(\delta)$, as a function of the number of beads $P$. Only PIMD and PIGSTA values are shown. PIGSTA values are obtained directly from the filtered coordinates. Error bars indicate 95% confidence intervals.