Activity2A_MD-simulations-of-alkanes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "9acb02e9",
   "metadata": {},
   "source": [
    "# Alkane Simulation with OpenMM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c3c6360d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# OpenMM imports\n",
    "import openmm as mm\n",
    "from openmm import app, unit"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a492e963",
   "metadata": {},
   "source": [
    "## Simulation Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a387c91e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read in our structure and force field\n",
    "pdb_file = \"data/ethane.pdb\"\n",
    "ff_file = \"data/ethane.gaff2.xml\"\n",
    "\n",
    "pdb = app.PDBFile(pdb_file)\n",
    "forcefield = app.ForceField(ff_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e819391",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create OpenMM representation of the system\n",
    "system = forcefield.createSystem(pdb.topology, nonbondedMethod=app.NoCutoff, constraints=app.HBonds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "377794a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Set up the integrator\n",
    "# run gas phase simulation using a Langevin thermostat (integrator)\n",
    "# at 298.15 K\n",
    "# with friction constant (or coupling constant) of 5.0 ps^-1\n",
    "# with 2 fs step size \n",
    "\n",
    "# Intergrator parameters\n",
    "temperature = 298.15*unit.kelvin # temperature in Kelvin\n",
    "friction_constant = 5.0/unit.picoseconds # coupling constant in picoseconds^-1 \n",
    "stepsize = 2.0*unit.femtoseconds # in femtoseconds\n",
    "\n",
    "# Set temperature, coupling constant, and time step\n",
    "integrator = mm.LangevinIntegrator(temperature, friction_constant, stepsize)\n",
    "\n",
    "integrator.setConstraintTolerance(1e-5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff6f185e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize the simulation\n",
    "platform = mm.Platform.getPlatformByName('Reference')\n",
    "simulation = app.Simulation(pdb.topology, system, integrator, platform)\n",
    "simulation.context.setPositions(pdb.positions)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "320e6fb6",
   "metadata": {},
   "source": [
    "# Energy Minimization\n",
    "\n",
    "The minimization will move the atoms slightly to get them in a more energetically favorable position."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c36fbdd",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Minimizing...')\n",
    "\n",
    "st_initial = simulation.context.getState(getPositions=True,getEnergy=True)\n",
    "print(F\"Potential energy before minimization is {st_initial.getPotentialEnergy()}\")\n",
    "\n",
    "simulation.minimizeEnergy(maxIterations=100)\n",
    "\n",
    "st_minimized = simulation.context.getState(getPositions=True,getEnergy=True)\n",
    "print(F\"Potential energy after minimization is {st_minimized.getPotentialEnergy()}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff3bbf1a",
   "metadata": {},
   "source": [
    "# Equilibration\n",
    "\n",
    "We'll run the simulation for a few time steps to make sure the system is equilibrated. We usually do this at the start of a simulation to make sure the data we collect later isn't influenced by the starting position of the atoms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "804044f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sys import stdout\n",
    "\n",
    "print('Equilibrating...')\n",
    "\n",
    "# Sets up what information OpenMM will report as the simulations runs\n",
    "simulation.reporters.append(app.StateDataReporter(stdout, 100, step=True, \n",
    "    potentialEnergy=True, temperature=True, separator='\\t'))\n",
    "\n",
    "# Set initial velocities\n",
    "simulation.context.setVelocitiesToTemperature(150.0*unit.kelvin)\n",
    "\n",
    "simulation.step(2500)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7ad08a8",
   "metadata": {},
   "source": [
    "# Production\n",
    "\n",
    "This is the portion of the simulation we will collect data from and analyze."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4521dba4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "\n",
    "print('Running Production...')\n",
    "\n",
    "# Begin timer\n",
    "tinit=time.time()\n",
    "\n",
    "# Clear simulation reporters\n",
    "simulation.reporters.clear()\n",
    "\n",
    "# Reinitialize simulation reporters. \n",
    "# We do this because we want different information printed \n",
    "# from the production run than the equilibration run.\n",
    "# Output basic simulation information below every 250000 steps (How many picoseconds is?) \n",
    "\n",
    "simulation.reporters.append(app.StateDataReporter(stdout, 250000, \n",
    "    step=True, time=True, potentialEnergy=True, temperature=True, \n",
    "    speed=True, separator=','))\n",
    "\n",
    "# write out a trajectory (i.e., coordinates vs. time) to a DCD\n",
    "# file every 100 steps = 0.2 ps\n",
    "simulation.reporters.append(app.DCDReporter(F'ethane_sim_{temperature}.dcd', 100))\n",
    "\n",
    "# run the simulation for 1.0x10^7 steps - 20 ns\n",
    "simulation.step(10000000)\n",
    "\n",
    "# End timer\n",
    "tfinal=time.time()\n",
    "print('Done!')\n",
    "print('Time required for simulation:', tfinal-tinit, 'seconds')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc79dd4b",
   "metadata": {},
   "source": [
    "### Summary of this worked example\n",
    "\n",
    "This workflow contains all of the necessary commands to use OpenMM to carry out a gas phase simulation of a single ethane molecule. \n",
    "\n",
    "As written, it requires the files:\n",
    "* ethane.gaff2.xml (force field information)\n",
    "* ethane.pdb (starting coordinates and topology)\n",
    "to be imported to the notebook\n",
    "\n",
    "Briefly, this notebook carried out: \n",
    "1. Up to 100 steps of energy minimization using the L-BFGS algorithm.\n",
    "2. A 5.0 ps MD simulation to bring the ethane molecule to an equilibrium temperature of 298 K (25 °C) in which output is printed every 0.2 ps (100 steps).\n",
    "3.\tA 20 ns MD simulation at 298 K in which output is printed every 500 ps and structures are saved every 0.2 ps into a file called ethane_sim.dcd.  (This trajectory file is written in a commonly used binary format called DCD.)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64061cc7",
   "metadata": {},
   "source": [
    "## Your Turn! - Simulation of Butane\n",
    "\n",
    "Create a copy of this notebook of the ethane simulation. \n",
    "Modify this notebook to:\n",
    "1. Read in the files **butane.gaff2.xml** and **butane.pdb**\n",
    "2. Carry out a 10 ps MD simulation to bring the butane molecule to an equilibrium temperature of 298 K in which output is printed every 0.5 ps (Leave the minimization portion beforehand unchanged.)\n",
    "3. Carry out a 40 ns MD simulation at 298 K in which output is printed every 1 ns and structures are (still) saved every 0.2 ps into a file called butane_sim.dcd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9429506",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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": 5
}