06-Molecular_dynamics.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.animation import FuncAnimation\n",
    "mpl.rcParams[\"animation.embed_limit\"] = 50\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "import chemiscope\n",
    "from widget_code_input import WidgetCodeInput\n",
    "from ipywidgets import Layout, Output, Text, HTML, HBox, VBox\n",
    "from scwidgets import (AnswerRegistry, TextareaAnswer, CodeDemo,\n",
    "                       ParametersBox, PyplotOutput, ClearedOutput,\n",
    "                       AnimationOutput,CheckRegistry,Answer)\n",
    "from collections.abc import Iterable\n",
    "\n",
    "import ase\n",
    "from ase.io import read, write\n",
    "from ase.calculators import lj, eam\n",
    "from ase.optimize import BFGS, LBFGS\n",
    "from tqdm.notebook import tqdm\n",
    "import math\n",
    "from IPython.display import HTML\n",
    "\n",
    "from warnings import filterwarnings\n",
    "filterwarnings(\"ignore\", message=\"chemiscope\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib widget"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(mpl.get_backend())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%%html\n",
    "<style>\n",
    ".jp-CodeCell.jp-mod-outputsScrolled .jp-Cell-outputArea  {  height:auto !important;\n",
    "    max-height: 5000px; overflow-y: hidden }\n",
    "</style>\n",
    "<style>\n",
    ".output_wrapper, .output {\n",
    "    height:auto !important;\n",
    "    max-height:4000px;  /* your desired max-height here */\n",
    "}\n",
    ".output_scroll {\n",
    "    box-shadow:none !important;\n",
    "    webkit-box-shadow:none !important;\n",
    "}\n",
    "</style>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "check_registry = CheckRegistry() \n",
    "answer_registry = AnswerRegistry(prefix=\"module_06\")\n",
    "display(answer_registry)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can write here general comments you may have on this module. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "module_summary = TextareaAnswer(\"general comments on this module\")\n",
    "answer_registry.register_answer_widget(\"module-summary\", module_summary)\n",
    "display(module_summary)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Reference textbook / figure credits: Allen, Tildesley, Computer simulations of liquids (2017): Chapter 3_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Classical mechanics: Newton, Hamilton and planetary motion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The basic idea behind molecular dynamics is to apply to atoms the same kinematic laws that are applied to macroscopic objects. This is as simple as Newton's second law: if $\\mathbf{x}$ corresponds to the position of a particle, its motion is governed by \n",
    "\n",
    "$$\n",
    "\\ddot{\\mathbf{x}} = \\mathbf{f}/m = -\\frac{1}{m}\\frac{\\partial V}{\\partial \\mathbf{x}} \n",
    "$$\n",
    "\n",
    "where we also use the definition of the force acting on each particle as the derivative of a potential energy. \n",
    "\n",
    "The equations of motion can also be written in several alternative forms. A particularly elegant one involves formulating _Hamilton's equations_\n",
    "\n",
    "$$\n",
    "\\dot{\\mathbf{x}} = \\mathbf{p}/m; \\quad \\dot{\\mathbf{p}} = \\mathbf{f}\n",
    "$$\n",
    "\n",
    "where $\\mathbf{p}$ is the _momentum_ of the particle. It's clear that the second law can be recovered by taking the time derivative of the first equation, and substituting the equation for the derivative of the momentum. \n",
    "Note also that these equations can be written in a symmetric form by defining the Hamiltonian \n",
    "\n",
    "$$\n",
    "H(\\mathbf{x},\\mathbf{p}) = V(\\mathbf{x}) + \\mathbf{p}^2/2m\n",
    "$$\n",
    "\n",
    "that corresponds to the total (potential+kinetic) energy. Then, one can write\n",
    "\n",
    "$$\n",
    "\\dot{\\mathbf{x}} = \\frac{\\partial H}{\\partial\\mathbf{p}}; \\quad \\dot{\\mathbf{p}} = -\\frac{\\partial H}{\\partial\\mathbf{x}}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span style=\"color:blue\">**01** Compute $d{H}/d{t}$ for a particle that follows Hamilton's equations. Sketch the key steps of the derivation. If at $t=0$ the Hamiltonian evaluates to $1$kJ, how large will it be at $t=10s$? </span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ex1_txt = TextareaAnswer(\"Enter your answer here\")\n",
    "answer_registry.register_answer_widget(\"ex1-answer\", ex1_txt)\n",
    "display(ex1_txt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The widget below shows the analytical solution for the motion of two planets with a given mass and initial velocities. The input is the initial position and velocity of the first body. The initial position and velocity of the other one are set such a way to make the center of the mass motionless and located at [0, 0] \n",
    "The initial parameters are set to match roughly those for the earth ($m_2$) and moon ($m_1$). Experiment with different parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# constants to be used for units conversion\n",
    "G = 6.67430e-11 # gravitational constant in SI units\n",
    "m_earth = 5.972e24 # mass of earth in kg\n",
    "d_em = 383.4e6 # earth-moon distance in m\n",
    "\n",
    "def get_theoretical_trajectories(x1_initial, x2_initial, v1_initial, v2_initial, \n",
    "                                 m1, m2, k=G, n_points = 200):\n",
    "    \"\"\"Computes the elliptical orbits for the Keplerian 2-body problem. \n",
    "    k is the universal gravitational constant. Assumes inputs in SI units\"\"\"\n",
    "        \n",
    "    delta_x_initial = x1_initial - x2_initial\n",
    "    delta_v_initial = v1_initial - v2_initial\n",
    "    mu = m1 * m2 / (m1 + m2)\n",
    "    r = np.sqrt(np.sum(delta_x_initial * delta_x_initial))\n",
    "    L = mu * np.abs(np.cross(delta_x_initial, delta_v_initial))\n",
    "    alpha = k*m1*m2\n",
    "    E = -alpha / r + 0.5 * mu * np.sum(delta_v_initial ** 2)\n",
    "    e = np.sqrt(1 + 2 * E * L * L / (alpha**2 * mu))\n",
    "    theta_initial = np.arctan2(delta_x_initial[1], delta_x_initial[0])\n",
    "    \n",
    "    if np.abs(e) > 1e-7:\n",
    "        cos = (L * L / (alpha * mu * r) - 1) / e\n",
    "        # fixes rounding errors\n",
    "        cos = min(1,max(-1,cos))\n",
    "        \n",
    "        delta_theta = np.arccos(cos)\n",
    "        first_theta_0 = theta_initial - delta_theta\n",
    "        second_theta_0 = theta_initial + delta_theta\n",
    "\n",
    "        def get_velocity_direction(theta_0):\n",
    "            r_p_theta = -L * L / (alpha * mu) * e * np.sin(theta_0 - theta_initial) / (e * np.cos(theta_0 - theta_initial) + 1) ** 2\n",
    "            vx = r_p_theta * np.cos(theta_initial) - np.sin(theta_initial)\n",
    "            vy = r_p_theta * np.sin(theta_initial) + np.cos(theta_initial)\n",
    "\n",
    "            length = np.sqrt(vx ** 2 + vy ** 2)\n",
    "            return  np.array([vx / length, vy / length])\n",
    "\n",
    "        first_direction = get_velocity_direction(first_theta_0)\n",
    "        second_direction = get_velocity_direction(second_theta_0)\n",
    "\n",
    "        if np.abs(np.dot(first_direction, delta_v_initial)) > np.abs(np.dot(second_direction, delta_v_initial)):\n",
    "            theta_0 = first_theta_0\n",
    "        else:\n",
    "            theta_0 = second_theta_0\n",
    "    else:\n",
    "        theta_0 = 0.0\n",
    "    \n",
    "    theoretical_1, theoretical_2 = [], []\n",
    "    has_break = False\n",
    "    for theta in np.linspace(0, 2 * np.pi, n_points):\n",
    "        r_now = L * L / (alpha * mu * (1 + e * np.cos(theta - theta_0)))\n",
    "        if (r_now > 0.0):\n",
    "            r_first = r_now * m2 / (m1 + m2)\n",
    "            r_second = r_now * m1 / (m1 + m2)\n",
    "\n",
    "            theoretical_1.append([r_first * np.cos(theta), r_first * np.sin(theta)])\n",
    "            theoretical_2.append([-r_second * np.cos(theta), -r_second * np.sin(theta)])\n",
    "        else:\n",
    "            theoretical_1.append(None)\n",
    "            theoretical_2.append(None)\n",
    "            has_break = True\n",
    "            \n",
    "            \n",
    "    if not has_break:\n",
    "        theoretical_1_aligned = np.array(theoretical_1)\n",
    "        theoretical_2_aligned = np.array(theoretical_2)\n",
    "    else:\n",
    "        if (theoretical_1[0] is None) and (theoretical_1[-1] is None):\n",
    "            theoretical_1_aligned, theoretical_2_aligned = [], []\n",
    "            for i in range(len(theoretical_1)):\n",
    "                if theoretical_1[i] is not None:\n",
    "                    theoretical_1_aligned.append(theoretical_1[i])\n",
    "                    theoretical_2_aligned.append(theoretical_2[i])\n",
    "            theoretical_1_aligned = np.array(theoretical_1_aligned)\n",
    "            theoretical_2_aligned = np.array(theoretical_2_aligned)\n",
    "        else:\n",
    "            theoretical_1_aligned, theoretical_2_aligned = [], []\n",
    "            for i in range(len(theoretical_1)):\n",
    "                if theoretical_1[i] is None:\n",
    "                    if (i + 1 == len(theoretical_1)) or (theoretical_1[i + 1] is not None):\n",
    "                        begining = i + 1\n",
    "                        break\n",
    "            for i in range(begining, len(theoretical_1)):\n",
    "                theoretical_1_aligned.append(theoretical_1[i])\n",
    "                theoretical_2_aligned.append(theoretical_2[i])\n",
    "            for i in range(len(theoretical_1)):\n",
    "                if theoretical_1[i] is None:\n",
    "                    break\n",
    "                theoretical_1_aligned.append(theoretical_1[i])\n",
    "                theoretical_2_aligned.append(theoretical_2[i])\n",
    "            theoretical_1_aligned = np.array(theoretical_1_aligned)\n",
    "        theoretical_2_aligned = np.array(theoretical_2_aligned)\n",
    "\n",
    "    return theoretical_1_aligned, theoretical_2_aligned, E, e"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     3
    ],
    "tags": []
   },
   "outputs": [],
   "source": [
    "def plot_theoretical(m1, m2, x1, y1, v1_x, v1_y, visualizers, k=G,\n",
    "                     n_points = 200, margin = 0.1):\n",
    "    \"\"\" Plots the theoretical gravitational 2B problem trajectory. \n",
    "    masses are expressed in earth masses, positions in units of the earth-moon distance\n",
    "    and velocities in km/s (which is roughly the mean orbital velocity of the moon)\n",
    "    \"\"\"\n",
    "    \n",
    "    ax = visualizers[1].figure.get_axes()[0]\n",
    "    # converts to SI units\n",
    "    m1 *= m_earth; m2*= m_earth\n",
    "    x1 *= d_em; y1 *= d_em\n",
    "    v1_x *= 1e3\n",
    "    v1_y *= 1e3\n",
    "    \n",
    "    x1_initial = np.array([x1, y1])\n",
    "    v1_initial = np.array([v1_x, v1_y])\n",
    "    if np.abs(np.cross(x1_initial, v1_initial)) < 1e-7:\n",
    "        ax.set_title(f\"angular momentum is zero, cannot evaluate analytical solution\")\n",
    "        return\n",
    "    \n",
    "    x2_initial = -m1 * x1_initial / m2\n",
    "    v2_initial = -m1 * v1_initial / m2\n",
    "        \n",
    "    theoretical_1, theoretical_2, E, e = get_theoretical_trajectories(x1_initial, x2_initial, v1_initial, v2_initial, \n",
    "                                 m1, m2, k, n_points = n_points)\n",
    "    \n",
    "    if (e < 1):\n",
    "        min_x = min(np.min(theoretical_1[:, 0]), np.min(theoretical_2[:, 0]))\n",
    "        max_x = max(np.max(theoretical_1[:, 0]), np.max(theoretical_2[:, 0]))\n",
    "        spread = max_x - min_x\n",
    "        min_x -= spread * margin\n",
    "        max_x += spread * margin\n",
    "\n",
    "        min_y = min(np.min(theoretical_1[:, 1]), np.min(theoretical_2[:, 1]))\n",
    "        max_y = max(np.max(theoretical_1[:, 1]), np.max(theoretical_2[:, 1]))\n",
    "        spread = max_y - min_y\n",
    "        min_y -= spread * margin\n",
    "        max_y += spread * margin\n",
    "    else:\n",
    "        first_d = np.min(np.sqrt(np.sum(theoretical_1 ** 2, axis = 1)))\n",
    "        second_d = np.min(np.sqrt(np.sum(theoretical_2 ** 2, axis = 1)))\n",
    "        first_d_initial = np.sqrt(np.sum(x1_initial ** 2))\n",
    "        second_d_initial = np.sqrt(np.sum(x2_initial ** 2))\n",
    "        \n",
    "        limit = 5 * max(first_d, second_d, first_d_initial, second_d_initial)\n",
    "        min_x, min_y = -limit, -limit\n",
    "        max_x, max_y = limit, limit\n",
    "        \n",
    "   \n",
    "    min_both = min(min_x, min_y)\n",
    "    max_both = max(max_x, max_y)\n",
    "    ax.set_xlim([min_both, max_both])\n",
    "    ax.set_ylim([min_both, max_both])\n",
    "    \n",
    "    max_spread = max(max_x - min_x, max_y - min_y)\n",
    "    \n",
    "    circle1 = plt.Circle((x1_initial[0], x1_initial[1]), 0.03 * max_spread, color='b',\n",
    "                             label = f'first; $m$ = {m1/m_earth} $m_E$')\n",
    "    ax.add_patch(circle1)\n",
    "    \n",
    "    \n",
    "    circle2 = plt.Circle((x2_initial[0], x2_initial[1]), 0.03 * max_spread, color='g', \n",
    "                             label = f'second; $m$ = {m2/m_earth} $m_E$')\n",
    "    ax.add_patch(circle2)\n",
    "    ax.plot(theoretical_1[:, 0], theoretical_1[:, 1], color = 'b')\n",
    "    ax.plot(theoretical_2[:, 0], theoretical_2[:, 1], color = 'g')\n",
    "    \n",
    "    \n",
    "    def get_center_mass():\n",
    "        return (x1_initial * m1 + x2_initial * m2) / (m1 + m2)\n",
    "    center_mass = get_center_mass()\n",
    "    \n",
    "    ax.plot([center_mass[0]], [center_mass[1]], 'x', color = 'red',\n",
    "                   label = 'center of mass')\n",
    "    \n",
    "    first_length = np.sqrt(np.sum(v1_initial ** 2))\n",
    "    second_length = np.sqrt(np.sum(v2_initial ** 2))\n",
    "    \n",
    "    max_length = max(first_length, second_length)\n",
    "   \n",
    "    multiplier = 0.3 * max_spread / max_length\n",
    "    \n",
    "    ax.arrow(x1_initial[0], x1_initial[1], v1_initial[0] * multiplier , v1_initial[1] * multiplier, fc = 'b',\n",
    "             head_width = 0.03 * max_spread)\n",
    "    ax.arrow(x2_initial[0], x2_initial[1], v2_initial[0] * multiplier , v2_initial[1] * multiplier, fc = 'g',\n",
    "             head_width = 0.03 * max_spread)\n",
    "    \n",
    "    ax.legend(loc = 'upper left')\n",
    "    if e < 1.0:\n",
    "        ax.set_title(f\"Energy = {(E/1e30):9.3f} $\\\\times 10^{{30}}$J, eccentricity = {e:9.3f}\")\n",
    "    else:\n",
    "        ax.set_title(f\"Energy = {(E/1e30):9.3f} $\\\\times 10^{{30}}$J, eccentricity = {e:9.3f}\\n trajectory is unbound\")\n",
    "                         \n",
    "    \n",
    "fig_ax_analytical, _ = plt.subplots(figsize=(6, 6))\n",
    "analytical_pyplot_output = PyplotOutput(fig_ax_analytical)\n",
    "analytical_pb = ParametersBox( m1 = (0.01, 0.01, 10, 0.01, r'$m_1$ / $m$(moon)'),\n",
    "                               m2 = (1.0, 0.01, 10, 0.01, r'$m_2$ / $m$(earth)'),\n",
    "                               x1 = (1.0, -10, 10, 0.01, r'$x_1(0)$ / $d$(earth-moon)'),\n",
    "                               y1 = (0.0, -10, 10, 0.01, r'$y_1(0)$ / $d$(earth-moon)'),\n",
    "                               v1_x = (0.0, -2, 2, 0.01, r\"$v_{1x}$ / km/s\"),\n",
    "                               v1_y = (1.0, -2, 2, 0.01, r'$v_{1y}$ / km/s'),\n",
    "                               )\n",
    "\n",
    "analytical_code_demo = CodeDemo(\n",
    "            input_parameters_box=analytical_pb,\n",
    "            visualizers = [ClearedOutput(),analytical_pyplot_output],\n",
    "            update_visualizers = plot_theoretical)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "display(analytical_code_demo)\n",
    "analytical_code_demo.run_demo()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Integrators: numerically solving the equations of motion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Even in the simple case of a central potential, analytical solutions are possible only for a two-body setup. Even just the [three-body case](https://en.wikipedia.org/wiki/Three-body_problem) (e.g. a Sun-Earth-Moon scenario) doesn't have a closed-form solution. \n",
    "\n",
    "It then becomes necessary to use numerical methods to integrate the equations of motion. The simplest method possible corresponds to a naïve discretization of the expression for the first-order derivatives in the Hamiltonian formulation:\n",
    "\n",
    "$$\n",
    "\\mathbf{x}(dt) = \\mathbf{x}(0) + dt\\  \\mathbf{p}(0)/m\\\\\n",
    "\\mathbf{p}(dt) = \\mathbf{p}(0) + dt\\ \\mathbf{f}(0)\\\\\n",
    "$$\n",
    "\n",
    "This is called _forward Euler_ method, and $dt$ is the _timestep_ used to integrate the equations of motion. The procedure is iterated many times, until the trajectory has evolved for a sufficiently long time to "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span style=\"color:blue\">**02a** Complete the function below to implement a forward Euler integrator for the equations of motion of two planets, given the masses and the initial condition. </span>\n",
    "\n",
    "_Most of the routine is already implemented, only modify the part indicated, but try to understand what the rest of the code is doing. The routine returns the time series of planet positions, potential and kinetic energy._"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def simulate(updater, m1, m2, x1_initial, x2_initial, v1_initial, v2_initial, dt, num):\n",
    "    x1, x2, v1, v2 = x1_initial, x2_initial, v1_initial, v2_initial\n",
    "    x1_log, x2_log, v1_log, v2_log = [x1], [x2], [v1], [v2]\n",
    "    for _ in tqdm(range(num), leave=True):\n",
    "        x1, x2, v1, v2 = updater(m1, m2, x1, x2, v1, v2, dt)\n",
    "        \n",
    "        x1_log.append(x1)\n",
    "        x2_log.append(x2)\n",
    "\n",
    "        v1_log.append(v1)\n",
    "        v2_log.append(v2)\n",
    "\n",
    "    return np.array(x1_log), np.array(x2_log), np.array(v1_log), np.array(v2_log)\n",
    "        \n",
    "def simulate_center_fix(updater, m1, m2, x1_initial, v1_initial, dt, T):\n",
    "    \n",
    "    x2_initial = -m1 * x1_initial / m2\n",
    "    v2_initial = -m1 * v1_initial / m2\n",
    "    num = math.ceil(T / dt)\n",
    "    return simulate(updater, m1, m2, x1_initial, x2_initial, v1_initial, v2_initial, dt, num)\n",
    "\n",
    "def make_anim_simulation(updater, base_figure, ax, m1, m2, x1, y1, v1_x, v1_y, dt, T, \n",
    "                         duration = 10, fps = 15, margin = 0.1, compress_to = 10000):\n",
    "    \n",
    "    \"\"\" Creates a matplotlib animation object visualizing a trajectory of a planetary body, \n",
    "    and returns it. \"\"\" \n",
    "    \n",
    "    m1 *= m_earth; m2*= m_earth\n",
    "    x1 *= d_em; y1 *= d_em\n",
    "    v1_x *= 1e3\n",
    "    v1_y *= 1e3\n",
    "    dt*=86400; T*=86400\n",
    "   \n",
    "    x1_initial = np.array([x1, y1])\n",
    "    v1_initial = np.array([v1_x, v1_y])\n",
    "    if np.abs(np.cross(x1_initial, v1_initial)) < 1e-7:\n",
    "        print(f\"angular momentum is zero \\n It is inappropriate initial condition for this model\")\n",
    "        return\n",
    "    \n",
    "    print(\"Performing simulation...\")\n",
    "    x1_log, x2_log, v1_log, v2_log = simulate_center_fix(updater, m1, m2, x1_initial, v1_initial, dt, T)\n",
    "        \n",
    "    print(\"Simulation finished.\\n\")\n",
    "    \n",
    "    print(\"Making animation...\")\n",
    "    \n",
    "    delta_x_initial = x1_log[0] - x2_log[0]\n",
    "    delta_v_initial = v1_log[0] - v2_log[0]\n",
    "\n",
    "    theoretical_1, theoretical_2, E, e = get_theoretical_trajectories(x1_log[0], x2_log[0], v1_log[0], v2_log[0], \n",
    "                                                                    m1, m2, G)\n",
    "\n",
    "    T = dt * x1_log.shape[0]\n",
    "    \n",
    "    total_draws = duration * fps    \n",
    "    if compress_to is not None:\n",
    "        step = math.floor(x1_log.shape[0] / compress_to)\n",
    "        if step > 1:\n",
    "            x1_log = x1_log[::step]\n",
    "            x2_log = x2_log[::step]\n",
    "            v1_log = v1_log[::step]\n",
    "            v2_log = v2_log[::step]\n",
    "            dt = dt * step\n",
    "\n",
    "    def get_energies():\n",
    "        deltas = x1_log - x2_log\n",
    "        distances = np.sqrt(np.sum(deltas ** 2, axis = 1))\n",
    "        potential_energies = -G*m1*m2 / distances\n",
    "        kinetic_energies = 0.5 * (m1 * np.sum(v1_log ** 2, axis = 1) + m2 * np.sum(v2_log ** 2, axis = 1))\n",
    "        return potential_energies, kinetic_energies, potential_energies + kinetic_energies\n",
    "\n",
    "    def get_center_mass():\n",
    "        return (x1_log * m1 + x2_log * m2) / (m1 + m2)\n",
    "\n",
    "    potential_energies, kinetic_energies, energies = get_energies()\n",
    "    center_mass = get_center_mass()\n",
    "    \n",
    "    if e < 1:\n",
    "        min_x = min(np.min(x1_log[:, 0]), np.min(x2_log[:, 0]), np.min(theoretical_1[:, 0]), np.min(theoretical_2[:, 0]))\n",
    "        max_x = max(np.max(x1_log[:, 0]), np.max(x2_log[:, 0]), np.max(theoretical_1[:, 0]), np.max(theoretical_2[:, 0]))\n",
    "    else:\n",
    "        min_x = min(np.min(x1_log[:, 0]), np.min(x2_log[:, 0]))\n",
    "        max_x = max(np.max(x1_log[:, 0]), np.max(x2_log[:, 0]))\n",
    "    \n",
    "    spread = max_x - min_x\n",
    "    min_x -= spread * margin\n",
    "    max_x += spread * margin\n",
    "    \n",
    "    if e < 1:\n",
    "        min_y = min(np.min(x1_log[:, 1]), np.min(x2_log[:, 1]), np.min(theoretical_1[:, 1]), np.min(theoretical_2[:, 1]))\n",
    "        max_y = max(np.max(x1_log[:, 1]), np.max(x2_log[:, 1]), np.max(theoretical_1[:, 1]), np.max(theoretical_2[:, 1]))\n",
    "    else:\n",
    "        min_y = min(np.min(x1_log[:, 1]), np.min(x2_log[:, 1]))\n",
    "        max_y = max(np.max(x1_log[:, 1]), np.max(x2_log[:, 1]))\n",
    "        \n",
    "    spread = max_y - min_y\n",
    "    min_y -= spread * margin\n",
    "    max_y += spread * margin\n",
    "\n",
    "    min_both = min(min_x, min_y)\n",
    "    max_both = max(max_x, max_y)\n",
    "    max_spread = max(max_x - min_x, max_y - min_y)\n",
    "    \n",
    "    min_e = min(np.min(potential_energies), np.min(kinetic_energies), np.min(energies))\n",
    "    max_e = max(np.max(potential_energies), np.max(kinetic_energies), np.max(energies))\n",
    "    delta = max_e - min_e\n",
    "    min_e -= delta * margin\n",
    "    max_e += delta * margin\n",
    "\n",
    "    def init():\n",
    "\n",
    "        ax[0].set_xlabel('$x$/m', fontsize = 12)\n",
    "        ax[0].set_ylabel('$y$/m', fontsize = 12)\n",
    "        ax[0].set_xlim([min_both, max_both])\n",
    "        ax[0].set_ylim([min_both, max_both])\n",
    "\n",
    "\n",
    "        line_center, = ax[0].plot([], 'x', color = 'red',\n",
    "                       label = 'center of mass')\n",
    "\n",
    "        line_theoretical_1, = ax[0].plot([], color = 'b', label = 'analytical traj. 1')\n",
    "        line_theoretical_2, = ax[0].plot([], color = 'g', label = 'analytical traj. 2')\n",
    "\n",
    "        line_1,  = ax[0].plot([], color = 'gray', alpha = 0.5)\n",
    "        line_2,  = ax[0].plot([], color = 'gray', alpha = 0.5)\n",
    "\n",
    "        circle_1 = plt.Circle((x1_log[0, 0], x1_log[0, 1]), 0.03 * max_spread, color='b',\n",
    "                                 label = f'first; $m$ = {m1/m_earth} $m_E$')\n",
    "        ax[0].add_patch(circle_1)\n",
    "\n",
    "        circle_2 = plt.Circle((x2_log[0, 0], x2_log[0, 1]), 0.03 * max_spread, color='g', \n",
    "                                 label = f'second; $m$ = {m2/m_earth} $m_E$')\n",
    "        ax[0].add_patch(circle_2)\n",
    "\n",
    "        ax[0].legend(loc = 'upper left')\n",
    "        if e < 1.0:\n",
    "            ax[0].set_title(f\"Energy = {(E/1e30):9.3f} $\\\\times 10^{{30}}$J, eccentricity = {e:9.3f}\")\n",
    "        else:\n",
    "            ax[0].set_title(f\"Energy = {(E/1e30):9.3f} $\\\\times 10^{{30}}$J, eccentricity = {e:9.3f}\\n trajectory is unbound\")\n",
    "\n",
    "\n",
    "        line_pot, = ax[1].plot([], label = 'potential energy')\n",
    "        line_kin, = ax[1].plot([], label = 'kinetic energy')\n",
    "        line_tot, = ax[1].plot([], label = 'total energy')\n",
    "\n",
    "\n",
    "        ax[1].set_xlim([0, T/86400])\n",
    "        ax[1].set_ylim([min_e, max_e])\n",
    "        ax[1].legend(loc = 'upper right')\n",
    "        ax[1].set_xlabel(\"Time / days\", fontsize = 12)\n",
    "        ax[1].set_ylabel(\"Energy / J\", fontsize = 12)\n",
    "        return base_figure, ax, line_pot, line_kin, line_tot, line_theoretical_1,\\\n",
    "            line_theoretical_2, line_center, circle_1, circle_2, line_1, line_2\n",
    "\n",
    "    base_figure, ax, line_pot, line_kin, line_tot, line_theoretical_1, line_theoretical_2, line_center, circle_1, circle_2, line_1, line_2 = init()\n",
    "\n",
    "    def animate(frame_num):\n",
    "        index = int(x1_log.shape[0] * frame_num / total_draws)\n",
    "        if index > x1_log.shape[0] - 1:\n",
    "            index = x1_log.shape[0] - 1\n",
    "        line_pot.set_data((np.arange(index) * dt/86400, potential_energies[0:index]))\n",
    "        line_kin.set_data((np.arange(index) * dt/86400, kinetic_energies[0:index]))\n",
    "        line_tot.set_data((np.arange(index) * dt/86400, energies[0:index]))\n",
    "\n",
    "        line_theoretical_1.set_data((theoretical_1[:, 0], theoretical_1[:, 1]))\n",
    "        line_theoretical_2.set_data((theoretical_2[:, 0], theoretical_2[:, 1]))\n",
    "        line_center.set_data((center_mass[index, 0], center_mass[index, 1]))\n",
    "        circle_1.center = x1_log[index, 0], x1_log[index, 1]\n",
    "        circle_2.center = x2_log[index, 0], x2_log[index, 1]\n",
    "        line_1.set_data((x1_log[0:index, 0], x1_log[0:index, 1]))\n",
    "        line_2.set_data((x2_log[0:index, 0], x2_log[0:index, 1]))\n",
    "        return line_pot, line_kin, line_tot, line_theoretical_1, line_theoretical_2, line_center, circle_1, circle_2, line_1, line_2\n",
    "\n",
    "    anim = FuncAnimation(base_figure, animate, frames=tqdm(range(total_draws), leave=True), \n",
    "                         interval=duration, blit = True)\n",
    "    \n",
    "    anim =  HTML(anim.to_jshtml(fps=fps))\n",
    "    print(\"Animation finished.\")\n",
    "    \n",
    "    return anim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ex02_wci = WidgetCodeInput(\n",
    "        function_name=\"euler_update\", \n",
    "        function_parameters=\"m1, m2, r1, r2, v1, v2, dt\",\n",
    "        docstring=\"\"\"Update the position of the particles in a two-body problem by\n",
    "        the forward Euler method described above.\n",
    "        Takes as arguments the masses of the two bodies (m1,m2), the initial positions\n",
    "        and velocities (r1,r2,v1,v2) and the time step (dt). Returns the new positions\n",
    "        and velocities (r1_new,r2_new,v1_new,v2_new). \n",
    "        >> All positions and velocities are numpy arrays. <<\n",
    "\"\"\",\n",
    "        function_body=\"\"\"\n",
    "\n",
    "import numpy as np\n",
    "return r1, r2, v1, v2 # <-- remove this line after having finished the implementation of the function\n",
    "G = 6.67430e-11 # gravitational constant in SI units\n",
    "# you can compute forces by calling this function with the positions of the two bodies\n",
    "def get_forces(r1, r2):\n",
    "    delta = r2 - r1\n",
    "    norm = np.sqrt(np.sum(delta ** 2))\n",
    "    direction = delta / norm\n",
    "\n",
    "    f1 = direction * G*m1*m2 / (norm ** 2)\n",
    "    f2 = -f1\n",
    "    return f1, f2\n",
    "\n",
    "# Modify this to implement a forward Euler integrator. \n",
    "r1_new = r1\n",
    "r2_new = r2\n",
    "\n",
    "v1_new = v1\n",
    "v2_new = v2\n",
    "\n",
    "return r1_new, r2_new, v1_new, v2_new\n",
    "\"\"\"\n",
    "        )\n",
    "\n",
    "def convert_celestial_input(m1, m2, r1, r2, v1, v2, dt):\n",
    "    \n",
    "    m1 *= m_earth; m2*= m_earth\n",
    "    r1 *= d_em; r2 *= d_em\n",
    "    v1 *= 1e3; v2 *= 1e3\n",
    "    dt*=86400\n",
    "    \n",
    "    return (m1, m2, r1, r2, v1, v2, dt)\n",
    "    \n",
    "first_input = (0.01, 1.0, np.array([1.0, 0.0]), np.array([-1.0, 0.0]), np.array([0.0, 1.0]), \n",
    "               np.array([0.0, -2.0]), 0.1)\n",
    "second_input = (0.01, 1.0, np.array([1.0, 0.3]), np.array([-1.0, -0.2]), np.array([0.0, 1.0]), \n",
    "               np.array([0.0, -2.0]), 0.1)\n",
    "third_input = (0.5, 1.0, np.array([1.0, 0.3]), np.array([-1.0, -0.2]), np.array([0.0, 1.0]), \n",
    "               np.array([0.0, -2.0]), 0.1)\n",
    "\n",
    "inputs = [convert_celestial_input(*inp) for inp in [first_input, second_input, third_input]]\n",
    "ref_input = [{\"m1\" : value[0], \"m2\" : value[1], \"r1\" : value[2], \"r2\" : value[3], \"v1\" : value[4], \"v2\" : value[5], \"dt\" : value[6]} for value in inputs]\n",
    "\n",
    "ref_output = [[np.array([3.834e+08, 8.640e+06]), np.array([-3.834e+08, -1.728e+07]), np.array([  -5.85699774, 1000.        ]), np.array([ 5.85699774e-02, -2.00000000e+03])],\n",
    "[np.array([3.8340e+08, 1.2366e+08]), np.array([-3.834e+08, -9.396e+07]), np.array([ -5.34787994, 998.66303002]), np.array([ 5.34787994e-02, -1.99998663e+03])],\n",
    "[np.array([3.8340e+08, 1.2366e+08]), np.array([-3.834e+08, -9.396e+07]), np.array([ -5.34787994, 998.66303002]), np.array([    2.67393997, -1999.33151501])]]\n",
    "\n",
    "ex02_pb = ParametersBox(    \n",
    "                            m1 = (0.01, 0.01, 10, 0.01, r'$m_1$ / $m$(moon)'),\n",
    "                            m2 = (1.0, 0.01, 10, 0.01, r'$m_2$ / $m$(earth)'),\n",
    "                            x1 = (1.0, -10, 10, 0.01, r'$x_1(0)$ / $d$(earth-moon)'),\n",
    "                            y1 = (0.0, -10, 10, 0.01, r'$y_1(0)$ / $d$(earth-moon)'),\n",
    "                            v1_x = (0.0, -2, 2, 0.01, r\"$v_{1x}$ / km/s\"),\n",
    "                            v1_y = (1.0, -2, 2, 0.01, r'$v_{1y}$ / km/s'),                               \n",
    "                            dt = (0.01, 0.001, 1, 0.001, r'$dt$ / days', dict(readout_format='.3f')),\n",
    "                            T = (100.0, 2, 1000, 1, r'total simulation / days'),\n",
    "                            refresh_mode =\"click\")\n",
    "with plt.ioff():\n",
    "    fig_ax_ex02,_ = plt.subplots(1, 2, figsize=(10, 5))\n",
    "\n",
    "ex02_pyplot_output = AnimationOutput(fig_ax_ex02)\n",
    "    \n",
    "anim_duration = 10\n",
    "anim_fps = 10\n",
    "def make_anim_euler_simulation(m1, m2, x1, y1, v1_x, v1_y, dt, T, code_input, visualizers):\n",
    "    \n",
    "    print_output = visualizers[0]\n",
    "    print_output.clear_output()\n",
    "    \n",
    "    figure = visualizers[1].figure\n",
    "    axes = figure.get_axes()\n",
    "    updater = code_input.get_function_object()\n",
    "    with print_output:\n",
    "        animation = make_anim_simulation(updater, figure, axes, m1, m2, x1, y1, v1_x, v1_y, dt, T,\n",
    "                                        duration=anim_duration, fps=anim_fps)    \n",
    "    # clears progress bars\n",
    "    print_output.clear_output()        \n",
    "    visualizers[1].animation = animation\n",
    "        \n",
    "ex02_code_demo = CodeDemo(\n",
    "            input_parameters_box=ex02_pb,\n",
    "            code_input= ex02_wci,\n",
    "            check_registry=check_registry,\n",
    "            visualizers = [ClearedOutput(),ex02_pyplot_output],\n",
    "            update_visualizers = make_anim_euler_simulation,\n",
    "            merge_check_and_update_buttons = False)\n",
    "\n",
    "check_registry.add_check(ex02_code_demo,\n",
    "                         inputs_parameters=ref_input,\n",
    "                         reference_outputs = ref_output,\n",
    "                         equal= np.allclose)\n",
    "\n",
    "answer_registry.register_answer_widget(\"ex02-function\", ex02_code_demo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "display(ex02_code_demo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# hack to avoid a mysterious figure disappearence bug\n",
    "ex02_code_demo.visualizers[1].clear_output()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The widget above allows you to run the integration with different starting conditions and integration time step. Experiment with a trajectory initialized with earth-moon settings (set sliders to $(m_1,m_2,x_1,y_1,v_{1x},v_{1y})=$ (0.01,1.0,1.0,0.0,0.0,1.0)) and varying the integration time step.\n",
    "\n",
    "<span style=\"color:blue\">**02b** Is the total energy conserved? What is the largest time step you can use before the trajectory becomes completely unstable? </span>\n",
    "\n",
    "_NB: you will see that for this problem the accuracy of the integrator degrades progressively as $dt$ is increased, so there isn't a clear-cut transition from a \"stable\" to \"unstable\" regime.  To give a consistent answer, look roughly for the value that leads to a doubling of the orbit radius after one revolution._"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ex2b_txt = TextareaAnswer(\"Enter your answer here\")\n",
    "answer_registry.register_answer_widget(\"ex2b-answer\", ex2b_txt)\n",
    "display(ex2b_txt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span style=\"color:blue\">**02c** Try with a more eccentric trajectory (set sliders to $(m_1,m_2,x_1,y_1,v_{1x},v_{1y})=$ (0.01,1.0,1.0,0.0,-0.5,1.0). At which point of the orbit is the violation of energy conservation more severe? Is the stability limit the same as for the earth-moon settings? </span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ex2c_txt = TextareaAnswer(\"Enter your answer here\")\n",
    "answer_registry.register_answer_widget(\"ex2c-answer\", ex2c_txt)\n",
    "display(ex2c_txt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The forward Euler integrator is not very performant, both in terms of theoretical accuracy, nor in terms of practical stability. Several [high-order integrators](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods) exist for generic differential equations. When integrating Hamilton's equations, there are considerations other than the asymptotic accuracy of the approximate integration. For example, Newtonian dynamics is time reversible, and [symplectic](https://en.wikipedia.org/wiki/Symplectic_integrator) (a property related to how a swarm of approximate trajectories started from a given volume of $\\mathbf{x},\\mathbf{p}$ space evolve in time). \n",
    "\n",
    "A simple integrator that preserves these two properties is the _velocity Verlet_ integrator,\n",
    "\n",
    "$$\n",
    "\\mathbf{p}(dt/2) = \\mathbf{p}(0) + \\frac{dt}{2} \\mathbf{f}(0)\\\\\n",
    "\\mathbf{x}(dt) = \\mathbf{x}(0) + dt \\mathbf{p}(dt/2)/m \\\\\n",
    "\\mathbf{p}(dt) = \\mathbf{p}(dt/2) + \\frac{dt}{2} \\mathbf{f}(dt) \\\\\n",
    "$$\n",
    "\n",
    "In practice, the momentum is evolved in two steps, first using the force evaluated at the starting position, and then using the force evaluated at the final position. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span style=\"color:blue\">**03a** Complete the function below to implement a velocity Verlet integrator for the equations of motion of two planets. </span>\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ex03_wci = WidgetCodeInput(\n",
    "        function_name=\"verlet_update\", \n",
    "        function_parameters=\"m1, m2, r1, r2, v1, v2, dt\",\n",
    "        docstring=\"\"\"Update the position of the particles in a two-body problem using\n",
    "        the Verlet integrator described above.\n",
    "        Takes as arguments the masses of the two bodies (m1,m2), the initial positions\n",
    "        and velocities (r1,r2,v1,v2) and the time step (dt). Returns the new positions\n",
    "        and velocities (r1_new,r2_new,v1_new,v2_new). \n",
    "        >> All positions and velocities are numpy arrays. <<\n",
    "\"\"\",\n",
    "        function_body=\"\"\"\n",
    "\n",
    "import numpy as np\n",
    "return r1, r2, v1, v2 # <-- remove this line after having finished the implementation of the function\n",
    "G = 6.67430e-11 # gravitational constant in SI units\n",
    "# you can compute forces by calling this function with the positions of the two bodies\n",
    "def get_forces(r1, r2):\n",
    "    ...\n",
    "    return f1, f2\n",
    "\n",
    "# Modify this to implement a velocity Verlet integrator. \n",
    "f1, f2 = get_forces(r1, r2)\n",
    "\n",
    "v1_mid = ...\n",
    "v2_mid = ...\n",
    "\n",
    "\n",
    "r1_new = ...\n",
    "r2_new = ...\n",
    "\n",
    "f1, f2 = get_forces(r1_new, r2_new)\n",
    "\n",
    "v1_new = ...\n",
    "v2_new = ...\n",
    "\n",
    "return r1_new, r2_new, v1_new, v2_new\n",
    "\"\"\"\n",
    "        )\n",
    "\n",
    "\n",
    "def reference_verlet_update(m1, m2, r1, r2, v1, v2, dt):\n",
    "    import numpy as np\n",
    "\n",
    "    def get_forces(r1, r2):\n",
    "        delta = r2 - r1\n",
    "        norm = np.sqrt(np.sum(delta ** 2))\n",
    "        direction = delta / norm\n",
    "\n",
    "        f1 = direction * G*m1*m2 / (norm ** 2)\n",
    "        f2 = -f1\n",
    "        return f1, f2\n",
    "\n",
    "    # edit these before updating the reference values...\n",
    "    f1, f2 = 0,0# \n",
    "\n",
    "    v1_mid = 0\n",
    "    v2_mid = 0\n",
    "\n",
    "\n",
    "    r1_new = 0\n",
    "    r2_new = 0\n",
    "\n",
    "    f1, f2 = 0,0 #\n",
    "\n",
    "    v1_new = 0\n",
    "    v2_new = 0\n",
    "\n",
    "    return r1_new, r2_new, v1_new, v2_new\n",
    "    \n",
    "first_input = (0.01, 1.0, np.array([1.0, 0.0]), np.array([-1.0, 0.0]), np.array([0.0, 1.0]), \n",
    "               np.array([0.0, -2.0]), 0.1)\n",
    "second_input = (0.01, 1.0, np.array([1.0, 0.3]), np.array([-1.0, -0.2]), np.array([0.0, 1.0]), \n",
    "               np.array([0.0, -2.0]), 0.1)\n",
    "third_input = (0.5, 1.0, np.array([1.0, 0.3]), np.array([-1.0, -0.2]), np.array([0.0, 1.0]), \n",
    "               np.array([0.0, -2.0]), 0.1)\n",
    "\n",
    "inputs = [convert_celestial_input(*inp) for inp in [first_input, second_input, third_input]]\n",
    "#ref_values = [[inp, reference_verlet_update(*inp)] for inp in inputs]\n",
    "ref_values = [[np.array([3.83374698e+08, 8.64000000e+06]), np.array([-3.83399747e+08, -1.72800000e+07]), np.array([ -5.85218015, 999.90116804]), np.array([ 5.85218015e-02, -1.99999901e+03])],\n",
    "              [np.array([3.83376897e+08, 1.23654224e+08]), np.array([-3.83399769e+08, -9.39599422e+07]), np.array([ -5.2813386 , 998.59152525]), np.array([ 5.28133860e-02, -1.99998592e+03])], \n",
    "              [np.array([3.83376897e+08, 1.23654224e+08]), np.array([-3.83388449e+08, -9.39571121e+07]), np.array([ -5.28141457, 998.59150239]), np.array([    2.64070728, -1999.29575119])]]\n",
    "                   \n",
    "ex03_pb = ParametersBox(    m1 = (0.01, 0.01, 10, 0.01, r'$m_1$ / $m$(moon)'),\n",
    "                            m2 = (1.0, 0.01, 10, 0.01, r'$m_2$ / $m$(earth)'),\n",
    "                            x1 = (1.0, -10, 10, 0.01, r'$x_1(0)$ / $d$(earth-moon)'),\n",
    "                            y1 = (0.0, -10, 10, 0.01, r'$y_1(0)$ / $d$(earth-moon)'),\n",
    "                            v1_x = (0.0, -2, 2, 0.01, r\"$v_{1x}$ / km/s\"),\n",
    "                            v1_y = (1.0, -2, 2, 0.01, r'$v_{1y}$ / km/s'),                                                           \n",
    "                            dt = (0.01, 0.001, 1, 0.001, r'$dt$ / days', dict(readout_format='.3f')),\n",
    "                            T = (100.0, 2, 1000, 1, r'total simulation / days'),\n",
    "                            refresh_mode =\"click\")\n",
    "\n",
    "fig_ax_ex03,_ = plt.subplots(1, 2, figsize=(10, 5))\n",
    "ex03_pyplot_output = AnimationOutput(fig_ax_ex03)\n",
    "\n",
    "def make_anim_verlet_simulation(m1, m2, x1, y1, v1_x, v1_y, dt, T, code_input, visualizers):\n",
    "    print_output = visualizers[0]\n",
    "    figure = visualizers[1].figure\n",
    "    axes = figure.get_axes()\n",
    "    print_output.clear_output()\n",
    "    updater = code_input.get_function_object()\n",
    "    with print_output:\n",
    "        animation = make_anim_simulation(updater, figure, axes, \n",
    "                                        m1, m2, x1, y1, v1_x, v1_y, dt, T,\n",
    "                                        duration=anim_duration, fps=anim_fps)    \n",
    "    print_output.clear_output()\n",
    "    visualizers[1].animation = animation \n",
    "        \n",
    "ex03_code_demo = CodeDemo(\n",
    "            input_parameters_box=ex03_pb,\n",
    "            code_input= ex03_wci,\n",
    "            check_registry=check_registry,\n",
    "            visualizers = [ClearedOutput(),ex03_pyplot_output],\n",
    "            update_visualizers = make_anim_verlet_simulation,\n",
    "            merge_check_and_update_buttons = False)\n",
    "\n",
    "check_registry.add_check(ex03_code_demo,\n",
    "                         inputs_parameters=ref_input,\n",
    "                         reference_outputs = ref_values,\n",
    "                         equal= np.allclose)\n",
    "\n",
    "answer_registry.register_answer_widget(\"ex03-function\", ex03_code_demo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "display(ex03_code_demo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# needed to avoid a mysterious output disappearance bug\n",
    "ex03_code_demo.visualizers[1].clear_output()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span style=\"color:blue\">**03b** Experiment with a trajectory initialized at earth-moon conditions and varying the integration time step. Try also the more eccentric trajectory (setting $v_{1x}=-0.5$ km/s). Find the stability limit and compare with the observations made for the forward Euler integrator. </span>\n",
    "\n",
    "_NB: as for the previous exercise, use the time step that leads to a doubling of the orbit radius after one period as the instability limit._"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ex3b_txt = TextareaAnswer(\"Enter your answer here\")\n",
    "answer_registry.register_answer_widget(\"ex3b-answer\", ex3b_txt)\n",
    "display(ex3b_txt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Classical molecular dynamics (MD) for molecules and crystals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This far we have looked at the use of integrators to predict the trajectories of celestial bodies. As long as we treat nuclei as classical particles (which works decently unless one considers cryogenic temperatures, or light atoms [such as hydrogen](http://arxiv.org/abs/1803.00600)) the motion of atoms follows the same Hamiltonian equations, driven by the interatomic potentials we have learned about in the [dedicated module](./04-Potentials.ipynb). \n",
    "\n",
    "Then, we use integrators implemented in `ase.md` to run some short trajectories. The usage is pretty simple\n",
    "\n",
    "```python\n",
    "from ase.md.velocitydistribution import MaxwellBoltzmannDistribution, Stationary\n",
    "from ase.md.verlet import VelocityVerlet\n",
    "from ase.calculators import eam\n",
    "from ase import units\n",
    "\n",
    "structure = ase.Atoms(...)\n",
    "structure.calc = eam.EAM(...)\n",
    "\n",
    "# set the initial velocities of the atoms from a Maxwell-Bolztmann distribution at 300K\n",
    "MaxwellBoltzmannDistribution(structure, temperature_K=300)\n",
    "# removes the center-of-mass velocity \n",
    "Stationary(structure)\n",
    "\n",
    "\n",
    "# initialize the VV integrator\n",
    "vv_integrator = VelocityVerlet(structure, 2 * units.fs)  # time step of 2fs\n",
    "\n",
    "# run for a given number of steps\n",
    "vv_integrator.run(1000)\n",
    "```\n",
    "\n",
    "You may note that we call `MaxwellBoltzmannDistribution` to set the initial velocities of the atoms. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The problem with this strategy is that the simulation will run to completion, and we are only able of checking what is happening after it is `VelocityVerlet.run` returns (it will modify `structure` in place). In order to follow what's going on, we need to store copies of the structure and/or its properties after short trajectory segments. This is called an _MD loop_ and allows to analyze the outcome of the simulation by _post processing_.\n",
    "\n",
    "```python\n",
    "trajectory = []\n",
    "for i in range(100):\n",
    "    vv_integrator.run(10)\n",
    "    print(\"V=%10.5f  K=%10.5f\" %(structure.get_potential_energy(), structure.get_kinetic_energy) )    \n",
    "    trajectory.append(structure.copy()) # <-- must make a copy otherwise the reference gets overwritten\n",
    "```\n",
    "\n",
    "_NB: the `\"format_string\" % (tuple_of_values)` syntax allows to format nicely the output. See [here](https://docs.python.org/3/library/stdtypes.html#old-string-formatting) for reference, or [here](https://docs.python.org/3/tutorial/inputoutput.html) for a more modern approach_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ex04_wci = WidgetCodeInput(\n",
    "        function_name=\"LJ_dimer_MD\", \n",
    "        function_parameters=\"time_step\",\n",
    "        docstring=\"\"\"\n",
    "Performs Molecular Dynamics for the dimer described by Lennard-Jones potential. \n",
    ":param time_step: time step for molecular dynamic in fs. \n",
    "\"\"\",\n",
    "        function_body=\"\"\"\n",
    "import numpy as np\n",
    "from ase.md.velocitydistribution import MaxwellBoltzmannDistribution, Stationary\n",
    "from ase.md.verlet import VelocityVerlet\n",
    "from ase.calculators import lj\n",
    "from ase import units\n",
    "import ase.io\n",
    "import math\n",
    "from tqdm.notebook import tqdm\n",
    "\n",
    "# LJ parameters fitted to yield roughly equilibrium distance and binding energy of H2\n",
    "e0 = 4.74 # in eV\n",
    "d0 = 0.74 # in Å\n",
    "\n",
    "structure = ase.Atoms(\"H2\", positions=[[0,-d0/2,0], [0,d0/2,0]])\n",
    "structure.calc = lj.LennardJones(sigma=d0/(2**(1/6)), epsilon=e0, rc = 4 * d0)\n",
    "\n",
    "np.random.seed(123456)\n",
    "MaxwellBoltzmannDistribution(structure, temperature_K=2000)\n",
    "Stationary(structure)\n",
    "vv_integrator = VelocityVerlet(structure, time_step * units.fs)  \n",
    "\n",
    "filename = 'module_06-ar_trajectory.extxyz'\n",
    "trajectory = []\n",
    "structure.info['time'] = 0\n",
    "structure.info['potential_energy'] = structure.get_potential_energy()\n",
    "structure.info['kinetic_energy'] = structure.get_kinetic_energy()\n",
    "structure.info['dimer_distance'] = structure.get_distance(0,1)\n",
    "trajectory.append(structure.copy())\n",
    "# run segments corresponding to 0.05 fs, and a total of 5fs\n",
    "nstep = math.ceil(0.05 / time_step)\n",
    "pbar = tqdm(range(int(5/(nstep*time_step))))\n",
    "for i in pbar:\n",
    "    vv_integrator.run(nstep)\n",
    "    pbar.set_description(f\"V={structure.get_potential_energy():10.5f} K={structure.get_kinetic_energy():10.5f}\")\n",
    "    structure.info['time'] = i*time_step*nstep    \n",
    "    structure.info['potential_energy'] = structure.get_potential_energy()\n",
    "    structure.info['kinetic_energy'] = structure.get_kinetic_energy()\n",
    "    structure.info['dimer_distance'] = structure.get_distance(0,1)\n",
    "    trajectory.append(structure.copy())\n",
    "\n",
    "ase.io.write(filename, trajectory)\n",
    "return filename\n",
    "\"\"\"\n",
    "        )\n",
    "\n",
    "def visualise_md_lj(time_step, code_input, visualizers):    \n",
    "    print_output = visualizers[0]\n",
    "    with print_output:\n",
    "        filename = code_input.get_function_object()(time_step)\n",
    "    if filename is None:\n",
    "        return\n",
    "    \n",
    "    frames = ase.io.read(filename, ':')\n",
    "    times = [structure.info['time'] for structure in frames]\n",
    "    potential_energies = [structure.info['potential_energy'] for structure in frames]\n",
    "    kinetic_energies = [structure.info['kinetic_energy'] for structure in frames]\n",
    "    dimer_distances = [structure.info['dimer_distance'] for structure in frames]\n",
    "    for f in frames:\n",
    "        f.info.clear()\n",
    "        f.arrays.pop('momenta')\n",
    "    \n",
    "    properties=dict(\n",
    "                time=times,\n",
    "                potential_energy=potential_energies,\n",
    "                kinetic_energy = kinetic_energies,\n",
    "                total_energy = np.asarray(potential_energies)+np.asarray(kinetic_energies),\n",
    "                dimer_distance = dimer_distances\n",
    "                )\n",
    "\n",
    "    chemiscope.write_input('module-06_H2_trajectory.chemiscope.json.gz', frames=frames, properties=properties)\n",
    "    mycs = chemiscope.show(frames, properties,\n",
    "                         settings={\"structure\":[{\"bonds\":False, \"unitCell\":False,\"supercell\":{\"0\":1,\"1\":1,\"2\":1},\n",
    "                                                \"keepOrientation\": True, 'playbackDelay': 100}],\n",
    "                                   'map': {\n",
    "                                       'x' : {'property':'time' },\n",
    "                                       'y' : {'property':'total_energy' },\n",
    "                                       'color' : {'property':'dimer_distance' }\n",
    "                                   } }\n",
    "                                   )\n",
    "    with print_output:\n",
    "         display( mycs )\n",
    "\n",
    "            \n",
    "ex04_pb = ParametersBox(time_step=(0.05,0.01,0.5,0.01, r\"time step, fs\"),refresh_mode=\"click\")\n",
    "ex04_code_demo = CodeDemo(\n",
    "            code_input= ex04_wci,\n",
    "            input_parameters_box=ex04_pb,\n",
    "            visualizers = [ClearedOutput()],\n",
    "            update_visualizers = visualise_md_lj,\n",
    "            merge_check_and_update_buttons = True)\n",
    "\n",
    "answer_registry.register_answer_widget(\"ex04-function\", ex04_code_demo)\n",
    "display(ex04_code_demo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Download chemiscope datafile](./module-06_H2_trajectory.chemiscope.json.gz)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span style=\"color:blue\">**04a** The widget above integrates MD for an H₂ dimer. Look at the trajectory, plotting potential, kinetic and total energy. How do the three quantities change over time? Comment in particular on the scale of the fluctuations in kinetic and potential energy in comparison with the fluctuations in total energy. Is the behavior during compression and during stretching the same? Why? </span>\n",
    "\n",
    "_Hint: to answer the last question, think of the shape of the LJ potential, and try to change the initialization temperature to lower and higher values._"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ex4a_txt = TextareaAnswer(\"Enter your answer here\")\n",
    "answer_registry.register_answer_widget(\"ex4a-answer\", ex4a_txt)\n",
    "display(ex4a_txt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span style=\"color:blue\">**04b** Increase progressively the time step and re-run the simulation. What happens to the different energies? What is the largest time step you can use before the simulation becomes unstable?  </span>\n",
    "\n",
    "_NB: similar to the planetary case, this is not a clear-cut transition, but ultimtely the molecule will break down. Give a rough estimate (+- 0.1fs) of the point where the error in the conservation of the total energy is comparable to the fluctuations in potential and kinetic energy_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ex4b_txt = TextareaAnswer(\"Enter your answer here\")\n",
    "answer_registry.register_answer_widget(\"ex4b-answer\", ex4b_txt)\n",
    "display(ex4b_txt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span style=\"color:blue\">**04c** Change the function to generate an He dimer, without changing the LJ parameters. In practice, you are just changing the masses to 4 atomic mass units. What do you observe in terms of the vibrational frequencies, and of the limiting time step for the dynamics? </span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ex4c_txt = TextareaAnswer(\"Enter your answer here\")\n",
    "answer_registry.register_answer_widget(\"ex4c-answer\", ex4c_txt)\n",
    "display(ex4c_txt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span style=\"color:blue\">**05a** Now, a slightly more serious exercise. Write a function that generates a 2x2x2 supercell of _fcc_ aluminum, instrument the structure with an EAM calculator, initialize the velocities at 400K and runs 1ps of molecular dynamics. The time step in fs is given as a parameter.  </span>\n",
    "\n",
    "_NB: most of the \"logging\" code is already implemented_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ex05_wci = WidgetCodeInput(\n",
    "        function_name=\"aluminum_MD\", \n",
    "        function_parameters=\"time_step\",\n",
    "        docstring=\"\"\"\n",
    "Performs Molecular Dynamics for the 2x2x2 supercell of fcc aluminum described by EAM potential. \n",
    ":param time_step: time step for molecular dynamic in fs. \n",
    "\"\"\",\n",
    "        function_body=\"\"\"\n",
    "from ase.md.velocitydistribution import MaxwellBoltzmannDistribution, Stationary\n",
    "from ase.md.verlet import VelocityVerlet\n",
    "from ase.calculators import eam\n",
    "from ase import units\n",
    "import ase.io\n",
    "import math\n",
    "from tqdm.notebook import tqdm\n",
    "filename = 'aluminum_trajectory.extxyz'\n",
    "return filename # <-- remove this after having completed the code below\n",
    "structure = ...\n",
    "\n",
    "structure.calc = eam.EAM(potential='data/Al99.eam.alloy')\n",
    "\n",
    "MaxwellBoltzmannDistribution(...)\n",
    "Stationary(structure)\n",
    "vv_integrator = VelocityVerlet(...)  \n",
    "\n",
    "\n",
    "trajectory = []\n",
    "structure.info['time'] = 0\n",
    "structure.info['potential_energy'] = structure.get_potential_energy()\n",
    "structure.info['kinetic_energy'] = structure.get_kinetic_energy()\n",
    "trajectory.append(structure.copy())\n",
    "\n",
    "# runs for 1ps storing a structure every 20fs (already set up)\n",
    "nstep = math.ceil(20 / time_step)\n",
    "pbar = tqdm(range(int(1000/(nstep*time_step))))\n",
    "for i in pbar:  # this is the MD loop\n",
    "\n",
    "    vv_integrator.run(nstep)\n",
    "    \n",
    "    # diagnostics\n",
    "    pbar.set_description(f\"V={structure.get_potential_energy():10.5f} K={structure.get_kinetic_energy():10.5f}\")\n",
    "    structure.info['time'] = i*time_step*nstep        \n",
    "    structure.info['potential_energy'] = structure.get_potential_energy()\n",
    "    structure.info['kinetic_energy'] = structure.get_kinetic_energy()\n",
    "    trajectory.append(structure.copy())\n",
    "\n",
    "# write out the trajectory \n",
    "ase.io.write(filename, trajectory)\n",
    "return filename\n",
    "\"\"\"\n",
    "        )\n",
    "\n",
    "\n",
    "ex05_pb = ParametersBox(time_step=(5.0,1,50,1, r\"time step, fs\"),refresh_mode=\"click\")\n",
    "\n",
    "def visualise_md_aluminum(time_step, code_input, visualizers):\n",
    "    print_output = visualizers[0]\n",
    "    with print_output:\n",
    "        filename = code_input.get_function_object()(time_step)\n",
    "    if filename is None:\n",
    "        return\n",
    "    \n",
    "    frames = ase.io.read(filename, ':')\n",
    "    \n",
    "    times = [structure.info['time'] for structure in frames]\n",
    "    potential_energies = [structure.info['potential_energy'] for structure in frames]\n",
    "    kinetic_energies = [structure.info['kinetic_energy'] for structure in frames]\n",
    "    for f in frames:\n",
    "        f.info.clear()\n",
    "        f.arrays.pop('momenta')\n",
    "        \n",
    "    properties=dict(\n",
    "                time=times,\n",
    "                potential_energy=potential_energies,\n",
    "                kinetic_energy = kinetic_energies,\n",
    "                total_energy = np.asarray(potential_energies)+np.asarray(kinetic_energies)\n",
    "                )\n",
    "\n",
    "    chemiscope.write_input(\"module_06-md_aluminum_nve.chemiscope.json.gz\", frames=frames, properties=properties)\n",
    "    with print_output:\n",
    "        display(chemiscope.show(frames, properties,\n",
    "                         settings={\"structure\":[{\"bonds\":False, \"unitCell\":True,\"supercell\":{\"0\":1,\"1\":1,\"2\":1},\n",
    "                                   \"environments\": {\"cutoff\": 40}}],\n",
    "                                   'map': {'y' : {'property':'total_energy' }} } )\n",
    "               )\n",
    "\n",
    "\n",
    "ex05_code_demo = CodeDemo(\n",
    "            code_input= ex05_wci,\n",
    "            input_parameters_box=ex05_pb,\n",
    "            visualizers = [ClearedOutput()],\n",
    "            update_visualizers = visualise_md_aluminum,\n",
    "            merge_check_and_update_buttons = True)\n",
    "\n",
    "answer_registry.register_answer_widget(\"ex05-function\", ex05_code_demo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display(ex05_code_demo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Download chemiscope datafile](./module_06-md_aluminum_nve.chemiscope.json.gz)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span style=\"color:blue\">**05b** What is the largest time step you can use before the simulation becomes unstable? What are the differences relative to the hydrogen molecule above that explain the different behavior?  </span>\n",
    "\n",
    "_NB: when you approach the stability limit some trajectories will blow up completely in the short time that is simulated by the demo, and some will not. You don't have to be very precise in quoting the point where things start go awry, just give a rough estimate to within a few fs._"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ex5b_txt = TextareaAnswer(\"Enter your answer here\")\n",
    "answer_registry.register_answer_widget(\"ex5b-answer\", ex5b_txt)\n",
    "display(ex5b_txt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Constant-temperature simulations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The simulations in the previous section describe the evolution of a material at a constant total energy $H=V(\\mathbf{x}) + \\mathbf{p}^2/2m$. In most practical situations, materials are not studied in isolation, but as part of a larger setup that is in thermal equilibrium with the environment. This is even more relevant because usually a simulation supercell is meant to describe a small portion of a macroscopic sample, and one must simulate the flow of energy between different parts of the system. \n",
    "\n",
    "There are dozens of schemes that have been proposed to modify Hamilton's equations to describe the energy exchange with a heat bath. Here we will use [Langevin dynamics](https://en.wikipedia.org/wiki/Langevin_equation), in which the time evolution for $\\mathbf{p}$ is written as \n",
    "\n",
    "$$\n",
    "\\dot{\\mathbf{p}} = -\\frac{\\partial V}{\\partial\\mathbf{x}} -\\gamma \\mathbf{p} + \\sqrt{2m\\gamma k_B T} \\boldsymbol{\\xi}\n",
    "$$\n",
    "\n",
    "$-\\gamma \\mathbf{p}$ corresponds to a viscous drag term, and $\\sqrt{2m\\gamma k_B T} \\boldsymbol{\\xi}$ describes stochastic collisions that compensate the energy dissipated by friction ($\\boldsymbol{\\xi}$ corresponds to uncorrelated Gaussian noise. The net result is that the energy fluctuates, in a way that is consistent with isothermal conditions at temperature $T$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A Langevin integrator is implemented in ASE in the `ase.md.langevin` module. Its usage is very similar to that for the NVE integrator\n",
    "\n",
    "```python\n",
    "from ase.md.langevin import Langevin\n",
    "from ase import units\n",
    "\n",
    "lan = Langevin(structure, timestep=2*units.fs, temperature_K=300, friction=1/(1000*units.fs))\n",
    "\n",
    "lan.run(1000)\n",
    "```\n",
    "\n",
    "The friction is expressed in inverse time units, and $1/\\gamma$ roughly indicates the time it takes for the momenta to equilibrate at the target temperature. Note that a too large friction would slow down diffusion, and so reduce the efficiency of the simulation. A value around $1/\\text{ps}$ is usually a decent guess. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span style=\"color:blue\">**06** Write a routine that runs constant-temperature molecular dynamics for a 2x2x2 supercell of Al. Use an aggressive thermostat time constant ($1/\\gamma$) of 200fs, and a time step of 10fs. Run first 500fs to equilibrate the structure, and a further 2ps to accumulate statistics.  The function has already a stub of a MD loop that will save the trajectories that you run, for further processing. </span>\n",
    "\n",
    "_NB_: once you get something that works, run at several temperatures, e.g. 100K, 200K, 400K, 700K - the resulting trajectories will be used further down. DON'T CHANGE THE NAMING CONVENTION. Each trajectory should take a few minutes, so you can start working on the next answer while you run these."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ex06_wci = WidgetCodeInput(\n",
    "        function_name=\"aluminum_langevin\", \n",
    "        function_parameters=\"T\",\n",
    "        docstring=\"\"\"\n",
    "Performs Langevin Dynamics for the 2x2x2 supercell of fcc aluminum described by EAM potential. \n",
    ":param time_step: time step for molecular dynamic in fs. \n",
    "\"\"\",\n",
    "        function_body=\"\"\"\n",
    "from ase.md.velocitydistribution import MaxwellBoltzmannDistribution, Stationary\n",
    "from ase.md.langevin import Langevin\n",
    "from ase.calculators import eam\n",
    "from ase import units\n",
    "import ase.io\n",
    "import math\n",
    "from tqdm.notebook import tqdm\n",
    "filename = f'aluminum_langevin_T_{T}.extxyz'\n",
    "\n",
    "return filename # <-- remove this after having completed the code below\n",
    "\n",
    "structure = ...\n",
    "\n",
    "# initialize potential (don't forget to instrument the cell)\n",
    "calc = eam.EAM(potential='data/Al99.eam.alloy')\n",
    "\n",
    "# initialize the velocities at 400K, and remove the center of mass\n",
    "MaxwellBoltzmannDistribution(...)\n",
    "Stationary(structure)\n",
    "lan = Langevin(...)  \n",
    "\n",
    "print(f\"Equilibrating at {T}K (be patient, it'll take ~1 minute)\")\n",
    "lan.run(50)\n",
    "\n",
    "\n",
    "trajectory = []\n",
    "\n",
    "structure.info['index'] = 0\n",
    "structure.info['potential energy'] = structure.get_potential_energy()\n",
    "structure.info['kinetic energy'] = structure.get_kinetic_energy()\n",
    "trajectory.append(structure.copy())\n",
    "pbar = tqdm(range(200))\n",
    "for i in pbar:\n",
    "    lan.run(1)  # save at every step\n",
    "    pbar.set_description(f\"V={structure.get_potential_energy():10.5f} K={structure.get_kinetic_energy():10.5f}\")\n",
    "\n",
    "    structure.info['index'] = i\n",
    "    structure.info['potential energy'] = structure.get_potential_energy()\n",
    "    structure.info['kinetic energy'] = structure.get_kinetic_energy()\n",
    "    trajectory.append(structure.copy())\n",
    "\n",
    "ase.io.write(filename, trajectory)\n",
    "return filename\n",
    "\"\"\"\n",
    "        )\n",
    "\n",
    "ex06_pb = ParametersBox(T = (100,10,1000,10,r'Target temp.'),refresh_mode = \"click\")\n",
    "def langevin_updater(T, code_input, visualizers):\n",
    "    print_output = visualizers[0]\n",
    "    with print_output:\n",
    "        filename = code_input.get_function_object()(T)\n",
    "        frames = read(filename, \"1:\")\n",
    "        properties = chemiscope.extract_properties(frames)\n",
    "        mycs = chemiscope.show(frames, properties,\n",
    "                     settings={\"structure\":[{\"bonds\":False, \"unitCell\":True,\"supercell\":{\"0\":1,\"1\":1,\"2\":1},\n",
    "                                            \"keepOrientation\": True, 'playbackDelay': 100}],\n",
    "                               'map': {\n",
    "                                   'x' : {'property':'index' },\n",
    "                                   'y' : {'property':'potential energy' }\n",
    "                               } }\n",
    "                               )\n",
    "        display(mycs)\n",
    "\n",
    "ex06_code_demo = CodeDemo(\n",
    "            code_input= ex06_wci,\n",
    "            input_parameters_box=ex06_pb,\n",
    "            visualizers = [ClearedOutput()],\n",
    "            update_visualizers = langevin_updater,\n",
    "            merge_check_and_update_buttons = True)\n",
    "\n",
    "answer_registry.register_answer_widget(\"ex06-function\", ex06_code_demo)\n",
    "display(ex06_code_demo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Termal fluctuations cause vibrations of the atoms in a crystal around their equilibrium positions. These vibrations have all sorts of effects, from thermal expansion to changes of the electrical conductivity. One of the most direct consequences of the atomic motion is the attenuation of the intensity of diffraction peaks, quantified by a term known as the [Debye-Waller factor](https://en.wikipedia.org/wiki/Debye%E2%80%93Waller_factor). \n",
    "\n",
    "The central quantity is the _mean square displacement_ (MSD)\n",
    "\n",
    "$$\n",
    "\\langle u^2 \\rangle = \\frac{1}{N} \\langle \\sum_i (\\mathbf{x}_i - \\mathbf{x}_i^{(0)})^2\\rangle\n",
    "$$\n",
    "\n",
    "where $\\mathbf{x}_i^{(0)}$ indicates the equilibrium lattice position for the atom $i$, and $\\mathbf{x}_i$ the position of the same atom in a snapshot of a finite-temperature simulation. $N$ is the number of atoms in the supercell.\n",
    "\n",
    "$\\langle \\cdot \\rangle$ indicates averaging over thermal fluctuations, but we can compute it (through the so-called [ergodic hypothesis](https://en.wikipedia.org/wiki/Ergodic_hypothesis)) running a finite-temperature molecular dynamics simulations, and collecting representative snapshots of the structure. In the previous demo we ran several of these finite-temperature trajectories, so now we can load them, and use them to compute the MSD at various temperatures. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span style=\"color:blue\">**07** Write a routine that, given a pre-computed trajectory, evaluates the mean-square displacement of the atoms and returns its value. \n",
    "Compute the mean of the positions at all time steps along the trajectory as an estimate of the equilibrium reference position $\\mathbf{x}^{(0)}_i$, so that the function can be self-contained, and does not require building the idealized minimum-energy structure. </span>\n",
    "\n",
    "_Hints:_\n",
    "* Note that the MSD is a scalar - this function should return a single value for each `filename` it loads\n",
    "* You can achieve the appropriate averaging with explicit summations, but you can also more easily (and quickly) achieve the same goal with a judicious use of [`numpy.mean`](https://numpy.org/doc/stable/reference/generated/numpy.ndarray.mean.html)\n",
    "* As you could show using statistical mechanics, for a harmonic crystal the MSD should grow linearly with $T$. Here you may observe small deviations due to anharmonicity, and noise due to the limited length of the simulations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ex07_wci = WidgetCodeInput(\n",
    "        function_name=\"compute_msd\", \n",
    "        function_parameters=\"filename\",\n",
    "        docstring=\"\"\"\n",
    "Computes the mean square displacement of the positions of the atoms over the trajectory. \n",
    ":param filename: path to the saved trajectory in the extxyz format. \n",
    "\"\"\",\n",
    "        function_body=\"\"\"\n",
    "\n",
    "import ase.io\n",
    "import numpy as np\n",
    "\n",
    "structures = ase.io.read(filename, ':')\n",
    "return 0.0 # <-- remove this after having completed the code below\n",
    "return ...\n",
    "\"\"\"\n",
    "        )\n",
    "\n",
    "def plot_msd(code_input,visualizers):\n",
    "    print_output = visualizers[0]\n",
    "    ax = visualizers[1].figure.get_axes()[0]\n",
    "    import os\n",
    "    files = os.listdir('.')\n",
    "    func = code_input.get_function_object()            \n",
    "    T_grid, msd_values = [], []\n",
    "    for file in files:\n",
    "        if file.startswith('aluminum_langevin_T_'):\n",
    "            try:\n",
    "                T = float(file.split('.')[0].split('_')[-1])\n",
    "            except:\n",
    "                pass            \n",
    "            \n",
    "            with print_output:\n",
    "                msd = func(file)\n",
    "            T_grid.append(T)\n",
    "            msd_values.append(msd)\n",
    "            \n",
    "    T_grid, msd_values = np.array(T_grid), np.array(msd_values)\n",
    "    indices = np.argsort(T_grid)\n",
    "    T_grid = T_grid[indices]\n",
    "    msd_values = msd_values[indices]\n",
    "    \n",
    "    ax.plot(T_grid, msd_values, 'o', color = 'blue')\n",
    "    ax.plot(T_grid, msd_values, color = 'blue')\n",
    "    ax.set_xlabel('$T$ / K')\n",
    "    ax.set_ylabel('\\n\\n\\nmsd/ Å$^2$')\n",
    "                     \n",
    "ex07_figure,_ = plt.subplots(1, figsize=(6,3.8), tight_layout=True)\n",
    "ex07_pyplot_output = PyplotOutput(ex07_figure)\n",
    " \n",
    "ex07_code_demo = CodeDemo(\n",
    "            code_input= ex07_wci,\n",
    "            check_registry=check_registry,\n",
    "            visualizers = [ClearedOutput(),ex07_pyplot_output],\n",
    "            update_visualizers = plot_msd,\n",
    "            merge_check_and_update_buttons=False\n",
    "           )\n",
    "\n",
    "check_registry.add_check(ex07_code_demo,\n",
    "                         inputs_parameters=[{\"filename\": \"module_06-ex07_ref.xyz\"}],\n",
    "                         reference_outputs = [0.008466126604490386*3], \n",
    "                         equal=np.allclose\n",
    "                        )     \n",
    "\n",
    "answer_registry.register_answer_widget(\"ex07-function\", ex07_code_demo)\n",
    "\n",
    "display(ex07_code_demo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Melting point calculations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we run something a bit more complicated, in fact pushing the limits of what can be achieved with ASE and pure Python code. We will run a simulation in which Al is heated above its melting point, and then quenched back to low temperature. \n",
    "\n",
    "This can be achieved with a loop that adjusts the target temperature of the thermostat between a low value $T_{\\mathrm{low}}$ and a high value $T_{\\mathrm{high}}$, running in between short stretches of MD simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_ramp(nramp, temp_lo, temp_hi,visualizers):\n",
    "    ax = visualizers[1].figure.get_axes()[0]\n",
    "    ax.plot(np.arange(nramp*2), \n",
    "            np.concatenate([\n",
    "            np.linspace(temp_lo, temp_hi, nramp),\n",
    "            np.linspace(temp_hi, temp_lo, nramp)\n",
    "            ]) )\n",
    "    ax.set_xlabel(\"MD segment\")\n",
    "    ax.set_ylabel(r'$T$/K')\n",
    "    \n",
    "ramp_pb = ParametersBox(nramp = (200,50,4000,50, r'$n_{\\mathrm{ramp}}$'),\n",
    "                       temp_lo = (350.,100,1000,50, r'$T_{\\mathrm{low}}/K$'),\n",
    "                       temp_hi = (3500.,1500,4000,500, r'$T_{\\mathrm{high}}/K$'),\n",
    "                      )\n",
    "\n",
    "ramp_figure,_ = plt.subplots(1, figsize=(8,3.8), tight_layout=True)\n",
    "ramp_pyplot_output = PyplotOutput(ramp_figure)\n",
    "\n",
    "ramp_demo = CodeDemo(\n",
    "            input_parameters_box=ramp_pb,\n",
    "            visualizers = [ClearedOutput(),ramp_pyplot_output],\n",
    "            update_visualizers = plot_ramp)\n",
    "display(ramp_demo)\n",
    "ramp_demo.run_demo()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function below runs a melt-and-quench simulation. Try to understand it as much as you can: it uses several non-trivial algorithms, such as constant-pressure integrator that allow to change the size of the simulation cell in response to the internal pressure, which is needed to accomodate thermal expansion and the latent volume of fusion. It is also _slow_, because everything is implemented in Python: start by running with the default parameters, keeping in mind it might already take up to one hour to complete the simulation.\n",
    "At the end of the simulation, the routine saves a `.xyz` file containing the simulation \"movie\". We will look at it further down. Later you may want to run more simulations, changing tre supercell size, the length of the ramp, and the temperature bounds: if you experiment with multiple parameters, change the `filename` so that you can load the trajectory later. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ex08_wci = WidgetCodeInput(\n",
    "        function_name=\"melt_and_quench\", \n",
    "        function_parameters=\"nrep, nramp, temp_lo, temp_hi\",\n",
    "        docstring=\"\"\"\n",
    "Creates a nrep×nrep×nrep Al fcc structure, and runs a MD simulation in which the temperature \n",
    "is raised from temp_low to temp_hi over nramp short trajectory segments. \n",
    ":return: name of the trajectory file\n",
    "Structures to be visualized and line energies estimated for the structures along the relaxation\n",
    "\"\"\",\n",
    "        function_body=\"\"\"\n",
    "import numpy as np\n",
    "from ase.io import write\n",
    "from ase import Atoms\n",
    "from ase.calculators import lj, eam\n",
    "from ase.md.nptberendsen import NPTBerendsen\n",
    "from ase.md.langevin import Langevin\n",
    "from ase import units\n",
    "from tqdm.notebook import tqdm\n",
    "\n",
    "# use a LJ potential with parameters fitted to match lattice parameter and cohesive energy of Al. aggressive cutoff to reduce cost\n",
    "calc = lj.LennardJones(sigma=2.62, epsilon=0.41, rc=2*2.62)\n",
    "#calc = eam.EAM(potential='data/Al99.eam.alloy')   #<<- this is way too slow, but you can try it\n",
    "a0 = 4.05 # lattice parameter of Al\n",
    "# constructs a unit cell\n",
    "fcc_cell = Atoms(\"Al4\", cell=np.eye(3)*a0, positions=a0*np.asarray([[0,0,0],[0.5,0.5,0],[0.5,0,0.5],[0,0.5,0.5]] ),\n",
    "                 pbc=True, calculator=calc ) \n",
    "\n",
    "# creates a supercell\n",
    "suxcell = fcc_cell.repeat(nrep)\n",
    "# attach calculator\n",
    "suxcell.calc = calc\n",
    "\n",
    "# this is a long story. ASE doesn't have a decent barostat, NPT is very unstable, and Berendsen \n",
    "# doesn't sample the correct ensemble. so, we use a combination of Berendsen with a loose coupling and \n",
    "# a Langevin integrator to enforce canonical sampling\n",
    "dyn = NPTBerendsen(suxcell, timestep=2*units.fs, temperature_K=temp_lo, pressure_au=2*units.bar, \n",
    "                   taut=1e4*units.fs, taup=1e2*units.fs, compressibility_au=2.2 )\n",
    "lan = Langevin(suxcell, timestep=2*units.fs, temperature_K=temp_lo, friction=0.02)\n",
    "\n",
    "filename = \"traj-exx.xyz\"\n",
    "\n",
    "# prints some stats\n",
    "from time import time\n",
    "start = [0.0]\n",
    "iramp = [1]\n",
    "traj = []\n",
    "\n",
    "# ugly but effective: store reference to globals in the default args\n",
    "def printenergy(atoms=suxcell, t=traj, md=dyn, start=start, nramp=nramp, iramp=iramp):  \n",
    "    elapsed = time() - start[0]\n",
    "    epot = atoms.get_potential_energy() / len(atoms)\n",
    "    ekin = atoms.get_kinetic_energy() / len(atoms)\n",
    "    pressure = np.trace(atoms.get_stress(voigt=False))/3\n",
    "    a = atoms.copy(); a.calc = atoms.calc\n",
    "    volume = np.linalg.det(a.cell)\n",
    "    a.arrays.pop('momenta')\n",
    "    a.info['time'] = md.dt*md.get_number_of_steps()/(1000*units.fs)\n",
    "    a.info['target_temperature'] = dyn.temperature\n",
    "    a.info['potential'] = epot\n",
    "    a.info['kinetic_temperature'] = ekin / (1.5 * units.kB)\n",
    "    a.info['stress'] = pressure\n",
    "    a.info['volume'] = volume\n",
    "    t.append(a)\n",
    "    print('Energy/at.: V = %.3feV  K = %.3feV (T=%3.0fK)  '\n",
    "          'Volume = %.3fÅ³, Time (elapsed/total): %.3fs/%.3fs' % (epot, ekin, ekin / (1.5 * units.kB), volume, elapsed, elapsed*(2*nramp/iramp[0]-1)))\n",
    "\n",
    "start[0] = time()\n",
    "print(\"Starting equilibration (be patient)\")\n",
    "lan.run(400)\n",
    "\n",
    "# attaches the function as a callback - it'll be called every 50 steps to store and print the trajectory\n",
    "dyn.attach(printenergy, interval=50)\n",
    "\n",
    "# temperature change between ramp steps\n",
    "framp = (temp_hi-temp_lo)*(1./nramp)\n",
    "# ramp temperature up\n",
    "print('Starting temperature ramp')\n",
    "for i in tqdm(range(nramp)):\n",
    "    iramp[0]+=1\n",
    "    lan.run(75)\n",
    "    dyn.run(25)\n",
    "    suxcell.wrap()\n",
    "    lan.set_temperature(temperature_K=lan.temp/units.kB+framp)\n",
    "    dyn.set_temperature(temperature_K=lan.temp/units.kB+framp)\n",
    "# and then quench\n",
    "print('Quench')\n",
    "for i in tqdm(range(nramp)):\n",
    "    iramp[0]+=1\n",
    "    lan.run(75)    \n",
    "    dyn.run(25)\n",
    "    suxcell.wrap()\n",
    "    lan.set_temperature(temperature_K=lan.temp/units.kB-framp)\n",
    "    dyn.set_temperature(temperature_K=lan.temp/units.kB-framp)   \n",
    "\n",
    "write(filename, traj)\n",
    "return filename\n",
    "\"\"\"\n",
    "        )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ex08_pb = ParametersBox(nrep = (2,1,3,1,r'$n_{\\mathrm{rep}}$'), \n",
    "                       nramp = (100,50,4000,50, r'$n_{\\mathrm{ramp}}$'),\n",
    "                       temp_lo = (400.,100,1000,50, r'$T_{\\mathrm{low}}/K$'),\n",
    "                       temp_hi = (4000.,1500,5000,500, r'$T_{\\mathrm{high}}/K$'),\n",
    "                       refresh_mode=\"click\"\n",
    "                      )\n",
    "def ex08_updater(nrep,nramp,temp_lo,temp_hi,code_input, visualizers):\n",
    "    print_output = visualizers[0]\n",
    "    with print_output:\n",
    "        filename = code_input.get_function_object()(nrep,nramp,temp_lo,temp_hi)    \n",
    "\n",
    "ex08_code_demo = CodeDemo(\n",
    "            code_input= ex08_wci,\n",
    "            input_parameters_box=ex08_pb,\n",
    "            visualizers = [ClearedOutput()],\n",
    "            update_visualizers = ex08_updater)\n",
    "\n",
    "answer_registry.register_answer_widget(\"ex08-function\", ex08_code_demo)\n",
    "display(ex08_code_demo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here you can load the output trajectory and visualize it. Time is expressed in picoseconds ($10^{-12}$s), temperatures in K, volumes in Å³ and energies in eV/atom. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_traj(filename,visualizers):\n",
    "    clear_output = visualizers[0]\n",
    "    \n",
    "    traj = read(filename, \":\")\n",
    "    for t in traj:\n",
    "        t.arrays.pop('forces', -1)\n",
    "        t.arrays.pop('momenta', -1)\n",
    "        t.info.pop('stress', -1)\n",
    "    properties = {'temperature / K': {'property' : 'target_temperature','target':'structure','values':[t.info[\"target_temperature\"] for t in traj]}, \n",
    "                              'energy / eV': {'property' : 'potential','target':'structure','values':[t.info[\"potential\"] for t in traj]},\n",
    "                              'time / ps':{'property': 'time','target':'structure','values':[t.info[\"time\"] for t in traj]}}\n",
    "    \n",
    "    chemiscope.write_input(\"./module-06_temperature_ramp.chemiscope.json.gz\", frames=traj,\n",
    "                          properties=properties)\n",
    "    clear_output.clear_output()  \n",
    "    with clear_output:              \n",
    "        cs = chemiscope.show(traj, properties=properties, \n",
    "                            settings={'structure': [{'bonds': False, 'unitCell': True, 'keepOrientation': True}],\n",
    "                                     'map': {'joinPoints': True} }\n",
    "                            )\n",
    "        display(cs)\n",
    "\n",
    "load_traj_demo = CodeDemo(\n",
    "            visualizers = [ClearedOutput()],\n",
    "            input_parameters_box=ParametersBox(filename = Text(\"traj-exx.xyz\"),refresh_mode = \"click\"),\n",
    "            update_visualizers = load_traj)\n",
    "display(load_traj_demo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Download chemiscope datafile](./module-06_temperature_ramp.chemiscope.json.gz)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Inspect the trajectory. Plot target temperature versus simulation time to visualize the temperature ramp. Then plot potential (a proxy for the enthalpy) versus time. It is also instructive to look at the potential energy as a function of the target temperature, which is the default visualization you get when you load a trajectory. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span style=\"color:blue\">**08a** What happens as the temperature approaches the maximum temperature along the ramp? And as the temperature is decreased? Is the potential curve symmetric in the heating and cooling directions? Compare the starting and final configurations: what differences do you observe? </span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ex8a_txt = TextareaAnswer(\"Enter your answer here\")\n",
    "answer_registry.register_answer_widget(\"ex8a-answer\", ex8a_txt)\n",
    "display(ex8a_txt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Besides the fact we are running the simulation with a very lousy model of the interatomic potential for Al (as we have seen in the [potentials module](./04-Potentials.ipynb)) there are two serious limitations to these simulations: they are too small (they suffer from _finite size effects_) and they are too fast (they are strongly out-of-equilibrium). Given that larger and longer simulations are too lengthy, we have prepared some for you to inspect. You can load a 100ps, 3×3×3 run inserting in the input box above the filename `data/traj-n3-r1000.xyz`, and an even longer, 400ps trajectory loading `data/traj-n3-r4000.xyz`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span style=\"color:blue\">**08b** In these trajectories you can observe a clearer discontinuity in the potential (and the volume) as the temperature increases and decreases. Are these discontinuities at the same temperature? How can you explain this observation? What would be your best estimate for the melting point of Al with this model potential? </span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ex8b_txt = TextareaAnswer(\"Enter your answer here\")\n",
    "answer_registry.register_answer_widget(\"ex8b-answer\", ex8b_txt)\n",
    "display(ex8b_txt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span style=\"color:blue\">**08c** Look carefully at the final configuration: can you give a better explanation for the difference in energy between the starting and the final states of the trajectory? Is this a completely unreasonable artefact? Can you think of realistic processing conditions that would lead to similar phenomena? </span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ex8c_txt = TextareaAnswer(\"Enter your answer here\")\n",
    "answer_registry.register_answer_widget(\"ex8c-answer\", ex8c_txt)\n",
    "display(ex8c_txt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A more challenging exercise (optional!)\n",
    "\n",
    "The problem of the celestial mechanics of two bodies can be solved analytically with comparatively elementary tools. \n",
    "However the problem gets _substantially_ harder when you consider the dynamics of three objects of comparable mass \n",
    "(the so-called [three-body problem](https://en.wikipedia.org/wiki/Three-body_problem)).\n",
    "\n",
    "Implement a simulation of three particles interacting through gravity, and compare the trajectories obtained for \n",
    "different time steps, and different integrators. The errors will be much more dramatic, as even with three particles\n",
    "you start seeing exponential deviation between trajectories that differ by an infinitesimal amount \n",
    "(read about it in terms of the [Lyuapunov exponent](https://en.wikipedia.org/wiki/Lyapunov_exponent)). \n",
    "The reason we can get away with simulating _thousands_ of atoms with an approximate integrator is that \n",
    "typically we don't care about the actual trajectories (which are kind of important if you want to predict\n",
    "if a meteor will crash into the Earth...) but only about their statistical properties, that are typically\n",
    "not affected by moderate integration errors. "
   ]
  }
 ],
 "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.12.3"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {
    "height": "160px",
    "width": "252px"
   },
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "461px",
    "left": "0px",
    "right": "927.667px",
    "top": "107px",
    "width": "139px"
   },
   "toc_section_display": "block",
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}