# 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
>>> 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 the `bond_length` and `bond_k` levers.

• `angle`: OpenMM::HarmonicAngleForce. This models all of the angles between atoms in the system. It uses parameters that are controlled by the `angle_size` and `angle_k` levers.

• `torsion`: OpenMM::PeriodicTorsionForce. This models all of the torsions (dihedrals and impropers) in the system. It uses parameters that are controlled by the `torsion_phase` and `torsion_k` levers.

• `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 the `charge`, `sigma`, `epsilon`, `charge_scale` and `lj_scale` levers.

• `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 the `charge`, `sigma`, `epsilon`, `alpha` and `kappa` levers.

• `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 the `charge`, `sigma`, `epsilon`, `alpha` and `kappa` levers.

• `ghost-14`: OpenMM::CustomBondForce. This models all of the 1-4 non-bonded interactions involving ghost atoms. It uses parameters that are controlled by the `charge`, `sigma`, `epsilon`, `alpha`, `kappa`, `charge_scale` and `lj_scale` levers.

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