Lambda schedules and levers#
In the last chapter we saw how
we could use a LambdaSchedule to control how the
forcefield parameters of perturbable molecules are morphed as a function
of λ.
For example, let’s load up a merged molecule system that represents
neopentane to methane in water, and then explore the
LambdaSchedule that is associated with this system.
>>> import sire as sr
>>> mols = sr.load_test_files("neo_meth_solv.bss")
>>> mols = sr.morph.link_to_reference(mols)
>>> print(mols)
System( name=BioSimSpace_System num_molecules=879 num_residues=879 num_atoms=2651 )
>>> d = mols.dynamics()
>>> s = d.get_schedule()
>>> print(s)
LambdaSchedule(
morph: (-λ + 1) * initial + λ * final
)
This shows that all perturbable parameters of all perturbable molecules in this system will be morphed in a single stage (called “morph”), using a linear interpolation between the initial and final values of the parameters.
As we saw in the last section, we can find the exact values of all of the perturbable parameters of a perturbable molecule via the perturbation object.
>>> p = mols[0].perturbation()
>>> print(p)
Perturbation( Molecule( Merged_Molecule:6 num_atoms=17 num_residues=1 ) )
>>> p_omm = p.to_openmm()
>>> print(p_omm.changed_bonds())
bond length0 length1 k0 k1
0 C2:2-C4:4 0.15375 0.10969 251793.12 276646.08
Note
In this case, we know that only the first molecule is perturbable,
hence why it is safe to use mols[0].perturbation(). In general,
you would need to find perturbable molecule(s), e.g. using
pert_mols = mols.molecules("property is_perturbable").
From this, we can see that the bond between atoms C2 and C4 is perturbable, and the above schedule will morph the bond length from 0.15375 nm to 0.10969 nm, and the force constant from 251793.12 kJ mol-1 nm-2 to 276646.08 kJ mol-1 nm-2, linearly with respect to λ.
Note
The parameters are directly as would be used in an OpenMM force, i.e. in OpenMM default units of nanometers and kilojoules per mole.
Controlling individual levers#
We can also control how individual parameters in individual forces are morphed. We call these individual morphing controls “levers”. So, the bond potential has two levers, its bond length, and its force constant.
We can list the levers that are available in the OpenMM context using the
get_levers() function.
>>> print(s.get_levers())
['charge', 'sigma', 'epsilon', 'charge_scale', 'lj_scale',
'bond_length', 'bond_k', 'angle_size', 'angle_k', 'torsion_phase',
'torsion_k', 'alpha', 'kappa']
The parameter representing bond length is connected to the bond_length lever,
while the parameter representing the bond force constant is
connected to the bond_k lever.
We can set the morphing equation for individual levers by naming the lever
in the set_equation() function.
>>> l = s.lam()
>>> init = s.initial()
>>> fin = s.final()
>>> s.set_equation(stage="morph", lever="bond_length",
... equation=(1-l**2)*init + l**2*fin)
>>> print(s)
LambdaSchedule(
morph: initial * (-λ + 1) + final * λ
bond_length: initial * (-λ^2 + 1) + final * λ^2
)
Note
We extracted the symbols representing λ and the initial and final
states to l, init and fin, just to make it easier to
write the equations.
We can see that the bond_length lever in the morph stage is now
interpolated from the initial to final value by λ^2, rather than λ.
All of the other levers continue to use the default equation for this stage, which is the linear interpolation between the initial and final values.
You can change the default equation used for a stage using the
set_default_equation() function, e.g.
>>> s.set_default_equation(stage="morph", equation=(1-l)*init + l*fin)
>>> print(s)
LambdaSchedule(
morph: initial * (-λ + 1) + final * λ
bond_length: initial * (-λ^2 + 1) + final * λ^2
)
Controlling individual levers in individual forces#
Multiple OpenMM Force objects are combined in the OpenMM context to
model the total force acting on each atom in the system. OpenMM is very
flexible, and supports the arbitrary combination of lots of different
Force objects. In sire, we use a simple collection of Force objects
that, when combined, model perturbable systems. You can list the names
of the Force objects used via the get_forces()
function.
>>> print(s.get_forces())
['clj', 'bond', 'angle', 'torsion', 'ghost/ghost',
'ghost/non-ghost', 'ghost-14']
In this case, as we have a perturbable system, the Force objects used are;
bond: OpenMM::HarmonicBondForce. This models all of the bonds between atoms in the system. It uses parameters that are controlled by thebond_lengthandbond_klevers.angle: OpenMM::HarmonicAngleForce. This models all of the angles between atoms in the system. It uses parameters that are controlled by theangle_sizeandangle_klevers.torsion: OpenMM::PeriodicTorsionForce. This models all of the torsions (dihedrals and impropers) in the system. It uses parameters that are controlled by thetorsion_phaseandtorsion_klevers.clj: OpenMM::NonbondedForce. This models all of the electrostatic (coulomb) and van der Waals (Lennard-Jones) interactions between non-ghost atoms in the system. Non-ghost atoms are any atoms that are not ghosts in either end state. It uses parameters that are controlled by thecharge,sigma,epsilon,charge_scaleandlj_scalelevers.ghost/ghost: OpenMM::CustomNonbondedForce. This models all of the electrostatic (coulomb) and van der Waals (Lennard-Jones) interactions between ghost atoms in the system. Ghost atoms are any atoms that are ghosts in either end state. It uses parameters that are controlled by thecharge,sigma,epsilon,alphaandkappalevers.ghost/non-ghost: OpenMM::CustomNonbondedForce. This models all of the electrostatic (coulomb) and van der Waals (Lennard-Jones) interactions between the ghost atoms and the non-ghost atoms in the system. It uses parameters that are controlled by thecharge,sigma,epsilon,alphaandkappalevers.ghost-14: OpenMM::CustomBondForce. This models all of the 1-4 non-bonded interactions involving ghost atoms. It uses parameters that are controlled by thecharge,sigma,epsilon,alpha,kappa,charge_scaleandlj_scalelevers.
Some levers, like bond_length, are used only by a single Force object.
However, others, like charge, are used by multiple Force objects.
By default, setting a lever will affect the parameters in all of the Force
objects that use that lever. However, you can limit which Force objects
are affected by specifying the force in the set_equation()
function.
>>> s.set_equation(stage="morph", force="ghost/ghost", lever="alpha",
equation=0.5*s.get_equation("morph"))
>>> print(s)
LambdaSchedule(
morph: initial * (-λ + 1) + final * λ
bond_length: (-λ^2 + 1) * initial + final * λ^2
ghost/ghost::alpha: 0.5 * (initial * (-λ + 1) + final * λ)
)
Here, we have set the alpha lever in the ghost/ghost Force object
to set the alpha parameter to equal half of its linearly interpolated
value.
Note
The alpha parameter controls the amount of softening used in the
soft-core potential for modelling ghost atoms. An alpha value of
0.0 means that the soft-core potential is not used, while an alpha
value of 1.0 means that the soft-core potential is on and strong.
Scaling up alpha will gradually soften any ghost atoms.
Controlling individual levers for individual molecules#
We can also control how individual levers for individual forces are morphed for individual perturbable molecules in the system. This is useful if you have multiple perturbable molecules, and you want to control how each one perturbs separately.
To do this, we use the set_molecule_schedule() function
to set the schedule for a specific perturbable molecule.
First, let’s get the original schedule for our simulation…
>>> orig_s = d.get_schedule()
>>> print(orig_s)
LambdaSchedule(
morph: initial * (-λ + 1) + final * λ
)
Now, let’s set the schedule to be used only for the first perturbable molecule in the system to the custom one we created earlier.
>>> orig_s.set_molecule_schedule(0, s)
>>> print(orig_s)
LambdaSchedule(
morph: initial * (-λ + 1) + final * λ
Molecule schedules:
0: LambdaSchedule(
morph: initial * (-λ + 1) + final * λ
bond_length: (-λ^2 + 1) * initial + final * λ^2
ghost/ghost::alpha: 0.5 * (initial * (-λ + 1) + final * λ)
)
)
This shows that the default for all perturbable molecules except the first is to use the default morph equation for all levers in all forces.
However, for the first perturbable molecule (which has index 0),
this uses our custom equation for the bond_length lever in the
morph stage, and our custom equation for the alpha lever in
the ghost/ghost force in the morph stage.
Once you are happy, we can set the schedule to be used for the simulaton
via the set_schedule() function.
>>> d.set_schedule(orig_s)
>>> print(d.get_schedule())
LambdaSchedule(
morph: initial * (-λ + 1) + final * λ
Molecule schedules:
0: LambdaSchedule(
morph: initial * (-λ + 1) + final * λ
bond_length: (-λ^2 + 1) * initial + final * λ^2
ghost/ghost::alpha: 0.5 * (initial * (-λ + 1) + final * λ)
)
)
Viewing the effect of levers on a merged molecule#
You can view the effect of the LambdaSchedule on a
the PerturbableOpenMMMolecule using
the get_lever_values()
function.
>>> df = p_omm.get_lever_values(schedule=orig_s)
>>> print(df)
clj-charge-1 clj-charge-2 clj-charge-4 ... bond-bond_length-3 angle-angle_k-19 angle-angle_size-19
λ ...
0.00 -0.085335 -0.060235 -0.085335 ... 0.10969 387.438400 1.916372
0.01 -0.084482 -0.059362 -0.085566 ... 0.10969 386.861008 1.915985
0.02 -0.083629 -0.058489 -0.085797 ... 0.10969 386.283616 1.915597
0.03 -0.082775 -0.057615 -0.086027 ... 0.10969 385.706224 1.915210
0.04 -0.081922 -0.056742 -0.086258 ... 0.10969 385.128832 1.914822
... ... ... ... ... ... ... ...
0.96 -0.003413 0.023607 -0.107477 ... 0.10969 332.008768 1.879176
0.97 -0.002560 0.024480 -0.107708 ... 0.10969 331.431376 1.878788
0.98 -0.001707 0.025353 -0.107939 ... 0.10969 330.853984 1.878401
0.99 -0.000853 0.026227 -0.108169 ... 0.10969 330.276592 1.878013
1.00 0.000000 0.027100 -0.108400 ... 0.10969 329.699200 1.877626
[101 rows x 19 columns]
It can be useful to plot these, to check that the morphing is as expected.
>>> ax = df.plot()
>>> ax.legend(bbox_to_anchor=(1.0, 1.0))
Note
The line ax.legend(bbox_to_anchor=(1.0, 1.0)) is used to move the
legend outside of the plot area, so that it doesn’t obscure the data.
Unfortunately, the large values of the bond force constant make it very
difficult to see the effect of the LambdaSchedule on the
other parameters. To fix this, lets filter out any columns that contain
values that have an absolute value greater than 5.
>>> skip_columns = df.abs().gt(5).apply(lambda x: x.index[x].tolist(), axis=1)[0]
>>> ax = df.loc[:, ~df.columns.isin(skip_columns)].plot()
>>> ax.legend(bbox_to_anchor=(1.0, 1.0))
