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
while for convention 2
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.