Run in Binder Download .ipynb Download .py

The Mesh class#

This tutorial will show you how to create custom meshes and winding paths around the the combined Brillouin zone (BZ) and parameter space.

from pythtb import Mesh
import numpy as np

Mesh accepts the following parameters:

  • axis_types: List of strings defining the type of each axis; options are "k" for crystal momentum (reduced units) and "l" for adiabatic parameters (\(\lambda\))

  • axis_names (optional): List of strings defining the names of each axis, primarily used for adiabatic \(\lambda\) cycles (see Three-site Thouless pump for example)

  • dim_k (optional): Integer defining the dimension of k-space; required if the number of ‘k’ axes is less than dim_k (e.g., for a 1D path in a 2D BZ)

To build a simple regular k-mesh, we will

  • initialize a Mesh object with the desired axis types (and names and k-dimension if needed)

  • use the build_grid() method to create a regular grid with a specified shape

  • (optional) specify whether the grid is \(\Gamma\)-centered

  • (optional) specify whether the k-points should include the endpoints or not

In this case we will uniformly sample a 2D BZ with a 10x10 grid. Each leading axis in mesh.points corresponds to one of the k-axes defined during initialization. The last axis corresponds to the Cartesian components of each k-point in reduced units.

mesh = Mesh(axis_types=["k", "k"])
mesh.build_grid(shape=(10, 10), gamma_centered=True, k_endpoints=False)

print(mesh.points.shape)

(10, 10, 2)

In the following cell, we create a mesh representating a 1-dimensional path in a 2D k-space. This means that there will only be a single k axis. We will use the method build_custom() to specifiy the mesh points directly. Here we create a path from the \(\Gamma\) = [0, 0] point through the \(X = \) [0.5, 0.5] point to the \(\Gamma = \) [1,1] point in the next cell.

mesh = Mesh(["k"], dim_k=2)
points = np.linspace(
    [0, 0], [1, 1], 10, endpoint=False
)  # path from (0, 0) to (0.5, 0.5) to (1, 1)
mesh.build_custom(points)

We access the mesh points via the points attribute. Notice that each point is a 2D vector corresponding to the two k-space dimensions.

print("Mesh shape:", mesh.shape)
print("Mesh points along the path from (0,0) to (1,1):")
print(mesh.points)

Mesh shape: (10, 2)
Mesh points along the path from (0,0) to (1,1):
[[0.  0. ]
 [0.1 0.1]
 [0.2 0.2]
 [0.3 0.3]
 [0.4 0.4]
 [0.5 0.5]
 [0.6 0.6]
 [0.7 0.7]
 [0.8 0.8]
 [0.9 0.9]]

We can see the mesh information by printing the Mesh object. This tells us

  • the type of mesh (either “path” or “grid”)

  • the number of dimensions in k-space and parameter space

  • The total number of mesh points

  • The shape of the mesh array

  • The k-axes (type (k), name, and number of points)

  • The parameter space axes (type (l), name, and number of points)

  • The axes that are looped (periodic) or not

  • The loops in the Mesh (if any)

print(mesh)

Mesh Summary
========================================
Type: path
Dimensionality: 2 k-dim(s) + 0 λ-dim(s)
Number of mesh points: 10
Full shape: (10, 2)
k-axes: [Axis(type=k, name=k_0, size=10)]
λ-axes: []
Loops: None

The code doesn’t inherently know which axes should wind the BZ or parameter space if the endpoints are not included, so we have to specify this manually. This is important for properly handling the periodic boundary conditions with the Bloch wavefunctions.

We can declare an axis as winding the BZ by calling the method loop. It takes the following parameters:

  • axis_idx: The index of the axis to wind (0-indexed)

  • component_idx: The index of the component of the combined (\(k\), \(\lambda\)) vector (0-indexed)

  • winds_bz: Boolean indicating whether the axis winds the BZ by a reciprocal lattice vector (True)

  • closed: Boolean indicating whether the axis includes the endpoints (is closed) (False).

We set axis_idx=0 to wind the first (and only) axis of our mesh, and component_idx=0 to wind around the first k-space dimension. We do the same thing for the second k-space dimension by calling loop again with component_idx=1. After calling loop, we can see that the mesh information now indicates that both axes wind the BZ.

Note

When using with k-points, we should set winds_bz=True if the axis winds around the BZ by a reciprocal lattice vector. This ensures that the Bloch wavefunctions are treated with the correct periodic boundary conditions. If the axis does not wind the BZ, we should set winds_bz=False. This would be the case for loops internal to the BZ, such as small circular paths around Dirac points, for example.

# Wind the first k-space dimension
mesh.loop(axis_idx=0, component_idx=0, winds_bz=True, closed=False)

# Now wind the second k-space dimension
mesh.loop(axis_idx=0, component_idx=1, winds_bz=True, closed=False)
print(mesh)

Mesh Summary
========================================
Type: path
Dimensionality: 2 k-dim(s) + 0 λ-dim(s)
Number of mesh points: 10
Full shape: (10, 2)
k-axes: [Axis(type=k, name=k_0, size=10)]
λ-axes: []
Loops: (axis 0, comp 0, winds_bz=yes, closed=no), (axis 0, comp 1, winds_bz=yes, closed=no)

We will do the same thing as we have done above, but this time include the endpoints of the BZ.

Note

The Mesh class will automatically detect an axis as winding the BZ if the start and end points are the same modulo 1 (i.e. the axis is closed) for the k-components. So even for \(\Gamma\)-centered meshes, where the k-component(s) range from \([-0.5, 0.5]\), the class will still detect that the axis winds the BZ if the endpoints (\(0.5\)) are included.

mesh = Mesh(dim_k=2, axis_types=["k"])

# Path from (-0.5,-0.5) to (0.5, 0.5) including endpoints (endpoint=True)
points = np.linspace([-0.5, -0.5], [0.5, 0.5], 10, endpoint=True)

mesh.build_custom(points)
print(mesh)

Mesh Summary
========================================
Type: path
Dimensionality: 2 k-dim(s) + 0 λ-dim(s)
Number of mesh points: 10
Full shape: (10, 2)
k-axes: [Axis(type=k, name=k_0, size=10)]
λ-axes: []
Loops: (axis 0, comp 0, winds_bz=yes, closed=yes), (axis 0, comp 1, winds_bz=yes, closed=yes)