Interface to OpenMM

This subpackage provides an interface to the OpenMM molecular simulation toolkit, expanding the structure capabilities of Biotite into the realm of molecular dynamics.

Like other biotite.interface subpackages, OpenMM objects are created with to_<x>() functions (to_topology(), to_system()) and converted back to an AtomArray with from_<x>() functions (from_topology(), from_system(), from_states(), etc.). Like in the RDKit interface tutorial, explaining all the functionalities of OpenMM would exceed the scope of this tutorial. Please refer to the OpenMM documentation for further information. This tutorial will only give a few examples of how OpenMM and Biotite work in tandem.

First example: MD simulation of lysozyme

Here the aim is to run the prominent ‘Lysozyme in water’ example in OpenMM. For structure preparation before and analysis afterwards, we will use Biotite.

import numpy as np
import openmm
import openmm.app as app
from openmm.unit import angstrom, kelvin, nanometer, picosecond
import biotite.database.rcsb as rcsb
import biotite.interface.openmm as openmm_interface
import biotite.interface.pymol as pymol_interface
import biotite.structure.io.pdbx as pdbx
import biotite.structure as struc

BUFFER = 10

atoms = pdbx.get_structure(
    pdbx.BinaryCIFFile.read(rcsb.fetch("1aki", "bcif")), model=1, include_bonds=True
)
# Remove solvent, as proper solvent addition is handled by OpenMM
molecule = atoms[struc.filter_amino_acids(atoms)]
# Create a box with some buffer around the protein
# This box will be used as starting simulation box
box_dim = np.max(molecule.coord, axis=0) - np.min(molecule.coord, axis=0)
molecule.box = np.diag(box_dim + BUFFER)

The AtomArray is converted to an openmm.topology.Topology simply by calling to_topology().

# Create an OpenMM Topology from the AtomArray
topology = openmm_interface.to_topology(molecule)

The lysozyme structure we use does not have hydrogen atoms, but we require them for the MD simulation. Hence, we use functionality from OpenMM to model the missing hydrogen atoms.

force_field = openmm.app.ForceField("amber14-all.xml", "amber14/tip3pfb.xml")
# Add hydrogen atoms and water
modeller = app.Modeller(topology, molecule.coord * angstrom)
modeller.addHydrogens(force_field)
modeller.addSolvent(force_field)
topology = modeller.topology

The last prerequisite for the MD simulation is the energy minimization to get a proper starting conformation. After that, we can finally run the simulation for the given number of steps. After a certain number of steps, we record the current state (i.e. conformation) of the system for later analysis.

TIME_STEP = 0.004
FRAME_STEP = 2.0
N_FRAMES = 60

system = force_field.createSystem(
    topology,
    nonbondedMethod=app.PME,
    nonbondedCutoff=1 * nanometer,
    constraints=app.HBonds,
)
integrator = openmm.LangevinMiddleIntegrator(
    300 * kelvin, 1 / picosecond, TIME_STEP * picosecond
)
simulation = app.Simulation(topology, system, integrator)
simulation.context.setPositions(modeller.positions)
simulation.minimizeEnergy()

# Run simulation and record the current state (the coordinates)
# every FRAME_STEP picoseconds
states = []
for i in range(N_FRAMES):
    simulation.step(FRAME_STEP // TIME_STEP)
    states.append(simulation.context.getState(getPositions=True))

While each state contains information about the current atom coordinates, box dimensions etc., it misses the topology of the system, i.e. residues, atoms and bonds. Hence, to get the trajectory as AtomArrayStack via from_states(), we need to provide an AtomArray as template. However, we cannot use the original input lysozyme structure from above, as hydrogen atoms were added meanwhile and thus the topology has changed. To solve this problem we obtain an updated template first using from_topology().

template = openmm_interface.from_topology(topology)
trajectory = openmm_interface.from_states(template, states)
# Center protein in box
trajectory.coord -= struc.centroid(
    trajectory.coord[:, struc.filter_amino_acids(trajectory)]
)[:, np.newaxis, :]
trajectory.coord += np.sum(trajectory.box / 2, axis=-2)[:, np.newaxis, :]
# Remove segmentation over periodic boundary
trajectory = struc.remove_pbc(trajectory)

Now we could analyze the trajectory using functionality from biotite.structure, as exemplified in the trajectory tutorial, or use some custom analysis. For this tutorial we will simply visualize the trajectory using biotite.interface.pymol.

pymol_object = pymol_interface.PyMOLObject.from_structure(trajectory)
pymol_object.color("biotite_lightorange", struc.filter_amino_acids(trajectory))
# Draw simulation box
# As the box does not extend much during the simulation, simply draw the mean size
box = pymol_interface.draw_box(
    np.mean(trajectory.box, axis=0), color=(0, 0, 0), width=1
)

# Create an isometric view
pymol_interface.cmd.turn("y", 45)
pymol_interface.cmd.turn("x", 45)
pymol_object.zoom(buffer=15)
pymol_interface.cmd.set("orthoscopic", 1)

pymol_interface.cmd.mset()
pymol_interface.play(
    (1400, 1400), format="mp4", html_attributes="loop autoplay width=700"
)