A second order Langevin sampler preserving positive volume for isothermal isobaric ensemble

TL;DR

The paper addresses the challenge of sampling the isothermal-isobaric ensemble with molecular dynamics while guaranteeing a positive simulation-volume. It derives a zero-mass limit and reformulates the dynamics using a log-volume variable $\epsilon=\log V$ with a carefully chosen friction $\gamma(V)=1/(\lambda V^2)$ to yield additive noise and a well-posed, invariant $NPT$ process. A second-order weak-splitting scheme is proposed to preserve positivity of $V$ and achieve second-order accuracy, validated through numerical tests on a free gas, an artificial interacting system, and Lennard–Jones fluids, demonstrating adherence to pressure-virial relations and superior sampling efficiency compared to first-order methods. The result is a robust, higher-order, positive-volume-preserving NPT sampler suitable for realistic MD applications and accurate free-energy and equation-of-state calculations.

Abstract

We propose in this work a second-order Langevin sampler for the isothermal-isobaric ensemble (the NPT ensemble), preserving a positive volume for the simulation box. We first derive the suitable equations of motion for particles to be coupled with the overdamped Langevin equation of volume by sending the artificial mass of the periodic box to zero in the work of Liang et. al. [J. Chem. Phys. 157(14)]. We prove the well-posedness of the new system of equations and show that its invariant measure is the desired ensemble. The new continuous time equations not only justify the previous cell-rescaling methods, but also allow us to choose a suitable friction coefficient so that one has additive noise after a change of variable by taking logarithm of the volume. This observation allows us to propose a second order weak scheme that guarantees the positivity of the volume. Various numerical experiments have been performed to demonstrate the efficacy of our method.

Figures (6)

• Figure 1: Marginal distribution of the stationary distribution concerning $V$. Left panel: the numerical results of the experiment with $N=1$ particle. Right panel: the numerical results of the experiment with $N=100$ particles. In both panels. the red circles are the empirical density obtained by $10^7$ iterations with a time step of $\Delta t=10^{-3}$ while the blue curve is the exact density \ref{['eq:rhoV_toy']}. The empirical density is perfectly overlapping with the exact one.
• Figure 2: Relative error verse time step with the left panel for Scheme \ref{['eq:EM_Vrp']} and the right panel for Algorithm \ref{['alg:2ndLangevinSampler']}. The numerical results obtained by taking the time step $\Delta t_{\text{ref}}=2^{-14}$ serve as the reference solution. The solid curves of various colors present the relative errors in different test functions, as shown in the legend. The magenta dotted curve indicates the second-order convergence rate.
• Figure 3: Relative error verse time step. Left panel: results obtained by Scheme \ref{['eq:EM_Vrp']}. The reference solution is given by the numerical result of taking time step$\Delta t_{\text{ref}} = 2^{-14}$. The right panel: results obtained by Algorithm \ref{['alg:2ndLangevinSampler']}. The reference solution is given by the numerical result of taking time step$\Delta t_{\text{ref}} = 2^{-12}$. For both panels, the solid curves of various colors present the relative errors in different test functions, as shown in the legend. The magenta dotted curve indicates the second-order convergence rate while the black dotted curve indicates the first-order convergence rate.
• Figure 4: Left panel: Ensemble average $\langle P\rangle$ under different pressure $P_0$. Right panel: Ensemble average $\langle PV\rangle$ under different pressure $P_0$. For both panels, the numerical results are from simulations of the system with $N=10$ particles. The ensembles averages (red squares) are finely overlapped with the theoretical values (blue dotted line), which are $P_0$ for $\langle P\rangle$ and $P_0\langle V\rangle-\beta^{-1}$ for $\langle PV\rangle$, respectively, according to the pressure virial theorems Eq. \ref{['eq:P-virial']}.
• Figure 5: The ensemble average of pressure $\langle P\rangle$ verse the ensemble average of density $\langle\rho\rangle=\langle N/V\rangle$. The blue triangles, orange circles, green diamonds and red stars represent the numerical results of the systems with $N=5, 10, 100$ and $200$ particles, respectively. The results of $N=100$ and $N=200$ particles are closely aligned with the reference solution given by the fitting curve in johnson1993lennard (blue curve), indicating that the relation between $\langle P\rangle$ and $\langle \rho\rangle$ converges as $N\rightarrow\infty$.

Theorems & Definitions (15)

• Proposition 2.1
• proof : Proof of Proposition \ref{['prop:convergence']}
• Lemma 2.1
• proof
• Theorem 2.1
• proof
• Definition 3.1
• Proposition 3.1
• Lemma 4.1
• proof