The 2023.4 release introduces support for GPU-accelerated molecular dynamics alchemical free energy simulations. This builds on the integration with OpenMM introduced in the 2023.2 release, and the support for trajectories and units introduced in 2023.3.
In this post, we will show how you can run a relative hydration free energy simulation, and then process the results using alchemlyb.
The concept of a merged molecule is central to the way that free energy calculations are implemented in sire. A merged molecule is one that represents both a “reference” state and a “perturbed” state. These are the two states that the free energy simulation will morph between, and for which the free energy difference will be calculated.
For example, here we have pre-prepared a merged molecule that represents the perturbation from ethane to methanol.
import sire as sr
mols = sr.load(sr.expand(sr.tutorial_url, "merged_molecule.s3"))
This system contains a single merged molecule in a box of water. Merged molecules are idenfitied by the molecule property is_perturbable, which will be True. We can extract the merged molecule from this system using
mol = mols["molecule property is_perturbable"]
print(mol)
Molecule( Merged_Molecule:6 num_atoms=8 num_residues=1 )
Merged molecules contain two sets of the molecular properties; one that represents the reference state, and one that represents the perturbed state. These are identified by the 0 and 1 suffixes.
For example, the reference state atomic charges are in the “charge0” property;
print(mol.property("charge0"))
SireMol::AtomCharges( size=8 0: -0.09435 |e| 1: -0.09435 |e| 2: 0.03145 |e| 3: 0.03145 |e| 4: 0.03145 |e| 5: 0.03145 |e| 6: 0.03145 |e| 7: 0.03145 |e| )
while the perturbed state atomic charges are in the “charge1” property:
print(mol.property("charge1"))
SireMol::AtomCharges( size=8 0: -0.5988 |e| 1: 0.1167 |e| 2: 0 |e| 3: 0.396 |e| 4: 0 |e| 5: 0.0287 |e| 6: 0.0287 |e| 7: 0.0287 |e| )
We can view the perturbed state by linking to its properties, e.g.
mol = mol.perturbation().link_to_perturbed().commit()
mol["not element Xx"].view()
(noting to only view non-dummy atoms via mol["not element Xx"])
We can view the reference state using
mol = mol.perturbation().link_to_reference().commit()
mol.view()
A λ-coordinate is used to morph from the reference state (at λ=0) to the perturbed state (at λ=1). We can run dynamics at any λ-value just by passing this in as an argument to the minimisation and dynamics functions.
mol = mol.minimisation(lambda_value=0.5).run().commit()
d = mol.dynamics(lambda_value=0.5, temperature="25oC")
d.run("10ps")
print(d)
Dynamics(completed=10 ps, energy=16.3839 kcal mol-1, speed=179.5 ns day-1)
The next step is to calculate and store the energy during the trajectory of the molecules as a function of λ. We do this by creating an EnergyTrajectory. We do this by telling the dynamics simulation to save the energy periodically, via the energy_frequency argument.
d = mol.dynamics(lambda_value=0.5, temperature="25oC")
d.run("10ps", energy_frequency="0.1ps")
print(d)
Dynamics(completed=10 ps, energy=16.1006 kcal mol-1, speed=177.8 ns day-1)
The energy trajectory is retrieved via the energy_trajectory() function.
t = d.energy_trajectory()
print(t)
EnergyTrajectory( size=100 time lambda 0.5 kinetic potential 0.1 0.5 2.76584 1.08741 2.76584 0.2 0.5 3.93688 1.25623 3.93688 0.3 0.5 4.28899 1.4863 4.28899 0.4 0.5 4.31009 2.95505 4.31009 0.5 0.5 4.67531 3.20987 4.67531 ... 9.6 0.5 6.27478 6.38104 6.27478 9.7 0.5 6.17797 7.88701 6.17797 9.8 0.5 9.33133 6.4523 9.33133 9.9 0.5 5.75732 10.3856 5.75732 10 0.5 8.44169 7.65892 8.44169 )
We calculate the free energy across λ by collecting and averaging the energy across many dynamics simulations run across λ. For this to work, the energy at neighbouring λ-values has to be evaluated in addition to the energy for the simulated λ-value. The lambda_windows argument lets us tell the simulation to evalute the energy at extra λ-values during the trajectory.
For example, let’s run the simulation at λ=0.5, while also calculating the energy at λ=0 and λ=1.
d = mol.dynamics(lambda_value=0.5, temperature="25oC")
d.run("10ps", energy_frequency="0.1ps", lambda_windows=[0, 1])
print(d)
t0_5 = d.energy_trajectory()
print(t0_5)
Dynamics(completed=10 ps, energy=14.9342 kcal mol-1, speed=146.8 ns day-1) EnergyTrajectory( size=100 time lambda 0 0.5 1 kinetic potential 0.1 0.5 4.4009 2.84375 6.48395 1.04331 2.84375 0.2 0.5 4.97366 3.27012 6.74341 2.64619 3.27012 0.3 0.5 6.87626 3.87773 5.86002 2.73733 3.87773 0.4 0.5 7.04107 4.16728 6.37846 3.85838 4.16728 0.5 0.5 7.40754 5.29728 8.20447 4.60734 5.29728 ... 9.6 0.5 14.6052 9.4928 9.89531 6.19703 9.4928 9.7 0.5 10.775 9.22934 12.5969 6.46628 9.22934 9.8 0.5 12.0069 8.82677 10.7986 7.36023 8.82677 9.9 0.5 13.6869 9.652 10.5415 7.06599 9.652 10 0.5 12.1086 10.1357 13.1283 4.79852 10.1357 )
To calculate a free energy, we would also need to run simulations at λ=0 and λ=1…
d = mol.dynamics(lambda_value=0, temperature="25oC")
d.run("10ps", energy_frequency="0.1ps", lambda_windows=[0.5, 1])
t0_0 = d.energy_trajectory()
d = mol.dynamics(lambda_value=1, temperature="25oC")
d.run("10ps", energy_frequency="0.1ps", lambda_windows=[0, 0.5])
t1_0 = d.energy_trajectory()
alchemlyb is a great tool that can process a table of energy values across λ, and compute a free energy with associated error.
First, we combine the three EnergyTrajectory objects to into a single pandas DataFrame in alchemlyb format.
df = sr.morph.to_alchemlyb([t0_0, t0_5, t1_0])
print(df)
0.0 0.5 1.0 time fep-lambda 0.1 0.0 2.685941 7.224478 17.266447 0.2 0.0 2.265553 6.088665 15.488273 0.3 0.0 3.428572 6.842974 15.979260 0.4 0.0 4.832400 10.716389 22.198475 0.5 0.0 4.741605 5.834284 12.575234 ... ... ... ... 9.6 1.0 14.166821 9.616899 10.351943 9.7 1.0 14.249995 8.776641 8.706912 9.8 1.0 18.603194 11.711414 9.567702 9.9 1.0 15.942173 10.188642 9.893304 10.0 1.0 15.318865 9.876553 10.382260 [300 rows x 3 columns]
This DataFrame can be passed directly into alchemlyb to calculate the relative free energy using the BAR method.
from alchemlyb.estimators import BAR
b = BAR()
b.fit(df)
print(b.delta_f_.loc[0.00, 1.00])
3.029635052005364
One of the nice things about alchemlyb is that it is easy to use different free energy estimators to calculate the free energy. Let’s try the more accurate MBAR method.
from alchemlyb.estimators import MBAR
b = MBAR()
b.fit(df)
print(b.delta_f_.loc[0.00, 1.00])
3.0231274550384573
Not bad – both methods agree that the free energy calculated from these simulations is 3.0 kcal mol-1. This is only a rough estimate of the relative free energy, as we only used three λ-windows, and only ran very short (10 ps) dynamics simulations. A better estimate could be made by running more λ-windows and running longer dynamics simulations at each window.
This script let’s you run longer simulations across more λ-windows, running the perturbation both for ethane to methanol in the gas phase and also in a box of water. Running 250 ps of dynamics across 21 evenly-spaced λ-windows gives a gas-phase relative free energy of 2.98 +/- 0.01 kcal mol-1. Our value calculated above using only three λ-windows and short simulations is remarkably close!
Combining this with the water-phase relative free energy of -3.30 +/- 0.02 kcal mol-1 gives a predicted relative hydration free energy of ethane and methanol of -6.28 +/- 0.02 kcal mol-1. This is in excellent agreement with values computed with other codes.
Now that we know that the free energy code works, our focus now is checking accuracy and performance for protein-ligand relative and absolute binding free energy simulations. To do this, we’ve implemented support for restraints. These restraints can be morphed using λ, via our new LambdaSchedule objects. This schedule gives you control over how forcefield parameters and restraints are changed as a function of λ. For example, here we show how you can set up a schedule that slowly turns on some restraints before morphing from the reference molecule to the perturbed molecule. Then, after the morph, these restraints are slowly turned off.
We think that this framework, and the concepts of merged molecules and λ-schedules gives a lot of power to users and developers of free energy simulations. We’re looking forward to exploring this new functionality ourselves, and also seeing what everyone else in the community creates.
As with all of our diaries, if you want to try this yourself, please feel free connect to try.openbiosim.org and starting a notebook. You can download the notebook used to generate this post onto the server by running this command in one of the notebook code cells.
! wget https://github.com/OpenBioSim/posts/raw/main/sire/004_alchemy/alchemy.ipynb
Have a play and let us know what you think.
