{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solvating an MDMC Universe  \n",
    "\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For simulating systems in solution, MDMC's inbuilt ``solvate`` method can be used to add solvent molecules of a desired density to your universe.\n",
    "\n",
    "\n",
    "\n",
    "## Create the Universe\n",
    "See the [_Building a Universe_](https://github.com/MDMCproject/MDMCv0.2_pilot/wiki/Building-a-Universe)  wiki for details on how to create a Universe of your own specifications. For the purposes of this tutorial, we will be solvating a simple universe that contains only 4 Hydrogen molecules:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import the Atom, Molecule, and Universe classes\n",
    "# Import the HarmonicPotential class (needed to create a Bond)\n",
    "from MDMC.MD import Atom, Bond, HarmonicPotential, Molecule, Universe\n",
    "\n",
    "# Initialise a Universe with dimensions in Ang\n",
    "universe = Universe([10.0, 15.0, 20.0])\n",
    "\n",
    "# Create a pair of Hydrogen Atoms\n",
    "H1 = Atom('H')\n",
    "H2 = H1.copy(position=[1., 1., 1.])\n",
    "\n",
    "# Initialise a H-H Bond\n",
    "HH_bond = Bond(H1, H2, function=HarmonicPotential(1., 100., interaction_type='bond'))\n",
    "\n",
    "# Make a H2 Molecule\n",
    "H2_mol_1 = Molecule(atoms=[H1, H2], interactions=[HH_bond])\n",
    "\n",
    "# Create 3 copies of the Molecule at different positions\n",
    "H2_mol_2 = H2_mol_1.copy(position=[3.5, 3.5, 3.5])\n",
    "H2_mol_3 = H2_mol_1.copy(position=[6.5, 6.5, 6.5])\n",
    "H2_mol_4 = H2_mol_1.copy(position=[9.5, 9.5, 9.5])\n",
    "\n",
    "# Add the 4 Hydrogen Molecules to the Universe\n",
    "for molecule in [H2_mol_1, H2_mol_2, H2_mol_3, H2_mol_4]:\n",
    "    universe.add_structure(molecule)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solvating the Universe\n",
    "The ```solvate``` method accepts 3 main parameters; ``density``, ``tolerance``, and ``solvent`` (as well as some ``**settings``). An example call to ```solvate``` the Universe created above would therefore be:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "universe.solvate(0.6, tolerance=1, solvent='SPCE')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parameters \n",
    "\n",
    "Explanations of each parameter passed to ``solvate`` can be found below.\n",
    "\n",
    "\n",
    "#### 1. ``density``\n",
    "The desired density of the bulk solvent in your universe, in MDMC units of **amu / Ang ^ 3**. See [Units in MDMC](units.ipynb) for instructions on how to convert your density into MDMC units.\n",
    "\n",
    "\n",
    "**Note**: with this parameter you are specifying the **bulk density of the solvent**. If you have any solute molecules already present in your universe, the total density of your universe after solvation will be higher than the desired density you pass to ``solvate`` (plus or minus the tolerance you pass, see below).\n",
    "\n",
    "In the above example, passing the density as ``0.6`` means that the universe will be solvated with SPCE water with a bulk density of 0.6 amu / Ang ^ 3 (+/- the tolerance). The density of the universe **in total will be greater** than 0.6 amu / Ang ^ 3, because of the 4 Hydrogen molecules.\n",
    "\n",
    "\n",
    "\n",
    "#### 2. ``tolerance``\n",
    "With this parameter you can specify the percentage tolerance of the ``density`` that you would like to be achieved for the bulk density of the solvent. \n",
    "\n",
    "For the above example, passing a density of ``0.6`` and setting ``tolerance=1`` will achieve a bulk solvent density of:\n",
    "\n",
    "* (0.6 amu / Ang ^ 3)  +/-  1 %\n",
    "* equivalent to (0.6  +/-  0.006) amu / Ang ^ 3\n",
    "\n",
    "**Note**: the tolerance has a default value of 1 %. Setting the tolerance to anything lower than increases the risk of ``solvate`` not converging to within the specified tolerance.\n",
    "\n",
    "\n",
    "\n",
    "#### 3. ``solvent``\n",
    "\n",
    "\n",
    "Using a solvent with pre-defined coordinates.\n",
    "\n",
    "MDMC has a few in-built solvents that can be used to solvate your Universe. These have pre-defined atomic coordinates and interactions. \n",
    "\n",
    "The in-built solvents are listed in the help method for solvate:  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "help(universe.solvate)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Example: solvating with SPCE water"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "universe.solvate(0.6, tolerance=1, solvent='SPCE')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 4. ``**settings``\n",
    "\n",
    "##### ``constraint_algorithm``\n",
    "\n",
    "You can specify a `ConstraintAlgorithm` which is applied to the Universe. If specifying ``solvent`` as a string representing one of the in-built solvents (i.e. 'SPCE'), then a ``Shake(1e-4, 100)`` constraint algorithm is applied by default.\n"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "add738a18a8c9080256f416f48e2f1975a0b0948846decec7f4bbc677823fdf5"
  },
  "kernelspec": {
   "display_name": "Python 3.7.4 32-bit ('mdmc-venv': virtualenv)",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": ""
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}