notebook/Lennard_Jones.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# **Lennard-Jones Potential**\n",
    "\n",
    "<i class=\"fa fa-home fa-2x\"></i><a href=\"./index.ipynb\" style=\"font-size: 20px\"> Go back to index</a>\n",
    "\n",
    "**Source code:** https://github.com/osscar-org/materials-science/blob/master/notebook/Lennard_Jones.ipynb\n",
    "\n",
    "<hr style=\"height:1px;border:none;color:#cccccc;background-color:#cccccc;\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **Goals**\n",
    "\n",
    "* We demonstrate the Lennard-Jones potential in this notebook.\n",
    "    * Understand the interactions between particles and concept of the potentials.\n",
    "    * Know the formula of the Lennard-Jones potential.\n",
    "    * Understand the application of the Lennard-Jones potential in molecular modelling.\n",
    "    * Know the insuffisance of the Lennard-Jones potential."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **Background theory**\n",
    "\n",
    "In the molecular modelling, it is crucial to construct mathematical models for\n",
    "the interactions of the simulated particles. Primarly, the simplest model should\n",
    "be developed for two-body systems. The Lennard-Jones potential was proposed by\n",
    "Sr Lennard Jones in 1931 \n",
    "[[J E Lennard-Jones 1931 Proc. Phys. Soc. 43 461]](https://iopscience.iop.org/article/10.1088/0959-5309/43/5/301). The Lennard-Jones\n",
    "potential is a simplifed and realistic model, which has proved successfully for \n",
    "a varties of systems. Therefore, it is always used to be the first potential \n",
    "model to teach students."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<details close>\n",
    "<summary style=\"font-size: 20px\">Lennard-Jones potential</summary>\n",
    "The mathematical formula of the Lennard-Jones potential is:\n",
    "    \n",
    "$$\\large E(r)= \\epsilon \\left[\\frac{a}{r^{12}}-\\frac{b}{r^6}\\right]$$\n",
    "    \n",
    "where, $r$ is the distance between two atoms. $a$ and $b$ are two parameters.\n",
    "$\\epsilon$ is the dispersion energy, which controls the depth of the\n",
    "well potential. For a two-body system like two nuclei, the particles repulse \n",
    "each other at short distance. While particles shows attraction at long distance.\n",
    "The interaction should trend to vanlish at infinite distance.\n",
    "The Lennard-Jones potential fulfills all these three conditions.\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<details close>\n",
    "<summary style=\"font-size: 20px\">Molecular simulations with LJ potential</summary>\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"height:1px;border:none;color:#cccccc;background-color:#cccccc;\" />\n",
    "\n",
    "## Interactive visualization\n",
    "(be patient, it might take a few seconds to load)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from bqplot import pyplot as plt\n",
    "from ipywidgets import Checkbox, HBox, VBox, link, jslink, FloatSlider, Button\n",
    "\n",
    "a = 1.0\n",
    "b = 2.0\n",
    "\n",
    "fig1 = plt.figure(1, title='Lennard Jones Potential', animation_duration=50)\n",
    "plt.clear()\n",
    "n = 500\n",
    "x = np.linspace(0.7, 1.6, n)\n",
    "y = (a/x**12 - b/x**6)\n",
    "\n",
    "y1 = a/x**12\n",
    "y2 = -b/x**6\n",
    "\n",
    "axes_options = {'x': {'label': 'Distance'}, 'y': {'label': 'Energy'}}\n",
    "\n",
    "line_LJ = plt.plot(x, y, 'r', axes_options=axes_options)\n",
    "line_r = plt.plot(x, y1, 'b', line_style = 'dash_dotted')\n",
    "line_a = plt.plot(x, y2, 'y', line_style = 'dash_dotted')\n",
    "\n",
    "px = pow(2*a/b, -1/6.0)\n",
    "py = (a/px**12 - b/px**6)\n",
    "p_dist = plt.scatter(x = [px], y = [py], colors =['orange'], stroke = 'black', default_size = 100)\n",
    "#p1 = plt.scatter(x = [1.0,], y = [-1.0], colors =['orange'], stroke = 'black', default_size = 100)\n",
    "\n",
    "\n",
    "plt.ylim(-3, 8.0)\n",
    "plt.xlim(0.75, 1.5)\n",
    "\n",
    "line_r.visible = False\n",
    "line_a.visible = False\n",
    "\n",
    "cb1 = Checkbox(value = False, description = 'Repulsion Energy', disabled = False, indent = True)\n",
    "cb2 = Checkbox(value = False, description = 'Attraction Energy', disabled = False, indent = True)\n",
    "\n",
    "dis = FloatSlider(value = 1.0, min = 0.8, max = 1.5, description = 'Distance: ', step = 0.01)\n",
    "\n",
    "a_slider = FloatSlider(value = 1.0, min = 0.5, max = 10.0, description = 'Parameter a: ', step = 0.1)\n",
    "b_slider = FloatSlider(value = 2.0, min = 0.5, max = 10.0, description = 'Parameter b: ', step = 0.1)\n",
    "\n",
    "energy_button = Button(description = 'Energy Mimimum')\n",
    "\n",
    "def on_distance_change(change):\n",
    "    a = a_slider.value\n",
    "    b = b_slider.value\n",
    "    p_dist.x = [dis.value]\n",
    "    p_dist.y = [a/dis.value**12-b/dis.value**6]\n",
    "    p_dist2.x = [dis.value]\n",
    "    p_dist2.y = [-12.0*a/dis.value**13+6.0*b/dis.value**7]\n",
    "\n",
    "def on_button_clicked(b):\n",
    "    #p_dist.x = [1.0]\n",
    "    #p_dist.y = [-1.0]\n",
    "    a = a_slider.value\n",
    "    b = b_slider.value\n",
    "    dis.value = pow(2.0*a/b, 1.0/6.0)\n",
    "    \n",
    "dis.observe(on_distance_change, names = 'value')\n",
    "energy_button.on_click(on_button_clicked)\n",
    "\n",
    "jslink((cb1, 'value'), (line_r, 'visible'))\n",
    "jslink((cb2, 'value'), (line_a, 'visible'))\n",
    "\n",
    "fig2 = plt.figure(2, title='Force', animation_duration=50)\n",
    "\n",
    "plt.clear()\n",
    "nx = 500\n",
    "xx = np.linspace(0.7, 1.6, n)\n",
    "\n",
    "yx = (-12.0*a/xx**13 + 6.0*b/xx**7)\n",
    "\n",
    "yx1 = -12.0*a/xx**13\n",
    "yx2 = 6.0*b/xx**7\n",
    "\n",
    "axes_options2 = {'x': {'label': 'Distance'}, 'y': {'label': 'Force'}}\n",
    "\n",
    "l_LJ = plt.plot(xx, yx, 'r', axes_options=axes_options2)\n",
    "l_r = plt.plot(xx, yx1, 'b', line_style = 'dash_dotted')\n",
    "l_a = plt.plot(xx, yx2, 'y', line_style = 'dash_dotted')\n",
    "\n",
    "plt.ylim(-30.0, 10.0)\n",
    "plt.xlim(0.75, 1.5)\n",
    "\n",
    "cb3 = Checkbox(value = False, description = 'Repulsion Force', disabled = False, indent = True)\n",
    "cb4 = Checkbox(value = False, description = 'Attraction Force', disabled = False, indent = True)\n",
    "\n",
    "jslink((cb3, 'value'), (l_r, 'visible'))\n",
    "jslink((cb4, 'value'), (l_a, 'visible'))\n",
    "\n",
    "def on_a_changed(b):\n",
    "    a = a_slider.value\n",
    "    b = b_slider.value\n",
    "    line_LJ.y = (a/line_LJ.x**12 - b/line_LJ.x**6)\n",
    "    line_r.y = a/line_r.x**12\n",
    "    line_a.y = -b/line_a.x**6\n",
    "    l_LJ.y = (-12*a/l_LJ.x**13 + 6.0*b/l_LJ.x**7)\n",
    "    l_r.y = -12*a/l_r.x**13\n",
    "    l_a.y = 6*b/l_a.x**7\n",
    "    px = pow(2*a/b, 1.0/6.0)\n",
    "    py = (a/px**12 - b/px**6)\n",
    "    p_dist.x = [px]\n",
    "    p_dist.y = [py]\n",
    "    px2 = pow(2*a/b, 1.0/6.0)\n",
    "    py2 = (-12*a/px2**13 + 6*b/px2**7)\n",
    "    p_dist2.x = [px2]\n",
    "    p_dist2.y = [py2]\n",
    "    dis.value = px2\n",
    "    \n",
    "def on_b_changed(b):\n",
    "    a = a_slider.value\n",
    "    b = b_slider.value\n",
    "    line_LJ.y = (a/line_LJ.x**12 - b/line_LJ.x**6)\n",
    "    line_r.y = a/line_r.x**12\n",
    "    line_a.y = -b/line_a.x**6\n",
    "    l_LJ.y = (-12*a/l_LJ.x**13 + 6.0*b/l_LJ.x**7)\n",
    "    l_r.y = -12*a/l_r.x**13\n",
    "    l_a.y = 6*b/l_a.x**7\n",
    "    px = pow(2*a/b, 1.0/6.0)\n",
    "    py = (a/px**12 - b/px**6)\n",
    "    p_dist.x = [px]\n",
    "    p_dist.y = [py]\n",
    "    px2 = pow(2*a/b, 1.0/6.0)\n",
    "    py2 = (-12*a/px2**13 + 6*b/px2**7)\n",
    "    p_dist2.x = [px2]\n",
    "    p_dist2.y = [py2]\n",
    "    dis.value = px2\n",
    "\n",
    "a_slider.observe(on_a_changed, names = 'value');\n",
    "b_slider.observe(on_b_changed, names = 'value');\n",
    "\n",
    "px2 = pow(2*a/b, 1.0/6.0)\n",
    "py2 = (-12*a/px2**13 + 6*b/px2**7)\n",
    "p_dist2 = plt.scatter(x = [px2], y = [py2], colors =['orange'], stroke = 'black', default_size = 100)\n",
    "#p2 = plt.scatter(x = [1.0,], y = [0.0], colors =['orange'], stroke = 'black', default_size = 100)\n",
    "\n",
    "fig1.layout.min_width = \"400px\"\n",
    "fig2.layout.min_width = \"400px\"\n",
    "\n",
    "display(HBox([VBox([fig1, cb1, cb2]), VBox([fig2, cb3, cb4])]))\n",
    "display(a_slider, b_slider, HBox([dis, energy_button]))"
   ]
  },
  {
   "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}