diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 3dd88e2b1..8b0262fe6 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -21,6 +21,13 @@ organisation on `GitHub `__. * 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 `__ - October 2024 -------------------------------------------------------------------------------------------- diff --git a/doc/source/tutorial/part06/03_restraints.rst b/doc/source/tutorial/part06/03_restraints.rst index 0b35d49fd..341dca762 100644 --- a/doc/source/tutorial/part06/03_restraints.rst +++ b/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,9 +250,9 @@ 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, @@ -253,13 +260,13 @@ the system supplied. For example, >>> 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", diff --git a/doc/source/tutorial/part08/01_intro.rst b/doc/source/tutorial/part08/01_intro.rst index 166b27e17..632cbd442 100644 --- a/doc/source/tutorial/part08/01_intro.rst +++ b/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") diff --git a/doc/source/tutorial/part08/02_emle.rst b/doc/source/tutorial/part08/02_emle.rst index 84305c4c5..be1442b37 100644 --- a/doc/source/tutorial/part08/02_emle.rst +++ b/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 diff --git a/doc/source/tutorial/part08/03_adp_pmf.rst b/doc/source/tutorial/part08/03_adp_pmf.rst index 27b2b9a55..f980c6c54 100644 --- a/doc/source/tutorial/part08/03_adp_pmf.rst +++ b/doc/source/tutorial/part08/03_adp_pmf.rst @@ -55,10 +55,10 @@ calculation. We can now create a ``dynamics`` that will create an ``OpenMM`` context for us and can be used to run a simulation: ->>> d = mols.dynamics( +>>> d = qm_mols.dynamics( ... timestep="1fs", ... constraint="none", -... engine=engine, +... qm_engine=engine, ... platform="cpu", ... ) diff --git a/doc/source/tutorial/part08/04_diels_alder.rst b/doc/source/tutorial/part08/04_diels_alder.rst index a32eb781e..0e3de5413 100644 --- a/doc/source/tutorial/part08/04_diels_alder.rst +++ b/doc/source/tutorial/part08/04_diels_alder.rst @@ -101,7 +101,7 @@ keyword argument. This specifies a selection for the atoms that should be fixed during simulation. Here we will fix all atoms more than 20 Å from the reaction site: ->>> d = mols.dynamics( +>>> d = qm_mols.dynamics( ... timestep="1fs", ... constraint="none", ... perturbable_constraint="none", diff --git a/src/sire/mol/_dynamics.py b/src/sire/mol/_dynamics.py index 40c8b125a..fe6eb7d2d 100644 --- a/src/sire/mol/_dynamics.py +++ b/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 diff --git a/src/sire/restraints/_restraints.py b/src/sire/restraints/_restraints.py index d4c64774e..077cb97c0 100644 --- a/src/sire/restraints/_restraints.py +++ b/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 diff --git a/src/sire/restraints/_standard_state_correction.py b/src/sire/restraints/_standard_state_correction.py index 3420b5eb2..5366c8344 100644 --- a/src/sire/restraints/_standard_state_correction.py +++ b/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 diff --git a/tests/convert/test_rdkit.py b/tests/convert/test_rdkit.py index 780f4e3a8..304630f15 100644 --- a/tests/convert/test_rdkit.py +++ b/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]"] diff --git a/tests/restraints/test_standard_state_correction.py b/tests/restraints/test_standard_state_correction.py index 977e064dc..dab0e1bac 100644 --- a/tests/restraints/test_standard_state_correction.py +++ b/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"], diff --git a/wrapper/Convert/SireOpenMM/openmmminimise.cpp b/wrapper/Convert/SireOpenMM/openmmminimise.cpp index dd5a5f536..500401b4b 100644 --- a/wrapper/Convert/SireOpenMM/openmmminimise.cpp +++ b/wrapper/Convert/SireOpenMM/openmmminimise.cpp @@ -50,6 +50,7 @@ #include // CHAR_BIT #include #include // uint64_t +#include inline auto is_ieee754_nan(double const x) -> bool @@ -91,7 +92,6 @@ inline auto is_ieee754_nan(double const x) #include #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::max(); @@ -629,8 +631,6 @@ namespace SireOpenMM timeout = std::numeric_limits::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"); } diff --git a/wrapper/Convert/SireOpenMM/torchqm.cpp b/wrapper/Convert/SireOpenMM/torchqm.cpp index 1332c91c7..c75456c75 100644 --- a/wrapper/Convert/SireOpenMM/torchqm.cpp +++ b/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(); diff --git a/wrapper/Convert/SireOpenMM/torchqm.h b/wrapper/Convert/SireOpenMM/torchqm.h index 839ecf071..97969d75c 100644 --- a/wrapper/Convert/SireOpenMM/torchqm.h +++ b/wrapper/Convert/SireOpenMM/torchqm.h @@ -296,6 +296,7 @@ namespace SireOpenMM double neighbour_list_cutoff; QSet neighbour_list; int max_num_mm=0; + c10::DeviceType device; }; #endif diff --git a/wrapper/Convert/SireRDKit/sire_rdkit.cpp b/wrapper/Convert/SireRDKit/sire_rdkit.cpp index f5574cdc4..c5b971203 100644 --- a/wrapper/Convert/SireRDKit/sire_rdkit.cpp +++ b/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().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)); }