foyer-overview.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MoSDeF - A Molecular Simulation Design Framework\n",
    "\n",
    "## Foyer Overview\n",
    "Foyer is designed to be a force field agnostic tool for performing atom-typing that can output to various file formats. That is, the goal is to create a single, flexible tool that can work with a range of force fields and simulation engines, rather than creating software specific to a single forcefield/simulation engine. \n",
    "\n",
    "A key goal of Foyer is to allow force field usage to be unambiguous and the abilication of force fields to be reproducible.  Foyer employs a force field file format that is:\n",
    "- simultaneously human and machine readable\n",
    "- can be automatically validated\n",
    "- easy to disseminate \n",
    "- easy to evolve and expand.\n",
    "\n",
    "### File Format\n",
    "Foyer force fields are defined within an XML file that contains both the 'rules' required for atomtyping as well as the force field parameters within a single file. \n",
    "\n",
    "The Foyer XML format is an extension of the [OpenMM forcefield XML format](http://docs.openmm.org/7.0.0/userguide/application.html#creating-force-fields). The only differences reside in the `AtomTypes` section, where several additional attributes are available (which we will examine in a moment) that allow for atomtyping.\n",
    "\n",
    "The `AtomTypes` section of the Foyer XML is similar to that used for OpenMM forcefield XMLs; however, each `Type` in Foyer XML supports four additional attributes not found in OpenMM:\n",
    "* `def` - SMARTS string describing the chemical substructure of this atom type, as discussed later (Follow [this link](https://github.com/mosdef-hub/foyer/blob/master/docs/smarts.md) for more on SMARTS-based atomtyping using Foyer.)\n",
    "* `desc` - Brief description of the atomtype\n",
    "* `doi` - DOI reference for parameters associated with this atomtype\n",
    "* `overrides` - One or more atomtypes to 'override', providing precedence to this atomtype (discussed later)\n",
    "\n",
    "\n",
    "Since all the usage rules are encoded in a text file outside of the software itself, rules can be evolved and corrected without changing the software itself; it also makes it easier to develop and distribute new force fields and subsets specific to a project/manuscript. \n",
    "\n",
    "\n",
    "Let's quickly look at a Foyer XML file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cat utils/oplsaa-alcohol.xml"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Defining Chemical Context (i.e., \"Rules\") using SMARTS\n",
    "\n",
    "Classical force fields are typically able to achieve high accuracy by creating sets of highly specific fitting parameters (i.e., atom types), in which each atom type describes an interaction site within a different chemical context. The chemical context is typically defined by the bonded environment of an interaction site (e.g., the number of bonds and the identity of the bonded neighbors) and may also consider, among other factors, the bonded environment of the neighbors, and/or the specific molecule/structure within which the interaction site is included.  Chemical context effectively defines the 'rules' for when an atomtype should apply. \n",
    "\n",
    "Foyer relies upon using [SMARTS](http://www.daylight.com/dayhtml/doc/theory/theory.smarts.html) to define the chemical context of an atom-type. SMARTS is a language for describing chemical structures and substructures.  In Foyer we tend focus on substructures, rather than trying to describe the entire molecule in a single string, given that the goal is to uniquely describe differentiate each interaction site/atom-type in the system.\n",
    "\n",
    "\n",
    "\n",
    "### Example, Methane: \n",
    "Consider defining the chemical context of the atoms (i.e., interactions sites) in a methane molecule ($CH_4$) using SMARTS.  The first thing to note is that there are multiple valid SMARTS that could describe the atom.  \n",
    "\n",
    "#### Carbon:\n",
    "- For example, the simplest way to describe the carbon atom in methane would be:\n",
    " - ```C``` \n",
    "- While in some cases such a simple definition may be sufficient, typically we wish to also define the number of bond, e.g.:\n",
    " - ```[C;X4]```\n",
    " - Note, ```;``` represents ```AND```, thus stating that the atom is a carbon AND it has 4 bonds. \n",
    "- Often it is useful to further define the identity of those bonds. In this case, we would typically describe methane as: \n",
    " - ```[C;X4](H)(H)(H)(H)``` \n",
    "\n",
    "#### Hydrogen:\n",
    "- Since hydrogen can only form a single bond, we could simply define it as:\n",
    " - ```H```\n",
    " - however, this may be insufficient as the parameters will likely be different based on what it is bonded to. \n",
    "- Additional specificity can be encoded by providing information about the neighbors, e.g.:\n",
    " - ```H(C)``` \n",
    " - this specifies this hydrogen must be bonded to a carbon.\n",
    "- This could be made further specific by further specifyin that this is a hydrogen that is bonded to a carbon with four bonds. \n",
    " -  ```H[C;X4]```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1:\n",
    "\n",
    "#### Setting up a force field\n",
    "Modify the \"mff_test1.xml\" file (in the \"mff_test1\" directory) defining SMARTS strings that appropriately describe the atom types for simple linear alkanes based on our fictious forcefield \"mosdef_ff\". \n",
    "\n",
    "Here, the atom types for \"mosdef_ff\" are:\n",
    "\n",
    "* `mff_0`: carbon in methane (CH_4)\n",
    "* `mff_1`: carbon in a terminal methyl group of a linear chain (-CH_3)\n",
    "* `mff_2`: carbon in a middle methyl group of a linear chain (-CH_2-)\n",
    "* `mff_3`: a generic hydrogen atom\n",
    "\n",
    "<img src=\"./utils/mff.png\" alt=\"Drawing\" style=\"width: 700px;\"/>\n",
    "\n",
    "\n",
    "After defining the appropriate SMARTS strings, run the test suite (below) to ensure SMARTS are correctly encoded.  This uses pytest to compare the output of your forcefield to the atom types listed in known mol2 files for methane, ethane, and propane. \n",
    "\n",
    "Note, don't worry about the bonds/angles/dihedrals listed in the file.  Those are just to avoid warnings about an underspecified force field (Foyer will provide warnings if connected pairs/triples/quartets are identified, but not parameterized as bond/angles/dihedrals)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%cd mff_test1/\n",
    "!py.test -v --tb=line\n",
    "%cd .."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Working with coarse-grained and united atom forcefields\n",
    "\n",
    "Foyer allows non-atomistic types to be defined within SMARTS, allowing coarse-grained and united atom forcefields to be handled as well. Non-elemental species can easily be defined by pre-pending the name of custom \"element\" with an underscore.\n",
    "\n",
    "For example, let us consider defining the SMARTS for a simple homopolymer composed of \"A\" beads, where terminal groups are of atom type \"cg_term\" and middle groups are \"cg_mid\". Here, the key differentiation between atom types is the number of bonds (1 for cg_term, 2 for cg_mid): \n",
    "\n",
    "- `cg_term` : `[_A;X1]`\n",
    "- `cg_mid` : `[_A;X2]`\n",
    "\n",
    "Similarly, the following lines could be used to describe beads representing $-CH_2-$ groups in a polymer using the TraPPE forcefield. \n",
    "\n",
    "`  <Type name=\"CH2_sp3\" class=\"CH2\" element=\"_CH2\" mass=\"14.02700\" \n",
    "   def=\"[_CH2;X2]([_CH3,_CH2])[_CH3,_CH2]\" \n",
    "   desc=\"Alkane CH2, united atom\" doi=\"10.1021/jp972543+\"/>`\n",
    "  \n",
    "Here, the SMARTS definition `[_CH2;X2]([_CH3,_CH2])[_CH3,_CH2]` states that for atom-type `CH2_sp3`\n",
    "\n",
    "- our bead is `_CH2` with 2 bonded neighbors, i.e., `[_CH2;X2]`\n",
    "- those neighbors can be either `_CH2` or `_CH3`, since, i.e., `([_CH3,_CH2])[_CH3,_CH2]` (note, `,`  represents an \"or\" statement).\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using  `overrides` to set rule precendence\n",
    "\n",
    "Force fields often contain atom types that ostensively have matching chemical context (at least in term of local bonded environment), but require different parameters.  For example, consider defining the a force field for both alkenes and benzene in a single file for the OPLS force field (note, for clarity only only 2 atom types are included):\n",
    "\n",
    "`\n",
    "<ForceField>\n",
    "  <AtomTypes>\n",
    "    <Type name=\"opls_142\" class=\"CM\" element=\"C\" mass=\"12.01100\"\\\\\n",
    "        def=\"[C;X3](C)(C)H\" desc=\"alkene C (RH-C=)\"/>\n",
    "    <Type name=\"opls_145\" class=\"CA\" element=\"C\" mass=\"12.01100\"\\\\\n",
    "        def=\"[C;X3;r6]1[C;X3;r6][C;X3;r6][C;X3;r6][C;X3;r6][C;X3;r6]1\"\\\\\n",
    "        overrides=\"opls_142\"/>\n",
    "  </AtomTypes>\n",
    "</ForceField>\n",
    "`\n",
    "\n",
    "When atom-typing a benzene molecule, the carbon atoms in the ring will match the SMARTS patterns for both `opls_142` (an alkene carbon) and `opls_145` (a benzene carbon). \n",
    "\n",
    "\n",
    "`Foyer` allows rule precedence to be explicitly stated via the use of the `overrides` attribute added to the XML file format. This allows atom type usage rules to be encoded in any order within the file, eliminating incorrectly placed rule order as a source of error. This approach also means that all rules are evaluated, helping to better ensure completeness and remove ambiguity.  Providing the `overrides` indicates that if the `opls_145` pattern matches, it will supersede `opls_142`.\n",
    "\n",
    "\n",
    "`Foyer` iteratively evaluates all rules on all interaction sites in the system, maintaining for each interaction site a \"whitelist\" consisting of rules that evaluate to `True` and a \"blacklist\" consisting  of rules that have been superseded by another rule (i.e., those that appear in the `overrides` attribute). The set difference between the white- and blacklists of an interaction site yields the correct atom type if the force field is implemented correctly. Thus, in this example the difference between the whitelist (containing `opls_142` and `opls_145`) and blacklist (containing only `opls_142`) would be `opls_145`.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 2\n",
    "Let us modify the prior fictious \"mosdef_ff\" to include rules that will require overrides. \n",
    "\n",
    "Here we will add a new atom type `mff_4` which is a generic carbon that should be used when more specific matches are not found. \n",
    "\n",
    "Here, the atom types for \"mosdef_ff\" are:\n",
    "\n",
    "* `mff_0`: carbon in methane (CH_4)\n",
    "* `mff_1`: carbon in a terminal methyl group of a linear chain (-CH_3)\n",
    "* `mff_2`: carbon in a middle methyl group of a linear chain (-CH_2-)\n",
    "* `mff_3`: a hydrogen atom \n",
    "* `mff_4`: a generic carbon\n",
    "\n",
    "Here, we have an additional test case that is simply a linear chain of connected carbons (\"linear.mol2\"). \n",
    "\n",
    "Make these edits to \"mff_test2.xml\" (in the \"mff_test2\" directory).  Here is an example of encoding an override.\n",
    "\n",
    "```\n",
    "<Type name=\"opls_138\" def=\"[C;X4](H)(H)(H)H\"\n",
    "    class=\"CT\" element=\"C\" mass=\"12.01100\" desc=\"alkane CH4\"\n",
    "    doi=\"10.1021/ja9621760\" overrides=\"opls_135,opls_136\"/>\n",
    "```\n",
    "\n",
    "When complete, run pytest below, as before.  \n",
    "#### Try to break the force field\n",
    "Assuming you did everything correct, modify the override statement of mff_4 such that you will cause the methane test to fail and examine the \"FAILURES\" reported by Foyer. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%cd mff_test2/\n",
    "!py.test -v --tb=line\n",
    "%cd .."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Force field application with `Foyer`. Generating data files.\n",
    "\n",
    "If we wanted to actually run a simulation of any of these systems we've built with mBuild, we would need to apply a force field and write the necessary data files. mBuild handles all of this through a single `save` command, where we can pass as arguments the name of the force field to apply (which uses `Foyer` under the hood) and the name of the file to create, which will be formatted based on the extension.\n",
    "\n",
    "\n",
    "First, let's consider how we would write to Gromacs `TOP` and `GRO` formats.\n",
    "\n",
    "The `GRO` format contains no force field information, so we do not have to pass a force field file to `save` when writing to this format. Let's start by using the \"mosdef_ff\" force field we set up previously to parameterize octane. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import mbuild as mb\n",
    "from mbuild.lib.atoms import H\n",
    "from mbuild.lib.moieties import CH2\n",
    "\n",
    "\n",
    "class Alkane(mb.Compound):\n",
    "    def __init__(self, chain_length=1):\n",
    "        super(Alkane, self).__init__()\n",
    "        hydrogen = H()\n",
    "        last_monomer = CH2()\n",
    "        # top capping CH2 -> CH3\n",
    "        mb.force_overlap(move_this=hydrogen,\n",
    "                         from_positions=hydrogen['up'],\n",
    "                         to_positions=last_monomer['up'])\n",
    "        # Add to our `self` compound\n",
    "        self.add(hydrogen)\n",
    "        self.add(last_monomer)\n",
    "        # loop over n = chain_length CH2's to add to Alkane\n",
    "        for _ in range(chain_length-1):\n",
    "            current_monomer = CH2()\n",
    "            mb.force_overlap(move_this=current_monomer,\n",
    "                             from_positions=current_monomer['up'],\n",
    "                             to_positions=last_monomer['down'])\n",
    "            self.add(current_monomer)\n",
    "            last_monomer=current_monomer\n",
    "        # bottom cap\n",
    "        hydrogen = H()\n",
    "        mb.force_overlap(move_this=hydrogen,\n",
    "                         from_positions=hydrogen['up'],\n",
    "                         to_positions=last_monomer['down'])\n",
    "        self.add(hydrogen)\n",
    "\n",
    "octane = Alkane(chain_length=8)\n",
    "\n",
    "\n",
    "octane.save('octane_test.gro', overwrite=True)\n",
    "octane.save('octane_test.top', forcefield_files='mff_test2/mff_test2.xml', overwrite=True, combining_rule='geometric')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can easily inspect the output:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cat octane_test.top"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we wish to setup an entire system, we can also specify a `residues` argument. In this case, we are saying to treat every `Compound` with the name `Octane` as a separate residue. This helps to speed up the calculation by avoiding duplicating the work, since the box is effectively filled with duplicates of the same molecule.  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With this force field XML file, Foyer will use the SMARTS strings to atom-type our system and will then apply the proper force field parameters. We'll execute the `save` method again, this time passing through our force field file and changing the desired file format from GRO to TOP. Additionally, as OPLS uses geometric mixing rules as opposed to Lorentz-Berthelot, we can feed this to `save` as well.\n",
    "\n",
    "**Note:** The warning message about unparameterized impropers can be safely ignored, as OPLS does not include any impropers for our system. By default, Foyer will warn the user if improper parameters are not specified for all possible impropers and will exit with an error if bond, angle, or dihedral parameters are not specified for all possible bonds, angles, dihedrals. (This behavior may be overridden if desired.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "octane_box = mb.fill_box(octane, n_compounds=10, box=[3, 3, 3], seed=2)\n",
    "\n",
    "octane_box.save('octane_system.gro', residues='Octane', overwrite=True)\n",
    "octane_box.save('octane_system.top', forcefield_files='mff_test2/mff_test2.xml', residues='Octane',\n",
    "                 combining_rule='geometric', overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, just to prove that these files were written correctly, we can take a quick peek."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cat octane_system.top"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also trivially change the force field we apply.  For example, instead of our made up force field, let us apply the OPLS_AA forcefield.  We have included a file that includes parameters for alkanes and alcohols which would apply to this system.  We viewed this file previously.  The only change in the save command below (aside from changing the name of the top file) is the forcefield_files argument we pass to the function.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cat utils/oplsaa-alcohol.xml"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "octane_box.save('octane_system_opls.top', forcefield_files='utils/oplsaa-alcohol.xml', residues='Octane',\n",
    "                 combining_rule='geometric', overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cat octane_system_opls.top"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.6"
  },
  "nbdime-conflicts": {
   "local_diff": [
    {
     "key": "language_info",
     "op": "add",
     "value": {
      "codemirror_mode": {
       "name": "ipython",
       "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.5.5"
     }
    }
   ],
   "remote_diff": [
    {
     "key": "language_info",
     "op": "add",
     "value": {
      "codemirror_mode": {
       "name": "ipython",
       "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.6.7"
     }
    }
   ]
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}