mbuild-overview-STUDENT.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MoSDeF - A Molecular Simulation Design Framework\n",
    "\n",
    "## Overview\n",
    "\n",
    "MoSDeF consists of two core Python package, [mBuild](https://github.com/mosdef-hub/mbuild) and [Foyer](https://github.com/mosdef-hub/foyer), 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 **reproducible** 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 and the application of a force field (atom-typing and parameter assignment). The two 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",
    "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",
    "---\n",
    "\n",
    "**Pre-requisites**\n",
    "\n",
    "We have designed this tutorial for users that have some knowledge of Python and object-oriented programming (OOP). However, we encourage all users to work through the notebook, even those new to the world of Python and OOP, in order to still obtain an idea of the general concept of our tools. The syntax can be picked up later.\n",
    "\n",
    "\n",
    "------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Primer on using Jupyter notebooks and Binder\n",
    "\n",
    "[Jupyter notebooks](https://jupyter-notebook.readthedocs.io/en/stable/) provide an interactive environment for \"developing, documenting, and executing code\". Several languages are supported; however here we will be using Python. \n",
    "\n",
    "Jupyter notebooks feature two primary types of cells:\n",
    "1. Markdown cells, like this cell, which contain explanatory text\n",
    "2. Code cells, that can be executed by either clicking on the \"run cell\" icon or by hitting SHIFT + ENTER.\n",
    "\n",
    "Cells do not have to be executed in order (however the cells in this tutorial are designed to be executed _sequentially_), and the order in which cells have been executed is recorded by the bracketed number to the left of the _code_ cell (e.g. [ 1 ]). When a cell is executed you will first see an asterisk (i.e. [ * ]) which means that the cell is still running. When the asterisk is replaced by a number this means the execution has completed.\n",
    "\n",
    "Markdown cells will _not_ have numbers to the left of their cell. These are text based and not meant to be considered executable code. Executing these cells will render the Markdown cells in HTML. More information can be found [here](https://www.markdownguide.org/getting-started)\n",
    "\n",
    "---\n",
    "\n",
    "[Binder](https://mybinder.readthedocs.io/en/latest/) provides the ability to deploy Jupyter notebooks in the cloud, such that users do not need to set up their own computing environment to execute the notebook cells.\n",
    "* We will not be using Binder during this session, but all of our notebooks are hosted on Binder as well.\n",
    "* Binder is a free service that is community supported, and can be slow to access with multiple users trying to access the same notebook at once.\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "### mBuild Units\n",
    "\n",
    "Within mBuild, units to describe various aspects of system initialization are kept constant within the package.\n",
    "This provides a controlled environment that limits possible Input/Output (IO) errors when reading in/saving your structure of interest to various simulation engines.\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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import mbuild as mb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The `Compound` class\n",
    "\n",
    "The base class of mBuild is the `Compound` class, which defines the primary building block used for constructing molecules. **Molecules are constructed hierarchically**; however, 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. mBuild `Compounds` feature [a variety of useful methods and attributes](http://mosdef-hub.github.io/mbuild/data_structures.html) to facilitate system construction.\n",
    "\n",
    "<img src=\"./utils/hierarchical_design_image.png\" alt=\"Drawing\" style=\"width: 700px;\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compounds can be created by generating and connecting particles one-by-one; however, it is typically more practical to load coordinates and bonds from a structure file (such as PDB, MOL2, etc.) or to import `Compound` class definitions that other users have already defined.\n",
    "\n",
    "We'll explore both approaches here for constructing a linear alkane chain.\n",
    "\n",
    "---\n",
    "\n",
    "## Creating an Alkane from CH2 building blocks\n",
    "\n",
    "In our first approach for creating an alkane we will set-up routines for connecting CH2 building blocks (and we will cap the ends of our chain with hydrogens).\n",
    "\n",
    "### Loading from a PDB structure file\n",
    "\n",
    "First, we'll load a CH2 moiety into an mBuild `Compound` by reading from a PDB structure file (created using [Avogadro](https://avogadro.cc/)). This will create an mBuild `Compound` containing three atoms (C, H, H), as well as two C-H bonds. The `visualize` method allows us to view our `Compound` directly within the notebook. This visualization is provided by [`nglview`](https://github.com/arose/nglview)."
   ]
  },
  {
   "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": [
    "### 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` \n",
    "within the notebook will provide us with a summary of the contents."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# simply call the compound to print a summary of the number particles and bonds\n",
    "ch2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can examine the coordinates in multiple ways as shown below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# view the coordinates of the atoms in the compound\n",
    "ch2.xyz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use the list function to iterate over the atoms and their positions in the compound\n",
    "list(ch2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To view the bonds, we can call the `bonds` function as part of the `Compound` taking advantage of `list`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#list the pairs of atoms that are bonded; each pair appears between parantheses, i.e., (atom1, atom2)\n",
    "list(ch2.bonds())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we can also format the output of bonds to simply list the pairs of bonded atoms by name alone\n",
    "for pair in ch2.bonds():\n",
    "    print('{}-{}'.format(pair[0].name, pair[1].name))\n",
    "\n",
    "# equivalent shorthand output using list comprehension\n",
    "['{}-{}'.format(pair[0].name, pair[1].name) for pair in ch2.bonds()]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Formatting the output and exploring `Compound`: more examples\n",
    "\n",
    "Lets take a look at all the underlying attributes and data structures that a `Compound` inherently contains.\n",
    "\n",
    "To do this, we can use the built-in Python function: `dir()`. When an object is passed to `dir()`, it will\n",
    "attempt to list all of the available attributes that object contains."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# call dir on ch2\n",
    "dir(ch2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Information overload! All attributes prefixed with an underscore `_attribute` are considered _private_. This means \n",
    "that you should not be interacting and changing values of an object by accessing their private attributes. \n",
    "Attributes without an underscore are public and can be accessed like so: `ch2.xyz`\n",
    "\n",
    "Double underscores as prefix and suffixes are usually special functions, you can access and edit these, but most \n",
    "cases probably will _not_ overwrite these. For more information check here: [https://docs.python.org/3.5/reference/datamodel.html#objects-values-and-types](https://docs.python.org/3.5/reference/datamodel.html#objects-values-and-types)\n",
    "\n",
    "Lets use some python methods and strip out all the underscore and doubleunderscore (also known as: _dunderscore_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# to look through all the available attributes of the Compound class\n",
    "# Call the __dir__ method on the class\n",
    "# accessing the objects \"__dir__\" method achieves the same results as using the dir() method\n",
    "for x in ch2.__dir__():\n",
    "    if '__' in str(x) or '_' in str(x)[0]:\n",
    "        continue # skip these attrs\n",
    "    else:\n",
    "        print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example\n",
    "Print out the **private** (single underscore at beginning) attributes of `ch2`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for x in ch2.__dir__():\n",
    "    if '_' in str(x)[0] and '__' not in str(x):\n",
    "        print(x)\n",
    "    else:\n",
    "        continue"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using the `help` python method\n",
    "\n",
    "Python has a built-in function that is able to provide documentation \n",
    "on the fly. As long as a function or object contains a documentation \n",
    "string (_docstring_). `help()` will return that documentation\n",
    "interactively. If you are familiar with the unix command line, this is similar to `man`.\n",
    "\n",
    "Lets check out the help information for the `Compound` class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "help(help)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "help(mb.Compound)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1\n",
    "* `Save` `ch2` out to an XYZ file with a unique name and read it back in as a new `Compound` with the variable name: new_ch2\n",
    "\n",
    "* List the coordinates and the bonds of this compound, compare to the original `ch2` compound"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# How many bonds do you expect to find\n",
    "ch2.save('ch2.xyz', overwrite=True)\n",
    "\n",
    "# list coordinates\n",
    "\n",
    "# list bonds\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_ch2.add_bond?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adding bonds to a `Compound`\n",
    "\n",
    "Can add bonds explicitly using `compound.add_bonds` method\n",
    "\n",
    "Cons:\n",
    "Not flexible\n",
    "Tedious"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class CH2(mb.Compound):\n",
    "    def __init__(self):\n",
    "        super(CH2, self).__init__()\n",
    "        \n",
    "        carbon = mb.Particle(pos=[0.0, 0.0, 0.0], name='C')\n",
    "        hydrogen0 = mb.Particle(pos=[0.1, 0.0, 0.0], name='H')\n",
    "        hydrogen1 = mb.Particle(pos=[-0.1, 0.0, 0.0], name='H')\n",
    "        self.add([carbon, hydrogen0, hydrogen1])\n",
    "        self.add_bond((carbon, hydrogen0))\n",
    "        self.add_bond((carbon, hydrogen1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_ch2 = CH2()\n",
    "print(test_ch2)\n",
    "test_ch2.visualize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(list(test_ch2.particles()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2\n",
    "Using the above `CH2` class as inspiration, and create a `CH4` `Compound` Class\n",
    "\n",
    "_Note: Do not worry about placing the hydrogens in chemically realistic locations_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class CH4(mb.Compound):\n",
    "    def __init__(self):\n",
    "        super(CH4, self).__init__()\n",
    "        # Fill in here\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ch4 = CH4()\n",
    "ch4.visualize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ch4.energy_minimize()\n",
    "ch4.visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading in Compounds\n",
    "\n",
    "We can create compounds in multiple ways, some of these methods have been shown above.\n",
    "\n",
    "One of the easiest ways is to create compounds by loading in from a file (`mb.load`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ch2 = mb.load('utils/ch2.pdb')\n",
    "ch2.visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Explicitly create `Particles` and `add` them to a `Compound`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ch2 = mb.Compound()\n",
    "carbon = mb.Particle(pos=[0.0, 0.0, 0.0], name='C')\n",
    "hydrogen0 = mb.Particle(pos=[0.1, 0.0, 0.0], name='H')\n",
    "hydrogen1 = mb.Particle(pos=[-0.1, 0.0, 0.0], name='H')\n",
    "ch2.add([carbon, hydrogen0, hydrogen1])\n",
    "ch2.add_bond((carbon, hydrogen0))\n",
    "ch2.add_bond((carbon, hydrogen1))\n",
    "ch2.visualize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class CH2(mb.Compound):\n",
    "    def __init__(self):\n",
    "        super(CH2, self).__init__()\n",
    "        carbon = mb.Particle(pos=[0.0, 0.0, 0.0], name='C')\n",
    "        hydrogen0 = mb.Particle(pos=[0.1, 0.0, 0.0], name='H')\n",
    "        hydrogen1 = mb.Particle(pos=[-0.1, 0.0, 0.0], name='H')\n",
    "        self.add([carbon, hydrogen0, hydrogen1])\n",
    "        self.add_bond((carbon, hydrogen0))\n",
    "        self.add_bond((carbon, hydrogen1))\n",
    "\n",
    "ch2 = CH2()\n",
    "ch2.visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use `load` within your classes! Always remember to add to `self`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class CH2(mb.Compound):\n",
    "    def __init__(self):\n",
    "        super(CH2, self).__init__()\n",
    "        ch2 = mb.load('utils/ch2.pdb')\n",
    "        self.add(ch2)\n",
    "\n",
    "ch2 = CH2()\n",
    "ch2.visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Same code as above, without creating a temporary compound to add. `load` can add directly to a compound as well!"
   ]
  },
  {
   "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",
    "\n",
    "ch2 = CH2()\n",
    "ch2.visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use `load` to return the compound, but add it using `self.add` since `self` is a `Compound`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class CH2(mb.Compound):\n",
    "    def __init__(self):\n",
    "        super(CH2, self).__init__()\n",
    "        self.add(mb.load('utils/ch2.pdb'))\n",
    "\n",
    "ch2 = CH2()\n",
    "ch2.visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adding `Ports` to the `Compound` classes\n",
    "\n",
    "In order to connect `Compounds` we need to define locations where bonds can be formed (we can also add bonds manually through the `add_bond` method, which was covered above). mBuild handles this by allowing users to define `Ports` on particles, which essentially act as dangling bonds.\n",
    "\n",
    "However, if one had to re-write the commands for loading a CH2 molecule and adding `Ports` each time they wanted to create a molecule that included a CH2 unit, the process would be quite cumbersome. Instead, we can create a reusable class that defines our CH2 `Compound`. This approach allows one to encapsulate the routines for creating a molecular moiety into an object that can be instantiated and in a manner that can easily be shared with others.\n",
    "\n",
    "Below is a class definition for a CH2 moiety that uses the _same_ command we used above to load coordinates and bonds from a PDB structure file and features a few lines that 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",
    "        # Note: carbon is just a reference to the carbon particle in the compound self, NOT a copy\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": [
    "We can also view a summary of the  ports associated with a `Compound`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ch2 = CH2()\n",
    "ch2.all_ports()"
   ]
  },
  {
   "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": [
    "### Connecting two `CH2` compounds with `Ports`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "first_ch2 = CH2()\n",
    "second_ch2 = CH2()\n",
    "second_ch2.translate([1, 1, 1])\n",
    "parent = mb.Compound()\n",
    "parent.add((first_ch2, second_ch2))\n",
    "parent.visualize(show_ports=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To connect 2 ports, you use the `force_overlap` method and tell which particles/compounds to connect.\n",
    "\n",
    "This _physically_ moves and orients compounds based on the location and orientation defined in the `Port` class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mb.force_overlap(move_this=first_ch2,\n",
    "                 from_positions=first_ch2['up'],\n",
    "                 to_positions=second_ch2['down'])\n",
    "\n",
    "parent.visualize(show_ports=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Connecting compounds can be quite confusing if they are not wrapped in a `Compound` class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Ethane(mb.Compound):\n",
    "    def __init__(self):\n",
    "        super(Ethane, self).__init__()\n",
    "        first_ch2 = CH2()\n",
    "        second_ch2 = CH2()\n",
    "        second_ch2.translate([1, 1, 1])\n",
    "\n",
    "        self.add((first_ch2, second_ch2))\n",
    "        \n",
    "        mb.force_overlap(move_this=first_ch2,\n",
    "                 from_positions=first_ch2['up'],\n",
    "                 to_positions=second_ch2['down'])\n",
    "\n",
    "Ethane().visualize(show_ports=True)"
   ]
  },
  {
   "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",
    "        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)"
   ]
  },
  {
   "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 Particles 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 function 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 function and control of this function. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "butane.energy_minimize()\n",
    "butane.visualize()"
   ]
  },
  {
   "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.energy_minimize()\n",
    "decane.visualize()"
   ]
  },
  {
   "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": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 3:\n",
    "Add a port to the oxygen atom, labeled `down`"
   ]
  },
  {
   "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', pos=[0.0, 0.0, 0.0]), label='O')\n",
    "        self.add(mb.Particle(name='H', pos=[0.0, 0.1, 0.0]), label='H')\n",
    "        self.add_bond((self['O'], self['H']))\n",
    "        # add the port to the oxygen atom along the [0,-1, 0] direction\n",
    "        \n",
    "        \n",
    "hydroxyl = OH()\n",
    "hydroxyl.visualize(show_ports=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise : 4\n",
    "Convert the previous `Alkane` class to create an alcohol where it can accept the `up_cap` as a variable.\n",
    "\n",
    "#### Part 2: \n",
    "If no variable is passed as `up_cap`, make the default behavior to add a Hydroxyl group.\n",
    "\n",
    "A skeleton has been laid out below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mbuild.lib.atoms import H\n",
    "\n",
    "class Alcohol(mb.Compound):\n",
    "    def __init__(self, chain_length=1, up_cap=None):\n",
    "        super(Alcohol, self).__init__()\n",
    "\n",
    "        # check if nothing was passed as up_cap\n",
    "        if up_cap is None:\n",
    "            # input default behavior here\n",
    "\n",
    "        # something was passed as up_cap,\n",
    "        else:\n",
    "            pass\n",
    "\n",
    "        last_monomer = CH2()\n",
    "        \n",
    "        # Add to our `self` compound\n",
    "        self.add(up_cap)\n",
    "        self.add(last_monomer)\n",
    "        \n",
    "        # top capping CH2\n",
    "        mb.force_overlap(move_this=up_cap,\n",
    "                         from_positions=up_cap['down'],\n",
    "                         to_positions=last_monomer['up'])\n",
    "\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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "butane = Alcohol(chain_length=4, up_cap=None)\n",
    "butane.visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using `mbuild.Polymer`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Build an alkane using `Polymer`\n",
    "\n",
    "Note the absence of terminal groups\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alkane_block = mb.recipes.Polymer(monomers=CH2(), n=3, port_labels=('up', 'down'))\n",
    "alkane_block.visualize(show_ports=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create a semifluorinated alkane with 2 distinct blocks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class CF2(mb.Compound):\n",
    "    def __init__(self):\n",
    "        super(CF2, self).__init__()\n",
    "        \n",
    "        mb.load('./utils/cf2.pdb', compound=self)\n",
    "        carbon = list(self.particles_by_name('C'))[0]\n",
    "        up_port = mb.Port(anchor=carbon, orientation=[0, 0, 1], separation=0.075)\n",
    "        down_port = mb.Port(anchor=carbon, orientation=[0, 0, -1], separation=0.075)\n",
    "        self.add(up_port, label='up')\n",
    "        self.add(down_port, label='down')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pfa_block = mb.recipes.Polymer(monomers=CF2(), n=3, port_labels=('up', 'down'))\n",
    "alkane_block = mb.recipes.Polymer(monomers=CH2(), n=3, port_labels=('up', 'down'))\n",
    "\n",
    "semifluorinated_hexane = mb.recipes.Polymer(monomers=(alkane_block, pfa_block), sequence='ABAB', n=1, port_labels=('up', 'down'))\n",
    "\n",
    "semifluorinated_hexane.visualize(show_ports=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 5:\n",
    "Create an alcohol using the `Polymer` Class.\n",
    "\n",
    "#### Part2:\n",
    "Create an `Alcohol` class where the `chain_length` can be an argument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mbuild.lib.atoms import H\n",
    "\n",
    "class Alcohol(mb.Compound):\n",
    "    def __init__(self, chain_length=1, up_cap=None):\n",
    "        super(Alcohol, self).__init__()\n",
    "\n",
    "        # check if nothing was passed as up_cap\n",
    "        if up_cap is None:\n",
    "            # input default behavior here\n",
    "            \n",
    "        # something was passed as up_cap,\n",
    "        else:\n",
    "            pass\n",
    "        \n",
    "        # Add to our `self` compound\n",
    "        self.add(up_cap)\n",
    "        \n",
    "        # create the backbone using polymer\n",
    "        internal_chain =  # input polymer code here\n",
    "\n",
    "        \n",
    "        # top capping CH2\n",
    "        mb.force_overlap(move_this=up_cap,\n",
    "                         from_positions=up_cap['down'],\n",
    "                         to_positions=internal_chain['up'])\n",
    "\n",
    "        self.add(internal_chain)\n",
    "\n",
    "        # bottom cap\n",
    "        hydrogen = H()\n",
    "        mb.force_overlap(move_this=hydrogen,\n",
    "                         from_positions=hydrogen['up'],\n",
    "                         to_positions=internal_chain['down'])\n",
    "        self.add(hydrogen)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "octanol = Alcohol(chain_length=8,)\n",
    "octanol.visualize(show_ports=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Copying `Compounds`\n",
    "\n",
    "Copying a compound requires using the `mb.clone` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ch2_1 = CH2()\n",
    "\n",
    "# This does NOT make a copy of ch2_1, but a reference to ch2_1\n",
    "ch2_2 = ch2_1\n",
    "\n",
    "\n",
    "print(id(ch2_1), id(ch2_2))\n",
    "\n",
    "ch2_2.rotate(around=[0,0,1], theta=1.5)\n",
    "\n",
    "# They are the exact same\n",
    "print(ch2_2.xyz)\n",
    "print(ch2_1.xyz)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ch2_1 = CH2()\n",
    "\n",
    "# This DOES make a unique copy of ch2_1\n",
    "ch2_2 = mb.clone(ch2_1)\n",
    "\n",
    "\n",
    "print(id(ch2_1), id(ch2_2))\n",
    "\n",
    "ch2_2.rotate(around=[0,0,1], theta=1.5)\n",
    "\n",
    "# They are NOT the exact same!\n",
    "print(ch2_2.xyz)\n",
    "print(ch2_1.xyz)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting up bulk systems\n",
    "\n",
    "\n",
    "Typically we aren't desiring to run simulations of a single molecule. Fortunately, mBuild offers several routines to help create more complex systems. \n",
    "\n",
    "For example, mBuild provides users with an interface to [PACKMOL](http://m3g.iqm.unicamp.br/packmol/home.shtml) to set up bulk systems through the `fill_box` function. 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, see [mBuild Tutorial 06: Setting Up Bulk Systems](https://github.com/mosdef-hub/mbuild_tutorials/blob/master/mBuild_06_Setting_Up_Bulk_Systems.ipynb).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "octanol = Alcohol(chain_length=8)\n",
    "\n",
    "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['H'])\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": [
    "Again, the above routines could be wrapped into a class. Here we'll load the `AlkaneMonolayer` class from mBuild's `examples` 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, while `chain_length` can be toggled to alter the lengths of the alkylsilane chains."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mbuild.examples import AlkaneMonolayer\n",
    "pattern = mb.Random2DPattern(n=50, seed=123)\n",
    "monolayer = AlkaneMonolayer(pattern=pattern, tile_x=2, tile_y=1, chain_length=12)\n",
    "monolayer.visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This concludes the mBuild overview. For more in-depth tutorials into mBuild and Foyer, refer to the [mosdef_tutorials repository](https://github.com/mosdef-hub/mosdef_tutorials) or use our [Binder link](https://mybinder.org/v2/gh/mosdef-hub/mosdef_tutorials/master)."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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": 1
}