Restraints¶
Automated selection and imposition of receptor-ligand restraints for absolute alchemical binding free energy calculations, along with computation of the standard state correction.
-
exception
yank.restraints.
RestraintStateError
[source]¶ Error raised by an
RestraintState
.-
with_traceback
()¶ Exception.with_traceback(tb) – set self.__traceback__ to tb and return self.
-
-
exception
yank.restraints.
RestraintParameterError
[source]¶ Error raised by a
ReceptorLigandRestraint
.-
with_traceback
()¶ Exception.with_traceback(tb) – set self.__traceback__ to tb and return self.
-
-
yank.restraints.
available_restraint_classes
()[source]¶ Return all available restraint classes.
Returns: - restraint_classes : dict of {str
restraint_classes[name]
is the class corresponding toname
-
yank.restraints.
available_restraint_types
()[source]¶ List all available restraint types.
Returns: - available_restraint_types : list of str
List of names of available restraint classes
-
yank.restraints.
create_restraint
(restraint_type, **kwargs)[source]¶ Factory of receptor-ligand restraint objects.
Parameters: - restraint_type : str
Restraint type name matching a register (imported) subclass of
ReceptorLigandRestraint
.- kwargs
Parameters to pass to the restraint constructor.
-
class
yank.restraints.
RestraintState
(lambda_restraints)[source]¶ The state of a restraint.
A
ComposableState
controlling the strength of a restraint through itslambda_restraints
property.Parameters: - lambda_restraints : float
The strength of the restraint. Must be between 0 and 1.
Examples
Create a system in a thermodynamic state.
>>> from openmmtools import testsystems, states >>> system_container = testsystems.LysozymeImplicit() >>> system, positions = system_container.system, system_container.positions >>> thermodynamic_state = states.ThermodynamicState(system, 300*unit.kelvin) >>> sampler_state = states.SamplerState(positions)
Identify ligand atoms. Topography automatically identify receptor atoms too.
>>> from yank.yank import Topography >>> topography = Topography(system_container.topology, ligand_atoms=range(2603, 2621))
Apply a Harmonic restraint between receptor and protein. Let the restraint automatically determine all the parameters.
>>> restraint = Harmonic() >>> restraint.determine_missing_parameters(thermodynamic_state, sampler_state, topography) >>> restraint.restrain_state(thermodynamic_state)
Create a
RestraintState
object to control the strength of the restraint.>>> restraint_state = RestraintState(lambda_restraints=1.0)
RestraintState
implements theIComposableState
interface, so it can be used withCompoundThermodynamicState
.>>> compound_state = states.CompoundThermodynamicState(thermodynamic_state=thermodynamic_state, ... composable_states=[restraint_state]) >>> compound_state.lambda_restraints 1.0 >>> integrator = openmm.VerletIntegrator(1.0*unit.femtosecond) >>> context = compound_state.create_context(integrator) >>> context.getParameter('lambda_restraints') 1.0
You can control the parameters in the OpenMM Context by setting the state’s attributes. To To deactivate the restraint, set lambda_restraints to 0.0.
>>> compound_state.lambda_restraints = 0.0 >>> compound_state.apply_to_context(context) >>> context.getParameter('lambda_restraints') 0.0
Attributes: lambda_restraints
Float: the strength of the applied restraint (between 0 and 1 inclusive).
-
lambda_restraints
¶ Float: the strength of the applied restraint (between 0 and 1 inclusive).
-
apply_to_system
(system)[source]¶ Set the strength of the system’s restraint to this.
System is updated in-place
Parameters: - system : simtk.openmm.System
The system to modify.
Raises: - RestraintStateError
If the system does not have any
CustomForce
with alambda_restraint
global parameter.
-
check_system_consistency
(system)[source]¶ Check if the system’s restraint is in this restraint state.
It raises a
RestraintStateError
if the restraint is not consistent with the state.Parameters: - system : simtk.openmm.System
The system with the restraint to test.
Raises: - RestraintStateError
If the system is not consistent with this state.
-
class
yank.restraints.
ReceptorLigandRestraint
[source]¶ A restraint preventing a ligand from drifting too far from its receptor.
With replica exchange simulations, keeping the ligand close to the binding pocket can enhance mixing between the interacting and the decoupled state. This should be always used in implicit simulation, where there are no periodic boundary conditions.
This restraint strength is controlled by a global context parameter called
lambda_restraints
. You can easily control this variable through theRestraintState
object.Notes
Creating a subclass requires the following:
1. Implement a constructor. Optionally this can leave all or a subset of the restraint parameters undefined. In this case, you need to provide an implementation of
determine_missing_parameters()
.2. Implement
restrain_state()
that add the restrainForce
to the state’s System.- Implement
get_standard_state_correction()
to return standard state correction.
4. Optionally, implement
determine_missing_parameters()
to fill in the parameters left undefined in the constructor.-
restrain_state
(thermodynamic_state)[source]¶ Add the restraint force to the state’s System.
Parameters: - thermodynamic_state : openmmtools.states.ThermodynamicState
The thermodynamic state holding the system to modify.
-
get_standard_state_correction
(thermodynamic_state)[source]¶ Return the standard state correction.
Parameters: - thermodynamic_state : openmmtools.states.ThermodynamicState
The thermodynamic state.
-
determine_missing_parameters
(thermodynamic_state, sampler_state, topography)[source]¶ Automatically choose undefined parameters.
Optionally, a
ReceptorLigandRestraint
can support the automatic determination of all or a subset of the parameters that can be left undefined in the constructor, making implementation of this method optional.Parameters: - thermodynamic_state : openmmtools.states.ThermodynamicState
The thermodynmaic state to inspect
- sampler_state : openmmtools.states.SamplerState
The sampler state holding the positions of all atoms.
- topography : yank.Topography
The topography with labeled receptor and ligand atoms.
- Implement
-
class
yank.restraints.
RadiallySymmetricRestraint
(restrained_receptor_atoms=None, restrained_ligand_atoms=None)[source]¶ Base class for radially-symmetric restraints between ligand and protein.
The restraint is applied between the centroids of two groups of atoms that belong to the receptor and the ligand respectively. The centroids are determined by a mass-weighted average of the group particles positions. The restraint strength is controlled by a global context parameter called ‘lambda_restraints’.
With OpenCL, groups with more than 1 atom are supported only on 64bit platforms.
The class allows the restrained atoms to be temporarily undefined, but in this case,
determine_missing_parameters()
must be called before using the restraint.Parameters: - restrained_receptor_atoms : iterable of int, int, or str, optional
The indices of the receptor atoms to restrain, an MDTraj DSL expression, any other
Topography
region name, orTopography Selection
. This can temporarily be left undefined, butdetermine_missing_parameters()
must be called before using the Restraint object. The same if a DSL expression or Topography selection is provided (default is None).- restrained_ligand_atoms : iterable of int, int, or str, optional
The indices of the ligand atoms to restrain, an MDTraj DSL expression, or a
Topography
region name, orTopography Selection
. This can temporarily be left undefined, butdetermine_missing_parameters()
must be called before using the Restraint object. The same if a DSL expression or Topography selection is provided (default is None).
Notes
To create a subclass, follow these steps:
1. Implement the
_create_restraint_force()
returning the force object modeling the restraint.- Implement the property
_are_restraint_parameters_defined()
.
3. Optionally, you can overwrite the
_determine_restraint_parameters()
method to automatically determine these parameters from the atoms positions.Attributes: - restrained_receptor_atoms : list of int, str, or None
The indices of the receptor atoms to restrain, an MDTraj selection string, or a Topography selection string.
- restrained_ligand_atoms : list of int, str, None
The indices of the receptor atoms to restrain, an MDTraj selection string, or a Topography selection string.
-
restrain_state
(thermodynamic_state)[source]¶ Add the restraining Force(s) to the thermodynamic state’s system.
All the parameters must be defined at this point. An exception is raised if they are not.
Parameters: - thermodynamic_state : openmmtools.states.ThermodynamicState
The thermodynamic state holding the system to modify.
Raises: - RestraintParameterError
If the restraint has undefined parameters.
-
get_standard_state_correction
(thermodynamic_state)[source]¶ Return the standard state correction.
Parameters: - thermodynamic_state : openmmtools.states.ThermodynamicState
The thermodynamic state.
Returns: - correction : float
The unit-less standard-state correction, in kT (at the temperature of the given thermodynamic state).
-
determine_missing_parameters
(thermodynamic_state, sampler_state, topography)[source]¶ Automatically determine missing parameters.
If some parameters have been left undefined (i.e. the atoms to restrain or the restraint force parameters) this attempts to find them using the information in the states and the topography.
Parameters: - thermodynamic_state : openmmtools.states.ThermodynamicState
The thermodynamic state.
- sampler_state : openmmtools.states.SamplerState, optional
The sampler state holding the positions of all atoms.
- topography : yank.Topography, optional
The topography with labeled receptor and ligand atoms.
-
class
yank.restraints.
Harmonic
(spring_constant=None, **kwargs)[source]¶ Impose a single harmonic restraint between ligand and protein.
This can be used to prevent the ligand from drifting too far from the protein in implicit solvent calculations or to keep the ligand close to the binding pocket in the decoupled states to increase mixing.
The restraint is applied between the centroids of two groups of atoms that belong to the receptor and the ligand respectively. The centroids are determined by a mass-weighted average of the group particles positions.
The energy expression of the restraint is given by
E = lambda_restraints * (K/2)*r^2
where K is the spring constant, r is the distance between the two group centroids, and lambda_restraints is a scale factor that can be used to control the strength of the restraint. You can control
lambda_restraints
throughRestraintState
class.The class supports automatic determination of the parameters left undefined or defined by strings in the constructor through
determine_missing_parameters()
.With OpenCL, groups with more than 1 atom are supported only on 64bit platforms.
Parameters: - spring_constant : simtk.unit.Quantity, optional
The spring constant K (see energy expression above) in units compatible with joule/nanometer**2/mole (default is None).
- restrained_receptor_atoms : iterable of int, int, or str, optional
The indices of the receptor atoms to restrain, an MDTraj DSL expression, or a
Topography
region name, orTopography Select String
. This can temporarily be left undefined, butdetermine_missing_parameters()
must be called before using the Restraint object. The same if a DSL expression or Topography region is provided (default is None).- restrained_ligand_atoms : iterable of int, int, or str, optional
The indices of the ligand atoms to restrain, an MDTraj DSL expression. or a
Topography
region name, orTopography Select String
. This can temporarily be left undefined, butdetermine_missing_parameters()
must be called before using the Restraint object. The same if a DSL expression or Topography region is provided (default is None).
Examples
Create the ThermodynamicState.
>>> from openmmtools import testsystems, states >>> system_container = testsystems.LysozymeImplicit() >>> system, positions = system_container.system, system_container.positions >>> thermodynamic_state = states.ThermodynamicState(system, 300*unit.kelvin) >>> sampler_state = states.SamplerState(positions)
Identify ligand atoms. Topography automatically identify receptor atoms too.
>>> from yank.yank import Topography >>> topography = Topography(system_container.topology, ligand_atoms=range(2603, 2621))
you can create a completely defined restraint
>>> restraint = Harmonic(spring_constant=8*unit.kilojoule_per_mole/unit.nanometers**2, ... restrained_receptor_atoms=[1644, 1650, 1678], ... restrained_ligand_atoms='resname TMP')
Or automatically identify the parameters. When trying to impose a restraint with undefined parameters, RestraintParameterError is raised.
>>> restraint = Harmonic() >>> try: ... restraint.restrain_state(thermodynamic_state) ... except RestraintParameterError: ... print('There are undefined parameters. Choosing restraint parameters automatically.') ... restraint.determine_missing_parameters(thermodynamic_state, sampler_state, topography) ... restraint.restrain_state(thermodynamic_state) ... There are undefined parameters. Choosing restraint parameters automatically.
Get standard state correction.
>>> correction = restraint.get_standard_state_correction(thermodynamic_state)
Attributes: - restrained_receptor_atoms : list of int, str, or None
The indices of the receptor atoms to restrain, an MDTraj selection string, or a Topography region selection string.
- restrained_ligand_atoms : list of int, str, or None
The indices of the ligand atoms to restrain, an MDTraj selection string, or a Topography region selection string.
-
determine_missing_parameters
(thermodynamic_state, sampler_state, topography)¶ Automatically determine missing parameters.
If some parameters have been left undefined (i.e. the atoms to restrain or the restraint force parameters) this attempts to find them using the information in the states and the topography.
Parameters: - thermodynamic_state : openmmtools.states.ThermodynamicState
The thermodynamic state.
- sampler_state : openmmtools.states.SamplerState, optional
The sampler state holding the positions of all atoms.
- topography : yank.Topography, optional
The topography with labeled receptor and ligand atoms.
-
get_standard_state_correction
(thermodynamic_state)¶ Return the standard state correction.
Parameters: - thermodynamic_state : openmmtools.states.ThermodynamicState
The thermodynamic state.
Returns: - correction : float
The unit-less standard-state correction, in kT (at the temperature of the given thermodynamic state).
-
restrain_state
(thermodynamic_state)¶ Add the restraining Force(s) to the thermodynamic state’s system.
All the parameters must be defined at this point. An exception is raised if they are not.
Parameters: - thermodynamic_state : openmmtools.states.ThermodynamicState
The thermodynamic state holding the system to modify.
Raises: - RestraintParameterError
If the restraint has undefined parameters.
-
class
yank.restraints.
FlatBottom
(spring_constant=None, well_radius=None, **kwargs)[source]¶ A receptor-ligand restraint using a flat potential well with harmonic walls.
An alternative choice to receptor-ligand restraints that uses a flat potential inside most of the protein volume with harmonic restraining walls outside of this. It can be used to prevent the ligand from drifting too far from protein in implicit solvent calculations while still exploring the surface of the protein for putative binding sites.
The restraint is applied between the centroids of two groups of atoms that belong to the receptor and the ligand respectively. The centroids are determined by a mass-weighted average of the group particles positions.
More precisely, the energy expression of the restraint is given by
E = lambda_restraints * step(r-r0) * (K/2)*(r-r0)^2
where
K
is the spring constant,r
is the distance between the restrained atoms,r0
is another parameter defining the distance at which the restraint is imposed, andlambda_restraints
is a scale factor that can be used to control the strength of the restraint. You can controllambda_restraints
through the classRestraintState
.The class supports automatic determination of the parameters left undefined in the constructor through
determine_missing_parameters()
.With OpenCL, groups with more than 1 atom are supported only on 64bit platforms.
Parameters: - spring_constant : simtk.unit.Quantity, optional
The spring constant K (see energy expression above) in units compatible with joule/nanometer**2/mole (default is None).
- well_radius : simtk.unit.Quantity, optional
The distance r0 (see energy expression above) at which the harmonic restraint is imposed in units of distance (default is None).
- restrained_receptor_atoms : iterable of int, int, or str, optional
The indices of the receptor atoms to restrain, an MDTraj DSL expression, or a
Topography
region name, orTopography Select String
. This can temporarily be left undefined, butdetermine_missing_parameters()
must be called before using the Restraint object. The same if a DSL expression or Topography region is provided (default is None).- restrained_ligand_atoms : iterable of int, int, or str, optional
The indices of the ligand atoms to restrain, an MDTraj DSL expression. or a
Topography
region name, orTopography Select String
. This can temporarily be left undefined, butdetermine_missing_parameters()
must be called before using the Restraint object. The same if a DSL expression or Topography region is provided (default is None).
Examples
Create the ThermodynamicState.
>>> from openmmtools import testsystems, states >>> system_container = testsystems.LysozymeImplicit() >>> system, positions = system_container.system, system_container.positions >>> thermodynamic_state = states.ThermodynamicState(system, 298*unit.kelvin) >>> sampler_state = states.SamplerState(positions)
Identify ligand atoms. Topography automatically identify receptor atoms too.
>>> from yank.yank import Topography >>> topography = Topography(system_container.topology, ligand_atoms=range(2603, 2621))
You can create a completely defined restraint
>>> restraint = FlatBottom(spring_constant=0.6*unit.kilocalorie_per_mole/unit.angstroms**2, ... well_radius=5.2*unit.nanometers, restrained_receptor_atoms=[1644, 1650, 1678], ... restrained_ligand_atoms='resname TMP')
or automatically identify the parameters. When trying to impose a restraint with undefined parameters, RestraintParameterError is raised.
>>> restraint = FlatBottom() >>> try: ... restraint.restrain_state(thermodynamic_state) ... except RestraintParameterError: ... print('There are undefined parameters. Choosing restraint parameters automatically.') ... restraint.determine_missing_parameters(thermodynamic_state, sampler_state, topography) ... restraint.restrain_state(thermodynamic_state) ... There are undefined parameters. Choosing restraint parameters automatically.
Get standard state correction.
>>> correction = restraint.get_standard_state_correction(thermodynamic_state)
Attributes: - restrained_receptor_atoms : list of int or None
The indices of the receptor atoms to restrain, an MDTraj selection string, or a Topography region selection string.
- restrained_ligand_atoms : list of int or None
The indices of the ligand atoms to restrain, an MDTraj selection string, or a Topography region selection string.
-
determine_missing_parameters
(thermodynamic_state, sampler_state, topography)¶ Automatically determine missing parameters.
If some parameters have been left undefined (i.e. the atoms to restrain or the restraint force parameters) this attempts to find them using the information in the states and the topography.
Parameters: - thermodynamic_state : openmmtools.states.ThermodynamicState
The thermodynamic state.
- sampler_state : openmmtools.states.SamplerState, optional
The sampler state holding the positions of all atoms.
- topography : yank.Topography, optional
The topography with labeled receptor and ligand atoms.
-
get_standard_state_correction
(thermodynamic_state)¶ Return the standard state correction.
Parameters: - thermodynamic_state : openmmtools.states.ThermodynamicState
The thermodynamic state.
Returns: - correction : float
The unit-less standard-state correction, in kT (at the temperature of the given thermodynamic state).
-
restrain_state
(thermodynamic_state)¶ Add the restraining Force(s) to the thermodynamic state’s system.
All the parameters must be defined at this point. An exception is raised if they are not.
Parameters: - thermodynamic_state : openmmtools.states.ThermodynamicState
The thermodynamic state holding the system to modify.
Raises: - RestraintParameterError
If the restraint has undefined parameters.
-
class
yank.restraints.
BoreschLike
(restrained_receptor_atoms=None, restrained_ligand_atoms=None, K_r=None, r_aA0=None, K_thetaA=None, theta_A0=None, K_thetaB=None, theta_B0=None, K_phiA=None, phi_A0=None, K_phiB=None, phi_B0=None, K_phiC=None, phi_C0=None)[source]¶ Abstract class to impose Boresch-like orientational restraints on protein-ligand system. Subclasses are specific implementations of the energy functions
This restraints the ligand binding mode by constraining 1 distance, 2 angles and 3 dihedrals between 3 atoms of the receptor and 3 atoms of the ligand.
More precisely, the energy expression of the restraint is given by
E = lambda_restraints * { K_r/2 * [|r3 - l1| - r_aA0]^2 + + K_thetaA/2 * [angle(r2,r3,l1) - theta_A0]^2 + + K_thetaB/2 * [angle(r3,l1,l2) - theta_B0]^2 + + K_phiA/2 * hav(dihedral(r1,r2,r3,l1) - phi_A0) * 2 + + K_phiB/2 * hav(dihedral(r2,r3,l1,l2) - phi_B0) * 2 + + K_phiC/2 * hav(dihedral(r3,l1,l2,l3) - phi_C0) * 2 }
, where
hav
is the Haversine function(1-cos(x))/2
and the parameters are:r1
,r2
,r3
: the coordinates of the 3 receptor atoms.l1
,l2
,l3
: the coordinates of the 3 ligand atoms.K_r
: the spring constant for the restrained distance|r3 - l1|
.r_aA0
: the equilibrium distance of|r3 - l1|
.K_thetaA
,K_thetaB
: the spring constants forangle(r2,r3,l1)
andangle(r3,l1,l2)
.theta_A0
,theta_B0
: the equilibrium angles ofangle(r2,r3,l1)
andangle(r3,l1,l2)
.K_phiA
,K_phiB
,K_phiC
: the spring constants fordihedral(r1,r2,r3,l1)
,dihedral(r2,r3,l1,l2)
,dihedral(r3,l1,l2,l3)
.phi_A0
,phi_B0
,phi_C0
: the equilibrium torsion ofdihedral(r1,r2,r3,l1)
,dihedral(r2,r3,l1,l2)
,dihedral(r3,l1,l2,l3)
.lambda_restraints
: a scale factor that can be used to control the strength of the restraint.You can control
lambda_restraints
through the classRestraintState
.The class supports automatic determination of the parameters left undefined in the constructor through
determine_missing_parameters()
.This function used to be based on the Boresch orientational restraints [1] and has similar form to its energy equation
E = lambda_restraints * { K_r/2 * [|r3 - l1| - r_aA0]^2 + + K_thetaA/2 * [angle(r2,r3,l1) - theta_A0]^2 + + K_thetaB/2 * [angle(r3,l1,l2) - theta_B0]^2 + + K_phiA/2 * [dihedral(r1,r2,r3,l1) - phi_A0]^2 + + K_phiB/2 * [dihedral(r2,r3,l1,l2) - phi_B0]^2 + + K_phiC/2 * [dihedral(r3,l1,l2,l3) - phi_C0]^2 }
However, the form at the top is periodic with the dihedral angle and imposes a more steep energy penalty while still maintaining the same Taylor series expanded force and energy near phi_X0. The
*2
on thehav()
functions in the energy expression are shown as the explicit correction to thehav()
function to make the leading spring constant force consistent with the original harmonic Boresch restraint. In practice, the1/2
from thehav()
function is omitted.Warning: Symmetry corrections for symmetric ligands are not automatically applied. See Ref [1] and [2] for more information on correcting for ligand symmetry.
Warning: Only heavy atoms can be restrained. Hydrogens will automatically be excluded.
Parameters: - restrained_receptor_atoms : iterable of int, str, or None; Optional
The indices of the receptor atoms to restrain, an MDTraj DSL expression, or a
Topography
region name, orTopography Select String
. If this is a list of three ints, the receptor atoms will be restrained in order, r1, r2, r3. If there are more than three entries or the selection string resolves more than three atoms, the three restrained atoms will be chosen at random from the selection. This can temporarily be left undefined, butdetermine_missing_parameters()
must be called before using the Restraint object. The same if a DSL expression or Topography region is provided (default is None).- restrained_ligand_atoms : iterable of int, str, or None; Optional
The indices of the ligand atoms to restrain, an MDTraj DSL expression, or a
Topography
region name, orTopography Select String
. If this is a list of three ints, the receptor atoms will be restrained in order, l1, l2, l3. If there are more than three entries or the selection string resolves more than three atoms, the three restrained atoms will be chosen at random from the selection. This can temporarily be left undefined, butdetermine_missing_parameters()
must be called before using the Restraint object. The same if a DSL expression or Topography region is provided (default is None).- K_r : simtk.unit.Quantity, optional
The spring constant for the restrained distance
|r3 - l1|
(units compatible with kilocalories_per_mole/angstrom**2).- r_aA0 : simtk.unit.Quantity, optional
The equilibrium distance between r3 and l1 (units of length).
- K_thetaA, K_thetaB : simtk.unit.Quantity, optional
The spring constants for
angle(r2, r3, l1)
andangle(r3, l1, l2)
(units compatible with kilocalories_per_mole/radians**2).- theta_A0, theta_B0 : simtk.unit.Quantity, optional
The equilibrium angles of
angle(r2, r3, l1)
andangle(r3, l1, l2)
(units compatible with radians).- K_phiA, K_phiB, K_phiC : simtk.unit.Quantity, optional
The spring constants for
dihedral(r1, r2, r3, l1)
,dihedral(r2, r3, l1, l2)
anddihedral(r3,l1,l2,l3)
(units compatible with kilocalories_per_mole/radians**2).- phi_A0, phi_B0, phi_C0 : simtk.unit.Quantity, optional
The equilibrium torsion of
dihedral(r1,r2,r3,l1)
,dihedral(r2,r3,l1,l2)
anddihedral(r3,l1,l2,l3)
(units compatible with radians).
References
- [1] Boresch S, Tettinger F, Leitgeb M, Karplus M. J Phys Chem B. 107:9535, 2003.
- http://dx.doi.org/10.1021/jp0217839
- [2] Mobley DL, Chodera JD, and Dill KA. J Chem Phys 125:084902, 2006.
- https://dx.doi.org/10.1063%2F1.2221683
Examples
Create the ThermodynamicState.
>>> from openmmtools import testsystems, states >>> system_container = testsystems.LysozymeImplicit() >>> system, positions = system_container.system, system_container.positions >>> thermodynamic_state = states.ThermodynamicState(system, 298*unit.kelvin) >>> sampler_state = states.SamplerState(positions)
Identify ligand atoms. Topography automatically identify receptor atoms too.
>>> from yank.yank import Topography >>> topography = Topography(system_container.topology, ligand_atoms=range(2603, 2621))
Create a partially defined restraint
>>> restraint = Boresch(restrained_receptor_atoms=[1335, 1339, 1397], ... restrained_ligand_atoms=[2609, 2607, 2606], ... K_r=20.0*unit.kilocalories_per_mole/unit.angstrom**2, ... r_aA0=0.35*unit.nanometer)
and automatically identify the other parameters. When trying to impose a restraint with undefined parameters, RestraintParameterError is raised.
>>> try: ... restraint.restrain_state(thermodynamic_state) ... except RestraintParameterError: ... print('There are undefined parameters. Choosing restraint parameters automatically.') ... restraint.determine_missing_parameters(thermodynamic_state, sampler_state, topography) ... restraint.restrain_state(thermodynamic_state) ... There are undefined parameters. Choosing restraint parameters automatically.
Get standard state correction.
>>> correction = restraint.get_standard_state_correction(thermodynamic_state)
Attributes: - restrained_receptor_atoms : list of int
The indices of the 3 receptor atoms to restrain [r1, r2, r3].
- restrained_ligand_atoms : list of int
The indices of the 3 ligand atoms to restrain [l1, l2, l3].
-
restrain_state
(thermodynamic_state)[source]¶ Add the restraint force to the state’s
System
.Parameters: - thermodynamic_state : openmmtools.states.ThermodynamicState
The thermodynamic state holding the system to modify.
-
get_standard_state_correction
(thermodynamic_state)[source]¶ Return the standard state correction.
Parameters: - thermodynamic_state : openmmtools.states.ThermodynamicState
The thermodynamic state.
Returns: - DeltaG : float
Computed standard-state correction in dimensionless units (kT).
-
determine_missing_parameters
(thermodynamic_state, sampler_state, topography)[source]¶ Determine parameters and restrained atoms automatically.
Currently, all equilibrium values are measured from the initial structure, while spring constants are set to 20 kcal/(mol A**2) or 20 kcal/(mol rad**2) as in Ref [1]. The restrained atoms are selected so that the analytical standard state correction will be valid.
Parameters that have been already specified are left untouched.
Future iterations of this feature will introduce the ability to extract equilibrium parameters and spring constants from a short simulation.
Parameters: - thermodynamic_state : openmmtools.states.ThermodynamicState
The thermodynamic state.
- sampler_state : openmmtools.states.SamplerState, optional
The sampler state holding the positions of all atoms.
- topography : yank.Topography, optional
The topography with labeled receptor and ligand atoms.
-
class
yank.restraints.
Boresch
(restrained_receptor_atoms=None, restrained_ligand_atoms=None, K_r=None, r_aA0=None, K_thetaA=None, theta_A0=None, K_thetaB=None, theta_B0=None, K_phiA=None, phi_A0=None, K_phiB=None, phi_B0=None, K_phiC=None, phi_C0=None)[source]¶ Impose Boresch-style orientational restraints on protein-ligand system.
This restraints the ligand binding mode by constraining 1 distance, 2 angles and 3 dihedrals between 3 atoms of the receptor and 3 atoms of the ligand.
More precisely, the energy expression of the restraint is given by
E = lambda_restraints * { K_r/2 * [|r3 - l1| - r_aA0]^2 + + K_thetaA/2 * [angle(r2,r3,l1) - theta_A0]^2 + + K_thetaB/2 * [angle(r3,l1,l2) - theta_B0]^2 + + K_phiA/2 * [dihedral(r1,r2,r3,l1) - phi_A0]^2 + + K_phiB/2 * [dihedral(r2,r3,l1,l2) - phi_B0]^2 + + K_phiC/2 * [dihedral(r3,l1,l2,l3) - phi_C0]^2 }
, where the parameters are:
r1
,r2
,r3
: the coordinates of the 3 receptor atoms.l1
,l2
,l3
: the coordinates of the 3 ligand atoms.K_r
: the spring constant for the restrained distance|r3 - l1|
.r_aA0
: the equilibrium distance of|r3 - l1|
.K_thetaA
,K_thetaB
: the spring constants forangle(r2,r3,l1)
andangle(r3,l1,l2)
.theta_A0
,theta_B0
: the equilibrium angles ofangle(r2,r3,l1)
andangle(r3,l1,l2)
.K_phiA
,K_phiB
,K_phiC
: the spring constants fordihedral(r1,r2,r3,l1)
,dihedral(r2,r3,l1,l2)
,dihedral(r3,l1,l2,l3)
.phi_A0
,phi_B0
,phi_C0
: the equilibrium torsion ofdihedral(r1,r2,r3,l1)
,dihedral(r2,r3,l1,l2)
,dihedral(r3,l1,l2,l3)
.lambda_restraints
: a scale factor that can be used to control the strength of the restraint.You can control
lambda_restraints
through the classRestraintState
.The class supports automatic determination of the parameters left undefined in the constructor through
determine_missing_parameters()
.Warning: Symmetry corrections for symmetric ligands are not automatically applied. See Ref [1] and [2] for more information on correcting for ligand symmetry.
Warning: Only heavy atoms can be restrained. Hydrogens will automatically be excluded.
Parameters: - restrained_receptor_atoms : iterable of int, str, or None; Optional
The indices of the receptor atoms to restrain, an MDTraj DSL expression, or a
Topography
region name, orTopography Select String
. If this is a list of three ints, the receptor atoms will be restrained in order, r1, r2, r3. If there are more than three entries or the selection string resolves more than three atoms, the three restrained atoms will be chosen at random from the selection. This can temporarily be left undefined, butdetermine_missing_parameters()
must be called before using the Restraint object. The same if a DSL expression or Topography region is provided (default is None).- restrained_ligand_atoms : iterable of int, str, or None; Optional
The indices of the ligand atoms to restrain, an MDTraj DSL expression, or a
Topography
region name, orTopography Select String
. If this is a list of three ints, the receptor atoms will be restrained in order, l1, l2, l3. If there are more than three entries or the selection string resolves more than three atoms, the three restrained atoms will be chosen at random from the selection. This can temporarily be left undefined, butdetermine_missing_parameters()
must be called before using the Restraint object. The same if a DSL expression or Topography region is provided (default is None).- K_r : simtk.unit.Quantity, optional
The spring constant for the restrained distance
|r3 - l1|
(units compatible with kilocalories_per_mole/angstrom**2).- r_aA0 : simtk.unit.Quantity, optional
The equilibrium distance between r3 and l1 (units of length).
- K_thetaA, K_thetaB : simtk.unit.Quantity, optional
The spring constants for
angle(r2, r3, l1)
andangle(r3, l1, l2)
(units compatible with kilocalories_per_mole/radians**2).- theta_A0, theta_B0 : simtk.unit.Quantity, optional
The equilibrium angles of
angle(r2, r3, l1)
andangle(r3, l1, l2)
(units compatible with radians).- K_phiA, K_phiB, K_phiC : simtk.unit.Quantity, optional
The spring constants for
dihedral(r1, r2, r3, l1)
,dihedral(r2, r3, l1, l2)
anddihedral(r3,l1,l2,l3)
(units compatible with kilocalories_per_mole/radians**2).- phi_A0, phi_B0, phi_C0 : simtk.unit.Quantity, optional
The equilibrium torsion of
dihedral(r1,r2,r3,l1)
,dihedral(r2,r3,l1,l2)
anddihedral(r3,l1,l2,l3)
(units compatible with radians).
References
- [1] Boresch S, Tettinger F, Leitgeb M, Karplus M. J Phys Chem B. 107:9535, 2003.
- http://dx.doi.org/10.1021/jp0217839
- [2] Mobley DL, Chodera JD, and Dill KA. J Chem Phys 125:084902, 2006.
- https://dx.doi.org/10.1063%2F1.2221683
Examples
Create the ThermodynamicState.
>>> from openmmtools import testsystems, states >>> system_container = testsystems.LysozymeImplicit() >>> system, positions = system_container.system, system_container.positions >>> thermodynamic_state = states.ThermodynamicState(system, 298*unit.kelvin) >>> sampler_state = states.SamplerState(positions)
Identify ligand atoms. Topography automatically identify receptor atoms too.
>>> from yank.yank import Topography >>> topography = Topography(system_container.topology, ligand_atoms=range(2603, 2621))
Create a partially defined restraint
>>> restraint = Boresch(restrained_receptor_atoms=[1335, 1339, 1397], ... restrained_ligand_atoms=[2609, 2607, 2606], ... K_r=20.0*unit.kilocalories_per_mole/unit.angstrom**2, ... r_aA0=0.35*unit.nanometer)
and automatically identify the other parameters. When trying to impose a restraint with undefined parameters, RestraintParameterError is raised.
>>> try: ... restraint.restrain_state(thermodynamic_state) ... except RestraintParameterError: ... print('There are undefined parameters. Choosing restraint parameters automatically.') ... restraint.determine_missing_parameters(thermodynamic_state, sampler_state, topography) ... restraint.restrain_state(thermodynamic_state) ... There are undefined parameters. Choosing restraint parameters automatically.
Get standard state correction.
>>> correction = restraint.get_standard_state_correction(thermodynamic_state)
Attributes: - restrained_receptor_atoms : list of int
The indices of the 3 receptor atoms to restrain [r1, r2, r3].
- restrained_ligand_atoms : list of int
The indices of the 3 ligand atoms to restrain [l1, l2, l3].
-
determine_missing_parameters
(thermodynamic_state, sampler_state, topography)¶ Determine parameters and restrained atoms automatically.
Currently, all equilibrium values are measured from the initial structure, while spring constants are set to 20 kcal/(mol A**2) or 20 kcal/(mol rad**2) as in Ref [1]. The restrained atoms are selected so that the analytical standard state correction will be valid.
Parameters that have been already specified are left untouched.
Future iterations of this feature will introduce the ability to extract equilibrium parameters and spring constants from a short simulation.
Parameters: - thermodynamic_state : openmmtools.states.ThermodynamicState
The thermodynamic state.
- sampler_state : openmmtools.states.SamplerState, optional
The sampler state holding the positions of all atoms.
- topography : yank.Topography, optional
The topography with labeled receptor and ligand atoms.
-
get_standard_state_correction
(thermodynamic_state)¶ Return the standard state correction.
Parameters: - thermodynamic_state : openmmtools.states.ThermodynamicState
The thermodynamic state.
Returns: - DeltaG : float
Computed standard-state correction in dimensionless units (kT).
-
restrain_state
(thermodynamic_state)¶ Add the restraint force to the state’s
System
.Parameters: - thermodynamic_state : openmmtools.states.ThermodynamicState
The thermodynamic state holding the system to modify.
-
class
yank.restraints.
PeriodicTorsionBoresch
(restrained_receptor_atoms=None, restrained_ligand_atoms=None, K_r=None, r_aA0=None, K_thetaA=None, theta_A0=None, K_thetaB=None, theta_B0=None, K_phiA=None, phi_A0=None, K_phiB=None, phi_B0=None, K_phiC=None, phi_C0=None)[source]¶ Impose Boresch-style orientational restraints on protein-ligand system where torsions are restrained by a periodic instead of harmonic force
This restraints the ligand binding mode by constraining 1 distance, 2 angles and 3 dihedrals between 3 atoms of the receptor and 3 atoms of the ligand.
More precisely, the energy expression of the restraint is given by
E = lambda_restraints * { K_r/2 * [|r3 - l1| - r_aA0]^2 + + K_thetaA/2 * [angle(r2,r3,l1) - theta_A0]^2 + + K_thetaB/2 * [angle(r3,l1,l2) - theta_B0]^2 + + K_phiA/2 * hav(dihedral(r1,r2,r3,l1) - phi_A0) * 2 + + K_phiB/2 * hav(dihedral(r2,r3,l1,l2) - phi_B0) * 2 + + K_phiC/2 * hav(dihedral(r3,l1,l2,l3) - phi_C0) * 2 }
, where
hav
is the Haversine function(1-cos(x))/2
and the parameters are:r1
,r2
,r3
: the coordinates of the 3 receptor atoms.l1
,l2
,l3
: the coordinates of the 3 ligand atoms.K_r
: the spring constant for the restrained distance|r3 - l1|
.r_aA0
: the equilibrium distance of|r3 - l1|
.K_thetaA
,K_thetaB
: the spring constants forangle(r2,r3,l1)
andangle(r3,l1,l2)
.theta_A0
,theta_B0
: the equilibrium angles ofangle(r2,r3,l1)
andangle(r3,l1,l2)
.K_phiA
,K_phiB
,K_phiC
: the spring constants fordihedral(r1,r2,r3,l1)
,dihedral(r2,r3,l1,l2)
,dihedral(r3,l1,l2,l3)
.phi_A0
,phi_B0
,phi_C0
: the equilibrium torsion ofdihedral(r1,r2,r3,l1)
,dihedral(r2,r3,l1,l2)
,dihedral(r3,l1,l2,l3)
.lambda_restraints
: a scale factor that can be used to control the strength of the restraint.You can control
lambda_restraints
through the classRestraintState
.The class supports automatic determination of the parameters left undefined in the constructor through
determine_missing_parameters()
.This function used to be based on the Boresch orientational restraints [1] and has similar form to its energy equation
E = lambda_restraints * { K_r/2 * [|r3 - l1| - r_aA0]^2 + + K_thetaA/2 * [angle(r2,r3,l1) - theta_A0]^2 + + K_thetaB/2 * [angle(r3,l1,l2) - theta_B0]^2 + + K_phiA/2 * [dihedral(r1,r2,r3,l1) - phi_A0]^2 + + K_phiB/2 * [dihedral(r2,r3,l1,l2) - phi_B0]^2 + + K_phiC/2 * [dihedral(r3,l1,l2,l3) - phi_C0]^2 }
However, the form at the top is periodic with the dihedral angle and imposes a more steep energy penalty while still maintaining the same Taylor series expanded force and energy near phi_X0. The
*2
on thehav()
functions in the energy expression are shown as the explicit correction to thehav()
function to make the leading spring constant force consistent with the original harmonic Boresch restraint. In practice, the1/2
from thehav()
function is omitted.Warning: Symmetry corrections for symmetric ligands are not automatically applied. See Ref [1] and [2] for more information on correcting for ligand symmetry.
Warning: Only heavy atoms can be restrained. Hydrogens will automatically be excluded.
Parameters: - restrained_receptor_atoms : iterable of int, str, or None; Optional
The indices of the receptor atoms to restrain, an MDTraj DSL expression, or a
Topography
region name, orTopography Select String
. If this is a list of three ints, the receptor atoms will be restrained in order, r1, r2, r3. If there are more than three entries or the selection string resolves more than three atoms, the three restrained atoms will be chosen at random from the selection. This can temporarily be left undefined, butdetermine_missing_parameters()
must be called before using the Restraint object. The same if a DSL expression or Topography region is provided (default is None).- restrained_ligand_atoms : iterable of int, str, or None; Optional
The indices of the ligand atoms to restrain, an MDTraj DSL expression, or a
Topography
region name, orTopography Select String
. If this is a list of three ints, the receptor atoms will be restrained in order, l1, l2, l3. If there are more than three entries or the selection string resolves more than three atoms, the three restrained atoms will be chosen at random from the selection. This can temporarily be left undefined, butdetermine_missing_parameters()
must be called before using the Restraint object. The same if a DSL expression or Topography region is provided (default is None).- K_r : simtk.unit.Quantity, optional
The spring constant for the restrained distance
|r3 - l1|
(units compatible with kilocalories_per_mole/angstrom**2).- r_aA0 : simtk.unit.Quantity, optional
The equilibrium distance between r3 and l1 (units of length).
- K_thetaA, K_thetaB : simtk.unit.Quantity, optional
The spring constants for
angle(r2, r3, l1)
andangle(r3, l1, l2)
(units compatible with kilocalories_per_mole/radians**2).- theta_A0, theta_B0 : simtk.unit.Quantity, optional
The equilibrium angles of
angle(r2, r3, l1)
andangle(r3, l1, l2)
(units compatible with radians).- K_phiA, K_phiB, K_phiC : simtk.unit.Quantity, optional
The spring constants for
dihedral(r1, r2, r3, l1)
,dihedral(r2, r3, l1, l2)
anddihedral(r3,l1,l2,l3)
(units compatible with kilocalories_per_mole/radians**2).- phi_A0, phi_B0, phi_C0 : simtk.unit.Quantity, optional
The equilibrium torsion of
dihedral(r1,r2,r3,l1)
,dihedral(r2,r3,l1,l2)
anddihedral(r3,l1,l2,l3)
(units compatible with radians).
References
- [1] Boresch S, Tettinger F, Leitgeb M, Karplus M. J Phys Chem B. 107:9535, 2003.
- http://dx.doi.org/10.1021/jp0217839
- [2] Mobley DL, Chodera JD, and Dill KA. J Chem Phys 125:084902, 2006.
- https://dx.doi.org/10.1063%2F1.2221683
Attributes: - restrained_receptor_atoms : list of int
The indices of the 3 receptor atoms to restrain [r1, r2, r3].
- restrained_ligand_atoms : list of int
The indices of the 3 ligand atoms to restrain [l1, l2, l3].
-
determine_missing_parameters
(thermodynamic_state, sampler_state, topography)¶ Determine parameters and restrained atoms automatically.
Currently, all equilibrium values are measured from the initial structure, while spring constants are set to 20 kcal/(mol A**2) or 20 kcal/(mol rad**2) as in Ref [1]. The restrained atoms are selected so that the analytical standard state correction will be valid.
Parameters that have been already specified are left untouched.
Future iterations of this feature will introduce the ability to extract equilibrium parameters and spring constants from a short simulation.
Parameters: - thermodynamic_state : openmmtools.states.ThermodynamicState
The thermodynamic state.
- sampler_state : openmmtools.states.SamplerState, optional
The sampler state holding the positions of all atoms.
- topography : yank.Topography, optional
The topography with labeled receptor and ligand atoms.
-
get_standard_state_correction
(thermodynamic_state)¶ Return the standard state correction.
Parameters: - thermodynamic_state : openmmtools.states.ThermodynamicState
The thermodynamic state.
Returns: - DeltaG : float
Computed standard-state correction in dimensionless units (kT).
-
restrain_state
(thermodynamic_state)¶ Add the restraint force to the state’s
System
.Parameters: - thermodynamic_state : openmmtools.states.ThermodynamicState
The thermodynamic state holding the system to modify.
-
class
yank.restraints.
RMSD
(restrained_receptor_atoms=None, restrained_ligand_atoms=None, K_RMSD=None, RMSD0=None, reference_sampler_state=None)[source]¶ Impose RMSD restraint on protein-ligand system.
This restrains both protein and ligand using a flat-bottom RMSD restraint of the form:
E = lambda_restraints * step(RMSD-RMSD0) * (K/2)*(RMSD-RMSD0)^2
Parameters: - restrained_receptor_atoms : iterable of int, str, or None; Optional
The indices of the receptor atoms to restrain, an MDTraj DSL expression, or a
Topography
region name, orTopography Select String
. Any number of receptor atoms can be selected. This can temporarily be left undefined, butdetermine_missing_parameters()
must be called before using the Restraint object. The same if a DSL expression or Topography region is provided. If no selection is given, all receptor atoms will be restrained. If an empty list is provided, no receptor atoms will be restrained. (default is None).- restrained_ligand_atoms : iterable of int, str, or None; Optional
The indices of the ligand atoms to restrain, an MDTraj DSL expression, or a
Topography
region name, orTopography Select String
. Any number of ligand atoms can be selected This can temporarily be left undefined, butdetermine_missing_parameters()
must be called before using the Restraint object. The same if a DSL expression or Topography region is provided If no selection is given, all ligand atoms will be restrained. If an empty list is provided, no receptor atoms will be restrained. (default is None).- K_RMSD : simtk.unit.Quantity, optional, default=0.6*kilocalories_per_mole/angstrom**2
The spring constant (units compatible with kilocalories_per_mole/angstrom**2).
- RMSD0 : simtk.unit.Quantity, optional, default=2.0*angstrom
The RMSD at which the restraint becomes nonzero.
- reference_sampler_state : openmmtools.states.SamplerState or None, Optional
Reference sampler state with coordinates to use as the structure to align the RMSD to. The sampler state must have the same number of particles as the thermodynamic state you will apply this force to. This can temporarily be left undefined, but
determine_missing_parameters()
must be called before using the Restraint object.
Examples
Create the ThermodynamicState.
>>> from openmmtools import testsystems, states >>> system_container = testsystems.LysozymeImplicit() >>> system, positions = system_container.system, system_container.positions >>> thermodynamic_state = states.ThermodynamicState(system, 298*unit.kelvin) >>> sampler_state = states.SamplerState(positions)
Identify ligand atoms. Topography automatically identify receptor atoms too.
>>> from yank.yank import Topography >>> topography = Topography(system_container.topology, ligand_atoms=range(2603, 2621))
Create a restraint
>>> restraint = RMSD(restrained_receptor_atoms=[1335, 1339, 1397], ... restrained_ligand_atoms=[2609, 2607, 2606], ... K_RMSD=1.0*unit.kilocalories_per_mole/unit.angstrom**2, ... RMSD0=1*unit.angstroms)
Find missing parameters
>>> restraint.determine_missing_parameters(thermodynamic_state, sampler_state, topography)
Get standard state correction.
>>> correction = restraint.get_standard_state_correction(thermodynamic_state)
Attributes: - restrained_receptor_atoms : list of int
The indices of the restrained receptor atoms
- restrained_ligand_atoms : list of int
The indices of the restrained_ligand_atoms
-
dev_validation
(wrapped_function)¶ Decorator function which will only execute the wrapped_function if
validate()
is true, else will do nothing.
-
restrain_state
(thermodynamic_state)[source]¶ Add the restraint force to the state’s
System
.Parameters: - thermodynamic_state : openmmtools.states.ThermodynamicState
The thermodynamic state holding the system to modify.
-
get_standard_state_correction
(thermodynamic_state)[source]¶ Return the standard state correction.
Parameters: - thermodynamic_state : openmmtools.states.ThermodynamicState
The thermodynamic state.
Returns: - DeltaG : float
Computed standard-state correction in dimensionless units (kT).
-
determine_missing_parameters
(thermodynamic_state, sampler_state, topography)[source]¶ Set reference positions for RMSD restraint.
Future iterations of this feature will introduce the ability to extract equilibrium parameters and spring constants from a short simulation.
Parameters: - thermodynamic_state : openmmtools.states.ThermodynamicState
The thermodynamic state.
- sampler_state : openmmtools.states.SamplerState, optional
The sampler state holding the positions of all atoms to be used as reference
- topography : yank.Topography, optional
The topography with labeled receptor and ligand atoms.