Fix issue #506

doc/source/changelog.rst

@@ -87,6 +87,9 @@ organisation on `GitHub <https://github.com/openbiosim/sire>`__.
 * Fixed ``sire.restraints.boresch_search()`` failing for charged ligands, since alchemical
   ions used to keep the charge constant were counted as additional perturbable molecules.
 
+* Fixed MM non-bonded interactions not being removed between QM atoms in different
+  molecules, e.g. when a QM water reacts with a QM ligand.
+
 `2026.2.0 <https://github.com/openbiosim/sire/compare/2026.1.0...2026.2.0>`__ - September 2026
 ----------------------------------------------------------------------------------------------
 

tests/qm/test_qm.py

@@ -660,6 +660,76 @@ def callback(numbers_qm, charges_mm, xyz_qm, xyz_mm, cell=None, idx_mm=None):
     assert lever.was_force_changed("qmff"), "qmff should be changed when lambda changes"
 
 
+def test_qm_cross_molecule_exceptions(ala_mols):
+    """
+    Make sure that the MM non-bonded interactions between QM atoms in
+    different molecules are scaled from fully on at lambda = 0 to fully
+    off at lambda = 1.
+    """
+
+    import math
+    import openmm
+
+    def callback(numbers_qm, charges_mm, xyz_qm, xyz_mm, cell=None, idx_mm=None):
+        return (0.0, xyz_qm, xyz_mm)
+
+    mols = ala_mols.clone()
+    qm_mols, engine = sr.qm.create_engine(
+        mols, mols["molidx 0 or molidx 1"].atoms(), callback, callback=None
+    )
+
+    d = qm_mols.dynamics(
+        timestep="1fs", constraint="none", qm_engine=engine, platform="cpu"
+    )
+
+    ctx = d._d._omm_mols
+    system = ctx.getSystem()
+
+    cljff = [f for f in system.getForces() if isinstance(f, openmm.NonbondedForce)][0]
+
+    num_atoms0 = qm_mols[0].num_atoms()
+    num_atoms1 = qm_mols[1].num_atoms()
+
+    cross_pairs = {
+        (i, j)
+        for i in range(num_atoms0)
+        for j in range(num_atoms0, num_atoms0 + num_atoms1)
+    }
+
+    def get_exceptions():
+        exceptions = {}
+        for k in range(cljff.getNumExceptions()):
+            i, j, q, sigma, epsilon = cljff.getExceptionParameters(k)
+            pair = (min(i, j), max(i, j))
+            if pair in cross_pairs:
+                exceptions[pair] = (
+                    q.value_in_unit(openmm.unit.elementary_charge**2),
+                    sigma.value_in_unit(openmm.unit.nanometer),
+                    epsilon.value_in_unit(openmm.unit.kilojoule_per_mole),
+                )
+        return exceptions
+
+    # The context is initialised at lambda = 1, so the pairs are switched off.
+    exceptions = get_exceptions()
+    assert set(exceptions.keys()) == cross_pairs
+    for q, sigma, epsilon in exceptions.values():
+        assert q == pytest.approx(0.0)
+        assert epsilon == pytest.approx(0.0, abs=1e-8)
+
+    # At lambda = 0 the pairs should have the full MM interaction.
+    ctx.set_lambda(0.0)
+    exceptions = get_exceptions()
+    for (i, j), (q, sigma, epsilon) in exceptions.items():
+        q0, sig0, eps0 = cljff.getParticleParameters(i)
+        q1, sig1, eps1 = cljff.getParticleParameters(j)
+        q0 = q0.value_in_unit(openmm.unit.elementary_charge)
+        q1 = q1.value_in_unit(openmm.unit.elementary_charge)
+        eps0 = eps0.value_in_unit(openmm.unit.kilojoule_per_mole)
+        eps1 = eps1.value_in_unit(openmm.unit.kilojoule_per_mole)
+        assert q == pytest.approx(q0 * q1)
+        assert epsilon == pytest.approx(math.sqrt(eps0 * eps1), abs=1e-8)
+
+
 def test_property_bases(ala_mols):
     """
     Make sure that the SireOpenMM property classes are exposed as

wrapper/Convert/SireOpenMM/lambdalever.cpp

@@ -279,7 +279,8 @@ LambdaLever::LambdaLever(const LambdaLever &other)
       perturbable_maps(other.perturbable_maps),
       lambda_cache(other.lambda_cache),
       last_rest2_scale(-1.0),
-      last_qmff_lam(-1.0)
+      last_qmff_lam(-1.0),
+      qm_exception_idxs(other.qm_exception_idxs)
 {
 }
 
@@ -299,6 +300,7 @@ LambdaLever &LambdaLever::operator=(const LambdaLever &other)
         start_indices = other.start_indices;
         perturbable_maps = other.perturbable_maps;
         lambda_cache = other.lambda_cache;
+        qm_exception_idxs = other.qm_exception_idxs;
         Property::operator=(other);
     }
 
@@ -312,7 +314,8 @@ bool LambdaLever::operator==(const LambdaLever &other) const
            lambda_schedule == other.lambda_schedule and
            perturbable_mols == other.perturbable_mols and
            start_indices == other.start_indices and
-           perturbable_maps == other.perturbable_maps;
+           perturbable_maps == other.perturbable_maps and
+           qm_exception_idxs == other.qm_exception_idxs;
 }
 
 bool LambdaLever::operator!=(const LambdaLever &other) const
@@ -2047,6 +2050,47 @@ double LambdaLever::setLambda(OpenMM::Context &context,
         }
     }
 
+    // scale the MM interactions between QM atoms in different molecules
+    // from fully on (MM, lambda = 0) to fully off (QM, lambda = 1)
+    if (cljff != nullptr and not qm_exception_idxs.isEmpty())
+    {
+        const double coul_scale = this->lambda_schedule.morph(
+            "clj", "charge_scale", 1.0, 0.0, lambda_value);
+        const double lj_scale = this->lambda_schedule.morph(
+            "clj", "lj_scale", 1.0, 0.0, lambda_value);
+
+        for (const auto &idx : qm_exception_idxs)
+        {
+            int atom0, atom1;
+            double old_charge, old_sigma, old_epsilon;
+            cljff->getExceptionParameters(idx, atom0, atom1,
+                                          old_charge, old_sigma, old_epsilon);
+
+            double q0, sig0, eps0;
+            double q1, sig1, eps1;
+            cljff->getParticleParameters(atom0, q0, sig0, eps0);
+            cljff->getParticleParameters(atom1, q1, sig1, eps1);
+
+            const double charge = coul_scale * q0 * q1;
+            double sigma = 0.5 * (sig0 + sig1);
+            double epsilon = lj_scale * std::sqrt(eps0 * eps1);
+
+            // stop OpenMM from turning this exception into an exclusion
+            if (charge == 0 and epsilon == 0)
+            {
+                sigma = 1e-9;
+                epsilon = 1e-9;
+            }
+
+            if (charge != old_charge or sigma != old_sigma or epsilon != old_epsilon)
+            {
+                cljff->setExceptionParameters(idx, atom0, atom1,
+                                              charge, sigma, epsilon);
+                has_changed_cljff = true;
+            }
+        }
+    }
+
     // update the parameters in the context for forces whose parameters changed
     if (has_changed_cljff and cljff)
         cljff->updateParametersInContext(context);
@@ -2683,6 +2727,14 @@ void LambdaLever::setConstraintIndicies(int mol_idx, const QVector<qint32> &cons
     this->perturbable_mols[mol_idx].setConstraintIndicies(constraint_idxs);
 }
 
+/** Set the indices of the NonbondedForce exceptions between QM atoms
+ *  in different molecules
+ */
+void LambdaLever::setQMExceptionIndicies(const QVector<qint32> &exception_idxs)
+{
+    this->qm_exception_idxs = exception_idxs;
+}
+
 /** Return all of the property maps used to find the perturbable properties
  *  of the perturbable molecules. This is indexed by molecule number
  */

wrapper/Convert/SireOpenMM/lambdalever.h

@@ -139,6 +139,8 @@ namespace SireOpenMM
 
         void setConstraintIndicies(int idx, const QVector<qint32> &constraint_idxs);
 
+        void setQMExceptionIndicies(const QVector<qint32> &exception_idxs);
+
         void setSchedule(const SireCAS::LambdaSchedule &schedule);
 
         void addLever(const QString &lever_name);
@@ -231,6 +233,10 @@ namespace SireOpenMM
          *  non-empty these atoms are excluded from background-lrc (their LRC
          *  is handled by the gcmc-lrc CustomVolumeForce instead). */
         QSet<int> gcmc_water_atoms;
+
+        /** Indices of the NonbondedForce exceptions between QM atoms
+         *  in different molecules */
+        QVector<qint32> qm_exception_idxs;
     };
 
 #ifndef SIRE_SKIP_INLINE_FUNCTION

wrapper/Convert/SireOpenMM/qmmm.h

@@ -32,6 +32,8 @@
 #include "OpenMM.h"
 #include "openmm/Force.h"
 
+#include <QVector>
+
 #include "sireglobal.h"
 
 #include "SireBase/property.h"
@@ -54,6 +56,9 @@ namespace SireOpenMM
         //! Set the lambda weighting factor.
         virtual void setLambda(double lambda) = 0;
 
+        //! Get the indices of the QM atoms.
+        virtual QVector<int> getAtoms() const = 0;
+
     protected:
         virtual OpenMM::ForceImpl *createImpl() const = 0;
     };
@@ -74,7 +79,7 @@ namespace SireOpenMM
         }
 
         //! Create a QM force object.
-        virtual QMForce* createForce() const=0;
+        virtual QMForce *createForce() const = 0;
 
         //! Get a null QM engine.
         static const NullQMEngine &null();
@@ -96,7 +101,7 @@ namespace SireOpenMM
         const char *what() const;
 
         //! Create a QM force object.
-        QMForce* createForce() const;
+        QMForce *createForce() const;
     };
 }
 

wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp

@@ -3034,6 +3034,67 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system,
         }
     }
 
+    // add exceptions for QM atoms in different molecules, since these can't
+    // be removed via each molecule's intrascale. These are scaled by the lever.
+    if (qmff != 0 and cljff != 0)
+    {
+        const auto qm_atoms = qmff->getAtoms();
+
+        QVector<int> qm_mol_idxs(qm_atoms.count());
+
+        for (int i = 0; i < qm_atoms.count(); ++i)
+        {
+            int mol_idx = 0;
+
+            while (mol_idx + 1 < nmols and start_indexes[mol_idx + 1] <= qm_atoms[i])
+                mol_idx++;
+
+            qm_mol_idxs[i] = mol_idx;
+        }
+
+        QVector<qint32> qm_exception_idxs;
+
+        for (int i = 0; i < qm_atoms.count(); ++i)
+        {
+            for (int j = i + 1; j < qm_atoms.count(); ++j)
+            {
+                if (qm_mol_idxs[i] == qm_mol_idxs[j])
+                    continue;
+
+                const int atom0 = qm_atoms[i];
+                const int atom1 = qm_atoms[j];
+
+                double q0, sig0, eps0;
+                double q1, sig1, eps1;
+                cljff->getParticleParameters(atom0, q0, sig0, eps0);
+                cljff->getParticleParameters(atom1, q1, sig1, eps1);
+
+                const double charge = q0 * q1;
+                double sigma = 0.5 * (sig0 + sig1);
+                double epsilon = std::sqrt(eps0 * eps1);
+
+                // stop OpenMM from turning this exception into an exclusion
+                if (charge == 0 and epsilon == 0)
+                {
+                    sigma = 1e-9;
+                    epsilon = 1e-9;
+                }
+
+                qm_exception_idxs.append(cljff->addException(atom0, atom1,
+                                                             charge, sigma,
+                                                             epsilon, false));
+
+                if (ghost_ghostff != 0)
+                {
+                    ghost_ghostff->addExclusion(atom0, atom1);
+                    ghost_nonghostff->addExclusion(atom0, atom1);
+                }
+            }
+        }
+
+        lambda_lever.setQMExceptionIndicies(qm_exception_idxs);
+    }
+
     // Stage 6 is complete. We have set up all of the exceptions. The
     // total energy / force calculated for the system should now be
     // correct.