On the importance of numerical integration details for homogeneous flow simulation

This paper presents a reversible, energy-conserving numerical integration scheme for the Sllod equations of motion implemented in LAMMPS, demonstrating that such details are critical for eliminating systematic errors in viscosity calculations and accurately simulating homogeneous flows, particularly at high rates.

The Gist

Imagine you are trying to simulate how a fluid, like honey or water, flows through a pipe. In the real world, the fluid moves because of pressure or gravity. In a computer simulation, we want to see how the tiny atoms inside that fluid behave when we "push" them to flow.

This paper is about fixing a very subtle, but very important, bug in the computer code used to simulate this flow. The authors, Stephen Sanderson and Debra Searles, discovered that the standard way computers calculate this flow was slightly "leaking" energy, leading to wrong answers, especially when the fluid is moving very fast.

Here is the breakdown of their discovery using simple analogies:

1. The Problem: The "Leaky Bucket" of Physics

To simulate flow, scientists use a set of rules called the Sllod equations. Think of these equations as the instruction manual for how atoms should move when you push the fluid.

However, when you translate these rules into computer code, you have to break time down into tiny little steps (like frames in a movie). The standard way of doing this in popular software (called LAMMPS) was like trying to walk a tightrope while juggling. It worked okay for slow, gentle flows, but as the flow got faster, the "juggling" got messy.

The Analogy: Imagine you are pushing a heavy box across a floor.

• The Real World: You push, the box moves, and the energy you put in is exactly equal to the movement plus the friction.
• The Old Computer Code: It was like pushing the box, but every time you took a step, the computer accidentally added a tiny bit of extra "phantom" energy or lost a tiny bit of real energy. Over a million steps, this tiny error added up. The box would end up moving faster or slower than it should, or the "friction" (viscosity) would be calculated incorrectly.

2. The Solution: A Perfectly Balanced Scale

The authors designed a new, more precise way to take those tiny time steps. They call it a reversible, energy-conserving integration scheme.

The Analogy: Think of a perfectly balanced scale.

• In the old method, the scale was slightly unbalanced. If you put a weight on one side, the scale would slowly tip over time, even if you didn't touch it.
• The new method is like a scale with a magical counter-weight. No matter how much the fluid flows or how fast the atoms move, the total "energy weight" on the scale stays exactly the same. If the fluid gains kinetic energy from the flow, the system automatically accounts for it so the total energy budget never changes.

3. The "Peculiar" vs. "Lab" Frame: The Moving Train

One of the tricky parts of this simulation is deciding where you are watching from.

• The Lab Frame: You are standing on the ground watching the train go by.
• The Peculiar Frame: You are sitting on the train, watching the other passengers.

The authors found that the old code was mixing these two perspectives up. It was like trying to calculate the speed of a passenger on the train by looking at their speed relative to the ground, but forgetting to subtract the speed of the train itself. This caused a "glitch" where the atoms seemed to gain or lose energy just because the computer changed its point of view.

Their new code keeps the perspective consistent, ensuring that the "passengers" (atoms) don't magically speed up or slow down just because the "train" (the simulation box) is stretching or rotating.

4. The Result: Why It Matters

Why should you care if a computer simulation gets the viscosity (thickness) of a fluid slightly wrong?

• High-Speed Flaws: The error was small at slow speeds, but at high speeds (like in industrial manufacturing or blood flow in arteries), the error became huge. The old code would tell engineers that a fluid was thicker or thinner than it actually was.
• The "Ghost" Error: The authors showed that this energy leak wasn't just a math error; it showed up as a systematic error in the pressure calculations. It was like a ghost in the machine, making the fluid feel "heavier" than it really is.

The Takeaway

The authors took a popular, widely used tool (LAMMPS) and fixed the "tightrope walking" instructions. They made the simulation:

1. More Accurate: It now correctly predicts how fluids behave, especially when they are moving fast or being twisted in complex ways.
2. More Stable: The energy doesn't leak out, so the simulation doesn't drift off into nonsense over long periods.
3. More Versatile: It can handle complex flows (like twisting and stretching at the same time) that the old code struggled with.

In short, they didn't just fix a bug; they polished the lens through which scientists view the microscopic world of flowing fluids, ensuring that what we see on the screen matches reality much more closely.

Technical Summary

1. Problem Statement

The paper addresses critical numerical inaccuracies in the simulation of homogeneous non-equilibrium fluid flows using the Sllod equations of motion. While these equations are the standard for modeling bulk rheological properties (like viscosity) in Molecular Dynamics (MD) without explicit walls, existing implementations in widely used software (specifically LAMMPS) suffer from subtle but significant integration errors.

The core problems identified are:

• Lack of Reversibility and Energy Conservation: Many existing codes do not implement a reversible numerical integration scheme. This leads to a systematic drift in the system's energy, violating the conservation of a specific extended Hamiltonian ($H'$).
• Frame of Reference Errors: Common implementations store velocities in the laboratory frame. When updating particle positions, failing to correctly account for the "peculiar" momentum (momentum relative to the streaming velocity) introduces non-reversible step changes.
• Boundary Condition Handling: In complex flows (mixed shear, elongational), the periodic unit cell deforms. Existing codes often update the unit cell shape at the wrong time step relative to force calculations or fail to handle the coupling between different flow tensor components (e.g., $xy$ and $yz$ shear) and lattice remapping correctly.
• Consequence: These errors manifest as systematic overestimations of shear stress and viscosity, particularly at high flow rates, and produce inaccurate transient responses.

2. Methodology

The authors developed a rigorous, reversible, and energy-conserving integration scheme for the Sllod equations, implemented within the LAMMPS MD package.

A. Theoretical Framework:

• Conserved Quantity ($H'$): They derived a conserved quantity for Sllod dynamics coupled with a Nosé-Hoover thermostat. This quantity ($H'$) includes the internal energy, thermostat variables, and a new variable $\lambda$ (or $\kappa$) representing the energy source/sink for the flow work.
  $$H' = \sum \frac{p_i^2}{2m_i} + \Phi + \frac{1}{2}Qp_\eta^2 + Ndk_BT\eta + \lambda$$
  Ensuring $\dot{H}' = 0$ serves as the primary test for integrator stability.
• Liouville Operator Factorization: The authors applied Trotter factorization to the Liouville operator ($iL$) associated with the equations of motion. They decomposed the operator into uncoupled sub-operators for position ($iL_q$), momentum ($iL_p$), flow gradient ($iL_{\nabla u}$), thermostat ($iL_\eta, iL_{p_\eta}$), and unit cell deformation ($iL_V$).

B. Numerical Scheme:

• Propagator: They constructed a velocity-Verlet-style propagator (Eq. 19) that splits the time step ($\delta t$) into specific sequences. Crucially, the "Sllod bookkeeping" (updating $\lambda$) and thermostat steps are interleaved with position and momentum updates.
• Frame Handling: The scheme explicitly handles the conversion between laboratory frame and peculiar frame velocities. To ensure reversibility, velocities are converted to the peculiar frame before position updates and converted back afterward, preventing artificial jumps in momentum due to streaming velocity changes.
• Unit Cell Integration: The lattice vectors ($a, b, c$) are integrated analytically or via splitting at every time step, immediately following position updates and preceding force calculations. This ensures the force calculation uses the correct, current box geometry.
• General Flow Tensors: Unlike previous work limited to planar shear, this scheme handles arbitrary triangular flow tensors (including mixed shear and elongational flows) and correctly manages the coupling between off-diagonal components during lattice remapping.

C. Validation:

• The new implementation was tested against the standard LAMMPS implementation (fix nvt/sllod).
• Simulations were run on a Lennard-Jones fluid (256 particles) under various conditions: planar shear, mixed shear, biaxial extension, and uniform expansion.
• Metrics included the conservation of $H'$, transient response accuracy, and steady-state viscosity calculations using both Direct Ensemble Averaging and the Transient Time Correlation Function (TTCF) formalism.

3. Key Contributions

1. Derivation of a General Conserved Quantity: Provided a simplified expression for the conserved quantity in thermostatted Sllod dynamics that explicitly accounts for deforming boundary conditions and arbitrary flow tensors.
2. Reversible Integration Algorithm: Developed a specific Trotter-splitting integrator that is reversible and conserves $H'$ to order $O(\delta t^3)$ (with global error $O(\delta t^2)$), suitable for general triangular flow tensors.
3. Correction of Frame-of-Reference Artifacts: Demonstrated that storing velocities in the laboratory frame without careful peculiar-frame conversion leads to systematic errors in pressure and viscosity.
4. Robust Unit Cell Handling: Resolved issues regarding the timing of unit cell updates and the handling of mixed flow components (e.g., $xz$ tilt generation from $xy$ and $yz$ shear) during lattice remapping.
5. LAMMPS Implementation: Released a modified version of LAMMPS (fix nvt/sllod and fix deform) that incorporates these corrections, making high-accuracy non-equilibrium MD accessible to the broader community.

4. Results

• Energy Conservation: The modified scheme conserved the extended Hamiltonian $H'$ to within numerical precision ($<0.002\%$ deviation) for all tested flow types (shear, extensional, mixed). In contrast, the standard LAMMPS implementation showed significant drift, especially in mixed flows and at high shear rates.
• Viscosity Accuracy:
  • Direct Average: The standard implementation overestimated shear viscosity at high shear rates due to the systematic error in the pressure tensor. The new scheme yielded lower, more accurate viscosity values.
  • TTCF Agreement: The Transient Time Correlation Function (TTCF) results, which rely on fluctuations rather than direct averages, agreed well between the old and new schemes. However, the new scheme's direct averages aligned with the TTCF results, confirming that the previous discrepancy was an artifact of the integration error, not a failure of the response theory.
• Transient Responses: The new scheme accurately captured transient responses and mixed flow behaviors, which were distorted in the standard implementation due to improper timing of unit cell updates.
• Time Step Scaling: The error in the rate of change of $H'$ scaled as $O(\delta t^2)$, consistent with the theoretical expectations of the second-order integrator.

5. Significance

This work is significant for the field of computational rheology and non-equilibrium molecular dynamics (NEMD) for several reasons:

• Correcting Systematic Bias: It identifies and fixes a long-standing source of error in widely used simulation codes. Researchers studying high-shear-rate rheology or complex flows may have been reporting inaccurate viscosity values due to these integration artifacts.
• Enabling Complex Flows: By providing a robust method for integrating arbitrary triangular flow tensors, the work enables accurate simulation of complex flow geometries (mixed shear, extensional flows) that were previously difficult or impossible to simulate reliably with standard Sllod implementations.
• Validation of Response Theory: The results reconcile discrepancies between direct ensemble averages and exact response theory (TTCF), reinforcing the validity of using NEMD to predict transport properties when the numerical integration is performed correctly.
• Open Source Impact: By implementing these corrections in LAMMPS, the authors provide the community with a tool to perform high-fidelity simulations of homogeneous flows, ensuring that future studies on fluid rheology are built on a numerically stable foundation.