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 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"],