In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
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
plot = ToyPlot()
pes = (
toys.OuterWalls([1.0, 1.0], [0.0, 0.0]) +
toys.Gaussian(-0.7, [12.0, 12.0], [-0.5, 0.5]) +
toys.Gaussian(-0.7, [12.0, 12.0], [-0.5, -0.5]) +
toys.Gaussian(-0.7, [12.0, 12.0], [0.5, -0.5])
)
plot.contour_range = np.arange(-1.5, 1.0, 0.1)
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.add_pes(pes)
plot.plot()
Out[3]:
In [4]:
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]
opX = paths.FunctionCV(name="opX", f=xval).with_diskcache()
opY = paths.FunctionCV(name="opY", f=yval).with_diskcache()
opXprime = paths.FunctionCV(name="opXprime", f=xprime).with_diskcache()
In [5]:
x_under_min = paths.CVDefinedVolume(opX, float("-inf"), -0.3)
x_over_max = paths.CVDefinedVolume(opX, 0.3, float("inf"))
y_under_min = paths.CVDefinedVolume(opY, float("-inf"), -0.3)
y_over_max = paths.CVDefinedVolume(opY, 0.3, 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 [6]:
#plot.add_states([stateA, stateB, stateC])
#plot.plot()
In [7]:
interfacesAB = paths.VolumeInterfaceSet(opX, float("-inf"), [-0.3, -0.2, -0.1])
interfacesAC = paths.VolumeInterfaceSet(opY, float("-inf"), [-0.3, -0.2, -0.1, 0.0])
interfacesBA = paths.VolumeInterfaceSet(opXprime, float("-inf"), [-0.3, -0.2, -0.1])
In [8]:
mistis = paths.MISTISNetwork(
[(stateA, interfacesAB, stateB),
(stateA, interfacesAC, stateC),
(stateB, interfacesBA, stateA)],
ms_outers=paths.MSOuterTISInterface.from_lambdas(
{iface: 0.0 for iface in [interfacesAB, interfacesBA]}
)
)
In [9]:
scheme = paths.DefaultScheme(mistis, engine=toy_eng)
In [10]:
tisAB = mistis.transitions[(stateA, stateB)]
tisAC = mistis.transitions[(stateA, stateC)]
tisBA = mistis.transitions[(stateB, stateA)]
In [11]:
import logging.config
logging.config.fileConfig("../resources/debug_logging.conf", disable_existing_loggers=False)
In [12]:
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=mistis.ms_outers
).run()
In [13]:
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 [14]:
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 [15]:
initial_trajectories = [s.trajectory for s in list(init_AB)+list(init_AC)+list(init_BA)]
In [16]:
plot.plot(initial_trajectories)
Out[16]:
In [17]:
sset = scheme.initial_conditions_from_trajectories(initial_trajectories)
print scheme.initial_conditions_report(sset)
In [18]:
plot.plot([s.trajectory for s in sset])
Out[18]:
In [19]:
minus_samples = []
for minus in mistis.minus_ensembles:
samp = minus.extend_sample_from_trajectories(
trajectories=sset,
replica=-mistis.minus_ensembles.index(minus)-1,
engine=toy_eng
)
minus_samples.append(samp)
sset = sset.apply_samples(minus_samples)
print scheme.initial_conditions_report(sset)
In [20]:
# 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 mistis.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 [21]:
#! skip
# tests should not run this, users should. Undoes the previous cell
sset = better_initial_conditions
In [22]:
plot.plot([s.trajectory for s in sset])
Out[22]:
In [23]:
#logging.config.fileConfig("../resources/debug_logging.conf", disable_existing_loggers=False)
storage = paths.Storage("mistis.nc", "w")
storage.save(template)
Out[23]:
In [24]:
mistis_calc = paths.PathSampling(
storage=storage,
move_scheme=scheme,
sample_set=sset
)
mistis_calc.save_frequency = 100
In [25]:
import logging.config
logging.config.fileConfig("../resources/logging.conf", disable_existing_loggers=False)
mistis_calc.run(100)
In [26]:
#! skip
# skip this during testing; leave for full calculation
mistis_calc.run_until(1000)
In [ ]: