hi

README

@@ -1,3 +1,3 @@
 Q: What do you call it when Batman skips church?
 
-Find the answer somewhere in the codes!
+Find the answer somewhere in the codes! (Joke credit: Lee-Ping Wang)

Trp-Cage_AnalysisAndVisualization.ipynb

@@ -0,0 +1,123 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 86,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import simtk.openmm as mm\n",
+    "from simtk.openmm import app\n",
+    "from simtk import unit\n",
+    "import matplotlib.pyplot as plt\n",
+    "import mdtraj as md\n",
+    "import nglview as nv\n",
+    "%matplotlib inline"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "First load the trajectory:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 113,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "traj = md.load_dcd('trajectory.dcd', top='trpcage.pdb')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "And here are a couple of different ways of visualizing the trajectory:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 114,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "09607cdcf9894ab598527a9dd375baca",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": [
+       "NGLWidget(count=200)"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "view = nv.show_mdtraj(traj)\n",
+    "view.clear_representations()\n",
+    "view.add_representation(repr_type='cartoon', selection='protein', color='forestgreen')\n",
+    "view.add_representation('licorice', selection='protein')\n",
+    "## the next two commented lines may be useful for those who are doing ##\n",
+    "## the first option in the activity (adding an extra force)           ## \n",
+    "#view.add_representation('surface', selection='1.CA or 20.CA', surfaceType='vws')\n",
+    "#view.add_distance(atom_pair=[['1.CA', '20.CA']],label_color='black')\n",
+    "view"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 115,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "09607cdcf9894ab598527a9dd375baca",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": [
+       "NGLWidget(count=200, n_components=1)"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "view.clear_representations()\n",
+    "view.add_representation(repr_type='cartoon', selection='protein', color='forestgreen')\n",
+    "view.add_representation('licorice', selection='water')\n",
+    "view"
+   ]
+  }
+ ],
+ "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.6.5"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

run_openmm.py

@@ -1,4 +1,3 @@
-from __future__ import print_function
 from simtk.openmm import app
 import simtk.openmm as mm
 from simtk import unit
@@ -27,21 +26,14 @@
 # Use the CPU platform
 platform = mm.Platform.getPlatformByName('CPU')
 
-### If you want to add any forces to your System or modify any
-### of the existing forces, you should do it here - after the
-### System has been created, but before the Simulation is created.
-
-# Implemented option B, to create a merge conflict.
-def neutralizeCharge():
-    atoms = list(pdb.topology.atoms())
-    for f in system.getForces():
-        if isinstance(f, mm.NonbondedForce):
-            print("Found nonbonded force")
-            for i in range(f.getNumParticles()):
-                if atoms[i].residue.name != 'HOH':
-                    chg, sig, eps = f.getParticleParameters(i)
-                    f.setParticleParameters(i, 0.0, sig, eps)    
-
+def ZeroProteinCharges(a):
+    list = mm.System.getForces(a)
+    for i in list:
+        if isinstance(i, mm.NonbondedForce):
+            for l in range(2363):
+                a = mm.NonbondedForce.getParticleParameters(i,l)
+                mm.NonbondedForce.setParticleParameters(i,l,0,a[1], a[2])
+ZeroProteinCharges(system)
 # Create a Simulation object by putting together the objects above
 simulation = app.Simulation(pdb.topology, system, integrator, platform)
 
@@ -55,18 +47,29 @@ def neutralizeCharge():
 # Initialize the random velocities of the system from a Maxwell-Boltzmann distribution
 simulation.context.setVelocitiesToTemperature(300*unit.kelvin)
 
-# Add reporters to the simulation object, which do things at regular intervals
-# while the simulation is running.
-# This reporter creates a DCD trajectory file
+# *Incredibly Short* NPT equilibration (5 ps)
+# First, add reporters to the simulation object, which do things at regular intervals
+# while the simulation is running. This reporter prints information to the terminal.
+simulation.reporters.append(app.StateDataReporter(stdout, 100, step=True,
+    potentialEnergy=True, temperature=True, density=True, progress=True,
+    remainingTime=True, speed=True, totalSteps=2500, separator='\t'))
+# Next, run the equilibration simulation itself.
+print('Running Equilibration...')
+simulation.step(2500)
+
+# Production (40 ps)
+# Before doing anything, remember to clear the previous reporters. In the code below
+# the first reporter creates a DCD trajectory file. The second reporter is otherwise
+# very similar to the one above in the equilibration section. The only difference --
+# other than ouput frequency -- is that the totalSteps parameter is modified to be
+# the number of production steps + equilibration steps. Run the code to see why...
 # A: Christian Bale
+simulation.reporters.clear()
 simulation.reporters.append(app.DCDReporter('trajectory.dcd', 100))
-
-# This reporter prints information to the terminal
-simulation.reporters.append(app.StateDataReporter(stdout, 100, step=True, 
+simulation.reporters.append(app.StateDataReporter(stdout, 500, step=True, 
     potentialEnergy=True, temperature=True, density=True, progress=True, 
-    remainingTime=True, speed=True, totalSteps=10000, separator='\t'))
-
-# Run the simulation itself
+    remainingTime=True, speed=True, totalSteps=22500, separator='\t'))
+# Run the production simulation!
 print('Running Production...')
-simulation.step(10000)
+simulation.step(20000)
 print('Done!')