Cross-engine validation: water-dimer dissociation

This page records a reproducibility check that the external molecular-dynamics engines driven by PyRETIS reproduce the same rare-event kinetics on one simple molecular system with a shared force field. The system is the dissociation (“splitting”) of a hydrogen-bonded water dimer, with the O–O distance as the order parameter. The case and the scripts below live in examples/validation/engines/water_dimer/.

The point of the test is that GROMACS, LAMMPS and OpenMM integrate the identical potential, so if their rate-constant estimates agree (within statistical error) the PyRETIS machinery that drives them is reproducible across engines.

The shared system

  • Force field: flexible SPC/Fw water (Wu, Tepper and Voth, J. Chem. Phys. 124, 024503, 2006) – harmonic O–H bonds and H–O–H angle, Lennard-Jones on oxygen, point charges. The model is flexible (no SETTLE/SHAKE constraints), so the forces are identical across engines with no constraint-solver divergence.
  • Conditions: gas-phase dimer in a 5 nm box with a bare (unshifted) 2 nm cutoff, so every engine evaluates the same \(1/r\) + LJ for all dimer pairs. Langevin dynamics at T = 250 K with matched friction \(\gamma = 1\)/ps and a 0.5 fs time step. The temperature was chosen from MD scans: at 300 K the dimer splits in < 10 ps (not a rare event); at 150 K it stays bound for > 500 ps (too rare to converge); at 250 K it dissociates in ~100 ps – a genuine but tractable rare event.
  • Order parameter: O–O distance. State A (bound) is O–O < 3.5 Å, state B (dissociated) is O–O > 6.5 Å, with eight interfaces between.

Every engine’s input is generated from a single source, wd_common.py, so each force field is derived from the same numbers.

The force-field gate

Before any rate is compared, check_forcefield.py evaluates identical configurations – including deliberately strained ones that exercise the bond/angle constants and the factor-of-two prefactor convention that differs between LAMMPS and GROMACS/OpenMM – in all three engines and compares the potential energy:

conda run -n pyretis python check_forcefield.py

The engines agree to \(\le 3\times10^{-5}\) relative (the residual is the documented difference between the codes’ electrostatic constants). This gate must pass before a rate comparison means anything: a green RETIS run with a mismatched force field would silently compare different physics.

Engines and samplers

  • GROMACS – native PyRETIS RETIS, canonical streaming engine (class = "gromacs" -> GromacsEngine).
  • OpenMM – native PyRETIS RETIS, in-process engine.
  • LAMMPS – the streaming engine via infinite swapping (class = "lammps"); the rate is obtained from WHAM. This is a different sampler from the other two, so agreement is additionally a cross-sampler check.

Results

Flux through \(\lambda_0\) and the dissociation rate constant, converted to a common unit (1/ps). RETIS, 3000 cycles per engine.

Table 60 Water-dimer dissociation: cross-engine rate comparison (250 K).
engine sampler cycles flux (1/ps) rate (1/ps)
GROMACS native RETIS 3000 0.161 (8 %) 0.0127 (22 %)
OpenMM native RETIS 3000 0.177 (7 %) 0.0076 (18 %)
LAMMPS infinite swapping 3000 0.221 0.0185

The bracketed numbers are relative block-error estimates.

Findings:

  • The two native engines agree on the flux through \(\lambda_0\) within statistical error (0.161 vs 0.177 /ps, ~1 \(\sigma\)), and on the rate within ~2 \(\sigma\). The corresponding escape time (~80–130 ps) matches the direct-MD dissociation timescale (~100 ps), an independent sanity check.
  • The streaming LAMMPS engine (a different engine and sampler) reproduces the same flux and rate to within ~1.4x / ~1.5x. The small systematic offset is consistent with the WHAM point-matching estimator and the infinite-swapping sampler; resolving whether it closes fully would need a longer run.
  • Overall the three engines agree to within a factor of ~1.4 on the flux and ~2.4 on the rate at 3000 cycles – the level expected for RETIS with ~20 % block errors, and well within an order of magnitude. The force-field gate above guarantees they are integrating the same potential, so the residual spread is statistical / estimator-level, not a physics difference.

Reproducing the result

On a node with gmx, lmp and openmm (conda env pyretis):

cd examples/validation/engines/water_dimer
python check_forcefield.py            # the force-field gate (must pass)
python build_inputs.py                # GROMACS + LAMMPS inputs from one source
python compare_engines.py run --steps 3000
python compare_engines.py compare     # flux + rate per engine, in 1/ps

The streaming-LAMMPS load paths are bootstrapped by lammps_stream/generate_water_loadpaths.py (adapted from generate_H2_loadpaths); compare_engines.py calls it automatically.

Note

The LAMMPS case here uses the canonical continuous engine (class = "lammps" -> LAMMPSEngine) via infinite swapping. For native PyRETIS RETIS with LAMMPS, the relaunch-per-step engine class = "lammps_steps" is used instead, but its kick-initiation is slow on this system: lambda_0 sits at the bound-state 99th percentile, so the velocity-kick walk needed to reach it requires many per-kick LAMMPS starts. Making lammps_steps kick-init performant here is a follow-up; the case directory lammps/ is wired for it. GROMACS and OpenMM use native RETIS. Because the rate is sampler-independent, mixing samplers is sound – and the agreement between them is an extra cross-sampler validation.