{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Applying a force field"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As well as defining interactions by parameters (e.g. the Lennard-Jones potential), it is also possible to define atomic interaction via a [force field](https://en.wikipedia.org/wiki/Force_field_(chemistry)). These define the interactions between atoms based on existing experimental or calculated data rather than calculating them during the simulation, and are an attractive choice for accurate simulation.\n",
    "\n",
    "To apply a `ForceField`, it is necessary to specify the `ForceField`'s `atom_type`.  It is important to note that this is different from the `Atom`'s `atom_type` which every atom possesses, regardless of whether or not a `ForceField` is used.\n",
    "\n",
    "This can be done by either passing the `ForceField` `atom_type` or the `ForceField` atom `name` as the `Atom.name`.\n",
    "\n",
    "It is also important to note that if you require the charge to be set from the `ForceField`, you should pass this when creating each `Atom`. Any float can be passed as this will be changed when the `ForceField` is applied.\n",
    "\n",
    "Below is an example using methanol and the [OPLSAA](https://en.wikipedia.org/wiki/OPLS) `ForceField`. Any other force field has an analogous procedure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from MDMC.MD import Atom, Bond, BondAngle, DihedralAngle, Dispersion, Molecule, Universe"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first create our methanol molecule. All bonded and nonbonded interactions need to be specified as normal, except the `function` for the `Dispersion` should not be passed as a parameter.  This is because the `ForceField` provides the interaction function as well as the parameter values.\n",
    "\n",
    "Note that we are setting the `name` parameter of the atoms to match the 'atom type' IDs used in OPLSAA for the same atoms. This is important, as MDMC will use these names to define the force field interactions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the atoms\n",
    "HC1 = Atom('H', position=[-0.7006,  0.3636,  0.8900], name='98', charge=0., atom_type=1)\n",
    "C = Atom('C', position=[-0.3366, -0.1504,  0.0000], name='99', charge=0., atom_type=2)\n",
    "O = Atom('O', position=[ 1.0849, -0.1713,  0.0000], name='96', charge=0., atom_type=3)\n",
    "HO = Atom('H', position=[ 1.3606,  0.7699,  0.0000], name='97', charge=0., atom_type=4)\n",
    "\n",
    "# Create the bonds with harmonic potentials\n",
    "CH_bond = Bond(C, HC1)\n",
    "CO_bond = Bond(C, O)\n",
    "OH_bond = Bond(O, HO)\n",
    "\n",
    "# Create the H-C-O and H-O-C bond angles\n",
    "HCO_angle = BondAngle((HC1, C, O))\n",
    "HOC_angle = BondAngle((HO, O, C))\n",
    "\n",
    "# Create the H-C-O-H dihedral\n",
    "HCOH_dihedral = DihedralAngle((HC1, C, O, HO))\n",
    "\n",
    "# Duplicate the HC1 atom\n",
    "HC2 = HC1.copy(position=[-0.7006,  0.3636, -0.8900])\n",
    "\n",
    "# Create an HCH bond angle\n",
    "HCH_angle = BondAngle((HC1, C, HC2))\n",
    "\n",
    "# Duplicate the HC1 atom again\n",
    "# This atom will have all existing bond (CH_bond) and bond angles (HCO_angle and HCH_angle) defined\n",
    "HC3 = HC1.copy(position=[-0.7076, -1.1754,  0.0000])\n",
    "\n",
    "# remember to add a H1-C-H3 angle! as HC3 is a copy of HC1, it only has HCH_angle\n",
    "# defined between itself and HC2\n",
    "H2CH3_angle = BondAngle((HC1, C, HC3))\n",
    "\n",
    "# Create the methanol Molecule\n",
    "methanol = Molecule(atoms=[HC1, HC2, HC3, C, O, HO])\n",
    "\n",
    "# Create a universe and add the methanol\n",
    "universe = Universe(10.)\n",
    "universe.add_structure(methanol)\n",
    "\n",
    "# Add dispersion interactions (no dispersion for HO)\n",
    "HC_disp = Dispersion(universe, (1, 1))\n",
    "C_disp = Dispersion(universe, (2, 2))\n",
    "O_disp = Dispersion(universe, (3, 3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that all of the interactions have been defined and the methanol has been added to a `Universe`, a `ForceField` can be applied to the `Universe`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "universe.add_force_field('OPLSAA')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This sets all of the `InteractionFunction` and `Parameter` values for each `Interaction`.  For example, the charges of all atoms are set:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('\\n'.join(['{0} {1} charge = {2}'.format(atom.element, atom.name, atom.charge) for atom in methanol.atoms]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly, the `InteractionFunction` and `Parameter` values are set for `BondedInteraction`. The `HCO_angle` now has a `HarmonicPotential` `function`, which has two `Parameter` values, `equilibrium,_state` set to `109.5 deg`, and `potential_strength` set to `146.44 kJ / mol deg ^ 2`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "HCO_angle"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is worth noting that applying a `ForceField` changes each `Atom.name` from the `ForceField` `atom_type` to the `ForceField` atom `name`.  This can be seen when the charges are printed above, but to emphasize the point:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H_atom = Atom('H', position=[-0.7006,  0.3636,  0.8900], name='98', charge=0., atom_type=1)\n",
    "print('Before OPLSAA applied, H_atom name is: {}'.format(H_atom.name))\n",
    "\n",
    "uni = Universe(10.)\n",
    "uni.add_structure(H_atom)\n",
    "uni.add_force_field('OPLSAA')\n",
    "\n",
    "print('After OPLSAA applied, H_atom name is: {}'.format(H_atom.name))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If these `Atom.name` were originally set to the `Forcefield` atom `name` there is no effect:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H_atom = Atom('H', position=[-0.7006,  0.3636,  0.8900], name='Methanol CH3-', charge=0., atom_type=1)\n",
    "print('Before OPLSAA applied, H_atom name is: {}'.format(H_atom.name))\n",
    "\n",
    "uni = Universe(10.)\n",
    "uni.add_structure(H_atom)\n",
    "uni.add_force_field('OPLSAA')\n",
    "\n",
    "print('After OPLSAA applied, H_atom name is: {}'.format(H_atom.name))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Determining the `ForceField` `atom_type` or atom `name`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To determine the correct `ForceField` `atom_type` or `name` for each `Atom`, there are two methods:\n",
    "\n",
    "- Search through the .dat file for the `ForceField` (MDMC/MD/force_fields/data/oplsaa.dat)\n",
    "- Import the `ForceField` and use `ForceField.filter_element`\n",
    "\n",
    "The latter of these methods is shown below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from MDMC.MD import OPLSAA\n",
    "oplsaa = OPLSAA()\n",
    "\n",
    "# If we wanted to determine the correct type of a Chlorine atom\n",
    "chlorines = oplsaa.filter_element('Cl')\n",
    "print(chlorines.to_string())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So if the Cl atom is part of a chloroalkene, the correct atom could be created with either of the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specifying the ForceField atom_type:\n",
    "Cl = Atom('Cl', name='168', charge=0.)\n",
    "\n",
    "# Is equivalent to specifying the ForceField atom name\n",
    "Cl = Atom('Cl', name='Chloroalkene Cl-CH=', charge=0.)"
   ]
  }
 ],
 "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}