{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Obtaining an equilibrated initial trajectory" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is the first file to run in the alanine dipeptide TPS example. This teaches you how to:\n", "\n", "* set up an engine using OpenMM\n", "* set up states using MDTraj-based collective variables\n", "* obtain a initial trajectory using high temperature MD\n", "* equilibrate by using shooting moves until the first decorrelated trajectory\n", "\n", "We assume at this point that you are familiar with the basic concepts of OPS. If you find this file confusing, we recommend working through the toy model examples." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from __future__ import print_function\n", "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "import openpathsampling as paths\n", "\n", "import openpathsampling.engines.openmm as omm\n", "from simtk.openmm import app\n", "import simtk.openmm as mm\n", "import simtk.unit as unit\n", "from openmmtools.integrators import VVVRIntegrator\n", "\n", "import mdtraj as md\n", "\n", "import numpy as np\n", "\n", "import os\n", "initial_pdb = os.path.join(\"..\", \"resources\", \"AD_initial_frame.pdb\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Setting up high-temperature engine" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we set things up for the OpenMM simulation. We will need a `openmm.System` object and an `openmm.Integrator` object.\n", "\n", "To learn more about OpenMM, read the [OpenMM documentation](http://docs.openmm.org). The code we use here is based on output from the convenient web-based [OpenMM builder](http://builder.openmm.org)." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# this cell is all OpenMM specific\n", "forcefield = app.ForceField('amber96.xml', 'tip3p.xml')\n", "pdb = app.PDBFile(initial_pdb)\n", "system = forcefield.createSystem(\n", " pdb.topology, \n", " nonbondedMethod=app.PME, \n", " nonbondedCutoff=1.0*unit.nanometers,\n", " constraints=app.HBonds, \n", " rigidWater=True,\n", " ewaldErrorTolerance=0.0005\n", ")\n", "hi_T_integrator = VVVRIntegrator(\n", " 500*unit.kelvin, \n", " 1.0/unit.picoseconds, \n", " 2.0*unit.femtoseconds)\n", "hi_T_integrator.setConstraintTolerance(0.00001)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The storage file will need a template snapshot. In addition, the OPS OpenMM-based `Engine` has a few properties and options that are set by these dictionaries." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "template = omm.snapshot_from_pdb(initial_pdb)\n", "openmm_properties = {'OpenCLPrecision': 'mixed'}\n", "engine_options = {\n", " 'n_steps_per_frame': 10,\n", " 'n_frames_max': 20000\n", "}" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "hi_T_engine = omm.Engine(\n", " template.topology, \n", " system, \n", " hi_T_integrator, \n", " openmm_properties=openmm_properties,\n", " options=engine_options\n", ")\n", "hi_T_engine.name = '500K'" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "hi_T_engine.current_snapshot = template\n", "hi_T_engine.minimize()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Defining states\n", "\n", "First we define the CVs using the `md.compute_dihedrals` function. Then we define our states using `PeriodicCVDefinedVolume` (since our CVs are periodic.)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "# define the CVs\n", "psi = paths.MDTrajFunctionCV(\n", " \"psi\", md.compute_dihedrals,\n", " topology=template.topology, indices=[[6,8,14,16]]\n", ").enable_diskcache()\n", "phi = paths.MDTrajFunctionCV(\n", " \"phi\", md.compute_dihedrals,\n", " topology=template.topology, indices=[[4,6,8,14]]\n", ").enable_diskcache()" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "# define the states\n", "deg = 180.0/np.pi\n", "C_7eq = (\n", " paths.PeriodicCVDefinedVolume(phi, lambda_min=-180/deg, lambda_max=0/deg, \n", " period_min=-np.pi, period_max=np.pi) &\n", " paths.PeriodicCVDefinedVolume(psi, lambda_min=100/deg, lambda_max=200/deg,\n", " period_min=-np.pi, period_max=np.pi)\n", ").named(\"C_7eq\")\n", "# similarly, without bothering with the labels:\n", "alpha_R = (\n", " paths.PeriodicCVDefinedVolume(phi, -180/deg, 0/deg, -np.pi, np.pi) &\n", " paths.PeriodicCVDefinedVolume(psi, -100/deg, 0/deg, -np.pi, np.pi)\n", ").named(\"alpha_R\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Getting a first trajectory\n", "\n", "We'll use the `VisitAllStatesEnsemble` to create a high-temperature trajectory that visits all of our states. Here we only have 2 states, but this approach generalizes to multiple states." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Ran 2143 frames. Found states [alpha_R,C_7eq]. Looking for [].\n" ] } ], "source": [ "visit_all = paths.VisitAllStatesEnsemble([C_7eq, alpha_R])\n", "trajectory = hi_T_engine.generate(hi_T_engine.current_snapshot, [visit_all.can_append])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plotting the trajectory" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "# create a network so we can use its ensemble to obtain an initial trajectory\n", "# use all-to-all because initial traj can be A->B or B->A; will be reversed\n", "tmp_network = paths.TPSNetwork.from_states_all_to_all([C_7eq, alpha_R])" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[Trajectory[33]]\n" ] } ], "source": [ "# take the subtrajectory matching the ensemble (only one ensemble, only one subtraj)\n", "subtrajectories = []\n", "for ens in tmp_network.analysis_ensembles:\n", " subtrajectories += ens.split(trajectory)\n", "print(subtrajectories)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA6GElEQVR4nO2df3Qc5XX3v3dmdlemhFeJcCswkVXATqzGYNlmwxZeWbx2FdwGssRp36RunQCNWAJJnDRZ6tImJpzIiekBpRCSNbUT60DbQ18TFwMuxsbrH9klRkYGExsISfgV19SIuDQBaaXd+/4hPcPM7MxqV7ur3ZHv55zn2Ds/n5ldfefOfe69DzEzBEEQBP+i1boDgiAIQnmIkAuCIPgcEXJBEASfI0IuCILgc0TIBUEQfI5Ri5OeeeaZ3NraWotTC4Ig+JaDBw++wcwznctrIuStra3o7++vxakFQRB8CxG97LZcXCuCIAg+R4RcEATB54iQC4Ig+BwRckEQBJ8jQi4IguBzRMgFQRB8zikp5Ol0GuvWrUM6nZ7U+krtU85+9ch0uhZB8BM1iSMvh3Q6jWQyic7OTkQikZK3TafTWLp0KTKZDILBIHp7ezE4OIjOzk4AwPr167Ft2zYwM0KhEHbt2gUABc/pdsyBgQEAwKpVq2znTiaTaGpqwuDgIJqamrB69WrXvkQikZKu1e261Xmsx+vr6wMAtLe32/oIwFyn+pxOp7F+/Xq88MILmDt3LpYvX247XqF7sGvXrrKuQRCEEmDmKW+LFi3iyZBKpXjGjBms6zrPmDGDU6lUSdumUinu6upiTdMYAGuaxrquMxGxYRgcCAQYgNmIiMPhMIdCIdZ1nQOBAIfDYU4kErZz9fT0mMckItZ13TyGruucSCTM/ljPHQgEXPsSDAY5kUhwKBRiIuJQKFTwWhOJBHd1dXEikeBEImE7PxFxIBDgeDzOoVDIdn2qBQIB27Wr8xuGkbetpmmu976np8c8r67r3NPTU9L3JQjTmVQqZf5NlAOAfnbR1LJFGUADgAMAngbwUwC3TLTPZIXcKhaapnFXV5fnjXFuO2/ePFMorcLpJmzFtJUrV5pfTCKRKLhtIBDgWCyWdz5N09gwDNZ1PW/de9/7XtvnWCzmep3Oc6sHQzmNiLirq8tzva7rHIvFbD9MN9G2fgdE5HkNgjAdUX8DiUSiYgZNNYWcAJw+/v8AgJ8AuLjQPuVa5FYr1uvGpFIpDgaDpnC7CVU0Gi1b8HRd5+bm5rxlzu3a2tpY1/W8vicSCe7p6ZmwL1YRTKVSHIvFOBaL8bx588oWbrfW2trquc76ENQ0jePxuNkvp7hb3wKCwaBY5cIpgdWwsb55q7fVyVI1IbcdDDgNwFMAPlxou8kKOTOb7hGrZe12Y5wiYm2GYZiuFuW+cLpVJts0TeOVK1fy7NmzXR8iSviUEHuJnvNBoLZVD6hqiHc5ra2tLc/lxMwci8XM+1Duj1gQ/ILTI6DevOvWIh87NnQAhwD8BsC3J9q+HCFnzncnKIvQivVGOpt1+0QiweFwuKAFWk5zE/NoNGp71UokEhyLxTgajXI0GvV8qASDQV6wYEHNRbtQ6+josP1QJ/KTV8p3KAj1hPN3r96869ZHbjsY0AhgN4APuazrBtAPoL+lpaWsi4nFYjbx0HXdVSDUjXQKejgcNq3bQn7yOXPmcENDQ8XFrrm52fa0tvYhEAhwNBrlpqammovyZJsa4LV+F24/YhkMFaYz1TBSpkTIx86DrwP4SqFtyrXInUJORJ7uFTXY4IzACAQCNt/2VLdAIOD6kKmn5vY2Ucr1TRRVZI0gEreLIEyMl5CXnRBERDOJqHH8/zMALAPwXLnHLcSqVasQCATMz4FAwIwDtxKJRLBmzRrMnz8fuVzOtm5kZATHjx+vZjcLMmvWLFxxxRV5/SqH1tZWRKNRaNq7X6uu62hsbCz5WESkHswFtwmHw67rRkdHsXbtWtx000348Ic/jKuuuspMFFIx5zt37kQul4OmaQgGg67foSAIE1OJhKCzAGwmIh1jmaL3M/NDFTiuJ5FIBHv27MlLYLFiTUTp6+urqGAWSyExHBoaMhOPKkVLSwvefvtt2zEvueQS7N+/v+RjldsvZsaOHTuwY8cOc9lDDz2EvXv3IplMIpPJmCK+bNkyrF27VhKGhGlNNZPjqJJCUiyLFy/mSs8QZL1Jhw8fxo033ohsNotQKISPfOQj2Lp1a0XPVwwzZ87EG2+84SqKxVi8VhobG3Hy5MmC2+i6jmw2W9Z5qk0sFgMAbNq0CdlsVrJAhVMCr8znUiGig8y82Lncdyn6blhvkmEYGB0dNQVteHgYx44dg6ZpU26VnzhxAoC7mHqJq2EYWLhwIQ4cOGBbfvbZZ08o5E4RL3SeWrFx40Zks1kQEa644grE43FTxCvxQxeEekS9hWazWWQyGSSTyYr+vn1XNMutMFMymcTw8DCy2az5ryKXy6G/v9/mN55qShHTbDaLa6+9FqFQyFwWCATw0Y9+tKbXYEXX9UnvOzIyglwuh2w2i4ceetcD5/wOk8lkBXoqCPVBZ2cngsEgdF2vyniQryxyN6sNAA4cOJBnbRMRgDERVetmz56Nl192nbu0bmBmbN++Hbt378b69etx7NgxdHZ24s4775zUG0U13kTcLP/JkMvlTMG2foe5XA5NTU0VOYcg1AORSAS7du2qmuvQV0LufD3p6+vD5s2bMTQ0ZNuOiHDRRRfhqaeewujoKIAxcXj11Vdd/cj1Rn9/P7Zu3YqHH34Yo6Oj6O/vn7SLpN5cKwoigqZpuPfee/G1r33N9p1omobBwcEa9k4QKk8kEqmau7A+3tWLpKmpyRSAYDAIAMhkMnliFQgEsHDhQttyIkIul/O0TstxF0yGcDhsvjU4ee2117B+/XqMjIyYbxTM7Lm9lfPPPx9tbW3mZ+e9aWpqwsyZMzFr1qzyLqBMiAijo6M4cuQIRkdHzX4SEXRdF4tcEErANxZ5Op3G6tWrTSFubW3FGWecgWAwaA5yLl++HM3NzWZ97c2bNyOTyUDXdVM4lKtBiSMwNsA4a9asKXO7NDY24uDBgyVby/PmzcPcuXNx7NgxPPnkk3n7ExF++ctfur5xqAdgvVi6bg9UwzDAzMhms1i9ejXmz58vA57CtKDaEVm+EXJr7DEAHD16FEePHsXKlStx4sQJrFixAt3d3bZ9rD6pw4cPY+PGjXjqqadM61YJYTabxauvvjpl1zJR9IkXR44cwXPPPYdAIIBAIGBa7Eqk1SCiGxdddFFeJEyt0DQNuq5jZGTEXBaNRtHc3Ix77rkHuVzOdJ1JOKLgJ4qZzKYqEVlu6Z7VbpNJ0Vd1OeCSRl6onC0ze06S4IemaVpeqryu6xyNRm0TOXR0dBRMta+nYltExPF4nKPRqG2iDmvpYcMwzAk9pA6L4Ae8age5TboyWVCtFP2pQo36RqNR23Ie9yF7haxt2LAB119/vTnoWc8o/7DTF+4MO9Q0DS+88IJpfWezWTNm3Q1mxqFDhyre38nCzLjtttvw4IMP4vDhwwCAdevW4fDhw+ab0ujoqBmOqOJuBaGecYsVB6ofegjAPxa5IpVKuZZ51TTNZtl5Fcuqp6ZpGnd0dHAoFDKnflOWqpqEYsaMGRyNRm1WuXWmI9Wi0ahpzZZT7KoW90DVavYq3ztRAS7rdy6Wu1ArClXzrPup3ibTyhFyrzrjRGTW/bXOzFGPonbWWWfxvHnzzMkY1ANHibkSaquwO6/VKXSxWIzj8XjdXrNXsz54vPodDocL/iakHK5QL1TboJg2Qm79ow2FQhwOh22lULu6umzzRLpZrwC4sbGx5iKmWkdHh+dsQm4PLU3TWNM0DoVCpiVunUbOj01Z5m7rnG9axUz8LAjTES8h903UikL5ylXlw/b2dhw+fNgMQRwaGhp7QuHdh5QbJ0+eNOuyWLFmhE6ErutmFI11+3POOQdvvPFGXqKSF3v37vVc5xaFksvlQET4xCc+gT/4gz/Atm3b6j7JaSKWLVuGBQsW4I477sDo6CiICO9///vR3t6O+fPnFxz5Vz5ItU7K4QpTTc0Lvrmpe7VbuRNLWCdhVu6HWCxW1/7wUpthGEVFmoTD4bqenKKYpmYUUj5+Xdc5Ho/b3CWxWKyg1S0+cqFWTKVrD36PWrGiCizlcjmMjIzg9ttvB1C5GiC1JhgMYu7cufjtb3874bYHDhyoyHVHo1F0dHQUlT1aaS655BJs377dzNLNZrN5EQAACo78q0lEJN5cmGq8olWmEt+5VoCxV2mrWyObzeL48ePQNG1aiHkmk8GRI0fKOoamaTbX0kR1yefOnYvzzjsP+/btK+u8k2Hv3r15IZbDw8PmNRjG2M+0t7cXAwMDU94/QSiEl2tvSt0tbmZ6tVu5rhVm5ng8nueKUIODhZJjpkObM2dOUdupQVEVGVIomsU6aDzZVo3BVmskTygUMgd2JTpFqCecrr1quVswnVwrAPLmoRwdHTUt9DfeeKMGPZo6vvKVryAej09Yn9xaU8b8wjXNdT9r+YPJUo2JO3K5nDlRSCaTwcjIiCQJCXWH07U31e4W3wq5V3U8TdPw9ttvT3FvqovTbz04OIhvf/vb2L9/P2KxmG0i6olwqwC5YMEC/OxnP6tIXyeLpmme16FqswSDQQQCgepmyAlCBZiSbE4LvvSRA/Cs4jcyMjKlBbCqxbx58/Dcc8+5hlCqolvq6b9p06ay5uac6vR9IsKZZ56JwcFB8y3h7rvvBgB87nOfs41zhEIhfPGLX8ShQ4ewYsUKzJ8/XwppCXVPtSeScOJbIe/s7PQUr2w2i+bmZhw/frwGPSsfTdOwZMkSvPjii7YKgYr77rsPwJgAn3baachms2Dmupto2Q3VR2ttGMMwTIFWaJqGZcuWYcWKFfjCF76ATCaDPXv2YPfu3VizZk0Nei4IpVHNiSSc+Na1AhSeDOLiiy+e8skiKkUgEMCqVatw1113mREbVn71q19h/fr12LFjB7Zu3Qrg3YJbgUCgJiGExaBpGj74wQ/mLR8ZGTEtF8MwQEQIBAJYu3YtBgYGMDw8DGbG8PCwmQgmCMK7+FbI+/r6bFmZCxYsMIU7GAwiHo9j0aJFtepeWVx99dWIRCLo7u7G3r17EQ6HC26vfN6lZKVWC7cHj5WZM2e6LldjHqrv6l/nW5Vf37IEoZqULeRE9H4i2k1ER4nop0T0xUp0bCKcf9DPPPMMmBmBQAB33nknAKChoWEqulI0mqbhfe9734Tb7dmzBxs2bAAw9nq2cOHCCfdhHiv96iw5oAYKq2Glux2zULngXC6Hffv25b0pMTNWr16Nvr4+002kkoKam5tt2zo/C4JQGR/5KIC/ZuaniOg9AA4S0WPMXF5GywQ4/6Ctg4K33norXn/99Tz/8plnnonTTz8dL730UjW75gkz49e//vWE2x09ehTXXXcdAKC7uxurVq3Cxo0bXf3lXV1d2LdvnzmlXTabtQ0WEhHuvvtuDAwM4J577pnyhClN08z5UtV35IyaYWYMDQ3hyJEj5mxHuq7jlVdeQXt7O4LBIEZGRkyXk9csLMlkEk1NTRgcHJTBUOHUwi24vJwG4N8B/FGhbcqtfqhqjYdCIbPcq/o/6iBhp9TW2trqWou7q6vLvO5YLJa3nojMe9HV1WWWxLUm5hARx2IxZh6bKanYe6QSiFQtm3A4nLdNMTVeotEop1IpjsViHAwGJ0waUrVWAoGArZaOSrZwS7RIJBK28r0TzRglCLWi3JpAmIoytgBaAbwC4AyXdd0A+gH0t7S0TOoinH/EiUTCvCmliFQ9NE3TeM6cOWb5WcMwuLW11baNKt+qrt1ZFCwYDNrqr1uFzSqyuq6bU6o5a5sXEnH1/2g0yrFYzHZ/1RR7zv3a2tpsy4PBoPmjjcfjRX1HzixUwzA8p83yKpbmLKwlRbWEWmD93VUi27PqQg7gdAAHAXx8om0na5EXqjvd09PjGyFX1f2cVriygMPhcJ4VGovFTKtT0zSORqOmJW6tx67uiVN4reeYTJ+d+3kdJ5FI2M6taZr51uBmjXd0dJiVLNX2Kg3f+tBT12X9YwgGgxwOh/P64rTIZeIJoRY4f3cTVfAsBi8hr0gcOREFAGwBcB8zP1CJY7pRqO50Z2cnGhoazFA15ZetVw4dOpQ3MMjj/uOFCxfiO9/5jun35vHBP4WmaWhubsbq1avNKpDKr3zgwAFcf/31aG9vR0NDA9555528c0wG635e8eqapmFwcBBnnHGGuT6Xy+HkyZP4+te/nvd9BINBfOtb3wKAPP/24cOHccMNNyCXyyEUCpnfdSQSQW9vLzZu3IiBgQH09/eb37eu6/jyl7+MxsZGm4/cmS7d19cnSUVC1fGq4FmVuvlu6l5KA0AA+gD0FrtPJXzkblaV1X9ej/NXqrkplVvIzS8eDAYndH+EQiHb013TNA6Hw7bjGYbB8Xi8KoWsvHzjbtdFRNzV1ZXnllG+81K/a2XlWI+naRp3dXV5Hs85q5Rb4S1xvQiVRL1Fh0Ih22+tbn3kAC4d/4N6BsCh8fbHhfapRPVD5uKm/qqnFovFzAdNT08Pd3V15W0TjUZdBxatIhiLxVxf25zbuk06EQgEOBqN8jnnnDPp67jwwgs9++R06agHinX7eDw+6e/c+d2quVoLibj1nru93orrRagkTvef+tuoBF5CXrZrhZn3Y8wqn1LS6TQ6OzvNsLRkMmmb+kuF4jlpbm7G66+/PmkXw+kAIgAem8S+atoy1W9nHzRNw/bt213DDIExl0ZDQwPa29uRTCbR29truiIA5IUXnn322Th8+LDpfiEiEBHi8Ti2bt2K9evX246/YMGCouquqNdExZ//+Z+jpaUFW7duxYYNG8zrIiLT1aFpmukCclauLAWre03XdVxzzTVYtWqVq4sknU7jsssuM19ld+/eDQDYvHmz7fXWrVKduFyEyWL9PQFAS0uL1CP3wmmBqqeesrC8BvvKsdRbAX4G4N8C/LsOi7QYF0YsFnO1nFVra2uzuUu6uro4Ho+zYRhMRGwYBq9cudKsve60HlW0ChFxKBQyX+WcA6JWq9R6X1TES1tbG8+ePdt0BTm3dVrkhSZ+Vu6WSlooxb6euv1G3PYXi1yoJF6/p0q47zAV4YfFtmoIeTQazQtNVJ+VuCmBLDYEzvr5UoD/C+A3AV7m2FbTNJ45c2ZZQm4YRl5svNreKu5WwbRGcygK+ZXVvXHzwbe1tbnup5Ypf3up4w7qwRGLxTgajeb5DKuJl5C7IT5yoZJUy1iYdkKeSqVM0XMb/FPWrPrXehNXrlxp+wM//fTTbVars10N8DDAzwE8x0OEzz///IKCZrWQQ6FQ3vpwOGz6tK0PHevgnNMytsZXq3sy0UCweltxu4ZCP7LJjDuo70b1X71JKIGfTPhVqb8RNehtjWcXhKmmUOh0KUw7IWd2D7ZXQmFNaFGJM1aRKyY5RQP4HwBmgB8F+H3jWYcTWe/A2HRs6iFiHWxTfZ1I+K1Wt7JordEgmqblJQwV+8RXAqf63tHRkfcQdL4OqjecYqxxZbnrus6zZ8+2LVfTtjkTupzfZ6UQS1uoFOX8lpyGp1jkBXD6gwu9Vvf09BQUozMAfnhcxL8DsI6xaItCPm6nhW3tl1v2ZbGWrRI863Upn7b1ekp54hd6CCr/uzNrVKXqF7LMne4fa9N1nVeuXOn6luSWoSoI9UK5rpFKvR16Cblvy9i6EYlEsGLFCtc5KY8fP450Oo1169YhnU57ThUHjFUS2wPgjwBcB+CLAFjT8NZbb+H48eMTzpUJwBbsn0wmMTQ0hGw2i6GhISSTSXR3dxc176ZKsnH2l5ltEzGUOrWUdY5BNZvJsmXLzOiSTCaDjRs3mv0eHh7GHXfcgYMHD8IwDESjUUSjUYTD4bwqiF410bPZLO677z7s3LkTt99+O4aHh81IkS1btkzpHIeCUArlzsGZTCYxOjoK5rEqpRX/fbupe7VbNS1yZVk6fcrO1/pYLOZpOX5y3BL/M4f7pJTkGmsxK2e8uDWOWlnGHR0drsdRRaOs9US8YqfLdSNYrQ6rPxuwR6ZYLX5nDRhrbRa3hCev70MscqGesWqLCkwoBecbeKn7KzCdXStKwJzJHm51OJSIhMNhcxAuGAzyhRdeaG57AOCjAFORou1shmHwggULPNerqoZW4XVW8LP21ZkAUyiLsdh75ZZlpgZCOzo68vrR2tpqirlX2KPaVr06quOpaBU31021feSCUCnU3+hE1TXV794aZmt1jbpFmxXLtBVyr/Rrleru5c9VESFWkdE0jVfOns0M8A0lprbrus6tra2m4BXatqury0ypV9avilZxi8m2iqrTNz7ZexUMBs3r1nWdOzo6bGJbyJJ2e6OIRqNmP9189M4MSxFrwW8UMw5lDSQA7NFq1Qw/9O3kywpnFtUVV1yBY8eOYWBgANu2bTOzGdU8kBdeeCGefPJJ5HI5jIyM4O2338bIyIhZ0Omql1/GmwB+aCnwpPzGVlpbW3H55Zfj+PHjaG5uRnt7O1avXo2hoaGxJ2QBduzY4bmOmXHppZdi7969tvOr6yvkU3ebcMFKX1+f2T810YPCej4rzgJZuVwOt99+O6LRKABg6dKlyGQyMAwDgUAA2WzW9NE7+yPZkoKfKVS0T5FMJm2Z2cqfvmbNGuzatatqxdp8L+TWm2sYBh555BHP9Pd//Md/BAAcOHAAwJgoWeeQbAUQBfAkgDUA3oOxAjI/1jQ0zJ+PQ08/bW67Zs0adHd3m5/XrVuHTCYzoYhPhK7raGtrw/79+820eusx1RRozh9COp02RTUYDGLXrl2IRCK2mXN+8IMfmMdyq6poxVpNsLe315aWn8vlzMEa60P0s5/9LFpaWswfuFt/BMGvqKCAQmLc2dmJQCBgq3ZordxZrb8B3wu59ea+8soruOeee1zFNJvNYmBgAC0tLaaFTUS4//77ze0vB6ADuBjAIgBDGBNzjI7itaefRhJj0SxzPvtZdH/2s7bjqweKs64JAJs1P2fOHPzsZz8zP8+ePRvHjx/HyMgINE3Dl770Jbz11lswDAPZbBZElDd1m5cl4DaqrsTUWtaXiPBXf/VXaG9vx+c+9zmbtX/llVdi7ty5OHToEFasWIHu7m6cd9556O3txfPPPw8AtrKyVgvFWvNEPdikfokwnZhIjCORCJLJJPr6+gDAsw5QxXHzt1S7VTtqRfmA29rabP5dZ9VAa2QGEfFpAP8fgM/FWDIQAP4AwNcB/C8A/yfGolkYYJ41i3nlSuZ77uF/ufVWnvfBD3Jra6s54UMsFuNwOGzzOavSrSqaIxAImOnv559/Pq9cudLsm3WCCWsmqJeP3M0HZ/XpWaNEVCastbaKGiNQvnO3aBK3OileA5RSv0QQKg88fOTEZboCJsPixYu5v7+/Kse2+mUB5FW/c7obVq9ebZvAQcV6Au/OEm+9Rx8A8JnWVvzVeefhzGefBV5/HQDwEsbizl8E0NHRgSeeeAKjo6O2SR++9KUvobGx0ZxAoampCdu3b8fWrVtdr0VVO7zoootsPuxYLIbvfe97Ba9dXafVvdHb24uBgQFs2rQJ2WwWhmGAmW2uKOs167qORYsW4cknnzQ/33rrrVizZk3J34VY44JQPkR0kJkX561wU/dqt2pZ5G5YJyZ2wxmO57Ti3aJIzPCjH/+Yr73kEv7RuJX+QZcoDxV/bo08sc616dze2VQYpXWZtfa3+n8hy9i63GqlOwtgqawzZZGHQiFbjHg5qcWCIJQPpmv4YSFKfb1PpVJ5ou2V3GIW5vrqV/kZgJ/wEHFVgta5btasWZ6hfioM0eresNZpUJ/V9oZhmKGDEyUrWO+JcrUol4rzoWCNy7dOHiHhg4JQG7yE3PeDnYUoZcKAdDqN1atX2wYmNU1DPB5HMpnMCz/M5XJ47LHH8Ns9e/BtANdb1mmahkAggKuvvhpnnHFG3gQOAPCrX/3K3FbTNMydOxcvvPACmNl0g6hJIyKRCObPn2+6KdR1KUZHR835P3O5HG688UbMnz/f81o//elP4/jx43jkkUdMl8mdd95pi8KJRCLYsGEDiAiapiEUCqG9vV0iUQShDpnWQl5M3Cfwbuiec6LiSy+91BQqdRxN00w/OjPjk8PDAIBbMDb79BtEWLZsGdauXYtIJIJ169a5xqEDYyK+bNkyrFixwvSZDwwM4Pjx4xgYGLCNeDtHy1WEDADT162iT5whis4xAbWf9Yk+MDCQd0/Ug03XdfPBIpEoglCHuJnp1W5T6SMvxhXgVWvbWu/bmsFo3eY5SyTLYwAHHVElzkiaaDRqm/zXGRVideF4+aStqe/W6BOV4m91rzijdLzK0DrPZZ1hSea2FIT6AKeij7xYlEA5RY6I8tJwnUJ+kzUkEeCDy5e7Ht/6MLEOwDoHHyc6fyExdasF4QxBLFRi1loIy5pmrOu67cEgPnJBqA0i5BOQSqXyokOsMdteVQp1In7IIeb80EMFz+NVh9tpkavp3grVHY/FYqawutWCcJ5PVVJ0zj5kfSiUOlO9IAhTg5eQT2sfeSlEIhH09vaaM9xrmoa77747Lx7bSQ7AKgBPAZgN4FUinPWpT8F45hmgtTVve+cA7ODgoC3tFxiriaIGI++55x5s3rzZHFh0liRQMeFqgNQ5JuBMKwaAF154AceOHcO1116Ln//853jggQfw8Y9/3PR3q3OouizMbPOJS3y4INQZbupeagOwCcB/AXi2mO3r1SLv6ekxZ8FRWZjMXJR7YhHAOwEOaxq/09DAvGgR8zvvuJ6nGD9zoUprXmV7rTXQ3UIQneVmneMCzqnjVPVI5aqxzhBU7pRVgiCUDqrpWgHQAWChX4XcOdUZLK6NaDRqzmmp3BFeEyYowTv6rW+N3drrrvM830R+5mIEv5CbRu1jLR9baKIH4N066c57YhhG3oTVqhWamV4QhMriJeQVca0w814iaq3EsWqBcnc4QwRHRkZs6fO6rmP58uXYtm2bbTvDMPDlL38ZjY2N6OzsxAcjEeDNN4H167HtzTdx5pe+ZHNBFFMFrZhKa85tnG6bvr4+bN682Sya5VXpULFgwQKsW7fOPJYqAAYAP/nJT1z3OX78eMFjCoIwBbip+2QaxqrATiuL3NlUdmOhIlLmMffu5T2axr8FeLHDBTHZyI+J9nNa6FbXiyqa5UzLh8Uat+7rtMC9LHJxrwjC1IFqF80at8gfYuYPeazvBtANAC0tLYtefvnlipy3UliTZlRSzkMPPWRmTAJjSTiqPKybpWw9xpYtW3D4scdwkBlvAbjlT/4E8y+5xFaoq5TsSK96417X4VYT/POf/zzuuOMOW2EwhSq7q4p8AbCVvf3mN79pXtfQ0BD27dsH5tILaQmCMHm8imZNWdQKM28AsAEYq344VectFjd3RzqdNiNIAKC5ubngtkuXLjXdEQ3jy/8vgMcBXPXww/jkI49AH68zrmaqLzY7sthyA86+9fb2YsuWLVixYgUGBgZcRVxhnX3I6mYyDMN8aHV3d+c9VLwyZgVBmCLczPTJNPjYtTIRXgOP1oHErq4um2vmMYDfAviHAD88Hl/+ReQXxCrWLTGZrEpnVql1sNOZ5RkMBjkej9smgFZZovF4PM+lI4lBgjD1oJqDnUT0LwA6AZxJRK8B+Dozb6zEseuBQrPvqFhr55RsawHcD+DTluPcBuBJZjypabj22mtLmj2kmMHPQv22WthqhqBVq1bZZjJRxcGY2VYHxs0VJHNwCkId4abu1W7TwSLv6emZcHD0dwHe7sj6/P54/PbHr7ySP3PppQVLzlay386ZfyY7y1BXV5dY4YJQIyCZnZPHyxrWdT2v7K318x9Go0jOmYPHb7sNPQB+DeBejFnETQ8+ODZgsH8/HnjpJZx1xRUVz5Z0y+p0nsOZpel2nda5SHfu3Il9+/ZJCVtBqCfc1L3azW8WuReJRML0NRuGwdFo1OZf7unp4Z6eHiYijgD8EsDDAK/7vd/jJxyW+tNE/LuWYlfM1fdDWy1w6zyebttZxwC8Mk3FUheE6gKxyCtPd3e3bcKHrVu34n8x4zoAQWa0b9+OdzIZ/DWAEQAbAXwDwN+Mz/Np5QJmvM6M9Dvv4NHPfQ5vLlmCdd//Pp4YGUEwFKqKBWz1oWezWSQSCbOui1qvLPO1a9di3759eZEqxYZFCoJQRdzUvdptuljkVlKpFBuGwR9xVkIss33AxQKuZJ+d5XtVRcVCUTqFqjFWo5+CIIwBD4tcKyzzgpN0Oo1169YhnU7blieTSWSzWewBcBLAPwNoIMJpABo1DX999dXovflm9G/bhvu+9S2cbxiYQ4QLQiEcvf12jJxxhuv5vgNgCTPObGwseP7JoHzi1113HYLBIHRdRzAYBADXKJ1IJII1a9bYLG5VKVHtKzHlglAD3NS92s2vFnkqleLzGxpMa/n5W25hHh5m5jF/Ocat2u8D/I6m8fvG6347rVprMSozamV0lLdcfPGEFvrrAM9xHLNS16as7clMWi0+ckGoPpCJJcqnp6eHf1fT8gU2GOT/uPJKbh13UUTGl1/rmHZNHaNQedqLAoEJxfyfLHVfqoWIsyDUH15CLq6VEujs7MT/hEIwdB2XhkJ4p6VlbEUmg488+CB+yYz/AfBnRPgfAH/JjFwuh8HBQdsxvFwRkUgE39mzB7fdcguevewyz35cCyDKjE2bNlXExeKGmxtFEIT6RIS8BJRP+dZbb8Vtu3djxssvA8PDwP33A+OifjqA1cx4D4AlAGbmcmhqanI9hluERyQSwVe/9jV86PHHceT22/HrUMi1L9cByGazpv9aEIRTl4pVPyyFxYsXc39//5Sft+q89BKwaRPe+Yd/wIx33gEA9AI4+clPInTBBZNL9nnjDQx+4hNo2rMHADD6O7+DztHRqoYlCoJQn3hVPxQhLwOvuSvT+/dj/dKl+HQmg49irMTkPgA/0HWcfvXV+NQ115QmvszAD38IfOELwG9+g5HGRmz6zGdwwZ/9mYi4IJxCiJBXmHQ6jRf/9//GH2SzaCDCuWefjQaiMVfL8DByQ0OgkRGQ4/7+E4AvzJgxOUv65z8H/vIvgXQawzNnYuNnPoP2q64SMReEU4Sa1yOfbiR378bcbBYLAYAZJ0ZG8NZFF+G1EydwVmsrzmptBUIhvPL66/jepk14e3QUQxirTV5KHXIb550H7N2LV6+/Hmf90z+h7bbbsPSuu8S9IginODLYOUk6L7sMf9nQgE8R4TiAphMn8G//8R9Y2t+P87ZtQzoaBb7xDbQkErhy715kYjFsDoXwy3ITZwwD9557Li7WNKwFbAk7giCcmoiQT5JIJIJdjz+OC775Tbzy6KMYCIdxfTaLZ3M5/MnwMJK7d9u2/d73vofdu3d7RquUQmdnJ46EQtgv2ZSCIEB85BUjnU7jby67DN8ZHsYCAL/+wz/Ee++7D2htrdr5Kl32VhCE+kYGO6eAdDqNPbt2YcnTT2PRgw9C1zTo3/gGsHo1EAjUunuCIPgcLyEX10oFiUQiWLJ0Kf72wQdx08gI9g8NAfE4sGgRUKUMTEEQBIlaqTDJZBLxTAbLrQsPHwYuuQTo7saBq67CrqeeEpeIIAgVQ4S8wnR2duKTuo6Xs1n7CmYgkcDsRAIHNA23SlamIAgVQlwrFSYSieDmu+/G7xuOZ+Rf/AV+esEF0AH8fS4nYYOCIFQMEfIq0N3djX/euxffjcffXXjvvTht2TK0NjTgYk2TsEFBECqGRK1Um1/8Yiwjc5zjV12FzQsXomPpUnGrCIJQElWNWiGiy4noeSJ6kYj+phLHnDacey5w5Ij5sflHP8JNu3cj8oEP1LBTgiBMJ8oWciLSAXwXwHIAbQA+RURt5R53WjFvHnDw4LufH38c7yxYADz/fM26JAjC9KESFnkYwIvM/AtmzgD4VwAfq8BxpxcLF+LZ737X/Jh59VWMXnQRsHNnDTslCMJ0oBJCPgvAq5bPr40vs0FE3UTUT0T9J06cqMBp/ce2//5v/ImmoRdARNPw5mmnAZdfDtx9d627JgiCj6mEkJPLsrwRVGbewMyLmXnxzJkzK3Ba/9HZ2YndoRC+out4KRTCL++9F1i+HLjhBuDGG4HR0Vp3URAEH1KJhKDXALzf8vkcAMcqcNxph5qvUxW7+nAkAlx2GbBmDXDbbTh54AB+cPnluHj5coloEQShaMoOPyQiA8ALAJYC+BWAJwH8OTP/1GufUyr8sEhevPlmtPT04BkASxoasPPxx0XMBUGwUbXwQ2YeBXAjgEcBHAVwfyERF9z5t9NPR5em4WsAhkdGJOtTEISiqUitFWZ+BMAjlTjWqUpnZyduDYWQyWQk61MQhJKQoll1gtN/Lm4VQRCKRYS8johEIiLggiCUjBTNEgRB8Dki5IIgCD5HhFwQBMHniJALgiD4HBFyQRAEnyNCLgiC4HNEyAVBEHyOCLkgCILPESEXBEHwOSLkgiAIPkeEXBAEweeIkAuCIPgcEXJBEASfI0IuCILgc0TIBUEQfI4IuSAIgs8RIRcEQfA5IuSCIAg+R4RcEATB54iQC4Ig+JyyhJyI/pSIfkpEOSJaXKlOCYIgCMVTrkX+LICPA9hbgb4IgiAIk8AoZ2dmPgoARFSZ3giCIAglIz5yQRAEnzOhRU5EOwE0u6y6mZn/vdgTEVE3gG4AaGlpKbqDgiAIQmEmFHJmXlaJEzHzBgAbAGDx4sVciWMKgiAI4loRBEHwPeWGH15FRK8BiAB4mIgerUy3BEEQhGIpN2rlRwB+VKG+CIIgCJNAXCuCIAg+R4RcEATB54iQC4Ig+BwRckEQBJ8jQi4IguBzRMgFQRB8jgi5IAiCzxEhFwRB8Dki5IIgCD5HhFwQBMHniJALgiD4HBFyQRAEnyNCLgiC4HNEyAVBEHyOCLkgCILPESEXBEHwOSLkgiAIPkeEXBAEweeIkAuCIPgcEXJBEASfI0IuCILgc0TIBUEQfI4IuSAIgs8pS8iJ6DYieo6IniGiHxFRY4X6JQiCIBRJuRb5YwA+xMwXAHgBwJryuyQIgiCUQllCzsw7mHl0/OMTAM4pv0uCIAhCKVTSR34NgO1eK4mom4j6iaj/xIkTFTytIAjCqY0x0QZEtBNAs8uqm5n538e3uRnAKID7vI7DzBsAbACAxYsX86R6KwiCIOQxoZAz87JC64no0wA+CmApM4tAC4IgTDETCnkhiOhyADcBWMLMb1emS4IgCEIplOsjvwvAewA8RkSHiOj7FeiTIAiCUAJlWeTMfH6lOiJUn3Q6jWQyic7OTkQikVp3RxCEClGWkAv+IZ1OY+nSpchkMggGg9i1a5eIuSBMEyRF/xQhmUwik8kgm80ik8kgmUzWukuCIFQIEfJThM7OTgSDQei6jmAwiKamJqxbtw7pdLrWXRMEoUzEtXKKEIlEsGvXLiSTSZw8eRI33HADcrkcQqGQuFkEweeIkNcp5QxMeu2r/r9kyRKMjo5VVhgeHkYymRQhFwQfI0Jeh5QzMDnRvslkEtls1vysaRo6OzsrfQmCIEwh4iOvQ5wDk319fUX7syca1Ozs7EQoFIKmaQgEAvjud78r1rgg+ByxyOsQNTCZyWRgGAY2bdqEbDZblHVu3TcYDJrWttXdonzlEk8uCNMDEfI6JBKJ4POf/zweeOABnH322fjxj39ss7DdxLeQULu5W9askdLxgjBdECGvQzZs2ID169cDAF588UUEAgEAsFnYVqxCbRgGrr76aqxatcoUfKu7ZXh4GGvXrsXatWvN9eohcPLkSRw6dAgrVqxAd3f31FysIAhlI0Jeh2zZssX2ub29HdFo1NMVYhXqbDaLRCKBzZs3m24Y5W4ZHh5GLpfDY489hj179mD37t0AgKVLl2JoaAiqeOWOHTsAQMRcEHyCDHbWmHQ6nTeQuWLFCts21157LdasWZMn4mrfpqYmBINBEBEAgJltA50qhnzx4sXm+uHhYfT19ZkPAWcFYufDRBCE+kUs8hKpZOEpr1BBZQlv2bLF083h3Le3txcDAwO2gVGrGyYSiWDhwoU4cOCA7TjKWrda5ED+w6Ra90AQhArAzFPeFi1axH4klUrxjBkzWNd1njFjBqdSqbKO19PTw7quMwDWdZ17enrK3jeVSnFPT49r31KpFAeDQSYiDgaD5jZqn3g8zl1dXZxIJDzPW+l7IAhC8QDoZxdNFYu8BNxitCeKIJlMqGAxeO0biUQ8zxmJRJBMJpFMJtHU1GRzvRRrWRd7DwRBmDpEyEvATTydom11eei6jmuuucYWQeLcfrIx3c59AWDdunUTHketW7p0KYaHh6HrOu66666iBzabmppARNA0reSHjyAI1YG4BtNsLl68mPv7+6f8vJXAKsQA8nzcyWQSf//3f2+mwRMRGhoaTB/2D37wA4yOjpp+7cHBQfNY1uMqq9ltvdugZykp/evWrcPf/d3fIZfLAQAMw8DevXsLnsN6nsk8AARBKB8iOsjMi/NWuPlbqt386iN34uanVj5kImIADIA1TWPDMPKWBQIB1nWdQ6EQB4NB1nWdg8Egh0Ihc1vlzw6FQua20WiUY7GY6Z+29kPTNO7q6uJUKuXqL0+lUhyLxVjTNLMvADgajRb0fadSKe7q6jL3s55HEISpAR4+chHyMvAa+FNiqcQ3EAjYhJOI2DAMcxkR2YTbKvheywBwKBQyBTsUCtm213WdDcMwxT8Wi3Eikch7yLidwznwqq7T2l8l5jLgKQhTh5eQi4+8CNz84IXqlqjBw1WrVpkuktWrV9v85u3t7eYywzDAzMhms9B13UzsURARiAi5XM4WIqgGGzs7O23L1bEUKklIHdu6rXUfAK6+bzXAmcvloGkazj33XPziF79ALpeTAU9BqAOmhZBPJq7ZbR+vZc54bSXAE9UtsUaDzJ8/P+/Y1mXAu/7pw4cP44YbbkA2m4WmaaaIa5oGZjZ920pwnaVp3VD7EZGrkANjD4xly5aZ6fvqfqiEI3XNX/3qV233QAY8BaHGuJnp1W6VdK0Ucm+4xVOnUimORqOs67rp4kgkEp7HicViNlfC+eefb7oYdF3nWCzmGbdd7nX19PRwLBaz+b+tPnkV723tu5svXrlYgsEgNzc357lVVPM6ZigU4o6ODg6Hw7b11bhuQRC8QTV85ABuBfAMgEMAdgA4u5j9KinkhQYcnf7hWCzGhmHkCVggELAJpvU4Vt+ztVkF0stXrHzl1oFJtbxYEfQSaU3T8vzYPT09nEgkzOtQ26k+BAIBTxFXTV2H9QFmbdZEIkEQppZqCfkZlv9/AcD3i9mv2ha5lwi5LbOKnRrQ03XdtEC99mlra7MNYLoJazAYzLN2J5MZac28tPYhGo3m7d/T05PXL3VPvK7d7V64PfDUPVQPDJUFKta5IEwNVRFy24GANQC+V8y2lY5asQqJU0CLaYFAgFOpVJ41W6jNnj07zy1hFTI34TQMg6PRaMHoEOt1OMXRKdJEZL5xWF1KTstbvY0Ucx2apvGCBQs8r9swjLwHiorAUW4qQRCqQ9WEHMA3AbwK4FkAMwts1w2gH0B/S0tL1S7U6mpRYXhuYX5WMYzH4xyNRnnWrFkFxdtp6VrXdXR0mMLr9UBQ/VGfnfVOrO4gFVfu9Pu7hQ+qWPNoNMrRaJQbGxtt68855xxOJBKubiKvNw6v1tbWxuecc86ED0VBECrPpIUcwM5xkXa2jzm2WwPglomOx1WwyK04XReqEJQ1mcXaFixYUJTvGABfeOGFtnhvp8hrmmYKsFMsVQKQ9aESi8XMfjsfQIWs9mL93U6BjUajfOGFF7qKeamCXuhh5/WWIQhCeXgJ+YThh8y8bKJtxvlnAA8D+HqR25eNW7igtQaJit9WEyq48cwzz+Sta2howNDQUN62mUwGgUDAtX63OsbIyEjefpqm4ayzzsKSJUvwwAMPmGF7q1atMrdxztPJzJ7laAFg06ZNBUMJnYyMjGDr1q2e64s9jpX3ve99ePPNN23LiAhNTU0A3EM3VckBiTsXhAripu7FNgBzLP//PID/V8x+lbDIixk0tFq5Xs3LGnWz3ufNm8fRaLTg8ZRrxMvCjcfjBcvMFvKRW6/LrX+VaCrE0m259V6GQiFOJBKubwbq+3CWDlBRPpINKgiTA1XK7PwWEX0AQA7AywBiZR6vaIopp+o1aYIVTdMwf/58HDp0yLbczYJ//vnn8fOf/xyapuWtJyIYhoG77roL8+fPR19fHzZt2oRMJmPb7tChQ3j00Udd++IsJ+tltTY1Nbme3+sai0FljwLAK6+8gkAggNHRUTAzNE1DKBQyC38BsFV0VMlL6vzWjFP1lqFpGrLZrGSDCkI1cFP3arepssjVdsqn7ObXLhTRAYt/2+qz7ujosG1jGEZerLg6t9OCr0RUh1v0ihoPKNZ37jZg6pbolEgkJvRxO2vLuCVmqTovYpELwuTBdCyaVWpijRIUZ5KOSvxRUSUq0sVZbMoqQolEgsPhsGsstxNrzHWlrtuZ8OQs2OXM8FQPLrW9M4Rw5cqVZQvtRN+HDHwKQnl4CbnUIx/HWWe8mDostWSi/ljrpAwODtpqm6vtN2zYYJsXtN6uURAEO171yEXIBUEQfIKXkGu16IwgCIJQOUTIBUEQfI4IuSAIgs8RIRcEQfA5IuSCIAg+R4RcEATB59Qk/JCITmAspb/SnAngjSocd6rwe/8B/1+D3/sPyDXUA9Xq/2xmnulcWBMhrxZE1O8WY+kX/N5/wP/X4Pf+A3IN9cBU919cK4IgCD5HhFwQBMHnTDch31DrDpSJ3/sP+P8a/N5/QK6hHpjS/k8rH7kgCMKpyHSzyAVBEE45RMgFQRB8zrQSciK6lYieIaJDRLSDiM6udZ9KhYhuI6Lnxq/jR0TUWOs+lQoR/SkR/ZSIckTkmxAyIrqciJ4noheJ6G9q3Z9SIaJNRPRfRPRsrfsyGYjo/US0m4iOjv9+vljrPpUKETUQ0QEienr8Gm6ZkvNOJx85EZ3BzG+N//8LANqYecrmEa0ERNQF4HFmHiWibwMAM99U426VBBHNw9g8rgkAX2Hmui8+T0Q6gBcA/BGA1wA8CeBTzHykph0rASLqAPAbAH3M/KFa96dUiOgsAGcx81NE9B4ABwFEffYdEIDfYebfEFEAwH4AX2TmJ6p53mllkSsRH+d3MDaNma9g5h3MPDr+8QkA59SyP5OBmY8y8/O17keJhAG8yMy/YOYMgH8F8LEa96kkmHkvgDdr3Y/Jwsz/ycxPjf//fwAcBTCrtr0qjfEZ2X4z/jEw3qquQ9NKyAGAiL5JRK8CWAnga7XuT5lcA2B7rTtxijALwKuWz6/BZyIynSCiVgDtAH5S466UDBHpRHQIwH8BeIyZq34NvhNyItpJRM+6tI8BADPfzMzvB3AfgBtr21t3JrqG8W1uBjCKseuoO4q5Bp9BLst890Y3HSCi0wFsAbDa8ZbtC5g5y8wLMPY2HSaiqru5jGqfoNIw87IiN/1nAA8D+HoVuzMpJroGIvo0gI8CWMp1OohRwvfgF14D8H7L53MAHKtRX05Zxv3KWwDcx8wP1Lo/5cDMJ4koCeByAFUdgPadRV4IIppj+XglgOdq1ZfJQkSXA7gJwJXM/Hat+3MK8SSAOUT0+0QUBPBJAA/WuE+nFOMDhRsBHGXm22vdn8lARDNVpBkRzQCwDFOgQ9MtamULgA9gLGLiZQAxZv5VbXtVGkT0IoAQgMHxRU/4MPLmKgB3ApgJ4CSAQ8z8kZp2qgiI6I8B9ALQAWxi5m/WtkelQUT/AqATYyVUXwfwdWbeWNNOlQARXQpgH4DDGPsbBoC/ZeZHater0iCiCwBsxthvSANwPzN/o+rnnU5CLgiCcCoyrVwrgiAIpyIi5IIgCD5HhFwQBMHniJALgiD4HBFyQRAEnyNCLgiC4HNEyAVBEHzO/wfsNBCdQvyqdQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.plot(phi(trajectory), psi(trajectory), 'k.')\n", "plt.plot(phi(subtrajectories[0]), psi(subtrajectories[0]), 'r')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Setting up the normal temperature engine\n", "\n", "We'll create another engine that uses a 300K integrator, and equilibrate to a 300K path from the 500K path." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "integrator = VVVRIntegrator(\n", " 300*unit.kelvin, \n", " 1.0/unit.picoseconds, \n", " 2.0*unit.femtoseconds\n", ")\n", "integrator.setConstraintTolerance(0.00001)\n", "engine = omm.Engine(\n", " template.topology, \n", " system, \n", " integrator, \n", " openmm_properties=openmm_properties,\n", " options=engine_options\n", ").named('300K')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Using TPS to equilibrate\n", "\n", "This is, again, a simple path sampling setup. We use the same `TPSNetwork` we'll use later, and only shooting moves. One the initial conditions are correctly set up, we run one step at a time until the initial trajectory is decorrelated.\n", "\n", "This setup of a path sampler always consists of defining a `network` and a `move_scheme`. See toy model notebooks for further discussion." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "network = paths.TPSNetwork(initial_states=C_7eq,\n", " final_states=alpha_R).named('tps_network')\n", "scheme = paths.OneWayShootingMoveScheme(network, \n", " selector=paths.UniformSelector(),\n", " engine=engine).named('equil_scheme')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The move scheme's `initial_conditions_from_trajectories` method will take whatever input you give it and try to generate valid initial conditions. This includes looking for subtrajectories that satisfy the ensembles you're sampling, as well as reversing the trajectory. If it succeeds, you'll see the message \"No missing ensembles. No extra ensembles.\" It's also a good practice to use `scheme.assert_initial_conditions` to ensure that your initial conditions are valid. This will halt a script if the initial conditions are not valid." ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "No missing ensembles.\n", "No extra ensembles.\n" ] } ], "source": [ "# make subtrajectories into initial conditions (trajectories become a sampleset)\n", "initial_conditions = scheme.initial_conditions_from_trajectories(subtrajectories)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "# check that initial conditions are valid and complete (raise AssertionError otherwise)\n", "scheme.assert_initial_conditions(initial_conditions)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can create a `StepVisualizer2D`, which plays the simulation as it is running. It requires CVs for the $x$ direction and the $y$ direction, as well as bounds of the plot. You can set the `background` attribute to an existing matplotlib `Figure`, and the simulation will plot on top of that. See the toy model MSTIS example for an example of that.\n", "\n", "This isn't necessary, and isn't even a useful for long simulations (which typically won't be in an interactive notebook), but it can be fun!" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "visualizer = paths.StepVisualizer2D(network, phi, psi, xlim=[-3.14, 3.14], \n", " ylim=[-3.14, 3.14])" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "sampler = paths.PathSampling(storage=paths.Storage(\"ad_tps_equil.nc\", \"w\", template),\n", " move_scheme=scheme,\n", " sample_set=initial_conditions)\n", "sampler.live_visualizer = visualizer" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# initially, these trajectories are correlated (actually, identical)\n", "# once decorrelated, we have a (somewhat) reasonable 300K trajectory\n", "initial_conditions[0].trajectory.is_correlated(sampler.sample_set[0].trajectory)" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAD4CAYAAAD2FnFTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAU1ElEQVR4nO3deZAdZbnH8d+TWbKcIKEwmJgEiSZgEFAg4IIiskgusgQtFheugDIFIhgtJCCLJVQoNpcIliTs12IXUiJLgSy5YEjEIbIkGYIRCQkCiSySmZCESZ77x3vmziS8M+ec6Z7TZ/l+qqbO6eV0P11Q7y/db7/d5u4CAGBLg7IuAABQmQgIAEAUAQEAiCIgAABRBAQAIKoxi50OGjTIhw4dmsWuAaBqrV271t29bP+wzyQghg4dqo6Ojix2DQBVy8zeLef+EieRmQ0xsyfN7BkzW2xmP0ujMABA8QaiLU7jDGK9pP3dvd3MmiT92czud/cFKWwbAFCc1NvixAHhYSh2e36yKf/H8GwAKKOBaItT6ewwswYze1rSKkl/cve/RNZpMbNWM2vt7OxMY7cAUG8au9rR/F9Lz4XFtMWlsDSfxWRmIyTNkXSauy/qbb1cLud0UgNAacxsrbvnilhvhIpoiwtJ9XYpd39b0lxJU9LcLgCgeGm1xWncxTQyn1Yys6GSDpT0fNLtAgCKNxBtcRp3MY2WdKOZNSgEzu3ufk8K2wUAFC/1tjjVPohi0QcBAKUrtg8iLTyLCQAQRUAAAKIICABAFAEBAIgiIAAAUQQEACCKgAAARBEQAIAoAgIAEEVAAACiCAgAQBQBAQCIIiAAAFEEBAAgioAAAEQREACAKAICABBFQAAAoggIAEAUAQEAiCIgAABRBAQAIIqAAABEERAAgCgCAgAQRUAAAKISB4SZjTOzR82szcwWm9kP0igMAFCatNtjc/ekBY2WNNrdF5rZVpKekjTV3Zf09ptcLucdHR2J9gsA9cbM1rp7ro/lJbfHfUl8BuHur7r7wvz3NZLaJI1Jul0AQGnSbo8b0ypMksxsB0m7S/pLZFmLpBZJam5uTnO3AFAvGs2stcf0bHefHVuxr/a4WIkvMfUoZrik/5U0w93v6mtdLjEBQOkKXWLqsV7R7XFfUrmLycyaJN0p6aYkxQAAkkmzPU7jLiaTdK2kNnf/RdLtAQD6J+32OI0ziH0kHSdpfzN7Ov93SArbBQCUJtX2OLU+iFLQBwEApSu2DyIttTGSev36rCsAgJpT3QGxYoV00EHSl76UdSUAUHOqOyCGD5cef1yaP1964onu+bNmSWPHSi+/nF1tAFDlqjsgttlGOv308P3MM8Onu3TyydIrr0hXXRX/3euvhz8AQK+qOyAkafp0qalJmjdPevLJcDbRpasj3F166inpvPOk8eOlUaPC57p12dQMAFUg1UdtZGLbbaXvfU+aOVM66yxpwoTuZQ8/LB13nHTPPdLbb2/+u85OafDgspYKANWkNm5zXbVKGjMmNPpDhsTPDEaMCGcaq1eH6fvvl6ZMSa8GABhg3ObaH9ttJ7W0hO/r1oWzir33lnbaSTr//HB5adq07nD4xjcIBwAooDbOICTp1VelceOkjRtDh/Ull4T569aFQJgzJ0xfcIF07rmSWbr7B4ABxhlEf40eHc4WRo2STjklzHvtNWmvvUI4NDWFz/POIxwAoAi1cwaxpSVLpH33ld54Q/rgB0OH9W67Dew+AWAAcQaRBnfp2GNDOOy2m7RoEeEAACWqzYB45BHpuefC96VLw2WnK6+U/vOfbOsCgCpSmwFx4YXd37se5HfaadLIkeHMYv78cJYBAOhVbQbE7rvH57/3nnTbbdLnPhdGUv/qV9Kbb5a1NACoFrXbSb1mTRj/sGCBdPXV0oc/LF1+uXTLLdK110rt7WG9hgbpsMPCOIl99+UOJwAVq9yd1LUbEH3ZsEG6+27p5z8PAdJlzBjp5ptDUEyZIi1cKL30kjRsWGalAkAXAqLc/vnP8NTXWbNCJ/bWW4dQ2GabsHyvvcJDAAEgYwREVjo7pT32CHc/7b9/eBz44sVh2YoV4f0SAJAhAiJLixdLn/pUCIueJk0KA+8AIEMMlMvSJz4RgmDq1M3nt7VJc+dmUREAZIaA2NLEieGZTYsXS4cc0j3/hBPCbbIAUCcIiJglS8LdTLNnh/6HbbcNHddXXJF1ZQBQNvRBxJx0knTNNeHlQ1dcEUZgT50a3kC3fLn0oQ9lXSGAOkQndSVYsiQ83G/jxjD92c+Glw0tWyYdc4x0663Z1gegLlVlJ7WZXWdmq8xsURrby9zOO0sXXdQ9PX9+CAcpPKqj5+A6AKgQabfFafVB3CCptt7hecYZ0p57xpcdf3z32QUAVI4blGJbnEpAuPtjkmrrqXeDBkl33CE1N4fpqVOlrbYK35culW68MbPSACAm7ba4bHcxmVmLmbWaWWvnlgPRKtX48eE9EpL0wAPS449LRx4Zpl9+Obu6ANSrxq52NP/XMpA7S62T2sx2kHSPu+9SaN2K76TuyT08emPuXGmffUJIrFoV7mwaxF3CAMqnmE7qUtriQmjhCjELjwjP5aR586QZM8K4CMIBQI2jlSvGqFHSddeF7+edF94tcdllvMIUQE1L6zbXWyTNl7STma00s++ksd2KcvTRYfDcqFFhTMSZZ4YBc6ecIr34YtbVAUDqbTED5Uq1caN0773SBReEN9Z1OeigcHbx+c/zVjoAA4KR1NXkb38LA+ruukvatCnM6+rIJiQApIyAqEavvirNnBn6JTZtklauDK8vBYAUVeWjNure6NHSxRdLEyaE6ZdeyrQcAEgDAZGmj340fBIQAGoAAZGmSZPCJ3c1AagBBESaJk4Mn21t2dYBACkgINL0kY+Ez65Hg3d57bUwdgIAqggBkaYddgifK1aEz3XrpB//OIy8/vjHpfb2zEoDgFIREGl46SVpl12kc88N06tXS488In3sY9Lll4cH/r35pvTrX2daJgCUgnEQaZg2LYyDiNl+e+nkk6Wf/EQaPjyMmRg+vKzlAagNjIOoNu7d76g+7LD3Lz/ySOlb3wqvMW1vDy8d+sc/ylsjAPQDAZHUokXS66+Hs4I5c6Tf/S7M73oc+MyZ4SyiZyhMmCBNny6tXVv+egGgSAREUrffHj4PP1xqaAhnC+5SZ6f02GPSV74Snsu0fv3mv7v00hAct98e1geAClO/ATFtmrT77tKaNcm2c9NN4fOb39x8vpn0hS9I99wjvfJKdwd2T2+8IR1zjPSZz0iLFyerAwBSVr+d1KNHh/EJM2dKp5/ev20sXx5ubTWT3npL2nrrvtd/7z2pubn35d/9brjrqdB2ANQlOqnLZfDg8Hnppd2XeFavDo/wvu++8Aa5GTOklhbp4IOlXXcNLwsaPDj8i9+9+5lL7tIJJxTeZ1OTdP754funPy2deOLmy6+5RhoxQrrwwu7HhwNARur3DGKnnaQXXgjfH3wwhMOWl4n6sueem78wqKVFmjWruN8uXSqNGycNGya98450443xs5jnngvjKwBAvA+ifD75SenZZ8P3/faTdtxRmj1b2nbbMMBt7Njw6Izttw9nDqNGhVtUJ0/u3kZjY3gV6RlnhP6MJNylRx+VDjige96hh0p//GOy7QKoGeUOiMZy7ajiDBnS/X3u3O67jC6+OPQFxGzaFM482tulH/5Q+s53wiWhNJhJ++8fgqKtTbr+eun449PZNgD0Q30GhHt3QGyzTehgnj8/TOf6COdBg6Tnnx/4+iZNCn0jAJCh+guIDRvC5aTly8P0W29tvryvgACAOlJ/dzE1NPS9/IQTpF/+MpxRrFtXnpoAoALVZyf16tXSdtsVXm/QoNBh/cUvSmef3f1KUQDIAJ3U5bBxY+/Lpk8Pl5/mzQvvdfj738Pfxo1hbAQA1InavcS0YUMY4GYWnofU2todDJ2dvf/uooukW24J75X+8pfDvGHDwq2sAFBHavcS05o10gc+8P75kyeHAXE33CA980z8tyNHhtte33knjJyeNy8MjAOADFXlozbMbIqZLTWzZWZ2VhrbTGyrrcJAuPHjN5/f2hrGMHSFw9ix7//t6tUhHIYNkx5+mHAAUBXSbosTB4SZNUj6jaT/krSzpK+b2c5Jt5uKXXcNl4o2bJB+//t4GGy/vTR06ObzBg+W7rhD+ve/pX32KU+tAJDAQLTFaZxB7C1pmbu/6O4bJN0q6YgUtpuepibpa18Lnc5Ll0qnntodCk88Ib37bhicdttt4e6m9eulo46SfvSjMGoaACpf6m1xGgExRtKKHtMr8/M2Y2YtZtZqZq2dfXUSD7Qdd5SuvDK8i+H668MzmfbbL4x7OProcMZx2mlh3auuCpeoDj00PKbbLHyeemp29QOoZ41d7Wj+r6XHsqLa4lIk7qQ2s6MkHezu381PHydpb3c/rbffZD4OohgLF0rHHhtucd2SmXTKKdJvflP+ugDUrb46qfvTFheSxhnESknjekyPlfSvFLabrT32kJYs6X63dE/u0tVXl78mAOhd6m1xGgHxV0kTzWy8mTVLOlbS3SlsN3uNjb2/uOe998pbCwD0LfW2OPFIanfvNLPvS3pAUoOk69y9dl6w3NQUD4OmpvLXAgC9GIi2OJVHbbj7fZLuS2NbFeekk6Tf/rb7taRS6IM46aTsagKAiLTb4vp8FlMpujqir746nEk0NYVwoIMaQI2r3UdtAECNqcpHbQAAag8BAQCIIiAAAFEEBAAgioAAAEQREACAKAICABBFQAAAoggIAEAUAQEAiCIgAABRBAQAIIqAAABEERAAgCgCAgAQRUAAAKIICABAFAEBAIgiIAAAUQQEACCKgAAARBEQAIAoAgIAEEVAAACiEgWEmR1lZovNbJOZTU6rKABAuvrTXic9g1gk6auSHku4HQDAwCq5vW5Msjd3b5MkM0uyGQDAAOtPe50oIEphZi2SWiSpubm5XLsFgFrSaGatPaZnu/vsAdtZoRXM7CFJoyKLznH3PxS7o/xBzJakXC7nRVcIAOjS6e699h+k1V53KRgQ7n5gqRsFAJRf2u01t7kCAKKS3uZ6pJmtlPRZSfea2QPplAUASFN/2mtzL393QC6X846OjrLvFwCqmZmtdfdcufbHJSYAQBQBAQCIIiAAAFEEBAAgioAAAEQREACAKAICABBFQAAAoggIAEAUAQEAiCIgAABRBAQAIIqAAABEERAAgCgCAgAQRUAAAKIICABAFAEBAIgiIAAAUQQEACCKgAAARBEQAIAoAgIAEEVAAACiCAgAQFSigDCzy8zseTN71szmmNmIlOoCAKSoP+110jOIP0naxd13k/SCpLMTbg8AMDBKbq8TBYS7P+junfnJBZLGJtkeAGBg9Ke9TrMP4kRJ9/e20MxazKzVzFo7Ozt7Ww0A0LvGrnY0/9fSz+302V53MXfvewWzhySNiiw6x93/kF/nHEmTJX3VC21QUi6X846OjkKrAQB6MLO17p7rY3mq7XVjoYLc/cACBX9b0qGSDigmHAAAAyPt9rpgQBTY2RRJ0yV90d3XJtkWAGDg9Ke9LniJqcAOl0kaLOmN/KwF7n5yod9xiQkASlfoElOB35bcXic6g3D3CUl+DwAoj/6014ykBgBEERAAgCgCAgAQRUAAAKIICABAFAEBAIgiIAAAUQQEACCKgAAARBEQAIAoAgIAEEVAAACiCAgAQBQBAQCIIiAAAFEEBAAgioAAAEQREACAKAICABBFQAAAoggIAEAUAQEAiCIgAABRBAQAIIqAAABEERAAgKhEAWFmF5rZs2b2tJk9aGYfTqswAEB6+tNem7sn2eEH3P2d/PfTJe3s7icX+l0ul/OOjo5+7xcA6pGZrXX3XD9/W3J7negMomtneTlJ/U8bAMCA6U973Zh0p2Y2Q9J/S/qPpC/1sV6LpBZJam5uTrpbAKhHjWbW2mN6trvPLvbHxbbX/79+oUtMZvaQpFGRRee4+x96rHe2pCHu/tNCO+USEwCUrtAlprTb60R9EFsU9hFJ97r7LoXWJSAAoHRJ+iC22E5R7XXSu5gm9pg8XNLzSbYHABgY/Wmvk/ZBXGxmO0naJGm5pIJ3MAEAMlFye53aJaZSmNkmSe/2srhRUmcZyymXWjwujql61OJx1eMxDXX3sg1wziQg+mJmre4+Oes60laLx8UxVY9aPC6OaeDxqA0AQBQBAQCIqsSAKHrQR5WpxePimKpHLR4XxzTAKq4PAgBQGSrxDAIAUAEICABAVEUGRC2+Z8LMLjOz5/PHNcfMRmRdUxrM7CgzW2xmm8ysYm7P6w8zm2JmS81smZmdlXU9aTCz68xslZktyrqWtJjZODN71Mza8v/v/SDrmpIysyFm9qSZPZM/pp9lXZNUoX0Q/X3PRCUzsy9LesTdO83sEkly9+kZl5WYmU1SGJk5S9IZ7t5a4CcVycwaJL0g6SBJKyX9VdLX3X1JpoUlZGb7SmqX9D/FPCetGpjZaEmj3X2hmW0l6SlJU6v5v5WZmaScu7ebWZOkP0v6gbsvyLKuijyDqMX3TLj7g+7eNUJygaSxWdaTFndvc/elWdeRgr0lLXP3F919g6RbJR2RcU2Juftjkt7Muo40ufur7r4w/32NpDZJY7KtKhkP2vOTTfm/zNu9igwIKTy33MxWSPqmpPOzridlJ0q6P+sisJkxklb0mF6pKm906oGZ7SBpd0l/ybiUxMyswcyelrRK0p/cPfNjyiwgzOwhM1sU+TtCktz9HHcfJ+kmSd/Pqs5SFDqm/DrnKDxr5absKi1NMcdVAywyL/N/waF3ZjZc0p2Spm1x1aEquftGd/+UwtWFvc0s80uCid8o11/ufmCRq94s6V5JBV9ElLVCx2Rm35Z0qKQDvBI7f3pRwn+rarZS0rge02Ml/SujWlBA/jr9nZJucve7sq4nTe7+tpnNlTRFUqY3F1TkJaZafM+EmU2RNF3S4e6+Nut68D5/lTTRzMabWbOkYyXdnXFNiMh36F4rqc3df5F1PWkws5Fddzaa2VBJB6oC2r1KvYvpTkmbPbfc3V/JtqpkzGyZpMGS3sjPWlDtd2ZJkpkdKekKSSMlvS3paXc/ONOi+snMDpH0K0kNkq5z9xnZVpScmd0iaT9JH5T0uqSfuvu1mRaVkJl9XtLjkp5TaCMk6Sfufl92VSVjZrtJulHh/71Bkm539wuyrapCAwIAkL2KvMQEAMgeAQEAiCIgAABRBAQAIIqAAABEERAAgCgCAgAQ9X90q1RBXSSh/AAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Step 6: All trajectories decorrelated!\n" ] } ], "source": [ "sampler.run_until_decorrelated()" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "False" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# now we have decorrelated: no frames are in common\n", "initial_conditions[0].trajectory.is_correlated(sampler.sample_set[0].trajectory)" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAD4CAYAAAD2FnFTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAW8UlEQVR4nO3deZBU1d3G8ecHM4M4CLglIoMbLokLFRMgmyLuS4xG4xqNGpeJBBJTWRWs+KqFvgmJppJYGlBL4lrugagxoriLOi5RkCVoNCIQokGFGQVafu8fZ+btGTgz3T19u2/3zPdT1dV9+96+91ylzjP3nHPPNXcXAAAb6pN2AQAAlYmAAABEERAAgCgCAgAQRUAAAKJq0jhonz59vH///mkcGgCqVktLi7t72f6wTyUg+vfvr+bm5jQODQBVy8w+Kufxik4iM9vEzJ4zs7+b2TwzuziJggEA8leKujiJK4g1kg5w99VmVivpSTN7wN3nJLBvAEB+Eq+Liw4ID7dir25drG19cXs2AJRRKeriRDo7zKyvmb0saYWkh9z92cg2jWbWZGZNmUwmicMCQG9T01aPtr4a26/Mpy4uhCU5F5OZDZZ0j6Tvu/vczrarr693OqkBoDBm1uLu9XlsN1h51MW5JDpcyt3fl/SopMOS3C8AIH9J1cVJjGLaujWtZGb9JR0kaUGx+wUA5K8UdXESo5iGSJpuZn0VAud2d/9LAvsFAOQv8bo40T6IfNEHAQCFy7cPIinMxQQAiCIgAABRBAQAIIqAAABEERAAgCgCAgAQRUAAAKIICABAFAEBAIgiIAAAUQQEACCKgAAARBEQAIAoAgIAEEVAAACiCAgAQBQBAQCIIiAAAFEEBAAgioAAAEQREACAKAICABBFQAAAoggIAEAUAQEAiCIgAABRRQeEmQ0zs9lmNt/M5pnZeUkUDABQmKTrY3P3Ygs0RNIQd3/RzDaT9IKkb7j7a539pr6+3pubm4s6LgD0NmbW4u71XawvuD7uStFXEO6+zN1fbP28StJ8SUOL3S8AoDBJ18c1SRVMksxsB0l7S3o2sq5RUqMk1dXVJXlYAOgtasysqd3yVHefGtuwq/o4X0U3MbUrzABJj0ma7O53d7UtTUwAULhcTUzttsu7Pu5KIqOYzKxW0l2Sbi6mMACA4iRZHycxiskkXSdpvrtfUez+AADdk3R9nMQVxFclfVvSAWb2cuvriAT2CwAoTKL1cWJ9EIWgDwIACpdvH0RSeuad1M3N0hlnSO+8k3ZJAKBq9cyAmDRJmj5d2m03aebMtEsDAFWpZwbEBRdI++wTriSOOkoaN05asybtUgFAVemZAfHpT0uPPSZNnhyWr7lG+uIXpU8+kZqapEsvlVatSreMAFDhen4ntVn287hx0tVXh8/Tp0unnVaeMgBAAuikLsZ770k/+5n0xBPZ7447Lvu5LRwkaeXK8pULAKpQ9QfE//yPtO++0vz50i23SFOmSGPGSIccIr30kjRwYPx3DLMFgC5VfxPTfvtJjz8u1dRIW20lLV+e3+9mz5bGjk2mDABQBuVuYqr+gHj+eWn06I7f7bqrtGhR179bty6ECgBUCfogCjVqVLgprr2uwmHLLaVTTiEcACCH6r+CkEKH83bbSatXx9fff38Y5rr55h1HNQFAFeEKojs231z6/e87X3/eeeGKgXAAgLz1jICQpNNPl/beO7v8m99kP//jH9KgQdJr3XosKwD0Sj0nIMxCh/WJJ4blW2+VTjih4zZ77CF95jPSM89IKTStAUA16Rl9EO3NnSvttVfu7bbfXvrhD8Pd1FtsEd/m6aelF1+Uvvc9qU/PyVIA1Ylhrkl46inphRekt96SXn9devLJcGPcxx9vvG3fvtLXvx7CYsyYjv0U224rLVsmnXRSuCIBgBQREKX0/vuhQ7vNl74kzZmTXR46NNyNPWZMWN5xR+nNN8Pnf/1LGjasXCUFgI0wiqmUBg+W1q+XZsyQXnkl9EW88UaYv2nQoPCAoaOOCkEiSTfemP3t6afTbwGgV+ldVxBdyWSkz39eevVV6XOfkxobpQMPDA8danPbbdlOcAAoM5qY0jRvXgiJtWvj6zfbLDQ5ddapDQAlRBNTmvbYIzQ5XXONdOihUv/+HdevWiWNH59O2QCgzLiC6Monn4Qpw595JvRRnHlm+G7WrND8BABlxBVEJenbVxo5Unr33dBJ3a9f+P7UU6WWlnTLBgAlRkAUoi0Uli+XJk5MtywAUGI0MeVj/nxp9907fmcmNTWFTm0AKIOqbGIys+vNbIWZzU1ifxXns5/N3iS3xx7h3V369rfTKxMAbCDpujipJqYbJB2W0L4q02mnhfe99pJmzgxDXleu7Hz7Dz4ILwAonxuUYF2cSEC4++OS/pvEvirWt74V3mfMCENgV66UFi+Ob7tqlbTNNuHO7W23la66KtylDQAllHRdXLZOajNrNLMmM2vKZDLlOmxydt89NDO1tEiPPBJGOG26aXzbAQOkNWvC52XLpAkTpIaGMNX45MnSggUb/2bNGkZGAcilpq0ebX01lvJgZQsId5/q7iPdfWRNuZ8H3dwsrVtX/H7ampluuKHr7cykDz/ceNrxhQulCy8MfRoNDdKPfyw991yYH2r//cO9FhMmSCtWFF9WAD1Rpq0ebX1NLeXBev4w13Xrwl/+n/qUdO+9xe2rfTNTrsAZMCBMOX7wwWF5k02kKVOkb34z3KH9zjvSFVeEZ2XvsktolspkQnNUQ4P005923ccBACXW8wOipkYaMiTM0HrMMWG21u7+hd6+mWnmzK63nTkzTCU+Y0aYVvzjj0Olf/fdYXjsww9LZ58tDRwYpvd4663sb9etk37969B/cdFFITwAoNzcveiXpFslLZO0TtISSWd1tf2mm27qZfXOO+5bbukeBqe619e733ij+/r1he/rwguz+zn4YPfnntt4m8suy25TW5v93Pa6447stgsXutfUdFz/i1+4jx7dsbyXX+7e3Nz9/wYAqp6kZk+wLs716j03yr36qjRqVLbzWJLGjg3PfGhoyH8/H38cHkH6pz+FeZmk0Ex0+eVh/0cfHTqxu3LyydKXvxw6uevrpalTpdmzs+uHDw9PxXvppfCsildfDd8PGiRdcon03e9mp/0A0Gsw3XcpPfRQGKLa/pz79ZOuvFI699yOjxvNZdky6Ze/lK6+uvPpwZNwzjnStGkdv9tyyxBIZ5wh1daW7tgAKgoBUWrTpoWHAW1o1Kjw3Onhwwvb38qV0j77SK+9tvG6wYPDM7Efflg64YTO91FfH0ZaFWrIEOlXvwpXJH37Fv57AFWlKqfaqCrnnBOabTb0/PPSzjuHv85nzgzzL8Uq7XvvlXbdVbrvvnDlcP758XCQQsf42WdLX/tauGpZsSI0LW2os3D40Y9C5d+ZZcvCdB/Dh0t33cUjUQEkqvddQUihIr3jDun++8NIo4ULO992003DENlPPpGWLs32O0jS1ltL//lPdnnsWOnRR7PLtbVhRNLw4dKdd0ojRkh9+oTgmDJFuuyyrst57bXSWWeFz5mMdP31of+hM7vsIv32t9LhhxfWXAagKtDElIZMJtzdfPHFoSLvjhEjQtPVhAnZ7yZODDfVLV0alvv1CzfPjR0rfeUr4RGn99wTnjUxalR4ENFNN0mLFoXtTz01dKJvWNabbw5XLsuXx8tyxBHhCgdAj0JApO2DD0Ilf9NN8fX9+oVmqLZKX5L23FN67DHpyCPD0+faW7kyjHr661/jN76NGhXupm7v7belJ56QDjggzOkU0xYUF1wQmpra22EH6Z//7PI0AVQfAqJSPPJI9rGibc99iDXbrFsXmpJaWsJNb+2boKTQJHTNNeHz0qWhr+Ppp8Nf+PPmhcD56KPuNwltGBS1tdKzz0p77929/QGoWAREJVm/PlT8Awbk3va++8IVRFu/Q3v77y89+GDHIalXXRWuVA45JKwrViYTyrDddoQD0EMxiqmS9OmTXzhIoS9ByoZD+87k2bOlujpp3LgwOkqSbrklvJ90UjJlrakJN+kRDgASwhVEEtzDSKd33w3LX/1qCIX99tu4T0IKfRjvvRc+33lnuHkv3yAC0GtxBVGNFi7MhsM224T7KCZOjIeDlA0HSTruuPB0uh/8oPTlBIACEBBJmDUrvNfWhs8zZ4bZWNu+a7P55p3vY8ORSACQMgIiCSNGSEOHhruZ339f+s53wvdXXtmxw3rDYa677hpGNL35pnT77eUqLQDkhT6IJC1YIH3hC2Hk05lnStddF4a4jhsX1m+xRQiJtv/mtbXS6tWhAxsAcih3H0SZn/1Z4d58M1TYa9eG15o18fe1a8Ow0sMPDzelSeEehzFjQjgccID0xz+G7889N1xdLFwo/eQn4beLFoV7IBoaCAcAFYsriDbnnx+m7y7E6NHhprQPPwxXDosXh7uqn302zOEEAAniCiIte+4Z7ntYvz77XUNDGJVUVxeeKd2vX3j/979D30GfPuGK4NBDQzhsu20Y3ko4AOgB6KRuc+qp4dkNxx+f/W75cmnffaW//CU80+H++8MzpU85JazfaSfpxBPDjLADB4b5k7baKp3yA0DCCIj2dtghjCaaO1c66KDQz3DlleHK4JJLss9teP318P7gg+H5EHV14cphp53SKjkAJI4+iM5kMtKOO0pLlnT8/qKLpBdeCFcVUphk74EHQjMTAJQQk/VVig8/lAYNiq8bMiR7Y9sNN4TnOQBAiTHVRqUYOFAaNix8njhR+sMfwvTfv/tdeIa0FJqdCAcAPRSjmGLeeiv0PaxYEZa33loaPz68JGmffaQ33pCOPTa9MgJAidHE1MY9dDRffnl2biVJ2m230N+w887plQ0ARBNT+bmHfoTttw9NSLNmhfsbjj46zMY6fz7hAKBX6t1NTEuWhAf2PPVUWB40SPr+90NTUmfPggaAXiKRKwgzO8zMFprZYjM7P4l9lpR7mEhv551DOPTvL117behzuPRSwgFAVUq6Li46IMysr6SrJB0uaXdJJ5vZ7sXuN1Hjx4eb2czCDKpDh0pnnx0m3zvwwHDj21lnMXEegKpViro4iSuI0ZIWu/sb7r5W0m2Sjk5gv8kYP166+urscxkymXAPQ9++0vTp0kMPhfsaAKC6JV4XJxEQQyW93W55Set3HZhZo5k1mVlTJpNJ4LB5mjYt+/yFjgWSTjstvANAdahpq0dbX43t1uVVFxd0sGJ+3CpWw25UI7v7VElTpTDMNYHj5qf9E93aK2dIAUAyMu4+spN1edXFhUjiCmKJpGHtlhskLU1gv8lo/0zofL4HgOqUeF2cREA8L2kXM9vRzOoknSRpRgL7TcY552zcjGQWvgeAniPxurjoJiZ3z5jZBEkPSuor6Xp3n1fsfhNz1VXhfdq00NxUWxvCoe17AOgBSlEXM9UGAFQJptoAAFQEAgIAEEVAAACiCAgAQBQBAQCIIiAAAFEEBAAgioAAAEQREACAKAICABBFQAAAoggIAEAUAQEAiCIgAABRBAQAIIqAAABEERAAgCgCAgAQRUAAAKIICABAFAEBAIgiIAAAUQQEACCKgAAARBEQAIAoAgIAEFVUQJjZ8WY2z8zWm9nIpAoFAEhWd+rrYq8g5ko6VtLjRe4HAFBaBdfXNcUczd3nS5KZFbMbAECJdae+LiogCmFmjZIaJamurq5chwWAnqTGzJraLU9196klO1iuDcxslqRtIqsmufuf8z1Q60lMlaT6+nrPu4QAgDYZd++0/yCp+rpNzoBw94MK3SkAoPySrq8Z5goAiCp2mOsxZrZE0pcl3WdmDyZTLABAkrpTX5t7+bsD6uvrvbm5uezHBYBqZmYt7l5fruPRxAQAiCIgAABRBAQAIIqAAABEERAAgCgCAgAQRUAAAKIICABAFAEBAIgiIAAAUQQEACCKgAAARBEQAIAoAgIAEEVAAACiCAgAQBQBAQCIIiAAAFEEBAAgioAAAEQREACAKAICABBFQAAAoggIAEAUAQEAiCoqIMxsipktMLNXzOweMxucULkAAAnqTn1d7BXEQ5L2dPcRkhZJuqDI/QEASqPg+rqogHD3v7l7pnVxjqSGYvYHACiN7tTXSfZBnCnpgc5WmlmjmTWZWVMmk+lsMwBA52ra6tHWV2M399Nlfd3G3L3rDcxmSdomsmqSu/+5dZtJkkZKOtZz7VBSfX29Nzc359oMANCOmbW4e30X6xOtr2tyFcjdD8pR4NMlHSnpwHzCAQBQGknX1zkDIsfBDpP0c0n7uXtLMfsCAJROd+rrnE1MOQ64WFI/Se+1fjXH3c/N9TuamACgcLmamHL8tuD6uqgrCHffuZjfAwDKozv1NXdSAwCiCAgAQBQBAQCIIiAAAFEEBAAgioAAAEQREACAKAICABBFQAAAoggIAEAUAQEAiCIgAABRBAQAIIqAAABEERAAgCgCAgAQRUAAAKIICABAFAEBAIgiIAAAUQQEACCKgAAARBEQAIAoAgIAEEVAAACiCAgAQFRRAWFml5rZK2b2spn9zcy2TapgAIDkdKe+Nncv5oAD3f3D1s8/kLS7u5+b63f19fXe3Nzc7eMCQG9kZi3uXt/N3xZcXxd1BdF2sFb1krqfNgCAkulOfV1T7EHNbLKk0yR9IGn/LrZrlNQoSXV1dcUeFgB6oxoza2q3PNXdp+b743zr6//fPlcTk5nNkrRNZNUkd/9zu+0ukLSJu1+U66A0MQFA4XI1MSVdXxfVB7FBwbaXdJ+775lrWwICAApXTB/EBvvJq74udhTTLu0Wj5K0oJj9AQBKozv1dbF9EP9rZrtJWi/pLUk5RzABAFJRcH2dWBNTIcxsvaSPOlldIylTxuKUS088L86pevTE8+qN59Tf3ct2g3MqAdEVM2ty95FplyNpPfG8OKfq0RPPi3MqPabaAABEERAAgKhKDIi8b/qoMj3xvDin6tETz4tzKrGK64MAAFSGSryCAABUAAICABBVkQHRE58zYWZTzGxB63ndY2aD0y5TEszseDObZ2brzaxihud1h5kdZmYLzWyxmZ2fdnmSYGbXm9kKM5ubdlmSYmbDzGy2mc1v/bd3XtplKpaZbWJmz5nZ31vP6eK0yyRVaB9Ed58zUcnM7BBJj7h7xsx+KUnu/vOUi1U0M/uswp2Zf5T0E3dvyvGTimRmfSUtknSwpCWSnpd0sru/lmrBimRmYyStlvSnfOZJqwZmNkTSEHd/0cw2k/SCpG9U8/8rMzNJ9e6+2sxqJT0p6Tx3n5NmuSryCqInPmfC3f/m7m13SM6R1JBmeZLi7vPdfWHa5UjAaEmL3f0Nd18r6TZJR6dcpqK5++OS/pt2OZLk7svc/cXWz6skzZc0NN1SFceD1a2Lta2v1Ou9igwIKcxbbmZvSzpF0i/SLk/CzpT0QNqFQAdDJb3dbnmJqrzS6Q3MbAdJe0t6NuWiFM3M+prZy5JWSHrI3VM/p9QCwsxmmdncyOtoSXL3Se4+TNLNkiakVc5C5Dqn1m0mKcy1cnN6JS1MPufVA1jku9T/gkPnzGyApLsk/XCDVoeq5O6fuPvnFFoXRptZ6k2CRT9Rrrvc/aA8N71F0n2Scj6IKG25zsnMTpd0pKQDvRI7fzpRwP+rarZE0rB2yw2SlqZUFuTQ2k5/l6Sb3f3utMuTJHd/38welXSYpFQHF1RkE1NPfM6EmR0m6eeSjnL3lrTLg408L2kXM9vRzOoknSRpRsplQkRrh+51kua7+xVplycJZrZ128hGM+sv6SBVQL1XqaOY7pLUYd5yd38n3VIVx8wWS+on6b3Wr+ZU+8gsSTKzYyT9XtLWkt6X9LK7H5pqobrJzI6Q9FtJfSVd7+6T0y1R8czsVkljJW0l6d+SLnL361ItVJHMbB9JT0h6VaGOkKSJ7n5/eqUqjpmNkDRd4d9eH0m3u/sl6ZaqQgMCAJC+imxiAgCkj4AAAEQREACAKAICABBFQAAAoggIAEAUAQEAiPo/0hfMUs2FgjQAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "DONE! Completed 15 Monte Carlo cycles.\n" ] } ], "source": [ "# run an extra 10 to decorrelate a little futher\n", "sampler.run(10)" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": [ "sampler.storage.close()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "From here, you can either extend this to a longer trajectory for the fixed length TPS in the `AD_tps_1b_trajectory.ipynb` notebook, or go straight to flexible length TPS in the `AD_tps_2a_run_flex.ipynb` notebook." ] } ], "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": {}, "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 }, "widgets": { "application/vnd.jupyter.widget-state+json": { "state": {}, "version_major": 2, "version_minor": 0 } } }, "nbformat": 4, "nbformat_minor": 4 }