Analytical Hamiltonian

NOTE: The diagonal terms of analytical Hamiltonian produced by print_hk is incorrect in TBPLaS v1.4. Make sure you have installed v1.4.1 or newer versions when running this example.

In this tutorial, we show how to get the analytical Hamiltonian and how to use it in calculations. The PrimitiveCell class has a print_hk method which can print the formula of Hamiltonian, while the DiagSolver class can utilize user-defined analytical Hamiltonian to calculate the eigenvalues and eigenstates. The script is located at examples/advanced/analy_ham.py. We begin with importing the necessary packages:

from math import cos, sin, pi
import numpy as np
import tbplas as tb

Derivation of Hamiltonian

The analytical Hamiltonian can be obtained by calling the print_hk method, after the primitive cell has been successfully created. We show this with the 2-band model of monolayer graphene:

1cell = tb.make_graphene_diamond()
2cell.print_hk(convention=1)
3cell.print_hk(convention=2)

The convention argument specifies the convention to evaluate the elements of Hamiltonian, and should be either 1 or 2. For convention 1 the elements are

\[H_{ij}(\mathbf{k}) = \sum_{\mathbf{R}} H_{ij}(\mathbf{R})\mathrm{e}^{\mathrm{i}\mathbf{k}\cdot(\mathbf{R}+\tau_j-\tau_i)}\]

while for convention 2

\[H_{ij}(\mathbf{k}) = \sum_{\mathbf{R}} H_{ij}(\mathbf{R})\mathrm{e}^{\mathrm{i}\mathbf{k}\cdot\mathbf{R}}\]

The output for convention 1 should look like

1ham[0, 0] = (0.0)
2ham[1, 1] = (0.0)
3ham[0, 1] = ((-2.7+0j) * exp_i(0.3333333333333333 * ka + 0.3333333333333333 * kb)
4 + (-2.7-0j) * exp_i(-0.6666666666666667 * ka + 0.3333333333333333 * kb)
5 + (-2.7-0j) * exp_i(0.3333333333333333 * ka - 0.6666666666666667 * kb))
6ham[1, 0] = ham[0, 1].conjugate()
7with exp_i(x) := cos(2 * pi * x) + 1j * sin(2 * pi * x)

and for convention 2

1ham[0, 0] = (0.0)
2ham[1, 1] = (0.0)
3ham[0, 1] = ((-2.7+0j) * 1
4 + (-2.7-0j) * exp_i(-ka)
5 + (-2.7-0j) * exp_i(-kb))
6ham[1, 0] = ham[0, 1].conjugate()
7with exp_i(x) := cos(2 * pi * x) + 1j * sin(2 * pi * x)

where the formula for each element of the Hamiltonian matrix is printed. Note that we have defined the function \(\exp(\mathrm{i} \cdot 2\pi \cdot x)\) using Euler formula since the exp function of numpy is rather slow if the input is not an array. We will utilize this function in the next section.

Usage of Hamiltonian

TBPLaS provides the FakePC class for holding the analytical Hamiltonian. Users should derive their own class from this class and implement the analytical Hamiltonian in the set_ham_dense and set_ham_csr methods. The set_ham_dense method should accept the fractional coordinate of k-point and the Hamiltonian matrix as input, and modifies the Hamiltonian in-place. On the contrary, the set_ham_csr function accepts only the fractional coordinate of k-point as input, and returns a new Hamiltonian matrix in CSR format.

We define the following function and class from the analytical Hamiltonians in previous section

 1def exp_i(x: float) -> complex:
 2    """
 3    Evaluate exp(i*2pi*x) using Euler formula.
 4
 5    :param x: incoming x
 6    :return: exp(i*2pi*x)
 7    """
 8    return cos(2 * pi * x) + 1j * sin(2 * pi * x)
 9
10
11class FakePC(tb.FakePC):
12    def set_ham_dense(self, kpt: np.ndarray,
13                      ham: np.ndarray,
14                      convention: int = 1) -> None:
15        ka, kb = kpt.item(0), kpt.item(1)
16        ham[0, 0] = 0.0
17        ham[1, 1] = 0.0
18        if convention == 1:
19            ham[0, 1] = -2.7 * (exp_i(1. / 3 * ka + 1. / 3 * kb) +
20                                exp_i(-2. / 3 * ka + 1. / 3 * kb) +
21                                exp_i(1. / 3 * ka - 2. / 3 * kb))
22            ham[1, 0] = ham[0, 1].conjugate()
23        else:
24            ham[0, 1] = -2.7 * (1.0 + exp_i(-ka) + exp_i(-kb))
25            ham[1, 0] = ham[0, 1].conjugate()

To demonstrate the usage of analytical Hamiltonian, we create a fake graphene primitive cell with 2 orbitals. Then we create a solver from the DiagSolver class using the fake primitive cell with analytical Hamiltonian

1# Create fake primitive cell, solver and visualizer
2vectors = tb.gen_lattice_vectors(a=0.246, b=0.246, c=1.0, gamma=60)
3cell = FakePC(num_orb=2, lat_vec=vectors, unit=tb.NM)
4solver = tb.DiagSolver(cell)
5vis = tb.Visualizer()

The band structure from the analytical Hamiltonian can be obtained by calling the calc_bands method of the DiagSolver class

 1# Evaluation of band structure
 2k_points = np.array([
 3    [0.0, 0.0, 0.0],
 4    [2. / 3, 1. / 3, 0.0],
 5    [1. / 2, 0.0, 0.0],
 6    [0.0, 0.0, 0.0],
 7])
 8k_label = ["G", "M", "K", "G"]
 9k_path, k_idx = tb.gen_kpath(k_points, [40, 40, 40])
10for convention in (1, 2):
11    k_len, bands = solver.calc_bands(k_path, convention)[:2]
12    vis.plot_bands(k_len, bands, k_idx, k_label)

The evaluation of DOS is similar

1# Evaluation of DOS
2k_mesh = tb.gen_kmesh((120, 120, 1))
3for convention in (1, 2):
4    energies, dos = solver.calc_dos(k_mesh, convention=convention)
5    vis.plot_dos(energies, dos)

Both conventions produce the same band structure and DOS as in Band structure and DOS.