notebook/statistical-mechanics/monte_carlo_parabolic.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4a374f26",
   "metadata": {},
   "source": [
    "# Monte Carlo Simulations to Sample the Canonical Distribution\n",
    "\n",
    "**Authors:** Dou Du, Taylor James Baird and Giovanni Pizzi\n",
    "\n",
    "<i class=\"fa fa-home fa-2x\"></i><a href=\"../index.ipynb\" style=\"font-size: 20px\"> Go back to index</a>\n",
    "\n",
    "**Source code:** https://github.com/osscar-org/quantum-mechanics/blob/master/notebook/statistical-mechanics/monte_carlo_parabolic.ipynb\n",
    "\n",
    "We demonstrate the Monte Carlo simulations sampling the canoical distribution in this notebook.\n",
    "\n",
    "<hr style=\"height:1px;border:none;color:#cccccc;background-color:#cccccc;\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4eb88314",
   "metadata": {},
   "source": [
    "## Goals\n",
    "* Understand how a Monte Carlo simulation samples the canonical distribution (with simple examples).\n",
    "* Appreciate the effect of temperature on the simulations.\n",
    "* Familiarize oneself with the role of numerical parameters (number of steps, step size).\n",
    "* Gain insight into the difficulties encountered when considering low temperatures, high barriers, and achieving ergodicity."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63440ccb",
   "metadata": {},
   "source": [
    "## **Background theory** \n",
    "\n",
    "[More on the background theory.](./theory/theory_monte_carlo_parabolic.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "180723bd",
   "metadata": {},
   "source": [
    "## Tasks and exercises\n",
    "\n",
    "1. Find the approximate <code>(x, y)</code> coordinates of the minima of the potential both from the potential energy surface plot and from the peaks of the probability distribution (after running the simulation with the default parameters).\n",
    "   Do this for both potentials (single and double well). Note that when you hover on the potential energy surface plot, you can also see some black isosurface lines, that can help you find where the minimum is.<br>\n",
    "\n",
    "    <details>\n",
    "    <summary style=\"color: red\">Solution</summary>\n",
    "\n",
    "    There is only one global minimum at <code>(x=0, y=0)</code> in the single parabolic well case.<br>\n",
    "\n",
    "    In the double well case, two global minima are located at around <code>(x=0, y=3.13)</code> and <code>(x=0, y=-3.13)</code>.<br>\n",
    "    \n",
    "    In the probability distribution plot, if you run with the default variables, you should see clear peaks around the minima.<br><br>\n",
    "    \n",
    "    </details>\n",
    "\n",
    "1. Select to the single-well case and use the slider to modify the temperature. What do you expect the converged average total energy to be? Is this the case?\n",
    "\n",
    "    <details>\n",
    "    <summary style=\"color: red\">Solution</summary>\n",
    "    For a single quantum well with two coordinates (x and y), we have two degrees of freedom. Therefore, by equipartition, we expect the average energy to be $2\\cdot \\frac {k_B T}{2}=k_B T$ (i.e., in the reduced units used here with $k_B=1$, the same numerical value as the temperature).\n",
    "    As you change the temperature (try e.g. 0.1, 1 and 5), you should see that the average energy tends to reach the same value as the temperature.\n",
    "    Note: in order to keep the equilibration period very short, it might be convenient to set the starting point to <code>(0, 0)</code>. Otherwise, note you can zoom in on the plot to better inspect the values.<br>\n",
    "    \n",
    "    You will not necessarily converge to the exact value: try to run the simulation more than once, or try to increase the simulation steps or tune the other numerical parameters (we will do this together in the tasks below).<br><br>\n",
    "    \n",
    "    </details>\n",
    "\n",
    "1. Select to the single well case. Set the starting position to, e.g., (10, 10) to be out of the minimum, and the temperature to 1. Change the max move size (for example, try values around 0.1, 1 and 5). Inspect the convergence of the total energy to the converged result, both in terms of how long the equilibration takes, and the statistical error.\n",
    "\n",
    "    <details>\n",
    "    <summary style=\"color: red\">Solution</summary>\n",
    "    You should notice that, with a smaller max move size, the simulation takes a longer time to reach convergence. You can for example check how many steps it takes to get within, say, 5% of the converged value (you can zoom in, and note that the y axis is automatically rescaled).\n",
    "    \n",
    "    On the other hand, while a bigger move size makes the simulation equilibrate quickly, the magnitude of the fluctuations remain larger. In addition, the acceptance rate drops significantly: most moves change the system coordinates by a large amount, increasing significantly the energy of the system, and only a few of these are statistically accepted. The acceptance rate is a typical parameter that one should monitor and tune in a simulation; if it is too large it indicates that we are accepting almost all moves (and probably then the move is small and we are far from the minimum); if it is too small, we are probably moving too much at every step; as a result the simulation becomes much more expensive (it has to perform many calculations and reject most of them) and the average quantities are affected by larger errors.<br><br>\n",
    "    \n",
    "    </details>\n",
    "\n",
    "1. Let us now investigate the double-well potential, with low temperature and small max move size. You can set the starting position to, e.g., (10, 10), the temperature to 0.5, and the max move size to 0.1. Run the simulation a few times. How many minima does the simulation typically samples? What happens if you increase the max move size? What happens for higher temperatures?\n",
    "\n",
    "    <details>\n",
    "    <summary style=\"color: red\">Solution</summary>\n",
    "    When the temperature is small, the probability to jump between the two wells across the barrier is very low, and might never happen in the simulation. In this case the simulation is not ergodic, because the distribution obtained weighting the phase-space points (that would occupy identically the two wells) is different from the time average over the finite simulation (that typically sees only one well).\n",
    "    \n",
    "    By starting at <code>(10, 10)</code> and using a small step, we typically end up almost always in the closest of the two wells. If we start at the top of the double well (starting coordinates <code>(0, 0)</code>), with a small temperature (1 or smaller) and small max move size, each simulation will randomly fall in one of the two wells and be trapped there.\n",
    "    \n",
    "    If you increase the move size, the fluctuations become larger and one gets a better representation of the canonical distribution (e.g. for a max move size of 10) at the expense of larger fluctuations in the averaged quantities). In addition, most moves are rejected because the majority of them bring the system to a much higher energy.\n",
    "    \n",
    "    As the temperature increases, the thermal fluctuations increase and there is a higher probability that they allow us to visit both wells, if a long-enough simulation is considered.<br><br>\n",
    "    \n",
    "    </details>\n",
    "\n",
    "<hr style=\"height:1px;border:none;color:#cccccc;background-color:#cccccc;\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "586c902d",
   "metadata": {},
   "source": [
    "## Interactive visualization\n",
    "(be patient, it might take a few seconds to load)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5a3bc18",
   "metadata": {},
   "outputs": [],
   "source": [
    "import plotly.graph_objects as go\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import random\n",
    "import math\n",
    "from plotly.subplots import make_subplots\n",
    "from ipywidgets import FloatSlider, IntSlider, FloatLogSlider, Button\n",
    "from ipywidgets import Output, HBox, VBox, IntProgress, Dropdown, Layout, Label"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f8782b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "#all the widget components\n",
    "\n",
    "description_width = '120px'\n",
    "\n",
    "temp_slider = FloatSlider(min=0.1, max=5, value=1., continuous_update=False, description=\"Temperature\", style={'description_width': description_width})\n",
    "move_size_slider = FloatLogSlider(min=-1, max=1, value=1., continuous_update=False, description=\"Max move size\", style={'description_width': description_width})\n",
    "num_iterations_slider = IntSlider(min=100000, max=1000000, value=100000, continuous_update=False, description=\"Num iterations\", style={'description_width': description_width})\n",
    "run_button = Button(description=\"Run Monte-Carlo\")\n",
    "run_button.style.button_color = 'green'\n",
    "traces_button = Button(description=\"Show Traces\")\n",
    "x_slider = FloatSlider(value=4.0, min=-10.0, max=10.0, continuous_update=False, description=\"Initial pos. x\", style={'description_width': description_width})\n",
    "y_slider = FloatSlider(value=4.0, min=-10.0, max=10.0, continuous_update=False, description=\"Initial pos. y\", style={'description_width': description_width})\n",
    "potential_dropdown = Dropdown(options=[('single parabolic well',1), ('double well',2)],\n",
    "                              value = 1, description='Potential:', layout=Layout(width='250px'))\n",
    "label1 = Label(value='(only shows 300 uniformly sampled steps in the whole simulation)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b88f2156",
   "metadata": {},
   "outputs": [],
   "source": [
    "A1 = 1. # potential curvature first axis\n",
    "A2 = 1.e-1 # potential curvature second axis\n",
    "\n",
    "# Starting point of our MC, for the two coordinates\n",
    "starting_x1 = 9.0\n",
    "starting_x2 = 9.0\n",
    "\n",
    "def get_energy(x, y):\n",
    "    \"\"\"Compute the potential energy according to the dropbox.\"\"\"\n",
    "    if potential_dropdown.value == 1:\n",
    "        return 0.5 * x**2 + 0.5 * 4.0 * y**2\n",
    "    elif potential_dropdown.value == 2:\n",
    "        return 0.5 * A1 * x**2 + 0.5 * A2 * (y**4 - 20 * y**2)\n",
    "\n",
    "def get_property(x1, x2):\n",
    "    # I can do something else, I decide to sample the energy of one of the two coordinates\n",
    "    #return 0.5 * A2 * (x2**4 - 4 * x2**2)\n",
    "    \n",
    "    # I just return the second coordinate, to see its average\n",
    "    return x2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6cfb4f23",
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = np.mgrid[-10:10:100j, -10:10:100j]\n",
    "z = get_energy(x, y)\n",
    "\n",
    "# figure one to show the potential energy surface and current position\n",
    "fig1 = go.FigureWidget(data=[go.Surface(z=z, x=x, y=y, opacity=0.5),\n",
    "                             go.Scatter3d(x=[starting_x1], \n",
    "                                          y=[starting_x2], \n",
    "                                          z=[get_energy(starting_x1, starting_x2)], \n",
    "                                          mode='markers+text',\n",
    "                                          text=['Current position'],\n",
    "                                          marker=dict(size=8, color='red'))])\n",
    "\n",
    "def on_potential_change(c):\n",
    "    \"\"\"Callback function for the change of the dropdown.\"\"\"\n",
    "    fig1.data[0].z = get_energy(x, y)\n",
    "    fig1.data[1].x = [x_slider.value]\n",
    "    fig1.data[1].y = [y_slider.value]\n",
    "    fig1.data[1].z = [get_energy(x_slider.value, y_slider.value)]\n",
    "    \n",
    "potential_dropdown.observe(on_potential_change)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68b8534b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# figure 2 is the histogram of probabilities for the (x,y) coordinates, displayed as a heatmap.\n",
    "fig2 = go.FigureWidget()\n",
    "\n",
    "fig2.add_trace(go.Scatter(\n",
    "    x=[starting_x1],\n",
    "    y=[starting_x2],\n",
    "    mode='lines+markers',\n",
    "    visible = False,\n",
    "    showlegend=False,\n",
    "    marker=dict(\n",
    "        symbol='circle',\n",
    "        opacity=0.7,\n",
    "        color='white',\n",
    "        size=5,\n",
    "        line=dict(width=1),\n",
    "    ),\n",
    "))\n",
    "\n",
    "fig2.add_trace(go.Histogram2d(\n",
    "    x=[],\n",
    "    y=[],\n",
    "    histnorm='probability',\n",
    "    autobinx = False,\n",
    "    xbins=dict(start=-10, end=10, size=0.1),\n",
    "    autobiny = False,\n",
    "    ybins=dict(start=-10, end=10, size=0.1),\n",
    "    colorscale=[[0, 'rgb(12,51,131)'], [0.25, 'rgb(10,136,186)'], [0.6, 'rgb(242,211,56)'], \n",
    "                [0.75, 'rgb(242,143,56)'], [1, 'rgb(217,30,30)']]\n",
    "))\n",
    "\n",
    "fig1.update_layout(title='Potential energy surface', autosize=False,\n",
    "                  width=420, height=420,\n",
    "                  margin=dict(l=10, r=10, b=30, t=30))\n",
    "\n",
    "fig2.update_layout(title='Probability histogram',\n",
    "    xaxis=dict( ticks='', showgrid=False, zeroline=False, nticks=20 ),\n",
    "    yaxis=dict( ticks='', showgrid=False, zeroline=False, nticks=20 ),\n",
    "    xaxis_title=\"x\",\n",
    "    yaxis_title=\"y\",\n",
    "    autosize=False,\n",
    "    height=420,\n",
    "    width=430,\n",
    "    margin=dict(l=30, r=30, b=60, t=60),\n",
    "    hovermode='closest',\n",
    "\n",
    ")\n",
    "\n",
    "fig3 = go.FigureWidget()\n",
    "\n",
    "fig3.add_trace(go.Scatter(x=[], y=[], mode='lines',\n",
    "    name='Total energy',\n",
    "    line=dict(color='red', width=2),\n",
    "    connectgaps=True,\n",
    "))\n",
    "\n",
    "fig3.update_layout(\n",
    "    height = 200,\n",
    "    width = 530,\n",
    "    xaxis_title = 'Monte-Carlo move',\n",
    "    yaxis_title = 'Total energy',\n",
    "    margin=dict(l=30, r=30, b=10, t=10),\n",
    ")\n",
    "\n",
    "display(HBox([fig1, fig2]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d886aac",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def run_mc(move_size, num_iterations, temp):\n",
    "    \"\"\"\n",
    "    :param move_size: max move size for the two axes x1 and x2.\n",
    "        In principle, the move size could/should be different for x1 and x2, here I choose the same for simplicity!\n",
    "        \n",
    "    :param num_iterations: total number of Monte-Carlo iterations\n",
    "    \"\"\"\n",
    "    x1 = starting_x1\n",
    "    x2 = starting_x2\n",
    "    energy = get_energy(x1, x2)\n",
    "\n",
    "    all_x1 = [x1]\n",
    "    all_x2 = [x2]\n",
    "    all_properties = [get_property(x1, x2)]\n",
    "    all_energies = [get_energy(x1, x2)]\n",
    "    count_accepted = 0\n",
    "    count_refused = 0\n",
    "\n",
    "    for iter_cnt in range(num_iterations):\n",
    "        # shift by a random value between -move_size and + move_size\n",
    "        new_x1 = x1 + ((random.random() - 0.5) * 2) * move_size\n",
    "        new_x2 = x2 + ((random.random() - 0.5) * 2) * move_size\n",
    "        new_energy = get_energy(new_x1, new_x2)\n",
    "\n",
    "        if new_energy < energy:\n",
    "            accepted = True\n",
    "        else:\n",
    "            # boltzmann_k = 1 in these units\n",
    "            probability = math.exp(-(new_energy - energy)/temp)\n",
    "            # random.random() is a random number between 0 and 1\n",
    "            # also probability is between 0 and 1\n",
    "            # so if I accept only if random.random() is < probability,\n",
    "            # I'm accepting the move with 'probability' probability\n",
    "            accepted = random.random() < probability\n",
    "\n",
    "        if accepted:\n",
    "            x1 = new_x1\n",
    "            x2 = new_x2\n",
    "            energy = new_energy\n",
    "            count_accepted += 1\n",
    "        else:\n",
    "            count_refused += 1\n",
    "        all_properties.append(get_property(x1, x2))\n",
    "        all_energies.append(energy)\n",
    "        all_x1.append(x1)\n",
    "        all_x2.append(x2)\n",
    "        \n",
    "    average_integrated_energy = np.cumsum(all_energies) / (1 + np.arange(len(all_energies)))\n",
    "    average_property = np.cumsum(all_properties) / (1 + np.arange(len(all_properties)))\n",
    "\n",
    "    return count_accepted, count_refused, average_integrated_energy, average_property, all_x1, all_x2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04718a3d",
   "metadata": {},
   "outputs": [],
   "source": [
    "controls = VBox([x_slider, y_slider, potential_dropdown, temp_slider, move_size_slider, num_iterations_slider])\n",
    "\n",
    "display(HBox([controls, fig3]))\n",
    "display(run_button, HBox([traces_button, label1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4af7ade",
   "metadata": {},
   "outputs": [],
   "source": [
    "def change_init_position(c):\n",
    "    \"\"\"Callback function for the change of the x, y sliders.\"\"\"\n",
    "    global starting_x1, starting_x2\n",
    "    starting_x1 = x_slider.value\n",
    "    starting_x2 = y_slider.value\n",
    "    fig1.data[1].x = [x_slider.value]\n",
    "    fig1.data[1].y = [y_slider.value]\n",
    "    fig1.data[1].z = [get_energy(x_slider.value, y_slider.value)]\n",
    "    \n",
    "x_slider.observe(change_init_position, names='value')\n",
    "y_slider.observe(change_init_position, names='value')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc104b7d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def interactive_plot(button):\n",
    "    \"\"\"Callback function for the run button.\"\"\"\n",
    "    global all_x1, all_x2\n",
    "    temp = temp_slider.value\n",
    "    move_size=move_size_slider.value\n",
    "    num_iterations=num_iterations_slider.value\n",
    "    \n",
    "    run_button.disabled = True\n",
    "    run_button.style.button_color = 'red'\n",
    "    run_button.description = \"Running...\"\n",
    "    traces_button.disabled = True\n",
    "    try:\n",
    "        count_accepted, count_refused, average_integrated_energy, average_property, all_x1, all_x2 = run_mc(\n",
    "            move_size=move_size, num_iterations=num_iterations, temp=temp)\n",
    "        \n",
    "        fig1.data[1].x = [all_x1[-1]]\n",
    "        fig1.data[1].y = [all_x2[-1]]\n",
    "        fig1.data[1].z = [get_energy(all_x1[-1], all_x2[-1])]\n",
    "    \n",
    "        snum = int(len(all_x1)/300)\n",
    "        fig2.data[0].x = all_x1[::snum]\n",
    "        fig2.data[0].y = all_x2[::snum]\n",
    "        \n",
    "        fig2.data[1].x = all_x1\n",
    "        fig2.data[1].y = all_x2\n",
    "        \n",
    "        skip_start = 10000\n",
    "        \n",
    "        fig3.data[0].x = np.arange(len(average_integrated_energy) - skip_start)\n",
    "        fig3.data[0].y = average_integrated_energy[skip_start:]\n",
    "    finally:\n",
    "        run_button.disabled = False\n",
    "        traces_button.disabled = False\n",
    "        run_button.style.button_color = 'green'\n",
    "        run_button.description=\"Run Monte-Carlo\"\n",
    "        \n",
    "run_button.on_click(interactive_plot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2efe9428",
   "metadata": {},
   "outputs": [],
   "source": [
    "def toggle_show_traces(button):\n",
    "    \"\"\"Callback function for the show traces button.\"\"\"\n",
    "    fig2.data[0].visible ^= True\n",
    "    if fig2.data[0].visible:\n",
    "        traces_button.description = \"Hide traces\"\n",
    "    else:\n",
    "        traces_button.description = \"Show traces\"\n",
    "    \n",
    "traces_button.on_click(toggle_show_traces)\n",
    "\n",
    "\n",
    "def show_current_point(trace, points, selector):\n",
    "    for i in points.point_inds:\n",
    "        fig1.data[1].x = [all_x1[i]]\n",
    "        fig1.data[1].y = [all_x2[i]]\n",
    "        fig1.data[1].z = [get_energy(all_x1[i], all_x2[i])]\n",
    "        \n",
    "\n",
    "        fig3.data[0].on_click(show_current_point)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e75a09c",
   "metadata": {},
   "source": [
    "<hr style=\"height:1px;border:none;color:#cccccc;background-color:#cccccc;\" />\n",
    "\n",
    "## Legend\n",
    "(How to use the interactive visualization)\n",
    "\n",
    "### Potential energy surface\n",
    "The left panel shows the potential energy surface (PES).\n",
    "\n",
    "You can rotate the surface by clicking and dragging.\n",
    "The `z` value indicates the value of the potential energy.\n",
    "\n",
    "When the mouse is on the potential energy surface, you will also see a black isoline, which you can use to find the global minimum.\n",
    "\n",
    "### Controls\n",
    "All the control widgets are on the bottom left.\n",
    "\n",
    "You can choose the initial `x` and `y` coordinates of the simulation.\n",
    "You can choose between two simple potentials: a single parabolic well, and a double well.\n",
    "The second can help understanding how the algorithm can overcome energy barriers.\n",
    "\n",
    "You can also tune a number of system and numerical parameters: the simulation temperature,\n",
    "the maximum move size, and the number of the simulation steps.\n",
    "\n",
    "### Running the simulation and resulting plots\n",
    "You can start the simulation by clicking the `Run Monte Carlo` button.\n",
    "This might take a few seconds before it completes. All plots will then be updated.\n",
    "\n",
    "The histogram of the canonical distribution as sampled by the simulation will be shown as a heatmap on the top right.\n",
    "By clicking the `Show traces` button, it will show 300 (evenly distributed) points on the heatmap.\n",
    "\n",
    "The total energy as a function of the step number is shown on the bottom right panel.\n",
    "10000 steps from the beginning of the simulation are skipped for the total energy figure as they are considered as an initial thermalization; note, however, that if you start in a point with very low probability (i.e., far away from the minimum), the simulation might require many more steps to thermalize."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68a20ed3",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}