Add OpenMM CMAP support

corelib/src/libs/SireSystem/merge.cpp

@@ -131,6 +131,7 @@ namespace SireSystem
                 "atomtype",
                 "bond",
                 "charge",
+                "cmap",
                 "connectivity",
                 "coordinates",
                 "dihedral",

doc/source/changelog.rst

@@ -25,6 +25,10 @@ organisation on `GitHub <https://github.com/openbiosim/sire>`__.
 
 * Fixed parsing of AMBER and GROMACS GLYCAM force field topologies.
 
+* 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 <https://github.com/openbiosim/sire/compare/2025.3.0...2025.4.0>`__ - February 2026
 ---------------------------------------------------------------------------------------------
 

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,

tests/convert/test_openmm_cmap.py

@@ -0,0 +1,187 @@
+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])
+    mols.add(multichain_cmap[2])
+
+    # 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 an 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)
+
+
+@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. 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.
+    """
+    import openmm
+
+    platform_name = openmm_platform or "CPU"
+
+    mol0 = multichain_cmap[0]
+
+    mols_pert = sr.load_test_files("merged_molecule_cmap.s3")
+    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"
+        )

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 > const & ( ::SireOpenMM::PerturbableOpenMMMolecule::*getCMAPGrids0_function_type)(  ) const;
+            getCMAPGrids0_function_type getCMAPGrids0_function_value( &::SireOpenMM::PerturbableOpenMMMolecule::getCMAPGrids0 );
+
+            PerturbableOpenMMMolecule_exposer.def(
+                "getCMAPGrids0"
+                , getCMAPGrids0_function_value
+                , 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 > const & ( ::SireOpenMM::PerturbableOpenMMMolecule::*getCMAPGrids1_function_type)(  ) const;
+            getCMAPGrids1_function_type getCMAPGrids1_function_value( &::SireOpenMM::PerturbableOpenMMMolecule::getCMAPGrids1 );
+
+            PerturbableOpenMMMolecule_exposer.def(
+                "getCMAPGrids1"
+                , getCMAPGrids1_function_value
+                , 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 > const & ( ::SireOpenMM::PerturbableOpenMMMolecule::*getCMAPGridSizes_function_type)(  ) const;
+            getCMAPGridSizes_function_type getCMAPGridSizes_function_value( &::SireOpenMM::PerturbableOpenMMMolecule::getCMAPGridSizes );
+
+            PerturbableOpenMMMolecule_exposer.def(
+                "getCMAPGridSizes"
+                , getCMAPGridSizes_function_value
+                , 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." );
+
+        }
+        { //::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 );
 

wrapper/Convert/SireOpenMM/lambdalever.cpp

@@ -1107,6 +1107,27 @@ PropertyList LambdaLever::getLeverValues(const QVector<double> &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<double>());
+            }
+        }
+
+        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<OpenMM::HarmonicBondForce>("bond", system);
     auto angff = this->getForce<OpenMM::HarmonicAngleForce>("angle", system);
     auto dihff = this->getForce<OpenMM::PeriodicTorsionForce>("torsion", system);
+    auto cmapff = this->getForce<OpenMM::CMAPTorsionForce>("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<double> 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