OpenMM examples

example_for_kevin.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "focal-glory",
   "metadata": {},
   "source": [
    "## Here's the API"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "advised-railway",
   "metadata": {},
   "source": [
    "http://docs.openmm.org/latest/api-python/index.html"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "sixth-columbia",
   "metadata": {},
   "source": [
    "## Here's the documentation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "unlike-latvia",
   "metadata": {},
   "source": [
    "http://docs.openmm.org/latest/userguide/index.html"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "plain-surprise",
   "metadata": {},
   "source": [
    "## Here's an example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "imposed-investing",
   "metadata": {},
   "outputs": [],
   "source": [
    "from simtk import openmm\n",
    "from simtk.openmm import app, unit\n",
    "from openmmforcefields.generators import SystemGenerator\n",
    "import openmmtools"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "experimental-mandate",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load PDB\n",
    "pdb = app.PDBFile(\"../../input/ala_vacuum_noh.pdb\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "prepared-water",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up modeller\n",
    "modeller = app.Modeller(pdb.topology, pdb.positions)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "infrared-killer",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create Forcefield object\n",
    "ffxmls = [\"amber14/protein.ff14SB.xml\", \"amber14/tip3p.xml\"]\n",
    "forcefield = app.ForceField(*ffxmls)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "thermal-taste",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[None, None, None]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Add hydrogens\n",
    "modeller.addHydrogens(forcefield, pH=7.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "fifty-slave",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add solvent\n",
    "modeller.addSolvent(forcefield, ionicStrength=150 * unit.millimolar, padding=1.0 * unit.nanometers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "subjective-favor",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save solvated PDB\n",
    "app.PDBFile.writeFile(modeller.topology, modeller.positions, open(\"solvated.pdb\", 'w'), keepIds=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "black-second",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Generate system and add barostat\n",
    "temperature = 300 * unit.kelvin\n",
    "system = forcefield.createSystem(modeller.topology, \n",
    "                                 nonbondedMethod=app.PME,\n",
    "                                 nonbondedCutoff=1 * unit.nanometer, \n",
    "                                 removeCMMotion=False,\n",
    "                                 constraints=app.HBonds, \n",
    "                                 hydrogenMass=4 * unit.amu)\n",
    "\n",
    "barostat = openmm.MonteCarloBarostat(1.0 * unit.atmosphere, temperature, 50)\n",
    "system.addForce(barostat)\n",
    "\n",
    "# If you already have prmtop and inpcrd files generated by tleap, you can load those in using:\n",
    "# http://docs.openmm.org/latest/api-python/generated/openmm.app.amberprmtopfile.AmberPrmtopFile.html#openmm.app.amberprmtopfile.AmberPrmtopFile\n",
    "# http://docs.openmm.org/latest/api-python/generated/openmm.app.amberinpcrdfile.AmberInpcrdFile.html#openmm.app.amberinpcrdfile.AmberInpcrdFile\n",
    "# and then call createSystem() (like above) on your loaded AmberPrmtopFile object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "brazilian-porcelain",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Serialize the system to an xml file. This is will be necessary for resuming a simulation\n",
    "with open('system.xml', 'w') as f:\n",
    "    f.write(openmm.XmlSerializer.serialize(system))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "adjusted-excuse",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up integrator\n",
    "integrator = openmmtools.integrators.LangevinIntegrator(temperature=300 * unit.kelvin, \n",
    "                                                        collision_rate=1 / unit.picoseconds, \n",
    "                                                        timestep=4 * unit.femtoseconds)\n",
    "\n",
    "# Note: we use the the openmmtools LangevinIntegrator rather than the OpenMM one because it is uses\n",
    "# BAOAB, wheras the OpenMM doesn't. Here are the OpenMM integrators: \n",
    "# http://docs.openmm.org/latest/api-python/library.html#integrators"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "extreme-austria",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configure platform\n",
    "platform = openmm.Platform.getPlatformByName('CUDA')\n",
    "platform.setPropertyDefaultValue('Precision', 'mixed')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "labeled-affect",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up simulation and minimize it\n",
    "simulation = app.Simulation(modeller.topology, system, integrator, platform)\n",
    "simulation.context.setPositions(modeller.positions)\n",
    "simulation.context.setPeriodicBoxVectors(*system.getDefaultPeriodicBoxVectors())\n",
    "simulation.context.setVelocitiesToTemperature(temperature)\n",
    "simulation.minimizeEnergy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "impaired-volleyball",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save energy minimized PDB\n",
    "positions_minimized = simulation.context.getState(getPositions=True, enforcePeriodicBox=False).getPositions()\n",
    "app.PDBFile.writeFile(modeller.topology, positions_minimized, open(\"minimized.pdb\", 'w'), keepIds=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "sonic-services",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up reporters for state data, checkpoint file, and trajectory, if desired\n",
    "simulation.reporters.append(app.StateDataReporter('ala.csv', 12500, time=True, step=True,\n",
    "        potentialEnergy=True, kineticEnergy=True, totalEnergy=True, volume=True, temperature=True))\n",
    "simulation.reporters.append(app.CheckpointReporter('ala.chk', 12500))\n",
    "simulation.reporters.append(app.DCDReporter('ala.dcd', 12500, enforcePeriodicBox=False))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "dental-teaching",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run the simulation\n",
    "n_steps = 100\n",
    "n_iterations = 10\n",
    "for _ in range(n_iterations):\n",
    "    simulation.step(n_steps)\n",
    "    \n",
    "    # Write the state every 100 steps. This is will be necessary for resuming a simulation\n",
    "    state = simulation.context.getState(getPositions=True, getVelocities=True)\n",
    "    with open('state.xml', 'w') as outfile:\n",
    "        state_xml = openmm.XmlSerializer.serialize(state)\n",
    "        outfile.write(state_xml)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dependent-certificate",
   "metadata": {},
   "source": [
    "## Example for resuming a simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "interior-construction",
   "metadata": {},
   "outputs": [],
   "source": [
    "from simtk import openmm\n",
    "from simtk.openmm import app\n",
    "import openmmtools"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "signal-supervision",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read in solvated PDB\n",
    "pdb = app.PDBFile(\"solvated.pdb\")\n",
    "\n",
    "# Load system\n",
    "with open('system.xml', 'r') as f:\n",
    "    system = openmm.XmlSerializer.deserialize(f.read())\n",
    "\n",
    "# Load state\n",
    "with open('state.xml', 'r') as infile:\n",
    "    state_xml = infile.read()\n",
    "    state = openmm.XmlSerializer.deserialize(state_xml)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "elementary-promise",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up integrator\n",
    "integrator = openmmtools.integrators.LangevinIntegrator(temperature=300 * unit.kelvin, \n",
    "                                                        collision_rate=1 / unit.picoseconds, \n",
    "                                                        timestep=4 * unit.femtoseconds)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "terminal-creek",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up platform\n",
    "platform = openmm.Platform.getPlatformByName('CUDA')\n",
    "platform.setPropertyDefaultValue('Precision', 'mixed')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "imposed-producer",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up simulation \n",
    "simulation = app.Simulation(pdb.topology, system, integrator, platform)\n",
    "simulation.context.setPeriodicBoxVectors(*state.getPeriodicBoxVectors())\n",
    "simulation.context.setPositions(state.getPositions())\n",
    "simulation.context.setVelocities(state.getVelocities())\n",
    "simulation.context.setTime(state.getTime())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "attached-settlement",
   "metadata": {},
   "source": [
    "... Everything else is the same as it was in the previous example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "based-lyric",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}