On the importance of numerical integration details for homogeneous flow simulation

TL;DR

This paper tackles the sensitivity of Sllod-based homogeneous flow simulations to numerical integration details by developing a reversible, energy-conserving integration scheme with $O(δt^3)$ accuracy and implementing it in LAMMPS. A conserved quantity $\mathcal{H}'$, built by augmenting the phase space with a Nosé-Hoover thermostat and an energy reservoir $\kappa$, guides the integration to maintain energy balance under general flow tensors, while a carefully designed operator-splitting propagator ensures reversibility. The authors demonstrate that energy-nonconserving schemes introduce systematic biases in direct averages of the pressure tensor and viscosity at high flow rates, which are mitigated when $\mathcal{H}'$ is preserved and TTCF-based viscosity calculations align with the corrected direct results. Practically, the work provides actionable implementation details for public MD codes, improving the accuracy of transient responses, mixed-flow states, and high-rate rheology in Sllod simulations.

Abstract

The Sllod equations of motion enable modeling of homogeneous flow at the atomic scale, and are commonly used to predict fluid properties such as viscosity. However, few publicly available codes support such simulations, and those that do often do not implement a reversible numerical integration scheme or have other subtle problems. Here, we demonstrate a reversible and energy-conserving integration scheme for the Sllod equations of motion with error on the order of $δt^3$, in line with typical operator splitting integrators used in standard molecular dynamics simulations. We discuss various implementation details, and implement the scheme in LAMMPS where we find that our changes enable more accurate simulation of transient responses, mixed flows, and steady states, especially at high rates of flow. Importantly, we show that a lack of energy conservation can manifest as a systematic error in the direct ensemble average of the pressure tensor, leading to an error in the calculated viscosity which becomes significant at high flow rates.

Figures (4)

• Figure 1: Check for conservation of $\mathcal{H}'$ under planar shear flow ($\dot{\gamma}_{xy} = 0.1$) with two different integration time steps, showing (a) the rate of change in the conserved quantity averaged over 40,000 trajectories and (b) the conserved quantity for a single trajectory. Note, the data from the peculiar and lab-frame implementations overlap. The inset in (a) additionally shows the rate of change in the conserved quantity at equilibrium for comparison (error envelopes omitted for clarity).
• Figure 2: Comparison of the shear rate dependence of the shear viscosity, $\eta = {-\langle P_{xy}\rangle}/{\dot{\gamma}_{xy}}$, using the time- and ensemble-averaged shear pressure over the final 0.5 time units of each simulation as given by the direct average (closed symbols) and the TTCF formalism (open symbols). TTCF results are offset on the $x$ axis for clarity, but correspond to the same shear rate values as those of the direct average results. To reduce uncertainty, 200,000 trajectories were used for $\dot{\gamma}_{xy}\in\{0.001, 0.5, 1\}$ rather than 40,000. To ensure sampling of the steady state, a trajectory length of 2.5 time units was used for $\dot{\gamma}_{xy} = 1$ rather than 1.5 as used for other simulations. Results from the lab-frame implementation overlapped exactly with those of the peculiar frame implementation, and were therefore omitted for clarity.
• Figure 3: Average percentage change in the conserved quantity under (a) mixed shear flow ($\dot{\gamma}_{xy} = 0.1$, $\dot{\gamma}_{xz} = 0.03$ and $\dot{\gamma}_{yz} = 0.5$), (b) biaxial extensional flow ($\dot{\epsilon}_{xx} = \dot{\epsilon}_{yy} = -\frac{1}{2}\dot{\epsilon}_{zz} = 0.05$), (c) mixed biaxial extensional and shear flow ($\dot{\epsilon}_{xx}=\dot{\epsilon}_{yy}=-\frac{1}{2}\dot{\epsilon}_{zz} = 0.05$, $\dot{\gamma}_{xy} = 0.1$, $\dot{\gamma}_{xz} = 0.03$, $\dot{\gamma}_{yz} = 0.5$), and (d) uniform expansion ($\dot{\epsilon}_{xx}=\dot{\epsilon}_{yy}=\dot{\epsilon}_{zz}=0.1$). The insets show only the 'Lab-Frame' and 'Peculiar' curves, which all remain within approximately $0.002\%$ change in $\mathcal{H}'$ over the duration simulated. Note, the 'LMP' implementation does not allow for remapping of the $yz$ tilt under the flows of (a) or (c), and hence the simulation time was reduced in those cases. The non-zero initial values of the 'LMP' curves with extensional flows are caused by a bug which results in the particles being integrated with $\dot{\epsilon}_{xx} = \dot{\epsilon}_{yy} = \dot{\epsilon}_{zz} = 0$ in the first time step.
• Figure 4: Ensemble average rate of change in the conserved quantity, $\langle\dot{\mathcal{H}}'\rangle$, averaged over the duration of the simulation, plotted as a function of the integration time step for (a) planar shear flow and (b) mixed shear and biaxial extensional flow. The solid line shows a least squares fit of the form $a\delta t^2$.