Add softcore ring-breaking nonbonded force

doc/source/changelog.rst

@@ -86,6 +86,8 @@ organisation on `GitHub <https://github.com/openbiosim/sire>`__.
 
 * Add support for using a switching function for QM/MM simulations.
 
+* Add softcore ``CustomBondForce`` for ring-breaking and ring-making pairs.
+
 `2025.4.0 <https://github.com/openbiosim/sire/compare/2025.3.0...2025.4.0>`__ - February 2026
 ---------------------------------------------------------------------------------------------
 

wrapper/Convert/SireOpenMM/lambdalever.cpp

@@ -1299,6 +1299,8 @@ double LambdaLever::setLambda(OpenMM::Context &context,
     auto ghost_ghostff = this->getForce<OpenMM::CustomNonbondedForce>("ghost/ghost", system);
     auto ghost_nonghostff = this->getForce<OpenMM::CustomNonbondedForce>("ghost/non-ghost", system);
     auto ghost_14ff = this->getForce<OpenMM::CustomBondForce>("ghost-14", system);
+    auto ring_breaking_ff = this->getForce<OpenMM::CustomBondForce>("ring-break", system);
+    auto ring_making_ff = this->getForce<OpenMM::CustomBondForce>("ring-make", system);
     auto bondff = this->getForce<OpenMM::HarmonicBondForce>("bond", system);
     auto angff = this->getForce<OpenMM::HarmonicAngleForce>("angle", system);
     auto dihff = this->getForce<OpenMM::PeriodicTorsionForce>("torsion", system);
@@ -1330,6 +1332,25 @@ double LambdaLever::setLambda(OpenMM::Context &context,
     bool has_changed_dihff = false;
     bool has_changed_cmap = false;
 
+    // Pre-compute ring-break/make kappa values so the per-mol exception update
+    // loop and the later global-parameter block both use the same values.
+    const double rb_kappa = (ring_breaking_ff != nullptr)
+                                ? std::max(0.0, std::min(1.0, this->lambda_schedule.morph(
+                                                                  "ring-break", "kappa", 0.0, 1.0, lambda_value)))
+                                : 0.0;
+    const double rb_alpha = (ring_breaking_ff != nullptr)
+                                ? std::max(0.0, std::min(1.0, this->lambda_schedule.morph(
+                                                                  "ring-break", "alpha", 1.0, 0.0, lambda_value)))
+                                : 1.0;
+    const double rm_kappa = (ring_making_ff != nullptr)
+                                ? std::max(0.0, std::min(1.0, this->lambda_schedule.morph(
+                                                                  "ring-make", "kappa", 1.0, 0.0, lambda_value)))
+                                : 1.0;
+    const double rm_alpha = (ring_making_ff != nullptr)
+                                ? std::max(0.0, std::min(1.0, this->lambda_schedule.morph(
+                                                                  "ring-make", "alpha", 0.0, 1.0, lambda_value)))
+                                : 0.0;
+
     // change the parameters for all of the perturbable molecules
     for (int i = 0; i < this->perturbable_mols.count(); ++i)
     {
@@ -1658,6 +1679,12 @@ double LambdaLever::setLambda(OpenMM::Context &context,
                         scale = rest2_scale;
                     }
 
+                    // ring-breaking/making pairs have idx=-1 (their CLJ
+                    // exception is fixed at 1e-9 and must not be updated;
+                    // the ring force handles the interaction via global params)
+                    if (boost::get<0>(idxs[j]) == -1)
+                        continue;
+
                     // don't set LJ terms for ghost atoms
                     if (atom0_is_ghost or atom1_is_ghost)
                     {
@@ -1725,6 +1752,55 @@ double LambdaLever::setLambda(OpenMM::Context &context,
             }
         }
 
+        // Update CLJ exceptions for ring-breaking pairs: the exception charge is
+        // rb_kappa * q_a0(λ) * q_a1(λ), where q_a0/q_a1 are read from the CLJ
+        // particle parameters that were already morphed (and REST2-scaled via
+        // sqrt_scale) earlier in this function. LJ stays at 1e-9; the ring-break
+        //  CustomBondForce provides the full LJ correction.
+        if (cljff != nullptr and ring_breaking_ff != nullptr)
+        {
+            const auto rb_exc = perturbable_mol.getExceptionIndicies("ring-break");
+            for (int j = 0; j < rb_exc.count(); ++j)
+            {
+                const int clj_idx = boost::get<0>(rb_exc[j]);
+                if (clj_idx < 0)
+                    continue;
+                const int bond_idx = boost::get<1>(rb_exc[j]);
+                int a0, a1;
+                std::vector<double> bp;
+                ring_breaking_ff->getBondParameters(bond_idx, a0, a1, bp);
+                double q_a0, sig_a0, eps_a0;
+                double q_a1, sig_a1, eps_a1;
+                cljff->getParticleParameters(a0, q_a0, sig_a0, eps_a0);
+                cljff->getParticleParameters(a1, q_a1, sig_a1, eps_a1);
+                cljff->setExceptionParameters(clj_idx, a0, a1,
+                                              rb_kappa * q_a0 * q_a1, 1e-9, 1e-9);
+                has_changed_cljff = true;
+            }
+        }
+
+        if (cljff != nullptr and ring_making_ff != nullptr)
+        {
+            const auto rm_exc = perturbable_mol.getExceptionIndicies("ring-make");
+            for (int j = 0; j < rm_exc.count(); ++j)
+            {
+                const int clj_idx = boost::get<0>(rm_exc[j]);
+                if (clj_idx < 0)
+                    continue;
+                const int bond_idx = boost::get<1>(rm_exc[j]);
+                int a0, a1;
+                std::vector<double> bp;
+                ring_making_ff->getBondParameters(bond_idx, a0, a1, bp);
+                double q_a0, sig_a0, eps_a0;
+                double q_a1, sig_a1, eps_a1;
+                cljff->getParticleParameters(a0, q_a0, sig_a0, eps_a0);
+                cljff->getParticleParameters(a1, q_a1, sig_a1, eps_a1);
+                cljff->setExceptionParameters(clj_idx, a0, a1,
+                                              rm_kappa * q_a0 * q_a1, 1e-9, 1e-9);
+                has_changed_cljff = true;
+            }
+        }
+
         // update all of the perturbable constraints
         if (update_constraints)
         {
@@ -2054,6 +2130,28 @@ double LambdaLever::setLambda(OpenMM::Context &context,
         context.setParameter("lrc_scale", lrc_scale);
     }
 
+    // Update ring-breaking/making softcore force global parameters using the
+    // values pre-computed before the per-mol loop (rb_alpha/kappa, rm_alpha/kappa).
+    bool has_changed_ring_breaking_ff = false;
+    if (ring_breaking_ff != nullptr)
+    {
+        has_changed_ring_breaking_ff =
+            (rb_alpha != context.getParameter("ring_break_alpha") ||
+             rb_kappa != context.getParameter("ring_break_kappa"));
+        context.setParameter("ring_break_alpha", rb_alpha);
+        context.setParameter("ring_break_kappa", rb_kappa);
+    }
+
+    bool has_changed_ring_making_ff = false;
+    if (ring_making_ff != nullptr)
+    {
+        has_changed_ring_making_ff =
+            (rm_alpha != context.getParameter("ring_make_alpha") ||
+             rm_kappa != context.getParameter("ring_make_kappa"));
+        context.setParameter("ring_make_alpha", rm_alpha);
+        context.setParameter("ring_make_kappa", rm_kappa);
+    }
+
     // Update the NonbondedForce (background) LRC via its own CustomVolumeForce.
     // Ghost atoms have epsilon=0 in cljff so they contribute nothing naturally.
     auto background_lrc_ff = this->getForce<OpenMM::CustomVolumeForce>("background-lrc", system);
@@ -2179,6 +2277,8 @@ double LambdaLever::setLambda(OpenMM::Context &context,
     last_changed_forces["background-lrc"] = has_changed_cljff;
     last_changed_forces["gcmc-lrc"] = false;
     last_changed_forces["ghost-14"] = has_changed_ghost14ff;
+    last_changed_forces["ring-break"] = has_changed_ring_breaking_ff;
+    last_changed_forces["ring-make"] = has_changed_ring_making_ff;
     last_changed_forces["bond"] = has_changed_bondff;
     last_changed_forces["angle"] = has_changed_angff;
     last_changed_forces["torsion"] = has_changed_dihff;

wrapper/Convert/SireOpenMM/openmmmolecule.cpp

@@ -294,6 +294,36 @@ OpenMMMolecule::OpenMMMolecule(const Molecule &mol,
                 std::swap(map0, map1);
             }
 
+            // Read ring-breaking/making bond pairs from molecule properties,
+            // swapping them if end states are inverted so that the members
+            // always reflect the λ=0/λ=1 convention of the (possibly swapped)
+            // end states.
+            auto read_ring_pairs = [&](const QString &prop_name)
+            {
+                QVector<QPair<qint32, qint32>> pairs;
+                if (mol.hasProperty(prop_name))
+                {
+                    const auto &flat = mol.property(prop_name)
+                                           .asA<SireBase::IntegerArrayProperty>()
+                                           .toVector();
+                    pairs.reserve(flat.count() / 2);
+                    for (int k = 0; k + 1 < flat.count(); k += 2)
+                        pairs.append(QPair<qint32, qint32>(flat[k], flat[k + 1]));
+                }
+                return pairs;
+            };
+
+            if (swap_end_states)
+            {
+                this->ring_breaking_pairs = read_ring_pairs("ring_making_bonds");
+                this->ring_making_pairs = read_ring_pairs("ring_breaking_bonds");
+            }
+            else
+            {
+                this->ring_breaking_pairs = read_ring_pairs("ring_breaking_bonds");
+                this->ring_making_pairs = read_ring_pairs("ring_making_bonds");
+            }
+
             // save this perturbable map - this will help us set
             // new properties from the results of dynamics, e.g.
             // updating coordinates after minimisation

wrapper/Convert/SireOpenMM/openmmmolecule.h

@@ -3,17 +3,17 @@
 
 #include <OpenMM.h>
 
-#include "SireMol/moleculeinfo.h"
-#include "SireMol/core.h"
 #include "SireMol/atom.h"
+#include "SireMol/core.h"
+#include "SireMol/moleculeinfo.h"
 #include "SireMol/selector.hpp"
 
-#include "SireMM/mmdetail.h"
-#include "SireMM/excludedpairs.h"
 #include "SireMM/amberparams.h"
-#include "SireMM/bond.h"
 #include "SireMM/angle.h"
+#include "SireMM/bond.h"
 #include "SireMM/dihedral.h"
+#include "SireMM/excludedpairs.h"
+#include "SireMM/mmdetail.h"
 
 #include <boost/tuple/tuple.hpp>
 
@@ -34,7 +34,7 @@ namespace SireOpenMM
         };
 
         Type type;
-        int vsite_idx;            // molecule-local index of the virtual site atom
+        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
@@ -208,6 +208,18 @@ namespace SireOpenMM
          */
         QSet<qint32> from_ghost_idxs;
 
+        /** Pairs of atom indices (molecule-local) whose bond is present
+         *  at λ=0 but absent at λ=1 (ring-breaking). Already swapped if
+         *  swap_end_states was active at construction time.
+         */
+        QVector<QPair<qint32, qint32>> ring_breaking_pairs;
+
+        /** Pairs of atom indices (molecule-local) whose bond is absent
+         *  at λ=0 but present at λ=1 (ring-making). Already swapped if
+         *  swap_end_states was active at construction time.
+         */
+        QVector<QPair<qint32, qint32>> ring_making_pairs;
+
         /** What type of constraint to use */
         qint32 constraint_type;