Fix issue #262

doc/source/changelog.rst

@@ -21,6 +21,13 @@ organisation on `GitHub <https://github.com/openbiosim/sire>`__.
 * Expose missing ``timeout`` kwarg in :meth:`dynamics.minimise()` method.
 * Expose missing ``include_constrained_energies`` kwarg in minimisation function.
 * Make minimisation function settings consistent across API.
+* Reload TorchQMForce module if OpenMM platform changes.
+* Fix calculation of non-equilibrium work when starting from QM state.
+* Fix stereochemistry in RDKit molecule conversion.
+* Fixed :func:`sire.restraints.get_standard_state_correction` to be consistent with the definition of
+  the "force constanst" as ``F = 2 ** k ** x`` (rather than ``F = k ** x``). Updated docstrings and
+  restraints documentation to make it clear how the force constants are defined.
+* Fix thread safety issue in Sire OpenMM minimiser.
 
 `2024.3.0 <https://github.com/openbiosim/sire/compare/2024.2.0...2024.3.0>`__ - October 2024
 --------------------------------------------------------------------------------------------

doc/source/tutorial/part06/03_restraints.rst

@@ -13,6 +13,13 @@ in the :mod:`sire.restraints` module. These functions return
 :class:`~sire.mm.Restraints` objects that contain all of the information that
 needs to be passed to OpenMM to add the restraints to the system.
 
+.. note::
+
+   For all harmonic restraints, the restraint energy is defined as ``E = k * x ** 2``,
+   and not as ``E = 0.5 * k * x ** 2``. Hence, the force is ``F = 2 * k * x`` and all
+   ``k`` values supplied to the restraints functions are half the value of the force
+   constants.
+
 Positional Restraints
 ---------------------
 
@@ -65,11 +72,11 @@ PositionalRestraints( size=3
 2: PositionalRestraint( 14 => ( 15.3698, 4.19397, 16.434 ), k=150 kcal mol-1 Å-2 : r0=0 Å )
 )
 
-The default force constant is 150 kcal mol-1 Å-2. By default, the
+The default half force constant is 150 kcal mol-1 Å-2. By default, the
 atoms are restrained to their current positions, and are held exactly in
 those positions (hence why the ``r0=0 Å`` in the output above).
 
-You can set the force constant and width of the half-harmonic restraint
+You can set the half force constant and width of the half-harmonic restraint
 used to hold the atoms in position using the ``k`` and ``r0`` arguments, e.g.
 
 >>> restraints = sr.restraints.positional(mols, atoms="resname ALA and element C",
@@ -177,7 +184,7 @@ BondRestraints( size=1
 
 creates a single harmonic distance (or bond) restraint that acts between
 atoms 0 and 1. By default, the equilibrium distance (r0) is the current
-distance between the atoms (1.09 Å), and the force constant (k) is
+distance between the atoms (1.09 Å), and the half force constant (k) is
 150 kcal mol-1 Å-2.
 
 You can set these via the ``k`` and ``r0`` arguments, e.g.
@@ -243,23 +250,23 @@ in the ligand. For more detail, please see J. Phys. Chem. B 2003, 107, 35, 9535
 To create a Boresch restraint, you need to specify the receptor and ligand anchor
 atoms (note that the order of the atoms is important). Like the distance restraints,
 the atoms can be specified using a search string, passing lists of atom indexes, or 
-molecule views holding the atoms. You can also specify the force constants and equilibrium
-values for the restraints. If not supplied, default force constants of 10 kcal mol-1 Å-2
-and 100 kcal mol-1 rad-2 are used for the distance and angle restraints, respectively,
+molecule views holding the atoms. You can also specify the half force constants and equilibrium
+values for the restraints. If not supplied, default half force constants of 5 kcal mol-1 Å-2
+and 50 kcal mol-1 rad-2 are used for the distance and angle restraints, respectively,
 and the equilibrium values are set to the current values of the distances and angles in
 the system supplied. For example,
 
 >>> boresch_restraints = sr.restraints.boresch(mols, receptor=[1574, 1554, 1576], ligand=[4,3,5])
 >>> boresch_restraint = boresch_restraints[0]
 >>> print(boresch_restraint)
 BoreschRestraint( [1574, 1554, 1576] => [4, 3, 5],
-                  k=[10 kcal mol-1 Å-2, 0.0304617 kcal mol-1 °-2, 0.0304617 kcal mol-1 °-2,
-                   0.0304617 kcal mol-1 °-2, 0.0304617 kcal mol-1 °-2, 0.0304617 kcal mol-1 °-2]
+                  k=[5 kcal mol-1 Å-2, 0.0152309 kcal mol-1 °-2, 0.0152309 kcal mol-1 °-2, 
+                  0.0152309 kcal mol-1 °-2, 0.0152309 kcal mol-1 °-2, 0.0152309 kcal mol-1 °-2]
                   r0=15.1197 Å, theta0=[80.5212°, 59.818°],
                   phi0=[170.562°Ⱐ128.435°Ⱐ192.21°] )
 
 creates a Boresch restraint where the receptor anchor atoms are r1 = 1574, r2 = 1554, and r3 = 1576,
-and the ligand anchor atoms are l1 = 4, l2 = 3, and l3 = 5. The default force constants have been set
+and the ligand anchor atoms are l1 = 4, l2 = 3, and l3 = 5. The default half force constants have been set
 and the equilibrium values have been set to the current values of the distances and angles in the
 system supplied.
 
@@ -268,10 +275,10 @@ system supplied.
    Boresch restraints can be subject to instabilities if any three contiguous anchor points
    approach collinearity (J. Chem. Theory Comput. 2023, 19, 12, 3686–3704). It is important to
    prevent this by ensuring the associated angles are sufficiently far from 0 or 180 degrees,
-   and that the `ktheta` force constants are high enough. Sire will raise a warning if the 
-   `theta0` values are too close to 0 or 180 degrees for the given temperature and force constants.
+   and that the `ktheta` half force constants are high enough. Sire will raise a warning if the 
+   `theta0` values are too close to 0 or 180 degrees for the given temperature and half force constants.
 
-Alternatively, we could have explicitly set the force constants and equilibrium values, e.g.
+Alternatively, we could have explicitly set the half force constants and equilibrium values, e.g.
 
 >>> boresch_restraints = sr.restraints.boresch(mols, receptor=[1574, 1554, 1576], ligand=[4,3,5],
                                                kr = "6.2012 kcal mol-1 A-2",

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

@@ -112,8 +112,12 @@ QM engine when creating a dynamics object, for example:
 ...     platform="cpu",
 ... )
 
-For QM/MM simulations it is recommended to use a 1 femtosecond timestep and no
-constraints. The simulation can then be run as usual:
+Next we will minimise the system ready for simulation:
+
+>>> d.minimise()
+
+When runinng QM/MM simulations it is recommended to use a 1 femtosecond timestep
+and no constraints. The simulation can then be run as usual:
 
 >>> d.run("100ps", energy_frequency="1ps", frame_frequency="1ps")
 

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

@@ -86,6 +86,10 @@ energy calculations.
 ...     platform="cpu",
 ... )
 
+Next we will minimise the system ready for simulation:
+
+>>> d.minimise()
+
 We can now run the simulation. The options below specify the run time, the
 frequency at which trajectory frames are saved, and the frequency at which
 energies are recorded. The ``energy_frequency`` also specifies the frequency

src/sire/mol/_dynamics.py

@@ -328,7 +328,10 @@ def _exit_dynamics_block(
             if self._is_interpolate:
                 nrg = nrgs["potential"]
 
-                if sim_lambda_value != 0.0:
+                if (
+                    len(self._lambda_interpolate) == 2
+                    and sim_lambda_value != self._lambda_interpolate[0]
+                ):
                     self._work += delta_lambda * (nrg - self._nrg_prev)
                 self._nrg_prev = nrg
                 nrgs["work"] = self._work

src/sire/restraints/_restraints.py

@@ -36,7 +36,9 @@ def boresch(
     Create a set of Boresch restraints that will restrain the 6
     external degrees of freedom of the ligand relative to the receptor.
     All of the atoms in both 'ligand' and 'receptor' must be contained in
-    'mols'.
+    'mols'. Note that restraint energies are defined as k*x**2 (so forces
+    are defined as 2*k*x) and hence the 'kr', 'ktheta' and 'kphi' values are
+    half the force constants for the distance, angle and torsion restraints.
 
     The BoreschRestraint will be a set of six restraints between
     three identified ligand atoms, and three identified receptor
@@ -47,7 +49,7 @@ def boresch(
     2. Two angle restraints, with specified force constants (ktheta)
        and equilibrium angles (theta0) parameters.
     3. Three torsion restraints, with specified force constants (kphi)
-         and equilibrium angles (phi0) parameters.
+       and equilibrium angles (phi0) parameters.
 
     This will create a single BoreschRestraint, which will be passed
     back in a BoreschRestraints object.
@@ -64,20 +66,20 @@ def boresch(
         The ligand atoms to restrain.
 
     kr : str or SireUnits::Dimension::GeneralUnit, optional
-        The force constant for the distance restraint. If None, this will
-        default to 10 kcal mol-1 A-2. Default is None.
+        Half the force constant for the distance restraint. If None, this will
+        default to 5 kcal mol-1 A-2. Default is None.
 
     ktheta : str or SireUnits::Dimension::GeneralUnit or list of str or SireUnits::Dimension::GeneralUnit, optional
-        The force constants for the angle restraints, in the order kthetaA,
-        kthetaB If None, this will default to 100 kcal mol-1 rad-2 for
+        Half the force constants for the angle restraints, in the order kthetaA,
+        kthetaB If None, this will default to 50 kcal mol-1 rad-2 for
         both angle restraints.  If a list, then this should be a list of
         length 2 containing the force constants for the two angle
         restraints. If a single value, then this will be used for both
         angle restraints. Default is None.
 
     kphi : str or SireUnits::Dimension::GeneralUnit or list of str or SireUnits::Dimension::GeneralUnit, optional
-        The force constants for the torsion restraints, in the order kthetaA,
-        kthetaB, kthetaC. If None, this will default to 100 kcal mol-1 rad-2
+        Half the force constants for the torsion restraints, in the order kthetaA,
+        kthetaB, kthetaC. If None, this will default to 50 kcal mol-1 rad-2
         for all three torsion restraints.  If a list, then this should be a
         list of length 3 containing the force constants for the three
         torsion restraints. If a single value, then this will be used for
@@ -173,8 +175,8 @@ def boresch(
 
     from .. import measure
 
-    default_distance_k = u("10 kcal mol-1 A-2")
-    default_angle_k = u("100 kcal mol-1 rad-2")
+    default_distance_k = u("5 kcal mol-1 A-2")
+    default_angle_k = u("50 kcal mol-1 rad-2")
 
     # Use the user-specified equilibrium values if they are provided.
     distance = [[ligand[0], receptor[0]]]
@@ -380,8 +382,10 @@ def distance(mols, atoms0, atoms1, r0=None, k=None, name=None, map=None):
     Create a set of distance restraints from all of the atoms in 'atoms0'
     to all of the atoms in 'atoms1' where all atoms are
     contained in the container 'mols', using the
-    passed values of the force constant 'k' and equilibrium
-    bond length r0.
+    passed values of 'k' and equilibrium bond length r0.
+    Note that 'k' corresponds to half the force constant, because
+    the restraint energy is defined as k*(r - r0)**2 (hence the force is
+    defined as 2*k*(r-r0)).
 
     These restraints will be per atom-atom distance. If a list of k and/or r0
     values are passed, then different values could be used for
@@ -489,8 +493,11 @@ def positional(mols, atoms, k=None, r0=None, position=None, name=None, map=None)
     """
     Create a set of position restraints for the atoms specified in
     'atoms' that are contained in the container 'mols', using the
-    passed values of the force constant 'k' and flat-bottom potential
-    well-width 'r0' for the restraints.
+    passed values of 'k' and flat-bottom potential
+    well-width 'r0' for the restraints. Note that 'k' values
+    correspond to half the force constants for the harmonic
+    restraints, because the harmonic restraint energy is defined as
+    k*(r - r0)**2 (hence the force is defined as 2*(r - r0)).
 
     These restraints will be per atom. If a list of k and/or r0
     values are passed, then different values could be used for

src/sire/restraints/_standard_state_correction.py

@@ -150,6 +150,10 @@ def _get_boresch_standard_state_correction(restraint, temperature):
             ).value()
             force_constants.append(force_const)
 
+    # Correct "force constants" by factor of two, because the force is defined as 2kx,
+    # rather than kx.
+    force_constants = [2 * k for k in force_constants]
+
     n_nonzero_k = len(force_constants)
 
     # Calculation

tests/convert/test_rdkit.py

@@ -11,7 +11,7 @@
     [
         "C1CCCCC1",
         "C",
-        "OCC(O)C(O)C(O)C(O)CO",
+        "OC[C@H](O)[C@H](O)[C@@H](O)[C@@H](O)CO",
         "C[C@H](N)C(=O)O",  # L-alanine
         "C[C@@H](N)C(=O)O",  # D-alanine
     ],
@@ -118,10 +118,13 @@ def test_rdkit_returns_null():
     "rdkit" not in sr.convert.supported_formats(),
     reason="rdkit support is not available",
 )
+@pytest.mark.xfail(reason="SMILES now mismatches since SDF stereochemistry is preserved")
 def test_rdkit_infer_bonds(ejm55_sdf, ejm55_gro):
     sdf = ejm55_sdf[0].molecule()
     gro = ejm55_gro["not (protein or water)"].molecule()
 
+    from rdkit import Chem
+
     assert sdf.smiles() == gro.smiles()
 
     match_sdf = sdf["smarts [NX3][CX3](=[OX1])[#6]"]

tests/restraints/test_standard_state_correction.py

@@ -1,16 +1,16 @@
 import pytest
 
 # Boresch parameters from old test_boresch_corr SOMD test so we can compare
-# to previous results.
+# to previous results. Note that the force constants are 
 BORESCH_PARAMS_DEFAULT = {
     "receptor_selection": "atomidx 1538 or atomidx 1518 or atomidx 1540",
     "ligand_selection": "atomidx 4 or atomidx 3 or atomidx 5",
-    "kr": "50.98 kcal mol-1 A-2",
-    "ktheta": ["133.48 kcal mol-1 rad-2", "76.78 kcal mol-1 rad-2"],
+    "kr": "25.49 kcal mol-1 A-2",
+    "ktheta": ["66.74 kcal mol-1 rad-2", "38.39 kcal mol-1 rad-2"],
     "kphi": [
-        "430.72 kcal mol-1 rad-2",
-        "98.46 kcal mol-1 rad-2",
-        "99.58 kcal mol-1 rad-2",
+        "215.36 kcal mol-1 rad-2",
+        "49.23 kcal mol-1 rad-2",
+        "49.79 kcal mol-1 rad-2",
     ],
     "r0": "5.92 A",
     "theta0": ["1.85 rad", "1.59 rad"],

wrapper/Convert/SireOpenMM/openmmminimise.cpp

@@ -50,6 +50,7 @@
 #include <limits.h> // CHAR_BIT
 #include <sstream>
 #include <stdint.h> // uint64_t
+#include <Python.h>
 
 inline auto is_ieee754_nan(double const x)
     -> bool
@@ -91,7 +92,6 @@ inline auto is_ieee754_nan(double const x)
 #include <algorithm>
 
 #include "SireError/errors.h"
-#include "SireBase/releasegil.h"
 
 #include "SireBase/progressbar.h"
 #include "SireUnits/units.h"
@@ -619,6 +619,8 @@ namespace SireOpenMM
                                     double starting_k, double ratchet_scale,
                                     double max_constraint_error, double timeout)
     {
+        PyThreadState *_save = PyEval_SaveThread();
+
         if (max_iterations < 0)
         {
             max_iterations = std::numeric_limits<int>::max();
@@ -629,8 +631,6 @@ namespace SireOpenMM
             timeout = std::numeric_limits<double>::max();
         }
 
-        auto gil = SireBase::release_gil();
-
         const OpenMM::System &system = context.getSystem();
 
         int num_particles = system.getNumParticles();
@@ -1105,6 +1105,8 @@ namespace SireOpenMM
                                            CODELOC);
         }
 
+        PyEval_RestoreThread(_save);
+
         return data.getLog().join("\n");
     }
 

wrapper/Convert/SireOpenMM/torchqm.cpp

@@ -391,12 +391,20 @@ double TorchQMForceImpl::computeForce(
         device = torch::kCPU;
     }
 
-    // If this is the first step, then setup information for the neighbour list.
-    if (this->step_count == 0)
+    // Move the Torch module to the correct device. Annoyingly, we have to
+    // re-load the module if the device has changed. This is because it
+    // appears that the overloaded .to() method isn't called via C++.
+    if (device != this->device)
     {
-        // Move the Torch module to the correct device.
+        this->torch_module = torch::jit::load(this->getOwner().getModulePath().toStdString());
+        this->torch_module.eval();
         this->torch_module.to(device);
+        this->device = device;
+    }
 
+    // If this is the first step, then setup information for the neighbour list.
+    if (this->step_count == 0)
+    {
         // Store the cutoff as a double in Angstom.
         this->cutoff = this->owner.getCutoff().value();
 

wrapper/Convert/SireOpenMM/torchqm.h

@@ -296,6 +296,7 @@ namespace SireOpenMM
         double neighbour_list_cutoff;
         QSet<int> neighbour_list;
         int max_num_mm=0;
+        c10::DeviceType device;
     };
 #endif
 

wrapper/Convert/SireRDKit/sire_rdkit.cpp

@@ -640,9 +640,6 @@ namespace SireRDKit
 
             a->setAtomicNum(element.nProtons());
 
-            // don't automatically add hydrogens to this atom
-            a->setNoImplicit(true);
-
             elements.append(element);
 
             try
@@ -711,7 +708,7 @@ namespace SireRDKit
                 bondtype = string_to_bondtype(bond.property(map["order"]).asA<SireMol::BondOrder>().toRDKit());
 
                 // one bond has bond info, so assume that all do
-                has_bond_info = false;
+                has_bond_info = true;
             }
             catch (...)
             {
@@ -748,6 +745,7 @@ namespace SireRDKit
             molecule.addBond(bond.atom0().index().value(),
                              bond.atom1().index().value(),
                              bondtype);
+            molecule.getBondWithIdx(i)->setStereo(stereo);
         }
 
         if (has_coords)
@@ -855,6 +853,19 @@ namespace SireRDKit
         {
         }
 
+        // try assigning stereochemistry from 3D coordinates as it is the most
+        // reliable way to do it
+        if (has_coords)
+        {
+            try
+            {
+                RDKit::MolOps::assignStereochemistryFrom3D(molecule);
+            }
+            catch (...)
+            {
+            }
+        }
+
         return ROMOL_SPTR(new RDKit::ROMol(molecule));
     }