diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 582c2eb0c..eb6e7eb4b 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -29,6 +29,8 @@ organisation on `GitHub `__. * 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 `__ - February 2026 --------------------------------------------------------------------------------------------- diff --git a/tests/convert/test_openmm_water.py b/tests/convert/test_openmm_water.py new file mode 100644 index 000000000..7d671dda4 --- /dev/null +++ b/tests/convert/test_openmm_water.py @@ -0,0 +1,197 @@ +import math + +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 not math.isnan(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 not math.isnan(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 not math.isnan(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 not math.isnan(e), "Potential energy is NaN" diff --git a/wrapper/Convert/SireOpenMM/openmmmolecule.cpp b/wrapper/Convert/SireOpenMM/openmmmolecule.cpp index 80842c4c8..d3e5ba161 100644 --- a/wrapper/Convert/SireOpenMM/openmmmolecule.cpp +++ b/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 vsite_names = {"EPW", "EP1", "EP2", "MW", "LP1", "LP2"}; + + QSet vsite_idxs; + + if (nats == 4 or 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,151 @@ 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 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 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 or h1_idx < 0 or h2_idx < 0) + { + throw SireError::incompatible_error( + QObject::tr("Molecule %1 contains virtual site atoms (%2) but the expected " + "parent atoms (O/OW, H1/HW1, H2/HW2) could not be found. " + "Cannot register virtual sites for this molecule.") + .arg(number.toString()) + .arg(vsite_idxs.count()), + CODELOC); + } + else + { + 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); } diff --git a/wrapper/Convert/SireOpenMM/openmmmolecule.h b/wrapper/Convert/SireOpenMM/openmmmolecule.h index 8d46bf5ea..dd5481efb 100644 --- a/wrapper/Convert/SireOpenMM/openmmmolecule.h +++ b/wrapper/Convert/SireOpenMM/openmmmolecule.h @@ -22,6 +22,28 @@ SIRE_BEGIN_HEADER namespace SireOpenMM { + /** Describes a virtual site (extra point) in an OpenMM molecule. + * Used for 4- and 5-point water models such as OPC, TIP4P, and TIP5P. + */ + struct VirtualSiteInfo + { + enum Type + { + ThreeParticleAverage, // used for 4-point water (OPC, TIP4P) + OutOfPlane // used for 5-point water (TIP5P) + }; + + Type type; + int vsite_idx; // molecule-local index of the virtual site atom + int p1_idx, p2_idx, p3_idx; // molecule-local parent indices (O, H1, H2) + + // Weights for ThreeParticleAverageSite: pos = w1*p1 + w2*p2 + w3*p3 + double w1, w2, w3; + + // Weights for OutOfPlaneSite: pos = p1 + w12*(p2-p1) + w13*(p3-p1) + wCross*(p2-p1)x(p3-p1) + double w12, w13, wCross; + }; + /** Internal class used to hold all of the extracted information * of an OpenMM Molecule. You should not use this outside * of the sire_to_openmm_system function. It holds lots of scratch @@ -142,6 +164,9 @@ namespace SireOpenMM /** All the CMAP parameters (atom0..4 indices, CMAPParameter) */ QVector> cmap_params; + /** Virtual site definitions for extra-point water models (OPC, TIP4P, TIP5P) */ + QVector virtual_sites; + /** All the constraints */ QVector> constraints; diff --git a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp index b3a6798ed..7a9978632 100644 --- a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp +++ b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp @@ -1807,6 +1807,30 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, } } + // Register virtual sites (OPC, TIP4P, TIP5P, …). + // setVirtualSite() must be called after all particles have been added + // to the System but before the Context is created. + for (const auto &vs : mol.virtual_sites) + { + const int vsite_atom = vs.vsite_idx + start_index; + const int p1 = vs.p1_idx + start_index; + const int p2 = vs.p2_idx + start_index; + const int p3 = vs.p3_idx + start_index; + + if (vs.type == VirtualSiteInfo::ThreeParticleAverage) + { + system.setVirtualSite(vsite_atom, + new OpenMM::ThreeParticleAverageSite( + p1, p2, p3, vs.w1, vs.w2, vs.w3)); + } + else + { + system.setVirtualSite(vsite_atom, + new OpenMM::OutOfPlaneSite( + p1, p2, p3, vs.w12, vs.w13, vs.wCross)); + } + } + // now add all of the bond parameters for (const auto &bond : mol.bond_params) {