{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# HTMD FFEvaluate - Easy MM force-field evaluation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "HTMD provides the FFEvaluate module which allows the easy evaluation of molecular mechanics force-fields. This serves multiple purposes:\n",
    "\n",
    " * __Application development__: Many computational biology applications require the evaluation of energies between molecules. This allows faster development of applications based on those energies.\n",
    " * __Educational__: Most molecular dynamics codes are either closed source or written in a highly optimized manner making them hard to understand. The FFEvaluate code has been written in python and [numba](https://numba.pydata.org/) which makes it easy to read, improve and extend.\n",
    " * __Debugging__: FFEvaluate can provide a detailed breakdown of all energies and visualization of forces which helps debug problems in a MD simulation or parameterization.\n",
    " \n",
    "We would like to acknowledge [parmed](http://parmed.github.io/ParmEd/html/index.html), an awesome python library which FFEvaluate uses to parse forcefield parameters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quick example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from htmd.ffevaluation.ffevaluate import FFEvaluate, loadParameters, viewForces\n",
    "from htmd.molecule.molecule import Molecule\n",
    "from htmd.home import home\n",
    "from os.path import join\n",
    "\n",
    "datadir = join(home(dataDir='test-ffevaluate'), 'thrombin-ligand-amber')\n",
    "\n",
    "# Load the parameters\n",
    "prm = loadParameters(join(datadir, 'structure.prmtop')) \n",
    "# Load the topology\n",
    "mol = Molecule(join(datadir, 'structure.prmtop'))\n",
    "# Create a FFEvaluate object\n",
    "ffev = FFEvaluate(mol, prm)\n",
    "\n",
    "# Load some coordinates\n",
    "mol.read(join(datadir, 'structure.pdb'))\n",
    "# Calculate energy and forces\n",
    "energies, forces, atom_energies = ffev.calculate(mol.coords)  \n",
    "# Visualize the forces\n",
    "viewForces(mol, forces * 0.1)\n",
    "# Calculate interaction energy between the protein and the ligand\n",
    "ffev = FFEvaluate(mol, prm, betweensets=('protein', 'resname MOL'))\n",
    "energies = ffev.calculateEnergies(mol.coords)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Detailed explanation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from htmd.ffevaluation.ffevaluate import FFEvaluate, loadParameters\n",
    "from htmd.molecule.molecule import Molecule\n",
    "from htmd.home import home\n",
    "from os.path import join\n",
    "datadir = join(home(dataDir='test-ffevaluate'), 'test-thrombin-ligand-amber')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `loadParameters` function serves as a utility function for loading `AMBER` and `CHARMM` forcefield parameters. It is simply a wrapper for the `parmed` loading functions as it can sometimes get a bit complex to obtain parameters from some files, but if you already know how to use `parmed` you can use it directly to obtain a `ParameterSet` object of your parameter files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prm = loadParameters(join(datadir, 'structure.prmtop')) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have loaded the parameters, we load the `prmtop` topology in a `Molecule` object to obtain the bond, angle, dihedral, improper and atomtype information. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mol = Molecule(join(datadir, 'structure.prmtop'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The main class for the calculations is the `FFEvaluate` class. When instantiated, the class constructor converts and stores the parameters in a format friendly to the `numba python` calculations. This is only done once for each system & parameterization combination as it incurrs a small computational overhead. The energy calculations are then done with the class methods `calculate` and `calculateEnergies` which can be evaluated for multiple different configurations of the system without having to incurr the same overhead by re-converting the parameters at every call."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ffev = FFEvaluate(mol, prm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's read in some coordinates into our Molecule to evaluate the energies and forces."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mol.read(join(datadir, 'structure.pdb'))\n",
    "energies, forces, atom_energies = ffev.calculate(mol.coords)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see, FFEvaluate returns:\n",
    "\n",
    "* The energy of the system broken down in bonded, angle, dihedral, improper, VdW and electrostatic terms in a (6, nframes) shaped array where the rows correspond to `0: bond 1: lennard-jones 2: electrostatic 3: angle 4: dihedral 5: improper`. You can sum over axis 0 to get the total energy.\n",
    "* The forces of all atoms in a `(natoms, 3, nframes)` shaped array\n",
    "* The energies of each individual atom in a `(natoms, 6, nframes)` shaped array. These energies are not physically meaningful as they are approximate atom energies calculated as the sum of all non-bonded energies on that atom plus the fraction of the bonded energy terms (i.e. for bonds we divide equally the bond energy between the two atoms, for angles between all three atoms and for dihedrals/impropers between all 4 atoms). However they can be useful for identifying issues in a parameterization or high-energy areas such as collisions.\n",
    "\n",
    "where `natoms` is the number of atoms in `mol` and `nframe` is the number of frames in the coordinates of `mol` (in this case only one).\n",
    "\n",
    "Warning: The first call to `calculate` in a python session can take some time as `numba` will compile the functions. Subsequent calls in the same session will be much faster."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`FFEvaluate` helps us visualize the forces using VMD. This can be very helpful for debugging cases where the simulation crashes due to atoms flying away because of too large forces. In case the vectors of the forces show up too large in VMD we can simply scale the forces down as much as we want for visualization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "viewForces(mol, forces * 0.1)\n",
    "# On systems with 1000s of atoms it can take a while to render all force vectors in VMD"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To just obtain energies in a nicely formatted fashion, the class also provides the `calculateEnergies` utility method which performs the same calculations but returns only the energies in a python dictionary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ffev.calculateEnergies(mol.coords)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lastly, we are often interested specifically in (non-bonded) interaction energies between two molecules or generally between two sets of atoms. FFEvaluate allows us to calculate these energies using the `betweensets` argument. This argument accepts a list or tuple of two `atomselect` strings which define the two atom sets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ffev = FFEvaluate(mol, prm, betweensets=('protein', 'resname MOL'))\n",
    "energies = ffev.calculateEnergies(mol.coords)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Disclaimers\n",
    "\n",
    "While FFEvaluate provides an easy and fast way of obtaining MM energies and forces, the user should be aware of some pitfalls.\n",
    "\n",
    " * The code is still in beta. There are some less-used or older force-field features which are not yet supported (i.e. Urey-Bradley terms, PRMTOP EXCLUDED_ATOMS_LIST)\n",
    " * It will never be as fast as professional MD codes. The purpose of `FFEvaluate` is to make an easily usable and understandable MM evaluation code. Therefore on systems with large amount of non-bonded interations such as solvated systems with many waters the code __will__ take a long time to evaluate the energies and forces.\n",
    " * Speed improvements such as cell-lists for non-bonded interations and particle-mesh Ewald are not implemented (but could be in the future). FFEvaluate will calculate all non-bonded interactions explicitly.\n",
    " * Implicit solvent calculations are still missing but might be added in the future.\n",
    " \n",
    "If you are interested in contributing code to the project, or want to report a bug, feel free to chat us up on the [HTMD github issues](https://github.com/acellera/htmd/issues/)!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}