{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reading atoms from configuration files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "MDMC enables reading atomic configurations and topologies (if included) from configuration files.\n",
    "\n",
    "The implemented file formats are:\n",
    "\n",
    "- [Crystallographic Information File](https://www.iucr.org/resources/cif/documentation) (CIF)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In theory, it is possible to read a configuration file by simply passing it to the `read` function, as MDMC will determine the file format from the extension, if it is a supported format:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from MDMC.readers.configurations import read\n",
    "atoms = read('data/Paracetamol.cif')\n",
    "\n",
    "# The return type of read is a list of atoms\n",
    "print('The return type is: {}'.format(type(atoms)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can use this to create a Molecule\n",
    "from MDMC.MD import Molecule\n",
    "paracetamol = Molecule(atoms=atoms)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The paracetamol `Molecule` read from the CIF file can be visualized using the `view` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from MDMC.gui import view\n",
    "view(paracetamol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From this it appears that the configuration and topology has been correctly read; however, MDMC does not automatically determine if certain atoms or bonded interactions are equivalent.  For example, all three hydrogens in the methyl group are equivalent, as are the bonds, bond angles and dihedrals these hydrogens are involved in.  In the case of paracetamol, there are 20 bonds in total, however only 12 of them are unique."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# To see all of the Bond interactions, we can filter the paracetamol interactions by name\n",
    "bonds = list(filter(lambda x: x.name == 'Bond', paracetamol.interactions))\n",
    "\n",
    "# There are 20 bonds in total\n",
    "print('Number of bonds: {}'.format(len(bonds)))\n",
    "\n",
    "# We can cast the list of bonds to a set to see the number of unique bonds - in this case it is still 20\n",
    "unique_bonds = set(bonds)\n",
    "print('Number of unique bonds: {}'.format(len(unique_bonds)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are two methods for specifying which atoms are equivalent; these can be seen by viewing the documentation relevant to the file format of the configuartion file, which can be accessed by passing `docstring=True` when calling `read`.  For example, with a CIF file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "read('Paracetamol.cif', docstring=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So it is possible to define equivalent atoms either by specifying the `name` of each `Atom` or the `atom_type` of each `Atom`.  **These names or atom types must be passed in the same order as the atoms are defined within the CIF file (i.e. in the \\_atom_site_label identifier).**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Equivalence specified by `atom_types`\n",
    "As with all `atom_types` specified in MDMC the integers used should be monotonically increasing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# from MDMC.readers.configurations import read\n",
    "atoms = read('data/Paracetamol.cif', atom_types=[1, 2, # Oxygens\n",
    "                                           3, # Nitrogen\n",
    "                                           4, 5, 6, 7, 6, 6, 6, 8, # Carbons\n",
    "                                           9, 9, 9, 10, 10, 10, 10, 11, 12]) # Hydrogens\n",
    "paracetamol = Molecule(atoms=atoms)\n",
    "\n",
    "# Print the total number of bonds and the number of unique bonds\n",
    "bonds = list(filter(lambda x: x.name == 'Bond', paracetamol.interactions))\n",
    "print('Number of bonds: {}'.format(len(bonds)))\n",
    "print('Number of unique bonds: {}'.format(len(set(bonds))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Equivalence specified by `names`\n",
    "The `names` must be strings which are only shared by equivalent atoms.  These set the `name` attribute for each `Atom`.  As covered in the tutorial [Applying a Forcefield](applying-a-forcefield.ipynb), if each `Atom` has an `Atom.name` which is defined in a `ForceField`, the `ForceField` can be applied to set the `Parameter` values for all interactions.  An example of this is shown below, where the `names` provided are the OPLSAA atom types:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# from MDMC.readers.configurations import read\n",
    "atoms = read('data/Paracetamol.cif', names=['109', '177', # Oxygens\n",
    "                                            '207', # Nitrogen\n",
    "                                            '208', '108', '90', '178', '90', '90', '90', '185', # Carbons\n",
    "                                            '85', '85', '85', '91', '91', '91', '91', '183', '110']) # Hydrogens\n",
    "paracetamol = Molecule(atoms=atoms)\n",
    "\n",
    "# Print the total number of bonds and the number of unique bonds\n",
    "bonds = list(filter(lambda x: x.name == 'Bond', paracetamol.interactions))\n",
    "print('Number of bonds: {}'.format(len(bonds)))\n",
    "print('Number of unique bonds: {}'.format(len(set(bonds))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The paracetamol can now be added to a `Universe` and the OPLSAA `ForceField` can be applied to set the `Parameter` values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from MDMC.MD import Universe\n",
    "universe = Universe(10.)\n",
    "universe.add_structure(paracetamol)\n",
    "universe.add_force_field('OPLSAA')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(universe)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "universe.nonbonded_interactions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The effect of this can be seen if we visualize the paracetamol molecule using the `'ASE'` viewer and select `View -> Charges` from the GUI menu bar."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "view(paracetamol, viewer='ASE')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Display for paracetamol with charges](../_static/images/gui_paracetamol_charges.png)"
   ]
  }
 ],
 "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.7.10"
  },
  "nbsphinx": {
   "execute": "never"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}