From 80b6c98b4672f9d7a7a37071edf9a5426776becb Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 19 Mar 2026 17:27:22 +0000 Subject: [PATCH 1/4] Add support for 4- and 5-point water models with OpenMM. --- doc/source/changelog.rst | 2 + tests/convert/test_openmm_water.py | 195 ++++++++++++++++++ wrapper/Convert/SireOpenMM/openmmmolecule.cpp | 176 ++++++++++++++++ wrapper/Convert/SireOpenMM/openmmmolecule.h | 25 +++ .../SireOpenMM/sire_to_openmm_system.cpp | 24 +++ 5 files changed, 422 insertions(+) create mode 100644 tests/convert/test_openmm_water.py 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..476533bdd --- /dev/null +++ b/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" 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) { From f002ba8a5159ff99e21d88bb41607ff3fcf9dfef Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Fri, 20 Mar 2026 09:29:10 +0000 Subject: [PATCH 2/4] Use math.isnan to make assertion explicit. [ci skip] --- tests/convert/test_openmm_water.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/tests/convert/test_openmm_water.py b/tests/convert/test_openmm_water.py index 476533bdd..7d671dda4 100644 --- a/tests/convert/test_openmm_water.py +++ b/tests/convert/test_openmm_water.py @@ -1,3 +1,5 @@ +import math + import pytest import sire as sr @@ -81,7 +83,7 @@ def test_tip3p_no_virtual_sites(tip3p_mols, openmm_platform): assert bad_constraints == 0 e = _potential_energy(omm) - assert e == e, "Potential energy is NaN" + assert not math.isnan(e), "Potential energy is NaN" @_skip_no_openmm @@ -118,7 +120,7 @@ def test_tip4p_virtual_sites(tip4p_mols, openmm_platform): ), f"Particle {i}: expected ThreeParticleAverageSite, got {type(vs).__name__}" e = _potential_energy(omm) - assert e == e, "Potential energy is NaN" + assert not math.isnan(e), "Potential energy is NaN" @_skip_no_openmm @@ -154,7 +156,7 @@ def test_opc_virtual_sites(opc_mols, openmm_platform): ), f"Particle {i}: expected ThreeParticleAverageSite, got {type(vs).__name__}" e = _potential_energy(omm) - assert e == e, "Potential energy is NaN" + assert not math.isnan(e), "Potential energy is NaN" @_skip_no_openmm @@ -192,4 +194,4 @@ def test_tip5p_virtual_sites(tip5p_mols, openmm_platform): ), f"Particle {i}: expected OutOfPlaneSite, got {type(vs).__name__}" e = _potential_energy(omm) - assert e == e, "Potential energy is NaN" + assert not math.isnan(e), "Potential energy is NaN" From 4315572e2b61d0f64e7905271d164cb8e612d9b3 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 23 Mar 2026 10:13:40 +0000 Subject: [PATCH 3/4] Catch edge cases and add functionality for coupling levers. --- corelib/src/libs/SireCAS/lambdaschedule.cpp | 79 ++++++++++++-- corelib/src/libs/SireCAS/lambdaschedule.h | 12 +++ doc/source/changelog.rst | 2 + doc/source/tutorial/part07/02_levers.rst | 25 +++++ tests/cas/test_lambdaschedule.py | 102 ++++++++++++++++++ wrapper/CAS/LambdaSchedule.pypp.cpp | 28 ++++- wrapper/Convert/SireOpenMM/openmmmolecule.cpp | 12 +++ 7 files changed, 253 insertions(+), 7 deletions(-) diff --git a/corelib/src/libs/SireCAS/lambdaschedule.cpp b/corelib/src/libs/SireCAS/lambdaschedule.cpp index 1c0b43d71..652b380e4 100644 --- a/corelib/src/libs/SireCAS/lambdaschedule.cpp +++ b/corelib/src/libs/SireCAS/lambdaschedule.cpp @@ -65,7 +65,7 @@ static RegisterMetaType r_schedule; QDataStream &operator<<(QDataStream &ds, const LambdaSchedule &schedule) { - writeHeader(ds, r_schedule, 3); + writeHeader(ds, r_schedule, 4); SharedDataStream sds(ds); @@ -75,6 +75,7 @@ QDataStream &operator<<(QDataStream &ds, const LambdaSchedule &schedule) << schedule.default_equations << schedule.stage_equations << schedule.mol_schedules + << schedule.coupled_levers << static_cast(schedule); return ds; @@ -91,21 +92,30 @@ QDataStream &operator>>(QDataStream &ds, LambdaSchedule &schedule) { VersionID v = readHeader(ds, r_schedule); - if (v == 1 or v == 2 or v == 3) + if (v == 1 or v == 2 or v == 3 or v == 4) { SharedDataStream sds(ds); sds >> schedule.constant_values; - if (v == 3) + if (v >= 3) sds >> schedule.force_names; sds >> schedule.lever_names >> schedule.stage_names >> schedule.default_equations >> schedule.stage_equations; - if (v == 2 or v == 3) + if (v >= 2) sds >> schedule.mol_schedules; + if (v >= 4) + sds >> schedule.coupled_levers; + else + { + // Populate the default coupling for streams written before v4 + schedule.coupled_levers[_get_lever_name("cmap", "cmap_grid")] = + _get_lever_name("torsion", "torsion_k"); + } + if (v < 3) { // need to make sure that the lever names are namespaced @@ -146,13 +156,17 @@ QDataStream &operator>>(QDataStream &ds, LambdaSchedule &schedule) } } else - throw version_error(v, "1, 2, 3", r_schedule, CODELOC); + throw version_error(v, "1, 2, 3, 4", r_schedule, CODELOC); return ds; } LambdaSchedule::LambdaSchedule() : ConcreteProperty() { + // By default, cmap_grid falls back to torsion_k so that customising + // torsion scaling automatically keeps the CMAP correction in sync. + coupled_levers[_get_lever_name("cmap", "cmap_grid")] = + _get_lever_name("torsion", "torsion_k"); } LambdaSchedule::LambdaSchedule(const LambdaSchedule &other) @@ -162,7 +176,8 @@ LambdaSchedule::LambdaSchedule(const LambdaSchedule &other) force_names(other.force_names), lever_names(other.lever_names), stage_names(other.stage_names), default_equations(other.default_equations), - stage_equations(other.stage_equations) + stage_equations(other.stage_equations), + coupled_levers(other.coupled_levers) { } @@ -1006,6 +1021,32 @@ void LambdaSchedule::removeEquation(const QString &stage, this->stage_equations[idx].remove(_get_lever_name(force, lever)); } +/** Couple the lever 'force::lever' to 'fallback_force::fallback_lever'. + * If no custom equation has been set for 'force::lever' at any stage, + * the equation for 'fallback_force::fallback_lever' will be used instead + * of the stage default. This is a single level of indirection — the + * fallback lever is not itself followed further. + * + * By default, 'cmap::cmap_grid' is coupled to 'torsion::torsion_k' so + * that custom torsion schedules automatically keep the CMAP correction + * in sync. + */ +void LambdaSchedule::coupleLever(const QString &force, const QString &lever, + const QString &fallback_force, + const QString &fallback_lever) +{ + coupled_levers[_get_lever_name(force, lever)] = + _get_lever_name(fallback_force, fallback_lever); +} + +/** Remove any coupling for the lever 'force::lever', reverting it to + * use the stage default equation when no custom equation is set. + */ +void LambdaSchedule::removeCoupledLever(const QString &force, const QString &lever) +{ + coupled_levers.remove(_get_lever_name(force, lever)); +} + /** Return whether or not the specified 'lever' in the specified 'force' * at the specified 'stage' has a custom equation set for it */ @@ -1077,6 +1118,32 @@ SireCAS::Expression LambdaSchedule::_getEquation(int stage, return it.value(); } + // Check coupled levers: if this lever is coupled to another, try to + // find a custom equation for the coupled lever before falling back to + // the stage default. This is a single level of indirection (no recursion) + // to prevent loops. + auto coupled_it = this->coupled_levers.find(lever_name); + + if (coupled_it != this->coupled_levers.end()) + { + const auto &coupled_name = coupled_it.value(); + const int sep = coupled_name.indexOf("::"); + const QString coupled_force = sep >= 0 ? coupled_name.left(sep) : "*"; + const QString coupled_lever_part = sep >= 0 ? coupled_name.mid(sep + 2) : coupled_name; + + it = equations.find(coupled_name); + if (it != equations.end()) + return it.value(); + + it = equations.find(_get_lever_name(coupled_force, "*")); + if (it != equations.end()) + return it.value(); + + it = equations.find(_get_lever_name("*", coupled_lever_part)); + if (it != equations.end()) + return it.value(); + } + // we don't have any match, so return the default equation for this stage return this->default_equations[stage]; } diff --git a/corelib/src/libs/SireCAS/lambdaschedule.h b/corelib/src/libs/SireCAS/lambdaschedule.h index 1ac9fa9f6..4d65f079f 100644 --- a/corelib/src/libs/SireCAS/lambdaschedule.h +++ b/corelib/src/libs/SireCAS/lambdaschedule.h @@ -164,6 +164,12 @@ namespace SireCAS const QString &force = "*", const QString &lever = "*") const; + void coupleLever(const QString &force, const QString &lever, + const QString &fallback_force, + const QString &fallback_lever); + + void removeCoupledLever(const QString &force, const QString &lever); + SireCAS::Expression getEquation(const QString &stage = "*", const QString &force = "*", const QString &lever = "*") const; @@ -248,6 +254,12 @@ namespace SireCAS particular stage */ QVector> stage_equations; + /** Coupled lever fallbacks: if a lever (force::lever) has no custom + equation set, use the equation for the paired lever instead of + falling straight back to the stage default. Stored as + _get_lever_name(force, lever) → _get_lever_name(fb_force, fb_lever). */ + QHash coupled_levers; + /** The symbol used to represent the lambda value */ static SireCAS::Symbol lambda_symbol; diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index eb6e7eb4b..472b773b9 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -31,6 +31,8 @@ organisation on `GitHub `__. * Add support for 4- and 5-point water models in the OpenMM conversion layer. +* Add functionality for coupling one lambda lever to another. + `2025.4.0 `__ - February 2026 --------------------------------------------------------------------------------------------- diff --git a/doc/source/tutorial/part07/02_levers.rst b/doc/source/tutorial/part07/02_levers.rst index 5d240fe34..9d1854653 100644 --- a/doc/source/tutorial/part07/02_levers.rst +++ b/doc/source/tutorial/part07/02_levers.rst @@ -159,6 +159,31 @@ In this case, as we have a perturbable system, the Force objects used are; It uses parameters that are controlled by the ``charge``, ``sigma``, ``epsilon``, ``alpha``, ``kappa``, ``charge_scale`` and ``lj_scale`` levers. +For systems that use a force field with CMAP backbone torsion corrections (such +as AMBER ff19SB), an additional Force object is present: + +* ``cmap``: `OpenMM::CMAPTorsionForce `__. + This models the cross-term energy correction (CMAP) for backbone φ/ψ dihedral + pairs. It uses a parameter controlled by the ``cmap_grid`` lever, which + interpolates the two-dimensional energy grid from its λ=0 value to its λ=1 value. + +.. note:: + + By default, ``cmap::cmap_grid`` is coupled to ``torsion::torsion_k``: + if you set a custom equation for ``torsion_k``, the same equation will + automatically be used for ``cmap_grid``, keeping the CMAP correction + in sync with the torsion perturbation. You can override this by + setting an explicit equation for ``cmap::cmap_grid``:: + + # freeze CMAP at its λ=0 grid while still perturbing torsions + s.set_equation(stage="morph", force="cmap", lever="cmap_grid", + equation=s.initial()) + + To remove the coupling entirely (so that ``cmap_grid`` falls back to + the stage default independently of ``torsion_k``):: + + s.remove_coupled_lever(force="cmap", lever="cmap_grid") + Some levers, like ``bond_length``, are used only by a single Force object. However, others, like ``charge``, are used by multiple Force objects. diff --git a/tests/cas/test_lambdaschedule.py b/tests/cas/test_lambdaschedule.py index 3b4213784..04476aea6 100644 --- a/tests/cas/test_lambdaschedule.py +++ b/tests/cas/test_lambdaschedule.py @@ -154,3 +154,105 @@ def test_lambdaschedule(): def test_has_force_specific_equation(force, lever, contained): l = sr.cas.LambdaSchedule.standard_decouple() assert l.has_force_specific_equation("decouple", force, lever) == contained + + +def test_coupled_lever_default(): + """cmap::cmap_grid should follow torsion::torsion_k by default.""" + l = sr.cas.LambdaSchedule.standard_morph() + morph_equation = (1 - l.lam()) * l.initial() + l.lam() * l.final() + + # With no custom equations, both should return the stage default. + _assert_same_equation( + l.lam(), + l.get_equation(stage="morph", force="torsion", lever="torsion_k"), + morph_equation, + ) + _assert_same_equation( + l.lam(), + l.get_equation(stage="morph", force="cmap", lever="cmap_grid"), + morph_equation, + ) + + +def test_coupled_lever_follows_torsion_k(): + """Setting a custom torsion_k equation should automatically apply to cmap_grid.""" + l = sr.cas.LambdaSchedule.standard_morph() + custom_eq = l.lam() ** 2 * l.final() + (1 - l.lam() ** 2) * l.initial() + + l.set_equation( + stage="morph", force="torsion", lever="torsion_k", equation=custom_eq + ) + + # cmap_grid should now follow the custom torsion_k equation via coupling. + _assert_same_equation( + l.lam(), + l.get_equation(stage="morph", force="cmap", lever="cmap_grid"), + custom_eq, + ) + + +def test_coupled_lever_explicit_override(): + """An explicit cmap_grid equation should take precedence over the coupling.""" + l = sr.cas.LambdaSchedule.standard_morph() + torsion_eq = l.lam() ** 2 * l.final() + (1 - l.lam() ** 2) * l.initial() + cmap_eq = l.initial() # freeze CMAP at λ=0 + + l.set_equation( + stage="morph", force="torsion", lever="torsion_k", equation=torsion_eq + ) + l.set_equation(stage="morph", force="cmap", lever="cmap_grid", equation=cmap_eq) + + # cmap_grid should use its own explicit equation, not torsion_k's. + _assert_same_equation( + l.lam(), + l.get_equation(stage="morph", force="cmap", lever="cmap_grid"), + cmap_eq, + ) + # torsion_k should be unaffected. + _assert_same_equation( + l.lam(), + l.get_equation(stage="morph", force="torsion", lever="torsion_k"), + torsion_eq, + ) + + +def test_remove_coupled_lever(): + """Removing the coupling makes cmap_grid fall back to the stage default.""" + l = sr.cas.LambdaSchedule.standard_morph() + morph_equation = (1 - l.lam()) * l.initial() + l.lam() * l.final() + custom_eq = l.lam() ** 2 * l.final() + (1 - l.lam() ** 2) * l.initial() + + l.set_equation( + stage="morph", force="torsion", lever="torsion_k", equation=custom_eq + ) + l.remove_coupled_lever(force="cmap", lever="cmap_grid") + + # cmap_grid should now use the stage default, not follow torsion_k. + _assert_same_equation( + l.lam(), + l.get_equation(stage="morph", force="cmap", lever="cmap_grid"), + morph_equation, + ) + + +def test_couple_lever_custom(): + """coupleLever can set an arbitrary coupling between levers.""" + l = sr.cas.LambdaSchedule.standard_morph() + custom_eq = l.lam() ** 2 * l.final() + (1 - l.lam() ** 2) * l.initial() + + # Couple bond_k to torsion_k (unusual, but should work). + l.couple_lever( + force="bond", + lever="bond_k", + fallback_force="torsion", + fallback_lever="torsion_k", + ) + l.set_equation( + stage="morph", force="torsion", lever="torsion_k", equation=custom_eq + ) + + _assert_same_equation( + l.lam(), + l.get_equation(stage="morph", force="bond", lever="bond_k"), + custom_eq, + ) diff --git a/wrapper/CAS/LambdaSchedule.pypp.cpp b/wrapper/CAS/LambdaSchedule.pypp.cpp index 34f75ae16..368102494 100644 --- a/wrapper/CAS/LambdaSchedule.pypp.cpp +++ b/wrapper/CAS/LambdaSchedule.pypp.cpp @@ -645,9 +645,35 @@ void register_LambdaSchedule_class(){ , ( bp::arg("stage")="*", bp::arg("force")="*", bp::arg("lever")="*" ) , "Remove the custom equation for the specified `lever` in the\n specified force at the specified `stage`.\n The lever will now use the equation specified for this\n lever for this stage, or the default lever for the stage\n if this isnt set\n" ); + } + { //::SireCAS::LambdaSchedule::coupleLever + + typedef void ( ::SireCAS::LambdaSchedule::*coupleLever_function_type)( ::QString const &,::QString const &,::QString const &,::QString const & ) ; + coupleLever_function_type coupleLever_function_value( &::SireCAS::LambdaSchedule::coupleLever ); + + LambdaSchedule_exposer.def( + "coupleLever" + , coupleLever_function_value + , ( bp::arg("force"), bp::arg("lever"), bp::arg("fallback_force"), bp::arg("fallback_lever") ) + , bp::release_gil_policy() + , "Couple force::lever to fallback_force::fallback_lever. When no\n custom equation is set for force::lever, the equation for\n fallback_force::fallback_lever will be used instead of the stage\n default. By default cmap::cmap_grid is coupled to torsion::torsion_k.\n" ); + + } + { //::SireCAS::LambdaSchedule::removeCoupledLever + + typedef void ( ::SireCAS::LambdaSchedule::*removeCoupledLever_function_type)( ::QString const &,::QString const & ) ; + removeCoupledLever_function_type removeCoupledLever_function_value( &::SireCAS::LambdaSchedule::removeCoupledLever ); + + LambdaSchedule_exposer.def( + "removeCoupledLever" + , removeCoupledLever_function_value + , ( bp::arg("force"), bp::arg("lever") ) + , bp::release_gil_policy() + , "Remove any coupling for force::lever, reverting it to use the\n stage default equation when no custom equation is set.\n" ); + } { //::SireCAS::LambdaSchedule::removeForce - + typedef void ( ::SireCAS::LambdaSchedule::*removeForce_function_type)( ::QString const & ) ; removeForce_function_type removeForce_function_value( &::SireCAS::LambdaSchedule::removeForce ); diff --git a/wrapper/Convert/SireOpenMM/openmmmolecule.cpp b/wrapper/Convert/SireOpenMM/openmmmolecule.cpp index d3e5ba161..23d3f77fc 100644 --- a/wrapper/Convert/SireOpenMM/openmmmolecule.cpp +++ b/wrapper/Convert/SireOpenMM/openmmmolecule.cpp @@ -1399,6 +1399,18 @@ void OpenMMMolecule::constructFromAmber(const Molecule &mol, } } } + + if (virtual_sites.size() != vsite_idxs.count()) + { + throw SireError::incompatible_error( + QObject::tr("Molecule %1: detected %2 virtual site atom(s) but only " + "registered %3. Check atom geometry (degenerate coordinates?) " + "or atom naming conventions.") + .arg(number.toString()) + .arg(vsite_idxs.count()) + .arg(virtual_sites.size()), + CODELOC); + } } this->buildExceptions(mol, constrained_pairs, map); From c81b8d4ab7850f3a9b16d9a165aab24f69147ee2 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 24 Mar 2026 12:28:59 +0000 Subject: [PATCH 4/4] Energy trajectory key is also required in else block. [ci skip] --- src/sire/mol/_dynamics.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/sire/mol/_dynamics.py b/src/sire/mol/_dynamics.py index 7f9ece5c5..ef9bf2e17 100644 --- a/src/sire/mol/_dynamics.py +++ b/src/sire/mol/_dynamics.py @@ -492,11 +492,11 @@ def _exit_dynamics_block( zip(lambda_windows, rest2_scale_factors) ): if lambda_value != sim_lambda_value: + key = f"{lambda_value:.5f}" if ( not has_lambda_index or abs(lambda_index - i) <= num_energy_neighbours ): - key = f"{lambda_value:.5f}" self._omm_mols.set_lambda( lambda_value, rest2_scale=rest2_scale,