Alanine-dipeptide conformational landscape

Note

The code in this tutorial was adapted from FastMBAR.

In a recent preprint we used the emle-engine interface to sander to compute free-energy surfaces for alanine-dipeptide as a function of the Φ and Ψ dihedral angles shown below. Compared to regular mechanical embedding, EMLE was found to be closer to the reference density-functional theory (DFT) surface. In this tutorial we will show how to use the sire-emle interface to set up and run the same calculations using OpenMM.

Creating a context with sire-emle

As in the previous section, we can first use sire-emle to create a QM/MM capable dynamics object for the alanine-dipeptide example system. We can then extract the underlying OpenMM context from this.

First we will create an EMLECalculator to compute the QM intramolecular interaction using ANI-2x along with the electrostatic embedding interaction.

>>> from emle import EMLECalculator
>>> calculator = EMLECalculator(device="cpu")

Next we will load the alanine-dipeptide system using sire:

>>> import sire as sr
>>> mols = sr.load_test_files("ala.crd", "ala.top")

We can then create an EMLEEngine that can be be used to perform QM/MM:

>>> qm_mols, engine = sr.qm.emle(mols, mols[0], calculator)

Here the first argument is the molecules that we are simulating, the second selection coresponding to the QM region (here this is the first molecule), and the third is calculator that was created above. The fourth and fifth arguments are optional, and specify the QM cutoff distance and the neigbour list update frequency respectively. (Shown are the default values.) The function returns a modified version of the molecules containing a “merged” dipeptide that can be interpolated between MM and QM levels of theory, along with an engine. The engine registers a Python callback that uses emle-engine to perform the QM calculation.

We can now create a dynamics that will create an OpenMM context for us and can be used to run a simulation:

>>> d = mols.dynamics(
...     timestep="1fs",
...     constraint="none",
...     engine=engine,
...     platform="cpu",
... )

Before extracting the context we will use the dynamics object to minimise the alanine-dipeptide system:

>>> d.minimise()

Setting up umbrella sampling with OpenMM

We can now extract the underlying OpenMM context from the dynamics object, then create a copy of the integrator and system.

>>> from copy import deepcopy
>>> context = d.context()
>>> omm_system = context.getSystem()
>>> integrator = deepcopy(context.getIntegrator())

In order to perform umbrella sampling we will need to add a biasing potentials for the two dihedral angles. Here we will use simple harmonic biasing potentials:

First the Φ dihedral, which is formed by atom indices 4, 6, 8, and 14:

>>> import openmm
>>> import openmm.app
>>> bias_torsion_phi = openmm.CustomTorsionForce(
...     "0.5*k_phi*dtheta^2; dtheta = min(tmp, 2*pi-tmp); tmp = abs(theta - phi)"
... )
>>> bias_torsion_phi.addGlobalParameter("pi", math.pi)
>>> bias_torsion_phi.addGlobalParameter("k_phi", 100.0)
>>> bias_torsion_phi.addGlobalParameter("phi", 0.0)
>>> bias_torsion_phi.addTorsion(4, 6, 8, 14)

Next the Ψ dihedral, which is formed by atom indices 6, 8, 14, and 16:

>>> bias_torsion_psi = openmm.CustomTorsionForce(
...     "0.5*k_psi*dtheta^2; dtheta = min(tmp, 2*pi-tmp); tmp = abs(theta - psi)"
... )
>>> bias_torsion_psi.addGlobalParameter("pi", math.pi)
>>> bias_torsion_psi.addGlobalParameter("k_psi", 100.0)
>>> bias_torsion_psi.addGlobalParameter("psi", 0.0)
>>> bias_torsion_psi.addTorsion(6, 8, 14, 16)

We can now add these forces to the system:

>>> omm_system.addForce(bias_torsion_phi)
>>> omm_system.addForce(bias_torsion_psi)

In order to run the simulation we will create a new context using the system and integrator and set the initial positions.

>>> new_context = openmm.Context(omm_system, integrator, context.getPlatform())
>>> new_context.setPositions(context.getState(getPositions=True).getPositions())

Running the simulation

We are almost ready to run an umbrella sampling simulation. In this example we will sample the Φ and Ψ dihedral angles on a 36x36 grid. We will first set the biasing potential centers:

>>> m = 36
>>> M = m * m
>>> phi = np.linspace(-math.pi, math.pi, m, endpoint=False)
>>> psi = np.linspace(-math.pi, math.pi, m, endpoint=False)

During the simulation we will save trajectories to disk which can later be post-processed to compute the dihedral angles. We will create a directory in which to store the files:

>>> os.makedirs("./output/traj", exist_ok=True)

The sampling is performed by looping over each of the umbrella windows sequentially. For each window we set the biasing potential center and run an initial equilibration of 5000 steps. We then run a production simulation of 100 cycles of 100 steps each, saving trajectory after each cycle:

>>> for idx in range(M):
... phi_index = idx // m
... psi_index = idx % m
...
... # Set the center of the biasing potentials.
... new_context.setParameter("phi", phi[phi_index])
... new_context.setParameter("psi", psi[psi_index])
...
... # Initial equilibrium.
... integrator.step(5000)
...
... # Production sampling.
... file_handle = open(f"./output/traj/phi_{phi_index}_phi_{psi_index}.dcd", "bw")
... dcd_file = DCDFile(file_handle, prm.topology, dt=integrator.getStepSize())
... for x in range(100):
...     integrator.step(100)
...     state = new_context.getState(getPositions=True)
...     positions = state.getPositions()
...     dcd_file.writeModel(positions)
... file_handle.close()

Note

This is not a particulary efficient way to perform the sampling. In practice, since it’s possible to get good single core performance it is better to run the windows in parallel, either individually, or in blocks.

Analysing the results

The trajectories saved to disk can be post-processed to compute the dihedral angles, for example using the approach here. The free-energy surface can then be compute using MBAR, or UWHAM. Example code is provided in the FastMBAR tutorial.

The resulting free-energy surface should look similar to the one shown below: