diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 834de8b..801cfd7 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -15,9 +15,9 @@ defaults: shell: bash -l {0} jobs: - tests: + intro-tutorial: runs-on: ubuntu-latest - name: "Test answer key" + name: "Intro tutorial" strategy: matrix: CONDA_PY: @@ -34,6 +34,7 @@ jobs: python-version: ${{ matrix.CONDA_PY }} environment-file: binder/environment.yml activate-environment: ops-tutorial + miniforge-variant: Mambaforge - name: "Install testing tools" run: python -m pip install pytest nbval - name: "Conda info" @@ -51,3 +52,42 @@ jobs: 2_tps_analysis_tutorial.ipynb \ 3_committor_analysis_tutorial.ipynb \ 4_mstis_sampling_tutorial.ipynb + + advanced-tutorial: + runs-on: ubuntu-latest + name: "Advanced tutorial" + strategy: + matrix: + CONDA_PY: + - 3.7 + - 3.8 + - 3.9 + + steps: + - uses: actions/checkout@v2 + - uses: actions/setup-python@v2 + - uses: conda-incubator/setup-miniconda@v2 + with: + auto-update-conda: true + python-version: ${{ matrix.CONDA_PY }} + environment-file: binder/environment.yml + activate-environment: ops-tutorial + miniforge-variant: Mambaforge + - name: "Install testing tools" + run: python -m pip install pytest nbval + - name: "Conda info" + run: | + conda info + conda list + - name: "Patch answers" + run: source devtools/patch-all + - name: "Run tests" + run: | + pytest --nbval-lax 5_custom_shooting_setup.ipynb + source devtools/run5.sh + pytest --nbval-lax \ + 6_custom_shooting_analysis.ipynb \ + 7_parallel_tis_setup.ipynb + source devtools/run7.sh + pytest --nbval-lax 8_parallel_tis_analysis.ipynb + diff --git a/1_tps_sampling_tutorial.ipynb b/1_tps_sampling_tutorial.ipynb index 44b5bd6..bf34b4e 100644 --- a/1_tps_sampling_tutorial.ipynb +++ b/1_tps_sampling_tutorial.ipynb @@ -326,7 +326,9 @@ "source": [ "## Running the simulation\n", "\n", - "Now all the parts are in place, and we can run the simulation! First, we open a file to store the results in, and then we build the simulation and run it." + "Now all the parts are in place, and we can run the simulation! First, we open a file to store everything in, and then we build the simulation and run it.\n", + "\n", + "In practice, for a long-running simulation you might save all the relevant information to a \"setup\" file, and then use the OPS Command Line Interface (CLI) to run it (perhaps on a remote computer). The advantage of that approach is that you can re-use the same exact objects, guaranteeing that they will behave the same way and making analysis easier. That approach is described in notebook 5 of this tutorial. In this notebook, we'll directly run the simulation." ] }, { @@ -398,7 +400,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.7.3" + "version": "3.7.10" }, "toc": { "base_numbering": 1, diff --git a/5_advanced_customize_shooting.ipynb b/5_advanced_customize_shooting.ipynb deleted file mode 100644 index 5e2819a..0000000 --- a/5_advanced_customize_shooting.ipynb +++ /dev/null @@ -1,1107 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Customizing shooting moves\n", - "\n", - "The tutorial focuses on how to customize shooting moves in the move scheme. We discuss three major reasons to do that:\n", - "\n", - "1. To use a different kind of shooting move, such as two-way shooting, when needed.\n", - "2. To use a different shooting point selector, such as a Gaussian biased selector, to get better efficiency.\n", - "3. To create move schemes that only sample part of the network, for example, to perform TIS with each ensemble sampled in parallel.\n", - "\n", - "In addition, you'll learn about:\n", - "\n", - "* Creating a \"setup\" file with engine and other simulation information.\n", - "* Using the OpenPathSampling command line interface (CLI) to run simulations from the setup file.\n", - "\n", - "Note that you'll need to install the CLI, which can be done with either `pip install openpathsampling-cli` or `conda install -c conda-forge openpathsampling-cli`.\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import matplotlib.pyplot as plt\n", - "\n", - "import openpathsampling as paths\n", - "from openpathsampling import strategies" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## A simple toy model example system\n", - "\n", - "In this notebook, we'll use a simple toy system. For the most part, the ideas here directly generalize to other engines, although we'll make a few comments where the units associated with OpenMM require special care.\n", - "\n", - "We'll start by loading a number things from the storage file:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "storage = paths.Storage(\"./inputs/2_state_toy.nc\", mode='r')\n", - "state_A = storage.volumes['A']\n", - "state_B = storage.volumes['B']\n", - "cv = storage.cvs['x']\n", - "engine = storage.engines['toy_engine']\n", - "initial_conditions = storage.tags['initial_conditions']" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The toy model here uses a simple 2D potential energy surface, described by:\n", - "\n", - "$$V(x, y) = x^6 + y^6 - e^{-12(x+0.6)^2 - 5 y^2} - e^{-12(x-0.6)^2 - 5 y^2}$$\n", - "\n", - "We can visualize the toy engine using the `toy_plot_helpers` included in this repository:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "%run toy_plot_helpers.py\n", - "\n", - "pes = engine.topology.pes\n", - "\n", - "plot = ToyPlot()\n", - "plot.contour_range = np.arange(-1.5, 1.0, 0.1)\n", - "plot.add_pes(pes)\n", - "fig = plot.plot() " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We run this simulation at a temperature of $T=0.1$, so the barrier between those wells is about $10\\ k_\\text{B} T$." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Changing the shooting move type (e.g., two-way shooting)\n", - "\n", - "TPS is usually introduced with two-way shooting. In that algorithm, you select a frame of the trajectory, make some modification (typically changing the velocities in some way consistent with the thermodynamic ensemble), and then integrate the equations of motion forward and backward.\n", - "\n", - "However, this is problematic if the velocity memory is short compared to your TPS trajectory length, as is often the case in condensed phase systems. In this case, each shot might as well be a committor shot. The ideal place to shoot from is the where the committor is 1/2, where each shot has only a 50% chance of landing in the correct state. But this gives you a maximum of a 25% acceptance rate.\n", - "\n", - "To get around this, we typically use the one-way shooting algorithm when working with condensed matter systems. One-way shooting works under the assumption that you are using a stochastic integrator. If that's the case, then instead of changing the velocities at the shooting point, you can use the fact that the integrator will generate a new sequence of random numbers. Therefore, you create a new trajectory by only running in one direction, and that trajectory is still a valid, physical trajectory.\n", - "\n", - "In general, OpenPathSampling defaults to using one-way shooting. However, there are circumstances where you might want to use the older two-way shooting algorithm. For example, if the velocity memory is shorter than the period between saved frames, it might be better to use two-way shooting. If you intend to use deterministic dynamics, then one-way shooting is not possible.\n", - "\n", - "This part of the tutorial will show you how to use two-way shooting instead of one-way shooting in OPS. There are two ways of doing this: either replace the existing one-way shooting strategy, or create a new move scheme from scratch." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tps_network = paths.TPSNetwork(initial_states=state_A, final_states=state_B)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Replacing parts of an existing move scheme\n", - "\n", - "Frequently, the easiest way to modify a move scheme (especially a complex move scheme) is to replace the parts that you want to change. In this case, we can start with the move scheme that only includes 1-way shooting, and work from there." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "modified_scheme = paths.OneWayShootingMoveScheme(\n", - " network=tps_network,\n", - " engine=engine\n", - ").named(\"2_way_from_1_way\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "For our two-way shooting, we'll need to create a `TwoWayShootingStrategy`. This will require a modifier; to keep things simple, we'll completely randomize velocities (consistent with a given temperature) using `paths.RandomVelocities`. In order to be consistent between engines, `RandomVelocities` takes its input as the inverse temperature, $\\beta = 1/(k_\\text{B}T)$.\n", - "\n", - "For the toy engine, we can obtain the temperature from `engine.integ.temperature`, and we work in units where $k_\\text{B}=1$, so it is easy to calculate $\\beta$. For other engines, you'll need to use the correct value of $k_\\text{B}$. For OpenMM, you also need to worry about units: use `simtk.unit.BOLTZMANN_CONSTANT_kB`.\n", - "\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "print(engine.integ.temperature)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# YOUR TURN: Set beta correctly (yes, this is as easy as you think).\n", - "# beta = ... # fill in the ellipsis and uncomment this line" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "modifier = paths.RandomVelocities(beta=beta, engine=engine)\n", - "shooting_strategy = strategies.TwoWayShootingStrategy(modifier=modifier, engine=engine)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "To modify a move scheme, just `append` new strategies. In this case, the new strategy for the group of moves called `'shooting'` will be overwritten. If you gave the `TwoWayShootingStrategy` a different string for its `group` argument, then you would create a second group of movers, and each move would have a 50/50 chance of using one-way shooting or of using two-way shooting." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "modified_scheme.append(shooting_strategy)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Building a move scheme from scratch\n", - "\n", - "In this particular case, our overall move scheme will not be too complicated. In fact, let's start by looking at the at the code for " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import inspect\n", - "from IPython.display import Code\n", - "\n", - "Code(inspect.getsource(paths.OneWayShootingMoveScheme))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The output above is the full source code for the `OneWayShootingMoveScheme`. You can see that all it does is pass the network up to the superclass's initialization, and then `append` two move strategies to itself. So we can recreate this:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "two_way_scheme = paths.MoveScheme(network=tps_network).named(\"2_way\")\n", - "global_strategy = strategies.OrganizeByMoveGroupStrategy()\n", - "two_way_scheme.append([shooting_strategy, global_strategy])" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Changing the shooting point selection\n", - "\n", - "A common problem in one-way shooting is that paths do not decorrelate quickly enough. This typically happens when there is a significant barrier withing the transition region, such that shooting points before the barrier always go back to the initial state (so only backward shots are accepted) and shooting points after the barrier always to do the final state (so only forward shots are accepted).\n", - "\n", - "The same issue manifests in two-way shooting as a significant decrease in the acceptance rate for the shooting move. One side of the barrier always creates $A\\to A$ trials, while the other side always creates $B\\to B$ trials. With fewer $A\\to B$ trials, the acceptance rate drops." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import numpy as np\n", - "xvals = np.arange(-1.0, 1.0, 0.01)\n", - "\n", - "def plot_gaussian_bias(selector):\n", - " alpha = selector.alpha\n", - " x_0 = selector.l_0\n", - " gaussians = np.exp(-alpha*(xvals - x_0)**2)\n", - " plt.plot(xvals, gaussians)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "biased_shooting_scheme = paths.MoveScheme(network=tps_network).named(\"biased_shooting\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "gaussian_sel = paths.GaussianBiasSelector(cv, alpha=100.0, l_0=0.0)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plot_gaussian_bias(gaussian_sel)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "paths.GaussianBiasSelector?" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# YOUR TURN: Create a move scheme using a two-way shooting strategy with this selector.\n", - "# 1. Create a two-way shooting strategy that uses this Gaussian selector. (Use the \n", - "# selector keyword argument of TwoWayShootingStrategy.)\n", - "# 2. Append that strategy and the global_strategy to the biased_shooting_scheme\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Next we verify that the initial conditions will work for the move schemes we've created." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "_ = two_way_scheme.initial_conditions_from_trajectories(initial_conditions)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "_ = biased_shooting_scheme.initial_conditions_from_trajectories(initial_conditions)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Run the simulations\n", - "\n", - "Now that we've created these various move schemes, we'll save them to a file, and use that as the input file for the OpenPathSampling command line interface." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "shooting_setup = paths.Storage(\"shooting_setup.nc\", mode='w')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# saving everything will take a few minutes\n", - "shooting_setup.tags['initial_conditions'] = initial_conditions\n", - "shooting_setup.save(two_way_scheme)\n", - "shooting_setup.save(biased_shooting_scheme)\n", - "shooting_setup.close()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Now, use the OpenPathSampling command line interface to run simulations with these. Run at least 500 steps (`-n 500`), or run more if you'd like.\n", - "\n", - "```\n", - "$ openpathsampling pathsampling shooting_setup.nc -o 2_way.nc --scheme 2_way -n 500\n", - "\n", - "$ openpathsampling pathsampling shooting_setup.nc -o biased.nc --scheme biased_shooting -n 500\n", - "```\n", - "\n", - "Those simulations will take a few minutes, so this is a good time to take a quick break. In the next section, we'll analyze the results." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Comparing two-way shooting to biased shooting\n", - "\n", - "Let's open the output file and analyze the results." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "std_two_way = paths.Storage(\"2_way.nc\", mode='r')\n", - "biased_two_way = paths.Storage(\"biased.nc\", mode='r')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "two_way_scheme = std_two_way.schemes['2_way']\n", - "biased_scheme = biased_two_way.schemes['biased_shooting']" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# each of these will take about a minute\n", - "two_way_scheme.move_summary(std_two_way.steps)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "biased_scheme.move_summary(biased_two_way.steps)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "As you can see, the acceptance rate for uniform 2-way shooting is pretty low. We improve this significantly by using the Gaussian biased shooting (keeping in mind that 25% is the theoretical maximum acceptance rate for the approach we use here.)\n", - "\n", - "Next we will plot the shooting points with uniform shooting point selection and with the Gaussian biased shooting point selection. First we write and use a little function to extract the shooting points." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "def get_shooting_points(steps):\n", - " \"\"\"Function to extract x,y positions of all shooting points\"\"\"\n", - " shooting_snaps = [step.change.canonical.details.shooting_snapshot for step in steps]\n", - " xy = [snap.xyz[0][:2] for snap in shooting_snaps] # get x and y positions\n", - " return tuple(zip(*xy)) # [[x1, y1], [x2, y2]] => ([x1, x2], [y1, y2])" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# the 0th step saves initial conditions, so we only have shooting moves as of step 1\n", - "std_x, std_y = get_shooting_points(std_two_way.steps[1:])\n", - "biased_x, biased_y = get_shooting_points(biased_two_way.steps[1:])" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The next two plots show the location of shooting points in the when using uniform shooting (first plot) and when using Gaussian biased shooting (second plot). Note that the shooting points fall in a much narrower range in the Gaussian biased plots. This is also why we get a better acceptance rate!" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plot.plot()\n", - "plt.scatter(std_x, std_y, c=range(len(std_x)), cmap='rainbow');" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plot.plot()\n", - "plt.scatter(biased_x, biased_y, c=range(len(std_x)), cmap='rainbow');" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## (Embarassingly) parallel TIS without replica exchange\n", - "\n", - "The OPS `DefaultScheme` is designed to provide reasonable default behaviors for TIS. These include replica exchange moves, path reversal moves, the minus interface move, as well as shooting. In general, replica exchange TIS is a much more efficient way to sample than TIS without replica exchange. However, replica exchange TIS is much harder to parallelize, because the duration of each trajectory is not known before running the trajectory.\n", - "\n", - "Therefore, in some cases you may want to sample each interface independently. This allows a naïve parallelization, since each interface is its own simulation.\n", - "\n", - "Here, we will make a RETIS scheme that is similar to the default scheme, except it doesn't include the minus interface move." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "interfaces = paths.VolumeInterfaceSet(cv=cv, minvals=float(\"-inf\"), \n", - " maxvals=[-0.60, -0.5, -0.4, -0.3, -0.2])\n", - "tis_network = paths.MISTISNetwork([(state_A, interfaces, state_B)]).named(\"tis\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# this is basically the DefaultScheme without the minus interface move\n", - "retis_scheme = paths.MoveScheme(network=tis_network).named(\"retis\")\n", - "retis_scheme.append([\n", - " strategies.OneWayShootingStrategy(engine=engine),\n", - " strategies.PathReversalStrategy(),\n", - " strategies.NearestNeighborRepExStrategy(),\n", - " strategies.OrganizeByMoveGroupStrategy()\n", - "])" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The `OneWayShootingStrategy` includes an `ensembles`, which selects specific ensembles. The (normal TIS) ensembles sampled by the TIS network are in the attribute `sampling_ensembles`. (Aside: other ensembles, such as the minus interface ensembles, are in the `special_ensembles` attribute.) This means that we can create a strategy " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "ens_0_strategy = strategies.OneWayShootingStrategy(\n", - " ensembles=[tis_network.sampling_ensembles[0]],\n", - " engine=engine\n", - ")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "From here, we can make a " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "scheme_0 = paths.MoveScheme(tis_network).named(\"scheme_0\")\n", - "scheme_1 = paths.MoveScheme(tis_network).named(\"scheme_1\")\n", - "scheme_2 = paths.MoveScheme(tis_network).named(\"scheme_2\")\n", - "scheme_3 = paths.MoveScheme(tis_network).named(\"scheme_3\")\n", - "scheme_4 = paths.MoveScheme(tis_network).named(\"scheme_4\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "global_strategy = strategies.OrganizeByMoveGroupStrategy()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# YOUR TURN: Make the correct strategies and append things to the scheme\n", - "# 1. Create a OneWayShootingStrategy for each ensemble\n", - "# 2. Append the global_strategy and the appropriate shooting strategy to each scheme\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Running TIS\n", - "\n", - "Again, we'll use the command line interface to run the TIS. So the first stage is to save the relevant things to a file, and then we can use the `--scheme` option in the `pathsampling` command to select which scheme to run.\n", - "\n", - "First, we check that all our schemes are correct:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "_ = retis_scheme.initial_conditions_from_trajectories(initial_conditions)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "_ = scheme_0.initial_conditions_from_trajectories(initial_conditions)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "_ = scheme_1.initial_conditions_from_trajectories(initial_conditions)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "_ = scheme_2.initial_conditions_from_trajectories(initial_conditions)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "_ = scheme_3.initial_conditions_from_trajectories(initial_conditions)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "_ = scheme_4.initial_conditions_from_trajectories(initial_conditions)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# saving everything will take a few minutes\n", - "parallel_setup = paths.Storage(\"parallel_setup.nc\", mode='w')\n", - "parallel_setup.tags['initial_conditions'] = initial_conditions\n", - "parallel_setup.save(retis_scheme)\n", - "parallel_setup.save(scheme_0)\n", - "parallel_setup.save(scheme_1)\n", - "parallel_setup.save(scheme_2)\n", - "parallel_setup.save(scheme_3)\n", - "parallel_setup.save(scheme_4)\n", - "parallel_setup.close()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "In OPS, each individual move, such as an attempt to swap a specific pair of replicas, counts as a Monte Carlo step. So in order to make a fair comparison of the approaches with an without replica exchange, we want to ensure that they both have about the same number of shooting moves for each ensemble.\n", - "\n", - "However, the `MoveScheme` can give us an estimate of how many total moves are required to get a certain number of moves of a certain mover. To get 250 trials of the 0th (only!) shooting mover in `scheme_0`, how many total steps do we need?" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "scheme_0.n_steps_for_trials(scheme_0.movers['shooting'][0], 250)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "That answer is probably pretty obvious. But what about our RETIS scheme? How many total steps to we need (on average) to get 250 trials of the 0th shooting mover from that move scheme?" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# YOUR TURN: Answer the question above\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "This should be a significantly larger number, and it is due to the many replica exchange and path reversal moves in that move scheme." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Now let's run the simulations. First, we equilibrate the initial condition. You can do this with:\n", - "\n", - "```\n", - "$ openpathsampling equilibrate parallel_setup.nc -o retis_equil.nc --scheme retis --extra-steps 50\n", - "```\n", - "\n", - "That will first run until the first decorrelated path (no frames in common with the initial trajectory), and then run an additional 50 MC steps. The results will be saved in the `scheme_0_equil.nc` file.\n", - "\n", - "Then you can run the full simulation with:\n", - "```\n", - "$ openpathsampling pathsampling retis_equil.nc -o retis.nc -n $NSTEPS > retis.out &\n", - "```\n", - "where you should replace `$NSTEPS` with the number of steps you found for RETIS above. \n", - "\n", - "If you're not familiar, the `> retis.out` redirects the output to the file `retis.out` (you can `cat retis.out` to see progress updates), and the `&` at the end of the command forces the command to run in the background, so that you can issue more commands from the same command line (i.e., run multiple things in parallel).\n", - "\n", - "* Why don't you need to specify a `--scheme` with the second command? (Hint: use the `openpathsampling contents` command on `retis_equil.nc` and `parallel_setup.nc`. How many move schemes are saved in each?) \n", - "\n", - "\n", - "Do the same for `scheme_0`, `scheme_1`, `scheme_2`, `scheme_3` and `scheme_4`, running the path sampling with 250 steps each. In this way, you will be running all 5 interfaces in parallel." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Analyzing the TIS simulations" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from openpathsampling.analysis import tis" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "storage_0 = paths.Storage('scheme_0.nc', mode='r')\n", - "storage_1 = paths.Storage('scheme_1.nc', mode='r')\n", - "storage_2 = paths.Storage('scheme_2.nc', mode='r')\n", - "storage_3 = paths.Storage('scheme_3.nc', mode='r')\n", - "storage_4 = paths.Storage('scheme_4.nc', mode='r')\n", - "storage_retis = paths.Storage('retis.nc', mode='r')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "scheme_0 = storage_0.schemes['scheme_0']\n", - "scheme_1 = storage_1.schemes['scheme_1']\n", - "scheme_2 = storage_2.schemes['scheme_2']\n", - "scheme_3 = storage_3.schemes['scheme_3']\n", - "scheme_4 = storage_4.schemes['scheme_4']\n", - "scheme_retis = storage_retis.schemes['retis']" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Comparing the move summaries" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": false - }, - "outputs": [], - "source": [ - "# takes a few minutes\n", - "scheme_retis.move_summary(storage_retis.steps)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "scheme_retis.move_summary(movers='shooting')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "scheme_retis.move_summary(movers='repex')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# each of these takes a minute or so\n", - "scheme_0.move_summary(storage_0.steps)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "scheme_1.move_summary(storage_1.steps)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "scheme_2.move_summary(storage_2.steps)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "scheme_3.move_summary(storage_3.steps)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "scheme_4.move_summary(storage_4.steps)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "scheme_0.move_summary(movers='shooting')" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### TIS analysis (crossing probabilities, etc.)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We don't actually have the flux here, so we can't calculate the actual rates. However, we can create a fake flux that says that the flux through the out of state $A$ and through the innermost interface is `1.0`. This allows us to use the rest of the `StandardTISAnalysis` object. It just means that the rate that gets reported is actually the total transition probability.\n", - "\n", - "You can get the actual flux either from including a minus interface move in your TIS simulation, or from using direct MD. The `paths.TrajectoryTransitionAnalysis` class will analyze existing MD trajectories, or the `paths.DirectSimulation` class can run MD and analyze the flux on the fly." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "network = storage_retis.networks[0]\n", - "state_A = storage_retis.volumes['A']\n", - "interface_0 = network.sampling_transitions[0].interfaces[0]\n", - "fake_flux = tis.DictFlux({(state_A, interface_0): 1.0})" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Finally, we assemble the `StandardTISAnalysis` and perform the analysis:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "%%time\n", - "# takes about 2 minutes\n", - "retis_analysis = tis.StandardTISAnalysis(\n", - " network=network,\n", - " flux_method=fake_flux,\n", - " max_lambda_calcs={t: {'bin_width': 0.025, 'bin_range': (-0.6, 0.6)}\n", - " for t in network.sampling_transitions},\n", - " steps=storage_retis.steps\n", - ")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Currently, the parallel analysis needs an extra step to run correctly. We need to create the `weighted_trajectories` object from the steps, and then perform the overall analysis using that as the input, instead of the steps themselves." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "%%time\n", - "# currently we need to manually join the weighted trajectories from each storage\n", - "# Future versions of OPS will simplify this\n", - "weighted_trajectories = {}\n", - "storages = [storage_0, storage_1, storage_2, storage_3, storage_4]\n", - "for storage, ensemble in zip(storages, network.sampling_ensembles):\n", - " weighted_trajectories.update(\n", - " tis.core.steps_to_weighted_trajectories(storage.steps, [ensemble])\n", - " )" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "%%time\n", - "# this will take a few minutes\n", - "state_A = storage_0.volumes['A']\n", - "interface_0 = network.sampling_transitions[0].interfaces[0]\n", - "fake_flux = tis.DictFlux({(state_A, interface_0): 1.0})\n", - "parallel_analysis = tis.StandardTISAnalysis(\n", - " network=network,\n", - " flux_method=fake_flux,\n", - " max_lambda_calcs={t: {'bin_width': 0.025, 'bin_range': (-0.6, 0.6)}\n", - " for t in network.sampling_transitions},\n", - " combiners={t.interfaces: paths.numerics.WHAM(cutoff=0.01,\n", - " interfaces=t.interfaces.lambdas)\n", - " for t in network.sampling_transitions}\n", - ")\n", - "parallel_analysis.results['flux'] = fake_flux.calculate('foo')\n", - "parallel_analysis.results = parallel_analysis.from_weighted_trajectories(weighted_trajectories)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Plotting the crossing probabilities\n", - "\n", - "One of the spot-checks to see if your simulation is converged is to plot the crossing probabilities functions. For each ensemble, the `StandardTISAnalysis` calculates a crossing probability along the order parameter, defined as the fraction of paths in that ensemble that reach at least the given value on the $x$ axis. As such, the crossing probability is always 1 for values less than the cutoff for the interface. Additionally, two ensemble crossing probabilities should never cross; the one from an outer interface should always be higher at a given value of the order parameter than one from an inner interface.\n", - "\n", - "There is also the *total* crossing probability, which is generated by using a histogram combining algorithm (usually WHAM) to combine the individual ensemble crossing probabilities into a good estimate for the true crossing probability (from the innermost interface). Like all crossing probabilities, this should be monotonically decreasing; if it is not, that is a sign of insufficient sampling.\n", - "\n", - "Since the y-axis is probability, and we're looking at rare events, we frequently plot crossing probabilities on a semi-log plot." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "for ensemble in network.transitions[(state_A, state_B)].ensembles:\n", - " crossing = retis_analysis.crossing_probability(ensemble)\n", - " label = \"Interface at $x$={:3.2f}\".format(ensemble.lambda_i)\n", - " plt.plot(crossing.x, crossing, label=label)\n", - "\n", - "tcp_AB = retis_analysis.total_crossing_probability[(state_A, state_B)]\n", - "plt.plot(tcp_AB.x, tcp_AB, lw=2, color='k', label=\"Total crossing probability\")\n", - "plt.legend()\n", - "plt.yscale('log')\n", - "plt.xlabel('$x$')\n", - "plt.ylabel('Crossing probability');" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "for ensemble in network.transitions[(state_A, state_B)].ensembles:\n", - " crossing = parallel_analysis.crossing_probability(ensemble)\n", - " label = \"x={:3.2f}\".format(ensemble.lambda_i)\n", - " plt.plot(crossing.x, crossing, label=label)\n", - "\n", - "tcp_AB = parallel_analysis.total_crossing_probability[(state_A, state_B)]\n", - "plt.plot(tcp_AB.x, tcp_AB, lw=2, color='k', label=\"Total crossing probability\")\n", - "plt.legend()\n", - "plt.yscale('log')\n", - "plt.xlabel('$x$')\n", - "plt.ylabel('Crossing probability');" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.7.3" - }, - "toc": { - "base_numbering": 1, - "nav_menu": {}, - "number_sections": true, - "sideBar": true, - "skip_h1_title": true, - "title_cell": "Table of Contents", - "title_sidebar": "Contents", - "toc_cell": false, - "toc_position": { - "height": "calc(100% - 180px)", - "left": "10px", - "top": "150px", - "width": "278.594px" - }, - "toc_section_display": true, - "toc_window_display": true - }, - "varInspector": { - "cols": { - "lenName": 16, - "lenType": 16, - "lenVar": 40 - }, - "kernels_config": { - "python": { - "delete_cmd_postfix": "", - "delete_cmd_prefix": "del ", - "library": "var_list.py", - "varRefreshCmd": "print(var_dic_list())" - }, - "r": { - "delete_cmd_postfix": ") ", - "delete_cmd_prefix": "rm(", - "library": "var_list.r", - "varRefreshCmd": "cat(var_dic_list()) " - } - }, - "types_to_exclude": [ - "module", - "function", - "builtin_function_or_method", - "instance", - "_Feature" - ], - "window_display": false - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} diff --git a/5_custom_shooting_setup.ipynb b/5_custom_shooting_setup.ipynb new file mode 100644 index 0000000..13f90dc --- /dev/null +++ b/5_custom_shooting_setup.ipynb @@ -0,0 +1,494 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Customizing shooting moves: Setup\n", + "\n", + "The tutorial focuses on how to customize shooting moves in the move scheme. We discuss three major reasons to do that:\n", + "\n", + "1. To use a different kind of shooting move, such as two-way shooting, when needed.\n", + "2. To use a different shooting point selector, such as a Gaussian biased selector, to get better efficiency.\n", + "3. To create move schemes that only sample part of the network, for example, to perform TIS with each ensemble sampled in parallel.\n", + "\n", + "In addition, you'll learn about:\n", + "\n", + "* Creating a \"setup\" file with engine and other simulation information.\n", + "* Using the OpenPathSampling command line interface (CLI) to run simulations from the setup file.\n", + "\n", + "Note that you'll need to install the CLI, which can be done with either `pip install openpathsampling-cli` or `conda install -c conda-forge openpathsampling-cli`. (Already installed in Binder.)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "\n", + "import openpathsampling as paths\n", + "from openpathsampling import strategies" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## A simple toy model example system\n", + "\n", + "In this notebook, we'll use a simple toy system. For the most part, the ideas here directly generalize to other engines, although we'll make a few comments where the units associated with OpenMM require special care.\n", + "\n", + "We'll start by loading a number things from the storage file:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import sys\n", + "if sys.version_info < (3, 8):\n", + " f_version = \"\"\n", + "else:\n", + " f_version = \"_38\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "storage = paths.Storage(f\"./inputs/2_state_toy{f_version}.nc\", mode='r')\n", + "state_A = storage.volumes['A']\n", + "state_B = storage.volumes['B']\n", + "cv = storage.cvs['x']\n", + "engine = storage.engines['toy_engine']\n", + "initial_conditions = storage.tags['initial_conditions']" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The toy model here uses a simple 2D potential energy surface, described by:\n", + "\n", + "$$V(x, y) = x^6 + y^6 - e^{-12(x+0.6)^2 - 5 y^2} - e^{-12(x-0.6)^2 - 5 y^2}$$\n", + "\n", + "We can visualize the toy engine using the `toy_plot_helpers` included in this repository:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%run toy_plot_helpers.py\n", + "\n", + "pes = engine.topology.pes\n", + "\n", + "plot = ToyPlot()\n", + "plot.contour_range = np.arange(-1.5, 1.0, 0.1)\n", + "plot.add_pes(pes)\n", + "fig = plot.plot() " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We run this simulation at a temperature of $T=0.1$, so the barrier between those wells is about $10\\ k_\\text{B} T$." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Changing the shooting move type (e.g., two-way shooting)\n", + "\n", + "TPS is usually introduced with two-way shooting. In that algorithm, you select a frame of the trajectory, make some modification (typically changing the velocities in some way consistent with the thermodynamic ensemble), and then integrate the equations of motion forward and backward.\n", + "\n", + "However, this is problematic if the velocity memory is short compared to your TPS trajectory length, as is often the case in condensed phase systems. In this case, each shot might as well be a committor shot. The ideal place to shoot from is the where the committor is 1/2, where each shot has only a 50% chance of landing in the correct state. But this gives you a maximum of a 25% acceptance rate.\n", + "\n", + "To get around this, we typically use the one-way shooting algorithm when working with condensed matter systems. One-way shooting works under the assumption that you are using a stochastic integrator. If that's the case, then instead of changing the velocities at the shooting point, you can use the fact that the integrator will generate a new sequence of random numbers. Therefore, you create a new trajectory by only running in one direction, and that trajectory is still a valid, physical trajectory.\n", + "\n", + "In general, OpenPathSampling defaults to using one-way shooting. However, there are circumstances where you might want to use the older two-way shooting algorithm. For example, if the velocity memory is shorter than the period between saved frames, it might be better to use two-way shooting. If you intend to use deterministic dynamics, then one-way shooting is not possible.\n", + "\n", + "This part of the tutorial will show you how to use two-way shooting instead of one-way shooting in OPS. There are two ways of doing this: either replace the existing one-way shooting strategy, or create a new move scheme from scratch." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "tps_network = paths.TPSNetwork(initial_states=state_A, final_states=state_B)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Replacing parts of an existing move scheme\n", + "\n", + "Frequently, the easiest way to modify a move scheme (especially a complex move scheme) is to replace the parts that you want to change. In this case, we can start with the move scheme that only includes 1-way shooting, and work from there." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "modified_scheme = paths.OneWayShootingMoveScheme(\n", + " network=tps_network,\n", + " engine=engine\n", + ").named(\"2_way_from_1_way\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For our two-way shooting, we'll need to create a `TwoWayShootingStrategy`. This will require a modifier; to keep things simple, we'll completely randomize velocities (consistent with a given temperature) using `paths.RandomVelocities`. In order to be consistent between engines, `RandomVelocities` takes its input as the inverse temperature, $\\beta = 1/(k_\\text{B}T)$.\n", + "\n", + "For the toy engine, we can obtain the temperature from `engine.integ.temperature`, and we work in units where $k_\\text{B}=1$, so it is easy to calculate $\\beta$. For other engines, you'll need to use the correct value of $k_\\text{B}$. For OpenMM, you also need to worry about units: use `simtk.unit.BOLTZMANN_CONSTANT_kB`.\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true, + "tags": [] + }, + "outputs": [], + "source": [ + "print(engine.integ.temperature)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# YOUR TURN: Set beta correctly (yes, this is as easy as you think).\n", + "# beta = ... # fill in the ellipsis and uncomment this line" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "modifier = paths.RandomVelocities(beta=beta, engine=engine)\n", + "shooting_strategy = strategies.TwoWayShootingStrategy(modifier=modifier, engine=engine)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To modify a move scheme, just `append` new strategies. In this case, the new strategy for the group of moves called `'shooting'` will be overwritten. If you gave the `TwoWayShootingStrategy` a different string for its `group` argument, then you would create a second group of movers, and each move would have a 50/50 chance of using one-way shooting or of using two-way shooting." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "modified_scheme.append(shooting_strategy)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Building a move scheme from scratch\n", + "\n", + "In this particular case, our overall move scheme will not be too complicated. In fact, let's start by looking at the at the code for " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import inspect\n", + "from IPython.display import Code\n", + "\n", + "Code(inspect.getsource(paths.OneWayShootingMoveScheme))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The output above is the full source code for the `OneWayShootingMoveScheme`. You can see that all it does is pass the network up to the superclass's initialization, and then `append` two move strategies to itself. So we can recreate this:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "two_way_scheme = paths.MoveScheme(network=tps_network).named(\"2_way\")\n", + "global_strategy = strategies.OrganizeByMoveGroupStrategy()\n", + "two_way_scheme.append([shooting_strategy, global_strategy])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Changing the shooting point selection\n", + "\n", + "A common problem in one-way shooting is that paths do not decorrelate quickly enough. This typically happens when there is a significant barrier withing the transition region, such that shooting points before the barrier always go back to the initial state (so only backward shots are accepted) and shooting points after the barrier always to do the final state (so only forward shots are accepted).\n", + "\n", + "The same issue manifests in two-way shooting as a significant decrease in the acceptance rate for the shooting move. One side of the barrier always creates $A\\to A$ trials, while the other side always creates $B\\to B$ trials. With fewer $A\\to B$ trials, the acceptance rate drops." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "xvals = np.arange(-1.0, 1.0, 0.01)\n", + "\n", + "def plot_gaussian_bias(selector):\n", + " alpha = selector.alpha\n", + " x_0 = selector.l_0\n", + " gaussians = np.exp(-alpha*(xvals - x_0)**2)\n", + " plt.plot(xvals, gaussians)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "biased_shooting_scheme = paths.MoveScheme(network=tps_network).named(\"biased_shooting\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "gaussian_sel = paths.GaussianBiasSelector(cv, alpha=100.0, l_0=0.0)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plot_gaussian_bias(gaussian_sel)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "paths.GaussianBiasSelector?" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "* The Gaussian selector has two parameters, $\\alpha$ and $l_0$. What do each of these parameters change?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# YOUR TURN: Create a move scheme using a two-way shooting strategy with this selector.\n", + "# 1. Create a two-way shooting strategy that uses this Gaussian selector. (Use the \n", + "# selector keyword argument of TwoWayShootingStrategy.)\n", + "# 2. Append that strategy and the global_strategy to the biased_shooting_scheme" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Comparing uniform vs Gaussian bias shooting point selection\n", + "\n", + "Now we have two different move schemes that sample the same TPS network. Both use two-way shooting with complete velocity randomization; the only difference is that one uses uniform shooting point selection whereas the other uses a Gaussian biased selection.\n", + "\n", + "Let's see how these two approaches compare with each other -- we'll run them for the same number of steps, and in the next notebook, we'll compare the acceptance rates.\n", + "\n", + "First, let's verify that the initial conditions will work for the move schemes we've created." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "_ = two_way_scheme.initial_conditions_from_trajectories(initial_conditions)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "_ = biased_shooting_scheme.initial_conditions_from_trajectories(initial_conditions)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Creating a setup file for the OpenPathSampling CLI\n", + "\n", + "The files that OPS stores to always include details of how the simulation is run, as well as storing the newly generated data. However, we can also store only the information needed to run the simulation. That's the idea behind a \"setup\" file for an OPS simulation.\n", + "\n", + "This is the same storage you've used before, but we're using it in a different way. One of the main reasons to use \"setup\" files is that path sampling is rarely just a single computational run. Often you run several copies of your simulation simultaneously. Or perhaps after running TPS you'll want to do a committor analysis with the same state definitions. Using a setup file ensures that the objects you use are identical (no accidental copy-paste errors). Because the OPS analysis tools can also verify that these objects are the same, this also makes analysis easier.\n", + "\n", + "In our case, we'll use a single setup file that stores our initial conditions and both of our move schemes (and therefore, the TPS network, state definitions, etc.) Then we'll use that one setup file to run both of our simulations." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "shooting_setup = paths.Storage(\"shooting_setup.nc\", mode='w')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# saving everything will take a few minutes\n", + "shooting_setup.tags['initial_conditions'] = initial_conditions\n", + "shooting_setup.save(two_way_scheme)\n", + "shooting_setup.save(biased_shooting_scheme)\n", + "shooting_setup.close()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Running simulations with the CLI\n", + "\n", + "If you've installed the OpenPathSampling CLI, you will have the command `openpathsampling` available in your shell. The `openpathsampling` command has several subcommands. You can see some of the common subcommands with `openpathsampling --help`, and you can also use the `--help` flag to get more information about a subcommand. For example, you can get more information about the `pathsampling` subcommand with `openpathsampling pathsampling --help`.\n", + "\n", + "* In a terminal, use the `openpathsampling contents` on the file `shooting_setup.nc` to list the contents of the file. How many snapshots are stored in that file? How many move schemes?" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now, use the `pathsampling` subcommand to run simulations with your two move schemes. Run at least 500 steps (`-n 500`), or run more if you'd like:\n", + "\n", + "```\n", + "$ openpathsampling pathsampling shooting_setup.nc -o 2_way.nc --scheme 2_way -n 500\n", + "\n", + "$ openpathsampling pathsampling shooting_setup.nc -o biased.nc --scheme biased_shooting -n 500\n", + "```\n", + "\n", + "Those simulations will take a few minutes, so this is a good time to take a quick break. In the next notebook, we'll analyze the results." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.10" + }, + "toc": { + "base_numbering": 1, + "nav_menu": {}, + "number_sections": true, + "sideBar": true, + "skip_h1_title": true, + "title_cell": "Table of Contents", + "title_sidebar": "Contents", + "toc_cell": false, + "toc_position": { + "height": "calc(100% - 180px)", + "left": "10px", + "top": "150px", + "width": "278.594px" + }, + "toc_section_display": true, + "toc_window_display": true + }, + "varInspector": { + "cols": { + "lenName": 16, + "lenType": 16, + "lenVar": 40 + }, + "kernels_config": { + "python": { + "delete_cmd_postfix": "", + "delete_cmd_prefix": "del ", + "library": "var_list.py", + "varRefreshCmd": "print(var_dic_list())" + }, + "r": { + "delete_cmd_postfix": ") ", + "delete_cmd_prefix": "rm(", + "library": "var_list.r", + "varRefreshCmd": "cat(var_dic_list()) " + } + }, + "types_to_exclude": [ + "module", + "function", + "builtin_function_or_method", + "instance", + "_Feature" + ], + "window_display": false + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/6_custom_shooting_analysis.ipynb b/6_custom_shooting_analysis.ipynb new file mode 100644 index 0000000..f3b0931 --- /dev/null +++ b/6_custom_shooting_analysis.ipynb @@ -0,0 +1,254 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "6d51955b-ba20-48b8-b32f-93fad50a316a", + "metadata": {}, + "source": [ + "# Customizing shooting moves: Analysis\n", + "\n", + "Now that you've run the simulation, let's compare the results of two-way shooting with uniform shooting point selection to those with biased shooting point selection.\n", + "\n", + "We'll open the output file and consider a few results. Although you can perform all the typical analyses you might do with TPS data (path densities, etc.), here we're just going to look at the acceptance probabilities and the shooting point distributions." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fb5cbd96-2a23-4598-ad0d-993a282f379e", + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "\n", + "import openpathsampling as paths" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2bb71581-0726-412e-aec3-962a3a819a50", + "metadata": {}, + "outputs": [], + "source": [ + "std_two_way = paths.Storage(\"2_way.nc\", mode='r')\n", + "biased_two_way = paths.Storage(\"biased.nc\", mode='r')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ff711075-4e90-43b0-9f14-c9cb5f41f5f5", + "metadata": {}, + "outputs": [], + "source": [ + "two_way_scheme = std_two_way.schemes['2_way']\n", + "biased_scheme = biased_two_way.schemes['biased_shooting']" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "643ba3ba-9b5f-4b73-8d1e-87de015e6e54", + "metadata": {}, + "outputs": [], + "source": [ + "two_way_scheme.move_summary(std_two_way.steps)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b550c1d3-f2a9-4e1b-b785-d49b0c9138ef", + "metadata": {}, + "outputs": [], + "source": [ + "biased_scheme.move_summary(biased_two_way.steps)" + ] + }, + { + "cell_type": "markdown", + "id": "36d9a273-bde7-4aa8-9108-122c892ec01e", + "metadata": {}, + "source": [ + "As you can see, the acceptance rate for uniform 2-way shooting with complete velocity randomization is pretty low. We improve this significantly by using the Gaussian biased shooting (keeping in mind that 25% is the theoretical maximum acceptance rate for the approach we use here.)\n", + "\n", + "Next we will plot the shooting points with uniform shooting point selection and with the Gaussian biased shooting point selection. First we write and use a little function to extract the shooting points." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "47bb27d3-4a0d-4063-ac8d-032f1275ae52", + "metadata": {}, + "outputs": [], + "source": [ + "def get_shooting_points(steps):\n", + " \"\"\"Function to extract x,y positions of all shooting points\"\"\"\n", + " shooting_snaps = [step.change.canonical.details.shooting_snapshot for step in steps]\n", + " xy = [snap.xyz[0][:2] for snap in shooting_snaps] # get x and y positions\n", + " return tuple(zip(*xy)) # [[x1, y1], [x2, y2]] => ([x1, x2], [y1, y2])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "56e2dc59-3cb9-4c3c-917b-670d5da84429", + "metadata": {}, + "outputs": [], + "source": [ + "# the 0th step saves initial conditions, so we only have shooting moves as of step 1\n", + "std_x, std_y = get_shooting_points(std_two_way.steps[1:])\n", + "biased_x, biased_y = get_shooting_points(biased_two_way.steps[1:])" + ] + }, + { + "cell_type": "markdown", + "id": "3aeae0e4-a8c9-461f-ab94-dbc8d2d2c555", + "metadata": {}, + "source": [ + "The next two plots show the location of shooting points in the when using uniform shooting (first plot) and when using Gaussian biased shooting (second plot). Note that the shooting points fall in a much narrower range in the Gaussian biased plots. This is also why we get a better acceptance rate!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a71c14cf-ae55-4f41-a3cd-52fa29c403e2", + "metadata": {}, + "outputs": [], + "source": [ + "%run toy_plot_helpers.py\n", + "# doesn't matter where we get the toy engine from; same in both\n", + "engine = std_two_way.engines['toy_engine']\n", + "pes = engine.topology.pes\n", + "plot = ToyPlot()\n", + "plot.contour_range = np.arange(-1.5, 1.0, 0.1)\n", + "plot.add_pes(pes)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "aefe10ec-c62a-474f-8052-8a94aac760f2", + "metadata": {}, + "outputs": [], + "source": [ + "# using c=range(len(std_x) and cmap='rainbow' here means that color corresponds to step order\n", + "plot.plot()\n", + "plt.scatter(std_x, std_y, c=range(len(std_x)), cmap='rainbow');" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "49071ee5-01be-4dc8-994f-493cbadbbca0", + "metadata": {}, + "outputs": [], + "source": [ + "plot.plot()\n", + "plt.scatter(biased_x, biased_y, c=range(len(std_x)), cmap='rainbow');" + ] + }, + { + "cell_type": "markdown", + "id": "345a5fb5-9bcb-4067-bcf2-3e26bc1ada85", + "metadata": {}, + "source": [ + "* How would you expect the shooting acceptance to change if we changed $l_0$, the center of the Gaussian? What useful information do we have about this potential energy surface that we might not have in a real system?\n", + "* What would be the effect of changing the Gaussian selector's $\\alpha$ value (Gaussian width parameter)?" + ] + }, + { + "cell_type": "markdown", + "id": "1a5b6c4c-2ac7-4778-b3ae-69800401a5e6", + "metadata": {}, + "source": [ + "Now let's make similar plots, but we'll plot the shooting points for accepted/rejected steps in different colors." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9807f40d-2882-4d53-8223-9ad8492f500c", + "metadata": {}, + "outputs": [], + "source": [ + "# create lists to plot the accepted shooting points\n", + "std_acc_steps = [step for step in std_two_way.steps[1:] if step.change.accepted]\n", + "std_acc_x, std_acc_y = get_shooting_points(std_acc_steps)\n", + "\n", + "biased_acc_steps = [step for step in biased_two_way.steps[1:] if step.change.accepted]\n", + "biased_acc_x, biased_acc_y = get_shooting_points(biased_acc_steps)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "522a3a56-c119-41b3-8a08-d0edb8e835c6", + "metadata": {}, + "outputs": [], + "source": [ + "# YOUR TURN: Do the same for the rejected shooting points.\n", + "\n", + "# std_rej_x, std_rej_y = ...\n", + "# biased_rej_x, biased_rej_y = ..." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "991d2207-c7c0-42df-9293-256e2f592594", + "metadata": {}, + "outputs": [], + "source": [ + "plot.plot()\n", + "plt.scatter(std_rej_x, std_rej_y, c='r', label='Rejected')\n", + "plt.scatter(std_acc_x, std_acc_y, c='b', label='Accepted')\n", + "plt.legend();" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b8c1a690-c668-4d9b-878c-cb3b6ab40a1c", + "metadata": {}, + "outputs": [], + "source": [ + "plot.plot()\n", + "plt.scatter(biased_rej_x, biased_rej_y, c='r', label='Rejected')\n", + "plt.scatter(biased_acc_x, biased_acc_y, c='b', label='Accepted')\n", + "plt.legend();" + ] + }, + { + "cell_type": "markdown", + "id": "ae380499-1341-447a-a523-53d365b8ec7b", + "metadata": {}, + "source": [ + "* Are the accepted are rejected shooting points distributed similarly with uniform shooting? What about with the biased selector?\n", + "* Given the plot of accepted and rejected shooting points for uniform shooting, what would you propose as a shooting point selection bias?" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.10" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/7_parallel_tis_setup.ipynb b/7_parallel_tis_setup.ipynb new file mode 100644 index 0000000..34d1977 --- /dev/null +++ b/7_parallel_tis_setup.ipynb @@ -0,0 +1,355 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "3df2ac2d-30f1-4010-ae44-d024bedb1544", + "metadata": {}, + "source": [ + "# (Embarrassingly) Parallel TIS vs. RETIS: Setup\n", + "\n", + "The OPS `DefaultScheme` is designed to provide reasonable default behaviors for TIS. These include replica exchange moves, path reversal moves, the minus interface move, as well as shooting. In general, replica exchange TIS is a more efficient way to sample than TIS without replica exchange. However, replica exchange TIS is much harder to parallelize, because the duration of each trajectory is not known before running the trajectory.\n", + "\n", + "Therefore, in some cases you may want to sample each interface independently. This allows a naïve parallelization, since each interface is its own simulation.\n", + "\n", + "In this notebook, we'll make move schemes for a custom version of RETIS (similar to the default scheme, but without the minus move) and for running each interface of a TIS simulation independently (in parallel). We'll use the OPS CLI to run it, and in the next notebook we'll compare the results.\n", + "\n", + "We'll use the same simple toy model as we used in notebooks 5 and 6." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5c1fb77d-673a-4b59-bfb9-5c048476b414", + "metadata": {}, + "outputs": [], + "source": [ + "import openpathsampling as paths\n", + "from openpathsampling import strategies" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1166098e", + "metadata": {}, + "outputs": [], + "source": [ + "import sys\n", + "if sys.version_info < (3, 8):\n", + " f_version = \"\"\n", + "else:\n", + " f_version = \"_38\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "de08a51c-cc74-4c8a-ba97-4049a17e8d65", + "metadata": {}, + "outputs": [], + "source": [ + "storage = paths.Storage(f\"./inputs/2_state_toy{f_version}.nc\", mode='r')\n", + "state_A = storage.volumes['A']\n", + "state_B = storage.volumes['B']\n", + "cv = storage.cvs['x']\n", + "engine = storage.engines['toy_engine']\n", + "initial_conditions = storage.tags['initial_conditions']" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "69347fa2-0532-413d-b2e3-adbfc6d5728c", + "metadata": {}, + "outputs": [], + "source": [ + "interfaces = paths.VolumeInterfaceSet(cv=cv, minvals=float(\"-inf\"), \n", + " maxvals=[-0.6, -0.45, -0.375, -0.3, -0.2])\n", + "tis_network = paths.MISTISNetwork([(state_A, interfaces, state_B)]).named(\"tis\")" + ] + }, + { + "cell_type": "markdown", + "id": "70bbba4c-d885-4d4c-a7ce-194f2c864447", + "metadata": {}, + "source": [ + "All move schemes need a \"global\" level strategy. This is usually the `OrganizeByMoveGroupStrategy`. We can re-use the same strategy in all our move schemes." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a3c99594-cb39-45a9-a9e8-feda08b49f4c", + "metadata": {}, + "outputs": [], + "source": [ + "global_strategy = strategies.OrganizeByMoveGroupStrategy()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f7399b89-9acd-4003-a961-05ecb0bd8882", + "metadata": {}, + "outputs": [], + "source": [ + "# this is basically the DefaultScheme without the minus interface move\n", + "retis_scheme = paths.MoveScheme(network=tis_network).named(\"retis\")\n", + "retis_scheme.append([\n", + " strategies.OneWayShootingStrategy(engine=engine),\n", + " strategies.PathReversalStrategy(),\n", + " strategies.NearestNeighborRepExStrategy(),\n", + " global_strategy\n", + "])" + ] + }, + { + "cell_type": "markdown", + "id": "0da16b81-da5b-4d89-995d-ec50514730cc", + "metadata": {}, + "source": [ + "The `OneWayShootingStrategy` includes an `ensembles` option, which selects specific ensembles to use. The (normal TIS) ensembles sampled by the TIS network are in the attribute `sampling_ensembles`. (Aside: other ensembles, such as the minus interface ensembles, are in the `special_ensembles` attribute.) This means that we can create a shooting strategy only samples a single ensemble (instead of the default, which is to sample all the `sampling_ensembles` in the network)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5fb5c252-17c7-4501-8832-f1bab84fbaf9", + "metadata": {}, + "outputs": [], + "source": [ + "ens_0_strategy = strategies.OneWayShootingStrategy(\n", + " ensembles=[tis_network.sampling_ensembles[0]],\n", + " engine=engine\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "9641b304-0a4a-45b0-8e5a-bb2934169086", + "metadata": {}, + "source": [ + "From here, we can make a moves scheme for each interface. I'll do that manually, although in practice you might come up with a loop to make this easier:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "184183fd-795e-4d4b-af1b-8714ba4700d3", + "metadata": {}, + "outputs": [], + "source": [ + "scheme_0 = paths.MoveScheme(tis_network).named(\"scheme_0\")\n", + "scheme_1 = paths.MoveScheme(tis_network).named(\"scheme_1\")\n", + "scheme_2 = paths.MoveScheme(tis_network).named(\"scheme_2\")\n", + "scheme_3 = paths.MoveScheme(tis_network).named(\"scheme_3\")\n", + "scheme_4 = paths.MoveScheme(tis_network).named(\"scheme_4\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1c46b83a-a25e-4cb8-8e80-683f33f14276", + "metadata": {}, + "outputs": [], + "source": [ + "# YOUR TURN: Make the correct strategies and append things to the scheme\n", + "# 1. Create a OneWayShootingStrategy for each ensemble\n", + "# 2. Append the global_strategy and the appropriate shooting strategy to each scheme" + ] + }, + { + "cell_type": "markdown", + "id": "48d17350-73c8-402b-b3ba-156cb7a611cc", + "metadata": {}, + "source": [ + "## Running TIS\n", + "\n", + "Again, we'll use the command line interface to run the TIS. So the first stage is to save the relevant things to a setup file, and then we can use the `--scheme` option in the `pathsampling` command to select which scheme to run.\n", + "\n", + "First, we check that all our schemes match our initial conditions:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5b44299f-23a4-4255-bfc2-12caef6ef913", + "metadata": {}, + "outputs": [], + "source": [ + "_ = retis_scheme.initial_conditions_from_trajectories(initial_conditions)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "95ed4cb2-4786-4822-a378-db3137263a74", + "metadata": {}, + "outputs": [], + "source": [ + "_ = scheme_0.initial_conditions_from_trajectories(initial_conditions)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "12e8d979-8a26-4b97-9f01-874c71ef58e7", + "metadata": {}, + "outputs": [], + "source": [ + "_ = scheme_1.initial_conditions_from_trajectories(initial_conditions)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f427e394-faba-45c3-8b54-da62990476f7", + "metadata": {}, + "outputs": [], + "source": [ + "_ = scheme_2.initial_conditions_from_trajectories(initial_conditions)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "48056996-0f5c-4c84-8058-7bd744c63093", + "metadata": {}, + "outputs": [], + "source": [ + "_ = scheme_3.initial_conditions_from_trajectories(initial_conditions)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c0cf2b3f-908e-4427-96d5-334cecc46717", + "metadata": {}, + "outputs": [], + "source": [ + "_ = scheme_4.initial_conditions_from_trajectories(initial_conditions)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "500beb8c-645c-4710-b38d-9120b163911e", + "metadata": {}, + "outputs": [], + "source": [ + "# saving everything will take a few minutes\n", + "parallel_setup = paths.Storage(\"parallel_setup.nc\", mode='w')\n", + "parallel_setup.tags['initial_conditions'] = initial_conditions\n", + "parallel_setup.save(retis_scheme)\n", + "parallel_setup.save(scheme_0)\n", + "parallel_setup.save(scheme_1)\n", + "parallel_setup.save(scheme_2)\n", + "parallel_setup.save(scheme_3)\n", + "parallel_setup.save(scheme_4)\n", + "parallel_setup.close()" + ] + }, + { + "cell_type": "markdown", + "id": "1bc1f3d5-94cd-4bf9-9695-3cccc24a5f7d", + "metadata": {}, + "source": [ + "In OPS, each individual move, such as an attempt to swap a specific pair of replicas, counts as a Monte Carlo step. So in order to make a fair comparison of the approaches with an without replica exchange, we want to ensure that they both have about the same number of shooting moves for each ensemble.\n", + "\n", + "However, the `MoveScheme` can give us an estimate of how many total moves are required to get a certain number of moves of a certain mover. To get 300 trials of the 0th (only!) shooting mover in `scheme_0`, how many total steps do we need?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cd514607-c420-47b1-8f69-d85ef9e7f1ac", + "metadata": {}, + "outputs": [], + "source": [ + "scheme_0.n_steps_for_trials(scheme_0.movers['shooting'][0], 300)" + ] + }, + { + "cell_type": "markdown", + "id": "4dcb1f5e-8fde-4763-b4a1-3779658fc8e7", + "metadata": {}, + "source": [ + "That answer is probably pretty obvious. But what about our RETIS scheme? How many total steps to we need (on average) to get 300 trials of the 0th shooting mover from that move scheme?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8df4dcf3-8712-4108-9865-e50b912d660d", + "metadata": {}, + "outputs": [], + "source": [ + "# YOUR TURN: Answer the question above" + ] + }, + { + "cell_type": "markdown", + "id": "b524c507-e7e6-4f03-af94-9f2585d1b8c3", + "metadata": {}, + "source": [ + "This should be a significantly larger number, and it is due to the many replica exchange and path reversal moves in that move scheme, as well as the fact that there are multiple ensembles to sample." + ] + }, + { + "cell_type": "markdown", + "id": "ce816c75-03c9-467c-b78b-5e1649d6e114", + "metadata": {}, + "source": [ + "Now let's run the simulations. First, we equilibrate the initial condition, since it is a transition trajectory, which is highly unlikely in inner interfaces. You can do this with:\n", + "\n", + "```\n", + "$ openpathsampling equilibrate parallel_setup.nc -o equil_retis.nc --scheme retis --extra-steps 50\n", + "```\n", + "\n", + "That will first run until the first decorrelated path (no frames in common with the initial trajectory), and then run an additional 50 MC steps. The results will be saved in the `equil_retis.nc` file.\n", + "\n", + "Then you can run the full simulation with:\n", + "```\n", + "$ openpathsampling pathsampling equil_retis.nc -o retis.nc -n $NSTEPS > retis.out &\n", + "```\n", + "where you should replace `$NSTEPS` with the number of steps you found for RETIS above. \n", + "\n", + "If you're not familiar, the `> retis.out` redirects the output to the file `retis.out` (you can `tail retis.out` to see progress updates), and the `&` at the end of the command forces the command to run in the background, so that you can issue more commands from the same command line (i.e., run multiple things in parallel).\n", + "\n", + "* Why don't you need to specify a `--scheme` with the second command? (Hint: use the `openpathsampling contents` command on `retis_equil.nc` and `parallel_setup.nc`. How many move schemes are saved in each?) \n", + "\n", + "\n", + "Do the same for `scheme_0`, `scheme_1`, `scheme_2`, `scheme_3` and `scheme_4`, running the path sampling with 300 steps each. In this way, you will be running all 5 interfaces in parallel." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ed256067", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.10" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/8_parallel_tis_analysis.ipynb b/8_parallel_tis_analysis.ipynb new file mode 100644 index 0000000..ff5fa7d --- /dev/null +++ b/8_parallel_tis_analysis.ipynb @@ -0,0 +1,424 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "d975f240-bdca-458f-acee-d9dd72e04ade", + "metadata": {}, + "source": [ + "# (Embarrassingly) Parallel TIS vs. RETIS: Analysis\n", + "\n", + "Now that you've run your embarrassingly parallel simulation and your RETIS simulation, let's compare the results. We'll check if the simulations seem consistent by looking at the acceptance of the shooting moves in each ensemble (you could also compare path length distributions or even do per-interface path density plots to give a more thorough check.) Then we'll look at the crossing probabilities to see if they seem to be converged (and perhaps whether one seems closer to converged than the other)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4504965b-df8b-4130-a19f-3f9ab1e5e93f", + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import openpathsampling as paths\n", + "from openpathsampling.analysis import tis" + ] + }, + { + "cell_type": "markdown", + "id": "15588f71-805d-4374-aa9b-580d8f1a4c65", + "metadata": {}, + "source": [ + "We'll start out by loading the output files we generated, and from each file we load in the move scheme that it used." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "73d2926f-7f26-44cd-85cf-b445156efeac", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "storage_0 = paths.Storage('scheme_0.nc', mode='r')\n", + "storage_1 = paths.Storage('scheme_1.nc', mode='r')\n", + "storage_2 = paths.Storage('scheme_2.nc', mode='r')\n", + "storage_3 = paths.Storage('scheme_3.nc', mode='r')\n", + "storage_4 = paths.Storage('scheme_4.nc', mode='r')\n", + "storage_retis = paths.Storage('retis.nc', mode='r')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "32cebb28-471c-467e-bcdf-89b25a9bb19a", + "metadata": {}, + "outputs": [], + "source": [ + "scheme_0 = storage_0.schemes['scheme_0']\n", + "scheme_1 = storage_1.schemes['scheme_1']\n", + "scheme_2 = storage_2.schemes['scheme_2']\n", + "scheme_3 = storage_3.schemes['scheme_3']\n", + "scheme_4 = storage_4.schemes['scheme_4']\n", + "scheme_retis = storage_retis.schemes['retis']" + ] + }, + { + "cell_type": "markdown", + "id": "71e4b708-a37f-4cf6-83d2-d3d5a80323fb", + "metadata": {}, + "source": [ + "## Comparing the acceptance rates\n", + "\n", + "Recall that we can get per-mover acceptance by passing the appropriate string to the `movers` keyword of the `move_summary`:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "822feb97-0d21-4bb7-bb85-77dbfd7b86ef", + "metadata": {}, + "outputs": [], + "source": [ + "scheme_retis.move_summary(storage_retis.steps)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "727e1c33-0679-422b-bb3c-5b9435e1fad7", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "scheme_retis.move_summary(movers='shooting')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "baf70277-3c76-440b-a1d1-1e6d6d3ff4db", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "scheme_retis.move_summary(movers='repex')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "eb6b9de7-35e7-419e-94ec-2b0928563196", + "metadata": {}, + "outputs": [], + "source": [ + "scheme_0.move_summary(storage_0.steps)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ef124f49-ac95-4149-a63f-d0e1309b610b", + "metadata": {}, + "outputs": [], + "source": [ + "scheme_1.move_summary(storage_1.steps)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "704453f9-8d5f-4a7b-8bd4-6d938dfacdea", + "metadata": {}, + "outputs": [], + "source": [ + "scheme_2.move_summary(storage_2.steps)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "135c7d76-2213-4280-8e8b-be842ecfa3e4", + "metadata": {}, + "outputs": [], + "source": [ + "scheme_3.move_summary(storage_3.steps)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2c769c6f-b98d-432c-a23d-d49f236589ea", + "metadata": {}, + "outputs": [], + "source": [ + "scheme_4.move_summary(storage_4.steps)" + ] + }, + { + "cell_type": "markdown", + "id": "36962444-6bd1-40ef-8751-b70dbf181198", + "metadata": {}, + "source": [ + "* Is there a major difference in the shooting move acceptance for any ensemble? Would you expect there to be one?" + ] + }, + { + "cell_type": "markdown", + "id": "043d97eb-17e6-4520-854d-34cc3d670a9f", + "metadata": {}, + "source": [ + "## TIS analysis (crossing probabilities, etc.)" + ] + }, + { + "cell_type": "markdown", + "id": "9e58c4ea-0ef5-43ad-8e6a-c96337bb2f71", + "metadata": {}, + "source": [ + "We don't have the flux here, so we can't calculate the actual rates. However, we can create a fake flux that says that the flux through the out of state $A$ and through the innermost interface is `1.0`. This allows us to use the rest of the `StandardTISAnalysis` object. It just means that the rate that gets reported is actually the total transition probability.\n", + "\n", + "You can get the actual flux either from including a minus interface move in your TIS simulation, or from using direct MD. The `paths.TrajectoryTransitionAnalysis` class will analyze existing MD trajectories, or the `paths.DirectSimulation` class can run MD and analyze the flux on the fly." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4ca161ea-a4a5-4307-84ab-c8c189e120a0", + "metadata": {}, + "outputs": [], + "source": [ + "# because we used a setup file, the netword/state/interface are the same in both\n", + "network = storage_retis.networks[0]\n", + "state_A = storage_retis.volumes['A']\n", + "state_B = storage_retis.volumes['B']\n", + "interface_0 = network.sampling_transitions[0].interfaces[0]\n", + "fake_flux = tis.DictFlux({(state_A, interface_0): 1.0})" + ] + }, + { + "cell_type": "markdown", + "id": "6ae2ed96-3d8a-439c-98dd-0165d0342ef7", + "metadata": {}, + "source": [ + "Finally, we assemble the `StandardTISAnalysis` and perform the analysis:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ff24a22f-c0b3-4135-8e67-c7c7d054879f", + "metadata": {}, + "outputs": [], + "source": [ + "%%time\n", + "retis_analysis = tis.StandardTISAnalysis(\n", + " network=network,\n", + " flux_method=fake_flux,\n", + " max_lambda_calcs={t: {'bin_width': 0.025, 'bin_range': (-0.6, 0.6)}\n", + " for t in network.sampling_transitions},\n", + " combiners={t.interfaces: paths.numerics.WHAM(cutoff=0.01, # lower cutoff, default is 0.05\n", + " interfaces=t.interfaces.lambdas)\n", + " for t in network.sampling_transitions},\n", + " steps=storage_retis.steps\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "8be276a9-82c2-45c0-8205-e35321eb5ae1", + "metadata": {}, + "source": [ + "Currently, the parallel analysis needs an extra step to run correctly. We need to create the `weighted_trajectories` object from the steps, and then perform the overall analysis using that as the input, instead of the steps themselves." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9a86866f-59e4-4ba1-8ab7-06643fc8a86d", + "metadata": {}, + "outputs": [], + "source": [ + "# currently we need to manually join the weighted trajectories from each storage\n", + "# Future versions of OPS will simplify this\n", + "weighted_trajectories = {}\n", + "storages = [storage_0, storage_1, storage_2, storage_3, storage_4]\n", + "for storage, ensemble in zip(storages, network.sampling_ensembles):\n", + " weighted_trajectories.update(\n", + " tis.core.steps_to_weighted_trajectories(storage.steps, [ensemble])\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7f974b4d-7f1a-4da4-8274-b1ae4656c54b", + "metadata": {}, + "outputs": [], + "source": [ + "#state_A = storage_0.volumes['A']\n", + "#interface_0 = network.sampling_transitions[0].interfaces[0]\n", + "#fake_flux = tis.DictFlux({(state_A, interface_0): 1.0})\n", + "parallel_analysis = tis.StandardTISAnalysis(\n", + " network=network,\n", + " flux_method=fake_flux,\n", + " max_lambda_calcs={t: {'bin_width': 0.025, 'bin_range': (-0.6, 0.6)}\n", + " for t in network.sampling_transitions},\n", + " combiners={t.interfaces: paths.numerics.WHAM(cutoff=0.01, # c cutoff, default is 0.05\n", + " interfaces=t.interfaces.lambdas)\n", + " for t in network.sampling_transitions}\n", + ")\n", + "parallel_analysis.results['flux'] = fake_flux.calculate('foo')\n", + "parallel_analysis.results = parallel_analysis.from_weighted_trajectories(weighted_trajectories)" + ] + }, + { + "cell_type": "markdown", + "id": "423f6e36-a97f-4fe3-9509-49a215ea5fd4", + "metadata": {}, + "source": [ + "### Plotting the crossing probabilities\n", + "\n", + "One of the spot-checks to see if your simulation is converged is to plot the crossing probabilities functions. For each ensemble, the `StandardTISAnalysis` calculates a crossing probability along the order parameter, defined as the fraction of paths in that ensemble that reach at least the given value on the $x$ axis. As such, the crossing probability is always 1 for values less than the cutoff for the interface. Additionally, two ensemble crossing probabilities should never cross; the one from an outer interface should always be higher at a given value of the order parameter than one from an inner interface.\n", + "\n", + "There is also the *total* crossing probability, which is generated by using a histogram combining algorithm (usually WHAM) to combine the individual ensemble crossing probabilities into a good estimate for the true crossing probability (from the innermost interface). Like all crossing probabilities, this should be monotonically decreasing; if it is not, that is a sign of insufficient sampling.\n", + "\n", + "Since the y-axis is probability, and we're looking at rare events, we frequently plot crossing probabilities on a semi-log plot." + ] + }, + { + "cell_type": "markdown", + "id": "eecd58a1", + "metadata": {}, + "source": [ + "Note: In some cases, you may get an error here about a lack of overlap between ensembles. Normally, this can be solved by continuing the simulation longer. For the purposes of this tutorial, you can probably just re-run the sampling stage and then re-run this notebook." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "372e1eb0-d3e2-49ed-9ad0-2ed5405a8a23", + "metadata": {}, + "outputs": [], + "source": [ + "for ensemble in network.transitions[(state_A, state_B)].ensembles:\n", + " crossing = retis_analysis.crossing_probability(ensemble)\n", + " label = \"Interface at $x$={:3.2f}\".format(ensemble.lambda_i)\n", + " plt.plot(crossing.x, crossing, label=label)\n", + "\n", + "tcp_AB = retis_analysis.total_crossing_probability[(state_A, state_B)]\n", + "plt.plot(tcp_AB.x, tcp_AB, lw=2, color='k', label=\"Total crossing probability\")\n", + "plt.legend(loc='upper right')\n", + "plt.yscale('log')\n", + "plt.xlabel('$x$')\n", + "plt.ylabel('Crossing probability')\n", + "plt.title(\"RETIS\");" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "98fc8567-14d5-4f03-9391-9f598eef0d7c", + "metadata": {}, + "outputs": [], + "source": [ + "for ensemble in network.transitions[(state_A, state_B)].ensembles:\n", + " crossing = parallel_analysis.crossing_probability(ensemble)\n", + " label = \"Interface at x={:3.2f}\".format(ensemble.lambda_i)\n", + " plt.plot(crossing.x, crossing, label=label)\n", + "\n", + "tcp_AB = parallel_analysis.total_crossing_probability[(state_A, state_B)]\n", + "plt.plot(tcp_AB.x, tcp_AB, lw=2, color='k', label=\"Total crossing probability\")\n", + "plt.legend(loc='upper right')\n", + "plt.yscale('log')\n", + "plt.xlabel('$x$')\n", + "plt.ylabel('Crossing probability')\n", + "plt.title('Parallel TIS');" + ] + }, + { + "cell_type": "markdown", + "id": "c4c66d93-fd79-4745-9ffb-45c97a01859a", + "metadata": {}, + "source": [ + "* Are your total crossing probabilities monotonically decreasing?\n", + "* Do your individual ensemble crossing probabilities cross each other?\n", + "* Based on the crossing probability plots, which approach seems more converged?" + ] + }, + { + "cell_type": "markdown", + "id": "d233219d-2d88-4f69-b98e-d5a9a3dc4545", + "metadata": {}, + "source": [ + "OPS keeps track of how long each sampling step took. We can look at that to get a rough comparison of the sampokling times (note that other aspects of the computing environment, such as other processes running at the same time, may have a significant effect on the timings here)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5721a2d6-3e8f-448c-8586-35d430e702e5", + "metadata": {}, + "outputs": [], + "source": [ + "def calculate_total_time(storage):\n", + " # step 0 (init conds) doesn't have timing data; all the others do\n", + " return sum(step.change.details.timing for step in storage.steps[1:])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "82efcac2-253b-4e53-be84-e41e803643e5", + "metadata": {}, + "outputs": [], + "source": [ + "%%time\n", + "calculate_total_time(storage_retis)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b7db6cf9-2a56-49ae-b29d-02ef7762a6e2", + "metadata": {}, + "outputs": [], + "source": [ + "total_times = [calculate_total_time(storage)\n", + " for storage in [storage_0, storage_1, storage_2, storage_3, storage_4]]\n", + "print(sum(total_times), total_times)" + ] + }, + { + "cell_type": "markdown", + "id": "0d1ea868-295d-410d-bdb4-40ae6f518aa5", + "metadata": {}, + "source": [ + "* The times reported here do not include the cost of storing results to disk. How do they compare to the time it actually took to run the simulation? Is storing to disk a significant overhead for simulations of toy models?\n", + "* The RETIS simulation includes replica exchange moves and path reversals, as well as the shooting moves used by both approaches. Do replica exchange and path reversal contribute significantly to the total simulation time?" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.10" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/binder/environment.yml b/binder/environment.yml index a9c949d..39a1464 100644 --- a/binder/environment.yml +++ b/binder/environment.yml @@ -9,6 +9,7 @@ dependencies: - openmm - tqdm - openpathsampling + - openpathsampling-cli - openmmtools - jupyterlab prefix: /Users/dwhs/miniconda3/envs/ops-tutorial diff --git a/devtools/a1.patch b/devtools/a1.patch index e009af8..7a15e84 100644 --- a/devtools/a1.patch +++ b/devtools/a1.patch @@ -1,8 +1,8 @@ diff --git a/1_tps_sampling_tutorial.ipynb b/1_tps_sampling_tutorial.ipynb -index 7814be1..555e20a 100644 +index bf34b4e..e7d75d3 100644 --- a/1_tps_sampling_tutorial.ipynb +++ b/1_tps_sampling_tutorial.ipynb -@@ -191,6 +191,16 @@ +@@ -187,6 +187,16 @@ "# Figure out which atoms make the dihedral, find their atom indices, then create the CV." ] }, @@ -19,7 +19,7 @@ index 7814be1..555e20a 100644 { "cell_type": "markdown", "metadata": {}, -@@ -232,6 +242,21 @@ +@@ -228,6 +238,21 @@ "# YOUR TURN: define the `alpha_R` state" ] }, diff --git a/devtools/a2.patch b/devtools/a2.patch index 1dc49b6..8623c0a 100644 --- a/devtools/a2.patch +++ b/devtools/a2.patch @@ -1,8 +1,8 @@ diff --git a/2_tps_analysis_tutorial.ipynb b/2_tps_analysis_tutorial.ipynb -index 6d25bc4..b335e13 100644 +index 5bd5f82..2a97107 100644 --- a/2_tps_analysis_tutorial.ipynb +++ b/2_tps_analysis_tutorial.ipynb -@@ -122,6 +122,16 @@ +@@ -126,6 +126,16 @@ "#traj = storage.steps[####].active[0].trajectory" ] }, @@ -19,7 +19,7 @@ index 6d25bc4..b335e13 100644 { "cell_type": "code", "execution_count": null, -@@ -165,6 +175,44 @@ +@@ -176,6 +186,44 @@ "# Did your selected snapshot come from an accepted or rejected trial?" ] }, @@ -64,7 +64,7 @@ index 6d25bc4..b335e13 100644 { "cell_type": "code", "execution_count": null, -@@ -422,6 +470,17 @@ +@@ -437,6 +485,17 @@ "# trajB = flexible.steps[####].active[0].trajectory" ] }, diff --git a/devtools/a3.patch b/devtools/a3.patch index 722e642..63e8636 100644 --- a/devtools/a3.patch +++ b/devtools/a3.patch @@ -1,5 +1,5 @@ diff --git a/3_committor_analysis_tutorial.ipynb b/3_committor_analysis_tutorial.ipynb -index 6cbc363..957f0ab 100644 +index ee1f4e3..abdd583 100644 --- a/3_committor_analysis_tutorial.ipynb +++ b/3_committor_analysis_tutorial.ipynb @@ -168,6 +168,22 @@ diff --git a/devtools/a5.patch b/devtools/a5.patch new file mode 100644 index 0000000..b5054e6 --- /dev/null +++ b/devtools/a5.patch @@ -0,0 +1,43 @@ +diff --git a/5_custom_shooting_setup.ipynb b/5_custom_shooting_setup.ipynb +index 362391b..a1df0ba 100644 +--- a/5_custom_shooting_setup.ipynb ++++ b/5_custom_shooting_setup.ipynb +@@ -170,6 +170,16 @@ + "# beta = ... # fill in the ellipsis and uncomment this line" + ] + }, ++ { ++ "cell_type": "code", ++ "execution_count": null, ++ "metadata": {}, ++ "outputs": [], ++ "source": [ ++ "### ANSWER:\n", ++ "beta = 1.0 / engine.integ.temperature # kB = 1" ++ ] ++ }, + { + "cell_type": "code", + "execution_count": null, +@@ -317,6 +327,21 @@ + "# 2. Append that strategy and the global_strategy to the biased_shooting_scheme" + ] + }, ++ { ++ "cell_type": "code", ++ "execution_count": null, ++ "metadata": {}, ++ "outputs": [], ++ "source": [ ++ "### ANSWER:\n", ++ "biased_strategy = strategies.TwoWayShootingStrategy(\n", ++ " modifier=modifier,\n", ++ " selector=gaussian_sel,\n", ++ " engine=engine\n", ++ ")\n", ++ "biased_shooting_scheme.append([global_strategy, biased_strategy])" ++ ] ++ }, + { + "cell_type": "markdown", + "metadata": {}, diff --git a/devtools/a6.patch b/devtools/a6.patch new file mode 100644 index 0000000..fa3c152 --- /dev/null +++ b/devtools/a6.patch @@ -0,0 +1,26 @@ +diff --git a/6_custom_shooting_analysis.ipynb b/6_custom_shooting_analysis.ipynb +index f3b0931..da7622e 100644 +--- a/6_custom_shooting_analysis.ipynb ++++ b/6_custom_shooting_analysis.ipynb +@@ -194,6 +194,21 @@ + "# biased_rej_x, biased_rej_y = ..." + ] + }, ++ { ++ "cell_type": "code", ++ "execution_count": null, ++ "id": "09dfccc9-bb21-47d4-8b9f-c915a1f6e05e", ++ "metadata": {}, ++ "outputs": [], ++ "source": [ ++ "### ANSWER:\n", ++ "std_rej_steps = [step for step in std_two_way.steps[1:] if not step.change.accepted]\n", ++ "std_rej_x, std_rej_y = get_shooting_points(std_rej_steps)\n", ++ "\n", ++ "biased_rej_steps = [step for step in biased_two_way.steps[1:] if not step.change.accepted]\n", ++ "biased_rej_x, biased_rej_y = get_shooting_points(biased_rej_steps)" ++ ] ++ }, + { + "cell_type": "code", + "execution_count": null, diff --git a/devtools/a7.patch b/devtools/a7.patch new file mode 100644 index 0000000..57ed949 --- /dev/null +++ b/devtools/a7.patch @@ -0,0 +1,45 @@ +diff --git a/7_parallel_tis_setup.ipynb b/7_parallel_tis_setup.ipynb +index 72dfb3c..2358afc 100644 +--- a/7_parallel_tis_setup.ipynb ++++ b/7_parallel_tis_setup.ipynb +@@ -144,6 +144,22 @@ + "# 2. Append the global_strategy and the appropriate shooting strategy to each scheme" + ] + }, ++ { ++ "cell_type": "code", ++ "execution_count": null, ++ "id": "3c79c94f-eeae-4e27-8828-0b3ef42bcbab", ++ "metadata": {}, ++ "outputs": [], ++ "source": [ ++ "### ANSWER:\n", ++ "for ens_num, scheme in enumerate([scheme_0, scheme_1, scheme_2, scheme_3, scheme_4]):\n", ++ " strat = strategies.OneWayShootingStrategy(\n", ++ " ensembles=tis_network.sampling_ensembles[ens_num],\n", ++ " engine=engine\n", ++ " )\n", ++ " scheme.append([strat, global_strategy])" ++ ] ++ }, + { + "cell_type": "markdown", + "id": "48d17350-73c8-402b-b3ba-156cb7a611cc", +@@ -273,6 +289,17 @@ + "# YOUR TURN: Answer the question above" + ] + }, ++ { ++ "cell_type": "code", ++ "execution_count": null, ++ "id": "3746cb59-74ee-455c-a0ad-bc0c971b5141", ++ "metadata": {}, ++ "outputs": [], ++ "source": [ ++ "### ANSWER:\n", ++ "retis_scheme.n_steps_for_trials(retis_scheme.movers['shooting'][0], 250)" ++ ] ++ }, + { + "cell_type": "markdown", + "id": "b524c507-e7e6-4f03-af94-9f2585d1b8c3", diff --git a/devtools/adv_tutorial_setup.py b/devtools/adv_tutorial_setup.py new file mode 100644 index 0000000..9ce334b --- /dev/null +++ b/devtools/adv_tutorial_setup.py @@ -0,0 +1,83 @@ +#!/usr/bin/env python + +import openpathsampling as paths +from openpathsampling.engines import toy as toys + +import numpy as np + +pes = ( + toys.OuterWalls(sigma=[1.0, 1.0], x0=[0.0, 0.0]) + + toys.Gaussian(A=-1.0, alpha=[12.0, 5.0], x0=[-0.6, 0.0]) + + toys.Gaussian(A=-1.0, alpha=[12.0, 5.0], x0=[0.6, 0.0]) +) + +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 +} + +engine = toys.Engine(options, topology).named("toy_engine") + +cv = paths.CoordinateFunctionCV("x", lambda snap: snap.xyz[0][0]) +state_A = paths.CVDefinedVolume(cv, float("-inf"), -0.6).named("A") +state_B = paths.CVDefinedVolume(cv, 0.6, float("inf")).named("B") + +snap = toys.Snapshot(velocities=np.array([[0.0, 0.0]]), + coordinates=np.array([[0.0, 0.0]]), + engine=engine) + +randomizer = paths.RandomVelocities(beta=10.0, engine=engine) + +committor_storage = paths.Storage("toy_committor.nc", mode='w') + +committor = paths.CommittorSimulation( + storage=committor_storage, + engine=engine, + states=[state_A, state_B], + randomizer=randomizer, + initial_snapshots=[snap], + direction=1 +) + +committor.run(10) + +committor_storage.close() + + + +committor_storage = paths.Storage("toy_committor.nc", mode='r') +spa = paths.ShootingPointAnalysis(committor_storage.steps, [state_A, state_B]) + +end_A = [step for step in committor_storage.steps if state_A(step.change.canonical.trials[0][-1])] + +end_B = [step for step in committor_storage.steps if state_B(step.change.canonical.trials[0][-1])] + +part_A = end_A[-1].change.canonical.trials[0].trajectory.reversed +part_B = end_B[-1].change.canonical.trials[0].trajectory + +traj = part_A + part_B[1:] + +equil_network = paths.TPSNetwork(state_A, state_B) +scheme = paths.OneWayShootingMoveScheme(network=equil_network, engine=engine) +init_conds = scheme.initial_conditions_from_trajectories(traj) + +sim = paths.PathSampling( + storage=None, + move_scheme=scheme, + sample_set=init_conds +) + +sim.run_until_decorrelated() + +sim.run(100) + +setup_nc = paths.Storage("2_state_toy.nc", mode='w') +setup_nc.tags['initial_conditions'] = sim.sample_set +#setup_nc.save(sim.sample_set[0].trajectory) +setup_nc.save(state_A) +setup_nc.save(state_B) +setup_nc.save(engine) +setup_nc.close() diff --git a/devtools/patch-all b/devtools/patch-all index 770ee20..8dbcda6 100644 --- a/devtools/patch-all +++ b/devtools/patch-all @@ -5,5 +5,8 @@ DIR="$( cd "$( dirname "${BASH_SOURCE[0]}" )" >/dev/null 2>&1 && pwd )" patch < ${DIR}/a1.patch patch < ${DIR}/a2.patch patch < ${DIR}/a3.patch +patch < ${DIR}/a5.patch +patch < ${DIR}/a6.patch +patch < ${DIR}/a7.patch #patch < ${DIR}/testing.patch diff --git a/devtools/run5.sh b/devtools/run5.sh new file mode 100644 index 0000000..97568bb --- /dev/null +++ b/devtools/run5.sh @@ -0,0 +1,6 @@ +#!/bin/bash + +openpathsampling pathsampling shooting_setup.nc -o 2_way.nc \ + --scheme 2_way -n 200 +openpathsampling pathsampling shooting_setup.nc -o biased.nc \ + --scheme biased_shooting -n 200 diff --git a/devtools/run7.sh b/devtools/run7.sh new file mode 100644 index 0000000..2f627a9 --- /dev/null +++ b/devtools/run7.sh @@ -0,0 +1,14 @@ +#!/bin/bash + +for iface in 0 1 2 3 4; do + echo "SAMPLING INTERFACE ${iface}" + openpathsampling equilibrate parallel_setup.nc -o equil_${iface}.nc \ + --scheme scheme_${iface} --extra-steps 5 # user does 50 + openpathsampling pathsampling equil_${iface}.nc -o scheme_${iface}.nc \ + -n 300 +done + +echo "SAAMPLING RETIS" +openpathsampling equilibrate parallel_setup.nc -o equil_retis.nc \ + --scheme retis --extra-steps 5 # user does 50 +openpathsampling pathsampling equil_retis.nc -o retis.nc -n 2850 diff --git a/inputs/2_state_toy_38.nc b/inputs/2_state_toy_38.nc new file mode 100644 index 0000000..e1d06e9 Binary files /dev/null and b/inputs/2_state_toy_38.nc differ