{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Filling With Packmol\n",
    "\n",
    "While the MDMC `universe.fill` method is serviceable for simpler systems, it isn't particularly sophisticated. For more complex filled universes, MDMC has an integration with packmol, a package for creating atomic configurations such as bilayers or mixtures. If you are not aware of how packmol works, then please read the user guide [here](https://m3g.github.io/packmol/userguide.shtml)."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note to run this tutorial, you need to have Packmol installed; this is already installed in the Docker container. The MDMC interface to packmol is contained in the objects `MDMC.MD.packmol.PackmolSetup` and `MDMC.MD.packmol.PackmolFiller`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from MDMC.MD.packmol import PackmolSetup, PackmolFiller"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "## Creating mixtures\n",
    "\n",
    "The first example shown is a simple system of water mixed in with ethanol. In a 1:1 molecular ratio."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "Firstly, we create the molecules needed for the system (See the notebook on [creating atomic configurations](./creating-atomic-configurations.ipynb) for information on this.)\n",
    "\n",
    "(NB: It is highly recommend to add any atom labels, dihedrals, bond angles, etc. to the `Molecule`s at this stage to ensure they are present in the final universe, rather than having to go in and add them all manually later.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from MDMC.MD import Atom, Molecule, Bond, BondAngle\n",
    "H1 = Atom('H')\n",
    "H2 = Atom('H', position=[0., 1.63298, 0.])\n",
    "O = Atom('O', position=[0., 0.81649, 0.57736])\n",
    "h2o_bonds = Bond((H1, O), (H2, O))\n",
    "HOH_angle = BondAngle(H1, O, H2)\n",
    "water_molecule = Molecule(atoms=[H1, H2, O], interactions=[HOH_angle, h2o_bonds], name=\"water\")\n",
    "\n",
    "h_1 = Atom(\"H\", position=[-1.9237, 0.3850, 0.0000])\n",
    "h_2 = Atom(\"H\", position=[2.0985, 0.2306, 0.0000])\n",
    "h_3 = Atom(\"H\", position=[1.1184, -1.0093, 0.8869])\n",
    "h_4 = Atom(\"H\", position=[1.1184, -1.0093, -0.8869])\n",
    "h_5 = Atom(\"H\", position=[-0.0227, 1.1812, 0.8852])\n",
    "h_6 = Atom(\"H\", position=[-0.0227, 1.1812, -0.8852])\n",
    "c_1 = Atom(\"C\", position=[1.1879, -0.3829, 0.0000])\n",
    "c_2 = Atom(\"C\", position=[0.0000, 0.5526, 0.0000])\n",
    "o_1 = Atom(\"O\", position=[-1.1867, -0.2472, 0.0000])\n",
    "ch_bond = Bond((h_2, c_1), (h_3, c_1), (h_4, c_1), (h_5, c_2), (h_6, c_2))\n",
    "co_bond = Bond((c_2, o_1))\n",
    "cc_bond = Bond((c_1, c_2))\n",
    "oh_bond = Bond((o_1, h_1))\n",
    "oh_bond_angle = BondAngle(h_1, o_1, c_2)\n",
    "ethanol_molecule = Molecule(atoms=[h_1,h_2,h_3,h_4,h_5,h_6,c_1,c_2,o_1], interactions=[ch_bond, co_bond, cc_bond, oh_bond, oh_bond_angle], name=\"ethanol\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "Next, put this into a `PackmolSetup` object.\n",
    "\n",
    "Currently, only single fixed molecules, or boxes, cubes and spheres of molecules are supported.\n",
    "\n",
    "Here we add a cube of molecules with the `setup.add_cube` method. We provide a desired density\n",
    "in $Ang^-3$; alternatively one could provide a parameter `n_structures` for the desired number\n",
    "of `Molecule`s one would like in the space. Later we will see `setup.add_box` for a non-cubic\n",
    "box; furthermore we have `setup.add_sphere` to add a sphere of given radius."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "setup = PackmolSetup()\n",
    "# Two identically overlapping cubes will create a mixture of water and ethanol in a 1:1 ratio of molecules\n",
    "setup.add_cube(structure=water_molecule, size=40., density=0.01)\n",
    "setup.add_cube(structure=ethanol_molecule, size=40., density=0.01)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "Next step is to pass this `PackmolSetup` object to a `PackmolFiller` object.\n",
    "\n",
    "Calling `fill_with_packmol` will setup & run packmol in the background and return a filled universe with the molecules specified."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "filler = PackmolFiller(setup_data=setup)\n",
    "universe = filler.fill_with_packmol()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can use `view` to look at the solved system. This may be slow - switch the viewer to ASE if it is not loading."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from MDMC.gui import view\n",
    "\n",
    "view(universe)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating bilayers"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`PackmolSetup` boxes, cubes, and spheres can be given a centre, for objects not centred on the origin.\n",
    "\n",
    "Here we will demonstrate how this can be used to create an interface between benzene and water. We get our benzene molecule from file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from MDMC.readers.configurations import read\n",
    "\n",
    "setup_2 = PackmolSetup()\n",
    "benzene_atoms = read(\"data/benzene.pdb\")\n",
    "benzene_molecule = Molecule(atoms=benzene_atoms, name='benzene')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define the centre of our boxes through the `origin` parameter;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "setup_2.add_box(structure=water_molecule, lengths=(20., 10., 20.,), origin=(0., 10., 0.), density=0.01)\n",
    "setup_2.add_box(structure=benzene_molecule, lengths=(20., 10., 20.,), origin=(0., 20., 0.), density=0.01)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and as before, we fill with the `PackmolFiller`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filler_2 = PackmolFiller(setup_data=setup_2)\n",
    "universe_2 = filler_2.fill_with_packmol()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "view(universe_2)"
   ]
  }
 ],
 "metadata": {
  "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}