{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Argon A-to-Z"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tutorial demonstrates a-to-z how to optimise Lennard Jones parameters for liquid argon, and without going into details. For details see other tutorials and wider MDMC documentation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Imports used for this tutorial\n",
    "import numpy as np\n",
    "import os\n",
    "from MDMC.control import Control\n",
    "from MDMC.MD import Atom, Dispersion, LennardJones, Simulation, Universe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Change the number of threads depending on the number of physical cores on your computer\n",
    "# as it was tested for LAMMPS\n",
    "os.environ[\"OMP_NUM_THREADS\"] = \"4\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Build universe with density 0.0176 atoms per AA^-3\n",
    "density = 0.0176\n",
    "# This means cubic universe of side:\n",
    "# 23.0668 A will contain 216 Ar atoms\n",
    "# 26.911 A will contain 343 Ar atoms\n",
    "# 30.7553 A will contain 512 Ar atoms\n",
    "# 38.4441 A will contain 1000 Ar atoms\n",
    "universe = Universe(dimensions=23.0668)\n",
    "Ar = Atom('Ar', charge=0., mass=36.0)\n",
    "# Calculating number of Ar atoms needed to obtain density\n",
    "n_ar_atoms = int(density * np.product(universe.dimensions))\n",
    "print(f'Number of argon atoms = {n_ar_atoms}')\n",
    "universe.fill(Ar, num_struc_units=(n_ar_atoms))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the Jupyter cell above, a box of Argon atoms is set up. However, at this point there is no interaction forces between the argon atoms! In the cell below an appropriate (for argon) force-field interaction potential is defined."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ar_dispersion = Dispersion(universe,\n",
    "                           (Ar.atom_type, Ar.atom_type),\n",
    "                           cutoff=8.,\n",
    "                           function=LennardJones(epsilon=1.0243, sigma=3.36))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case the interaction potential chosen is the humble Lennard Jones (to get info see doc or type `help(LennardJones)`).\n",
    "\n",
    "Also, a `cutoff` value is chosen (see `help(Dispersion)` for more info). A [rule of thumb for Lennard-Jones](https://en.wikipedia.org/wiki/Lennard-Jones_potential) is to pick `cutoff=2.5*sigma`. The value for argon is recommended to be between 8 and 12 ang. `cutoff` is not a force-field parameter and therefore will not be refined. Ideally, and for any system you want to pick at value of the `cutoff` which is small while not compromising accuracy. For this system picking a value between 8 and 12 ang is found to give near identifical results.\n",
    "\n",
    "Next (and before starting the refinement), we set up the MD engine and equilibrate the system. Note with MDMC the equilibration only needs to be done once. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MD Engine setup\n",
    "simulation = Simulation(universe,\n",
    "                        engine=\"lammps\",\n",
    "                        time_step=10.18893,\n",
    "                        temperature=120.,\n",
    "                        traj_step=15)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Energy Minimization and equilibration\n",
    "simulation.minimize(n_steps=50)\n",
    "simulation.run(n_steps=10000, equilibration=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK; time to set up the actual refinement of the force-field parameters. \n",
    "\n",
    "First we need some data to refine against:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# exp_datasets is a list of dictionaries with one dictionary per experimental\n",
    "# dataset\n",
    "# Dataset from: van Well et al. (1985). Physical Review A, 31(5), 3391-3414\n",
    "# resolution is None as the original author already accounted for instrument resolution\n",
    "exp_datasets = [{'file_name':'data/Well_s_q_omega_Ar_data.xml',\n",
    "                 'type':'SQw',\n",
    "                 'reader':'xml_SQw',\n",
    "                 'weight':1.,\n",
    "                 'auto_scale':True,\n",
    "                 'resolution':800}]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The number of `MD_steps` specified must be large enough to allow for successful calculation of all observables. This depends the `type` of the dataset provided and the value of the `traj_step` (specified when creating the `Simulation`). If a value for `MD_steps` is not provided, then the minimum number needed will be used automatically.\n",
    "\n",
    "Additionally, some observables will have an upper limit on the number of MD_steps that can be used in calculating their dependent variable(s). In these cases, the number of `MD_steps` is rounded down to a multiple of this upper limit so that we only run steps that will be useful. For example, if we use 1000 `MD_steps` in calculation, but a value of 2500 is provided, then we will run 2000 steps and use this to calculate the variable twice, without wasting time performing an additional 500 steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fit_parameters = universe.parameters\n",
    "fit_parameters['sigma'].constraints = [2.8,3.8]\n",
    "fit_parameters['epsilon'].constraints = [0.6, 1.4]\n",
    "\n",
    "\n",
    "control = Control(simulation=simulation,\n",
    "                  exp_datasets=exp_datasets,\n",
    "                  fit_parameters=fit_parameters,\n",
    "                  minimizer_type=\"GPO\",\n",
    "                  reset_config=True,\n",
    "                  MD_steps=4000, \n",
    "                  equilibration_steps=4000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And finally start the refinement! Bump up `n_steps` from 3 when you are ready."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run the refinement, i.e. refine the FF parameters against the data\n",
    "control.refine(n_steps=25)\n",
    "control.plot_results();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.6 64-bit",
   "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.2"
  },
  "vscode": {
   "interpreter": {
    "hash": "949777d72b0d2535278d3dc13498b2535136f6dfe0678499012e853ee9abcab1"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}