Foyer_01_SMARTS_and_Overrides.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Foyer Tutorial 01: SMARTS and Overrides\n",
    "\n",
    "This tutorial focuses on the basic usage of SMARTS to define chemical context and the use of overrides to set rule precedence, and demonstrates how to use a force field file to atomtype a system."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As shown in the prior tutorial, atom-type usage (i.e., chemical context) is defined via a SMARTS string contained in the `def` tag.  Below is the AtomTypes subsection of the XML file for a forcefield that describes alkanes using the OPLS all-atom force field (as detailed in [Jorgensen, et al.](http://dx.doi.org/10.1021/ja9621760)). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "! sed -n 2,18p utils/OPLSaa_alkanes.xml"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Defining SMARTS\n",
    "Focusing first on atom type `opls_140`, the SMARTS string, `def=\"H[C;X4]\"`, states that this atom-type applies when:\n",
    "- The element is hydrogen, i.e., `H`\n",
    "- When that hydrogen is connected to a carbon that has 4 neighbors, i.e., `[C;X4]`\n",
    "\n",
    "Similarly, for atom type `opls_138`, the SMARTS string, `def=\"[C;X4](H)(H)(H)H\"`, states that this atom-type applies when:\n",
    "- The element is carbon, with 4 neighbors, i.e., `[C;X4]`\n",
    "- 4 of those neighbors are hydrogens, i.e., `(H)(H)(H)H`\n",
    "\n",
    "For atom type `opls_136`, the SMARTS string, `def=\"[C;X4](H)H\"`, states that this atom-type applies when:\n",
    "- The element is carbon, with 4 neighbors, i.e., `[C;X4]`\n",
    "- At least 2 of those neighbors are hydrogens, i.e., `(H)H`\n",
    "\n",
    "For atom type `opls_135`, the SMARTS string, `def=\"[C;X4](H)(H)H\"`, states that this atom-type applies when:\n",
    "- The element is carbon, with 4 neighbors, i.e., `[C;X4]`\n",
    "- At least 3 of those neighbors are hydrogens, i.e., `(H)(H)H`\n",
    "\n",
    "### Atomtyping using SMARTS\n",
    "Let us now consider using these rules to atom-type the carbon in methane (i.e., CH4).\n",
    "\n",
    "- `opls_138` would obviously evaluate to `True`, as it is defined for a carbon, with 4 hydrogen neighbors. \n",
    "- `opls_135` and `opls_136` would also evaulate to `True`.  In the case of opls_135, our definition only states that at least 3 of the neighbors are hydrogen, and have not made any specific claims about the identity of the 4th neighbor; similarly, for opls_136, we have only stated that 2 neighbors must be hydrogen and not specified the identity of the other 2 neighbors. \n",
    "\n",
    "This is an important feature of Foyer to take note of.  Foyer will evaluate all the rules in a forcefield file, rather than just stopping at the first rule that evaluates to `True`. This allows rules to be defined in any order.  Furthermore, Foyer iterates over the rules, which allows recursive definitions of usage, i.e., referring to specific atom-types in the SMARTS string. \n",
    "\n",
    "### Overrides\n",
    "We will discuss two ways to address this. One approach is to employ `overrides`.  `overrides` provide a means to dictate rule precedence (i.e., which rules are more specific than others).  In the force field file above, `opls_138` has defined: `overrides=\"opls_135,opls_136\"`.  That is, if `opls_138` evaluates to `True`, then it takes precedence over `opls_135` and `opls_136`, even if they evaluate to `True`. \n",
    "\n",
    "Similarly, if opls_136 evaluates to `True`, it `overrides=\"opls_136\"`, thus taking precedence. \n",
    "\n",
    "`overrides` are especially useful if the chemical context of two different atom-types are effectively the same.  E.g., the terminal methyl group in an alkane has the same first neighbor environment as the methyl group in toluene, however the parameters (and thus atom-type) are different. Thus `overrides` can be used to state that if the more specific toluene rule evaluates to `True` it should take precedences over the more general alkane rule (as shown below):\n",
    "\n",
    "<img src=\"utils/ch3-toluene.png\" alt=\"Drawing\" style=\"width: 700px;\"/>\n",
    "\n",
    "### Better SMARTS definitions\n",
    "In many cases, `overrides` can be avoided by simply providing more specific SMARTS strings.  For example, the rules for `opls_135`, `opls_136`, and `opls_138` above can be made more specific by stating the identify of the other neighbors besides carbon and thus eliminate the need for `overrides`, as shown below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "! sed -n 2,18p utils/OPLSaa_alkanes2.xml"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that SMARTS can be used to define more than just the immediate local pattern. Recall the definition for `opls_140`, `def='H[C;X4]'`. This definition states not only that the element is hydrogen and is bonded to a carbon, but that the carbon it is bonded to has 4 total bonds. If necessary, the identify of those bonded neighbors could  be specified. \n",
    "\n",
    "Additionally, we note that since this force field only has a single hydrogen atom-type, the definition could actually be made less specific, i.e., `def='H'`, and still produce the correct output.  \n",
    "\n",
    "Thus, it is important to keep in mind that there are multiple valid SMARTS definitions that can be provided for a given atom-type; the specificity of the SMARTS definition and whether or not `overrides` are necessary will depend on the chemical context of the parameters themselves, the number of atom-types included in a forcefield, the intended usage of the forcefield, and the personal preferences of the individual(s) defining the forcefield.   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Applying a force field to atomtype a system\n",
    "With a force field file defined, we will apply it to an alkane system.  Here, we will use the routines previously discussed in the mBuild tutorial to construct a hexane, then perform atom-typing as part of saving our configuration to a desired file format.\n",
    "\n",
    "Note, the following command sequence may provide a long list of Deprecation Warnings (which will be suppressed in the next mBuild release)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import mbuild as mb\n",
    "\n",

    "hexane = mb.recipes.Alkane(6)\n",
    "\n",
    "hexane.save('hexane.top', forcefield_files='utils/OPLSaa_alkanes.xml', overwrite=True)\n",
    "hexane.save('hexane.lammps', forcefield_files='utils/OPLSaa_alkanes.xml', overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can easily view the output files. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# view the gromacs formatted file\n",
    "!cat hexane.top"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# view the lammps formatted file\n",
    "!cat hexane.lammps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Forcefield files allow a DOI for the original parameter source to be defined, thus enabling Foyer to output a bibtex formatted references file, that clearly states in the notes which atom-types are associated with a given reference. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hexane.save('hexane.top', forcefield_files='utils/OPLSaa_alkanes.xml', foyer_kwargs={'references_file':'hexane.bib'}, overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cat hexane.bib"
   ]
  }
 ],
 "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.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}