Merge pull request #7 from dwhswenson/molsim2021-updates

Updates based on feedback from MolSim 2021

Author: dwhswenson

2_tps_analysis_tutorial.ipynb

@@ -98,9 +98,13 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
+    "From looking at the path tree, we can get a sense of how well our sampling is going. Some questions to consider:\n",
+    "\n",
     "* Do you see acceptance of both forward and backward shots?\n",
     "* How many decorrelated paths do you see? \n",
-    "* How many accepted moves between decorrelated paths? How many total moves?"
+    "* How many accepted moves between decorrelated paths? How many total moves?\n",
+    "\n",
+    "If you don't see acceptance of both forward and backward shots, that might indicate a problem in either the state definitions or in the sampling (which might be remedied by using a more advanced shooting point selection algorithm). The number MC steps between decorrelated paths gives a sense of how quickly your simulation is exploring path space. If you see very few decorrelated paths (or none at all), that may indicate that you should switch to a more advanced algorithm such as spring shooting."
    ]
   },
   {
@@ -152,6 +156,13 @@
     "```"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Next let's visualize the shooting point from one of your MC steps. There are two reasons we do this: First, this gives you more practice visualizing OPS data with NGLView. Second, this shows you how to extract the shooting point snapshot from a MC step. Analyzing the shooting points can give ideas on how to troubleshoot problems with sampling (e.g., if you had very few accepted backward shots)."
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -319,14 +330,16 @@
     "path_lengths = [len(step.active[0].trajectory) for step in flexible.steps]\n",
     "plt.hist(path_lengths, bins=40, alpha=0.5);\n",
     "print(\"Maximum:\", max(path_lengths), \"(\"+str(max(path_lengths)*engine.snapshot_timestep)+\")\")\n",
-    "print(\"Average:\", \"{0:.2f}\".format(np.mean(path_lengths)), \"(\"+(np.mean(path_lengths)*engine.snapshot_timestep).format(\"%.3f\")+\")\")"
+    "print(\"Average:\", \"{0:.2f}\".format(np.mean(path_lengths)), \"(\"+(np.mean(path_lengths)*engine.snapshot_timestep).format(\"%.3f\")+\")\")\n",
+    "plt.ylabel(\"Count\")\n",
+    "plt.xlabel(\"Path length (Frames)\");"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "* Estimate (by eye) what the path length should be if you ensured that fixed-length TPS would include 99% of the paths. How much longer is that than the average path length here?"
+    "* Estimate (by eye) what the path length should be if you ensured that fixed-length TPS would include 99% of the paths. How much longer is that than the average path length here? What does that mean about the simulation time required for fixed-length TPS as compared the flexible-length TPS?"
    ]
   },
   {
@@ -395,7 +408,9 @@
     "                             xticklabels=tick_labels,\n",
     "                             yticklabels=tick_labels, \n",
     "                             label_format=\"{:4.2f}\")\n",
-    "ax = plotter.plot(cmap=\"Blues\")"
+    "ax = plotter.plot(cmap=\"Blues\")\n",
+    "plt.xlabel(\"$\\phi$\")\n",
+    "plt.ylabel(\"$\\psi$\");"
    ]
   },
   {
@@ -407,7 +422,7 @@
     "* Which one appears to be more frequent? \n",
     "* Do they look like completely separate channels, or does it look like some paths involves parts of both?\n",
     "\n",
-    "Next, you'll find a representative trajectory for each channel. This will just involve trial and error until you find one."
+    "Next, you'll find a representative trajectory for each channel. This will just involve trial and error until you find one. You're likely to also find trajectories that at least partly cross into the other channel -- these two channels are not very strictly separated."
    ]
   },
   {
@@ -431,14 +446,16 @@
     "ax = plotter.plot(xlim=(-np.pi, 0), ylim=(-np.pi/2, np.pi), cmap=\"Blues\")\n",
     "plotter.plot_trajectory(trajA, '-r', lw=1)\n",
     "plotter.plot_trajectory(trajB, '-k', lw=1)\n",
+    "plt.xlabel(\"$\\phi$\")\n",
+    "plt.ylabel(\"$\\psi$\");\n",
     "#plt.savefig(\"AD_tps_pathdensity.pdf\")  # uncomment this line to save as PDF"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "If you have the full simulation data, you could visualize these trajectories with NGLView. (Full simulation data not provided online, due to size), using the same approach as above."
+    "If you have the full simulation data, you could visualize these trajectories with NGLView, using the same approach as above. (Full simulation data not provided online, due to size.)"
    ]
   },
   {
@@ -465,7 +482,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.3"
+   "version": "3.7.10"
   },
   "toc": {
    "base_numbering": 1,

3_committor_analysis_tutorial.ipynb

@@ -236,7 +236,9 @@
    "source": [
     "You'll notice that this is *not* very sharply peaked in $\\psi$. This means that $\\psi$ alone probably doesn't define the transition state.\n",
     "\n",
-    "If $\\psi$ did define the transition state, we'd probably expect the transition state to be around $\\psi = 0.6$ (\\~35 degrees) to $\\psi = 1.0$ (\\~60 degrees). Now, instead of keeping the committor probability fixed and histogramming the values of $\\psi$, let's see what the values of the committor are for a fixed value of $\\psi = \\psi^*$. If $\\psi^*$ defined the transition state, then the committor would always be 0.5 here."
+    "If $\\psi$ did define the transition state, we'd probably expect the transition state to be around $\\psi = 0.6$ (\\~35 degrees) to $\\psi = 1.0$ (\\~60 degrees). Now, instead of keeping the committor probability fixed and histogramming the values of $\\psi$, let's see what the values of the committor are for a fixed value of $\\psi = \\psi^*$. If $\\psi^*$ defined the transition state, then the committor would always be 0.5 here.\n",
+    "\n",
+    "Previous, we tried to show that frames with a committor around 0.5 had a certain value of $\\psi$. Now we're trying to see if that value of $\\psi$ always implies that frames have a committor around 0.5. We need to show both directions before we can claim that the transition state is defined by that value of $\\psi$."
    ]
   },
   {
@@ -320,7 +322,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.3"
+   "version": "3.7.10"
   },
   "toc": {
    "base_numbering": 1,