From be8aa3a5c0bcb7b3fb2598ccc2048e305937942e Mon Sep 17 00:00:00 2001 From: finlayclark Date: Fri, 1 Nov 2024 13:34:23 +0000 Subject: [PATCH 01/13] Make Boresch restraints correction consistent with F=2kx Previously, the correction was implemented assuming F=kx. This lead to erroneously favourable absolute free energies of binding by 1.23 kcal mol-1. --- src/sire/restraints/_standard_state_correction.py | 4 ++++ tests/restraints/test_standard_state_correction.py | 12 ++++++------ 2 files changed, 10 insertions(+), 6 deletions(-) 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/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"], From 63e76ea91ac369a9388466c1af87d00b7b6575dd Mon Sep 17 00:00:00 2001 From: finlayclark Date: Fri, 1 Nov 2024 13:55:21 +0000 Subject: [PATCH 02/13] Update docs to make clear that F = 2kx for restraints Update the restraint function docstrings, restraints tutorial, and changelog. --- doc/source/changelog.rst | 3 ++ doc/source/tutorial/part06/03_restraints.rst | 31 ++++++++++------- src/sire/restraints/_restraints.py | 35 ++++++++++++-------- 3 files changed, 43 insertions(+), 26 deletions(-) diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 3dd88e2b1..424e8dcf7 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -16,6 +16,9 @@ organisation on `GitHub `__. --------------------------------------------------------------------------------------------- * Please add an item to this CHANGELOG for any new features or bug fixes when creating a PR. +* 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. * Fixed instantiaton of ``QByteArray`` in ``Sire::Mol::Frame::toByteArray`` and count bytes with ``QByteArray::size``. * Increase timeout before terminating ``QThread`` objects during ``PageCache`` cleanup. * Expose missing ``timeout`` kwarg in :meth:`dynamics.minimise()` method. 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/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 From 1bc4e46d6e552f77f2641269a6e7e8402535484a Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Fri, 29 Nov 2024 11:19:40 +0000 Subject: [PATCH 03/13] Workaround for device change during minimisation. [ref #262] --- wrapper/Convert/SireOpenMM/torchqm.cpp | 14 +++++++++++--- wrapper/Convert/SireOpenMM/torchqm.h | 1 + 2 files changed, 12 insertions(+), 3 deletions(-) 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 From 3e866dee1d3967c7e2b7bceaa843235c8e9ccebc Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Fri, 29 Nov 2024 11:22:45 +0000 Subject: [PATCH 04/13] Add minimisation step to QM/MM tutorials. [ref #262] --- doc/source/tutorial/part08/01_intro.rst | 8 ++++++-- doc/source/tutorial/part08/02_emle.rst | 4 ++++ 2 files changed, 10 insertions(+), 2 deletions(-) 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 From c064358843a8d38a1be788c3ff1e6a70d03d30a5 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Fri, 29 Nov 2024 11:43:24 +0000 Subject: [PATCH 05/13] Fix QM-to-MM non-equilibrium work. [ref #262] --- src/sire/mol/_dynamics.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) 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 From 217906a09fd85319f5273af285bfce0a67dc0969 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 2 Dec 2024 09:07:57 +0000 Subject: [PATCH 06/13] Update CHANGELOG. [ci skip] --- doc/source/changelog.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 3dd88e2b1..a23999dbb 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -21,6 +21,8 @@ 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. `2024.3.0 `__ - October 2024 -------------------------------------------------------------------------------------------- From 9c3aa12407b115355d6f231febb9935f6a54ae4c Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 2 Dec 2024 09:59:47 +0000 Subject: [PATCH 07/13] Fix setting of RDKit bond orders. [ref #261] --- wrapper/Convert/SireRDKit/sire_rdkit.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/wrapper/Convert/SireRDKit/sire_rdkit.cpp b/wrapper/Convert/SireRDKit/sire_rdkit.cpp index f5574cdc4..ab52f7452 100644 --- a/wrapper/Convert/SireRDKit/sire_rdkit.cpp +++ b/wrapper/Convert/SireRDKit/sire_rdkit.cpp @@ -711,7 +711,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 (...) { From 41d1f07817f978310b191002ae7021c8ef461bd0 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 2 Dec 2024 10:30:47 +0000 Subject: [PATCH 08/13] Set bond stereo information. [ref #261] --- wrapper/Convert/SireRDKit/sire_rdkit.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/wrapper/Convert/SireRDKit/sire_rdkit.cpp b/wrapper/Convert/SireRDKit/sire_rdkit.cpp index ab52f7452..412b45bbe 100644 --- a/wrapper/Convert/SireRDKit/sire_rdkit.cpp +++ b/wrapper/Convert/SireRDKit/sire_rdkit.cpp @@ -748,6 +748,7 @@ namespace SireRDKit molecule.addBond(bond.atom0().index().value(), bond.atom1().index().value(), bondtype); + molecule.getBondWithIdx(i)->setStereo(stereo); } if (has_coords) From b4a1c3b0a462e34e9d557fd7d1f5bd0379cf5802 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 2 Dec 2024 15:01:22 +0000 Subject: [PATCH 09/13] Allow implicit hydrogens. [ref #261] --- wrapper/Convert/SireRDKit/sire_rdkit.cpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/wrapper/Convert/SireRDKit/sire_rdkit.cpp b/wrapper/Convert/SireRDKit/sire_rdkit.cpp index 412b45bbe..8270583fc 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 From 3720ab91b63fe30949e38a52e3ef5729f5ee37eb Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 2 Dec 2024 15:01:56 +0000 Subject: [PATCH 10/13] Assign sterochemistry from 3D coordinates. [ref #261] [ci skip] --- wrapper/Convert/SireRDKit/sire_rdkit.cpp | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/wrapper/Convert/SireRDKit/sire_rdkit.cpp b/wrapper/Convert/SireRDKit/sire_rdkit.cpp index 8270583fc..c5b971203 100644 --- a/wrapper/Convert/SireRDKit/sire_rdkit.cpp +++ b/wrapper/Convert/SireRDKit/sire_rdkit.cpp @@ -853,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)); } From 88646ff71f05de814886b3fec62436b636ab3e55 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 2 Dec 2024 16:27:24 +0000 Subject: [PATCH 11/13] XFail test since SDF stereochemistry is now preserved. [ref #261] --- tests/convert/test_rdkit.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) 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]"] From 918690d76adc19554ae124f1529a8fb043f12e1a Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 2 Dec 2024 16:29:08 +0000 Subject: [PATCH 12/13] Update CHANGELOG. --- doc/source/changelog.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index a23999dbb..2e256e704 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -23,6 +23,7 @@ organisation on `GitHub `__. * 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. `2024.3.0 `__ - October 2024 -------------------------------------------------------------------------------------------- From 1e253b1d79d604de141dc802bb40b00b2484ab7e Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 3 Dec 2024 09:18:31 +0000 Subject: [PATCH 13/13] Fix thread safety issue in OpenMM minimiser. [closes #259] --- doc/source/changelog.rst | 7 ++++--- wrapper/Convert/SireOpenMM/openmmminimise.cpp | 8 +++++--- 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 9b07baad9..8b0262fe6 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -16,9 +16,6 @@ organisation on `GitHub `__. --------------------------------------------------------------------------------------------- * Please add an item to this CHANGELOG for any new features or bug fixes when creating a PR. -* 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. * Fixed instantiaton of ``QByteArray`` in ``Sire::Mol::Frame::toByteArray`` and count bytes with ``QByteArray::size``. * Increase timeout before terminating ``QThread`` objects during ``PageCache`` cleanup. * Expose missing ``timeout`` kwarg in :meth:`dynamics.minimise()` method. @@ -27,6 +24,10 @@ organisation on `GitHub `__. * 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/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"); }