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)
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)
In [ ]: