Add support for 4- and 5-point water models with OpenMM.

Author: lohedges

doc/source/changelog.rst

@@ -29,6 +29,8 @@ organisation on `GitHub <https://github.com/openbiosim/sire>`__.
 
 * Fix hang in ``sire.load`` function when shared GROMACS topology path is missing.
 
+* Add support for 4- and 5-point water models in the OpenMM conversion layer.
+
 `2025.4.0 <https://github.com/openbiosim/sire/compare/2025.3.0...2025.4.0>`__ - February 2026
 ---------------------------------------------------------------------------------------------
 

tests/convert/test_openmm_water.py

@@ -0,0 +1,195 @@
+import pytest
+import sire as sr
+
+_skip_no_openmm = pytest.mark.skipif(
+    "openmm" not in sr.convert.supported_formats(),
+    reason="openmm support is not available",
+)
+
+
+@pytest.fixture(scope="module")
+def tip3p_mols():
+    return sr.load_test_files("tip3p.s3")
+
+
+@pytest.fixture(scope="module")
+def tip4p_mols():
+    return sr.load_test_files("tip4p.s3")
+
+
+@pytest.fixture(scope="module")
+def tip5p_mols():
+    return sr.load_test_files("tip5p.s3")
+
+
+@pytest.fixture(scope="module")
+def opc_mols():
+    return sr.load_test_files("opc.s3")
+
+
+def _get_vsite_info(omm_context):
+    """Return (n_vsites, n_zero_mass, n_bad_constraints) for an OpenMM context."""
+    import openmm
+
+    sys = omm_context.getSystem()
+    natoms = sys.getNumParticles()
+
+    zero_mass = {
+        i
+        for i in range(natoms)
+        if sys.getParticleMass(i).value_in_unit(openmm.unit.dalton) == 0.0
+    }
+    n_vsites = sum(1 for i in zero_mass if sys.isVirtualSite(i))
+
+    bad_constraints = sum(
+        1
+        for k in range(sys.getNumConstraints())
+        if sys.getConstraintParameters(k)[0] in zero_mass
+        or sys.getConstraintParameters(k)[1] in zero_mass
+    )
+
+    return n_vsites, len(zero_mass), bad_constraints
+
+
+def _potential_energy(omm_context):
+    import openmm
+
+    return (
+        omm_context.getState(getEnergy=True)
+        .getPotentialEnergy()
+        .value_in_unit(openmm.unit.kilojoules_per_mole)
+    )
+
+
+@_skip_no_openmm
+def test_tip3p_no_virtual_sites(tip3p_mols, openmm_platform):
+    """TIP3P has no EP atoms — no virtual sites should be registered."""
+    omm_map = {
+        "constraint": "h-bonds",
+        "cutoff_type": "none",
+        "cutoff": "none",
+        "platform": openmm_platform or "CPU",
+    }
+
+    omm = sr.convert.to(tip3p_mols, "openmm", map=omm_map)
+    n_vsites, n_zero_mass, bad_constraints = _get_vsite_info(omm)
+
+    assert (
+        n_zero_mass == 0
+    ), f"TIP3P should have no zero-mass particles, got {n_zero_mass}"
+    assert n_vsites == 0, f"TIP3P should have no virtual sites, got {n_vsites}"
+    assert bad_constraints == 0
+
+    e = _potential_energy(omm)
+    assert e == e, "Potential energy is NaN"
+
+
+@_skip_no_openmm
+def test_tip4p_virtual_sites(tip4p_mols, openmm_platform):
+    """TIP4P: one ThreeParticleAverageSite per water, no bad constraints, finite energy."""
+    import openmm
+
+    omm_map = {
+        "constraint": "h-bonds",
+        "cutoff_type": "none",
+        "cutoff": "none",
+        "platform": openmm_platform or "CPU",
+    }
+
+    omm = sr.convert.to(tip4p_mols, "openmm", map=omm_map)
+    n_vsites, n_zero_mass, bad_constraints = _get_vsite_info(omm)
+
+    n_waters = tip4p_mols.num_molecules()
+
+    assert (
+        n_zero_mass == n_waters
+    ), f"Expected {n_waters} zero-mass EP atoms, got {n_zero_mass}"
+    assert n_vsites == n_waters, f"Expected {n_waters} virtual sites, got {n_vsites}"
+    assert bad_constraints == 0, f"Constraints on virtual sites: {bad_constraints}"
+
+    # All virtual sites should be ThreeParticleAverageSite
+    sys = omm.getSystem()
+    natoms = sys.getNumParticles()
+    for i in range(natoms):
+        if sys.isVirtualSite(i):
+            vs = sys.getVirtualSite(i)
+            assert isinstance(
+                vs, openmm.ThreeParticleAverageSite
+            ), f"Particle {i}: expected ThreeParticleAverageSite, got {type(vs).__name__}"
+
+    e = _potential_energy(omm)
+    assert e == e, "Potential energy is NaN"
+
+
+@_skip_no_openmm
+def test_opc_virtual_sites(opc_mols, openmm_platform):
+    """OPC: one ThreeParticleAverageSite per water, no bad constraints, finite energy."""
+    import openmm
+
+    omm_map = {
+        "constraint": "h-bonds",
+        "cutoff_type": "none",
+        "cutoff": "none",
+        "platform": openmm_platform or "CPU",
+    }
+
+    omm = sr.convert.to(opc_mols, "openmm", map=omm_map)
+    n_vsites, n_zero_mass, bad_constraints = _get_vsite_info(omm)
+
+    n_waters = opc_mols.num_molecules()
+
+    assert (
+        n_zero_mass == n_waters
+    ), f"Expected {n_waters} zero-mass EP atoms, got {n_zero_mass}"
+    assert n_vsites == n_waters, f"Expected {n_waters} virtual sites, got {n_vsites}"
+    assert bad_constraints == 0, f"Constraints on virtual sites: {bad_constraints}"
+
+    sys = omm.getSystem()
+    natoms = sys.getNumParticles()
+    for i in range(natoms):
+        if sys.isVirtualSite(i):
+            vs = sys.getVirtualSite(i)
+            assert isinstance(
+                vs, openmm.ThreeParticleAverageSite
+            ), f"Particle {i}: expected ThreeParticleAverageSite, got {type(vs).__name__}"
+
+    e = _potential_energy(omm)
+    assert e == e, "Potential energy is NaN"
+
+
+@_skip_no_openmm
+def test_tip5p_virtual_sites(tip5p_mols, openmm_platform):
+    """TIP5P: two OutOfPlaneSite virtual sites per water, no bad constraints, finite energy."""
+    import openmm
+
+    omm_map = {
+        "constraint": "h-bonds",
+        "cutoff_type": "none",
+        "cutoff": "none",
+        "platform": openmm_platform or "CPU",
+    }
+
+    omm = sr.convert.to(tip5p_mols, "openmm", map=omm_map)
+    n_vsites, n_zero_mass, bad_constraints = _get_vsite_info(omm)
+
+    n_waters = tip5p_mols.num_molecules()
+
+    assert (
+        n_zero_mass == 2 * n_waters
+    ), f"Expected {2 * n_waters} zero-mass EP atoms, got {n_zero_mass}"
+    assert (
+        n_vsites == 2 * n_waters
+    ), f"Expected {2 * n_waters} virtual sites, got {n_vsites}"
+    assert bad_constraints == 0, f"Constraints on virtual sites: {bad_constraints}"
+
+    sys = omm.getSystem()
+    natoms = sys.getNumParticles()
+    for i in range(natoms):
+        if sys.isVirtualSite(i):
+            vs = sys.getVirtualSite(i)
+            assert isinstance(
+                vs, openmm.OutOfPlaneSite
+            ), f"Particle {i}: expected OutOfPlaneSite, got {type(vs).__name__}"
+
+    e = _potential_energy(omm)
+    assert e == e, "Potential energy is NaN"

wrapper/Convert/SireOpenMM/openmmmolecule.cpp

@@ -742,6 +742,28 @@ void OpenMMMolecule::constructFromAmber(const Molecule &mol,
         this->unbonded_atoms.insert(i);
     }
 
+    // Detect virtual site (extra point) atoms by name.
+    // Only 4- and 5-atom molecules can have EP atoms (TIP4P, OPC, TIP5P).
+    // Supported names: AMBER EPW / EP1 / EP2, GROMACS MW / LP1 / LP2.
+    // We remove them from unbonded_atoms so that buildExceptions does
+    // not try to add an erroneous constraint for them.
+    static const QSet<QString> vsite_names = {"EPW", "EP1", "EP2", "MW", "LP1", "LP2"};
+
+    QSet<int> vsite_idxs;
+
+    if (nats >= 4 and nats <= 5)
+    {
+        for (int i = 0; i < nats; ++i)
+        {
+            if (masses_data[i] < 0.05 and
+                vsite_names.contains(mol.atom(SireMol::AtomIdx(i)).name().value()))
+            {
+                vsite_idxs.insert(i);
+                unbonded_atoms.remove(i);
+            }
+        }
+    }
+
     // now the bonds
     const double bond_k_to_openmm = 2.0 * (SireUnits::kcal_per_mol / (SireUnits::angstrom * SireUnits::angstrom)).to(SireUnits::kJ_per_mol / (SireUnits::nanometer * SireUnits::nanometer));
     const double bond_r0_to_openmm = SireUnits::angstrom.to(SireUnits::nanometer);
@@ -825,6 +847,11 @@ void OpenMMMolecule::constructFromAmber(const Molecule &mol,
         if (atom0 > atom1)
             std::swap(atom0, atom1);
 
+        // Skip bonds involving virtual site atoms: their positions are
+        // determined by OpenMM's VirtualSite machinery, not a HarmonicBondForce.
+        if (vsite_idxs.contains(atom0) or vsite_idxs.contains(atom1))
+            continue;
+
         const double k = bondparam.k() * bond_k_to_openmm;
         const double r0 = bondparam.r0() * bond_r0_to_openmm;
         double r0_1 = r0;
@@ -972,6 +999,10 @@ void OpenMMMolecule::constructFromAmber(const Molecule &mol,
         if (atom0 > atom2)
             std::swap(atom0, atom2);
 
+        // Skip angles involving virtual site atoms.
+        if (vsite_idxs.contains(atom0) or vsite_idxs.contains(atom1) or vsite_idxs.contains(atom2))
+            continue;
+
         const double k = angparam.k() * angle_k_to_openmm;
         const double theta0 = angparam.theta0(); // already in radians
 
@@ -1225,6 +1256,141 @@ void OpenMMMolecule::constructFromAmber(const Molecule &mol,
         }
     }
 
+    // Build virtual site definitions for any extra-point atoms detected above.
+    // Weights are derived from the actual atomic coordinates so that any 4- or
+    // 5-point water model (OPC, TIP4P, TIP5P, …) is handled automatically.
+    virtual_sites.clear();
+
+    if (not vsite_idxs.isEmpty())
+    {
+        // Map atom name → molecule-local index (first match wins).
+        QHash<QString, int> name_to_idx;
+
+        for (int i = 0; i < nats; ++i)
+        {
+            const QString name = mol.atom(SireMol::AtomIdx(i)).name().value();
+
+            if (not name_to_idx.contains(name))
+                name_to_idx.insert(name, i);
+        }
+
+        auto find_idx = [&](std::initializer_list<const char *> candidates) -> int
+        {
+            for (const char *n : candidates)
+            {
+                auto it = name_to_idx.find(QString(n));
+
+                if (it != name_to_idx.end())
+                    return it.value();
+            }
+
+            return -1;
+        };
+
+        const int o_idx = find_idx({"O", "OW"});
+        const int h1_idx = find_idx({"H1", "HW1"});
+        const int h2_idx = find_idx({"H2", "HW2"});
+
+        if (o_idx >= 0 and h1_idx >= 0 and h2_idx >= 0)
+        {
+            const auto &O = coords[o_idx];
+            const auto &H1 = coords[h1_idx];
+            const auto &H2 = coords[h2_idx];
+
+            // Bisector vector in the molecular plane: (H1-O) + (H2-O)
+            const double bx = (H1[0] - O[0]) + (H2[0] - O[0]);
+            const double by = (H1[1] - O[1]) + (H2[1] - O[1]);
+            const double bz = (H1[2] - O[2]) + (H2[2] - O[2]);
+            const double bisect_sq = bx * bx + by * by + bz * bz;
+
+            // Cross product (H1-O) x (H2-O) — needed for out-of-plane sites (TIP5P)
+            const double d1x = H1[0] - O[0], d1y = H1[1] - O[1], d1z = H1[2] - O[2];
+            const double d2x = H2[0] - O[0], d2y = H2[1] - O[1], d2z = H2[2] - O[2];
+            const double cx = d1y * d2z - d1z * d2y;
+            const double cy = d1z * d2x - d1x * d2z;
+            const double cz = d1x * d2y - d1y * d2x;
+            const double cross_sq = cx * cx + cy * cy + cz * cz;
+
+            // ── 4-point water (OPC, TIP4P): single EP on the bisector ──────
+            const int ep_idx = find_idx({"EPW", "MW"});
+
+            if (ep_idx >= 0 and vsite_idxs.contains(ep_idx) and bisect_sq > 1e-20)
+            {
+                const auto &EP = coords[ep_idx];
+                const double ex = EP[0] - O[0];
+                const double ey = EP[1] - O[1];
+                const double ez = EP[2] - O[2];
+                const double a = (ex * bx + ey * by + ez * bz) / bisect_sq;
+
+                VirtualSiteInfo vs;
+                vs.type = VirtualSiteInfo::ThreeParticleAverage;
+                vs.vsite_idx = ep_idx;
+                vs.p1_idx = o_idx;
+                vs.p2_idx = h1_idx;
+                vs.p3_idx = h2_idx;
+                vs.w1 = 1.0 - 2.0 * a;
+                vs.w2 = a;
+                vs.w3 = a;
+                vs.w12 = vs.w13 = vs.wCross = 0.0;
+                virtual_sites.append(vs);
+            }
+
+            // ── 5-point water (TIP5P): two out-of-plane EPs ─────────────────
+            const int ep1_idx = find_idx({"EP1", "LP1"});
+            const int ep2_idx = find_idx({"EP2", "LP2"});
+
+            if (ep1_idx >= 0 and ep2_idx >= 0 and
+                vsite_idxs.contains(ep1_idx) and vsite_idxs.contains(ep2_idx) and
+                bisect_sq > 1e-20 and cross_sq > 1e-20)
+            {
+                const auto &EP1 = coords[ep1_idx];
+                const auto &EP2 = coords[ep2_idx];
+
+                // a = dot((EP1+EP2)/2 - O, bisect) / bisect_sq
+                const double mx = 0.5 * (EP1[0] + EP2[0]) - O[0];
+                const double my = 0.5 * (EP1[1] + EP2[1]) - O[1];
+                const double mz = 0.5 * (EP1[2] + EP2[2]) - O[2];
+                const double a = (mx * bx + my * by + mz * bz) / bisect_sq;
+
+                // b = dot(EP1-EP2, cross) / (2 * cross_sq)
+                const double dx = EP1[0] - EP2[0];
+                const double dy = EP1[1] - EP2[1];
+                const double dz = EP1[2] - EP2[2];
+                const double b = (dx * cx + dy * cy + dz * cz) / (2.0 * cross_sq);
+
+                // EP1: wCross = +b
+                {
+                    VirtualSiteInfo vs;
+                    vs.type = VirtualSiteInfo::OutOfPlane;
+                    vs.vsite_idx = ep1_idx;
+                    vs.p1_idx = o_idx;
+                    vs.p2_idx = h1_idx;
+                    vs.p3_idx = h2_idx;
+                    vs.w1 = vs.w2 = vs.w3 = 0.0;
+                    vs.w12 = a;
+                    vs.w13 = a;
+                    vs.wCross = b;
+                    virtual_sites.append(vs);
+                }
+
+                // EP2: wCross = -b
+                {
+                    VirtualSiteInfo vs;
+                    vs.type = VirtualSiteInfo::OutOfPlane;
+                    vs.vsite_idx = ep2_idx;
+                    vs.p1_idx = o_idx;
+                    vs.p2_idx = h1_idx;
+                    vs.p3_idx = h2_idx;
+                    vs.w1 = vs.w2 = vs.w3 = 0.0;
+                    vs.w12 = a;
+                    vs.w13 = a;
+                    vs.wCross = -b;
+                    virtual_sites.append(vs);
+                }
+            }
+        }
+    }
+
     this->buildExceptions(mol, constrained_pairs, map);
 }