{
 "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, Molecule, Universe\n",
    "from MDMC.MD.force_fields.OPLSAA import add_opls_force_field"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first create our methanol molecule.\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\")\n",
    "HC2 = Atom(\"H\", position=[-0.7006,  0.3636, -0.8900], name=\"98\")\n",
    "HC3 = Atom(\"H\", position=[-0.7076, -1.1754,  0.0000], name=\"98\")\n",
    "C = Atom(\"C\", position=[-0.3366, -0.1504,  0.0000], name=\"99\")\n",
    "O = Atom(\"O\", position=[ 1.0849, -0.1713,  0.0000], name=\"96\")\n",
    "HO = Atom(\"H\", position=[ 1.3606,  0.7699,  0.0000], name=\"97\")\n",
    "\n",
    "# Create the methanol Molecule\n",
    "methanol = Molecule(\n",
    "    atoms=[HC1, HC2, HC3, C, O, HO],\n",
    "    interactions=[\n",
    "        # Create the bonds with harmonic potentials\n",
    "        Bond((C, HC1), (C, HC2), (C, HC3)),\n",
    "        Bond((O, HO)),\n",
    "        Bond((C, O), constrained=True),\n",
    "        # Create the H-C-O bond angles\n",
    "        BondAngle((HC1, C, O), (HC2, C, O), (HC3, C, O)),\n",
    "        # Create an HCH bond angle\n",
    "        BondAngle((HC1, C, HC2), (HC2, C, HC3), (HC3, C, HC1)),\n",
    "        # Create the H-O-C bond angle\n",
    "        BondAngle((HO, O, C)),\n",
    "        # Create the H-C-O-H dihedral\n",
    "        DihedralAngle((HC1, C, O, HO), (HC2, C, O, HO), (HC3, C, O, HO))\n",
    "    ]\n",
    ")\n",
    "\n",
    "# Create a universe and add the methanol\n",
    "universe = Universe(dimensions=15.0)\n",
    "universe.fill(methanol, num_density=0.01)"
   ]
  },
  {
   "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": [
    "add_opls_force_field(universe, cutoff=6.0, ewald=1e-4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This sets all of the `InteractionFunction` and `Parameter` values for each `Interaction`."
   ]
  },
  {
   "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.force_fields.OPLSAA 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.11.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}