{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Applying a `ForceField`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As well as defining `Parameter` values manually, it is also possible to apply a `ForceField` (e.g. OPLSAA).\n",
    "\n",
    "To apply a `ForceField`, it is necessary to specify the `ForceField` `atom_type`.  **It is important to note that this is different from the `Atom.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 `ForceField`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from MDMC.MD import Atom, Bond, BondAngle, DihedralAngle, Dispersion, Molecule, Universe\n",
    "\n",
    "# Alternatively, could simply import all (as shown below).\n",
    "#from MDMC.MD import *\n",
    "\n",
    "# Define the unique atoms using the ForceField atom_type\n",
    "# These can be seen in the oplsaa.dat file (MDMC/MD/force_fields/data/oplsaa.dat)\n",
    "# The H1 atom will be copied after the bond and bond angles have been defined\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)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All bonded and nonbonded interactions need to be specified as normal, except the `function` should not be passed as a parameter.  This is because the `ForceField` provides the `InteractionFunction` as well as the `Parameter` values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 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 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",
    "# 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))"
   ]
  },
  {
   "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.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}