overview.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MoSDeF - A Molecular Simulation Design Framework\n",
    "\n",
    "## Overview\n",
    "\n",
    "[MoSDeF](mosdef.org) consists of three core Python packages, [mBuild](https://github.com/mosdef-hub/mbuild), [Foyer](https://github.com/mosdef-hub/foyer) and [GMSO](https://github.com/mosdef-hub/gmso), that, when combined with tools for workflow and analysis management (such as [Signac and Signac-Flow](https://signac.io)) provide the means to perform complex molecular simulations in a **T**ransparent, **R**eproducible, **U**sable-by-other, and **E**xtensible (**TRUE**) manner. Reproducibility in this case is achieved by making all aspects of the simulation (system initialization, simulation execution, and analysis) scriptable, such that other researchers could execute your same scripts to achieve the same results. MoSDeF is also designed such that systems can be generated in a programmatic manner, facilitating screening of large structural/chemical parameter spaces.\n",
    "\n",
    "In this overview, we will be focusing specifically on the tools we've developed to address the issue of **_system initialization_**, including the creation of a molecular model, the application of a force field (atom-typing and parameter assignment), and output molecular files. The three tools contained within MoSDeF to address system initialization are:\n",
    "\n",
    "  - [**mBuild**](https://github.com/mosdef-hub/mbuild): A hierarchical, component-based molecule builder\n",
    "  \n",
    "  - [**Foyer**](https://github.com/mosdef-hub/foyer): A package for atom-typing as well as applying and disseminating forcefields\n",
    "\n",
    "  - [**GMSO**](https://github.com/mosdef-hub/gmso): A flexible data structure to store chemical topology for molecular simulation\n",
    "\n",
    "This overview is designed to introduce you to these tools in a general manner; however, more in-depth tutorials are also available from within the [mosdef_tutorials repository](https://github.com/mosdef-hub/mosdef_tutorials).\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "---\n",
    "\n",
    "### mBuild Units\n",
    "\n",
    "mBuild implicitly uses the below unit system:\n",
    "\n",
    "**Length**\n",
    "* nanometers [nm]\n",
    "\n",
    "**Angles**\n",
    "* Radians for all `Compound` operations\n",
    "* Degrees when building `Lattices`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Importing mBuild\n",
    "\n",
    "To begin, we need to import the `mbuild` package, here using the alias `mb`. This will give us access to all of the data structures and functions within mBuild. You can check to make sure the installed version is 0.10.1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import mbuild as mb\n",
    "print(mb.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The `Compound` class\n",
    "\n",
    "The building block of mBuild is the `Compound` class. A `Compound` can represent anything in the hierarchy of a chemical topology: an atom, a non-atomistic bead, a molecule, a chemical moeity, a section of a molecule, a collection of molecules, and even the entire system. Each level of the hierarchy inherits from the `Compound` class. This means that `Compounds` may contain other `Compounds`, and that the same methods and attributes are present for molecule components at any level of the hierarchy.\n",
    "\n",
    "<img src=\"./utils/hierarchical_design_image.png\" alt=\"Drawing\" style=\"width: 700px;\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's make a compound by creating an instance of the `Compound` class. The default arguments are fine for now, so we don't need to specify anything."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mycompound = mb.Compound()\n",
    "mycompound"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `Compound` class has a `visualize` method that enables in-notebook rendering."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mycompound.visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we've built a single-particle system, let's look at how we can build molecular structures. `Compounds` can be created by assembling systems piece-by-piece from constituent particles; however, it is typically more practical to load coordinates and bonds from a structure file (such as PDB, GRO, MOL2, XYZ, etc.) and operate on that structure.\n",
    "\n",
    "Let's look at how we can build a linear alkane chain using each approach.\n",
    "\n",
    "---\n",
    "\n",
    "## Creating an Alkane from CH2 building blocks\n",
    "\n",
    "In our first approach, we will generate a methylene (CH2) building block and use it to build an alkane (including capping the ends with hydrogens).\n",
    "\n",
    "### Loading from a PDB structure file\n",
    "\n",
    "First, let's load a CH2 moiety into an mBuild `Compound` by reading from a PDB structure file (created using [Avogadro](https://avogadro.cc/)). This creates an mBuild `Compound` containing three atoms (C, H, H), as well as two C-H bonds. Then, once it is loaded, we can visualize it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ch2 = mb.load('utils/ch2.pdb')\n",
    "ch2.visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note, formats such as PDB include bonding information.  One could presumably load other formats without bonding information, and specify these bonds manually.  Additionally, one can explicitly define atom locations and bonds; for example, see  [mBuild Tutorial 01: Basic Functionality](https://github.com/mosdef-hub/mbuild_tutorials/blob/master/mBuild_01_Basic_Functionality.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Ports` and `Compound` classes\n",
    "\n",
    "Connectivity in `mBuild` is handled with the `Port` class. Instances of this class define where bonds can potentially be formed. These can be thought of as a rough equivalent to dangling bonds.\n",
    "\n",
    "The CH2 moiety we just loaded does not have `Ports`, so we will need to add them ourselves.\n",
    "\n",
    "To make a `Port`, we will need to specify three things:\n",
    "- `anchor`: the `Compound` the `Port` is directly associated with\n",
    "- `orientation`: a vector describing the direction the `Port` points\n",
    "- `separation`: the distance from the anchor the `Port` will exist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select the carbon `Compound`, which we want to use as an achor\n",
    "carbon = list(ch2.particles_by_name('C'))[0]\n",
    "\n",
    "up_port = mb.Port(anchor=carbon, orientation=[0, 1, 0], separation=0.075)\n",
    "down_port = mb.Port(anchor=carbon, orientation=[0, -1, 0], separation=0.075)\n",
    "\n",
    "ch2.add(up_port, 'up')\n",
    "ch2.add(down_port, 'down')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To visualize our CH2 moiety with the `Ports` we just added, we can use the `show_ports` argument"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ch2.visualize(show_ports=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, if we wanted to do this for a number `N` of CH2 `Compounds` in an alkane of lengths `N`, the process would be quite cumbersome and include a lot of repeated code. Instead, we can create a reusable class that defines our CH2 `Compound`. This allows us to encapsulate the routines for creating this (or any!) molecular moiety into an object that can be instantiated and in a manner that is _reproducible_ and can easily be _shared with others_.\n",
    "\n",
    "Below is class definition that includes the _same_ commands we used above to load coordinates and bonds from a PDB structure file and add `Ports` to the carbon atom.\n",
    "\n",
    "For additional information, see [mBuild Tutorial 02: Reusing Components](https://github.com/mosdef-hub/mbuild_tutorials/blob/master/mBuild_02_Reusing_Components.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class CH2(mb.Compound):\n",
    "    def __init__(self):\n",
    "        super(CH2, self).__init__()\n",
    "        mb.load('utils/ch2.pdb', compound=self)\n",
    "        carbon = list(self.particles_by_name('C'))[0]\n",
    "        self.add(mb.Port(anchor=carbon, orientation=[0, 1, 0], separation=0.075), 'up')\n",
    "        self.add(mb.Port(anchor=carbon, orientation=[0, -1, 0], separation=0.075), 'down')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we instantiate this class and visualize we should see the same result we obtained earlier. We can pass the `show_ports=True` argument to `visualize` to see the `Ports` we've added to the carbon atom."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ch2 = CH2()\n",
    "ch2.visualize(show_ports=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Examining the `Compound` data structure\n",
    "\n",
    "Now that we have created a `Compound` we can examine the contents.  For example, simply calling the `Compound` will provide us with a summary of the contents."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ch2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use the `save` method to write this `Compound` to disk. A number of file formats are supported, each with their own uses; here we will use an `.xyz` file to show we can write out the coordinates of a `Compound`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ch2.save('ch2.xyz', overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's peek at the file we just wrote to see what we created:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cat ch2.xyz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use the `load` method to load this file back in:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ch2_xyz = mb.load('./ch2.xyz')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can examine the 3-D coordinates with the `.xyz` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ch2.xyz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use the list function to iterate over the atoms and their positions in the compound:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "list(ch2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To view the bonds, we can call the `bonds` method as part of the `Compound` taking advantage of `list`, including some ways of making the output prettier:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "list(ch2.bonds())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for pair in ch2.bonds():\n",
    "    print(pair[0].name, '-', pair[1].name)\n",
    "\n",
    "['{}-{}'.format(pair[0].name, pair[1].name) for pair in ch2.bonds()]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also view a summary of the `Ports` associated with a `Compound`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ch2.all_ports()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### More complex classes: connecting CH2 moieties into an alkane\n",
    "\n",
    "Now we'll create a more complex class that defines the instructions for connecting CH2 moieties into a linear alkane chain.\n",
    "\n",
    "The code below instantiates CH2 moieties inside of a `for` loop, where the number of iterations is dependent on the desired length of the chain. The length of the chain can be toggled through the `chain_length` argument provided to the class constructor. We also import a hydrogen `Compound` from mBuild's `atoms` library to cap the ends of our chain.\n",
    "\n",
    "This is shown pictorially below.\n",
    "<img src=\"./utils/figure_connecting.png\" alt=\"Drawing\" style=\"width: 700px;\"/>\n",
    "\n",
    "**Note:** For this general overview, we do not intend for users (particularly those new to Python and object-oriented programming) to get too bogged down in the syntax. Instead, the emphasis should be that with mBuild we can encapsulate a series of routines (a \"recipe\") into a class, and that these routines can be defined in a manner that gives the class structural/chemical flexibility.\n",
    "\n",
    "For additional examples, see tutorials: \n",
    "- [mBuild Tutorial 03: Connecting Components with Ports](https://github.com/mosdef-hub/mbuild_tutorials/blob/master/mBuild_03_Connecting_Components_with_Ports.ipynb) \n",
    "- [mBuild Tutorial 04: Constructing Larger Compounds](https://github.com/mosdef-hub/mbuild_tutorials/blob/master/mBuild_04_Constructing_Larger_Compounds.ipynb)\n",
    "- [mBuild Tutorial 05: Creating Flexible Classes](https://github.com/mosdef-hub/mbuild_tutorials/blob/master/mBuild_05_Creating_Flexible_Classes.ipynb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from mbuild.lib.atoms import H\n",
    "\n",
    "class Alkane(mb.Compound):\n",
    "    def __init__(self, chain_length=1):\n",
    "        super(Alkane, self).__init__()\n",
    "        last_monomer = CH2()\n",
    "        hydrogen = H()\n",
    "        mb.force_overlap(move_this=hydrogen,\n",
    "                         from_positions=hydrogen['up'],\n",
    "                         to_positions=last_monomer['up'])\n",
    "        self.add(last_monomer)\n",
    "        self.add(hydrogen)\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",
    "        hydrogen = H()\n",
    "        mb.force_overlap(move_this=hydrogen,\n",
    "                         from_positions=hydrogen['up'],\n",
    "                         to_positions=last_monomer['down'])\n",
    "        self.add(hydrogen)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because we've defined our class to take `chain_length` as an argument, we can toggle the chemistry of our system (in this case the number of carbons in a linear alkane) by changing the value we provide for this argument upon instantiation.\n",
    "\n",
    "For example, let's create a butane molecule."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "butane = Alkane(chain_length=4)\n",
    "butane.visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The geometry of this molecule is not entirely realistic as all backbone atoms featuring 180° angles in the alkane molecules, all hydrogen atoms in plane. This can be addressed by placing `Compounds` and `Ports` in more realistic locations, either manually or by using energy minimized inputs.\n",
    "\n",
    "Alternatively, a `Compound` can be constructed and then energy minimized, either through a simulation engine or using the energy_minimization method in `mBuild`, which uses the [Open Babel toolkit](http://openbabel.org/dev-api/). See tutorial [mBuild Tutorial 07: Energy Minimization](https://github.com/mosdef-hub/mbuild_tutorials/blob/master/mBuild_07_Energy_Minimization.ipynb) for more information about the use of this method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "butane.energy_minimize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's change the value of `chain_length` to create a decane."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "decane = Alkane(chain_length=10)\n",
    "decane.visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating an Alkane by importing a pre-defined class\n",
    "\n",
    "\n",
    "In addition to allowing chemistry to be easily tuned through the presence of arguments, by encapsulating `Compound` recipes into a class, we then never have to define it again(!) and can easily share our recipe with other researchers. In other words, if Researcher A at University Y creates a class defining routines for generating some complex molecular system, Researcher B at University Z can simply import this class and use it. We are currently at work towards the implementation of a plugin architecture for mBuild that will make it easier for researchers to disseminate the `Compound` classes they've created."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, an `Alkane` class already exists within the `utils` directory of this repository. We can import that class and in two lines of code (excluding the visualization) we'll be able to construct a linear alkane of arbitrary length."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from utils import Alkane\n",
    "octane = Alkane(chain_length = 8)\n",
    "octane.visualize()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One immediate concern you may have is that many of these molecules we've generated feature non-realistic configurations (such as the alkane above with all angles at 180 degrees). However, this can typically be resolved through an energy minimization of the system you create in whichever simulation package you like to use. Additionally, mBuild also offers a built-in energy minimization routine (built around OpenBabel and OpenMM) that will attempt to optimize a `Compound`'s geometry."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Altering `Compounds`: Creating an alcohol\n",
    "\n",
    "\n",
    "mBuild contains routines for the addition and removal of particles. Here, we'll explore this functionality by changing our hexane molecule into _hexanol_.\n",
    "Note, we could do this by manually changing the class itself, or simply by removing the temrinal hydration and adding a hydroxyl in its place. \n",
    "\n",
    "First, we'll define a class for a hydroxyl group featuring a single `Port` on the oxygen to represent the dangling bond.\n",
    "\n",
    "We'll also go ahead and instantiate this class and visualize the resulting `Compound`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class OH(mb.Compound):\n",
    "    def __init__(self):\n",
    "        super(OH, self).__init__()\n",
    "        self.add(mb.Particle(name='O', element=\"O\", pos=[0.0, 0.0, 0.0]), label='O')\n",
    "        self.add(mb.Particle(name='H', element=\"H\", pos=[0.0, 0.1, 0.0]), label='H')\n",
    "        self.add_bond((self['O'], self['H']))\n",
    "        self.add(mb.Port(anchor=self['O'], orientation=[0, -1, 0], separation=0.075), label='down')\n",
    "\n",
    "hydroxyl = OH()\n",
    "hydroxyl.visualize(show_ports=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, we will use the `label` assigned to the up cap to differentiate it from other hydogen atoms in the system. This label can be determined by examining the class source code itself included in util, or by quering the `Compound` instance itself. Below we can see that we have 3 labels, chain, 'up_cap', and 'down_cap'. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "list(octane.labels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll remove a hydrogen from one end of our octane. This will create a `Port` in it's place, representing the dangling bond on the carbon.\n",
    "\n",
    "In the class definition, the hydrogen atom we will remove was provided with a label, `up_cap`. We can use this label to refer to this hydrogen, and remove it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "octane.remove(octane[-1])\n",
    "octane.visualize(show_ports=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can use the `force_overlap` method to attach the hydroxyl cap to the octane with the dangling bond."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "available_ports = octane.all_ports()\n",
    "mb.force_overlap(move_this=hydroxyl,\n",
    "                 from_positions=hydroxyl['down'],\n",
    "                 to_positions=available_ports[0])\n",
    "octanol = mb.Compound(name='Octanol')\n",
    "octanol.add(octane, label='alkane')\n",
    "octanol.add(hydroxyl, label='hydroxyl')\n",
    "octanol.visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting up bulk systems\n",
    "\n",
    "\n",
    "We usually want to run simulations of many molecules, not just one. To accomplish this, mBuild includes a module that interfaces with [PACKMOL](http://m3g.iqm.unicamp.br/packmol/home.shtml) to initialize bulk systems. This module includes functions that allow users to access PACKMOL functionality through Python functions, bypassing the need to write a PACKMOL input script. Here we'll use `fill_box` to place ten octanol molecules into a 3nm x 3nm x 3nm box. We can provide a seed for PACKMOL's random number generator to ensure the configuration is reproducible.\n",
    "\n",
    "For additional information on other packing functions, see [mBuild Tutorial 06: Setting Up Bulk Systems](https://github.com/mosdef-hub/mbuild_tutorials/blob/master/mBuild_06_Setting_Up_Bulk_Systems.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "octanol_box = mb.fill_box(octanol, n_compounds=10, box=[3, 3, 3], seed=2)\n",
    "octanol_box.visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Surface functionalization\n",
    "\n",
    "mBuild also provides routines for functionalization surfaces. There are a few surfaces available within mBuild's `surfaces` library; however, in the future we are hoping to feature a more comprehensive `surfaces` plugin.\n",
    "\n",
    "Here, we'll load a surface of $\\beta$-cristobalite silica.\n",
    "\n",
    "**Note:** Harmless warning messages are currently generated by one of the packages mBuild depends on. To reduce clutter, we are filtering those here, so you can safely ignore the warnings filter.\n",
    "\n",
    "For additional information, see [mBuild Tutorial 09: Surface Functionalization](https://github.com/mosdef-hub/mbuild_tutorials/blob/master/mBuild_09_Surface_Functionalization.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings(action=\"ignore\", category=FutureWarning)\n",
    "\n",
    "from mbuild.lib.surfaces import Betacristobalite\n",
    "surface = Betacristobalite()\n",
    "surface.visualize(show_ports=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use the `TiledCompound` class to expand our surface if we wanted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tiled_surface = mb.recipes.TiledCompound(surface, n_tiles=(2, 1, 1))\n",
    "tiled_surface.visualize(show_ports=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now remove the end-hydrogen from our octanol molecule to generate a dangling bond/`Port` that we can use to attach copies to the surface."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "octanol.remove(octanol['alkane'][-1])\n",
    "octanol.add(octanol.all_ports()[0], 'down', containment=False)\n",
    "octanol.visualize(show_ports=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use mBuild's `Pattern` class to create a pattern to define the arrangement of molecules on the surface. Here, we will create a `Random2DPattern` of 50 points in *xy* space. The `apply_to_compound` method can be used to attach copies of a molecules to a surface at locations designated by the `Pattern`. We can also provide a backfill `Compound` (in this case hydrogen) to fill vacant sites."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pattern = mb.Random2DPattern(10)\n",
    "hydrogen = H()\n",
    "chains, backfills = pattern.apply_to_compound(guest=octanol, host=tiled_surface, backfill=hydrogen)\n",
    "functionalized_surface = mb.Compound(subcompounds=[tiled_surface, chains, backfills])\n",
    "functionalized_surface.visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now lets use a pre-built recipe that is available with mBuild. The `Monolayer` recipe from mBuild's `recipes` module. You can change `n` in the `Random2DPattern` instance to toggle the number of chains in the system. `tile_x` and `tile_y` can be toggled to expand the surface as desired."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mbuild.lib.surfaces import Betacristobalite\n",
    "surface = Betacristobalite()\n",
    "surface.visualize(show_ports=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pattern = mb.pattern.Random2DPattern(n=50, seed=123)\n",
    "monolayer = mb.recipes.Monolayer(pattern=pattern, surface=surface, chains=mb.recipes.Alkane(n=12,cap_end=False), tile_x=2, tile_y=1,)\n",
    "monolayer.visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Converting to GMSO and Applying force fields with `Foyer`\n",
    "\n",
    "In order to run a simulation of any system we have built with mBuild, we need to apply a force field and write the necessary data files. Here we will consider two of the input files needed to run a GROMACS simulation, the `.gro` and `.top` files. The `.gro` file contains atomic coordinates and box information. The `.top` file contains the relevant force field parameters. Other molecular simulation engines use different files but their contents overlap with these GROMACS files.\n",
    "\n",
    "\n",
    "Here, we willl utilize GMSO (General Simulation Molecular Object) library to write out these molecular files. The GMSO provides necessary data structure to store details related to molecular topology necessary to run a molecular simulation (e.g., atom positions, simulation box, mixing rules, potential parameters, etc.). This library also has a collection of writers, allow for outputting to multiple file format, including `.gro`, `.top`, `.lammps`, and more. A `mBuild.Compound` objcet can be converted to a `gmso.Topology` via the `mbuild.Compound.to_gmso()` method. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "octanol_box_top = octanol_box.to_gmso()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "octanol_box_top.save('system.gro', overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When first converted from a `mBuild.Compound`, the `gmso.Topoloy` object (in our case, `octanol_box_top`) only contains atoms (and their positions), bonds, and simulation box simulation. The other information, such as connection objects (angles, dihedrals, impropers) and potentials (atom type, bond type, etc.,) needed to be added. This process can either be done manually (e.g., if one wish to only include certain angles/dihedrals/impropers) or automatically using `Topology.identify_connections()`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Before identifying connections:\")\n",
    "print(octanol_box_top)\n",
    "\n",
    "octanol_box_top.identify_connections()\n",
    "print(\"After identifying connections:\")\n",
    "print(octanol_box_top)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the `.top` file contains the force field information, we need to apply a force field to our system. Here, we will use the OPLS all-atom force field. To do this, we will utilize the `foyer` library to perform the atom typing.\n",
    "\n",
    "`Foyer` force fields are defined within an XML file that contains both the \"rules\" required for atom-typing as well as the force field parameters within a single file. `Foyer` atom-types based on the local environment (chemical substructure) of an atom . The [SMARTS](http://www.daylight.com/dayhtml/doc/theory/theory.smarts.html) language can be used to define these chemical substructures. For example, the SMARTS string `\"[C;X4](F)(F)(F)(C)\"`, which prescribes the OPLS-AA atom type `opls_961`, states that this atom type should be used when:\n",
    "- the atom of interest is element carbon, i.e., `[C]`\n",
    "- the carbon has 4 connected neighbors, i.e., `[C;X4]`\n",
    "- three neighbors are F, i.e., `(F)(F)(F)`\n",
    "- one neighbor is C, i.e., `(C)`\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), which is an established file format for storing force field data. `foyer` extends this format to include information necessary for atom-typing, which is the process of specifying atom types in a system according to their definitions.\n",
    "\n",
    "The `AtomTypes` section of a `foyer` XML includes in each `Type` four additional attributes not found in OpenMM:\n",
    "* `def` - SMARTS string describing the chemical substructure of this atom type (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 atom type (optional)\n",
    "* `doi` - DOI reference for parameters associated with this atom type (optional)\n",
    "* `overrides` - One or more atom types to 'override', providing precedence to this atom type\n",
    "\n",
    "Let's first take a quick look at this file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "!cat utils/oplsaa-alcohol.xml"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With this force field XML file, `foyer` will use the SMARTS strings to atom-type our system and return a `type_map`. This `type_map` is a dictionary mapping each atom to the corresponding `atom_type` defined in the XML, and can be used by different data structures to continue subsequent steps in the simulation workflow. Within the MoSDeF ecosystem, the `type_map` can be taken up by `GMSO` to perform the final parameterization step, i.e., assigning connection types/potentials based on the assigned atom types/classes. This process is wrapped in the `gmso.parameterization` in `GMSO` as demonstrated below. \n",
    "\n",
    "**Note:** The warning message about unparameterized impropers can be safely ignored, as OPLS-AA 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": [
    "from gmso import ForceField\n",
    "from gmso.parameterization import apply\n",
    "\n",
    "oplsaa_alc = ForceField('utils/oplsaa-alcohol.xml')\n",
    "apply(\n",
    "    octanol_box_top,\n",
    "    oplsaa_alc\n",
    ")\n",
    "octanol_box_top.save('system.top', overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, just to check that these files were written correctly, we can look at their contents"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cat system.gro"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cat system.top"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Working with coarse-grained and united atom forcefields\n",
    "\n",
    "Foyer supports SMARTS definitions involving non-atomistic components such as coarse-grained or united atom beads. Non-elemental species can easily be defined by prepending the name of custom \"element\" with an underscore `_`.\n",
    "\n",
    "For example, the following lines could be used to describe beads representing `_CH2` 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",
    "- the bead of interest is `[_CH2]`\n",
    "- it has 2 connected neighbors, i.e., `[_CH2;X2]`\n",
    "- those neighbors can be either `_CH2` or `_CH3`, i.e., `([_CH3,_CH2])[_CH3,_CH2]`\n",
    "\n",
    "This uses the \"or\" operator in SMARTS, which is a comma. In this case, these neighbors can be _either_ `_CH2` or `_CH3`.\n",
    "\n",
    "\n",
    "For more information on nano-atomistic forcefields, see [Foyer Tutorial 02: SMARTS for Non-Atomistic Systems](https://github.com/mosdef-hub/foyer_tutorials/blob/master/Foyer_02_SMARTS_for_Non-Atomistic_Systems.ipynb).\n",
    "\n",
    "For additional information on Foyer itself, see the [Foyer tutorial repository](https://github.com/mosdef-hub/foyer_tutorials) and [github page](https://github.com/mosdef-hub/foyer)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This concludes the general MoSDeF overview. For more in-depth tutorials into mBuild and Foyer and some of their use cases, refer to any of these links:\n",
    "- [mosdef.org](mosdef.org)\n",
    "- [mosdef_tutorials](https://github.com/mosdef-hub/mosdef_tutorials)"
   ]
  },
  {
   "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}