Perturbation#
As seen in chapter 6, the concept of a
“merged molecule” is central to how sire implements the
perturbations (morphing) needed for free energy calculations.
A merged molecule is one that has been created with a single set of atoms, but with properties that represent either a “reference” (λ=0) or “perturbed” (λ=1) state.
For example, let’s load a merged molecule that represents neopentane in the reference state, and methane in the perturbed state.
>>> import sire as sr
>>> mols = sr.load_test_files("neo_meth_scratch.bss")
>>> print(mols)
System( name=BioSimSpace_System num_molecules=1 num_residues=1 num_atoms=17 )
>>> print(mols[0])
Molecule( Merged_Molecule:6 num_atoms=17 num_residues=1 )
Note
This molecule was created using BioSimSpace. We will cover how to create merged molecules in a later section in this tutorial.
Examining the perturbation#
Normally, a sire molecule will have one set of forcefield parameters.
For example, the atomic charges would be stored in the charge property,
the Lennard-Jones parameters in the LJ property, the bond parameters in
the bond property, etc.
However, a merged molecule has two sets of parameters, one for the reference
state (representing λ=0) and one for the perturbed state (representing λ=1).
The reference state parameters are stored in the 0 properties, e.g.
charge0, LJ0, bond0, etc. The perturbed state parameters are
stored in the 1 properties, e.g. charge1, LJ1, bond1, etc.
We can see this for our merged molecule above by printing out all of the property keys.
>>> print(mols[0].property_keys())
['gb_radii0', 'gb_radii1', 'atomtype0', 'element0', 'atomtype1', 'element1',
'time', 'intrascale0', 'intrascale1', 'improper0', 'mass0', 'improper1',
'mass1', 'ambertype0', 'ambertype1', 'treechain0', 'treechain1',
'parameters0', 'parameters1', 'space', 'charge0', 'charge1',
'connectivity', 'dihedral0', 'forcefield0', 'dihedral1', 'coordinates0',
'forcefield1', 'coordinates1', 'LJ0', 'name0', 'LJ1', 'name1', 'angle0',
'angle1', 'bond0', 'bond1', 'molecule0', 'molecule1', 'gb_screening0',
'gb_screening1', 'is_perturbable']
Notice how there is a 0 and 1 version of nearly every property.
These represent the molecule at the reference and perturbed end states, e.g.
>>> print(mols[0].property("charge0"))
SireMol::AtomCharges( size=17
0: -0.0853353 |e|
1: -0.0602353 |e|
2: -0.0853353 |e|
3: -0.0853353 |e|
4: -0.0853353 |e|
...
12: 0.0334647 |e|
13: 0.0334647 |e|
14: 0.0334647 |e|
15: 0.0334647 |e|
16: 0.0334647 |e|
)
>>> print(mols[0].property("charge1"))
SireMol::AtomCharges( size=17
0: 0 |e|
1: 0.0271 |e|
2: 0 |e|
3: -0.1084 |e|
4: 0 |e|
...
12: 0.0271 |e|
13: 0.0271 |e|
14: 0 |e|
15: 0 |e|
16: 0 |e|
)
prints the charges for neopentane, and then methane. Note how many of the
atom charges for methane are zero. This is because methane has fewer
atoms than neopentane, and so the charges (and other parameters) for the
extra atoms are “switched off” for the perturbed state. These extra atoms
are often called “dummy atoms”. In sire, we prefer to call them
“ghost atoms”. They are “ghosts” because they fade away to nothing as
λ increases from 0 to 1.
In addition to the 0 and 1 properties, a perturbable molecule must
also contain the is_perturbable property. This is a boolean that signals
whether or not a molecule is perturbable. It should be True if the molecule
can be perturbed.
>>> print(mols[0].property("is_perturbable"))
True
Perturbation objects#
Examining the perturbable properties directly can be a little cumbersome.
To make things easier, there are a number of helper classes and functions.
These are accessed via the sire.morph module, and the
perturbation() function of a Molecule.
Let’s create the Perturbation object for our molecule.
>>> pert = mols[0].perturbation()
>>> print(pert)
Perturbation( Molecule( Merged_Molecule:6 num_atoms=17 num_residues=1 ) )
We can extract a molecule that contains only the reference or perturbed
parameters using the extract_reference() and
extract_perturbed() functions.
>>> ref = pert.extract_reference()
>>> print(ref.property_keys())
['bond', 'atomtype', 'time', 'intrascale', 'element', 'dihedral',
'parameters', 'ambertype', 'angle', 'gb_radii', 'treechain', 'space',
'forcefield', 'name', 'connectivity', 'gb_screening', 'molecule',
'coordinates', 'LJ', 'charge', 'mass', 'improper']
>>> print(ref.property("charge"))
SireMol::AtomCharges( size=17
0: -0.0853353 |e|
1: -0.0602353 |e|
2: -0.0853353 |e|
3: -0.0853353 |e|
4: -0.0853353 |e|
...
12: 0.0334647 |e|
13: 0.0334647 |e|
14: 0.0334647 |e|
15: 0.0334647 |e|
16: 0.0334647 |e|
)
Note
You can extract the reference and perturbed molecules in a collection
using the sire.morph.extract_reference() and
sire.morph.extract_perturbed() functions. For example,
mols = sire.morph.extract_reference(mols) would extract the
reference state of all molecules in mols.
Note
By default, ghost atoms will be removed when you extract an end state.
You can retain ghost atoms by passing remove_ghosts=False to the
above extract_reference and extract_perturbed functions.
Extracting the reference or perturbed states can be useful if you want to
go back to the unmerged molecule, e.g. for visualisation via the
view() function. However, normally we would want to
keep the properties of the two end states, and then choose one end state
as the “current” state. We can do this by creating links from the
“standard” property names (e.g. charge, LJ etc.) to the equivalent
properties of the chosen end state. You could do this manually, but it is
much easier to use the link_to_reference()
and link_to_perturbed() functions.
For example, here we will link to the perturbed state.
>>> mol = pert.link_to_perturbed()
>>> print(mol.get_links())
{'improper': 'improper1', 'gb_screening': 'gb_screening1',
'mass': 'mass1', 'dihedral': 'dihedral1', 'parameters': 'parameters1',
'treechain': 'treechain1', 'bond': 'bond1', 'ambertype': 'ambertype1',
'molecule': 'molecule1', 'atomtype': 'atomtype1', 'charge': 'charge1',
'angle': 'angle1', 'forcefield': 'forcefield1',
'coordinates': 'coordinates1', 'intrascale': 'intrascale1',
'name': 'name1', 'LJ': 'LJ1', 'element': 'element1',
'gb_radii': 'gb_radii1'}
>>> print(mol.property("charge"))
SireMol::AtomCharges( size=17
0: 0 |e|
1: 0.0271 |e|
2: 0 |e|
3: -0.1084 |e|
4: 0 |e|
...
12: 0.0271 |e|
13: 0.0271 |e|
14: 0 |e|
15: 0 |e|
16: 0 |e|
)
Note
You can link the reference or perturbed states in a collection of
molecules using the sire.morph.link_to_reference() and
sire.morph.link_to_perturbed() functions. For example,
mols = sire.morph.link_to_perturbed(mols) would link all molecules
in mols to the perturbed state.
Note
It is a good idea when loading a system containing one or more merged
molecules to decide on which end state you want to use as the “current”.
After loading, simply call either
mols = sire.morph.link_to_reference(mols) or
mols = sire.morph.link_to_perturbed(mols) to update the molecules
with your chosen state.
Inspecting the changing parameters#
Under the hood, the properties of the two end states are converted into
parameters for the underlying OpenMM system used for the
free energy simulations. You can access these parameters by converting
the above perturbation into a
PerturbableOpenMMMolecule.
>>> pert_omm = pert.to_openmm()
>>> print(pert_omm)
PerturbableOpenMMMolecule()
This object contains all of the parameters needed to represent both
end states of this molecule in the OpenMM forces. For example,
get_charges0()
returns a list of the charges for the reference
state, while
get_charges0()
returns a list of the charges for the perturbed state.
>>> print(pert_omm.get_charges0())
[-0.08533529411764705, -0.06023529411764705, -0.08533529411764705,
-0.08533529411764705, -0.08533529411764705, 0.03346470588235294,
0.03346470588235294, 0.03346470588235294, 0.03346470588235294,
0.03346470588235294, 0.03346470588235294, 0.03346470588235294,
0.03346470588235294, 0.03346470588235294, 0.03346470588235294,
0.03346470588235294, 0.03346470588235294]
At this level, we are most interested in the parameters that change as
we morph from the reference to the perturbed state (i.e. as we move
from λ=0 to λ=1). You can access these via the
changed_atoms(),
changed_bonds(),
changed_angles(),
changed_torsions(),
changed_exceptions() and
changed_constraints()
functions.
>>> print(pert_omm.changed_bonds())
bond length0 length1 k0 k1
0 C2:2-C4:4 0.15375 0.10969 251793.12 276646.08
>>> print(pert_omm.changed_angles())
angle size0 size1 k0 k1
0 C2:2-C4:4-H12:12 1.916372 1.877626 387.4384 329.6992
1 C2:2-C4:4-H14:14 1.916372 1.877626 387.4384 329.6992
2 C2:2-C4:4-H13:13 1.916372 1.877626 387.4384 329.6992
Note
The parameters are directly as would be used in an OpenMM force, i.e. in OpenMM default units.
Here we see that the perturbation involves the C2-C4 bond changing from 0.15375 nm to 0.10969 nm, with an associated change in the force constant from 251793.12 kJ mol^-1 nm^-2 to 276646.08 kJ mol^-1 nm^-2.
Similarly, the C2-C4-H12 angle changes from 1.916372 radians to 1.877626 radians, with an associated change in the force constant from 387.4384 kJ mol^-1 rad^-2 to 329.6992 kJ mol^-1 rad^-2.
These functions return their output as pandas dataframes. You can get a raw
output by passing in to_pandas=False.
>>> print(pert_omm.changed_bonds(to_pandas=False))
[(Bond( C2:2 => C4:4 ), 0.15375000000000003, 0.10969000000000001, 251793.11999999997, 276646.08)]
As well as forcefield parameters, you can also access any changes in
constraint parameters caused by constraining perturbable bonds. For example,
here we can create the PerturbableOpenMMMolecule
used when the bonds constraint algorithm is used.
>>> pert_omm = pert.to_openmm(constraint="bonds")
Now, we can see how the constraint parameters will change across λ using
the changed_constraints()
function.
>>> print(pert_omm.changed_constraints())
atompair length0 length1
0 C2:2-C4:4 0.15375 0.10969
In this case, we see that the perturbing C2-C4 bond is constrained, with a constraint length of 0.15375 nm in the reference state, and 0.10969 nm in the perturbed state.
Visualising the perturbation#
Perturbations can involve changes in bond lengths, or angle / torsion sizes.
These can be difficult to visualise from the raw numbers. To help with this,
the view_reference() and
view_perturbed() functions can be used to
view a 3D movie of the perturbation from either the reference or perturbed
states.
>>> pert.view_reference()
Note
The movie loops from λ=0 to λ=1, and then back to λ=0. You can pass in
any of the visualisation options as used in the standard
view() function. The viewer may show some
bonds as broken - this is just because they are too long to be shown
when calculated at λ=0.