diff --git a/notebook/band-theory/density_of_states.ipynb b/notebook/band-theory/density_of_states.ipynb index 2b9c917..7f855c0 100644 --- a/notebook/band-theory/density_of_states.ipynb +++ b/notebook/band-theory/density_of_states.ipynb @@ -17,7 +17,7 @@ "\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", + "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", "
" ] @@ -28,8 +28,8 @@ "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." + "* 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." ] }, { @@ -50,36 +50,25 @@ "1. Investigate the influence of the number of k-points on the resulting DOS.\n", "
\n", " Solution\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", + " 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", "
\n", "\n", "2. Which method gives most accurate results? Which method is fastest and why?\n", "\n", "
\n", " Solution\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", + " 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 background theory notebook). 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", "
\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", + "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", "
\n", " Solution\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", + " 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", "
" ] }, @@ -95,7 +84,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 16, "metadata": {}, "outputs": [], "source": [ @@ -112,7 +101,7 @@ "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 ipywidgets import Button, RadioButtons, Layout, IntSlider, HBox, Text, VBox, Checkbox, Label, FloatSlider, Output, HTML\n", "from datetime import datetime\n", "\n", "%matplotlib widget" @@ -120,46 +109,7 @@ }, { "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, + "execution_count": 17, "metadata": {}, "outputs": [], "source": [ @@ -214,25 +164,182 @@ }, { "cell_type": "code", - "execution_count": null, + "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", - "lattices = np.zeros((3, 3, 3));\n", + "style = {'description_width': 'initial'}\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", + "output_bs = Output()\n", "\n", - "#Choose the cubic lattice for using the linear tetrahadron method (ASE)\n", - "real_lattice_bohr = lattices[0]" + "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\"Note that there may be a delay while the DOS is computed for G-vectors > 0.\")\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": [] }, @@ -240,14 +347,21 @@ "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", + "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=\"Tetrahedra\", style = {'button_color':'green'})\n", - "bthist = Button(description=\"Histogram\", style = {'button_color':'green'})\n", - "btgas = Button(description=\"Smearing\", style = {'button_color':'green'})\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", @@ -265,7 +379,6 @@ " btlti.style = {'button_color':'red'}\n", " btlti.description = \"Running\"\n", " \n", - " log_button('Tetrahedra')\n", " \n", " try:\n", " llti.remove()\n", @@ -276,20 +389,31 @@ " 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", + " eigs = _compute_dos(kpts, G, grange.value)\n", "\n", - " dosx = np.linspace(0, 10, 500)\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", - " llti, = ax.plot(dosy, dosx, 'r-', label='LTI')\n", - " ax.legend(loc=1, bbox_to_anchor=(1.3, 1.0))\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=\"Tetrahedron\"\n", + " btlti.description=\"Add plot\"\n", "\n", "btlti.on_click(compute_dos_lti)\n", "\n", @@ -303,9 +427,7 @@ " btgas.disabled = True\n", " btclear.disabled = True\n", " bthist.style = {'button_color':'green'}\n", - " \n", - " log_button('Histogram')\n", - " \n", + " \n", " try:\n", " lhist.remove()\n", " except:\n", @@ -314,14 +436,27 @@ " 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", + " 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.legend(loc=1, bbox_to_anchor=(1.3, 1.0))\n", + " ax_dos.legend(loc=1, bbox_to_anchor=(1.3, 1.0))\n", "\n", " btlti.disabled = False\n", " bthist.disabled = False\n", @@ -342,9 +477,7 @@ " btclear.disabled = True\n", " btgas.style = {'button_color':'red'}\n", " btgas.description = \"Running\"\n", - " \n", - " log_button('Smearing')\n", - " \n", + " \n", " try:\n", " lgas.remove()\n", " except:\n", @@ -353,24 +486,49 @@ " 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", + " eigs = _compute_dos(kpts, G, grange.value)\n", + " num_bins=nbin.value\n", + " \n", + " knum=nkpt.value\n", + " Gnum=grange.value\n", " \n", - " gx = np.linspace(-0.03, 5, 500)\n", - " gy = 0.0*gx\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", - " 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", + " 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 = \"Smearing\"\n", + " btgas.description = \"Add plot\"\n", " \n", "btgas.on_click(compute_dos_gaussian)\n", "\n", @@ -378,54 +536,45 @@ "def init_dos_plot():\n", " \"\"\"Init the DOS plot.\n", " \"\"\"\n", - " global hline, ann\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, 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", + " 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", - "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", + " \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", + " 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 = \"380px\"\n", + " fig.canvas.layout.width = \"850px\"\n", " fig.tight_layout()\n", " init_dos_plot()\n", - " cid = fig.canvas.mpl_connect('button_press_event', onclick)\n", - " plt.show()\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.clear()\n", + " ax_dos.clear()\n", " init_dos_plot()\n", " \n", "btclear.on_click(clear_plot)" @@ -433,7 +582,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 21, "metadata": {}, "outputs": [], "source": [ @@ -447,71 +596,45 @@ "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", - " +'
energy isosurface (isovalue can be set by
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", + " \"\"\"Update the k-points plot when tuning the k-points 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, + "execution_count": 22, "metadata": {}, - "outputs": [], + "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 = 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", + "method1 = HBox([HBox([HTML(value=f\"Simple histogram of the eigenvalues\"), nbin]), bthist])\n", + "method2 = HBox([HBox([HTML(value=f\"Gaussian smearing method\"), gstd]),btgas ])\n", + "method3 = HBox([HTML(value=\"Linear tetrahedron interpolation method\"),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\"Choose a method to calculate the DOS:\")\n", - "\n", - "display(HBox([VBox([figkpts, nkpt, label1, method1, method2, method3]), VBox([output, method4])]))" + "display(VBox([output,VBox([ nkpt_box,VBox([grange,grange_hint]), label1, method1, method2, method3])]))" ] }, { @@ -526,21 +649,16 @@ "\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", + "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 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", + "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 Gaussian functions with the slider next to \n", - "the \"Smearing\" button." + "deviation, $\\sigma$, of the Gaussian functions by adjusting the \"Gaussian $\\sigma$\" slider." ] }, { @@ -567,7 +685,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.12" + "version": "3.10.12" }, "voila": { "authors": "Dou Du and Giovanni Pizzi",