Seeding OpenMM Simulations with PyRosetta and ParmEd

In this example, we’ll go over starting a small protein (dodeca-alanine) simulation from scratch, using PyRosetta and OpenMM.

You can also follow along with the script called simulate_ala12_mutant.py in examples/rosetta/.

Generate pose from primary structure

To start, let’s create a PyRosetta Pose directly from a sequence of twelve alanines:

from rosetta import init, pose_from_sequence

init()

seq = 12*'A'
pose = pose_from_sequence(seq)

Mutate a residue

Want to mutate the fifth residue to a lysine? PyRosetta makes it as easy as:

from toolbox import mutate_residue

mutant = mutate_residue(pose, 5, 'K')

Load Pose into ParmEd

The next step is to use ParmEd’s load_rosetta() function to convert our mutant into a Structure:

from parmed import load_rosetta

struct = load_rosetta(mutant)

Start an OpenMM simulation

Finally, we can use the following code to solvate and start simulating our ParmEd Structure:

from simtk.openmm import *
from simtk.openmm.app import *
from simtk.unit import *

ff = ForceField('amber99sbildn.xml', 'tip3p.xml')

# Solvate Structure
mod = Modeller(struct.topology, struct.positions)
mod.addSolvent(ff, model='tip3p', boxSize=Vec3(4, 4, 4)*nanometers,
               positiveIon='Na+', negativeIon='Cl-',
               ionicStrength=0.1*molar)

# Create OpenMM System
system = ff.createSystem(mod.topology, nonbondedMethod=PME,
                         nonbondedCutoff=1*nanometers,
                         constraints=HBonds)
integrator = LangevinIntegrator(300*kelvin, 1/picoseconds,
                                2*femtoseconds)

system.addForce(MonteCarloBarostat(1*bar, 300*kelvin))

simulation = Simulation(mod.topology, system, integrator)
simulation.context.setPositions(mod.positions)

# Minimize System
simulation.minimizeEnergy(maxIterations=1000)

# Equilibrate System
simulation.reporters.append(
    PDBReporter('dodecaalanine.solv.pdb', 50000))

simulation.context.setVelocitiesToTemperature(300)

simulation.step(50000)