notebook/molecular-dynamics/verlet_integration.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# **Verlet Integration**\n",
    "**Authors**: Dou Du\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/quantum-mechanics/blob/master/notebook/molecular-dynamics/verlet_integration.ipynb\n",
    "\n",
    "<hr style=\"height:1px;border:none;color:#cccccc;background-color:#cccccc;\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Goals\n",
    "* Background theory: familiarize yourself with the main features of the Verlet integration scheme employed in molecular dynamics simulations.\n",
    "* Background theory: acquaint yourself with the finite-differences approach to solving differential equations governing motion of a system.\n",
    "* Appreciate the dependence of the accuracy of the Verlet integration scheme on the size of the timestep used in the method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Background theory\n",
    "[More on the background theory.](./theory/theory_verlet_integration.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tasks and exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Choose a small time step (say $\\Delta t=1$), observe how the earth orbits the sun, and check how \n",
    "the total energy and angular momentum evolves in time.\n",
    "\n",
    "    <details>\n",
    "    <summary style=\"color: red\">Solution</summary>\n",
    "    The earth is moving in an elliptical orbit. The earth moves slowly at the end of the ellipse farthest from the sun, whilst it moves quickest at the end of the ellipse nearest the sun. However, the total energy \n",
    "    and angular momentum remain constant. With a sufficiently small timestep, the Earth shall remain in this orbit for a very long time.\n",
    "    </details>\n",
    "    \n",
    "2. Does the earth still remain in the same orbit when a large time step is used? Try to explain your observations.\n",
    "\n",
    "    <details>\n",
    "    <summary style=\"color: red\">Solution</summary>\n",
    "    When choosing a large time step, the earth no longer maintains the same \n",
    "    elliptical orbital. You can observe that the ellipses traced out by the Earth shift during each circuit around the sun (subsequent orbits are no longer superimposed on the first). The reason for this discrepancy can intuitively be traced back to a deficiency of the Verlet method that becomes more pronounced the larger the timestep becomes. In the Verlet method, we approximate the velocities as being constant over a short time interval $\\Delta t$. However, this is of course not true in reality. When \n",
    "    using an overly large time step, significant accumulative errors will be introduced into the \n",
    "    calculations.\n",
    "    </details>\n",
    "    \n",
    "3. In molecular dynamics simulations, what are the problems associated with choosing (1) a very small timestep, (2) a very large timestep.\n",
    "\n",
    "    <details>\n",
    "    <summary style=\"color: red\">Solution</summary>\n",
    "    As we have seen in the previous exercises, choosing a small timestep in MD simulations can improve the \n",
    "    accuracy of the results. However, a small timestep means taking a \n",
    "    larger amount of computing time to obtain the same total simulation time. Practically, we cannot \n",
    "    afford a very small timestep for the MD simulations.\n",
    "    On the other hand, a large time step will introduce large errors. \n",
    "    For example, for too large a timestep, particles can assume configurations in which they are too close to each other (in the \n",
    "    real world, the nuclear-nuclear repulsion energy is considerable at short distances and prohibits this). Hence, it is crucial to choose a reasonable compromise for our \n",
    "    timestep in MD simulations so that they are sufficiently accurate but simultaneously do not require a prohibitively large amount of computer time.\n",
    "    </details>\n",
    "    \n",
    "<hr style=\"height:1px;border:none;color:#cccccc;background-color:#cccccc;\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interactive visualization\n",
    "(be patient, it might take a few seconds to load)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib widget\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.animation import FuncAnimation\n",
    "from ipywidgets import Button, IntSlider, Output, VBox, HBox, Layout\n",
    "from numpy.linalg import norm\n",
    "import matplotlib.gridspec as gridspec\n",
    "from datetime import datetime"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_force(r1, r2, m1, m2):\n",
    "    \"\"\"\"Computer the force for the earth.\"\"\"\n",
    "    x = G*m1*m2/np.linalg.norm(r1-r2)**2\n",
    "    v = (r2-r1)/np.linalg.norm(r2-r1)\n",
    "    return x*v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_acceleration(f, m):\n",
    "    \"\"\"Compute the acceleration.\"\"\"\n",
    "    return f/m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def verlet_update(r1, r2, v, m1, m2, dt, a):\n",
    "    \"\"\"Update the position, velocity and acceleration by the Verlet algorithm.\"\"\"\n",
    "    \n",
    "    rn = 2.0*r1[-1] - r1[-2]  + dt*dt*a\n",
    "    vn = (rn - r1[-2])/(2*dt)\n",
    "    \n",
    "    r1[-2] = r1[-1]\n",
    "    r1[-1] = rn\n",
    "    \n",
    "    f = compute_force(rn, r2, m1, m2)\n",
    "    a = compute_acceleration(f, m1)\n",
    "    \n",
    "    return r1, vn, a\n",
    "\n",
    "def compute_energies(r1, r2, m1, m2, v):\n",
    "    \"\"\"Compute the potential energy.\"\"\"\n",
    "    pot = -G*m1*m2/np.linalg.norm(r1-r2)\n",
    "    kint = 0.5*m1*np.linalg.norm(v)**2\n",
    "    angmom = np.cross(r1-r2, m1*v)\n",
    "    \n",
    "    return pot, kint, angmom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m1 = 1.0\n",
    "m2 = 1000.0\n",
    "G = 0.2\n",
    "v = [0.0, 3.0]\n",
    "dt = 0.01\n",
    "\n",
    "# initial position for the sun and earth, the sun is fix in the centre\n",
    "r1 = np.array([[10, 0],[10.0, 3.0*dt]])\n",
    "r2 = np.array([0, 0])\n",
    "\n",
    "# initial force and acceleration\n",
    "f = compute_force(r1[-1], r2, m1, m2)\n",
    "a = compute_acceleration(f, m1)\n",
    "\n",
    "# the initial total energy\n",
    "etot0 = -G*m1*m2/np.linalg.norm(r1[0]-r2) + 0.5*m1*np.linalg.norm(v)**2\n",
    "\n",
    "output = Output()\n",
    "\n",
    "with output:\n",
    "    \"\"\"Output the figures\"\"\"\n",
    "    global fig, ax1, ax2, ax3, sun, earth, origin, avector, fvector, trace, title\n",
    "    fig = plt.figure(tight_layout=True, figsize=(7.8, 5))\n",
    "    fig.canvas.header_visible = False\n",
    "    gs = gridspec.GridSpec(9, 5)\n",
    "    \n",
    "    ax1 = fig.add_subplot(gs[:, 0:3])\n",
    "    ax2 = fig.add_subplot(gs[0:4, 3:5])\n",
    "    ax3 = fig.add_subplot(gs[5:9, 3:5])\n",
    "\n",
    "    sun = ax1.scatter([0.0],[0.0], s=23**2, marker=\"o\", linewidth=0, color='gold', label = 'Sun')\n",
    "    earth, = ax1.plot([10.0],[0.0], 'r.', alpha = 1.0, markersize=15, label = 'Earth')\n",
    "    origin = ax1.plot([10.0],[0.0], 'kx', alpha = 1.0, markersize=15, label = 'Origin point')\n",
    "    avector = ax1.quiver([10.0], [0], [0], [0], color = ['blue'], units='inches', scale = 10.0, label=\"a\")\n",
    "    fvector = ax1.quiver([10.0], [0], [0], [0], color = ['gold'], units='inches', scale = 10.0, label=\"v\")\n",
    "    trace, = ax1.plot([10.0], [0.0], \"k\", linewidth=0.5)\n",
    "    title = ax1.set_title(\"\")\n",
    "\n",
    "    ax1.legend(loc=2, fontsize=10)\n",
    "\n",
    "    ax1.set_xlim(-11.0, 11.0)\n",
    "    ax1.set_ylim(-11.0, 11.0)\n",
    "    ax1.grid()\n",
    "    \n",
    "    lpot, = ax2.plot([], [], 'k-', label='pot')\n",
    "    lkint, = ax2.plot([], [], 'c-', label='kint')\n",
    "    ltot, = ax2.plot([],[], 'r-', label='tot')\n",
    "    ax2.axhline(etot0, color='blue', ls='--', linewidth=0.5)\n",
    "    \n",
    "    ax2.set_xlim(0, 20.0)\n",
    "    ax2.set_ylim(-95, 60)\n",
    "    ax2.legend(loc=3, ncol=3)\n",
    "    ax2.set_title(\"Energies as a function of time\", fontsize = 10)\n",
    "    \n",
    "    lang, = ax3.plot([],[], 'r-')\n",
    "    ax3.axhline(30.0, color='blue', ls='--', linewidth=0.5)\n",
    "\n",
    "    ax3.set_title(\"Angular momentum as a function of time\", fontsize = 10)\n",
    "    ax3.set_xlabel(\"Time\")\n",
    "    ax3.set_xlim(0, 20.0)\n",
    "    ax3.set_ylim(20, 40)\n",
    "\n",
    "\n",
    "xt = [];\n",
    "yt = [];\n",
    "\n",
    "lim1 = 0;\n",
    "lim2 = 20.0;\n",
    "\n",
    "t0 = 0.0;\n",
    "te = [];\n",
    "epot = [];\n",
    "ekint = [];\n",
    "etot = [];\n",
    "am = [];\n",
    "\n",
    "def init():\n",
    "    return (earth, avector, fvector, trace)\n",
    "\n",
    "    \n",
    "def update(i):\n",
    "    \"\"\"the function to update the frame\"\"\"\n",
    "    global r1, v, a, xt, yt, epot, ekint, etot, te, t0, lim1, lim2, am\n",
    "    current_time = datetime.now()\n",
    "    delta_time=(current_time-init_time).total_seconds()\n",
    "    if (delta_time>=60.0):\n",
    "        button_pause.description = \"Play\"\n",
    "        anim.event_source.stop()\n",
    "        return \n",
    "    r1, v, a = verlet_update(r1, r2, v, m1, m2, dt, a)\n",
    "    e1, e2, e3 = compute_energies(r1[-1,:], r2, m1, m2, v)\n",
    "\n",
    "    earth.set_data([r1[-1, 0]], [r1[-1, 1]])\n",
    "    xt.append(r1[-1, 0])\n",
    "    yt.append(r1[-1, 1])\n",
    "    \n",
    "    epot.append(e1)\n",
    "    ekint.append(e2)\n",
    "    etot.append(e1+e2)\n",
    "    am.append(e3)\n",
    "    \n",
    "    t0 = t0 + dt\n",
    "    te.append(t0)\n",
    "    \n",
    "    if te[-1] > 20.0:\n",
    "        lim1 = lim1 + dt\n",
    "        lim2 = lim2 + dt\n",
    "        ax2.set_xlim(lim1, lim2)\n",
    "        ax3.set_xlim(lim1, lim2)\n",
    "\n",
    "    lpot.set_data([te, epot])\n",
    "    lkint.set_data([te, ekint])\n",
    "    ltot.set_data([te, etot])\n",
    "    lang.set_data([te, am])\n",
    "    \n",
    "    trace.set_data(xt, yt)\n",
    "    avector.set_offsets(np.array([r1[-1, 0], r1[-1, 1]]))\n",
    "    fvector.set_offsets(np.array([r1[-1, 0], r1[-1, 1]]))\n",
    "    avector.set_UVC(a[0], a[1])\n",
    "    fvector.set_UVC(v[0], v[1])\n",
    "    \n",
    "    title.set_text(\"v(t) = %-5.3f\" % norm(v))\n",
    "\n",
    "    return (earth, avector, fvector, trace, lpot, lang)\n",
    "\n",
    "# create animation for the matplotlib figure\n",
    "anim = FuncAnimation(fig, update, frames = 20000, interval = 10,\n",
    "                    init_func=init, blit=True, repeat = True)\n",
    "\n",
    "\n",
    "def onClick(event):\n",
    "    \"\"\"callback function for the play and pause button\"\"\"\n",
    "    global init_time\n",
    "    if button_pause.description == \"Pause\":\n",
    "        button_pause.description = \"Play\"\n",
    "\n",
    "        anim.event_source.stop()\n",
    "    else:\n",
    "        button_pause.description = \"Pause\"\n",
    "        init_time=datetime.now()\n",
    "\n",
    "        anim.event_source.start()\n",
    "        \n",
    "\n",
    "def clear_data(b):\n",
    "    \"\"\"clear function for the Clear trace button\"\"\"\n",
    "    global xt, yt\n",
    "    \n",
    "    xt = [];\n",
    "    yt = [];\n",
    "\n",
    "button_pause = Button(description=\"Play\")\n",
    "button_pause.on_click(onClick)\n",
    "button_clear = Button(description=\"Clear trace\")\n",
    "button_clear.on_click(clear_data)\n",
    "\n",
    "plt.show()\n",
    "\n",
    "def on_dt_change(c):\n",
    "    \"\"\"callback for the change event of the slider\"\"\"\n",
    "    global dt, r1, xt, yt\n",
    "    \n",
    "    anim.event_source.stop()\n",
    "    \n",
    "    dt = w_dt.value*0.01\n",
    "    r1 = np.array([[10, 0],[10.0, 3.0*dt]])\n",
    "    \n",
    "    xt = []\n",
    "    yt = []\n",
    "    \n",
    "    anim.event_source.start()\n",
    "\n",
    "# the slider to change the time step\n",
    "w_dt = IntSlider(value = 1, min = 1, max = 50, step=1, description=\"$\\Delta t$: \", orientation='horizontal',\n",
    "                layout = Layout(width=\"400px\"))\n",
    "w_dt.observe(on_dt_change, names='value')\n",
    "\n",
    "anim.pause()\n",
    "display(output, w_dt, HBox([button_pause, button_clear]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"height:1px;border:none;color:#cccccc;background-color:#cccccc;\" />\n",
    "\n",
    "## Legend\n",
    "(How to use the interactive visualization)\n",
    "\n",
    "### Interactive figures\n",
    "The left figure shows the sun and earth system in two dimensions. The \n",
    "sun is in the center of the figure as shown in the yellow color. The \n",
    "earth is represented as a red dot. The trace of the earth's orbit is \n",
    "shown in the plot as the black line. The directions of the velocity \n",
    "and acceleration of the earth are shown as two vectors. The lengths\n",
    "of the vectors indicate the values of the velocity and acceleration.\n",
    "\n",
    "The top right figure shows how the energies change with time. The total \n",
    "energy is a summation of kinetic energy and potential energy. The plot \n",
    "is used to monitor the energy conservation of the system.\n",
    "\n",
    "The bottom right figure shows angular momentum as a function of time.\n",
    "\n",
    "### Controls\n",
    "You can click the \"Play & Pause\" button to run or stop the simulation.\n",
    "The slider is used to choose the time step $\\Delta t$. The \"Clear trace\"\n",
    "button can be used to remove the trace of the earth's orbit."
   ]
  },
  {
   "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.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}