Session1-MoSDeF_Tutorial/Building_an_Alkane.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# mBuild tutorial: Building an alkane"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To begin, we need to provide some import statements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import mbuild as mb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now 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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ch2 = mb.load('ch2.pdb')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "mBuild has a built-in `Compound.visualize` method for use in Jupyter notebooks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ch2.visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "mBuild Compound's have a variety of attributes and methods (http://mosdef-hub.github.io/mbuild/data_structures.html). These can be queried using `dir(ch2)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dir(ch2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One of the methods of mBuild `Compounds` is `particles()`, which is a generator that returns all of the `Particles` within a `Compound`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "list(ch2.particles())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another method, `particles_by_name` returns only `Particles` with a specified name. Here we will use the `particles_by_name` method to grab the carbon atom in our `Compound`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "carbon = list(ch2.particles_by_name('C'))[0]\n",
    "carbon"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `Port` class (http://mosdef-hub.github.io/mbuild/data_structures.html#mbuild.port.Port) is used to create ports, essentially dangling bonds, that can be used to connect `Compound`'s together and can accept several helpful arguments upon instantiation to orient the port in the desired direction.\n",
    "\n",
    "Here we will create two ports to add to our CH2 moiety to represent the two dangling bonds. We will orient these in opposite directions (+*y* and -*y*) and shift these from the carbon position by half of a C-C bond length\n",
    "- **Note**: mBuild uses _nanometers_ as the distance unit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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",
    "up_port"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can use the `add` function to add the ports we've just created to our `Compound`. `add` can be used to add `Compound`'s to other `Compound`'s."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ch2.add(up_port, label='up')\n",
    "ch2.add(down_port, label='down')\n",
    "ch2.available_ports()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Any part added to a `Compound` can be given an optional, descriptive string `label` for easy referencing. In the above code block, we added the labels `up` and `down` to the two `Port`'s that we added to our CH2 moiety."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ch2['up']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's visualize the ch2 `Compound` again, this time showing the two `Ports` we've added."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ch2.visualize(show_ports=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The mBuild `clone` function is an alternative to `deepcopy` and allows us to create copies of `Compound`'s. (If you are unfamiliar with the concept of a deep copy, here is a description copied from the Python manual: *A deep copy constructs a new compound object and then, recursively, inserts copies into it of the objects found in the original.*)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ch2_copy = mb.clone(ch2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `force_overlap` function takes a `Compound` and then rotates and translates it such that two other `Compounds` overlap. Typically, as in the following case, those two other `Compounds` are `Ports` - in our case, `ch2_copy['up']` and `ch2['down']`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mb.force_overlap(move_this=ch2_copy,\n",
    "                 from_positions=ch2_copy['up'],\n",
    "                 to_positions=ch2['down'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now have two CH2 moieties bonded together. Next, we'll add these two `Compounds` to a parent `Compound` that we will call `alkane`, which we can visualize or otherwise manipulate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alkane = mb.Compound()\n",
    "alkane.add(ch2, label='ch2[$]')\n",
    "alkane.add(ch2_copy, label='ch2[$]')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alkane.visualize(show_ports=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alkane"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It would be cumbersome to have to create `Compounds` by hand each time we wanted to use them. Instead, mBuild is intended to be used in such a way that reusable classes are defined for various `Compounds`. Here, we'll define a class for our CH2 moiety."
   ]
  },
  {
   "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('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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ch2 = CH2()\n",
    "ch2.visualize(show_ports=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have a class defined for our CH2 moiety, we can utilize this to create an alkane chain. Here we'll define a class for an alkane chain that can take `chain_length` as an argument. We will achieve this by successively adding CH2 moieties, and capping the first and last moieties with hydrogen atoms."
   ]
  },
  {
   "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": [
    "Now that we've created an `Alkane` class, we can instantiate this and provide any chain length (>0) that we want."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alkane = Alkane(chain_length=10)\n",
    "alkane.visualize(show_ports=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One downside to creating compounds in this manner is that you often end up with non-realistic configurations (such as the alkane above with all angles at 180 degrees). This is typically resolved through an energy minimization of the system you create in whichever simulation package you are using. However, mBuild also offers a built-in energy minimization routine that will attempt to optimize a `Compound`'s geometry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "alkane.energy_minimization(steps=5000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alkane.visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "mBuild also provides a `Polymer` recipe that essentially performs the same operations as the `Alkane` class we've defined above, but provides additional flexibility to support arbitrary monomer units as well as copolymers. Here we'll modify our `Alkane` class to use the `Polymer` recipe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Alkane(mb.Compound):\n",
    "    def __init__(self, chain_length=1):\n",
    "        super(Alkane, self).__init__()\n",
    "        chain = mb.Polymer(CH2(), chain_length)\n",
    "        self.add(chain, label='chain')\n",
    "        hydrogen = H()\n",
    "        mb.force_overlap(move_this=hydrogen,\n",
    "                         from_positions=hydrogen['up'],\n",
    "                         to_positions=chain['up'])\n",
    "        self.add(hydrogen, label='up-cap')\n",
    "        hydrogen = H()\n",
    "        mb.force_overlap(move_this=hydrogen,\n",
    "                         from_positions=hydrogen['up'],\n",
    "                         to_positions=chain['down'])\n",
    "        self.add(hydrogen, label='down-cap')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alkane = Alkane(6)\n",
    "alkane.visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Particles can be easily added and removed in mBuild. Here, we'll explore this functionality by changing our hexane molecule into _hexanol_.\n",
    "\n",
    "First, we'll create a class for a hydroxyl group."
   ]
  },
  {
   "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",
    "        self.add(mb.Port(anchor=self['O'], orientation=[0, -1, 0], separation=0.075), label='down')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's visualize to see if everything looks good."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hydroxyl = OH()\n",
    "hydroxyl.visualize(show_ports=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll remove a hydrogen from one end of our hexane and add a hydroxyl.\n",
    "\n",
    "Below, we remove the hydrogen compound that has the label `up-cap`.\n",
    "We see that there is a single port available."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alkane.remove(alkane['up-cap'])\n",
    "ports = alkane.all_ports()\n",
    "ports"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we force the single available port to overlap with the `down` port of the hydroxyl compound."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mb.force_overlap(move_this=hydroxyl,\n",
    "                 from_positions=hydroxyl['down'],\n",
    "                 to_positions=ports[0])\n",
    "hexanol = mb.Compound()\n",
    "hexanol.add(alkane, label='alkane')\n",
    "hexanol.add(hydroxyl, label='hydroxyl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now visualize to make sure our hexanol was successfully created."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hexanol.visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Running simulations of a single small molecule would likely not be very interesting. mBuild offers several routines to help create more complex systems. One of these routines is the `fill_box` function, which can be used to fill a box of a user-defined size with a specified number of molecules. We'll check this out now by placing 10 hexanols in a 3nm x 3nm x 3nm box."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "box = mb.fill_box(hexanol, n_compounds=10, box=[3, 3, 3], seed=2)\n",
    "box.visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now if we wanted to actually run a simulation of this system 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. Here we will save in Gromacs `TOP` and `GRO` formats."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "box.save('system.top', forcefield_name='oplsaa', overwrite=True)\n",
    "box.save('system.gro', overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The last bit of functionality we will observe is attaching molecules to a surface. In our case here, we first need to remove the hydrogen atom opposite of the hydroxyl so that we have an available port on our hexanol to attach to a surface."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hexanol.remove(hexanol['alkane']['down-cap'])\n",
    "hexanol.add(hexanol.all_ports()[0], 'down', containment=False)\n",
    "hexanol.visualize(show_ports=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The surface we will attach our hexanes to is crystalline silica, which we will import from mBuild's compound library. In many cases it is desirable to expand `Compounds` with periodicity, like surfaces, in one or more dimensions. mBuild features a built-in routine, `TiledCompound` to handle this functionality. Here we will use `TiledCompound` to expand our surface in the *x* dimension."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mbuild.lib.surfaces import Betacristobalite\n",
    "surface = Betacristobalite()\n",
    "tiled_surface = mb.TiledCompound(surface, n_tiles=(2, 1, 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "mBuild contains a `Pattern` class which can be used to create patterns to place molecules at desired locations in space. This is useful for many cases, including defining the arrangement of chains on a surface. Here, we will create a `Random2DPattern` of 10 points in *xy* space.\n",
    "\n",
    "The `Pattern.apply_to_compound` method can be used to attach `Compounds` to another `Compound` (e.g. a surface) by finding the vacant ports in the `host` `Compounds` closest to those defined by the pattern. Additionally, a `backfill` can be defined that will fill any leftover ports. Here, we will define a hydrogen atom as our backfill."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pattern = mb.Random2DPattern(10)\n",
    "hydrogen = H()\n",
    "chains, backfills = pattern.apply_to_compound(guest=hexanol, host=tiled_surface, backfill=hydrogen)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we'll add our surface, chains, and backfills to a parent `Compound` and visualize."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "functionalized_surface = mb.Compound()\n",
    "for part in [tiled_surface, chains, backfills]:\n",
    "    functionalized_surface.add(part)\n",
    "functionalized_surface.visualize()"
   ]
  }
 ],
 "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.5.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}