Activity3A_protein-MD.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Molecular Dynamics Simulation of a Protein\n",
    "\n",
    "In this activity, we will replicate the first published simulation of a protein from 1977.\n",
    "\n",
    "*McCammon, J. A.; Gelin, B. R.; Karplus, M. Nature 1977, 267 (5612), 585–590.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from openmm import app\n",
    "import openmm as mm\n",
    "from simtk import unit\n",
    "from sys import stdout\n",
    "import time as time\n",
    "import os "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulation Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#specify the names of the parameter/topology file and starting coordinate file\n",
    "topology_file = \n",
    "initial_coordinate_file = \n",
    "\n",
    "#specify the temperature at which you will run your simulation\n",
    "temperature = \n",
    "\n",
    "prmtop = app.AmberPrmtopFile(topology_file)\n",
    "inpcrd = app.AmberInpcrdFile(initial_coordinate_file)\n",
    "\n",
    "system = prmtop.createSystem(nonbondedMethod=app.NoCutoff, constraints=None)\n",
    "\n",
    "#set up the integrator specifying the temperature, friction coefficient, and timestep\n",
    "integrator = mm.LangevinIntegrator()\n",
    "\n",
    "platform = mm.Platform.getPlatformByName('Reference')\n",
    "simulation = app.Simulation(prmtop.topology, system, integrator, platform)\n",
    "simulation.context.setPositions(inpcrd.positions)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Minimization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "st = simulation.context.getState(getPositions=True,getEnergy=True)\n",
    "print(F\"Potential energy before minimization is {st.getPotentialEnergy()}\")\n",
    "\n",
    "print('Minimizing...')\n",
    "#set the number of iterations to 100\n",
    "simulation.minimizeEnergy()\n",
    "\n",
    "st = simulation.context.getState(getPositions=True,getEnergy=True)\n",
    "print(F\"Potential energy after minimization is {st.getPotentialEnergy()}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Equilibration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#set the initial velocities to your desired temperature\n",
    "simulation.context.setVelocitiesToTemperature()\n",
    "print('Equilibrating...')\n",
    "\n",
    "#equilibrate for 1 ps\n",
    "simulation.step()\n",
    "print('Done!')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Production Simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#specify the name of the trajectory file that will be generated\n",
    "trajectory_file = \n",
    "\n",
    "simulation.reporters.append(app.DCDReporter(trajectory_file, 1))\n",
    "\n",
    "#specify how often to print the energy\n",
    "print_freq = \n",
    "\n",
    "simulation.reporters.append(app.StateDataReporter(stdout, print_freq, step=True, time=True,\n",
    "    potentialEnergy=True, temperature=True, speed=True, separator='\\t'))\n",
    "\n",
    "tinit=time.time()\n",
    "\n",
    "print('Running Production...')\n",
    "\n",
    "#run simulation for 9 ps\n",
    "simulation.step()\n",
    "tfinal=time.time()\n",
    "\n",
    "print('Done!')\n",
    "print('Time required for simulation:', tfinal-tinit, 'seconds')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}