notebook/band-theory/density_of_states.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# **The Electronic Density of States (DOS)**\n",
    "\n",
    "**Authors:** Dou Du, Taylor James Baird and Giovanni Pizzi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<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/band-theory/density_of_states.ipynb\n",
    "\n",
    "This notebook demonstrates various approaches for the numerical calculation of the density of states (DOS) for a 3D free-electron model with periodic boundary conditions.\n",
    "\n",
    "<hr style=\"height:1px;border:none;color:#cccccc;background-color:#cccccc;\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **Goals**\n",
    "   \n",
    "* Familiarize yourself with various numerical methods employed to calculate the electronic density of states.\n",
    "* Examine the resulting DOS and compare the accuracy and computational cost of the various methods."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **Background theory**\n",
    "\n",
    "[More on the background theory.](./theory/theory_density_of_states.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **Tasks and exercises**\n",
    "\n",
    "1. Investigate the influence of the number of k-points on the resulting DOS.\n",
    "    <details>\n",
    "    <summary style=\"color: red\">Solution</summary>\n",
    "    In the right panel, the plotted blue line is the analytical result for the DOS of the free electron model. \n",
    "    By choosing different numbers of k-points via the \"Number of k-points slider\", we can investigate how \n",
    "    the quality of the calculated results varies with the density of the k-point mesh. You will observe that the numerical results converge to the analytical result with increasing number of k-points. This can be attributed to the fact that the DOS can be interpreted as a probability density of electronic states as a function of energy. Since energy is generally related to the k-vector magnitude, the quality with which we resolve the range of energy eigenvalues is in turn directly controlled by how fine our sampling of the k-point mesh is.\n",
    "    </details>   \n",
    "\n",
    "2. Which method gives most accurate results? Which method is fastest and why?\n",
    "\n",
    "    <details>\n",
    "    <summary style=\"color: red\">Solution</summary>\n",
    "    The linear tetrahedra interpolation (LTI) method is an accurate numerical approach which linearly interpolates the 3D k-points grid within a set of tetrahedra into which reciprocal space has been subdivided (see the <a href=\"./theory/theory_density_of_states.ipynb\">background theory notebook</a>). The LTI method can yield much better results compared to a simple histogram. Gaussian smearing renders the \n",
    "    histogram plot considerably smoother, and closer to the analytical \n",
    "    solution. The histogram method is a simple statistical representation of the eigenvalues (one simply represents the frequency with which a given range of eigenenergies occur). This latter approach is typically the fastest of those discussed, but shall also generally give results with the poorest level of resolution.\n",
    "    </details>\n",
    "\n",
    "3. Set the number of k-points to the maximum value and consider Gaussian smearing with a reasonable value for the smearing parameter (say $\\sigma=0.07$). Looking at the calculated DOS with the G-vector range set to $0$, how does it compare with the analytical result? Can you explain any discrepancies you notice? Increasing the G-vector range to $1$, how does the calculated DOS now compare?\n",
    "\n",
    "    <details>\n",
    "    <summary style=\"color: red\">Solution</summary>\n",
    "    You will notice that even with a high number of k-points and an appropriate smearing parameter, the calculated DOS starts to deviate significantly from the analytical result above an     energy of around . This discrepancy can be attributed to the fact that when we only consider the G-vector $G=0$, we are neglecting to include several bands in the electronic energy       dispersion that still contribute energy eigenvalues that fall within the range we are interested in. That is, and you can verify the claim by increasing the G-vector slider to the value 1, several bands are folded back into the first Brillouin zone. You can see which bands undergo folding by tracking those which change from black to red upon incrementing the G-     vector slider.\n",
    "    </details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"height:1px;border:none;color:#cccccc;background-color:#cccccc;\" />\n",
    "\n",
    "## Interactive visualization\n",
    "(be patient, it might take a few seconds to load)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import seekpath\n",
    "import re\n",
    "import os\n",
    "import matplotlib\n",
    "from ase.dft.dos import linear_tetrahedron_integration as lti\n",
    "from ase.dft.kpoints import monkhorst_pack\n",
    "from ase.cell import Cell\n",
    "from scipy.stats import norm\n",
    "import plotly.graph_objects as go\n",
    "import plotly.express as px\n",
    "import time\n",
    "import matplotlib.pyplot as plt\n",
    "from ipywidgets import Button, RadioButtons, Layout, IntSlider, HBox, Text, VBox, Checkbox, Label, FloatSlider, Output, HTML\n",
    "from datetime import datetime\n",
    "\n",
    "%matplotlib widget"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def _compute_dos(kpts, G, ranges):\n",
    "    \"\"\"initial all the engienvalue according to the kpoints\n",
    "    \n",
    "    Args: \n",
    "        kpts: a array of kpts (kx, ky, kz)\n",
    "        G: the reciprocal lattice vectors (3x3)\n",
    "        ranges: the range of the reciprocal lattice \n",
    "        \n",
    "    Returns:\n",
    "        The eigenvalues of the free electron model.\n",
    "    \"\"\"\n",
    "    eigs = []\n",
    "    n = ranges\n",
    "    \n",
    "    for i in range(-n, n+1):\n",
    "        for j in range(-n, n+1):\n",
    "            for k in range(-n, n+1):\n",
    "                g_vector = i*G[0] + j*G[1] + k*G[2]\n",
    "                eigs.append(np.sum(0.5*(kpts + g_vector)**2, axis=3))\n",
    "\n",
    "    eigs = np.moveaxis(eigs, 0, -1)\n",
    "    return eigs\n",
    "   \n",
    "def _compute_total_kpts(G, grange=0):\n",
    "    \"\"\"Get all the kpoints \n",
    "    \n",
    "    Args:\n",
    "        G: the reciprocal lattice vectors (3x3)\n",
    "        grange: the range of the reciprocal lattice\n",
    "        \n",
    "    Returns:\n",
    "        The kpoints (kx, ky, kz) as a array\n",
    "    \"\"\"\n",
    "    tot_kpts = []\n",
    "    n = grange\n",
    "    \n",
    "    shape = (nkpt.value, nkpt.value, nkpt.value)\n",
    "    kpts = np.dot(monkhorst_pack(shape), G).reshape(shape + (3,))\n",
    "    kpts = kpts.reshape(nkpt.value**3, 3)\n",
    "\n",
    "    for i in range(-n, n+1):\n",
    "        for j in range(-n, n+1):\n",
    "            for k in range(-n, n+1):\n",
    "                g_vector = i*G[0] + j*G[1] + k*G[2]\n",
    "                tot_kpts.extend(kpts+g_vector)\n",
    "    return np.array(tot_kpts)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "alat_bohr = 7.72\n",
    "emax=1.5 # consider DOS up to 1.5 eV\n",
    "\n",
    "#Choose the cubic lattice for using the linear tetrahadron method (ASE)\n",
    "real_lattice_bohr = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]) * alat_bohr / 2.0;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# modifications\n",
    "\n",
    "style = {'description_width': 'initial'}\n",
    "\n",
    "output_bs = Output()\n",
    "\n",
    "line1 = Label(\n",
    "    value=r'Number of k-points per dimension, $N_k$'\n",
    ")\n",
    "\n",
    "line2 = Label(\n",
    "    value=r'(total number $=N_k^3$)'\n",
    ")\n",
    "nkpt = IntSlider(value=10, min=5, max=40, description=\"\", style={'description_width': 'initial'}, continuous_update=False)\n",
    "\n",
    "grange = IntSlider(value=1, min=0, max=1, description=\"G-vector range:\", style=style)\n",
    "grange_hint =HTML(value=f\"<b>Note that there may be a delay while the DOS is computed for G-vectors > 0.</b>\")\n",
    "gcov = FloatSlider(value=0.5, min=0.1, max=1.0, description=\"Guassian covariance:\", style=style)\n",
    "\n",
    "def prettify(label):\n",
    "    \"\"\"\n",
    "    Prettifier for matplotlib, using LaTeX syntax\n",
    "    :param label: a string to prettify\n",
    "    \"\"\"\n",
    "\n",
    "    label = (\n",
    "        label\n",
    "            .replace('GAMMA', r'$\\Gamma$')\n",
    "            .replace('DELTA', r'$\\Delta$')\n",
    "            .replace('LAMBDA', r'$\\Lambda$')\n",
    "            .replace('SIGMA', r'$\\Sigma$')\n",
    "    )\n",
    "    label = re.sub(r'_(.?)', r'$_{\\1}$', label)\n",
    "\n",
    "    return label\n",
    "def _get_band_energies(kpoints_list, b1, b2, b3, g_vectors_range):\n",
    "    energy_data_curves = np.zeros(((2*g_vectors_range+1)**3, len(kpoints_list)), dtype=np.float_)\n",
    "\n",
    "    cnt = 0\n",
    "    for g_i in range(-g_vectors_range,g_vectors_range+1):\n",
    "        for g_j in range(-g_vectors_range,g_vectors_range+1):\n",
    "            for g_k in range(-g_vectors_range,g_vectors_range+1):\n",
    "                g_vector = b1 * g_i + b2*g_j + b3 * g_k\n",
    "                energy_data_curves[cnt] = np.sum(0.5*(kpoints_list + g_vector)**2, axis=1)# This is k^2 - NOTE: units to be double checked!\n",
    "                cnt += 1\n",
    "\n",
    "\n",
    "    # bands are ordered as follows: first band, second band, ...\n",
    "    return energy_data_curves\n",
    "\n",
    "def get_bands(real_lattice_bohr, reference_distance = 0.025, g_vectors_range = 3):\n",
    "\n",
    "    # Simple way to get automatically the band path:\n",
    "    # I go back to real space, just put a single atom at the origin,\n",
    "    # then go back with seekpath.\n",
    "    # NOTE! This might not give the most general path, as e.g. there are two\n",
    "    # options for cubic FCC (cF1 and cF2 in seekpath).\n",
    "    # But this should be general enough for this tool.\n",
    "\n",
    "    structure = (real_lattice_bohr, [[0., 0., 0.]], [1])\n",
    "    # Use a H atom at the origin\n",
    "    seekpath_path = seekpath.get_explicit_k_path(structure, reference_distance=reference_distance)\n",
    "    b1, b2, b3 = np.array(seekpath_path['reciprocal_primitive_lattice'])\n",
    "\n",
    "    all_kpoints_x = np.array(seekpath_path['explicit_kpoints_linearcoord'])\n",
    "    all_kpoints_list = np.array(seekpath_path['explicit_kpoints_abs'])\n",
    "\n",
    "    segments_data = []\n",
    "    for segment_indices in seekpath_path['explicit_segments']:\n",
    "        start_label = seekpath_path['explicit_kpoints_labels'][segment_indices[0]]\n",
    "        end_label = seekpath_path['explicit_kpoints_labels'][segment_indices[1]-1]\n",
    "\n",
    "        kpoints_x = all_kpoints_x[slice(*segment_indices)]\n",
    "        kpoints_list = all_kpoints_list[slice(*segment_indices)]\n",
    "\n",
    "        energy_bands = _get_band_energies(kpoints_list, b1, b2, b3, g_vectors_range)\n",
    "\n",
    "        segments_data.append({\n",
    "            'start_label': start_label,\n",
    "            'end_label': end_label,\n",
    "            'kpoints_list': kpoints_list,\n",
    "            'kpoints_x': kpoints_x,\n",
    "            'energy_bands': energy_bands,\n",
    "            'b1': b1,\n",
    "            'b2': b2,\n",
    "            'b3': b3,\n",
    "        })\n",
    "\n",
    "    return segments_data\n",
    "\n",
    "\n",
    "def plot_bandstructure(c):\n",
    "    global G, segments_data, lbands\n",
    "    \n",
    "    segments_data = get_bands(real_lattice_bohr)\n",
    "    G = np.array([segments_data[0]['b1'], segments_data[0]['b2'], segments_data[0]['b3']])\n",
    "    \n",
    "    x_ticks = []\n",
    "    x_labels = []\n",
    "    lbands = []\n",
    "\n",
    "    for segment_data in segments_data:\n",
    "        if not x_labels:\n",
    "            x_labels.append(prettify(segment_data['start_label']))\n",
    "            x_ticks.append(segment_data['kpoints_x'][0])\n",
    "        else:\n",
    "            if x_labels[-1] != prettify(segment_data['start_label']):\n",
    "                x_labels[-1] += \"|\" + prettify(segment_data['start_label'])\n",
    "        x_labels.append(prettify(segment_data['end_label']))\n",
    "        x_ticks.append(segment_data['kpoints_x'][-1])\n",
    "\n",
    "        for energy_band in segment_data['energy_bands']:\n",
    "            line, = ax_bs.plot(segment_data['kpoints_x'], energy_band, 'k')\n",
    "            lbands.append(line)\n",
    "\n",
    "    ax_bs.set_ylim([0, 1.0])\n",
    "    ax_bs.yaxis.tick_right()\n",
    "    ax_bs.yaxis.set_label_position(\"right\")\n",
    "    ax_bs.set_ylabel('Free-electron energy (eV)')\n",
    "    ax_bs.set_xlim([np.min(x_ticks), np.max(x_ticks)])\n",
    "    ax_bs.set_xticks(x_ticks)\n",
    "    ax_bs.set_xticklabels(x_labels)\n",
    "    ax_bs.grid(axis='x', color='red', linestyle='-', linewidth=0.5)\n",
    "    fig.tight_layout()\n",
    "    \n",
    "    update_bands_color('bands')\n",
    "    \n",
    "def update_bands_color(c):\n",
    "    n = 3\n",
    "    \n",
    "    shape = (nkpt.value, nkpt.value, nkpt.value)\n",
    "    kpts = np.dot(monkhorst_pack(shape), G).reshape(shape + (3,))\n",
    "    eigs = _compute_dos(kpts, G, grange.value)\n",
    "\n",
    "    index = 0\n",
    "    \n",
    "    for segment_data in segments_data:\n",
    "        for i in range(-n, n+1):\n",
    "            for j in range(-n, n+1):\n",
    "                for k in range(-n, n+1): \n",
    "                    if abs(i) <= grange.value and abs(j) <= grange.value and abs(k) <=grange.value:\n",
    "                        lbands[index].set_color('r')\n",
    "                    else:\n",
    "                        lbands[index].set_color('k')\n",
    "                    index+=1\n",
    "\n",
    "grange.observe(update_bands_color, names=\"value\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "style = {'description_width': 'initial'}\n",
    "\n",
    "line1 = Label(\n",
    "    value=r'Number of kpoints per dimension, $N_k$'\n",
    ")\n",
    "\n",
    "line2 = Label(\n",
    "    value=r'(total number $=N_k^3$)'\n",
    ")\n",
    "nkpt_box = HBox([VBox([line1, line2]), nkpt])\n",
    "nbin = IntSlider(value=50, min=20, max=100, description=\"Number of bins:\", layout=Layout(width=\"400px\"), style={'description_width': 'initial'})\n",
    "gstd = FloatSlider(value=0.01, min=0.01, max=0.1, step=0.01, description=\"Gaussian $\\sigma$ (eV):\", layout=Layout(width=\"300px\"), style={'description_width': 'initial'})\n",
    "\n",
    "#All buttons\n",
    "btlti = Button(description=\"Add plot\", style = {'button_color':'green'})\n",
    "bthist = Button(description=\"Add plot\", style = {'button_color':'green'})\n",
    "btgas = Button(description=\"Add plot\", style = {'button_color':'green'})\n",
    "btclear = Button(description=\"Clear plot\", style = {'button_color':'green'})\n",
    "\n",
    "#Ouput for the DOS figure\n",
    "output = Output()\n",
    "\n",
    "def compute_dos_lti(c):\n",
    "    \"\"\"Compute the DOS uing the ASE linear tetrahedron interpolation method.\n",
    "    \"\"\"\n",
    "    global llti\n",
    "    \n",
    "    btlti.disabled = True\n",
    "    bthist.disabled = True\n",
    "    btgas.disabled = True\n",
    "    btclear.disabled = True\n",
    "    btlti.style = {'button_color':'red'}\n",
    "    btlti.description = \"Running\"\n",
    "    \n",
    "    \n",
    "    try:\n",
    "        llti.remove()\n",
    "    except:\n",
    "        pass\n",
    "    \n",
    "    shape = (nkpt.value, nkpt.value, nkpt.value)\n",
    "    G = Cell(real_lattice_bohr).reciprocal()*2*np.pi\n",
    "    kpts = np.dot(monkhorst_pack(shape), G).reshape(shape + (3,))\n",
    "\n",
    "    eigs = _compute_dos(kpts, G, grange.value)\n",
    "\n",
    "    dosx = np.linspace(0, emax, 500)\n",
    "    dosy = lti(real_lattice_bohr, eigs, dosx)\n",
    "\n",
    "    # normalize lti dos for comparison with analytical result\n",
    "\n",
    "    norm_lti=(dosx[1]-dosx[0])*np.sum(dosy)\n",
    "    \n",
    "    vol=(alat_bohr/2.)**3\n",
    "    norm_an=(2./3.)*(vol/np.pi**2)*2**(1./2)*emax**(3./2) # analytical value of integrated DOS\n",
    "\n",
    "    dosy/=norm_lti\n",
    "    dosy*=norm_an\n",
    "    \n",
    "    \n",
    "    llti, = ax_dos.plot(dosy, dosx, 'r-', label='LTI')\n",
    "    ax_dos.legend(loc=1, bbox_to_anchor=(1.3, 1.0))\n",
    "    \n",
    "    btlti.disabled = False\n",
    "    bthist.disabled = False\n",
    "    btgas.disabled = False\n",
    "    btclear.disabled = False\n",
    "    btlti.style = {'button_color':'green'}\n",
    "    btlti.description=\"Add plot\"\n",
    "\n",
    "btlti.on_click(compute_dos_lti)\n",
    "\n",
    "def compute_dos_histogram(c):\n",
    "    \"\"\"Compute the DOS as a histogram.\n",
    "    \"\"\"\n",
    "    global lhist\n",
    "    \n",
    "    btlti.disabled = True\n",
    "    bthist.disabled = True\n",
    "    btgas.disabled = True\n",
    "    btclear.disabled = True\n",
    "    bthist.style = {'button_color':'green'}\n",
    "        \n",
    "    try:\n",
    "        lhist.remove()\n",
    "    except:\n",
    "        pass\n",
    "    \n",
    "    shape = (nkpt.value, nkpt.value, nkpt.value)\n",
    "    G = Cell(real_lattice_bohr).reciprocal()*2*np.pi\n",
    "    kpts = np.dot(monkhorst_pack(shape), G).reshape(shape + (3,))\n",
    "    eigs = _compute_dos(kpts, G, grange.value)\n",
    "    num_bins=nbin.value\n",
    "\n",
    "    vol=(alat_bohr/2.)**3\n",
    "\n",
    "    analy_x = np.linspace(0, emax, 500);\n",
    "    analy_y = (1.0/(2.0*np.pi**2))*(2.0)**(3./2.)*analy_x**0.5*(alat_bohr / 2.0)**3.0;\n",
    "\n",
    "    norm=np.sum(analy_y*(analy_x[1]-analy_x[0]))\n",
    "    norm_an=(2./3.)*(vol/np.pi**2)*2**(1./2)*emax**(3./2) # analytical value of integrated DOS\n",
    "  \n",
    "\n",
    "    hy, hx = np.histogram(eigs.ravel(),range=[0,emax],bins=num_bins )\n",
    "   \n",
    "    hy = hy/np.sum(hy*np.diff(hx)) # normalize histogram to make total area 1\n",
    "    hy*=norm_an # renormalize histogram for comparison with analytical DOS (norm gives integrated DOS of analytical curve)\n",
    "  \n",
    "\n",
    "    lhist = ax_dos.barh(hx[:-1]+np.diff(hx)[0], hy, color='yellow', edgecolor='black', \n",
    "                       height=np.diff(hx), label=\"Histogram\")\n",
    "    ax_dos.legend(loc=1, bbox_to_anchor=(1.3, 1.0))\n",
    "\n",
    "    btlti.disabled = False\n",
    "    bthist.disabled = False\n",
    "    btgas.disabled = False\n",
    "    btclear.disabled = False\n",
    "    bthist.style = {'button_color':'green'}\n",
    "    \n",
    "bthist.on_click(compute_dos_histogram)\n",
    "\n",
    "def compute_dos_gaussian(c):\n",
    "    \"\"\"Computing the DOS using Gaussian smearing method.\n",
    "    \"\"\"\n",
    "    global lgas\n",
    "    \n",
    "    btlti.disabled = True\n",
    "    bthist.disabled = True\n",
    "    btgas.disabled = True\n",
    "    btclear.disabled = True\n",
    "    btgas.style = {'button_color':'red'}\n",
    "    btgas.description = \"Running\"\n",
    "        \n",
    "    try:\n",
    "        lgas.remove()\n",
    "    except:\n",
    "        pass\n",
    "    \n",
    "    shape = (nkpt.value, nkpt.value, nkpt.value)\n",
    "    G = Cell(real_lattice_bohr).reciprocal()*2*np.pi\n",
    "    kpts = np.dot(monkhorst_pack(shape), G).reshape(shape + (3,))\n",
    "    eigs = _compute_dos(kpts, G, grange.value)\n",
    "    num_bins=nbin.value\n",
    "    \n",
    "    knum=nkpt.value\n",
    "    Gnum=grange.value\n",
    "    \n",
    "    # Accelerated approach to smearing: we smear out a histogram of evals instead of deltas centered on each eval\n",
    "    # instead of summing smeared gaussians centered on eigs, histogram and then smear\n",
    "\n",
    "    vol=(alat_bohr/2.)**3\n",
    "    norm_an=(2./3.)*(vol/np.pi**2)*2**(1./2)*emax**(3./2) # analytical value of integrated DOS\n",
    "\n",
    "    \n",
    "    hy, hx = np.histogram(eigs.ravel(),range=[0,emax],bins=num_bins ) # further right bin edge included in hx\n",
    "    hx=hx[:-1]\n",
    "    gx=hx\n",
    "    gy=0.0*hx\n",
    "\n",
    "    for i,Ei in enumerate(hx):\n",
    "        gy += hy*norm(Ei, gstd.value).pdf(hx)\n",
    "\n",
    "    gy = gy/np.sum(gy*(hx[1]-hx[0])) # normalize histogram to make total area 1\n",
    "    gy*=norm_an # renormalize histogram for comparison with analytical DOS (norm gives integrated DOS of analytical curve)\n",
    "    \n",
    "      \n",
    "    \n",
    "    # Standard approach to Gaussian smearing\n",
    "#     for eig in eigs.ravel():\n",
    "#         gy += norm(eig, gstd.value).pdf(gx)\n",
    "\n",
    "#     norm_g=np.sum((gx[1]-gx[0])*gy)\n",
    "#     gy = norm_an*gy/norm_g\n",
    "\n",
    "\n",
    "    lgas, = ax_dos.plot(gy, gx, 'k--', label=\"Gaussian smearing\")\n",
    "    ax_dos.legend(loc=1, bbox_to_anchor=(1.3, 1.0))\n",
    "    \n",
    "    btlti.disabled = False\n",
    "    bthist.disabled = False\n",
    "    btgas.disabled = False\n",
    "    btclear.disabled = False\n",
    "    btgas.style = {'button_color':'green'}\n",
    "    btgas.description = \"Add plot\"\n",
    "    \n",
    "btgas.on_click(compute_dos_gaussian)\n",
    "\n",
    "\n",
    "def init_dos_plot():\n",
    "    \"\"\"Init the DOS plot.\n",
    "    \"\"\"\n",
    "    global hline, ann,ax_dos\n",
    "    btlti.disabled = True\n",
    "    bthist.disabled = True\n",
    "    btgas.disabled = True\n",
    "    \n",
    "    analy_x = np.linspace(0, emax, 500);\n",
    "    analy_y = (1.0/(2.0*np.pi**2))*(2.0)**(3./2.)*analy_x**0.5*(alat_bohr / 2.0)**3.0;\n",
    "    lanaly, = ax_dos.plot(analy_y,analy_x, 'b', label='Analytical solution')\n",
    "    \n",
    "    ax_dos.legend(loc=1, bbox_to_anchor=(1.3, 1.0))\n",
    "    ax_dos.yaxis.tick_right()\n",
    "    ax_dos.yaxis.set_label_position(\"right\")\n",
    "    ax_dos.set_xlabel('Density of States')\n",
    "\n",
    "    ax_dos.set_ylabel('Energy (eV)')\n",
    "    fig.tight_layout()\n",
    "    \n",
    "    btlti.disabled = False\n",
    "    bthist.disabled = False\n",
    "    btgas.disabled = False\n",
    "     \n",
    "    \n",
    "with output:\n",
    "    \"\"\"Set the figure for the DOS\n",
    "    \"\"\"\n",
    "    global fig, ax_dos, ax_bs\n",
    "    fig, (ax_bs, ax_dos) = plt.subplots(1, 2, sharey=True)\n",
    "    fig.set_size_inches(8.5, 5.0)\n",
    "    fig.canvas.header_visible = False\n",
    "    fig.canvas.layout.width = \"850px\"\n",
    "    fig.tight_layout()\n",
    "    init_dos_plot()\n",
    "    plot_bandstructure('bands')\n",
    "    plt.show()    \n",
    "    \n",
    "def clear_plot(c):\n",
    "    \"\"\"Clear the DOS calculated results when the \"Clear\" button is clicked.\n",
    "    \"\"\"\n",
    "    ax_dos.clear()\n",
    "    init_dos_plot()\n",
    "    \n",
    "btclear.on_click(clear_plot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = px.data.gapminder()\n",
    "\n",
    "X, Y, Z = np.mgrid[-2:2:40j, -2:2:40j, -2:2:40j]\n",
    "\n",
    "# Fermi surface\n",
    "values = 0.5*(X * X + Y * Y + Z * Z)\n",
    "\n",
    "G = Cell(real_lattice_bohr).reciprocal()*2*np.pi\n",
    "kpts = _compute_total_kpts(G)\n",
    "\n",
    "def update_kpts_fig(c):\n",
    "    \"\"\"Update the k-points plot when tuning the k-points slider.\n",
    "    \"\"\"\n",
    "    kpts = _compute_total_kpts(G)\n",
    "    \n",
    "\n",
    "nkpt.observe(update_kpts_fig, names=\"value\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "00b43cbb9859412e905577e11aa76e09",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(Output(), VBox(children=(HBox(children=(VBox(children=(Label(value='Number of kpoints per dimen…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Group buttons with descriptions as labels\n",
    "method1 = HBox([HBox([HTML(value=f\"<b>Simple histogram of the eigenvalues</b>\"), nbin]), bthist])\n",
    "method2 = HBox([HBox([HTML(value=f\"<b>Gaussian smearing method</b>\"), gstd]),btgas ])\n",
    "method3 = HBox([HTML(value=\"<b>Linear tetrahedron interpolation method</b>\"),btlti ])\n",
    "method4 = HBox([btclear, Label(value=\"(Clear the calculated results)\")])\n",
    "nkpt_box = HBox([VBox([line1, line2]), nkpt,method4])\n",
    "\n",
    "label1 = HTML(value = f\"<b><font color='red'>Choose a method to calculate the DOS:</b>\")\n",
    "display(VBox([output,VBox([ nkpt_box,VBox([grange,grange_hint]), label1, method1, method2, method3])]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"height:1px;border:none;color:#cccccc;background-color:#cccccc;\" />\n",
    "\n",
    "## **Legend**\n",
    "\n",
    "(How to use the interactive visualization)\n",
    "\n",
    "## **Interactive figures**\n",
    "\n",
    "The left panel shows the electronic bandstructure of a free electron gas in 3 dimensions while the right panel displays the corresponding density of states. You can choose the number of k-points used in the computation of the DOS by manipulating the k-points slider. Similarly, the number of G-vectors employed in the calculation can be varied with the \"G-vector range\" slider (the value is capped at 1 to avoid excessive computational time in the case of the LTI method, however one should still expect long waiting times in in this case when $N_k$ is large).\n",
    "\n",
    "## **Controls**\n",
    "\n",
    "Three buttons allow one to compute the DOS with the three methods discussed \n",
    "earlier. The calculated DOS will appear in the figure on the right, superimposed on the analytic curve.\n",
    "Computing the DOS with a large number of k-points may take several seconds.\n",
    "When using the simple histogramming method, you can adjust the number of bins employed by varying the \"Number of bins\" slider.\n",
    "For the Gaussian smearing method, you can also tune the standard \n",
    "deviation, $\\sigma$, of the Gaussian functions by adjusting the \"Gaussian $\\sigma$\" slider."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  },
  "voila": {
   "authors": "Dou Du and Giovanni Pizzi",
   "notification": "true"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}