water3_fails.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:f61f20de6b2b22682e4319b130c3cb2a5fe4de85a3220dc3ced94a90236ccf9d"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from __future__ import print_function"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from simtk.openmm import app\n",
      "import simtk.openmm as mm\n",
      "from simtk import unit\n",
      "import mbpol"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Input system in pdb format"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# pdb = app.PDBFile(\"/home/zonca/Paesani/OpenMM/openmm/plugins/mbpol/python/water3.pdb\")\n",
      "pdb = app.PDBFile(\"/home/zonca/Paesani/OpenMM/openmm/plugins/mbpol/python/water3_fails.pdb\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Define the type of potential, first file defines all elements, only the water model is in the second xml file"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "forcefield = app.ForceField(\"/home/zonca/Paesani/OpenMM/openmm/plugins/mbpol/python/mbpol.xml\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Create the System, define an integrator, define the Simulation"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "system = forcefield.createSystem(pdb.topology, nonbondedMethod=app.CutoffNonPeriodic, nonBondedCutoff=1e3*unit.nanometer)\n",
      "integrator = mm.VerletIntegrator(0.00001*unit.femtoseconds)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "platform = mm.Platform.getPlatformByName('Reference')\n",
      "simulation = app.Simulation(pdb.topology, system, integrator, platform)\n",
      "simulation.context.setPositions(pdb.positions)\n",
      "simulation.context.computeVirtualSites()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Compute initial energy and forces with getState"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "state = simulation.context.getState(getForces=True, getEnergy=True)\n",
      "potential_energy = state.getPotentialEnergy()\n",
      "potential_energy.in_units_of(unit.kilocalorie_per_mole)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 7,
       "text": [
        "Quantity(value=-4.570117788168624, unit=kilocalorie/mole)"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "kilocalorie_per_mole_per_angstrom = unit.kilocalorie_per_mole/unit.angstrom\n",
      "for f in state.getForces():\n",
      "    print(f.in_units_of(kilocalorie_per_mole_per_angstrom))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(11.435647081051835, 17.894377982738252, 50.661575670007934) kcal/(A mol)\n",
        "(14.36783303635445, -21.50724325536891, -18.936391548521346) kcal/(A mol)\n",
        "(-16.231045906576846, 4.539194201166463, -28.323461551816848) kcal/(A mol)\n",
        "(-33.797609149637275, 5.65910479962967, -13.029774480663233) kcal/(A mol)\n",
        "(39.93871107810727, -3.496462044287423, -12.091543593613812) kcal/(A mol)\n",
        "(-57.93541548009736, 1.4975938780893399, 6.27053082023035) kcal/(A mol)\n",
        "(9.850156265979694, 2.0744000766332293, 7.147042103757614) kcal/(A mol)\n",
        "(-23.756732808508467, 0.6286318950020897, 2.061826844414644) kcal/(A mol)\n",
        "(-27.239800471924905, 25.609418415205067, 21.99425072451048) kcal/(A mol)\n",
        "(30.54710855724314, -34.66275158685176, -19.530109026050013) kcal/(A mol)\n",
        "(-4.733194160137288, 8.051472332675749, -7.191893598504327) kcal/(A mol)\n",
        "(0.6226899459526399, -4.8097093630309775, -15.417929066125732) kcal/(A mol)\n"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Local energy minimization"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from simtk.openmm import LocalEnergyMinimizer"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "LocalEnergyMinimizer.minimize(simulation.context)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "state = simulation.context.getState(getForces=True, getEnergy=True)\n",
      "potential_energy = state.getPotentialEnergy()\n",
      "potential_energy.in_units_of(unit.kilocalorie_per_mole)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 11,
       "text": [
        "Quantity(value=-17122.946663088784, unit=kilocalorie/mole)"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 11
    }
   ],
   "metadata": {}
  }
 ]
}