From 81f6fe4ff6cd1b99e3bce35f8fd12747d95f2f2b Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 18 Mar 2026 09:29:53 +0000 Subject: [PATCH 1/9] Add support for OpenMM CMAP. --- tests/convert/test_openmm_cmap.py | 81 +++++++ .../PerturbableOpenMMMolecule.pypp.cpp | 57 ++++- wrapper/Convert/SireOpenMM/lambdalever.cpp | 75 ++++++ wrapper/Convert/SireOpenMM/openmmmolecule.cpp | 227 ++++++++++++++++++ wrapper/Convert/SireOpenMM/openmmmolecule.h | 24 ++ .../SireOpenMM/sire_to_openmm_system.cpp | 104 +++++++- wrapper/Qt/sireqt_containers.cpp | 2 + 7 files changed, 568 insertions(+), 2 deletions(-) create mode 100644 tests/convert/test_openmm_cmap.py diff --git a/tests/convert/test_openmm_cmap.py b/tests/convert/test_openmm_cmap.py new file mode 100644 index 000000000..19ebe9f70 --- /dev/null +++ b/tests/convert/test_openmm_cmap.py @@ -0,0 +1,81 @@ +import sire as sr +import pytest + + +@pytest.mark.skipif( + "openmm" not in sr.convert.supported_formats(), + reason="openmm support is not available", +) +def test_openmm_cmap_energy(tmpdir, multichain_cmap, openmm_platform): + """ + Verify that Sire correctly adds CMAPTorsionForce to the OpenMM context by + comparing the total potential energy against a context built directly via + the OpenMM Python API from the same input files. + + The multichain_cmap fixture is a periodic solvated system with three protein + chains, each carrying CMAP backbone correction terms. Using it exercises the + multi-molecule CMAP code path in the conversion layer. + """ + import openmm + import openmm.app + import openmm.unit + + mols = sr.system.System() + mols.add(multichain_cmap[0]) + mols.add(multichain_cmap[1]) + + # Sanity-check: at least two molecules must carry CMAP so that the + # multi-chain code path is exercised. + cmap_mol_count = sum(1 for mol in mols.molecules() if mol.has_property("cmap")) + assert ( + cmap_mol_count >= 2 + ), "Expected at least two molecules with CMAP terms in multichain_cmap" + + # Save the Sire system to AMBER files so the direct OpenMM path reads the + # same topology and coordinates that Sire uses internally. + dir_path = str(tmpdir.mkdir("cmap_omm")) + prm7 = str(sr.save(mols, f"{dir_path}/system.prm7")[0]) + rst7 = str(sr.save(mols, f"{dir_path}/system.rst7")[0]) + + platform_name = openmm_platform or "CPU" + + # Create and OpenMM context via Sire's conversion layer, then get the + # potential energy. + sire_map = { + "constraint": "none", + "cutoff": "none", + "cutoff_type": "none", + "platform": platform_name, + } + omm_sire = sr.convert.to(mols, "openmm", map=sire_map) + sire_energy = ( + omm_sire.getState(getEnergy=True) + .getPotentialEnergy() + .value_in_unit(openmm.unit.kilojoules_per_mole) + ) + + # Create an OpenMM context directly from the AMBER files and get the + # potential energy. + prmtop = openmm.app.AmberPrmtopFile(prm7) + inpcrd = openmm.app.AmberInpcrdFile(rst7) + + omm_system = prmtop.createSystem( + nonbondedMethod=openmm.app.NoCutoff, + constraints=None, + ) + + integrator = openmm.VerletIntegrator(0.001) + platform = openmm.Platform.getPlatformByName(platform_name) + omm_context = openmm.Context(omm_system, integrator, platform) + omm_context.setPositions(inpcrd.positions) + if inpcrd.boxVectors is not None: + omm_context.setPeriodicBoxVectors(*inpcrd.boxVectors) + + direct_energy = ( + omm_context.getState(getEnergy=True) + .getPotentialEnergy() + .value_in_unit(openmm.unit.kilojoules_per_mole) + ) + + # Energies should agree to within 1 kJ/mol. + assert sire_energy == pytest.approx(direct_energy, abs=1.0) diff --git a/wrapper/Convert/SireOpenMM/PerturbableOpenMMMolecule.pypp.cpp b/wrapper/Convert/SireOpenMM/PerturbableOpenMMMolecule.pypp.cpp index 95cd479f1..807a90004 100644 --- a/wrapper/Convert/SireOpenMM/PerturbableOpenMMMolecule.pypp.cpp +++ b/wrapper/Convert/SireOpenMM/PerturbableOpenMMMolecule.pypp.cpp @@ -584,9 +584,64 @@ void register_PerturbableOpenMMMolecule_class(){ , bp::release_gil_policy() , "Return the torsion phase parameters of the perturbed state" ); + } + { //::SireOpenMM::PerturbableOpenMMMolecule::getCMAPGrids0 + + typedef ::QVector< double > ( ::SireOpenMM::PerturbableOpenMMMolecule::*getCMAPGrids0_function_type)( ) const; + getCMAPGrids0_function_type getCMAPGrids0_function_value( &::SireOpenMM::PerturbableOpenMMMolecule::getCMAPGrids0 ); + + PerturbableOpenMMMolecule_exposer.def( + "getCMAPGrids0" + , getCMAPGrids0_function_value + , bp::release_gil_policy() + , "Return the flat concatenated CMAP grid values (column-major, kJ/mol) " + "for the reference state. Grid k has getCMAPGridSizes()[k]^2 entries." ); + + } + { //::SireOpenMM::PerturbableOpenMMMolecule::getCMAPGrids1 + + typedef ::QVector< double > ( ::SireOpenMM::PerturbableOpenMMMolecule::*getCMAPGrids1_function_type)( ) const; + getCMAPGrids1_function_type getCMAPGrids1_function_value( &::SireOpenMM::PerturbableOpenMMMolecule::getCMAPGrids1 ); + + PerturbableOpenMMMolecule_exposer.def( + "getCMAPGrids1" + , getCMAPGrids1_function_value + , bp::release_gil_policy() + , "Return the flat concatenated CMAP grid values (column-major, kJ/mol) " + "for the perturbed state. Grid k has getCMAPGridSizes()[k]^2 entries." ); + + } + { //::SireOpenMM::PerturbableOpenMMMolecule::getCMAPGridSizes + + typedef ::QVector< int > ( ::SireOpenMM::PerturbableOpenMMMolecule::*getCMAPGridSizes_function_type)( ) const; + getCMAPGridSizes_function_type getCMAPGridSizes_function_value( &::SireOpenMM::PerturbableOpenMMMolecule::getCMAPGridSizes ); + + PerturbableOpenMMMolecule_exposer.def( + "getCMAPGridSizes" + , getCMAPGridSizes_function_value + , bp::release_gil_policy() + , "Return the grid dimension N for each CMAP torsion (grid is N x N). " + "Entries correspond to the grids in getCMAPGrids0/1." ); + + } + { //::SireOpenMM::PerturbableOpenMMMolecule::getCMAPAtoms + + typedef ::QVector< boost::tuples::tuple< int, int, int, int, int, + boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, + boost::tuples::null_type, boost::tuples::null_type > > + ( ::SireOpenMM::PerturbableOpenMMMolecule::*getCMAPAtoms_function_type)( ) const; + getCMAPAtoms_function_type getCMAPAtoms_function_value( &::SireOpenMM::PerturbableOpenMMMolecule::getCMAPAtoms ); + + PerturbableOpenMMMolecule_exposer.def( + "getCMAPAtoms" + , getCMAPAtoms_function_value + , bp::release_gil_policy() + , "Return the molecule-local 5-atom indices for each CMAP torsion, " + "in the same order as getCMAPGridSizes(). Used for REST2 scaling." ); + } { //::SireOpenMM::PerturbableOpenMMMolecule::isGhostAtom - + typedef bool ( ::SireOpenMM::PerturbableOpenMMMolecule::*isGhostAtom_function_type)( int ) const; isGhostAtom_function_type isGhostAtom_function_value( &::SireOpenMM::PerturbableOpenMMMolecule::isGhostAtom ); diff --git a/wrapper/Convert/SireOpenMM/lambdalever.cpp b/wrapper/Convert/SireOpenMM/lambdalever.cpp index 8651e2951..1c701c85a 100644 --- a/wrapper/Convert/SireOpenMM/lambdalever.cpp +++ b/wrapper/Convert/SireOpenMM/lambdalever.cpp @@ -1107,6 +1107,27 @@ PropertyList LambdaLever::getLeverValues(const QVector &lambda_values, idx += 1; } + const auto morphed_cmap_grid = cache.morph( + schedule, + "cmap", "cmap_grid", + mol.getCMAPGrids0(), + mol.getCMAPGrids1()); + + if (is_first) + { + for (int i = 0; i < morphed_cmap_grid.count(); ++i) + { + column_names.append(QString("cmap-cmap_grid-%1").arg(i + 1)); + lever_values.append(QVector()); + } + } + + for (const auto &val : morphed_cmap_grid) + { + lever_values[idx].append(val); + idx += 1; + } + is_first = false; } @@ -1156,6 +1177,7 @@ double LambdaLever::setLambda(OpenMM::Context &context, auto bondff = this->getForce("bond", system); auto angff = this->getForce("angle", system); auto dihff = this->getForce("torsion", system); + auto cmapff = this->getForce("cmap", system); // we know if we have peturbable ghost atoms if we have the ghost forcefields const bool have_ghost_atoms = (ghost_ghostff != 0 or ghost_nonghostff != 0); @@ -1781,6 +1803,59 @@ double LambdaLever::setLambda(OpenMM::Context &context, morphed_torsion_k[j] * scale); } } + + // update CMAP parameters for this perturbable molecule + start_index = start_idxs.value("cmap", -1); + + if (start_index != -1 and cmapff != 0) + { + const auto &grid0 = perturbable_mol.getCMAPGrids0(); + const auto &grid1 = perturbable_mol.getCMAPGrids1(); + const auto &sizes = perturbable_mol.getCMAPGridSizes(); + const auto &atoms = perturbable_mol.getCMAPAtoms(); + + // morph all grid values together using the lambda schedule + const auto morphed_grids = cache.morph( + schedule, + "cmap", "cmap_grid", + grid0, + grid1); + + int offset = 0; + + for (int j = 0; j < sizes.count(); ++j) + { + const int N = sizes[j]; + const int map_size = N * N; + + // CMAP is always a proper backbone dihedral pair, never improper. + // Apply REST2 scaling if all 5 atoms are within the REST2 region. + double scale = 1.0; + + const auto &cmap_atms = atoms[j]; + + if (perturbable_mol.isRest2(boost::get<0>(cmap_atms)) and + perturbable_mol.isRest2(boost::get<1>(cmap_atms)) and + perturbable_mol.isRest2(boost::get<2>(cmap_atms)) and + perturbable_mol.isRest2(boost::get<3>(cmap_atms)) and + perturbable_mol.isRest2(boost::get<4>(cmap_atms))) + { + scale = rest2_scale; + } + + std::vector energy(map_size); + + for (int k = 0; k < map_size; ++k) + { + energy[k] = morphed_grids[offset + k] * scale; + } + + cmapff->setMapParameters(start_index + j, N, energy); + offset += map_size; + } + + cmapff->updateParametersInContext(context); + } } // update the parameters in the context diff --git a/wrapper/Convert/SireOpenMM/openmmmolecule.cpp b/wrapper/Convert/SireOpenMM/openmmmolecule.cpp index 4552afb66..e39bce583 100644 --- a/wrapper/Convert/SireOpenMM/openmmmolecule.cpp +++ b/wrapper/Convert/SireOpenMM/openmmmolecule.cpp @@ -1166,6 +1166,65 @@ void OpenMMMolecule::constructFromAmber(const Molecule &mol, } } + // now the CMAP terms + cmap_params.clear(); + + const auto cmap_funcs = params.cmapFunctions(); + + if (not cmap_funcs.isEmpty()) + { + for (const auto &func : cmap_funcs.parameters()) + { + const int atom0 = molinfo.atomIdx(func.atom0()).value(); + const int atom1 = molinfo.atomIdx(func.atom1()).value(); + const int atom2 = molinfo.atomIdx(func.atom2()).value(); + const int atom3 = molinfo.atomIdx(func.atom3()).value(); + const int atom4 = molinfo.atomIdx(func.atom4()).value(); + + cmap_params.append(boost::make_tuple(atom0, atom1, atom2, atom3, atom4, + func.parameter())); + } + } + + if (is_perturbable) + { + // Add any CMAP terms that exist only in the other state, + // using the same atoms but a zero grid for this (missing) state. + const auto cmap_funcs1 = params1.cmapFunctions(); + + for (const auto &func1 : cmap_funcs1.parameters()) + { + const int atom0 = molinfo.atomIdx(func1.atom0()).value(); + const int atom1 = molinfo.atomIdx(func1.atom1()).value(); + const int atom2 = molinfo.atomIdx(func1.atom2()).value(); + const int atom3 = molinfo.atomIdx(func1.atom3()).value(); + const int atom4 = molinfo.atomIdx(func1.atom4()).value(); + + bool found = false; + + for (const auto &existing : cmap_params) + { + if (boost::get<0>(existing) == atom0 and + boost::get<1>(existing) == atom1 and + boost::get<2>(existing) == atom2 and + boost::get<3>(existing) == atom3 and + boost::get<4>(existing) == atom4) + { + found = true; + break; + } + } + + if (not found) + { + const int N = func1.parameter().nRows(); + const SireBase::Array2D null_grid(N, N, 0.0); + cmap_params.append(boost::make_tuple(atom0, atom1, atom2, atom3, atom4, + SireMM::CMAPParameter(null_grid))); + } + } + } + this->buildExceptions(mol, constrained_pairs, map); } @@ -1553,6 +1612,89 @@ void OpenMMMolecule::alignInternals(const PropertyMap &map) .arg(perturbed->exception_params.count()), CODELOC); } + + // align the CMAP parameters between the two end states + QVector> cmap_params_1; + cmap_params_1.reserve(cmap_params.count()); + + found_index_0 = QVector(cmap_params.count(), false); + found_index_1 = QVector(perturbed->cmap_params.count(), false); + + for (int i = 0; i < cmap_params.count(); ++i) + { + const auto &cmap0 = cmap_params.at(i); + + const int atom0 = boost::get<0>(cmap0); + const int atom1 = boost::get<1>(cmap0); + const int atom2 = boost::get<2>(cmap0); + const int atom3 = boost::get<3>(cmap0); + const int atom4 = boost::get<4>(cmap0); + + bool found = false; + + for (int j = 0; j < perturbed->cmap_params.count(); ++j) + { + if (not found_index_1[j]) + { + const auto &cmap1 = perturbed->cmap_params.at(j); + + if (boost::get<0>(cmap1) == atom0 and + boost::get<1>(cmap1) == atom1 and + boost::get<2>(cmap1) == atom2 and + boost::get<3>(cmap1) == atom3 and + boost::get<4>(cmap1) == atom4) + { + cmap_params_1.append(cmap1); + found_index_0[i] = true; + found_index_1[j] = true; + found = true; + break; + } + } + } + + if (not found) + { + // add a null CMAP (zero grid) for the missing perturbed state + found_index_0[i] = true; + const int N = boost::get<5>(cmap0).nRows(); + const SireBase::Array2D null_grid(N, N, 0.0); + cmap_params_1.append(boost::make_tuple(atom0, atom1, atom2, atom3, atom4, + SireMM::CMAPParameter(null_grid))); + } + } + + for (int j = 0; j < perturbed->cmap_params.count(); ++j) + { + if (not found_index_1[j]) + { + // add a CMAP term missing in the reference state + const auto &cmap1 = perturbed->cmap_params.at(j); + + const int atom0 = boost::get<0>(cmap1); + const int atom1 = boost::get<1>(cmap1); + const int atom2 = boost::get<2>(cmap1); + const int atom3 = boost::get<3>(cmap1); + const int atom4 = boost::get<4>(cmap1); + + // add a null CMAP to the reference state + const int N = boost::get<5>(cmap1).nRows(); + const SireBase::Array2D null_grid(N, N, 0.0); + cmap_params.append(boost::make_tuple(atom0, atom1, atom2, atom3, atom4, + SireMM::CMAPParameter(null_grid))); + cmap_params_1.append(cmap1); + found_index_1[j] = true; + } + } + + if (found_index_0.indexOf(false) != -1 or found_index_1.indexOf(false) != -1) + { + throw SireError::program_bug(QObject::tr( + "Failed to align the CMAP parameters!"), + CODELOC); + } + + perturbed->cmap_params = cmap_params_1; } /** Internal function that builds all of the exceptions for all of the @@ -2050,6 +2192,52 @@ QVector OpenMMMolecule::getTorsionKs() const return dih_ks; } +/** Return all CMAP grid values (column-major order, kJ/mol) for all CMAP terms + * concatenated. Grid k has getCMAPGridSizes()[k]^2 entries. + */ +QVector OpenMMMolecule::getCMAPGrids() const +{ + const double cmap_k_to_openmm = (SireUnits::kcal_per_mol).to(SireUnits::kJ_per_mol); + + QVector grids; + + for (const auto &cmap : this->cmap_params) + { + const auto ¶m = boost::get<5>(cmap); + // Apply the same AMBER→OpenMM grid transformation used by OpenMM's + // amber_file_parser.py: cyclic N/2 shift in both axes with phi↔psi swap. + // idx = ngrid*((j+ngrid//2)%ngrid)+((i+ngrid//2)%ngrid) + // where i=phi (outer/slow) and j=psi (inner/fast) in the output. + const auto flat_in = param.grid().toColumnMajorVector(); + const int N = param.nRows(); + + for (int phi = 0; phi < N; ++phi) + { + for (int psi = 0; psi < N; ++psi) + { + const int src = ((psi + N / 2) % N) * N + ((phi + N / 2) % N); + grids.append(flat_in[src] * cmap_k_to_openmm); + } + } + } + + return grids; +} + +/** Return the grid dimension N for each CMAP torsion (grid is N x N) */ +QVector OpenMMMolecule::getCMAPGridSizes() const +{ + QVector sizes; + sizes.reserve(this->cmap_params.count()); + + for (const auto &cmap : this->cmap_params) + { + sizes.append(boost::get<5>(cmap).nRows()); + } + + return sizes; +} + /** Return the atom indexes of the atoms in the exceptions, in * exception order for this molecule */ @@ -2237,6 +2425,18 @@ PerturbableOpenMMMolecule::PerturbableOpenMMMolecule(const OpenMMMolecule &mol, is_improper = mol.is_improper; is_rest2 = mol.is_rest2; + // populate the CMAP grid data and atom indices for both end states + cmap_grid0 = mol.getCMAPGrids(); + cmap_grid1 = mol.perturbed->getCMAPGrids(); + cmap_grid_sizes = mol.getCMAPGridSizes(); + + for (const auto &cmap : mol.cmap_params) + { + cmap_atoms.append(boost::make_tuple(boost::get<0>(cmap), boost::get<1>(cmap), + boost::get<2>(cmap), boost::get<3>(cmap), + boost::get<4>(cmap))); + } + bool fix_perturbable_zero_sigmas = false; if (map.specified("fix_perturbable_zero_sigmas")) @@ -2756,6 +2956,33 @@ QVector PerturbableOpenMMMolecule::getIsImproper() const return is_improper; } +/** Return the 5-atom indices (molecule-local) for each CMAP torsion, + * in the same order as getCMAPGridSizes(). Used for REST2 scaling. + */ +QVector> +PerturbableOpenMMMolecule::getCMAPAtoms() const +{ + return cmap_atoms; +} + +/** Return flat concatenated CMAP grid values for state 0 (column-major, kJ/mol) */ +QVector PerturbableOpenMMMolecule::getCMAPGrids0() const +{ + return cmap_grid0; +} + +/** Return flat concatenated CMAP grid values for state 1 (column-major, kJ/mol) */ +QVector PerturbableOpenMMMolecule::getCMAPGrids1() const +{ + return cmap_grid1; +} + +/** Return the grid dimension N for each CMAP torsion (grid is N x N) */ +QVector PerturbableOpenMMMolecule::getCMAPGridSizes() const +{ + return cmap_grid_sizes; +} + /** Return the index of the first atom in the Sire system */ int PerturbableOpenMMMolecule::getStartAtomIdx() const { diff --git a/wrapper/Convert/SireOpenMM/openmmmolecule.h b/wrapper/Convert/SireOpenMM/openmmmolecule.h index a903048cb..5edc96090 100644 --- a/wrapper/Convert/SireOpenMM/openmmmolecule.h +++ b/wrapper/Convert/SireOpenMM/openmmmolecule.h @@ -73,6 +73,9 @@ namespace SireOpenMM QVector getTorsionPhases() const; QVector getTorsionKs() const; + QVector getCMAPGrids() const; + QVector getCMAPGridSizes() const; + QVector> getExceptionAtoms() const; QVector getChargeScales() const; @@ -136,6 +139,9 @@ namespace SireOpenMM /** All the dihedral and improper parameters */ QVector> dih_params; + /** All the CMAP parameters (atom0..4 indices, CMAPParameter) */ + QVector> cmap_params; + /** All the constraints */ QVector> constraints; @@ -274,6 +280,14 @@ namespace SireOpenMM QVector getTorsionPhases0() const; QVector getTorsionPhases1() const; + QVector getCMAPGrids0() const; + QVector getCMAPGrids1() const; + QVector getCMAPGridSizes() const; + + /** Return the 5-atom indices (molecule-local) for each CMAP torsion, + * in the same order as getCMAPGridSizes(). Used for REST2 scaling. */ + QVector> getCMAPAtoms() const; + QVector getChargeScales0() const; QVector getChargeScales1() const; QVector getLJScales0() const; @@ -348,6 +362,16 @@ namespace SireOpenMM QVector charge_scl0, charge_scl1; QVector lj_scl0, lj_scl1; + /** Flat concatenation of all CMAP grid values (column-major, kJ/mol) for + * state 0 and state 1. Grid k occupies cmap_grid_sizes[k]^2 entries. */ + QVector cmap_grid0, cmap_grid1; + + /** The grid dimension N for each CMAP torsion (grid is N x N) */ + QVector cmap_grid_sizes; + + /** Molecule-local 5-atom indices for each CMAP torsion (for REST2 checks) */ + QVector> cmap_atoms; + /** The indexes of atoms that become ghosts in the * perturbed state */ diff --git a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp index 95275dbac..34f0ffe68 100644 --- a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp +++ b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp @@ -21,6 +21,7 @@ #include "SireMol/moleditor.h" #include "SireMM/amberparams.h" +#include "SireMM/cmapparameter.h" #include "SireMM/anglerestraints.h" #include "SireMM/atomljs.h" #include "SireMM/bondrestraints.h" @@ -1158,10 +1159,11 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, // the infomation in ffinfo _set_clj_cutoff(*cljff, ffinfo); - // now create the base bond, angle and torsion forcefields + // now create the base bond, angle, torsion and CMAP forcefields OpenMM::HarmonicBondForce *bondff = new OpenMM::HarmonicBondForce(); OpenMM::HarmonicAngleForce *angff = new OpenMM::HarmonicAngleForce(); OpenMM::PeriodicTorsionForce *dihff = new OpenMM::PeriodicTorsionForce(); + OpenMM::CMAPTorsionForce *cmapff = new OpenMM::CMAPTorsionForce(); // now create the engine for computing QM forces on atoms QMForce *qmff = 0; @@ -1242,6 +1244,9 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, lambda_lever.addLever("torsion_phase"); lambda_lever.addLever("torsion_k"); + lambda_lever.setForceIndex("cmap", system.addForce(cmapff)); + lambda_lever.addLever("cmap_grid"); + /// /// Stage 4 - define the forces for ghost atoms /// @@ -1584,6 +1589,12 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, from_ghost_idxs.reserve(n_ghost_atoms); to_ghost_idxs.reserve(n_ghost_atoms); + // map from CMAPParameter to OpenMM map index for non-perturbable molecules + // (enables sharing identical grids across molecules) + QHash shared_cmap_map_indices; + + const double cmap_k_to_openmm = (SireUnits::kcal_per_mol).to(SireUnits::kJ_per_mol); + // loop over every molecule and add them one by one for (int i = 0; i < nmols; ++i) { @@ -1632,6 +1643,7 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, start_indicies.insert("bond", bondff->getNumBonds()); start_indicies.insert("angle", angff->getNumAngles()); start_indicies.insert("torsion", dihff->getNumTorsions()); + start_indicies.insert("cmap", cmapff->getNumMaps()); // we can now record this as a perturbable molecule // in the lambda lever. The returned index is the @@ -1816,6 +1828,96 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, boost::get<4>(dih), boost::get<5>(dih), boost::get<6>(dih)); } + // now add all of the CMAP torsions + // CMAP defines two consecutive dihedrals sharing atoms 1-3: + // phi: atom0-atom1-atom2-atom3 + // psi: atom1-atom2-atom3-atom4 + if (any_perturbable and mol.isPerturbable()) + { + // perturbable molecules get a dedicated map per CMAP torsion + // (so the map values can be updated at each lambda step) + for (const auto &cmap : mol.cmap_params) + { + const auto ¶m = boost::get<5>(cmap); + // Apply the same AMBER→OpenMM grid transformation used by OpenMM's + // amber_file_parser.py: cyclic N/2 shift in both axes with phi↔psi swap. + // idx = ngrid*((j+ngrid//2)%ngrid)+((i+ngrid//2)%ngrid) + // where i=phi (outer/slow) and j=psi (inner/fast) in the output. + const auto flat_in = param.grid().toColumnMajorVector(); + const int N = param.nRows(); + + std::vector energy(N * N); + + for (int phi = 0; phi < N; ++phi) + { + for (int psi = 0; psi < N; ++psi) + { + const int src = ((psi + N / 2) % N) * N + ((phi + N / 2) % N); + energy[phi * N + psi] = flat_in[src] * cmap_k_to_openmm; + } + } + + const int map_index = cmapff->addMap(N, energy); + + cmapff->addTorsion(map_index, + boost::get<0>(cmap) + start_index, + boost::get<1>(cmap) + start_index, + boost::get<2>(cmap) + start_index, + boost::get<3>(cmap) + start_index, + boost::get<1>(cmap) + start_index, + boost::get<2>(cmap) + start_index, + boost::get<3>(cmap) + start_index, + boost::get<4>(cmap) + start_index); + } + } + else + { + // non-perturbable molecules share maps (deduplicated by CMAPParameter) + for (const auto &cmap : mol.cmap_params) + { + const auto ¶m = boost::get<5>(cmap); + int map_index; + + if (shared_cmap_map_indices.contains(param)) + { + map_index = shared_cmap_map_indices[param]; + } + else + { + // Apply the same AMBER→OpenMM grid transformation used by OpenMM's + // amber_file_parser.py: cyclic N/2 shift in both axes with phi↔psi swap. + // idx = ngrid*((j+ngrid//2)%ngrid)+((i+ngrid//2)%ngrid) + // where i=phi (outer/slow) and j=psi (inner/fast) in the output. + const auto flat_in = param.grid().toColumnMajorVector(); + const int N = param.nRows(); + + std::vector energy(N * N); + + for (int phi = 0; phi < N; ++phi) + { + for (int psi = 0; psi < N; ++psi) + { + const int src = ((psi + N / 2) % N) * N + ((phi + N / 2) % N); + energy[phi * N + psi] = flat_in[src] * cmap_k_to_openmm; + } + } + + map_index = cmapff->addMap(N, energy); + shared_cmap_map_indices.insert(param, map_index); + } + + cmapff->addTorsion(map_index, + boost::get<0>(cmap) + start_index, + boost::get<1>(cmap) + start_index, + boost::get<2>(cmap) + start_index, + boost::get<3>(cmap) + start_index, + boost::get<1>(cmap) + start_index, + boost::get<2>(cmap) + start_index, + boost::get<3>(cmap) + start_index, + boost::get<4>(cmap) + start_index); + } + } + for (const auto &constraint : mol.constraints) { const auto atom0 = boost::get<0>(constraint); diff --git a/wrapper/Qt/sireqt_containers.cpp b/wrapper/Qt/sireqt_containers.cpp index 8e7c5c4d2..7bbdf83ca 100644 --- a/wrapper/Qt/sireqt_containers.cpp +++ b/wrapper/Qt/sireqt_containers.cpp @@ -92,10 +92,12 @@ void register_SireQt_containers() register_tuple>(); register_tuple>(); + register_tuple>(); register_tuple, QVector, QVector>>(); register_list>>(); register_list>>(); + register_list>>(); register_dict>>(); From 3f66444d4d0826c2348e78b9d4b8fd9a7accab9c Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 18 Mar 2026 12:29:33 +0000 Subject: [PATCH 2/9] Expose CMAP merge. --- corelib/src/libs/SireSystem/merge.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/corelib/src/libs/SireSystem/merge.cpp b/corelib/src/libs/SireSystem/merge.cpp index d01703f76..e62d41d5c 100644 --- a/corelib/src/libs/SireSystem/merge.cpp +++ b/corelib/src/libs/SireSystem/merge.cpp @@ -131,6 +131,7 @@ namespace SireSystem "atomtype", "bond", "charge", + "cmap", "connectivity", "coordinates", "dihedral", From d688651fd9af96c0a3c04ec3a26da5f48fb33d3d Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 18 Mar 2026 12:33:55 +0000 Subject: [PATCH 3/9] Add test for perturbable CMAPs. --- tests/convert/test_openmm_cmap.py | 115 +++++++++++++++++++++++++++++- 1 file changed, 114 insertions(+), 1 deletion(-) diff --git a/tests/convert/test_openmm_cmap.py b/tests/convert/test_openmm_cmap.py index 19ebe9f70..c3cc04092 100644 --- a/tests/convert/test_openmm_cmap.py +++ b/tests/convert/test_openmm_cmap.py @@ -39,7 +39,7 @@ def test_openmm_cmap_energy(tmpdir, multichain_cmap, openmm_platform): platform_name = openmm_platform or "CPU" - # Create and OpenMM context via Sire's conversion layer, then get the + # Create an OpenMM context via Sire's conversion layer, then get the # potential energy. sire_map = { "constraint": "none", @@ -79,3 +79,116 @@ def test_openmm_cmap_energy(tmpdir, multichain_cmap, openmm_platform): # Energies should agree to within 1 kJ/mol. assert sire_energy == pytest.approx(direct_energy, abs=1.0) + + +@pytest.mark.slow +@pytest.mark.skipif( + "openmm" not in sr.convert.supported_formats(), + reason="openmm support is not available", +) +def test_openmm_cmap_perturbable(multichain_cmap, openmm_platform): + """ + Verify that CMAPTorsionForce grids are correctly handled for a perturbable + molecule. In practice, CMAP backbone terms are not perturbed in protein FEP, + so both end states carry identical CMAP. The test checks that the perturbable + code path correctly applies the same grids at all lambda values and that + set_lambda does not corrupt the force parameters. + + A region-of-interest merge is used to avoid the O(n²) intrascale merge cost + for a large protein. + """ + import openmm + import BioSimSpace as BSS + + platform_name = openmm_platform or "CPU" + + mol0 = multichain_cmap[0] + + # Both end states are identical: same molecule, same CMAP everywhere. + mol0_bss = BSS._SireWrappers.Molecule(mol0) + mol1_bss = BSS._SireWrappers.Molecule(mol0.clone()) + mapping = {i: i for i in range(mol0.num_atoms())} + pert_mol = BSS.Align.merge(mol0_bss, mol1_bss, mapping=mapping, roi=[1, 2, 3]) + + mols_pert = sr.system.System() + mols_pert.add(pert_mol._sire_object) + mols_pert = sr.morph.link_to_reference(mols_pert) + + omm_map = { + "constraint": "none", + "cutoff": "none", + "cutoff_type": "none", + "platform": platform_name, + } + + def get_cmap_torsion_grids(context): + """ + Return list of (size, grid) for each CMAP torsion, dereferencing the + map index. Grid values are returned as plain floats (kJ/mol) so that + pytest.approx can compare them. This is map-count-agnostic: the + non-perturbable path deduplicates maps while the perturbable path + allocates one map per torsion, but the per-torsion grid values must + agree. + """ + system = context.getSystem() + for force in system.getForces(): + if isinstance(force, openmm.CMAPTorsionForce): + maps = [] + for i in range(force.getNumMaps()): + size, grid = force.getMapParameters(i) + grid_floats = [ + v.value_in_unit(openmm.unit.kilojoules_per_mole) for v in grid + ] + maps.append((size, grid_floats)) + result = [] + for t in range(force.getNumTorsions()): + map_idx = force.getTorsionParameters(t)[0] + result.append(maps[map_idx]) + return result + return [] + + def unique_grids(torsion_grids, decimals=3): + """Return the sorted set of unique (size, rounded-grid) tuples. + + Torsion ordering can differ between the perturbable and non-perturbable + code paths, so we compare the sets of unique grid shapes rather than + comparing torsion-by-torsion.""" + seen = set() + result = [] + for size, grid in torsion_grids: + key = (size, tuple(round(v, decimals) for v in grid)) + if key not in seen: + seen.add(key) + result.append(key) + return sorted(result) + + # Reference: non-perturbable molecule. + mols_ref = sr.system.System() + mols_ref.add(mol0) + omm_ref = sr.convert.to(mols_ref, "openmm", map=omm_map) + ref_torsion_grids = get_cmap_torsion_grids(omm_ref) + + assert len(ref_torsion_grids) > 0, "Reference context has no CMAP torsions" + ref_unique = unique_grids(ref_torsion_grids) + + # Perturbable context — one context, lambda changed in place. + omm_pert = sr.convert.to(mols_pert, "openmm", map=omm_map) + + # At both lambda=0 and lambda=1 the set of unique CMAP grids must match the + # non-perturbable reference (cmap0 == cmap1 for an identity perturbation). + # We compare sets of unique grids because the perturbable and non-perturbable + # code paths may order torsions differently. + for lam in (0.0, 1.0): + omm_pert.set_lambda(lam) + pert_torsion_grids = get_cmap_torsion_grids(omm_pert) + + assert len(pert_torsion_grids) == len(ref_torsion_grids), ( + f"Wrong number of CMAP torsions at lambda={lam}: " + f"{len(pert_torsion_grids)} != {len(ref_torsion_grids)}" + ) + + pert_unique = unique_grids(pert_torsion_grids) + assert pert_unique == ref_unique, ( + f"Set of unique CMAP grids differs from reference at lambda={lam}: " + f"{len(pert_unique)} unique grids vs {len(ref_unique)} in reference" + ) From c518f9a5b82f1f774ff3cd9cf0d2c9dfe79c5604 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 18 Mar 2026 12:35:19 +0000 Subject: [PATCH 4/9] Update CHANGELOG. --- doc/source/changelog.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index bfe5a1bff..920a974cb 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -25,6 +25,8 @@ organisation on `GitHub `__. * Fixed parsing of AMBER and GROMACS GLYCAM force field topologies. +* Add support for CMAP terms in the OpenMM conversion layer. + `2025.4.0 `__ - February 2026 --------------------------------------------------------------------------------------------- From 1342b2b2cb0db66d09b599ff20378cd6c17a50c9 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 18 Mar 2026 12:40:14 +0000 Subject: [PATCH 5/9] Add timeout to avoid hang when there is no internet connection. --- doc/source/changelog.rst | 2 ++ src/sire/_load.py | 13 +++++++++---- 2 files changed, 11 insertions(+), 4 deletions(-) diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 920a974cb..582c2eb0c 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -27,6 +27,8 @@ organisation on `GitHub `__. * Add support for CMAP terms in the OpenMM conversion layer. +* Fix hang in ``sire.load`` function when shared GROMACS topology path is missing. + `2025.4.0 `__ - February 2026 --------------------------------------------------------------------------------------------- diff --git a/src/sire/_load.py b/src/sire/_load.py index 27891d69d..66334ce8a 100644 --- a/src/sire/_load.py +++ b/src/sire/_load.py @@ -77,11 +77,16 @@ def _get_gromacs_dir(): if not os.path.exists(gromacs_tbz2): try: + import socket import urllib.request - urllib.request.urlretrieve(f"{tutorial_url}/gromacs.tar.bz2", gromacs_tbz2) - except Exception: - # we cannot download - just give up + with urllib.request.urlopen( + f"{tutorial_url}/gromacs.tar.bz2", timeout=5 + ) as response: + with open(gromacs_tbz2, "wb") as f: + f.write(response.read()) + except (Exception, socket.timeout): + # we cannot download - continue without GROMACS return None if not os.path.exists(gromacs_tbz2): @@ -443,7 +448,7 @@ def load( gromacs_path = _get_gromacs_dir() m = { - "GROMACS_PATH": _get_gromacs_dir(), + "GROMACS_PATH": gromacs_path, "show_warnings": show_warnings, "parallel": parallel, "ignore_topology_frame": ignore_topology_frame, From 7660ee133848fb3f0a97c9e84a9b0a03da0e28eb Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 18 Mar 2026 12:57:37 +0000 Subject: [PATCH 6/9] Use pre-merge system avoid BioSimSpace dependency and speed up test. --- tests/convert/test_openmm_cmap.py | 22 ++++++---------------- 1 file changed, 6 insertions(+), 16 deletions(-) diff --git a/tests/convert/test_openmm_cmap.py b/tests/convert/test_openmm_cmap.py index c3cc04092..b54c86309 100644 --- a/tests/convert/test_openmm_cmap.py +++ b/tests/convert/test_openmm_cmap.py @@ -81,7 +81,6 @@ def test_openmm_cmap_energy(tmpdir, multichain_cmap, openmm_platform): assert sire_energy == pytest.approx(direct_energy, abs=1.0) -@pytest.mark.slow @pytest.mark.skipif( "openmm" not in sr.convert.supported_formats(), reason="openmm support is not available", @@ -89,29 +88,20 @@ def test_openmm_cmap_energy(tmpdir, multichain_cmap, openmm_platform): def test_openmm_cmap_perturbable(multichain_cmap, openmm_platform): """ Verify that CMAPTorsionForce grids are correctly handled for a perturbable - molecule. In practice, CMAP backbone terms are not perturbed in protein FEP, - so both end states carry identical CMAP. The test checks that the perturbable - code path correctly applies the same grids at all lambda values and that + molecule. The pre-merged stream file merged_molecule_cmap.s3 contains a + perturbable molecule whose two end states are identical (an identity + perturbation of a CHARMM protein chain), so both end states carry the same + CMAP backbone correction terms. The test checks that the perturbable code + path correctly applies the same grids at all lambda values and that set_lambda does not corrupt the force parameters. - - A region-of-interest merge is used to avoid the O(n²) intrascale merge cost - for a large protein. """ import openmm - import BioSimSpace as BSS platform_name = openmm_platform or "CPU" mol0 = multichain_cmap[0] - # Both end states are identical: same molecule, same CMAP everywhere. - mol0_bss = BSS._SireWrappers.Molecule(mol0) - mol1_bss = BSS._SireWrappers.Molecule(mol0.clone()) - mapping = {i: i for i in range(mol0.num_atoms())} - pert_mol = BSS.Align.merge(mol0_bss, mol1_bss, mapping=mapping, roi=[1, 2, 3]) - - mols_pert = sr.system.System() - mols_pert.add(pert_mol._sire_object) + mols_pert = sr.load_test_files("merged_molecule_cmap.s3") mols_pert = sr.morph.link_to_reference(mols_pert) omm_map = { From 7d63fabe794d49b3c5fb45c9680a3baa3ba7afce Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 18 Mar 2026 13:58:14 +0000 Subject: [PATCH 7/9] Add third chain. --- tests/convert/test_openmm_cmap.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/convert/test_openmm_cmap.py b/tests/convert/test_openmm_cmap.py index b54c86309..2b9b8813d 100644 --- a/tests/convert/test_openmm_cmap.py +++ b/tests/convert/test_openmm_cmap.py @@ -23,6 +23,7 @@ def test_openmm_cmap_energy(tmpdir, multichain_cmap, openmm_platform): mols = sr.system.System() mols.add(multichain_cmap[0]) mols.add(multichain_cmap[1]) + mols.add(multichain_cmap[2]) # Sanity-check: at least two molecules must carry CMAP so that the # multi-chain code path is exercised. From 78480ce7d7ea9f981e87afa2ec861d7ac23bdb80 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 18 Mar 2026 13:59:37 +0000 Subject: [PATCH 8/9] Formatting tweak. --- tests/convert/test_openmm_cmap.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/tests/convert/test_openmm_cmap.py b/tests/convert/test_openmm_cmap.py index 2b9b8813d..7670728a4 100644 --- a/tests/convert/test_openmm_cmap.py +++ b/tests/convert/test_openmm_cmap.py @@ -139,11 +139,13 @@ def get_cmap_torsion_grids(context): return [] def unique_grids(torsion_grids, decimals=3): - """Return the sorted set of unique (size, rounded-grid) tuples. + """ + Return the sorted set of unique (size, rounded-grid) tuples. Torsion ordering can differ between the perturbable and non-perturbable code paths, so we compare the sets of unique grid shapes rather than - comparing torsion-by-torsion.""" + comparing torsion-by-torsion. + """ seen = set() result = [] for size, grid in torsion_grids: From 710561114ebc5fdc0710353956576ec6245fb6dd Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 18 Mar 2026 14:06:42 +0000 Subject: [PATCH 9/9] Return grid values by const reference. --- .../SireOpenMM/PerturbableOpenMMMolecule.pypp.cpp | 12 ++++++------ wrapper/Convert/SireOpenMM/openmmmolecule.cpp | 6 +++--- wrapper/Convert/SireOpenMM/openmmmolecule.h | 6 +++--- 3 files changed, 12 insertions(+), 12 deletions(-) diff --git a/wrapper/Convert/SireOpenMM/PerturbableOpenMMMolecule.pypp.cpp b/wrapper/Convert/SireOpenMM/PerturbableOpenMMMolecule.pypp.cpp index 807a90004..7622c75f3 100644 --- a/wrapper/Convert/SireOpenMM/PerturbableOpenMMMolecule.pypp.cpp +++ b/wrapper/Convert/SireOpenMM/PerturbableOpenMMMolecule.pypp.cpp @@ -587,39 +587,39 @@ void register_PerturbableOpenMMMolecule_class(){ } { //::SireOpenMM::PerturbableOpenMMMolecule::getCMAPGrids0 - typedef ::QVector< double > ( ::SireOpenMM::PerturbableOpenMMMolecule::*getCMAPGrids0_function_type)( ) const; + typedef ::QVector< double > const & ( ::SireOpenMM::PerturbableOpenMMMolecule::*getCMAPGrids0_function_type)( ) const; getCMAPGrids0_function_type getCMAPGrids0_function_value( &::SireOpenMM::PerturbableOpenMMMolecule::getCMAPGrids0 ); PerturbableOpenMMMolecule_exposer.def( "getCMAPGrids0" , getCMAPGrids0_function_value - , bp::release_gil_policy() + , bp::return_value_policy< bp::copy_const_reference >() , "Return the flat concatenated CMAP grid values (column-major, kJ/mol) " "for the reference state. Grid k has getCMAPGridSizes()[k]^2 entries." ); } { //::SireOpenMM::PerturbableOpenMMMolecule::getCMAPGrids1 - typedef ::QVector< double > ( ::SireOpenMM::PerturbableOpenMMMolecule::*getCMAPGrids1_function_type)( ) const; + typedef ::QVector< double > const & ( ::SireOpenMM::PerturbableOpenMMMolecule::*getCMAPGrids1_function_type)( ) const; getCMAPGrids1_function_type getCMAPGrids1_function_value( &::SireOpenMM::PerturbableOpenMMMolecule::getCMAPGrids1 ); PerturbableOpenMMMolecule_exposer.def( "getCMAPGrids1" , getCMAPGrids1_function_value - , bp::release_gil_policy() + , bp::return_value_policy< bp::copy_const_reference >() , "Return the flat concatenated CMAP grid values (column-major, kJ/mol) " "for the perturbed state. Grid k has getCMAPGridSizes()[k]^2 entries." ); } { //::SireOpenMM::PerturbableOpenMMMolecule::getCMAPGridSizes - typedef ::QVector< int > ( ::SireOpenMM::PerturbableOpenMMMolecule::*getCMAPGridSizes_function_type)( ) const; + typedef ::QVector< int > const & ( ::SireOpenMM::PerturbableOpenMMMolecule::*getCMAPGridSizes_function_type)( ) const; getCMAPGridSizes_function_type getCMAPGridSizes_function_value( &::SireOpenMM::PerturbableOpenMMMolecule::getCMAPGridSizes ); PerturbableOpenMMMolecule_exposer.def( "getCMAPGridSizes" , getCMAPGridSizes_function_value - , bp::release_gil_policy() + , bp::return_value_policy< bp::copy_const_reference >() , "Return the grid dimension N for each CMAP torsion (grid is N x N). " "Entries correspond to the grids in getCMAPGrids0/1." ); diff --git a/wrapper/Convert/SireOpenMM/openmmmolecule.cpp b/wrapper/Convert/SireOpenMM/openmmmolecule.cpp index e39bce583..015e61b1c 100644 --- a/wrapper/Convert/SireOpenMM/openmmmolecule.cpp +++ b/wrapper/Convert/SireOpenMM/openmmmolecule.cpp @@ -2966,19 +2966,19 @@ PerturbableOpenMMMolecule::getCMAPAtoms() const } /** Return flat concatenated CMAP grid values for state 0 (column-major, kJ/mol) */ -QVector PerturbableOpenMMMolecule::getCMAPGrids0() const +const QVector &PerturbableOpenMMMolecule::getCMAPGrids0() const { return cmap_grid0; } /** Return flat concatenated CMAP grid values for state 1 (column-major, kJ/mol) */ -QVector PerturbableOpenMMMolecule::getCMAPGrids1() const +const QVector &PerturbableOpenMMMolecule::getCMAPGrids1() const { return cmap_grid1; } /** Return the grid dimension N for each CMAP torsion (grid is N x N) */ -QVector PerturbableOpenMMMolecule::getCMAPGridSizes() const +const QVector &PerturbableOpenMMMolecule::getCMAPGridSizes() const { return cmap_grid_sizes; } diff --git a/wrapper/Convert/SireOpenMM/openmmmolecule.h b/wrapper/Convert/SireOpenMM/openmmmolecule.h index 5edc96090..8d46bf5ea 100644 --- a/wrapper/Convert/SireOpenMM/openmmmolecule.h +++ b/wrapper/Convert/SireOpenMM/openmmmolecule.h @@ -280,9 +280,9 @@ namespace SireOpenMM QVector getTorsionPhases0() const; QVector getTorsionPhases1() const; - QVector getCMAPGrids0() const; - QVector getCMAPGrids1() const; - QVector getCMAPGridSizes() const; + const QVector &getCMAPGrids0() const; + const QVector &getCMAPGrids1() const; + const QVector &getCMAPGridSizes() const; /** Return the 5-atom indices (molecule-local) for each CMAP torsion, * in the same order as getCMAPGridSizes(). Used for REST2 scaling. */