docs/tutorials/tutorial_polymers.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Building a Simple Alkane"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The purpose of this tutorial is to demonstrate the construction of an alkane polymer and provide familiarity with many of the underlying functions in mBuild. Note that a robust polymer construction recipe already exists in mBuild, which will also be demonstrated at the end of the tutorial."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting up the monomer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first step is to construct the basic repeat unit for the alkane, i.e., a $CH_2$ group, similar to the construction of the $CH_3$ monomer in the prior methane tutorial.  Rather than importing the coordinates from a pdb file, as in the previous example, we will instead explicitly define them in the class. Recall that distance units are nm in mBuild."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import mbuild as mb\n",
    "\n",
    "class CH2(mb.Compound):\n",
    "    def __init__(self):\n",
    "        super(CH2, self).__init__()\n",
    "        # Add carbon        \n",
    "        self.add(mb.Particle(name='C', pos=[0,0,0]), label='C[$]')\n",
    "        \n",
    "        # Add hydrogens \n",
    "        self.add(mb.Particle(name='H', pos=[-0.109, 0, 0.0]), label='HC[$]')    \n",
    "        self.add(mb.Particle(name='H', pos=[0.109, 0, 0.0]), label='HC[$]')\n",
    "        \n",
    "        # Add bonds between the atoms\n",
    "        self.add_bond((self['C'][0], self['HC'][0]))\n",
    "        self.add_bond((self['C'][0], self['HC'][1]))\n",
    "        \n",
    "        # Add ports anchored to the carbon\n",
    "        self.add(mb.Port(anchor=self[0]), label='up')\n",
    "        self.add(mb.Port(anchor=self[0]), label='down')\n",
    "        \n",
    "        # Move the ports approximately half a C-C bond length away from the carbon\n",
    "        self['up'].translate([0, -0.154/2, 0]) \n",
    "        self['down'].translate([0, 0.154/2, 0]) \n",
    "\n",
    "monomer = CH2()\n",
    "monomer.visualize(show_ports=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This configuration of the monomer is not a particularly realistic conformation. One could use this monomer to construct a polymer and then apply an energy minimization scheme, or, as we will demonstrate here, we can use mBuild's rotation commands to provide a more realistic starting point. \n",
    "\n",
    "Below, we use the same basic script, but now apply a rotation to the hydrogen atoms.  Since the hydrogens start 180&deg; apart and we know they should be ~109.5&deg; apart, each should be rotated half of the difference closer to each other around the y-axis. Note that the rotation angle is given in radians. Similarly, the ports should be rotated around the x-axis by the same amount so that atoms can be added in a realistic orientation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import mbuild as mb\n",
    "\n",
    "class CH2(mb.Compound):\n",
    "    def __init__(self):\n",
    "        super(CH2, self).__init__()\n",
    "        # Add carbon\n",
    "        self.add(mb.Particle(name='C', pos=[0,0,0]), label='C[$]')\n",
    "        \n",
    "        # Add hydrogens \n",
    "        self.add(mb.Particle(name='H', pos=[-0.109, 0, 0.0]), label='HC[$]')    \n",
    "        self.add(mb.Particle(name='H', pos=[0.109, 0, 0.0]), label='HC[$]')\n",
    "        \n",
    "        # Rotate the hydrogens\n",
    "        theta = 0.5 * (180 - 109.5) * np.pi / 180\n",
    "        #mb.rotate(self['HC'][0], theta, around=[0, 1, 0])\n",
    "        #mb.rotate(self['HC'][1], -theta, around=[0, 1, 0])\n",
    "        self['HC'][0].rotate( theta, around=[0, 1, 0])\n",
    "        self['HC'][1].rotate(-theta, around=[0, 1, 0])\n",
    "        \n",
    "        # Add bonds between the atoms\n",
    "        self.add_bond((self['C'][0], self['HC'][0]))\n",
    "        self.add_bond((self['C'][0], self['HC'][1]))\n",
    "        \n",
    "        # Add the ports and appropriately rotate them\n",
    "        self.add(mb.Port(anchor=self[0]), label='up')\n",
    "        self['up'].translate([0, -0.154/2, 0]) \n",
    "        self['up'].rotate(theta, around=[1, 0, 0])\n",
    "        \n",
    "        self.add(mb.Port(anchor=self[0]), label='down')\n",
    "        self['down'].translate([0, 0.154/2, 0]) \n",
    "        self['down'].rotate(-theta, around=[1, 0, 0]) \n",
    "\n",
    "monomer = CH2()\n",
    "monomer.visualize(show_ports=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Defining the polymerization class"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With a basic monomer construct, we can now construct a polymer by connecting the ports together. Here, we first instantiate one instance of the CH2 class as ```1ast_monomer```, then use the clone function to make a copy. The ```force_overlap()``` function is used to connect the ```'up'``` port from ```current_monomer``` to the ```'down'``` port of ```last_mononer```."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class AlkanePolymer(mb.Compound):\n",
    "    def __init__(self):\n",
    "        super(AlkanePolymer, self).__init__()\n",
    "        last_monomer = CH2()\n",
    "        self.add(last_monomer)\n",
    "        for i in range(3):\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",
    "\n",
    "polymer = AlkanePolymer()\n",
    "polymer.visualize(show_ports=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualization of this structure demonstrates a problem; the polymer curls up on itself.  This is a result of the fact that ports not only define the location in space, but also an orientation. This can be trivially fixed, by  rotating the down port 180&deg; around the y-axis.  \n",
    "\n",
    "We can also add a variable ```chain_length``` both to the for loop and ```init```  that will allow the length of the polymer to be adjusted when the class is instantiated.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import mbuild as mb\n",
    "\n",
    "class CH2(mb.Compound):\n",
    "    def __init__(self):\n",
    "        super(CH2, self).__init__()\n",
    "         # Add carbons and hydrogens\n",
    "        self.add(mb.Particle(name='C', pos=[0,0,0]), label='C[$]')\n",
    "        self.add(mb.Particle(name='H', pos=[-0.109, 0, 0.0]), label='HC[$]')    \n",
    "        self.add(mb.Particle(name='H', pos=[0.109, 0, 0.0]), label='HC[$]')\n",
    "  \n",
    "        # rotate hydrogens\n",
    "        theta = 0.5 * (180 - 109.5) * np.pi / 180\n",
    "        self['HC'][0].rotate(theta, around=[0, 1, 0])\n",
    "        self['HC'][1].rotate(-theta, around=[0, 1, 0])\n",
    "        \n",
    "        # Add bonds between the atoms\n",
    "        self.add_bond((self['C'][0], self['HC'][0]))\n",
    "        self.add_bond((self['C'][0], self['HC'][1]))\n",
    "        \n",
    "        # Add ports\n",
    "        self.add(mb.Port(anchor=self[0]), label='up')\n",
    "        self['up'].translate([0, -0.154/2, 0]) \n",
    "        self['up'].rotate(theta, around=[1, 0, 0])\n",
    "        \n",
    "        self.add(mb.Port(anchor=self[0]), label='down')\n",
    "        self['down'].translate([0, 0.154/2, 0])\n",
    "        self['down'].rotate(np.pi, [0, 1, 0])\n",
    "        self['down'].rotate(-theta, around=[1, 0, 0]) \n",
    "\n",
    "\n",
    "class AlkanePolymer(mb.Compound):\n",
    "    def __init__(self, chain_length=1):\n",
    "        super(AlkanePolymer, self).__init__()\n",
    "        last_monomer = CH2()\n",
    "        self.add(last_monomer)\n",
    "        for i in range (chain_length-1):\n",
    "            current_monomer = CH2()\n",
    "    \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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "polymer = AlkanePolymer(chain_length=10)\n",
    "polymer.visualize(show_ports=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using mBuild's Polymer Class"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```mBuild``` provides a prebuilt class to perform this basic functionality. Since it is designed to be more general, it takes as an argument not just the chain length, but also the monomer and the port labels (e.g., 'up' and 'down', since these labels are user defined). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "polymer = mb.lib.recipes.Polymer(CH2(), 10, port_labels=('up', 'down'))\n",
    "polymer.visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Building a System of Alkanes "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A system of alkanes can be constructed by simply cloning the polymer constructed above and translating and/or rotating the alkanes in space. ```mBuild``` provides many routines that can be used to create different patterns, to which the polymers can be shifted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# create the polymer\n",
    "polymer = mb.lib.recipes.Polymer(CH2(), 10, port_labels=('up', 'down'))\n",
    "\n",
    "# the pattern we generate puts points in the xy-plane, so we'll rotate the polymer\n",
    "# so that it is oriented normal to the xy-plane\n",
    "polymer.rotate(np.pi/2, [1, 0, 0])\n",
    "\n",
    "# define a compound to hold all the polymers\n",
    "system = mb.Compound()\n",
    "\n",
    "# create a pattern of points to fill a disk\n",
    "# patterns are generated between 0 and 1, \n",
    "# and thus need to be scaled to provide appropriate spacing\n",
    "pattern_disk = mb.DiskPattern(50)\n",
    "pattern_disk.scale(5) \n",
    "\n",
    "# now clone the polymer and move it to the points in the pattern\n",
    "for pos in pattern_disk:\n",
    "    current_polymer = mb.clone(polymer)\n",
    "    current_polymer.translate(pos)\n",
    "    system.add(current_polymer)\n",
    "    \n",
    "system.visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Other patterns can be used, e.g., the ```Grid3DPattern```.  We can also use the rotation commands to randomize the orientation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import random\n",
    "\n",
    "polymer = mb.lib.recipes.Polymer(CH2(), 10, port_labels=('up', 'down'))\n",
    "system = mb.Compound()\n",
    "polymer.rotate(np.pi/2, [1, 0, 0])\n",
    "\n",
    "pattern_disk = mb.Grid3DPattern(5, 5, 5)\n",
    "pattern_disk.scale(8.0)\n",
    "    \n",
    "for pos in pattern_disk:\n",
    "    current_polymer = mb.clone(polymer)\n",
    "    for around in [(1, 0, 0), (0, 1, 0), (0, 0, 1)]:  # rotate around x, y, and z\n",
    "        current_polymer.rotate(random.uniform(0, np.pi), around)\n",
    "    current_polymer.translate(pos)\n",
    "    system.add(current_polymer)\n",
    "\n",
    "system.visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```mBuild``` also provides an interface to ```PACKMOL```, allowing the creation of a randomized configuration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "polymer = mb.lib.recipes.Polymer(CH2(), 5, port_labels=('up', 'down'))\n",
    "system = mb.fill_box(polymer, n_compounds=100, overlap=1.5, box=[10,10,10]) \n",
    "system.visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Variations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rather than a linear chain, the ```Polymer``` class we wrote can be easily changed such that small perturbations are given to each port. To avoid accumulation of deviations from the equilibrium angle, we will clone an unperturbed monomer each time (i.e., ```monomer_proto```) before applying a random variation.  \n",
    "\n",
    "We also define a variable ```delta```, which will control the maximum amount of perturbation. Note that large values of ```delta``` may result in the chain overlapping itself, as ```mBuild``` does not currently include routines to exclude such overlaps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import mbuild as mb\n",
    "\n",
    "import random\n",
    "\n",
    "class AlkanePolymer(mb.Compound):\n",
    "    def __init__(self, chain_length=1, delta=0):\n",
    "        super(AlkanePolymer, self).__init__()\n",
    "        monomer_proto = CH2()\n",
    "        last_monomer = CH2()\n",
    "        last_monomer['down'].rotate(random.uniform(-delta,delta), [1, 0, 0])\n",
    "        last_monomer['down'].rotate(random.uniform(-delta,delta), [0, 1, 0])\n",
    "        self.add(last_monomer)\n",
    "        for i in range(chain_length-1):\n",
    "            current_monomer = mb.clone(monomer_proto)\n",
    "            current_monomer['down'].rotate(random.uniform(-delta,delta), [1, 0, 0])\n",
    "            current_monomer['down'].rotate(random.uniform(-delta,delta), [0, 1, 0])\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",
    "\n",
    "polymer = AlkanePolymer(chain_length = 200, delta=0.4)\n",
    "polymer.visualize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  },
  "widgets": {
   "state": {
    "13dc75f8d2fc4f80a7f629a5a9f5dbd9": {
     "views": [
      {
       "cell_index": 25
      }
     ]
    },
    "14ffa5d91a0e4c38a3a404ae065a1313": {
     "views": [
      {
       "cell_index": 19
      }
     ]
    },
    "19b723d1ef3744af910bbc736f243f64": {
     "views": [
      {
       "cell_index": 18
      }
     ]
    },
    "1cc365c4b67740e990902533135c9b89": {
     "views": [
      {
       "cell_index": 25
      }
     ]
    },
    "30dadb6ffdc44e2ebf9003c2da261a7c": {
     "views": [
      {
       "cell_index": 9
      }
     ]
    },
    "348e47c083b94464b81f1f304a7df962": {
     "views": [
      {
       "cell_index": 15
      }
     ]
    },
    "46c25311de2f4c7690d28e8c373796c3": {
     "views": [
      {
       "cell_index": 9
      }
     ]
    },
    "4a65c610c7f9491982b4a40ad35e8ae9": {
     "views": [
      {
       "cell_index": 6
      }
     ]
    },
    "51619bcac3464aa8abd249aa50336912": {
     "views": [
      {
       "cell_index": 25
      }
     ]
    },
    "5d457c30ccd24c4ba920e443cc018b35": {
     "views": [
      {
       "cell_index": 12
      }
     ]
    },
    "6816d19b76bc4f73bd795c18eba0cf8c": {
     "views": [
      {
       "cell_index": 25
      }
     ]
    },
    "7ea9eec81642475f87d132172174eabb": {
     "views": [
      {
       "cell_index": 9
      }
     ]
    },
    "871689eb9fca49bdadab6309fd6bae37": {
     "views": [
      {
       "cell_index": 9
      }
     ]
    },
    "8d6c39fc9b354759a08d5e00f8326398": {
     "views": [
      {
       "cell_index": 9
      }
     ]
    },
    "9337b5e5d1b341eabfb486e91550c3f0": {
     "views": [
      {
       "cell_index": 18
      }
     ]
    },
    "94755327703e425db12abd2d59537621": {
     "views": [
      {
       "cell_index": 18
      }
     ]
    },
    "a660a5aa7cee45fa9caa13474e980ca4": {
     "views": [
      {
       "cell_index": 20
      }
     ]
    },
    "b4f9e409c7a94367a156187e2a3afaa2": {
     "views": [
      {
       "cell_index": 20
      }
     ]
    },
    "d0ffdeb89a254e05a06c5ae482957e8d": {
     "views": [
      {
       "cell_index": 18
      }
     ]
    },
    "da5fba9f177b4c109f92dad5f60ea85a": {
     "views": [
      {
       "cell_index": 22
      }
     ]
    },
    "da69e7285ef6409ca58d87622db87bfa": {
     "views": [
      {
       "cell_index": 9
      }
     ]
    },
    "eb96c7804aec442a9f201303b7fd905e": {
     "views": [
      {
       "cell_index": 4
      }
     ]
    },
    "f3f6f8cd7879436e9d46404efdc10ca5": {
     "views": [
      {
       "cell_index": 4
      }
     ]
    },
    "fea6d8d2746b46b2bd1edda7d91bbf6e": {
     "views": [
      {
       "cell_index": 18
      }
     ]
    }
   },
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}