Using ParmEd with MDAnalysis and OpenMM to simulate a selection of atoms

Here we use MDAnalysis to convert a ParmEd structure to an MDAnalysis Universe, select a subset of atoms, and convert it back to ParmEd to simulate with OpenMM.

Last updated: December 2022 with MDAnalysis 2.4.0-dev0

Last updated: December 2022

Minimum version of MDAnalysis: 1.0.0

Packages required:

• MDAnalysis [1, 2]

• MDAnalysisTests

• ParmEd

• OpenMM [3]

import parmed as pmd
import MDAnalysis as mda
from MDAnalysis.tests.datafiles import PRM7_ala2, RST7_ala2

import warnings
# suppress some MDAnalysis warnings when writing PDB files
warnings.filterwarnings('ignore')

Loading files: the difference between ParmEd and MDAnalysis

Both ParmEd and MDAnalysis read a number of file formats. However, while MDAnalysis is typically used to analyse simulations, ParmEd is often used to set them up. This requires ParmEd to read topology parameter information that MDAnalysis typically ignores, such as the equilibrium length and force constants of bonds in the system. For example, the ParmEd structure below.

pprm = pmd.load_file(PRM7_ala2, RST7_ala2)
pprm

<AmberParm 3026 atoms; 1003 residues; 3025 bonds; PBC (orthogonal); parameterized>

pprm.bonds[0]

<Bond <Atom C [10]; In ALA 0>--<Atom O [11]; In ALA 0>; type=<BondType; k=570.000, req=1.229>>

When MDAnalysis reads these files in, it does not include that information.

mprm = mda.Universe(PRM7_ala2, RST7_ala2, format='RESTRT')
mprm

<Universe with 3026 atoms>

The bond type simply shows the atom types involved in the connection.

mprm.atoms.bonds[0].type

('N3', 'H')

If you then convert this Universe to ParmEd, you can see that the resulting Structure is not parametrized.

mprm_converted = mprm.atoms.convert_to('PARMED')
mprm_converted

<Structure 3026 atoms; 1003 residues; 3025 bonds; parameterized>

While the bonds are present, there is no type information associated.

mprm_converted.bonds[0]

<Bond <Atom N [0]; In ALA 0>--<Atom H1 [1]; In ALA 0>; type=None>

Therefore, if you wish to use ParmEd functionality that requires parametrization on a MDAnalysis Universe, you need to create that Universe from a ParmEd structure in order to convert it back to something useable in ParmEd.

mprm_from_parmed = mda.Universe(pprm)
mprm_from_parmed

<Universe with 3026 atoms>

Now the bond type is actually a ParmEd Bond object.

mprm_from_parmed.bonds[0].type

<Bond <Atom N [0]; In ALA 0>--<Atom H1 [1]; In ALA 0>; type=<BondType; k=434.000, req=1.010>>

Using MDAnalysis to select atoms

One reason we might want to convert a ParmEd structure into MDAnalysis is to use its sophisticated atom selection syntax. While ParmEd has its own ways to select atoms, MDAnalysis allows you to select atoms based on geometric distance.

water = mprm_from_parmed.select_atoms('around 5 protein').residues.atoms
protein_shell = mprm_from_parmed.select_atoms('protein') + water
prm_protein_shell = protein_shell.convert_to('PARMED')

prm_protein_shell

<Structure 155 atoms; 46 residues; 154 bonds; PBC (orthogonal); parameterized>

Using ParmEd and OpenMM to create a simulation system

import sys
import openmm as mm
import openmm.app as app
from parmed import unit as u
from parmed.openmm import StateDataReporter, MdcrdReporter

You can create an OpenMM simulation system directly from a ParmEd structure, providing that it is parametrized.

system = prm_protein_shell.createSystem(nonbondedMethod=app.NoCutoff,
                                        constraints=app.HBonds,
                                        implicitSolvent=app.GBn2)

Here we set the integrator to do Langevin dynamics.

integrator = mm.LangevinIntegrator(
                        300*u.kelvin,       # Temperature of heat bath
                        1.0/u.picoseconds,  # Friction coefficient
                        2.0*u.femtoseconds, # Time step
                        )

We create the Simulation object and set particle positions.

sim = app.Simulation(prm_protein_shell.topology, system, integrator)
sim.context.setPositions(prm_protein_shell.positions)

We now minimise the energy.

sim.minimizeEnergy(maxIterations=500)

The reporter below reports energies and coordinates every 100 steps to stdout, but every 10 steps to the ala2_shell.nc file.

sim.reporters.append(
        StateDataReporter(sys.stdout, 100, step=True, potentialEnergy=True,
                          kineticEnergy=True, temperature=True, volume=True,
                          density=True)
)
sim.reporters.append(MdcrdReporter('ala2_shell.trj', 10, crds=True))

We can run dynamics for 500 steps (1 picosecond).

sim.step(500)

#"Step","Time (ps)","Potential Energy (kilocalorie/mole)","Kinetic Energy (kilocalorie/mole)","Total Energy (kilocalorie/mole)","Temperature (K)","Box Volume (angstrom**3)","Density (gram/(item*milliliter))"
100,0.20000000000000015,-623.6779995219885,20.140631869613383,-603.5373676523751,63.74314071570579,45325.8064191062,0.034909350700361955
200,0.4000000000000003,-614.1849904397706,38.40737137186695,-575.7776190679035,121.55559436896063,45325.8064191062,0.034909350700361955
300,0.6000000000000004,-606.5526783580306,48.61919832973248,-557.933480028298,153.87503334950912,45325.8064191062,0.034909350700361955
400,0.8000000000000006,-600.0374380078872,57.988937528818184,-542.0485004790689,183.52934648642113,45325.8064191062,0.034909350700361955
500,1.0000000000000007,-603.2854886173518,79.46589388029852,-523.8195947370533,251.50182419815255,45325.8064191062,0.034909350700361955

If we write a topology file out from our former protein_shell atomgroup, we can load the trajectory in for further analysis.

protein_shell.write('ala2_shell.pdb')

u = mda.Universe('ala2_shell.pdb', 'ala2_shell.trj')

u.trajectory

<TRJReader ala2_shell.trj with 50 frames of 155 atoms>

References

[1] R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler, D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein. MDAnalysis: A Python package for the rapid analysis of molecular dynamics simulations. In S. Benthall and S. Rostrup, editors, Proceedings of the 15th Python in Science Conference, pages 98-105, Austin, TX, 2016. SciPy, doi: 10.25080/majora-629e541a-00e.

[2] N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein. MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. J. Comput. Chem. 32 (2011), 2319-2327, doi:10.1002/jcc.21787. PMCID:PMC3144279

[3] Peter Eastman, Jason Swails, John D. Chodera, Robert T. McGibbon, Yutong Zhao, Kyle A. Beauchamp, Lee-Ping Wang, Andrew C. Simmonett, Matthew P. Harrigan, Chaya D. Stern, Rafal P. Wiewiora, Bernard R. Brooks, Vijay S. Pande. OpenMM 7: Rapid Development of High Performance Algorithms for Molecular Dynamics. PLoS Comput. Biol. 13:e1005659, 2017.

[4] Hai Nguyen, David A Case, Alexander S Rose. NGLview - Interactive molecular graphics for Jupyter notebooks. Bioinformatics. 34 (2018), 1241–1242, doi:10.1093/bioinformatics/btx789