Add support for using a switching function with QM/MM

doc/source/changelog.rst

@@ -84,6 +84,8 @@ organisation on `GitHub <https://github.com/openbiosim/sire>`__.
   ``reversed.morph(force, lever, initial, final, λ)`` equals
   ``original.morph(force, lever, final, initial, 1-λ)``.
 
+* Add support for using a switching function for QM/MM simulations.
+
 `2025.4.0 <https://github.com/openbiosim/sire/compare/2025.3.0...2025.4.0>`__ - February 2026
 ---------------------------------------------------------------------------------------------
 

doc/source/tutorial/part08/01_intro.rst

@@ -37,6 +37,7 @@ an engine to perform the calculation:
 ...     cutoff="7.5A",
 ...     neighbour_list_frequency=20,
 ...     mechanical_embedding=False,
+...     switch_width=0.2,
 ... )
 
 Here the first argument is the molecules that we are simulating, the second
@@ -65,7 +66,7 @@ signature:
         xyz_mm: List[List[float]],
         cell: Optional[List[List[float]]] = None,
         idx_mm: Optional[List[int]] = None,
-    ) -> Tuple[float, List[List[float]], List[List[float]]]:
+    ) -> Tuple[float, List[List[float]], List[List[float]], Optional[List[float]]]:
 
 The function takes the atomic numbers of the QM atoms, the charges of the MM
 atoms in mod electron charge, the coordinates of the QM atoms in Angstrom, and
@@ -75,10 +76,18 @@ QM/MM region. This is useful for obtaining any additional atomic properties
 that may be required by the callback. (Note that link atoms and virtual charges
 are always placed last in the list of MM charges and positions.) The function
 should return the calculated energy in kJ/mol, the forces on the QM atoms in
-kJ/mol/nm, and the forces on the MM atoms in kJ/mol/nm. The remaining arguments
-are optional and specify the QM cutoff distance, the neighbour list update
-frequency, and whether the electrostatics should be treated with mechanical
-embedding. When mechanical embedding is used, the electrostatics are treated
+kJ/mol/nm, and the forces on the MM atoms in kJ/mol/nm. Optionally, it may also
+return a fourth element: the gradient of the energy with respect to the effective
+MM charges (dE/dq) in kJ/mol/e. When present, this is used to apply a chain-rule
+force correction that accounts for the positional dependence of the charge switching
+function (see below), which is important for systems with a charged ML region. The
+remaining arguments are optional and specify the QM cutoff distance, the neighbour
+list update frequency, whether the electrostatics should be treated with mechanical
+embedding, and the width of the switching region as a fraction of the cutoff
+(``switch_width``, default 0.2). The switching function smoothly scales the MM
+charges to zero between ``(1 - switch_width) * cutoff`` and ``cutoff``, which
+avoids force discontinuities as MM atoms enter or leave the cutoff sphere. Setting
+``switch_width=0`` disables switching (not recommended for charged ML regions). When mechanical embedding is used, the electrostatics are treated
 at the MM level by ``OpenMM``. Note that this doesn't change the signature of
 the callback function, i.e. it will be passed empty lists for the MM specific
 arguments and should return an empty list for the MM forces. Atomic positions

doc/source/tutorial/part08/02_emle.rst

@@ -55,14 +55,20 @@ an engine to perform the calculation:
 ...     mols[0],
 ...     calculator,
 ...     cutoff="7.5A",
-...     neighbour_list_frequency=20
+...     neighbour_list_frequency=20,
+...     switch_width=0.2,
 ... )
 
 Here the first argument is the molecules that we are simulating, the second
 selection coresponding to the QM region (here this is the first molecule), and
-the third is calculator that was created above. The fourth and fifth arguments
-are optional, and specify the QM cutoff distance and the neighbour list update
-frequency respectively. (Shown are the default values.) The function returns a
+the third is calculator that was created above. The remaining arguments are
+optional and specify the QM cutoff distance, the neighbour list update frequency,
+and the width of the switching region as a fraction of the cutoff (``switch_width``,
+default 0.2). (Shown are the default values.) The switching function smoothly
+scales the MM charges to zero between ``(1 - switch_width) * cutoff`` and
+``cutoff``, which avoids force discontinuities as MM atoms enter or leave the
+cutoff sphere. This is particularly important for systems with a charged ML region.
+Setting ``switch_width=0`` disables switching. The function returns a
 modified version of the molecules containing a "merged" dipeptide that can be
 interpolated between MM and QM levels of theory, along with an engine. The
 engine registers a Python callback that uses ``emle-engine`` to perform the QM
@@ -182,7 +188,7 @@ computes the electrostatic embedding:
 Next we create a new engine bound to the calculator:
 
 >>> _, engine = sr.qm.emle(
->>> ... mols, mols[0], calculator, cutoff="7.5A", neighbour_list_frequency=20
+>>> ... mols, mols[0], calculator, cutoff="7.5A", neighbour_list_frequency=20, switch_width=0.2
 >>> ... )
 
 .. note::
@@ -423,7 +429,7 @@ The model is serialisable, so can be saved and loaded using the standard
 It is also possible to use the model with Sire when performing QM/MM dynamics:
 
 >>> qm_mols, engine = sr.qm.emle(
-...     mols, mols[0], model, cutoff="7.5A", neighbour_list_frequency=20
+...     mols, mols[0], model, cutoff="7.5A", neighbour_list_frequency=20, switch_width=0.2
 ... )
 
 The model will be serialised and loaded into a C++ ``TorchQMEngine`` object,

pixi.toml

@@ -133,6 +133,9 @@ loguru = "*"
 pygit2 = "*"
 pyyaml = "*"
 
+[feature.emle.system-requirements]
+cuda = "12"
+
 [feature.emle.target.linux-64.dependencies]
 ambertools = ">=22"
 deepmd-kit = "*"

src/sire/qm/__init__.py

@@ -19,6 +19,7 @@ def create_engine(
     neighbour_list_frequency=0,
     mechanical_embedding=False,
     redistribute_charge=True,
+    switch_width=0.2,
     map=None,
 ):
     """
@@ -74,6 +75,13 @@ def create_engine(
         molecule, the excess charge is redistributed over the MM atoms within
         the residues of the QM region.
 
+    switch_width : float, optional, default=0.2
+        The width of the switching region as a fraction of the cutoff (0 to 1).
+        A quintic switching function is applied to the MM charges over the last
+        ``switch_width * cutoff`` angstroms before the cutoff, smoothly scaling
+        them to zero. Set to 0 or None to disable switching (not recommended
+        for production MD with charged ML regions).
+
     Returns
     -------
 
@@ -133,6 +141,15 @@ def create_engine(
     if not isinstance(redistribute_charge, bool):
         raise TypeError("'redistribute_charge' must be of type 'bool'")
 
+    if switch_width is None:
+        switch_width = 0.0
+    if not isinstance(switch_width, (int, float)):
+        raise TypeError("'switch_width' must be of type 'int' or 'float'")
+    switch_width = float(switch_width)
+    if switch_width < 0.0 or switch_width > 1.0:
+        raise ValueError("'switch_width' must be between 0 and 1")
+    use_switch = switch_width > 0.0
+
     if map is not None:
         if not isinstance(map, dict):
             raise TypeError("'map' must be of type 'dict'")
@@ -147,6 +164,10 @@ def create_engine(
         mechanical_embedding,
     )
 
+    # Configure the switching function.
+    engine.set_switch_width(switch_width)
+    engine.set_use_switch(use_switch)
+
     from ._utils import (
         _check_charge,
         _create_qm_mol_to_atoms,

src/sire/qm/_emle.py

@@ -108,6 +108,7 @@ def emle(
     cutoff="7.5A",
     neighbour_list_frequency=0,
     redistribute_charge=True,
+    switch_width=0.2,
     map=None,
 ):
     """
@@ -143,6 +144,13 @@ def emle(
         molecule, the excess charge is redistributed over the MM atoms within
         the residues of the QM region.
 
+    switch_width : float, optional, default=0.2
+        The width of the switching region as a fraction of the cutoff (0 to 1).
+        A quintic switching function is applied to the MM charges over the last
+        ``switch_width * cutoff`` angstroms before the cutoff, smoothly scaling
+        them to zero. Set to 0 or None to disable switching (not recommended
+        for production MD with charged ML regions).
+
     Returns
     -------
 
@@ -159,7 +167,6 @@ def emle(
 
     try:
         import torch as _torch
-        from emle.models import EMLE as _EMLE
 
         has_model = True
     except:
@@ -227,6 +234,15 @@ def emle(
     if not isinstance(redistribute_charge, bool):
         raise TypeError("'redistribute_charge' must be of type 'bool'")
 
+    if switch_width is None:
+        switch_width = 0.0
+    if not isinstance(switch_width, (int, float)):
+        raise TypeError("'switch_width' must be of type 'int' or 'float'")
+    switch_width = float(switch_width)
+    if switch_width < 0.0 or switch_width > 1.0:
+        raise ValueError("'switch_width' must be between 0 and 1")
+    use_switch = switch_width > 0.0
+
     if map is not None:
         if not isinstance(map, dict):
             raise TypeError("'map' must be of type 'dict'")
@@ -246,7 +262,7 @@ def emle(
     # Create an engine from an EMLE model.
     else:
         try:
-            from emle.models import EMLE as _EMLE
+            pass
         except:
             raise ImportError(
                 "Could not import emle.models. Please reinstall emle-engine and try again."
@@ -276,6 +292,10 @@ def emle(
         except Exception as e:
             raise ValueError("Unable to create a TorchEMLEEngine: " + str(e))
 
+    # Configure the switching function.
+    engine.set_switch_width(switch_width)
+    engine.set_use_switch(use_switch)
+
     from ._utils import (
         _check_charge,
         _create_qm_mol_to_atoms,

wrapper/Convert/SireOpenMM/PyQMCallback.pypp.cpp

@@ -2,8 +2,8 @@
 
 // (C) Christopher Woods, GPL >= 3 License
 
-#include "boost/python.hpp"
 #include "PyQMCallback.pypp.hpp"
+#include "boost/python.hpp"
 
 namespace bp = boost::python;
 
@@ -51,69 +51,56 @@ namespace bp = boost::python;
 
 #include <mutex>
 
-SireOpenMM::PyQMCallback __copy__(const SireOpenMM::PyQMCallback &other){ return SireOpenMM::PyQMCallback(other); }
+SireOpenMM::PyQMCallback __copy__(const SireOpenMM::PyQMCallback &other) { return SireOpenMM::PyQMCallback(other); }
 
 #include "Qt/qdatastream.hpp"
 
-const char* pvt_get_name(const SireOpenMM::PyQMCallback&){ return "SireOpenMM::PyQMCallback";}
+const char *pvt_get_name(const SireOpenMM::PyQMCallback &) { return "SireOpenMM::PyQMCallback"; }
 
 #include "Helpers/release_gil_policy.hpp"
 
-void register_PyQMCallback_class(){
+void register_PyQMCallback_class()
+{
 
     { //::SireOpenMM::PyQMCallback
-        typedef bp::class_< SireOpenMM::PyQMCallback > PyQMCallback_exposer_t;
-        PyQMCallback_exposer_t PyQMCallback_exposer = PyQMCallback_exposer_t( "PyQMCallback", "A callback wrapper class to interface with external QM engines\nvia the CustomCPPForceImpl.", bp::init< >("Default constructor.") );
-        bp::scope PyQMCallback_scope( PyQMCallback_exposer );
-        PyQMCallback_exposer.def( bp::init< bp::api::object, bp::optional< QString > >(( bp::arg("arg0"), bp::arg("name")="" ), "Constructor\nPar:am py_object\nA Python object that contains the callback function.\n\nPar:am name\nThe name of a callback method that take the following arguments:\n- numbers_qm: A list of atomic numbers for the atoms in the ML region.\n- charges_mm: A list of the MM charges in mod electron charge.\n- xyz_qm: A list of positions for the atoms in the ML region in Angstrom.\n- xyz_mm: A list of positions for the atoms in the MM region in Angstrom.\n- cell: A list of cell vectors in Angstrom.\n-  idx_mm: A list of indices for the MM atoms in the QM/MM region.\nThe callback should return a tuple containing:\n- The energy in kJmol.\n- A list of forces for the QM atoms in kJmolnm.\n- A list of forces for the MM atoms in kJmolnm.\nIf empty, then the object is assumed to be a callable.\n") );
+        typedef bp::class_<SireOpenMM::PyQMCallback> PyQMCallback_exposer_t;
+        PyQMCallback_exposer_t PyQMCallback_exposer = PyQMCallback_exposer_t("PyQMCallback", "A callback wrapper class to interface with external QM engines\nvia the CustomCPPForceImpl.", bp::init<>("Default constructor."));
+        bp::scope PyQMCallback_scope(PyQMCallback_exposer);
+        PyQMCallback_exposer.def(bp::init<bp::api::object, bp::optional<QString>>((bp::arg("arg0"), bp::arg("name") = ""), "Constructor\nPar:am py_object\nA Python object that contains the callback function.\n\nPar:am name\nThe name of a callback method that take the following arguments:\n- numbers_qm: A list of atomic numbers for the atoms in the ML region.\n- charges_mm: A list of the MM charges in mod electron charge.\n- xyz_qm: A list of positions for the atoms in the ML region in Angstrom.\n- xyz_mm: A list of positions for the atoms in the MM region in Angstrom.\n- cell: A list of cell vectors in Angstrom.\n-  idx_mm: A list of indices for the MM atoms in the QM/MM region.\nThe callback should return a tuple containing:\n- The energy in kJmol.\n- A list of forces for the QM atoms in kJmolnm.\n- A list of forces for the MM atoms in kJmolnm.\n- (Optional) The gradient of the energy w.r.t. the effective MM charges in kJmol/e, used for the chain-rule switching correction.\nIf empty, then the object is assumed to be a callable.\n"));
         { //::SireOpenMM::PyQMCallback::call
-
-            typedef ::boost::tuples::tuple< double, QVector< QVector< double > >, QVector< QVector< double > >, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type > ( ::SireOpenMM::PyQMCallback::*call_function_type)( ::QVector< int >,::QVector< double >,::QVector< QVector< double > >,::QVector< QVector< double > >,::QVector<QVector< double > >, ::QVector< int > ) const;
-            call_function_type call_function_value( &::SireOpenMM::PyQMCallback::call );
-
-            PyQMCallback_exposer.def( 
-                "call"
-                , call_function_value
-                , ( bp::arg("numbers_qm"), bp::arg("charges_mm"), bp::arg("xyz_qm"), bp::arg("xyz_mm"), bp::arg("cell"), bp::arg("idx_mm") )
-                , bp::release_gil_policy()
-                , "Call the callback function.\nPar:am numbers_qm\nA vector of atomic numbers for the atoms in the ML region.\n\nPar:am charges_mm\nA vector of the charges on the MM atoms in mod electron charge.\n\nPar:am xyz_qm\nA vector of positions for the atoms in the ML region in Angstrom.\n\nPar:am xyz_mm\nA vector of positions for the atoms in the MM region in Angstrom.\n\nPar:am cell A list of cell vectors in Angstrom.\n\nPar:am idx_mm A vector of indices for the MM atoms in the QM/MM region. Note that len(idx_mm) <= len(charges_mm) since it only contains the indices of true MM atoms, not link atoms or virtual charges.\n\nReturn:s\nA tuple containing:\n- The energy in kJmol.\n- A vector of forces for the QM atoms in kJmolnm.\n- A vector of forces for the MM atoms in kJmolnm.\n" );
-
+
+            typedef ::boost::tuples::tuple<double, QVector<QVector<double>>, QVector<QVector<double>>, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type> (::SireOpenMM::PyQMCallback::*call_function_type)(::QVector<int>, ::QVector<double>, ::QVector<QVector<double>>, ::QVector<QVector<double>>, ::QVector<QVector<double>>, ::QVector<int>) const;
+            call_function_type call_function_value(&::SireOpenMM::PyQMCallback::call);
+
+            PyQMCallback_exposer.def(
+                "call", call_function_value, (bp::arg("numbers_qm"), bp::arg("charges_mm"), bp::arg("xyz_qm"), bp::arg("xyz_mm"), bp::arg("cell"), bp::arg("idx_mm")), bp::release_gil_policy(), "Call the callback function.\nPar:am numbers_qm\nA vector of atomic numbers for the atoms in the ML region.\n\nPar:am charges_mm\nA vector of the charges on the MM atoms in mod electron charge.\n\nPar:am xyz_qm\nA vector of positions for the atoms in the ML region in Angstrom.\n\nPar:am xyz_mm\nA vector of positions for the atoms in the MM region in Angstrom.\n\nPar:am cell A list of cell vectors in Angstrom.\n\nPar:am idx_mm A vector of indices for the MM atoms in the QM/MM region. Note that len(idx_mm) <= len(charges_mm) since it only contains the indices of true MM atoms, not link atoms or virtual charges.\n\nReturn:s\nA tuple containing:\n- The energy in kJmol.\n- A vector of forces for the QM atoms in kJmolnm.\n- A vector of forces for the MM atoms in kJmolnm.\n- (Optional) The gradient of the energy w.r.t. the effective MM charges in kJmol/e, used for the chain-rule switching correction.\n");
         }
         { //::SireOpenMM::PyQMCallback::typeName
-
-            typedef char const * ( *typeName_function_type )(  );
-            typeName_function_type typeName_function_value( &::SireOpenMM::PyQMCallback::typeName );
-
-            PyQMCallback_exposer.def( 
-                "typeName"
-                , typeName_function_value
-                , bp::release_gil_policy()
-                , "Return the C++ name for this class." );
-
+
+            typedef char const *(*typeName_function_type)();
+            typeName_function_type typeName_function_value(&::SireOpenMM::PyQMCallback::typeName);
+
+            PyQMCallback_exposer.def(
+                "typeName", typeName_function_value, bp::release_gil_policy(), "Return the C++ name for this class.");
         }
         { //::SireOpenMM::PyQMCallback::what
-
-            typedef char const * ( ::SireOpenMM::PyQMCallback::*what_function_type)(  ) const;
-            what_function_type what_function_value( &::SireOpenMM::PyQMCallback::what );
-
-            PyQMCallback_exposer.def( 
-                "what"
-                , what_function_value
-                , bp::release_gil_policy()
-                , "Return the C++ name for this class." );
-
+
+            typedef char const *(::SireOpenMM::PyQMCallback::*what_function_type)() const;
+            what_function_type what_function_value(&::SireOpenMM::PyQMCallback::what);
+
+            PyQMCallback_exposer.def(
+                "what", what_function_value, bp::release_gil_policy(), "Return the C++ name for this class.");
         }
-        PyQMCallback_exposer.staticmethod( "typeName" );
-        PyQMCallback_exposer.def( "__copy__", &__copy__);
-        PyQMCallback_exposer.def( "__deepcopy__", &__copy__);
-        PyQMCallback_exposer.def( "clone", &__copy__);
-        PyQMCallback_exposer.def( "__rlshift__", &__rlshift__QDataStream< ::SireOpenMM::PyQMCallback >,
-                            bp::return_internal_reference<1, bp::with_custodian_and_ward<1,2> >() );
-        PyQMCallback_exposer.def( "__rrshift__", &__rrshift__QDataStream< ::SireOpenMM::PyQMCallback >,
-                            bp::return_internal_reference<1, bp::with_custodian_and_ward<1,2> >() );
-        PyQMCallback_exposer.def_pickle(sire_pickle_suite< ::SireOpenMM::PyQMCallback >());
-        PyQMCallback_exposer.def( "__str__", &pvt_get_name);
-        PyQMCallback_exposer.def( "__repr__", &pvt_get_name);
+        PyQMCallback_exposer.staticmethod("typeName");
+        PyQMCallback_exposer.def("__copy__", &__copy__);
+        PyQMCallback_exposer.def("__deepcopy__", &__copy__);
+        PyQMCallback_exposer.def("clone", &__copy__);
+        PyQMCallback_exposer.def("__rlshift__", &__rlshift__QDataStream<::SireOpenMM::PyQMCallback>,
+                                 bp::return_internal_reference<1, bp::with_custodian_and_ward<1, 2>>());
+        PyQMCallback_exposer.def("__rrshift__", &__rrshift__QDataStream<::SireOpenMM::PyQMCallback>,
+                                 bp::return_internal_reference<1, bp::with_custodian_and_ward<1, 2>>());
+        PyQMCallback_exposer.def_pickle(sire_pickle_suite<::SireOpenMM::PyQMCallback>());
+        PyQMCallback_exposer.def("__str__", &pvt_get_name);
+        PyQMCallback_exposer.def("__repr__", &pvt_get_name);
     }
-
 }