Created
September 11, 2012 20:43
-
-
Save rmcgibbo/3701873 to your computer and use it in GitHub Desktop.
Milestoning With MSMBuilder/OpenMM :: Initial Sketch
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#!/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 | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment