{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Analyzing the flexible path length simulation"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"from __future__ import print_function\n",
"%matplotlib inline\n",
"import openpathsampling as paths\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import os\n",
"import openpathsampling.visualize as ops_vis\n",
"from IPython.display import SVG"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Load the file, and from the file pull our the engine (which tells us what the timestep was) and the move scheme (which gives us a starting point for much of the analysis)."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# note that this log will overwrite the log from the previous notebook\n",
"#import logging.config\n",
"#logging.config.fileConfig(\"logging.conf\", disable_existing_loggers=False)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"CPU times: user 1min 2s, sys: 3.52 s, total: 1min 6s\n",
"Wall time: 1min 6s\n"
]
}
],
"source": [
"%%time\n",
"flexible = paths.AnalysisStorage(\"ad_tps.nc\")\n",
"# opening as AnalysisStorage is a little slower, but speeds up the move_summary"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"engine = flexible.engines[0]\n",
"flex_scheme = flexible.schemes[0]"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"File size: 18.65GB for 10001 steps, 985686 snapshots\n"
]
}
],
"source": [
"print(\"File size: {0} for {1} steps, {2} snapshots\".format(\n",
" flexible.file_size_str,\n",
" len(flexible.steps),\n",
" len(flexible.snapshots)\n",
"))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"That tell us a little about the file we're dealing with. Now we'll start analyzing the contents of that file. We used a very simple move scheme (only shooting), so the main information that the `move_summary` gives us is the acceptance of the only kind of move in that scheme. See the MSTIS examples for more complicated move schemes, where you want to make sure that frequency at which the move runs is close to what was expected."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "6d9ea6ac33644ce4a0bd7cd57a00a0c1",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"HBox(children=(HTML(value=''), FloatProgress(value=0.0, max=10001.0), HTML(value='')))"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"shooting ran 100.000% (expected 100.00%) of the cycles with acceptance 5529/10000 (55.29%)\n"
]
}
],
"source": [
"flex_scheme.move_summary(flexible.steps)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Replica history tree and decorrelated trajectories\n",
"\n",
"The `ReplicaHistoryTree` object gives us both the history tree (often called the \"move tree\") and the number of decorrelated trajectories.\n",
"\n",
"A `ReplicaHistoryTree` is made for a certain set of Monte Carlo steps. First, we make a tree of only the first 25 steps in order to visualize it. (All 10000 steps would be unwieldy.) \n",
"\n",
"After the visualization, we make a second `PathTree` of all the steps. We won't visualize that; instead we use it to count the number of decorrelated trajectories."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"replica_history = ops_vis.ReplicaEvolution(replica=0)\n",
"tree = ops_vis.PathTree(\n",
" flexible.steps[0:25],\n",
" replica_history\n",
")\n",
"tree.options.css['scale_x'] = 3\n",
"\n",
"SVG(tree.svg())"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"# can write to svg file and open with programs that can read SVG\n",
"with open(\"flex_tps_tree.svg\", 'w') as f:\n",
" f.write(tree.svg())"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Decorrelated trajectories: 3\n"
]
}
],
"source": [
"print(\"Decorrelated trajectories:\", len(tree.generator.decorrelated_trajectories))"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"All decorrelated trajectories: 846\n",
"CPU times: user 1min 22s, sys: 321 ms, total: 1min 22s\n",
"Wall time: 1min 22s\n"
]
}
],
"source": [
"%%time\n",
"full_history = ops_vis.PathTree(\n",
" flexible.steps,\n",
" ops_vis.ReplicaEvolution(\n",
" replica=0\n",
" )\n",
")\n",
"\n",
"n_decorrelated = len(full_history.generator.decorrelated_trajectories)\n",
"\n",
"print(\"All decorrelated trajectories:\", n_decorrelated)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Path length distribution\n",
"\n",
"Flexible length TPS gives a distribution of path lengths. Here we calculate the length of every accepted trajectory, then histogram those lengths, and calculate the maximum and average path lengths.\n",
"\n",
"We also use `engine.snapshot_timestep` to convert the count of frames to time, including correct units."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Maximum: 449 (8.980 ps)\n",
"Average: 84.56 (1.691 ps)\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAQnUlEQVR4nO3df6jdd33H8edraY3VTmzX25IlYYmQbkuLm3rJujlEVkszFdN/Chlzhq0QJt2mdeCSCSv7I9D9wKqwCqF2RuwagjoaBDdDVGSgjbe2atOYNpqtuWvWXCfOukE09b0/zrd6vDm5yT3n3nNz83k+4HC+3/f38z3fz/m0vM433/M9n5uqQpLUhp9b6g5IksbH0Jekhhj6ktQQQ1+SGmLoS1JDLlvqDpzPNddcU+vWrVvqbkjSsvLoo49+p6omZtcv+tBft24dU1NTS90NSVpWkvzHoPp5L+8keSDJqSRP9NX+Lsk3k3w9yT8neWXftp1JjiU5muTWvvrrknyj2/ahJBnxPUmS5ulCrul/FNg8q3YAuLGqXg08BewESLIR2Arc0O1zX5IV3T4fBrYDG7rH7NeUJC2y84Z+VX0R+O6s2mer6ky3+mVgTbe8BdhbVaer6jhwDNiUZBXwiqr6UvV+Avwx4LYFeg+SpAu0EHfv/BHwmW55NXCib9t0V1vdLc+uS5LGaKTQT/I+4Azw4IulAc1qjvq5Xnd7kqkkUzMzM6N0UZLUZ+jQT7INeCvw+/XTWdumgbV9zdYAz3b1NQPqA1XV7qqarKrJiYmz7jiSJA1pqNBPshn4C+BtVfV/fZv2A1uTrEyynt4Xtoeq6iTwfJKburt23gE8PGLfJUnzdN779JM8BLwRuCbJNHA3vbt1VgIHujsvv1xVf1xVh5PsA56kd9nnzqp6oXupd9K7E+gKet8BfAZJ0ljlYp9Pf3JysvxxliTNT5JHq2pydv2i/0XuxereA0+dc9tdt1w/xp5I0oUz9M9hrlCXpOXKWTYlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQ/zD6IjjfH1W/65brx9QTSfpZnulLUkMMfUlqiKEvSQ0x9CWpIecN/SQPJDmV5Im+2tVJDiR5unu+qm/bziTHkhxNcmtf/XVJvtFt+1CSLPzbkSTN5ULO9D8KbJ5V2wEcrKoNwMFunSQbga3ADd0+9yVZ0e3zYWA7sKF7zH5NSdIiO2/oV9UXge/OKm8B9nTLe4Db+up7q+p0VR0HjgGbkqwCXlFVX6qqAj7Wt48kaUyGvaZ/XVWdBOier+3qq4ETfe2mu9rqbnl2faAk25NMJZmamZkZsouSpNkW+sdZg67T1xz1gapqN7AbYHJy8pztlit/vCVpqQx7pv9cd8mG7vlUV58G1va1WwM829XXDKhLksZo2NDfD2zrlrcBD/fVtyZZmWQ9vS9sD3WXgJ5PclN31847+vaRJI3JeS/vJHkIeCNwTZJp4G7gHmBfkjuAZ4DbAarqcJJ9wJPAGeDOqnqhe6l30rsT6ArgM91DkjRG5w39qvq9c2y6+RztdwG7BtSngBvn1TtJ0oLyF7mS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDblsqTuwlO498NRSd0GSxsozfUlqiKEvSQ0x9CWpIYa+JDVkpNBPcleSw0meSPJQkpcmuTrJgSRPd89X9bXfmeRYkqNJbh29+5Kk+Rg69JOsBv4MmKyqG4EVwFZgB3CwqjYAB7t1kmzstt8AbAbuS7JitO5LkuZj1Ms7lwFXJLkMeBnwLLAF2NNt3wPc1i1vAfZW1emqOg4cAzaNeHxJ0jwMHfpV9Z/A3wPPACeB/6mqzwLXVdXJrs1J4Npul9XAib6XmO5qZ0myPclUkqmZmZlhuyhJmmWUyztX0Tt7Xw/8IvDyJG+fa5cBtRrUsKp2V9VkVU1OTEwM20VJ0iyjXN55E3C8qmaq6kfAp4DfAp5Lsgqgez7VtZ8G1vbtv4be5SBJ0piMEvrPADcleVmSADcDR4D9wLauzTbg4W55P7A1ycok64ENwKERji9Jmqeh596pqkeSfAL4KnAGeAzYDVwJ7EtyB70Phtu79oeT7AOe7NrfWVUvjNh/SdI8jDThWlXdDdw9q3ya3ln/oPa7gF2jHFOSNDx/kStJDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqyEg/ztLiuPfAU+fcdtct14+xJ5IuNZ7pS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIZ4y+YyM9ftnOAtnZLm5pm+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDvHvnEuNkbZLm4pm+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDRgr9JK9M8okk30xyJMlvJrk6yYEkT3fPV/W135nkWJKjSW4dvfuSpPkY9Uz/g8C/VNWvAL8GHAF2AAeragNwsFsnyUZgK3ADsBm4L8mKEY8vSZqHoUM/ySuANwAfAaiqH1bV94AtwJ6u2R7gtm55C7C3qk5X1XHgGLBp2ONLkuZvlDP9VwEzwD8meSzJ/UleDlxXVScBuudru/argRN9+093tbMk2Z5kKsnUzMzMCF2UJPUbJfQvA14LfLiqXgP8L92lnHPIgFoNalhVu6tqsqomJyYmRuiiJKnfKKE/DUxX1SPd+ifofQg8l2QVQPd8qq/92r791wDPjnB8SdI8DR36VfVfwIkkv9yVbgaeBPYD27raNuDhbnk/sDXJyiTrgQ3AoWGPL0mav1EnXPtT4MEkLwG+DfwhvQ+SfUnuAJ4BbgeoqsNJ9tH7YDgD3FlVL4x4fEnSPIwU+lX1ODA5YNPN52i/C9g1yjElScPzF7mS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUkFH/XOJF7d4DTy11FyTpouKZviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0JakhI4d+khVJHkvy6W796iQHkjzdPV/V13ZnkmNJjia5ddRjS5LmZyHO9N8FHOlb3wEcrKoNwMFunSQbga3ADcBm4L4kKxbg+JKkCzRS6CdZA7wFuL+vvAXY0y3vAW7rq++tqtNVdRw4Bmwa5fiSpPkZ9Uz/A8B7gR/31a6rqpMA3fO1XX01cKKv3XRXO0uS7UmmkkzNzMyM2EVJ0ouGDv0kbwVOVdWjF7rLgFoNalhVu6tqsqomJyYmhu2iJGmWUWbZfD3wtiRvBl4KvCLJx4HnkqyqqpNJVgGnuvbTwNq+/dcAz45wfEnSPA19pl9VO6tqTVWto/cF7eeq6u3AfmBb12wb8HC3vB/YmmRlkvXABuDQ0D2XJM3bYsynfw+wL8kdwDPA7QBVdTjJPuBJ4AxwZ1W9sAjHlySdw4KEflV9AfhCt/zfwM3naLcL2LUQx5QkzZ+/yJWkhhj6ktQQQ1+SGmLoS1JDDH1Jashi3LKpi9S9B56ac/tdt1w/pp5IWiqGvn7CDwXp0mfo64LN9aHgB4K0PHhNX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNccI1LQhn6JSWB8/0Jakhhr4kNcTQl6SGGPqS1BBDX5IaMnToJ1mb5PNJjiQ5nORdXf3qJAeSPN09X9W3z84kx5IcTXLrQrwBSdKFG+VM/wzw51X1q8BNwJ1JNgI7gINVtQE42K3TbdsK3ABsBu5LsmKUzkuS5mfo0K+qk1X11W75eeAIsBrYAuzpmu0BbuuWtwB7q+p0VR0HjgGbhj2+JGn+FuSafpJ1wGuAR4Drquok9D4YgGu7ZquBE327TXe1Qa+3PclUkqmZmZmF6KIkiQUI/SRXAp8E3l1V35+r6YBaDWpYVburarKqJicmJkbtoiSpM1LoJ7mcXuA/WFWf6srPJVnVbV8FnOrq08Davt3XAM+OcnxJ0vyMcvdOgI8AR6rq/X2b9gPbuuVtwMN99a1JViZZD2wADg17fEnS/I0y4drrgT8AvpHk8a72l8A9wL4kdwDPALcDVNXhJPuAJ+nd+XNnVb0wwvElSfM0dOhX1b8x+Do9wM3n2GcXsGvYY0qSRuPUyhqLuaZedtplaXwMfV30nKtfWjjOvSNJDTH0Jakhhr4kNcRr+lpy57tmL2nheKYvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDvE9fy56TuUkXzjN9SWqIZ/pqmv9KUGs805ekhhj6ktQQQ1+SGmLoS1JD/CJXOoel/DON/olILRZDX5c05+qXfpaXdySpIYa+JDXEyzvSkLx0pOXIM31Jaohn+pJ+wruGLn1jD/0km4EPAiuA+6vqnnH3QbqUGdyay1hDP8kK4B+AW4Bp4CtJ9lfVk+Psh9SyUb6LWKoPFD/IFs64z/Q3Aceq6tsASfYCWwBDX5qHi/VL5FFmLV3MD6O5jPKBsZj/HRbrg2zcob8aONG3Pg38xuxGSbYD27vVHyQ5OuC1rgG+s+A9XP4cl7M5JoONdVzeM64DzdOAfl0U/78swHj90qDiuEM/A2p1VqFqN7B7zhdKpqpqcqE6dqlwXM7mmAzmuAx2qY/LuG/ZnAbW9q2vAZ4dcx8kqVnjDv2vABuSrE/yEmArsH/MfZCkZo318k5VnUnyJ8C/0rtl84GqOjzky815+adhjsvZHJPBHJfBLulxSdVZl9QlSZcop2GQpIYY+pLUkGUX+kk2Jzma5FiSHUvdn3FK8kCSU0me6KtdneRAkqe756v6tu3sxulokluXpteLL8naJJ9PciTJ4STv6urNjk2SlyY5lORr3Zj8dVdvdkz6JVmR5LEkn+7W2xmXqlo2D3pf/n4LeBXwEuBrwMal7tcY3/8bgNcCT/TV/hbY0S3vAP6mW97Yjc9KYH03biuW+j0s0risAl7bLf888FT3/psdG3q/ibmyW74ceAS4qeUxmTU+7wH+Cfh0t97MuCy3M/2fTONQVT8EXpzGoQlV9UXgu7PKW4A93fIe4La++t6qOl1Vx4Fj9MbvklNVJ6vqq93y88ARer/+bnZsqucH3erl3aNoeExelGQN8Bbg/r5yM+Oy3EJ/0DQOq5eoLxeL66rqJPTCD7i2qzc5VknWAa+hd2bb9Nh0lzAeB04BB6qq+THpfAB4L/Djvloz47LcQv+CpnEQ0OBYJbkS+CTw7qr6/lxNB9QuubGpqheq6tfp/fJ9U5Ib52jexJgkeStwqqoevdBdBtSW9bgst9B3GoezPZdkFUD3fKqrNzVWSS6nF/gPVtWnurJjA1TV94AvAJtxTF4PvC3Jv9O7PPw7ST5OQ+Oy3ELfaRzOth/Y1i1vAx7uq29NsjLJemADcGgJ+rfokgT4CHCkqt7ft6nZsUkykeSV3fIVwJuAb9LwmABU1c6qWlNV6+jlx+eq6u20NC5L/U3yfB/Am+ndnfEt4H1L3Z8xv/eHgJPAj+idgdwB/AJwEHi6e766r/37unE6CvzuUvd/Ecflt+n9k/vrwOPd480tjw3wauCxbkyeAP6qqzc7JgPG6I389O6dZsbFaRgkqSHL7fKOJGkEhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqyP8DJedGOQpBqqgAAAAASUVORK5CYII=\n",
"text/plain": [
"