Multiple Interface Set TIS on a Toy Model
This example demostrates multiple interface set TIS on a two dimensional toy model. It consists of thre notebooks: one to setup and run the sampling, one to calculate the flux, and one to analyze the results. In MISTIS, the flux must be performed as a separate calculation.
In [1]:
from __future__ import print_function
%matplotlib inline
import openpathsampling as paths
import numpy as np
# we use the %run magic because this isn't in a package
%run ../resources/toy_plot_helpers.py
In [2]:
import openpathsampling.engines.toy as toys
pes = (
toys.OuterWalls([1.0, 1.0], [0.0, 0.0]) +
toys.Gaussian(-1.0, [12.0, 12.0], [-0.5, 0.5]) +
toys.Gaussian(-1.0, [12.0, 12.0], [-0.5, -0.5]) +
toys.Gaussian(-1.0, [12.0, 12.0], [0.5, -0.5])
)
topology=toys.Topology(
n_spatial = 2,
masses =[1.0, 1.0],
pes = pes
)
integ = toys.LangevinBAOABIntegrator(dt=0.02, temperature=0.1, gamma=2.5)
options={
'integ' : integ,
'n_frames_max' : 5000,
'n_steps_per_frame' : 1
}
toy_eng = toys.Engine(
options=options,
topology=topology
)
toy_eng.initialized = True
template = toys.Snapshot(
coordinates=np.array([[-0.5, -0.5]]),
velocities=np.array([[0.0,0.0]]),
engine=toy_eng
)
toy_eng.current_snapshot = template
paths.PathMover.engine = toy_eng
In [3]:
plot = ToyPlot()
plot.contour_range = np.arange(-1.5, 1.0, 0.1)
plot.add_pes(pes)
plot.plot()
Out[3]:
In [ ]:
def xval(snapshot):
return snapshot.xyz[0][0]
def xprime(snapshot):
# this only exists until we set up the ability for the order parameter to decrease
return -snapshot.xyz[0][0]
def yval(snapshot):
return snapshot.xyz[0][1]
cvX = paths.FunctionCV(name="cvX", f=xval).with_diskcache()
cvY = paths.FunctionCV(name="cvY", f=yval).with_diskcache()
cvXprime = paths.FunctionCV(name="cvXprime", f=xprime).with_diskcache()
In [ ]:
x_under_min = paths.CVDefinedVolume(cvX, float("-inf"), -0.35)
x_over_max = paths.CVDefinedVolume(cvX, 0.35, float("inf"))
y_under_min = paths.CVDefinedVolume(cvY, float("-inf"), -0.35)
y_over_max = paths.CVDefinedVolume(cvY, 0.35, float("inf"))
stateA = (x_under_min & y_under_min).named("A")
stateB = (x_over_max & y_under_min).named("B")
stateC = (x_under_min & y_over_max).named("C")
In [ ]:
#plot.add_states([stateA, stateB, stateC])
#plot.plot()
In [ ]:
interfacesAB = paths.VolumeInterfaceSet(cvX, float("-inf"), [-0.35, -0.3, -0.27, -0.24, -0.2, -0.1])
interfacesAC = paths.VolumeInterfaceSet(cvY, float("-inf"), [-0.35, -0.3, -0.27, -0.24, -0.2, -0.1, 0.0])
interfacesBA = paths.VolumeInterfaceSet(cvXprime, float("-inf"), [-0.35, -0.3, -0.27, -0.24, -0.2, -0.1])
In [ ]:
ms_outer = paths.MSOuterTISInterface.from_lambdas(
{iface: 0.0 for iface in [interfacesAB, interfacesBA]}
)
network = paths.MISTISNetwork(
[(stateA, interfacesAB, stateB),
(stateA, interfacesAC, stateC),
(stateB, interfacesBA, stateA)],
ms_outers=ms_outer,
strict_sampling=True
)
In [ ]:
scheme = paths.DefaultScheme(network, engine=toy_eng)
In [ ]:
tisAB = network.transitions[(stateA, stateB)]
tisAC = network.transitions[(stateA, stateC)]
tisBA = network.transitions[(stateB, stateA)]
In [ ]:
import logging.config
logging.config.fileConfig("../resources/debug_logging.conf", disable_existing_loggers=False)
In [ ]:
snapA = toys.Snapshot(
coordinates=np.array([[-0.5, -0.5]]),
velocities=np.array([[0.5, 0.0]])
)
init_AB = paths.FullBootstrapping(
transition=tisAB,
snapshot=snapA,
engine=toy_eng,
forbidden_states=[stateC],
extra_ensembles=network.ms_outers
).run()
In [ ]:
snapA = toys.Snapshot(
coordinates=np.array([[-0.5, -0.5]]),
velocities=np.array([[0.0, 0.5]])
)
init_AC = paths.FullBootstrapping(
transition=tisAC,
snapshot=snapA,
engine=toy_eng,
forbidden_states=[stateB]
).run()
In [ ]:
snapB = toys.Snapshot(
coordinates=np.array([[0.5, -0.5]]),
velocities=np.array([[-0.5, 0.0]])
)
init_BA = paths.FullBootstrapping(
transition=tisBA,
snapshot=snapB,
engine=toy_eng,
forbidden_states=[stateC]
).run()
In [ ]:
initial_trajectories = [s.trajectory for s in list(init_AB)+list(init_AC)+list(init_BA)]
In [ ]:
plot.plot(initial_trajectories)
In [ ]:
sset = scheme.initial_conditions_from_trajectories(initial_trajectories)
print(scheme.initial_conditions_report(sset))
In [ ]:
plot.plot([s.trajectory for s in sset])
In [ ]:
minus_samples = []
for minus in network.minus_ensembles:
samp = minus.extend_sample_from_trajectories(
trajectories=sset,
replica=-network.minus_ensembles.index(minus)-1,
engine=toy_eng
)
minus_samples.append(samp)
print(minus_samples)
print(network.minus_ensembles)
sset = sset.apply_samples(minus_samples)
print(scheme.initial_conditions_report(sset))
In [ ]:
# The next two cells are for tests.
# This cell creates initial conditions that will pass analysis on low data.
# The next cell undoes that to use a better initial condition in practice.
better_initial_conditions = sset
for transition in network.sampling_transitions:
outermost_traj = sset[transition.ensembles[-1]].trajectory
for ensemble in transition.ensembles:
original = sset[ensemble]
sample = paths.Sample(replica=original.replica,
ensemble=ensemble,
trajectory=outermost_traj)
sset = sset.apply_samples(sample)
In [ ]:
# NBVAL_SKIP
# tests should not run this, users should. Undoes the previous cell
sset = better_initial_conditions
In [ ]:
plot.plot([s.trajectory for s in sset])
In [ ]:
#logging.config.fileConfig("../resources/debug_logging.conf", disable_existing_loggers=False)
storage = paths.Storage("mistis.nc", "w")
storage.save(template)
In [ ]:
mistis_calc = paths.PathSampling(
storage=storage,
move_scheme=scheme,
sample_set=sset
)
mistis_calc.save_frequency = 100
In [ ]:
#import logging.config
#logging.config.fileConfig("../resources/logging.conf", disable_existing_loggers=False)
mistis_calc.run(100)
In [ ]:
# NBVAL_SKIP
# skip this during testing; leave for full calculation
mistis_calc.run_until(100000)
In [ ]:
(toy_mistis_1_setup_run.ipynb; toy_mistis_1_setup_run.py)
In [ ]:
%matplotlib inline
import openpathsampling as paths
import numpy as np
As always, we load things from files so we don't have to set them up again.
In [ ]:
old = paths.Storage("mistis.nc", 'r')
engine = old.engines[0]
network = old.networks[0]
states = set(network.initial_states + network.final_states)
In [ ]:
# must ensure that the diskcache is disabled in order to save,
# otherwise it looks for things that aren't there!
cvs = old.cvs[:]
for cv in cvs:
cv.disable_diskcache()
The flux_pairs variable is a list of 2-tuples, where the first element is the state we're calculating the flux out of, and the second element is the interface we're calculating the flux through.
In [ ]:
flux_pairs = [(t.stateA, t.interfaces[0]) for t in network.transitions.values()]
Set up the simulation and run it!
In [ ]:
sim = paths.DirectSimulation(
storage=None,
engine=engine,
states=states,
flux_pairs=flux_pairs,
initial_snapshot=old.snapshots[0]
)
In [ ]:
%%time
sim.run(150000) # 30 sec
#sim.run(1500000) # 6 min
#sim.run(15000000) # 60 min
#sim.run(150000000) # 10 hr
#sim.run(800000000) # 2 days
Now we move on to the analysis.
In [ ]:
sim.rate_matrix
In [ ]:
sim.n_transitions
In [ ]:
fluxes = sim.fluxes
for f in fluxes:
print f[0].name, f[1].name, fluxes[f]
In [ ]:
sim.n_flux_events
In [ ]:
sim.results
In [ ]:
output = paths.Storage("direct_simulation.nc", 'w')
output.save(old.snapshots[0])
output.save(sim)
output.tag['direct_results'] = sim.results
output.close()
(toy_mistis_2_flux.ipynb; toy_mistis_2_flux.py)
In [1]:
from __future__ import print_function
# if our large test file is available, use it. Otherwise, use file generated from toy_mistis_1_setup_run.ipynb
import os
test_file = "../toy_mistis_1k_OPS1.nc"
filename = test_file if os.path.isfile(test_file) else "mistis.nc"
In [2]:
print(filename)
In [3]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import openpathsampling as paths
In [4]:
%%time
storage = paths.AnalysisStorage(filename)
In [5]:
mistis = storage.networks.load(0)
In [6]:
mistis.hist_args['max_lambda'] = { 'bin_width' : 0.01, 'bin_range' : (-0.35, 0.5) }
mistis.hist_args['pathlength'] = { 'bin_width' : 5, 'bin_range' : (0, 150) }
In [7]:
%%time
scheme = storage.schemes[0]
scheme.move_summary(storage.steps)
In [8]:
scheme.move_summary(storage.steps, 'shooting')
In [9]:
scheme.move_summary(storage.steps, 'minus')
In [10]:
scheme.move_summary(storage.steps, 'repex')
In [11]:
# we need to load the states and the innermost interface for each transition
stateA = storage.volumes['A']
stateB = storage.volumes['B']
stateC = storage.volumes['C']
inner_AB = mistis.transitions[(stateA, stateB)].interfaces[0]
inner_AC = mistis.transitions[(stateA, stateC)].interfaces[0]
inner_BA = mistis.transitions[(stateB, stateA)].interfaces[0]
In [12]:
# got these from mistis_flux.ipynb
fluxes = {(stateA, inner_AB): 0.0916199741819,
(stateA, inner_AC): 0.0915271110694,
(stateB, inner_BA): 0.0916882528979}
mistis.set_fluxes(fluxes)
In [13]:
%%time
rate = mistis.rate_matrix(storage.steps, force=True)
rate
In [35]:
import pandas as pd
pd.options.display.float_format = '{:.3e}'.format
rate
Out[35]:
In [46]:
# this can be copy-pasted into an article
print(rate.to_latex(float_format='{:.3e}'.format))
In [14]:
trans = list(mistis.transitions.values())[2]
trans_hists = trans.histograms['max_lambda']
print(trans)
In [15]:
for hist in trans_hists:
cross_prob = trans_hists[hist].reverse_cumulative()
plt.plot(cross_prob.x, np.log(cross_prob))
plt.plot(trans.tcp.x, np.log(trans.tcp), '-k', lw=2)
Out[15]:
In [16]:
len(storage.steps)
Out[16]:
In [17]:
#import logging.config
#logging.config.fileConfig("../resources/debug_logging.conf", disable_existing_loggers=False)
In [18]:
n_blocks = 1 # for testing code
In [19]:
# NBVAL_SKIP
n_blocks = 5 # for real examples
In [20]:
resampling = paths.numerics.BlockResampling(storage.steps, n_blocks=n_blocks)
In [21]:
rate_df_func = lambda steps: mistis.rate_matrix(steps, force=True)
In [22]:
%%time
rates = paths.numerics.ResamplingStatistics(function=rate_df_func, inputs=resampling.blocks)
In [37]:
rates.mean
Out[37]:
In [38]:
rates.std
Out[38]:
In [39]:
rates.percentile(0)
Out[39]:
In [40]:
rates.percentile(25)
Out[40]:
In [41]:
rates.percentile(50)
Out[41]:
In [42]:
rates.percentile(75)
Out[42]:
In [43]:
rates.percentile(100)
Out[43]:
In [ ]: