Committor Simulation


In [1]:
import openpathsampling as paths
import mdtraj as md
import numpy as np
from simtk import unit

In [2]:
input_storage = paths.Storage("snapshots.nc", "r")

In [3]:
print input_storage.file_size_str
print len(input_storage.snapshots)


42.89MB
2000

In [4]:
# 2000 instead of 1000 because reversed snapshots are automatically counted
snapshots_to_run = input_storage.snapshots[::2]

In [5]:
# conveniently, we saved the engine; have to re-do everything else
engine = input_storage.engines[0]

# set up the collective variables for our states 
phi = paths.MDTrajFunctionCV(name="phi",
                             f=md.compute_dihedrals,
                             topology=engine.topology,
                             indices=[[4, 6, 8, 14]])
psi = paths.MDTrajFunctionCV(name="psi",
                             f=md.compute_dihedrals,
                             topology=engine.topology,
                             indices=[[6, 8, 14, 16]])

In [6]:
# define our states
deg = 180.0/np.pi
C_7eq = (paths.PeriodicCVDefinedVolume(phi, lambda_min=-180/deg, lambda_max=0/deg, 
                                     period_min=-np.pi, period_max=np.pi) &
         paths.PeriodicCVDefinedVolume(psi, lambda_min=100/deg, lambda_max=200/deg,
                                     period_min=-np.pi, period_max=np.pi)
        ).named("C_7eq")
# similarly, without bothering with the labels:
alpha_R = (paths.PeriodicCVDefinedVolume(phi, -180/deg, 0/deg, -np.pi, np.pi) &
           paths.PeriodicCVDefinedVolume(psi, -100/deg, 0/deg, -np.pi, np.pi)).named("alpha_R")

In [7]:
# OpenMM requires everything to have units
# beta = 1.0 / (k_B T)
temperature = 300.0 * unit.kelvin
beta = 1.0 / (temperature * unit.BOLTZMANN_CONSTANT_kB)

In [8]:
randomizer = paths.RandomVelocities(beta=beta, engine=engine)

In [9]:
output_storage = paths.Storage("committor_simulation.nc", "w")

In [10]:
simulation = paths.CommittorSimulation(storage=output_storage,
                                       engine=engine,
                                       states=[C_7eq, alpha_R],
                                       randomizer=randomizer,
                                       initial_snapshots=snapshots_to_run)

In [11]:
simulation.run(n_per_snapshot=10)


Working on snapshot 1000 / 1000; shot 10 / 10

In [ ]: