.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "api/gallery/tight_binding/plot_tetrahedron_method.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_api_gallery_tight_binding_plot_tetrahedron_method.py: The Tetrahedron Method for the Density of States ===================================================== The density of states :math:`\rho(E) = \sum_n\int\frac{d\mathbf{k}}{V_{BZ}}\, \delta(E - E_n(\mathbf{k}))` needs an integral over the Brillouin zone, but a calculation only knows the bands on a finite mesh. Broadening each eigenvalue into a Gaussian of width :math:`\sigma` turns the sum into a smooth curve, at the price of smearing exactly the features that matter: band edges, and the van Hove singularities where :math:`\nabla_k E = 0`. Jepsen and Andersen (1971) and Lehmann and Taut (1972) found the way out: split the zone into tetrahedra with mesh points at their corners, interpolate each band *linearly* inside every tetrahedron, and integrate the density of states of that interpolation exactly -- it is a piecewise polynomial in :math:`E`, known in closed form. Blochl, Jepsen and Andersen (1994) gave the standard formulas. There is no :math:`\sigma` to choose, and the error falls with the mesh instead of with a width. :func:`tbkit.dos.tetrahedron_dos` implements the linear method in 1, 2 and 3 dimensions (segments, triangles, tetrahedra), and :meth:`~tbkit.kspace.KSpace.plot_dos` uses it with ``kernel='tetrahedron'``. .. GENERATED FROM PYTHON SOURCE LINES 25-38 .. code-block:: Python import numpy as np import matplotlib.pyplot as plt from scipy.special import ellipk from mpl_toolkits.mplot3d.art3d import Line3DCollection import tbkit.dos as dos import tbkit.lattices as lattices from tbkit.lattice import Lattice from tbkit.bridges import finite_system from tbkit.kspace import KSpace from tbkit.plot import Plot from tbkit.system import System .. GENERATED FROM PYTHON SOURCE LINES 39-47 The square lattice: a logarithmic van Hove singularity ----------------------------------------------------------- :math:`E = -2(\cos k_x + \cos k_y)` has saddle points at X, where the density of states diverges logarithmically. The exact result is :math:`\rho(E) = K(1 - E^2/16)/2\pi^2` per site, with :math:`K` the complete elliptic integral. On the same 80 x 80 mesh, the tetrahedron method is within 0.3% of it (away from the divergence itself), the Gaussian of width one mesh spacing in energy off by 6%. .. GENERATED FROM PYTHON SOURCE LINES 47-51 .. code-block:: Python sq = KSpace(lattices.square()) sq.set_hopping([{'i': 0, 'j': 0, 'R': (1, 0), 't': -1.}, {'i': 0, 'j': 0, 'R': (0, 1), 't': -1.}]) .. GENERATED FROM PYTHON SOURCE LINES 52-54 The square lattice, one site per cell with nearest-neighbour hopping, on :math:`6\times6` cells (:func:`~tbkit.bridges.finite_system`): .. GENERATED FROM PYTHON SOURCE LINES 54-57 .. code-block:: Python fig_lat = Plot(finite_system(sq, (6, 6))).lattice(plt_hop=True, ms=10, figsize=(4, 4)) .. image-sg:: /api/gallery/tight_binding/images/sphx_glr_plot_tetrahedron_method_001.png :alt: plot tetrahedron method :srcset: /api/gallery/tight_binding/images/sphx_glr_plot_tetrahedron_method_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 58-59 Its density of states, exact and on a mesh: .. GENERATED FROM PYTHON SOURCE LINES 59-94 .. code-block:: Python e_grid = np.linspace(-3.9, 3.9, 157) exact = ellipk(1 - e_grid ** 2 / 16) / (2 * np.pi ** 2) away = np.abs(e_grid) > 0.3 nk = 80 en = sq.mesh_bands(nk) _, rho_tet = dos.tetrahedron_dos(en.reshape(nk, nk, 1), e_grid) _, rho_gauss = dos.density_of_states(en, e_grid, broadening=8. / nk) rho_tet, rho_gauss = rho_tet / nk ** 2, rho_gauss / nk ** 2 # per site err_tet = np.abs(rho_tet - exact)[away].max() / exact[away].max() err_gauss = np.abs(rho_gauss - exact)[away].max() / exact[away].max() print('Square lattice, {0} x {0} mesh: tetrahedron error {1:.2%}, Gaussian error {2:.2%}'.format( nk, err_tet, err_gauss)) assert err_tet < 5e-3 and err_gauss > 5e-2 # the peak keeps growing with the mesh, like the logarithm it is peaks = [] for n in (20, 40, 80): en_n = sq.mesh_bands(n).reshape(n, n, 1) peaks.append(dos.tetrahedron_dos(en_n, np.array([-0.01, 0., 0.01]))[1][1] / n ** 2) print('rho(0) on 20, 40, 80 meshes:', np.round(peaks, 3)) assert peaks[0] < peaks[1] < peaks[2] fig, ax = plt.subplots(figsize=(7, 4.5)) ax.plot(e_grid, exact, 'k', lw=3, label='exact') ax.plot(e_grid, rho_gauss, 'r--', label=r'Gaussian, $\sigma = 0.1$') ax.plot(e_grid, rho_tet, 'b', label='tetrahedron') ax.set_ylim(0, 0.5) ax.set_xlabel('$E/t$') ax.set_ylabel(r'$\rho(E)$ per site') ax.set_title('Square lattice, {0} x {0} k-points'.format(nk)) ax.legend() fig.set_layout_engine('tight') .. image-sg:: /api/gallery/tight_binding/images/sphx_glr_plot_tetrahedron_method_002.png :alt: Square lattice, 80 x 80 k-points :srcset: /api/gallery/tight_binding/images/sphx_glr_plot_tetrahedron_method_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Square lattice, 80 x 80 mesh: tetrahedron error 0.26%, Gaussian error 6.21% rho(0) on 20, 40, 80 meshes: [0.384 0.442 0.467] .. GENERATED FROM PYTHON SOURCE LINES 96-103 The simple cubic lattice: band edges ------------------------------------------- In 3D the density of states starts as :math:`\sqrt{E - E_0}` at a band edge and has kinks (van Hove points) inside the band. The Gaussian leaks states below the band bottom :math:`E_0 = -6t`; the tetrahedron method puts none there, and converges to a 64 x 64 x 64 reference about ten times faster. .. GENERATED FROM PYTHON SOURCE LINES 103-108 .. code-block:: Python cubic = KSpace(Lattice(unit_cell=[{'tag': 'a', 'r0': (0., 0., 0.)}], prim_vec=[(1., 0., 0.), (0., 1., 0.), (0., 0., 1.)])) cubic.set_hopping([{'i': 0, 'j': 0, 'R': R, 't': -1.} for R in [(1, 0, 0), (0, 1, 0), (0, 0, 1)]]) .. GENERATED FROM PYTHON SOURCE LINES 109-111 The simple cubic lattice, one site per cell and six nearest neighbours, here :math:`4\times4\times4` sites of it: .. GENERATED FROM PYTHON SOURCE LINES 111-125 .. code-block:: Python lat_vis = Lattice(unit_cell=cubic.lat.unit_cell, prim_vec=cubic.lat.prim_vec) lat_vis.get_lattice(4, 4, 4) vis = System(lat_vis) vis.set_hopping([{'n': 1, 't': 1.}]) pos = np.column_stack([lat_vis.coor[f] for f in 'xyz']) bonds = np.stack([pos[vis.hop['i'].astype(int)], pos[vis.hop['j'].astype(int)]], axis=1) fig_lat = plt.figure(figsize=(4.5, 4.5)) ax_lat = fig_lat.add_subplot(projection='3d') ax_lat.add_collection3d(Line3DCollection(bonds, colors='k', linewidths=1)) ax_lat.scatter(*pos.T, s=30, c='b', depthshade=False) ax_lat.set_aspect('equal') ax_lat.set_axis_off() .. image-sg:: /api/gallery/tight_binding/images/sphx_glr_plot_tetrahedron_method_003.png :alt: plot tetrahedron method :srcset: /api/gallery/tight_binding/images/sphx_glr_plot_tetrahedron_method_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 126-127 Its density of states near the band bottom: .. GENERATED FROM PYTHON SOURCE LINES 127-146 .. code-block:: Python e_grid = np.linspace(-6.5, 6.5, 131) ref = dos.tetrahedron_dos(cubic.mesh_bands(64).reshape(64, 64, 64, 1), e_grid)[1] / 64 ** 3 nk = 32 en = cubic.mesh_bands(nk) rho_tet = dos.tetrahedron_dos(en.reshape(nk, nk, nk, 1), e_grid)[1] / nk ** 3 rho_gauss = dos.density_of_states(en, e_grid, broadening=12. / nk)[1] / nk ** 3 below = e_grid < -6.05 err_tet = np.abs(rho_tet - ref).max() / ref.max() err_gauss = np.abs(rho_gauss - ref).max() / ref.max() print('Simple cubic, 32^3 mesh: error {:.2%} (tetrahedron), {:.2%} (Gaussian)'.format(err_tet, err_gauss)) print('States below the band bottom: {:.1e} (tetrahedron), {:.1e} (Gaussian)'.format( rho_tet[below].max(), rho_gauss[below].max())) assert err_tet < 1e-2 and err_gauss > 5e-2 assert rho_tet[below].max() == 0. and rho_gauss[below].max() > 1e-3 # one state per site in all assert np.isclose(np.sum(rho_tet) * (e_grid[1] - e_grid[0]), 1.) .. rst-class:: sphx-glr-script-out .. code-block:: none Simple cubic, 32^3 mesh: error 0.53% (tetrahedron), 9.71% (Gaussian) States below the band bottom: 0.0e+00 (tetrahedron), 4.9e-03 (Gaussian) .. GENERATED FROM PYTHON SOURCE LINES 147-149 The same in one call: :meth:`~tbkit.kspace.KSpace.plot_dos` with ``kernel='tetrahedron'`` (no broadening to choose). .. GENERATED FROM PYTHON SOURCE LINES 149-151 .. code-block:: Python fig2 = cubic.plot_dos(nk=32, kernel='tetrahedron', fs=14, figsize=(7, 4.5)) .. image-sg:: /api/gallery/tight_binding/images/sphx_glr_plot_tetrahedron_method_004.png :alt: Density of states :srcset: /api/gallery/tight_binding/images/sphx_glr_plot_tetrahedron_method_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.891 seconds) .. _sphx_glr_download_api_gallery_tight_binding_plot_tetrahedron_method.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_tetrahedron_method.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_tetrahedron_method.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_tetrahedron_method.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_