Milestoning With MSMBuilder/OpenMM :: Initial Sketch

milestone.py

#!/usr/bin/env python
from simtk.openmm.app import *
from simtk.openmm import *
from simtk.unit import *

from msmbuilder import Trajectory
import msmbuilder.metrics

#interface to the msmbuilder fast RMSD calculator
rmsd_calculator = msmbuilder.metrics.RMSD()

# load in the path file with MSMBuilder's DCD reader
# note that this DCD/PDB should only have the reduced coordinates
path = Trajectory.LoadFromDCD('PATH_FILE.dcd', PDBFilename='PATH_FILE.pdb')
test_conf = path[0] # this is going to be overwritten later, but it's just going
                    # to be a place to dump the coordinates that we get out from
                    # OpenMM

#this precalculates some quantities needed for the Theobald RMSD calc.
rmsd_prepared_path = rmsd_calculator.prepare_trajectory(path)

# start up OpenMM to run the simulation
n_steps = 1000

pdb = PDBFile('input.pdb')

# these should be indices of the C-alpha atoms (the ones in PATH_FILE.pdb), as they
# are numbered in input.pdb
c_alpha_indices = [1,6,10,20] # only an example


forcefield = ForceField('amber99sb.xml', 'tip3p.xml')
system = forcefield.createSystem(pdb.topology, nonbondedMethod=PME, nonbondedCutoff=1*nanometer, constraints=HBonds)
integrator = LangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
simulation = Simulation(pdb.topology, system, integrator)
simulation.context.setPositions(pdb.positions)
simulation.minimizeEnergy()
simulation.reporters.append(PDBReporter('output.pdb', n_steps))
simulation.reporters.append(StateDataReporter(stdout, n_steps, step=True, potentialEnergy=True, temperature=True))

total_steps_taken = 0
while True:
    simulation.step(n_steps)
    total_steps_taken += n_steps

    # extract the coordinates, and prepare them for the RMSD calculation
    coordinates = simulation.context.getPositions(as_numpy=True)._value

    # extract only the c_alpha_indices,
    reduced_coordinates = coordinates[c_alpha_indices, :]
    
    test_conf['XYZList'] = np.array([reduced_coordinates])
    
    #this precalculates some quantities needed for the Theobald RMSD calc.
    rmsd_prepared_conf = rmsd_calculator.prepare_trajectory(test_conf)
    
    #calculate the rmsd from our current conf to each of the confs in the "path"
    distances = rmsd_calculator.one_to_all(rmsd_prepared_conf, rmsd_prepared_path, 0)
    
    # check if test_conf is closer to any of the other PATH conformations than it
    # is to the very first (index 0) path conf, which would indicate barrier crossing
    if np.any(distances[1:] < distances[0]):
        # break out of the loop
        print 'Barrier crossing detected'
        print 'total_steps_taken so far = %s' % total_steps_taken
        break