{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Building an MDMC Universe\n",
    "An MDMC simulation requires a configuration and a topology defined within a `Universe` object. Although MDMC is somewhat flexible with regards to the order in which these are created, a suggested approach is as follows:\n",
    "\n",
    "*Restriction on creating a* `Universe` *object*: MDMC currently only supports orthorhombic boxed Universes; the `dimensions` must be specified when creating (initialising) the `Universe` object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import the Universe class\n",
    "from MDMC.MD import Universe\n",
    "# Initialise a Universe with dimensions in Ang\n",
    "universe = Universe([10.0, 15.0, 20.0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To create a cubic `Universe`, only a single float needs to be specified for the `dimensions`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "universe = Universe(10.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create an atomic configuration\n",
    "Configurations can either be specified by the user or read from a CIF file.\n",
    "\n",
    "### Create an atom\n",
    "Each atom is specified using an `Atom` object, which possesses (amongst other things) a `position`, `velocity`, `element`, `mass`, `charge`, and `atom_type`.  At a minimum the `element` must be specified when creating an `Atom` object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import the Atom class\n",
    "from MDMC.MD import Atom\n",
    "# Create a hydrogen atom\n",
    "H1 = Atom('H')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This will create a Hydrogen atom with `mass` determined from an elemental lookup table, `position` set to the origin, `velocity` set to (0., 0., 0.), no `charge` and no `atom_type`. Note that if set, the velocity of atoms in the MD engine 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 passing the position of the new atom:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H2 = H1.copy(position=[1., 1., 1.])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The copied atom will have identical properties to the original attribute, except with a different ID (which is unique for all structural units) and 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'**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Atom type\n",
    "Each atom has an `atom_type`, which is used for applying non-bonded interactions between atoms; all atoms of the same `atom_type` will have the same non-bonded interactions applying to them.  The `atom_type` can be specified when the atom is created:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "O1 = Atom('O', position=[8.0, 8.0, 8.0], atom_type=1)\n",
    "# and to add the created atom to universe do\n",
    "universe.add_structure(O1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, the `atom_type` can be inferred by MDMC.  This occurs when an atom with no defined `atom_type` is added to a Universe:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "O2 = Atom('O', position=[9.0, 9.0, 9.0])\n",
    "universe.add_structure(O2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assuming no other atoms have been added to the universe, this sets ```O2.atom_type``` equal to 1.  MDMC infers the `atom_type` of each atom based on its element and its interactions: all atoms with the same element and the same interactions **when they are added to the universe** will have the same `atom_type`.  **Once an atom's** `atom_type` **is set, it is immutable i.e. it cannot be changed.**\n",
    "It is recommended that either all atom types are specified when atoms are created or none are (i.e. either MDMC is left to infer and assign all atom types or no atom types.  Only specifying some atom types could result in unexpected behaviour; if it is absolutely necessary then ensure that all atoms with specified atom types are added to the universe first).\n",
    "\n",
    "To see what atoms have been added to the universe:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "universe.structure_list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see other attributes of the universe (or any other MDMC containers) with:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "help(universe)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The same applies for methods, for example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "help(universe.add_structure)"
   ]
  },
  {
   "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`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Bond`\n",
    "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`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A `BondAngle` 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`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All `Bond` and `BondAngle` can be constrained; 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)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a constraint to be applied during a simulation, the `Universe` must have a `ConstraintAlgorithm`."
   ]
  },
  {
   "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."
   ]
  },
  {
   "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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create a molecule\n",
    "A molecule consists of two or more atoms and at least one bonded interaction:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import the Atom, Molecule, Bond, BondAngle and HarmonicPotential classes\n",
    "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)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Creating molecules using configuration files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition to creating molecules manually, it is also possible to create them from atomic configuration files (e.g CIF files).  Please see the tutorial [Reading atoms from configuration files](read-configurations.ipynb) for a detailed description on how to do this."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Add structural units to a universe\n",
    "There are two methods for adding a structural unit to a universe:"
   ]
  },
  {
   "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)"
   ]
  },
  {
   "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.\n",
    "\n",
    "### Create non-bonded interactions\n",
    "Non-bonded interactions are applied to atoms based on their `atom_type`, rather than to individual atoms.  They must also have a `Universe` specified, so that they know to which atoms they apply (atoms must have an `atom_type` once they have been added to a universe).  For example, to create a `Dispersion` interaction with a `LennardJones` interaction function between two atom types:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import the Dispersion interaction and Lennard-Jones function\n",
    "from MDMC.MD import Dispersion\n",
    "from MDMC.MD import LennardJones\n",
    "# Create a dispersion interaction with a Lennard-Jones function between atoms with atom_type 1 (O) and atom_type 2 (H)\n",
    "# The first LJ parameter is epsilon and the second is sigma\n",
    "# The cutoff for the dispersion interaction can also be set (in Ang)\n",
    "LJ_HO = Dispersion(universe, [1, 2], function=LennardJones(0.65, 3.), cutoff=10.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The exception to this is `Coulombic` interactions, which can be applied either to a list of atoms or to a list of atom types.  If the Coulombic interaction is applied to a list of atoms, the universe does not need to be specified:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import the Coulombic interaction\n",
    "from MDMC.MD import Coulombic\n",
    "# Create a Coulombic interaction with a Coulomb potential and a charge of 0.42\n",
    "c_H = Coulombic(atoms=[H1, H2], charge=0.42)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As `Coulombic` interactions typically have the same interaction function (i.e. a `Coulomb` function, where the force results in Coulomb's law), `Coulombic` interactions do not need to be specified with an interaction function (although other functions can be provided); to set the function as `Coulomb`, a value for the charge can be passed.  The warning highlights the `Coulomb` interaction function has been automatically set.  This is equivalent to:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from MDMC.MD import Coulomb\n",
    "c_H = Coulombic(atoms=[H1, H2], function=Coulomb(0.42))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As with all non-bonded interactions, a Coulombic interaction can also be created by specifying the atom types:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c_O = Coulombic(universe, atom_types=[1], charge=0.42)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case a universe must be provided."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Applying a `ForceField`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As well as manually setting the `Parameter` values for each interaction (e.g. `Coulombic(atoms=[H1, H2], function=Coulomb(0.42))` to set the `charge`), it is also possible to apply a `ForceField` to a `Universe`, in order to set the `Parameter` values.  Please see the tutorial [Applying a ForceField](applying-a-forcefield.ipynb) for a detailed description on how to do this."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adding kspace solvers to a universe\n",
    "There are several solvers for determining the long range energy contribution for non-bonded interactions, including Ewald, particle-particle particle-mesh (PPPM), and particle-mesh Ewald (PME).  If the long range contribution to the non-bonded energy (i.e. > cutoff) is to be calculated during a simulation, a kspace solver has to be added to the universe.  A solver can be specified for either electrostatic or dispersive interactions, or both. This can either be during universe initialisation or afterwards:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import Ewald and PPPM kspace solvers\n",
    "from MDMC.MD import Ewald, PPPM\n",
    "# Create an ewald solver\n",
    "ewald = Ewald(accuracy=1e-5)\n",
    "# Initialise a universe with an Ewald solver for both electrostatics and dispersive interactions\n",
    "uni1 = Universe(10., kspace_solver=ewald)\n",
    "# Initialise a universe and then add a PPPM solver for electrostatic interactions\n",
    "uni2 = Universe(10.)\n",
    "pppm = PPPM(accuracy=1e-4)\n",
    "uni2.electrostatic_solver = pppm\n",
    "# Initialise a universe with a PPPM solver for dispersive interactions\n",
    "uni3 = Universe(10., dispersive_solver=PPPM)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Not all kspace solvers are implemented for all MD engines, and they may also require different parameters to be specified - see the MD engine documentation for more information.\n",
    "\n",
    "### Adding constraint algorithms to a universe\n",
    "In a similar manner to kspace solvers, constraint algorithms (e.g. SHAKE, RATTLE) can also be passed to a universe.  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"
   ]
  },
  {
   "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.\n",
    "\n",
    "## Example Universe filled with water\n",
    "For water, SPCE, SPC, TIP3P and TIP3PFB forcefields are predefined, so parameters do not need to be set when interactions are defined; instead parameters are set by adding a forcefield to the universe.\n",
    "\n",
    "Note that all water models have (harmonic) potential strengths defined for their `Bond`s and `BondAngle`s. In order to create a rigid water molecule in accordance with the models, a `constraint_algorithm` should be passed to the `Universe` and `constrained=True` passed to the relevant interactions, as shown below.\n",
    "\n",
    "\"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](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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from MDMC.MD import *\n",
    "\n",
    "universe = Universe(dimensions=21.75, constraint_algorithm=Shake(1e-4, 100), electrostatic_solver=PPPM(accuracy=1e-5))\n",
    "H1 = Atom('H')\n",
    "H2 = Atom('H', position=(0., 1.63298, 0.))\n",
    "O = Atom('O', position=(0., 0.81649, 0.57736))\n",
    "H_coulombic = Coulombic(atoms=[H1, H2], cutoff=10.)\n",
    "O_coulombic = Coulombic(atoms=O, cutoff=10.)\n",
    "water_mol = Molecule(position=(0, 0, 0),\n",
    "                     velocity=(0, 0, 0),\n",
    "                     atoms=[H1, H2, O],\n",
    "                     interactions=[Bond((H1, O), (H2, O), constrained=True),\n",
    "                                   BondAngle(H1, O, H2, constrained=True)],\n",
    "                     name='water')\n",
    "universe.fill(water_mol, num_density=0.03356718472021752)\n",
    "O_dispersion = Dispersion(universe, [O.atom_type, O.atom_type], cutoff=10., vdw_tail_correction=True)\n",
    "universe.add_force_field('SPCE')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Accessing information\n",
    "After creating a `Universe`, the various properties and attributes of it can be accessed as follows:\n",
    "\n",
    "### Geometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('The dimension(s) of the Universe are {0}, giving a volume of {1}'\n",
    "      ''.format(universe.dimensions, universe.volume))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Forcefield"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('The current force field applied to the Universe is:\\n{}'\n",
    "      ''.format(universe.force_fields))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Configuration of atoms and molecules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('There are {0} atoms in the Universe\\n'\n",
    "      ''.format(universe.n_atoms))\n",
    "\n",
    "# For a complete list run `print(universe.atoms)`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('There are {0} molecules in the Universe\\n'\n",
    "      ''.format(universe.n_molecules))\n",
    "\n",
    "# For a complete list run `print(universe.molecules)`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Configuration of interactions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('There are {0} total interactions in the Universe\\n'\n",
    "      ''.format(universe.n_interactions))\n",
    "\n",
    "# For a complete list run `print(universe.interactions)`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('There are {0} bonded interactions in the Universe\\n'\n",
    "      ''.format(universe.n_bonded))\n",
    "\n",
    "# For a complete list run `print(universe.bonded_interactions)`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('There are {0} nonbonded interactions in the Universe\\n'\n",
    "      ''.format(universe.n_nonbonded))\n",
    "\n",
    "# For a complete list run `print(universe.nonbonded_interactions)`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Constraints"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('The current constraint algorithm applied to the Universe is {}'\n",
    "      ''.format(universe.constraint_algorithm.name))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solvers\n",
    "Note that a `kspace_solver` is mutually excusive with the other solver types."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('The kspace_solver is {}'.format(universe.kspace_solver and universe.kspace_solver.name))\n",
    "print('The electrostatic_solver is {}'.format(universe.electrostatic_solver and universe.electrostatic_solver.name))\n",
    "print('The dispersive_solver is {}'.format(universe.dispersive_solver and universe.dispersive_solver.name))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('The MDMC parameters of the Universe are:\\n{}'\n",
    "      ''.format(universe.parameters))"
   ]
  }
 ],
 "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
}