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

Error raised by an RestraintState.

with_traceback()

Exception.with_traceback(tb) – set self.__traceback__ to tb and return self.

exception yank.restraints.RestraintParameterError

Error raised by a ReceptorLigandRestraint.

with_traceback()

Exception.with_traceback(tb) – set self.__traceback__ to tb and return self.

yank.restraints.available_restraint_classes()

Return all available restraint classes.

Returns:

restraint_classes : dict of {str

restraint_classes[name] is the class corresponding to name

yank.restraints.available_restraint_types()

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)

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)

The state of a restraint.

A ComposableState controlling the strength of a restraint through its lambda_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 the IComposableState interface, so it can be used with CompoundThermodynamicState.

>>> 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)

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 a lambda_restraint global parameter.

check_system_consistency(system)

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.

apply_to_context(context)

Put the restraint in the Context into this state.

Parameters:

context : simtk.openmm.Context

The context to set.

Raises:

RestraintStateError

If the context does not have the required lambda global variables.

class yank.restraints.ReceptorLigandRestraint

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 the RestraintState 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 restrain Force to the state’s System.

3. 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)

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)

Return the standard state correction.

Parameters:

thermodynamic_state : openmmtools.states.ThermodynamicState

The thermodynamic state.

determine_missing_parameters(thermodynamic_state, sampler_state, topography)

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.

class yank.restraints.RadiallySymmetricRestraint(restrained_receptor_atoms=None, restrained_ligand_atoms=None)

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, or Topography Selection. This can temporarily be left undefined, but determine_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, or Topography Selection. This can temporarily be left undefined, but determine_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.

2. 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)

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)

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)

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)

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 through RestraintState 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, or Topography Select String. This can temporarily be left undefined, but determine_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, or Topography Select String. This can temporarily be left undefined, but determine_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)

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, and lambda_restraints is a scale factor that can be used to control the strength of the restraint. You can control lambda_restraints through the class RestraintState.

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, or Topography Select String. This can temporarily be left undefined, but determine_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, or Topography Select String. This can temporarily be left undefined, but determine_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)

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 for angle(r2,r3,l1) and angle(r3,l1,l2).

theta_A0, theta_B0: the equilibrium angles of angle(r2,r3,l1) and angle(r3,l1,l2).

K_phiA, K_phiB, K_phiC: the spring constants for dihedral(r1,r2,r3,l1), dihedral(r2,r3,l1,l2), dihedral(r3,l1,l2,l3).

phi_A0, phi_B0, phi_C0: the equilibrium torsion of dihedral(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 class RestraintState.

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 the hav() functions in the energy expression are shown as the explicit correction to the hav() function to make the leading spring constant force consistent with the original harmonic Boresch restraint. In practice, the 1/2 from the hav() 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, or Topography 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, but determine_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, or Topography 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, but determine_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) and angle(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) and angle(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) and dihedral(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) and dihedral(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)

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)

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)

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)

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 for angle(r2,r3,l1) and angle(r3,l1,l2).

theta_A0, theta_B0: the equilibrium angles of angle(r2,r3,l1) and angle(r3,l1,l2).

K_phiA, K_phiB, K_phiC: the spring constants for dihedral(r1,r2,r3,l1), dihedral(r2,r3,l1,l2), dihedral(r3,l1,l2,l3).

phi_A0, phi_B0, phi_C0: the equilibrium torsion of dihedral(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 class RestraintState.

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, or Topography 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, but determine_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, or Topography 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, but determine_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) and angle(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) and angle(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) and dihedral(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) and dihedral(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)

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 for angle(r2,r3,l1) and angle(r3,l1,l2).

theta_A0, theta_B0: the equilibrium angles of angle(r2,r3,l1) and angle(r3,l1,l2).

K_phiA, K_phiB, K_phiC: the spring constants for dihedral(r1,r2,r3,l1), dihedral(r2,r3,l1,l2), dihedral(r3,l1,l2,l3).

phi_A0, phi_B0, phi_C0: the equilibrium torsion of dihedral(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 class RestraintState.

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 the hav() functions in the energy expression are shown as the explicit correction to the hav() function to make the leading spring constant force consistent with the original harmonic Boresch restraint. In practice, the 1/2 from the hav() 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, or Topography 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, but determine_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, or Topography 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, but determine_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) and angle(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) and angle(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) and dihedral(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) and dihedral(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)

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, or Topography Select String. Any number of receptor atoms can be selected. This can temporarily be left undefined, but determine_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, or Topography Select String. Any number of ligand atoms can be selected This can temporarily be left undefined, but determine_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)

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)

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)

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.