Restraints¶
It can be useful to add restraints to the system to hold things in place. For example, you might want to hold the protein fixed while allowing the ligand to move around. Or you may want to hold a ligand in place while it is being decoupled from the simulation. Or you may want to use restraints to control the distances between atoms or pull a molecule into a particular conformation.
You can specify the restraints you want to add to a system via the functions
in the sire.restraints module. These functions return
Restraints objects that contain all of the information that
needs to be passed to OpenMM to add the restraints to the system.
Positional Restraints¶
The simplest type of restraint is a positional restraint, which holds a set of
atoms at a fixed position. You can create a positional restraint using the
positional() function. This function takes the
set of molecules or system you want to simulate, together with a search
string, list of atom indexes, or molecule views holding the atoms that
you want to restrain. For example,
>>> import sire as sr
>>> mols = sr.load(sr.expand(sr.tutorial_url, "ala.top", "ala.crd"))
>>> restraints = sr.restraints.positional(mols, atoms="resname ALA and element C")
>>> print(restraints)
PositionalRestraints( size=3
0: PositionalRestraint( 8 => ( 16.5371, 5.02707, 15.812 ), k=150 kcal mol-1 Å-2 : r0=0 Å )
1: PositionalRestraint( 10 => ( 16.0464, 6.38937, 15.2588 ), k=150 kcal mol-1 Å-2 : r0=0 Å )
2: PositionalRestraint( 14 => ( 15.3698, 4.19397, 16.434 ), k=150 kcal mol-1 Å-2 : r0=0 Å )
)
will create positional restraints for the three carbon atoms in the alanine
residue. These carbon atoms are at indexes 8, 10 and 14 in mols.atoms(),
e.g.
>>> print(mols.atoms([8, 10, 14]))
SireMol::SelectorM<SireMol::Atom>( size=3
0: MolNum(6) Atom( CA:9 [ 16.54, 5.03, 15.81] )
1: MolNum(6) Atom( CB:11 [ 16.05, 6.39, 15.26] )
2: MolNum(6) Atom( C:15 [ 15.37, 4.19, 16.43] )
)
You could have specified these indicies directly, e.g.
>>> restraints = sr.restraints.positional(mols, atoms=[8, 10, 14])
>>> print(restraints)
PositionalRestraints( size=3
0: PositionalRestraint( 8 => ( 16.5371, 5.02707, 15.812 ), k=150 kcal mol-1 Å-2 : r0=0 Å )
1: PositionalRestraint( 10 => ( 16.0464, 6.38937, 15.2588 ), k=150 kcal mol-1 Å-2 : r0=0 Å )
2: PositionalRestraint( 14 => ( 15.3698, 4.19397, 16.434 ), k=150 kcal mol-1 Å-2 : r0=0 Å )
)
or just passed in the atoms directly, e.g.
>>> atoms = mols.atoms([8, 10, 14])
>>> restraints = sr.restraints.positional(mols, atoms=atoms)
>>> print(restraints)
PositionalRestraints( size=3
0: PositionalRestraint( 8 => ( 16.5371, 5.02707, 15.812 ), k=150 kcal mol-1 Å-2 : r0=0 Å )
1: PositionalRestraint( 10 => ( 16.0464, 6.38937, 15.2588 ), k=150 kcal mol-1 Å-2 : r0=0 Å )
2: PositionalRestraint( 14 => ( 15.3698, 4.19397, 16.434 ), k=150 kcal mol-1 Å-2 : r0=0 Å )
)
The default force constant is 150 kcal mol-1 Å-2. By default, the
atoms are restrained to their current positions, and are held exactly in
those positions (hence why the r0=0 Å in the output above).
You can set the force constant and width of the half-harmonic restraint
used to hold the atoms in position using the k and r0 arguments, e.g.
>>> restraints = sr.restraints.positional(mols, atoms="resname ALA and element C",
... k="100 kcal mol-1 A-2", r0="1.5 A")
>>> print(restraints)
PositionalRestraints( size=3
0: PositionalRestraint( 8 => ( 16.5371, 5.02707, 15.812 ), k=100 kcal mol-1 Å-2 : r0=1.5 Å )
1: PositionalRestraint( 10 => ( 16.0464, 6.38937, 15.2588 ), k=100 kcal mol-1 Å-2 : r0=1.5 Å )
2: PositionalRestraint( 14 => ( 15.3698, 4.19397, 16.434 ), k=100 kcal mol-1 Å-2 : r0=1.5 Å )
)
By default, the atoms are restrained to their current positions. You can
specify a different position using the position argument, e.g.
>>> restraints = sr.restraints.positional(mols, atoms="resname ALA and element C",
... position=[(0,0,0), (1,1,1), (2,2,2)])
>>> print(restraints)
PositionalRestraints( size=3
0: PositionalRestraint( 8 => ( 0, 0, 0 ), k=150 kcal mol-1 Å-2 : r0=0 Å )
1: PositionalRestraint( 10 => ( 1, 1, 1 ), k=150 kcal mol-1 Å-2 : r0=0 Å )
2: PositionalRestraint( 14 => ( 2, 2, 2 ), k=150 kcal mol-1 Å-2 : r0=0 Å )
)
Note
The number of positions must equal the number of atoms to be restrained, or equal to 1. If you only pass in a single position then this will be used for all of the restraints.
Sometimes it can be useful to use the same position for all of the restraints, e.g. in the case of a spherical half-harmonic restraint that holds molecules within a spherical bubble, e.g.
>>> center = mols[0].coordinates()
>>> restraints = sr.restraints.positional(mols,
... atoms=mols[f"atoms within 10 of {center}"],
... position=center,
... r0="10 A")
>>> print(restraints)
PositionalRestraints( size=358
0: PositionalRestraint( 0 => ( 16.5471, 4.50102, 15.6589 ), k=150 kcal mol-1 Å-2 : r0=10 Å )
1: PositionalRestraint( 1 => ( 16.5471, 4.50102, 15.6589 ), k=150 kcal mol-1 Å-2 : r0=10 Å )
2: PositionalRestraint( 2 => ( 16.5471, 4.50102, 15.6589 ), k=150 kcal mol-1 Å-2 : r0=10 Å )
3: PositionalRestraint( 3 => ( 16.5471, 4.50102, 15.6589 ), k=150 kcal mol-1 Å-2 : r0=10 Å )
4: PositionalRestraint( 4 => ( 16.5471, 4.50102, 15.6589 ), k=150 kcal mol-1 Å-2 : r0=10 Å )
...
353: PositionalRestraint( 1892 => ( 16.5471, 4.50102, 15.6589 ), k=150 kcal mol-1 Å-2 : r0=10 Å )
354: PositionalRestraint( 1893 => ( 16.5471, 4.50102, 15.6589 ), k=150 kcal mol-1 Å-2 : r0=10 Å )
355: PositionalRestraint( 1906 => ( 16.5471, 4.50102, 15.6589 ), k=150 kcal mol-1 Å-2 : r0=10 Å )
356: PositionalRestraint( 1907 => ( 16.5471, 4.50102, 15.6589 ), k=150 kcal mol-1 Å-2 : r0=10 Å )
357: PositionalRestraint( 1908 => ( 16.5471, 4.50102, 15.6589 ), k=150 kcal mol-1 Å-2 : r0=10 Å )
)
has created half-harmonic restraints that affect all atoms within 10 Å of the center of the first molecule, and that restrain those atoms to be within a sphere of radius 10 Å centered on that molecule.
Note that the restraints only contain the indicies of the atoms that are
restrained, not the atoms themselves. You can get the atoms by looking up
those indicies from the mols molecular container from which the
atoms were selected. For example, here we get all of the atoms that are
subject to that spherical half-harmonic restraint;
>>> print(mols.atoms()[[restraint.atom() for restraint in restraints]])
SireMol::SelectorM<SireMol::Atom>( size=358
0: MolNum(6) Atom( HH31:1 [ 18.45, 3.49, 12.44] )
1: MolNum(6) Atom( CH3:2 [ 18.98, 3.45, 13.39] )
2: MolNum(6) Atom( HH32:3 [ 20.05, 3.63, 13.29] )
3: MolNum(6) Atom( HH33:4 [ 18.80, 2.43, 13.73] )
4: MolNum(6) Atom( C:5 [ 18.48, 4.55, 14.35] )
...
353: MolNum(630) Atom( H1:1893 [ 14.92, 1.28, 17.05] )
354: MolNum(630) Atom( H2:1894 [ 15.19, -0.21, 17.07] )
355: MolNum(623) Atom( O:1907 [ 21.65, 7.88, 9.79] )
356: MolNum(623) Atom( H1:1908 [ 22.33, 8.56, 9.83] )
357: MolNum(623) Atom( H2:1909 [ 21.07, 8.08, 10.53] )
)
Distance or Bond Restraints¶
The sire.restraints.distance() or sire.restraints.bond() functions
are used to create bond or distance restraints. These are identical restraints,
so the functions are just synonyms of each other (they are the same function
with different names).
These functions take the set of molecules or system you want to simulate, plus two search strings, lists of atom indexes, or molecule views holding the pairs of atoms that you want to restrain, e.g.
>>> restraints = sr.restraints.distance(mols, atoms0=0, atoms1=1)
>>> print(restraints)
BondRestraints( size=1
0: BondRestraint( 0 <=> 1, k=150 kcal mol-1 Å-2 : r0=1.09 Å )
)
or
>>> restraints = sr.restraints.bond(mols, atoms0=0, atoms1=1)
>>> print(restraints)
BondRestraints( size=1
0: BondRestraint( 0 <=> 1, k=150 kcal mol-1 Å-2 : r0=1.09 Å )
)
creates a single harmonic distance (or bond) restraint that acts between atoms 0 and 1. By default, the equilibrium distance (r0) is the current distance between the atoms (1.09 Å), and the force constant (k) is 150 kcal mol-1 Å-2.
You can set these via the k and r0 arguments, e.g.
>>> restraints = sr.restraints.bond(mols, atoms0=0, atoms1=1,
... k="100 kcal mol-1 A-2", r0="1.5 A")
>>> print(restraints)
BondRestraints( size=1
0: BondRestraint( 0 <=> 1, k=100 kcal mol-1 Å-2 : r0=1.5 Å )
)
You can specify as many atoms pairs as you like, e.g.
>>> restraints = sr.restraints.bond(mols, atoms0=[0, 1, 2], atoms1=[3, 4, 5])
>>> print(restraints)
BondRestraints( size=3
0: BondRestraint( 0 <=> 3, k=150 kcal mol-1 Å-2 : r0=1.71245 Å )
1: BondRestraint( 1 <=> 4, k=150 kcal mol-1 Å-2 : r0=1.54643 Å )
2: BondRestraint( 2 <=> 5, k=150 kcal mol-1 Å-2 : r0=2.48642 Å )
)
You can specify the atoms using a search string, passing the atoms themselves, or using the atom index, just as you could for the positional restraints.
>>> restraints = sr.restraints.bond(mols,
... atoms0=mols[0][0],
... atoms1=mols["water and element O"])
>>> print(restraints)
BondRestraints( size=630
0: BondRestraint( 0 <=> 22, k=150 kcal mol-1 Å-2 : r0=13.2847 Å )
1: BondRestraint( 0 <=> 25, k=150 kcal mol-1 Å-2 : r0=10.8445 Å )
2: BondRestraint( 0 <=> 28, k=150 kcal mol-1 Å-2 : r0=15.9183 Å )
3: BondRestraint( 0 <=> 31, k=150 kcal mol-1 Å-2 : r0=13.5108 Å )
4: BondRestraint( 0 <=> 34, k=150 kcal mol-1 Å-2 : r0=9.18423 Å )
...
625: BondRestraint( 0 <=> 1897, k=150 kcal mol-1 Å-2 : r0=9.52934 Å )
626: BondRestraint( 0 <=> 1900, k=150 kcal mol-1 Å-2 : r0=12.7207 Å )
627: BondRestraint( 0 <=> 1903, k=150 kcal mol-1 Å-2 : r0=10.8053 Å )
628: BondRestraint( 0 <=> 1906, k=150 kcal mol-1 Å-2 : r0=6.04142 Å )
629: BondRestraint( 0 <=> 1909, k=150 kcal mol-1 Å-2 : r0=17.1035 Å )
)
Note
If the number of atoms in atoms0 and atoms1 are not equal, then
the smaller list will be extended to match by appending multiple copies
of the last atom. In this case, this duplicates the single atom in
atoms0, meaning that this has created distance restraints between
the first atom in the first molecule and the oxygen atoms in all of
the water molecules.
Boresch Restraints¶
The sire.restraints.boresch() function is used to create Boresch restraints,
which are commonly used for restraining the relative positions and orientations
of a receptor and ligand during alchemical absolute binding free energy calculations.
They restrain the six relative external degrees of freedom of the receptor and ligand
by restraining one distance, two angles, and three dihedrals angles which are
defined according to three anchor points in the receptor and three anchor points
in the ligand. For more detail, please see J. Phys. Chem. B 2003, 107, 35, 9535–9551.
To create a Boresch restraint, you need to specify the receptor and ligand anchor atoms (note that the order of the atoms is important). Like the distance restraints, the atoms can be specified using a search string, passing lists of atom indexes, or molecule views holding the atoms. You can also specify the force constants and equilibrium values for the restraints. If not supplied, default force constants of 10 kcal mol-1 Å-2 and 100 kcal mol-1 rad-2 are used for the distance and angle restraints, respectively, and the equilibrium values are set to the current values of the distances and angles in the system supplied. For example,
>>> boresch_restraints = sr.restraints.boresch(mols, receptor=[1574, 1554, 1576], ligand=[4,3,5])
>>> boresch_restraint = boresch_restraints[0]
>>> print(boresch_restraint)
BoreschRestraint( [1574, 1554, 1576] => [4, 3, 5],
k=[10 kcal mol-1 Å-2, 0.0304617 kcal mol-1 °-2, 0.0304617 kcal mol-1 °-2,
0.0304617 kcal mol-1 °-2, 0.0304617 kcal mol-1 °-2, 0.0304617 kcal mol-1 °-2]
r0=15.1197 Å, theta0=[80.5212°, 59.818°],
phi0=[170.562°Ⱐ128.435°Ⱐ192.21°] )
creates a Boresch restraint where the receptor anchor atoms are r1 = 1574, r2 = 1554, and r3 = 1576, and the ligand anchor atoms are l1 = 4, l2 = 3, and l3 = 5. The default force constants have been set and the equilibrium values have been set to the current values of the distances and angles in the system supplied.
Note
Boresch restraints can be subject to instabilities if any three contiguous anchor points approach collinearity (J. Chem. Theory Comput. 2023, 19, 12, 3686–3704). It is important to prevent this by ensuring the associated angles are sufficiently far from 0 or 180 degrees, and that the ktheta force constants are high enough. Sire will raise a warning if the theta0 values are too close to 0 or 180 degrees for the given temperature and force constants.
Alternatively, we could have explicitly set the force constants and equilibrium values, e.g.
>>> boresch_restraints = sr.restraints.boresch(mols, receptor=[1574, 1554, 1576], ligand=[4,3,5],
kr = "6.2012 kcal mol-1 A-2",
ktheta = ["28.7685 kcal mol-1 rad-2", "24.8204 kcal mol-1 rad-2"],
kphi = ["59.8626 kcal mol-1 rad-2", "0.7923 kcal mol-1 rad-2",
"55.1775 kcal mol-1 rad-2"],
r0 = "16 A",
theta0 = ["1.2 rad", "1.3 rad"],
phi0 = ["2.2 rad", "2.5 rad", "1.5 rad"],
name = "boresch_restraint_1")
>>> boresch_restraint = boresch_restraints[0]
>>> print(boresch_restraint)
BoreschRestraint( [1574, 1554, 1576] => [4, 3, 5],
k=[6.2012 kcal mol-1 Å-2, 0.00876339 kcal mol-1 °-2, 0.00756073 kcal mol-1 °-2, 0.0182352 kcal mol-1 °-2, 0.000241348 kcal mol-1 °-2, 0.016808 kcal mol-1 °-2]
r0=16 Å, theta0=[68.7549°, 74.4845°],
phi0=[126.051°Ⱐ143.239°Ⱐ85.9437°] )
Note
sire.restraints.boresch() returns a list of Boresch restraints. If you are only
interested in a single Boresch restraint, you can extract it with the index, e.g.
boresch_restraint = boresch_restraints[0].
When performing an alchemical absolute binding free energy calculation, it is necessary to
calculate the free energy of releasing the decoupled ligand to the standard state volume.
The analytical Boresch correction is almost always accurate if stable restraints have been
selected (see 10.26434/chemrxiv-2023-8s9dz-v2). This can be calculated with
sire.restraints.standard_state_correction():
>>> boresch_restraint = boresch_restraints[0]
>>> from sire import u
>>> correction = sr.restraints.standard_state_correction(boresch_restraint, temperature=u("298 K"))
>>> print(correction)
-6.2399 kcal mol-1
Using restraints in minimisation or dynamics¶
You can use restraints in a minimisation or dynamics simulation by
passing them in via the restraints argument, e.g.
>>> restraints = sr.restraints.positional(mols, atoms="resname ALA and element C")
>>> d = mols.dynamics(timestep="4fs", temperature="25oC",
... restraints=restraints)
>>> d.run("10ps")
>>> mols = d.commit()
runs a dynamics simulation using positional restraints on the carbon atoms
of the ALA residue, while
>>> restraints = sr.restraints.distance(mols, atoms0=mols[0][0],
... atoms1=mols[1][0], r0="5 A")
>>> mols = mols.minimisation(restraints=restraints).run().commit()
performs a minimisation where a distance restraint is applied between the first atoms of the first two molecules, pulling them to be 5 Å apart.
You can pass in a list of restraints, meaning that you can use as many as you wish during a simulation.
>>> pos_rests = sr.restraints.positional(mols, atoms="resname ALA and element C")
>>> dst_rests = sr.restraints.distance(mols, atoms0=mols[0][0], atoms1=mols[1][0])
>>> mols = mols.minimisation(restraints=[pos_rests, dst_rests]).run().commit()
>>> d = mols.dynamics(timestep="4fs", temperature="25oC",
... restraints=[pos_rests, dst_rests])
>>> d.run("10ps")
>>> mols = d.commit()
Note
It is a good idea to run minimisation before dynamics whenever you add
restraints to a system. This is because the restraints could put a lot
of energy into the system, causing blow-ups or NaN errors.
Using restraints in the low-level API¶
You can use restraints in the low-level API by passing them in as a
map parameter using the key restraints, e.g.
>>> omm = sr.convert.to(mols, "openmm",
... map={"restraints": [pos_rests, dst_rests]})
>>> print(omm)
openmm::Context( num_atoms=1915 integrator=VerletIntegrator timestep=1.0 fs platform=HIP )
More information about the mappable options for the low-level API can be found in the detailed guide.
Fixing atoms in place¶
As well as restraining atoms, you also have the option of fixing atoms in space during the simulation. Atoms which are fixed are not moved by minimisation or the dynamics integrator. This can be useful if you, e.g. want to freeze atoms outside a radius of a ligand binding site, or if you want to minimise the solvent while keeping the solute fixed.
You can fix atoms by passing in the set of atoms to fix as the
fixed argument to the minimisation() or
dynamics() functions, e.g.
>>> mols = sr.load(sr.expand(sr.tutorial_url, "ala.top", "ala.crd"))
>>> d = mols.dynamics(timestep="4fs", temperature="25oC",
... fixed=mols[0])
>>> d.run("10ps")
>>> mols = d.commit()
will run dynamics where all of the atoms in mols[0] (the solute)
are fixed. You can pass in a molecular container containing all of the
atoms that should be fixed, or a search string, or the atom indexes
of the atoms. Here,
>>> d = mols.dynamics(timestep="1fs", temperature="25oC",
... fixed="element C")
>>> d.run("10ps")
>>> mols = d.commit()
we have run dynamics where all of the carbon atoms are fixed, while
>>> d = mols.dynamics(timestep="1fs", temperature="25oC",
... fixed=[0, 2, 4, 6, 8])
>>> d.run("10ps")
>>> mols = d.commit()
would run dynamics where atoms at indicies 0, 2, 4, 6 and 8 are fixed.
You could even use search strings to fix atoms by distances. If you do this, it is best to also add half-harmonic positional restraints that hold nearby molecules within the sphere of mobile atoms, e.g.
>>> center = mols[0].coordinates()
>>> radius = "10 A"
>>> restraints = sr.restraints.positional(mols,
... atoms=mols[f"molecules within 10 of {center}"],
... position=center,
... r0=radius)
>>> d = mols.dynamics(timestep="1fs", temperature="25oC",
... fixed=f"not molecules within {radius} of {center}")
>>> d.run("10ps")
>>> mols = d.commit()
This implements the restraints in OpenMM by setting the masses of the fixed atoms to zero. This signals the OpenMM integrator to skip over these atoms and not move them.
Note
Be careful using constraints with fixed atoms. If a constraint involves
a fixed atom and a mobile atom, then there is a high risk that the
constraint won’t be able to be satisfied during dynamics, and
Particle coordinate is NaN error or similar will occur.
It it safest to use a small timestep (e.g. 1 fs) when constraining
only parts of molecules.
Note
While the fixed atoms are not moved by the integrator, they are still included in the energy calculation. This means that the computational cost of the simulation will still be high, and also that energies will not be conserved. Make sure you use an integrator that provides a heat bath (e.g. NVT or NPT integrators).