{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Creating an atomic configuration\n",
    "An MDMC simulation requires a configuration and a topology defined within a `Universe` object. These are orthorhombic (cubic or cuboid) boxes containing a configuration of atoms and molecules. This how-to guide will explain how to create a `Universe` object and an atomic configuration inside it.\n",
    "\n",
    "In addition to creating molecules manually, it is also possible to create them from atomic configuration files (e.g CIF files).  Please see the guide on [Reading atoms from configuration files](read-configurations.ipynb) for a detailed description on how to do this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from MDMC.MD import Universe"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Creating the actual box that contains the atoms is trivial; specify a list or tuple of 3 floats, representing the x, y, and z lengths of the box respectively; these lengths are taken to be in angstroms (Å)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "box_universe = Universe((10.0, 15.0, 20.0))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To create a cubic `Universe`, only a single float needs to be specified for the `dimensions`; MDMC will infer that this is a cube of side length 10Å. This cube will be used as our base universe for the rest of the guide."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "universe = Universe(10.0)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create an atomic configuration\n",
    "Each atom is specified using an `Atom` object. At a minimum these can be specified just from an element symbol."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from MDMC.MD import Atom\n",
    "\n",
    "# Create some atoms!\n",
    "H1 = Atom('H')\n",
    "O1 = Atom('O', position=(8.0, 8.0, 8.0), atom_type=1)\n",
    "O2 = Atom(element='O', position=(9.0, 9.0, 9.0), velocity=(0., 0., 0.), charge=None)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Atom`s have a few parameters here! Let's break down some of the particularly important ones:\n",
    "- `element`: the first parameter gives the element of the atom. If atom `mass` is not given, it will be determined from an elemental lookup table.\n",
    "- `position`: the (x,y,z) coordinates of the atom. If not given, defaults to the origin.\n",
    "- `velocity`: the velocity of the atom along the (x,y,z) axes. Defaults to (0,0,0).\n",
    "- `charge`: the atomic charge of the atom. Defaults to None (no charge).\n",
    "- `atom_type`: an ID for the type of atom; if not given, is inferred based on the element and interactions of the atom, such that all atoms with the same element and interactions will have the same atom types. This is used to keep track of atom interactions. \n",
    "\n",
    "It's recommended you either define no atom types (and let MDMC do it) or define all atom types. Mixing the two may create unexpected behaviour.\n",
    "\n",
    " Note that if set, the `velocity` of atoms will be scaled when creating a `Simulation` in order to ensure the temperature is accurate. Otherwise, if the velocities of all `Atom` objects in a `Simulation` are 0, then the velocities of the atoms in the MD engine will be randomly chosen in order to provide an accurate temperature. For more details see [Running a Simulation](running-a-simulation.ipynb).\n",
    "\n",
    "Atoms can also be created by copying another atom and providing the position of the new atom:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H2 = H1.copy(position=[1., 1., 1.])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The copied atom will have identical properties to the original attribute, except for a different `position`.  This includes interactions, which will apply to the copied atom in the same way as the original atom e.g. if H1 is bonded to an atom O, H2 will also be bonded to atom O.  An example of this is shown in the section **'Create bonded interactions'**."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create bonded interactions\n",
    "There are three bonded interaction types within MDMC: `Bond`, `BondAngle` and `Dihedral`.  Each `Interaction` must have an `InteractionFunction` which describes the `Interaction`."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Bond`\n",
    "This interaction represents a bond between two atoms. Below this is demonstrated where a `Bond` with a `HarmonicPotential` function is created."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import Bond and HarmonicPotential\n",
    "from MDMC.MD import Bond\n",
    "from MDMC.MD import HarmonicPotential\n",
    "# Create a Bond with a HarmonicPotential\n",
    "# The first argument in the HarmonicPotential is the equilibrium state and the second is the potential strength\n",
    "HH_bond = Bond(H1, H2, function=HarmonicPotential(1., 100., interaction_type='bond'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see which units are supported:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from MDMC.common.units import SYSTEM\n",
    "SYSTEM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `BondAngle`"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A `BondAngle` represents bonds at a rotation around a central atom, and is created in the same manner except it requires a minimum of three atoms. The central atom should be the 2nd atom. For example, a water molecule would have:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from MDMC.MD import BondAngle\n",
    "\n",
    "H1 = Atom('H')\n",
    "H2 = Atom('H', position=[0., 1.63298, 0.])\n",
    "O = Atom('O', position=[0., 0.81649, 0.57736])\n",
    "HOH_angle = BondAngle(H1, O, H2)\n",
    "\n",
    "# The following is equivalent\n",
    "HOH_angle = BondAngle(H2, O, H1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Currently `HarmonicPotential` is the only `InteractionFunction` that can be applied to either `Bond` or `BondAngle`."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All `Bond` and `BondAngle` interactions can have [constraints](https://en.wikipedia.org/wiki/Constraint_(computational_chemistry)) imposed on them; this can either be set when creating the `Bond` (or `BondAngle`) or afterwards:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "HH_bond.constrained = True\n",
    "HOH_angle = BondAngle(H1, O, H2, constrained=True)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a constraint to be applied during a simulation, the `Universe` must have a `ConstraintAlgorithm`; these are described at the end of this guide.\n",
    "\n",
    "*NB: \"Constraining\" the ability of the bonds to oscillate during MD in this way should not be confused with constraining the value of an MDMC `Parameter` to a certain numerical range during refinement, as described in [Running a Refinement](../../../tutorials/running-a-refinement.ipynb). It is entirely possible to have a rigid bond but allow the length of that bond to change between refinement steps, or conversely have a bond that is free to oscillate during MD but the equilibrium length is not altered as part of the refinement.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `DihedralAngle`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A `DihedralAngle` is also created in the same manner, except it requires four atoms.  A `DihedralAngle` can be proper or improper, as specified by `DihedralAngle.improper`.  By default a `DihedralAngle` is proper.\n",
    "\n",
    "The atoms in a proper `DihedralAngle` must be specified so that the 2nd and 3rd atoms are the central two atoms:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from MDMC.MD import DihedralAngle\n",
    "\n",
    "C1 = Atom('C', position=[5.12033922, 4.63847287, 4.94075943])\n",
    "N = Atom('N', position=[5.12991894, 3.78609704, 3.79996577])\n",
    "C2 = Atom('C', position=[4.91462725, 4.08992816, 2.5091264])\n",
    "O = Atom('O', position=[4.67405373, 5.24130678, 2.1180462])\n",
    "proper = DihedralAngle(atoms=[C1, N, C2, O])\n",
    "\n",
    "# For a proper DihedralAngle, the equivalent atom order is:\n",
    "proper = DihedralAngle(atoms=[O, C2, N, C1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The atoms in an improper `DihedralAngle` must be specific so that the 1st atom is the central atom."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = Atom('N', position=[4.97080909, 2.91075722, 1.57280005])\n",
    "C1 = Atom('C', position=[5.12033922, 4.63847287, 4.94075943])\n",
    "C2 = Atom('C', position=[5.12991894, 3.78609704, 3.79996577])\n",
    "C3 = Atom('C', position=[4.91462725, 4.08992816, 2.5091264 ])\n",
    "improper = DihedralAngle(atoms=[N, C1, C2, C3], improper=True)\n",
    "\n",
    "# The following are some of the equivalent permutations\n",
    "# The only unique atom location is the first location (central atom)\n",
    "improper = DihedralAngle(atoms=[N, C2, C1, C3], improper=True)\n",
    "improper = DihedralAngle(atoms=[N, C3, C1, C2], improper=True)\n",
    "improper = DihedralAngle(atoms=[N, C1, C3, C2], improper=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Currently `Periodic` is the only `InteractionFunction` that can be applied to `DihedralAngle` interactions."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Copying bonded atoms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As mentioned above in the **'Create an atom'** section, if an `Atom` with a `BondedInteraction` is copied, the new atom will also have the same `BondedInteraction` (and be bonded to the same atom or atoms as the original).  For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wH1 = Atom('H', position=(5., 5., 5.))\n",
    "wO = Atom('O', position=(5., 6.63298, 5.))\n",
    "wBond = Bond((wH1, wO), function=HarmonicPotential(1., 100., interaction_type='bond'))\n",
    "wH2 = wH1.copy(position=(0., 0.81649, 0.57736))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Both atoms `wH1` and `wH2` are bonded to `wO`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wBond.atoms"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a molecule\n",
    "Bonded atoms can be combined into a `Molecule` object, which keeps track of the atoms in the molecule and the bonds between them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from MDMC.MD import Molecule, BondAngle\n",
    "# Create a H2 molecule\n",
    "H1.position = [0., 0., 0.]\n",
    "H2.position = [1., 1., 1.]\n",
    "H_mol = Molecule(position=[2., 1.5, 1.], atoms=[H1, H2], interactions=[HH_bond])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When a `Molecule` is created, the position of the atoms relative to one another is fixed.  The atoms are then moved so that the position of the molecular centre of mass is what was passed when creating the molecule.  In the example above, the atoms were at `[0., 0., 0.]` and `[1., 1., 1.]` before `H_mol` was created; therefore they will always be separated by `1.0 Ang` in each dimension, no matter where the molecule is moved to.  The molecular centre of mass is set to `[2., 1.5, 1.]`, so the atom positions are changed to `[2., 1.5, 1.]` and `[3., 2.5, 2.]` respectively.\n",
    "It is also possible to `copy` molecules:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H_mol2 = H_mol.copy(position=[5., 5., 5.])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When a `Molecule` is copied, each `Atom` is copied, as are all of the bonded interactions between these atoms (and all of the non-bonded interactions).\n",
    "\n",
    "One method for building molecules is to copy atoms, as the interactions are also copied. For example, to build methane:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m_C = Atom('C')\n",
    "# Define the H atom positions relative to a C at [0,0,0]\n",
    "d = 0.629118\n",
    "H_pos = [[d, d, d], [-d, -d, d], [d, -d, -d], [-d, d, -d]]\n",
    "m_H = Atom('H', position=H_pos[0])\n",
    "# Add a bond between \n",
    "CH_bonds = Bond(m_C, m_H, function=HarmonicPotential(1.09, 100., interaction_type='bond'))\n",
    "# Make three H atom copies (which are therefore all bonded to m_C)\n",
    "H_atoms = [m_H]\n",
    "for pos in H_pos[1:]:\n",
    "    H_atoms.append(m_H.copy(pos))\n",
    "# The next two lines simply create a list of unique HCH triplets\n",
    "# e.g. [(m_H1, m_C, m_H2), (m_H1, m_C, m_H3), ..., (m_H3, m_C, m_H4)].\n",
    "from itertools import combinations\n",
    "HCH_triplets = [(i[0], m_C, i[1]) for i in combinations(H_atoms, 2)]\n",
    "# Unpack the list of triplets with * notation i.e. [(...), (...), (...)] becomes (...), (...), (...)\n",
    "HCH_angles = BondAngle(*HCH_triplets, function=HarmonicPotential(109.5, 10., interaction_type='angle'))\n",
    "# Create a methane molecule by adding C atom to list of H atoms\n",
    "methane = Molecule(atoms=[m_C]+H_atoms)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Add structures to a universe\n",
    "Now we have created our structure objects, we need to add them to our universe box. There are two methods for adding a structure to a universe; either `add_structure`, which adds an individual structure to the universe, or `fill`, which fills the universe with copies of the structure. To `fill`, either specify a desired density of the structure with `num_density` or a desired number of structures with `num_struc_units`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add an individual structural unit to the universe\n",
    "universe.add_structure(H_mol)\n",
    "# Fill the universe with the structural unit repeated on a cubic lattice with a specific number density\n",
    "universe.fill(H_mol, num_density = 0.01)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Currently the `fill` command cannot be used in conjunction with `add_structure`. When using it with a cubic universe the density will be isotropic, however the exact number of structural units added by `fill` may be lower than expected as it will always add a cube number. Using `fill` with a non-cubic universe is not recommended as the density may or may not be isotropic depending on the dimensions of the universe and the number of units. A list of all structures in the universe can be viewed with `universe.structure_list` (try it!)."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solving bond constraints\n",
    "When constraining bonds, [constraint algorithms](https://en.wikipedia.org/wiki/Constraint_(computational_chemistry)) such as SHAKE or RATTLE can also be passed to a universe to solve constraints imposed on bonded interactions. A constraint algorithm is required if any of the bonded interactions are constrained."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import Shake and Rattle\n",
    "from MDMC.MD import Shake, Rattle\n",
    "# Initialise a universe with the Shake algorithm\n",
    "# The first Shake parameter is the accuracy and the second is the maximum number of iterations used for any constraint calculation\n",
    "shake = Shake(1e-4, 100)\n",
    "uni4 = Universe(10., constraint_algorithm=shake)\n",
    "# Add Rattle after universe initialisation\n",
    "rattle = Rattle(1e-5, 1000)\n",
    "uni5 = Universe(10.)\n",
    "uni5.constraint_algorithm = rattle"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Not all constraint algorithms are implemented for all MD engines, and they may also required different parameters to be specified - see the MD engine documentation for more information."
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "add738a18a8c9080256f416f48e2f1975a0b0948846decec7f4bbc677823fdf5"
  },
  "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.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}