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 in periodic boundary conditions.\n",
    "\n",
    "<hr style=\"height:1px;border:none;color:#cccccc;background-color:#cccccc;\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **Goals**\n",
    "   \n",
    "* Learn the methods to calculate the density of states.\n",
    "* Examine the resulting DOS and evaluate the accuracy of various methods to compute the DOS."
   ]
  },
  {
   "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 blue line is the analytical solution for the DOS. \n",
    "    By choosing the different number of kpoints from the slider, we can compare \n",
    "    the calculated results with the analytical solution. The density of states \n",
    "    is a probability distribution. The kpoints sampling is shown as red dots \n",
    "    in the left panel. The more kpoints the better results we can obtain. \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",
    "    Linear tetrahedra interpolation (LTI) is an accurate numerical method, \n",
    "    which interpolates the 3D kpoints grid. LTI method can give much better \n",
    "    results rather than a simple histogram. Gaussian smearing makes the \n",
    "    histogram plot much smoother, which is closer to the analytical \n",
    "    solution. The histogram method is a simple statistic of the eigenvalues, \n",
    "    which should be the fastest to compute.\n",
    "    </details>\n",
    "\n",
    "3. Why do the calculated results start to get diverge when the energy level is \n",
    "higher than a certain value? Could you explain it with the k-space plot?\n",
    "\n",
    "    <details>\n",
    "    <summary style=\"color: red\">Solution</summary>\n",
    "    In the free electron model, the energy isosurface is a sphere shown in \n",
    "    the left panel. The kpoints grid must be larger than the energy \n",
    "    isosurface to obtain the correct DOS at the energy level. Here, we \n",
    "    have a fixed length of the kpoints grid. When the energy is larger than \n",
    "    about 0.31, the kpoints grid cannot include the whole sphere (check it by \n",
    "    clicking on the right panel to move the isovalue above and below 0.31). \n",
    "    </details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"height:1px;border:none;color:#cccccc;background-color:#cccccc;\" />\n",
    "\n",
    "## Interactive visualization\n",
    "(be patient, it might take a few seconds to load)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "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, VBox, Checkbox, Label, FloatSlider, Output, HTML\n",
    "from datetime import datetime\n",
    "\n",
    "%matplotlib widget"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_kernel_id():\n",
    "    \"\"\"Get the current kernel ID, to distinguish different users.\n",
    "    \n",
    "    Call this only from within python jupyter notebooks.\n",
    "    \"\"\"\n",
    "    from IPython.lib import kernel\n",
    "    connection_file_path = kernel.get_connection_file()\n",
    "    connection_file = os.path.basename(connection_file_path)\n",
    "    kernel_id = connection_file.split('-', 1)[1].split('.')[0]\n",
    "    return kernel_id"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def log_button(name):\n",
    "    try:\n",
    "        allow_datacollection\n",
    "    except:\n",
    "        pass\n",
    "    else:\n",
    "        if allow_datacollection:\n",
    "            log_file = open('../log.dat', 'a+');\n",
    "            log_file.write(datetime.now().strftime(\"%Y-%m-%d %H:%M:%S\") + ' ')\n",
    "            log_file.write(get_kernel_id() + ' ')\n",
    "            log_file.write(name + ' ')\n",
    "            log_file.write(str(nkpt.value) + '\\n')\n",
    "            log_file.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alat_bohr = 7.72\n",
    "\n",
    "lattices = np.zeros((3, 3, 3));\n",
    "\n",
    "lattices[0] = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]) * alat_bohr / 2.0;\n",
    "lattices[1] = np.array([[0, 1, 1], [1, 0, 1], [1, 1, 0]]) * alat_bohr / 2.0;\n",
    "lattices[2] = np.array([[-1, 1, 1], [1, -1, 1], [1, 1, -1]]) * alat_bohr / 2.0;\n",
    "\n",
    "#Choose the cubic lattice for using the linear tetrahadron method (ASE)\n",
    "real_lattice_bohr = lattices[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "style = {'description_width': 'initial'}\n",
    "\n",
    "nkpt = IntSlider(value=4, min=4, max=25, description=\"Number of kpoints:\", style={'description_width': 'initial'}, continuous_update=False)\n",
    "nbin = IntSlider(value=30, min=30, max=500, description=\"Number of bins:\", layout=Layout(width=\"300px\"), 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=\"Tetrahedra\", style = {'button_color':'green'})\n",
    "bthist = Button(description=\"Histogram\", style = {'button_color':'green'})\n",
    "btgas = Button(description=\"Smearing\", 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",
    "    log_button('Tetrahedra')\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, 0)\n",
    "\n",
    "    dosx = np.linspace(0, 10, 500)\n",
    "    dosy = lti(real_lattice_bohr, eigs, dosx)\n",
    "    \n",
    "    llti, = ax.plot(dosy, dosx, 'r-', label='LTI')\n",
    "    ax.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=\"Tetrahedron\"\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",
    "    log_button('Histogram')\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, 0)\n",
    "    \n",
    "    hy, hx = np.histogram(eigs.ravel(), bins=nbin.value, range=(0.0, 3.0))\n",
    "    hy = hy/np.sum(hy*np.diff(hx))*np.shape(eigs)[-1]\n",
    "    \n",
    "    lhist = ax.barh(hx[:-1]+np.diff(hx)[0], hy, color='yellow', edgecolor='black', \n",
    "                       height=np.diff(hx), label=\"Histogram\")\n",
    "    ax.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",
    "    log_button('Smearing')\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, 0)\n",
    "    \n",
    "    gx = np.linspace(-0.03, 5, 500)\n",
    "    gy = 0.0*gx\n",
    "    \n",
    "    for eig in eigs.ravel():\n",
    "        gy += norm(eig, gstd.value).pdf(gx)\n",
    " \n",
    "    gy = gy/np.size(eigs)*np.shape(eigs)[-1]\n",
    "    lgas, = ax.plot(gy, gx, 'k--', label=\"Gaussian smearing\")\n",
    "    ax.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 = \"Smearing\"\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\n",
    "    btlti.disabled = True\n",
    "    bthist.disabled = True\n",
    "    btgas.disabled = True\n",
    "    \n",
    "    analy_x = np.linspace(0, 0.5, 500);\n",
    "    analy_y = 1.0/(2.0*np.pi**2)*2.0**0.5*analy_x**0.5*(alat_bohr / 2.0)**3.0;\n",
    "    lanaly, = ax.plot(analy_y, analy_x, 'b', label='Analytical solution')\n",
    "    \n",
    "    ax.set_ylim([-0.03, 0.5])\n",
    "    ax.set_xlim([0, analy_y.max() + 3.1])\n",
    "    ax.legend(loc=1, bbox_to_anchor=(1.3, 1.0))\n",
    "    ax.yaxis.tick_right()\n",
    "    ax.yaxis.set_label_position(\"right\")\n",
    "    ax.set_ylabel('Density of States (eV)')\n",
    "    hline = ax.axhline(0.3, color=\"red\")\n",
    "    ann = ax.annotate(r\"$\\frac{\\hbar^2k^2}{2m}$ isosurf. (click to move)\", xy=(0.2, 0.31), fontsize=10)\n",
    "    fig.tight_layout()\n",
    "    \n",
    "    btlti.disabled = False\n",
    "    bthist.disabled = False\n",
    "    btgas.disabled = False\n",
    "\n",
    "def onclick(event):\n",
    "    \"\"\"Click to move the isovalue line (red horizontal line) and update the kpoints plot.\n",
    "    \"\"\"\n",
    "    hline.set_ydata(event.ydata)\n",
    "    figkpts.data[0].isomin = event.ydata\n",
    "    figkpts.data[0].isomax = event.ydata\n",
    "    ann.set_position((0.5, event.ydata + 0.01))\n",
    "        \n",
    "with output:\n",
    "    \"\"\"Set the figure for the DOS\n",
    "    \"\"\"\n",
    "    global fig, ax\n",
    "    fig, ax = plt.subplots()\n",
    "    fig.set_size_inches(3.2, 5.0)\n",
    "    fig.canvas.header_visible = False\n",
    "    fig.canvas.layout.width = \"380px\"\n",
    "    fig.tight_layout()\n",
    "    init_dos_plot()\n",
    "    cid = fig.canvas.mpl_connect('button_press_event', onclick)\n",
    "    plt.show()\n",
    "    \n",
    "def clear_plot(c):\n",
    "    \"\"\"Clear the DOS calculated results when the \"Clear\" button is clicked.\n",
    "    \"\"\"\n",
    "    ax.clear()\n",
    "    init_dos_plot()\n",
    "    \n",
    "btclear.on_click(clear_plot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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",
    "#Init the kpoints plot with the plotly package\n",
    "figkpts = go.FigureWidget(data=[go.Isosurface(\n",
    "    x=X.flatten(),\n",
    "    y=Y.flatten(),\n",
    "    z=Z.flatten(),\n",
    "    value=values.flatten(),\n",
    "    opacity=0.5,\n",
    "    isomin=0.3,\n",
    "    isomax=0.3,\n",
    "    surface_count=1,\n",
    "    caps=dict(x_show=False, y_show=False)),\n",
    "    go.Scatter3d(x=kpts[:,0], y=kpts[:,1], z=kpts[:,2], mode='markers',\n",
    "                marker=dict(size=1.5, color='red'))],\n",
    "    layout=go.Layout(width=450, title='Kpoints (red dots) in reciprocal space and'\n",
    "                     +'<br>energy isosurface (isovalue can be set by <br> clicking on the left figure)',\n",
    "                     scene=dict(bgcolor = 'rgb(20, 24, 54)',\n",
    "                                xaxis = dict(title=r'kx', titlefont_color='white'),\n",
    "                                yaxis = dict(title=r'ky', titlefont_color='white'),\n",
    "                                zaxis = dict(title=r'kz', titlefont_color='white'))))\n",
    "\n",
    "\n",
    "def update_kpts_fig(c):\n",
    "    \"\"\"Update the kpoints plot when tuning the kpoints slider.\n",
    "    \"\"\"\n",
    "    kpts = _compute_total_kpts(G)\n",
    "    \n",
    "    with figkpts.batch_update():\n",
    "        figkpts.data[1].x = kpts[:, 0]\n",
    "        figkpts.data[1].y = kpts[:, 1]\n",
    "        figkpts.data[1].z = kpts[:, 2]\n",
    "        \n",
    "        if nkpt.value >= 8:\n",
    "            figkpts.data[1].marker['size'] = 1.0\n",
    "        else:\n",
    "            figkpts.data[1].marker['size'] = 1.5\n",
    "\n",
    "def half_sphere():\n",
    "    \"\"\"Only show half of the isosurface.\n",
    "    \"\"\"\n",
    "    X, Y, Z = np.mgrid[-6:6:40j, 0:6:40j, -6:6:40j]\n",
    "    values = 0.5*(X * X + Y * Y + Z * Z)\n",
    "    figkpts.data[0].x = X.flatten()\n",
    "    figkpts.data[0].y = Y.flatten()\n",
    "    figkpts.data[0].z = Z.flatten()\n",
    "    figkpts.data[0].value = values.flatten()\n",
    "    \n",
    "\n",
    "nkpt.observe(update_kpts_fig, names=\"value\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Group buttons with descriptions as labels\n",
    "method1 = VBox([HBox([bthist, nbin]), Label(value=\"(Simple histogram of the eigenvalues)\")])\n",
    "method2 = VBox([HBox([btgas, gstd]), Label(value=\"(Gaussian smearing method)\")])\n",
    "method3 = HBox([btlti, Label(value=\"(Linear tetrahedron interpolation method)\")])\n",
    "method4 = HBox([btclear, Label(value=\"(Clear the calculated results)\")])\n",
    "\n",
    "label1 = HTML(value = f\"<b><font color='red'>Choose a method to calculate the DOS:</b>\")\n",
    "\n",
    "display(HBox([VBox([figkpts, nkpt, label1, method1, method2, method3]), VBox([output, method4])]))"
   ]
  },
  {
   "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 kpoints (red dots) in the reciprocal space (k-space).\n",
    "A transparent sphere is employed to represent the energy isosurface \n",
    "(Fermi surface). The value of the isosurface can be set by clicking on the \n",
    "right panel (red horizontal line). Choose the number of kpoints in each \n",
    "dimension from the kpoints slider. The left panel will update accordingly \n",
    "when the number of the kpoints is changed.\n",
    "\n",
    "## **Controls**\n",
    "\n",
    "Three buttons allow to compute the DOS with the three methods discussed \n",
    "earlier. The DOS results will appear in the figure on the right.\n",
    "Get results with many kpoints might take several seconds.\n",
    "For the Gaussian smearing method, you can also tune the standard \n",
    "deviation $\\sigma$ of Gaussian functions with the slider next to \n",
    "the \"Smearing\" button."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  },
  "voila": {
   "authors": "Dou Du and Giovanni Pizzi",
   "notification": "true"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}