.. _developer-engine-validation:

###################################################
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 :math:`1/r` + LJ for all dimer
  pairs. **Langevin** dynamics at **T = 250 K** with **matched friction**
  :math:`\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:

.. code-block:: pyretis

    conda run -n pyretis python check_forcefield.py

The engines agree to :math:`\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 :math:`\lambda_0` and the dissociation rate constant, converted
to a common unit (1/ps). RETIS, 3000 cycles per engine.

.. table:: 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** :math:`\lambda_0`
  within statistical error (0.161 vs 0.177 /ps, ~1 :math:`\sigma`), and on the
  rate within ~2 :math:`\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``):

.. code-block:: 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.