Hopping generator

The @hopping_generator can be used to create new hoppings independent of the main lattice definition. It’s especially useful for creating additional local hoppings, e.g. to model defects. Here, we present a way create twisted bilayer graphene with an arbitrary rotation angle \(\theta\).

We start with two unconnected layers of graphene. A @site_position_modifier is applied to rotate just one layer. Then, a @hopping_generator finds and connects the layers via site pairs which satisfy the given criteria. The newly created hoppings all have identical energy at first. Finally, a @hopping_energy_modifier to applied to set the new interlayer hopping energy to the desired distance-dependent value.

This is an experimental feature, presented as is, without any additional support.

Source code

"""Construct a circular flake of twisted bilayer graphene (arbitrary angle)"""
import math
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import cKDTree

import pybinding as pb

c0 = 0.335  # [nm] graphene interlayer spacing


def two_graphene_monolayers():
    """Two individual layers of monolayer graphene without any interlayer hopping"""
    from pybinding.repository.graphene.constants import a_cc, a, t

    lat = pb.Lattice(a1=[a/2, a/2 * math.sqrt(3)], a2=[-a/2, a/2 * math.sqrt(3)])
    lat.add_sublattices(('A1', [0,  -a_cc,   0]),
                        ('B1', [0,      0,   0]),
                        ('A2', [0,      0, -c0]),
                        ('B2', [0,   a_cc, -c0]))
    lat.register_hopping_energies({'gamma0': t})
    lat.add_hoppings(
        # layer 1
        ([ 0,  0], 'A1', 'B1', 'gamma0'),
        ([ 0, -1], 'A1', 'B1', 'gamma0'),
        ([-1,  0], 'A1', 'B1', 'gamma0'),
        # layer 2
        ([ 0,  0], 'A2', 'B2', 'gamma0'),
        ([ 0, -1], 'A2', 'B2', 'gamma0'),
        ([-1,  0], 'A2', 'B2', 'gamma0'),
        # not interlayer hopping
    )
    lat.min_neighbors = 2
    return lat


def twist_layers(theta):
    """Rotate one layer and then a generate hopping between the rotated layers"""
    theta = theta / 180 * math.pi  # from degrees to radians

    @pb.site_position_modifier
    def rotate(x, y, z):
        """Rotate layer 2 by the given angle `theta`"""
        layer2 = (z < 0)
        x0 = x[layer2]
        y0 = y[layer2]
        x[layer2] = x0 * math.cos(theta) - y0 * math.sin(theta)
        y[layer2] = y0 * math.cos(theta) + x0 * math.sin(theta)
        return x, y, z

    @pb.hopping_generator('interlayer', energy=0.1)  # eV
    def interlayer_generator(x, y, z):
        """Generate hoppings for site pairs which have distance `d_min < d < d_max`"""
        positions = np.stack([x, y, z], axis=1)
        layer1 = (z == 0)
        layer2 = (z != 0)

        d_min = c0 * 0.98
        d_max = c0 * 1.1
        kdtree1 = cKDTree(positions[layer1])
        kdtree2 = cKDTree(positions[layer2])
        coo = kdtree1.sparse_distance_matrix(kdtree2, d_max, output_type='coo_matrix')

        idx = coo.data > d_min
        abs_idx1 = np.flatnonzero(layer1)
        abs_idx2 = np.flatnonzero(layer2)
        row, col = abs_idx1[coo.row[idx]], abs_idx2[coo.col[idx]]
        return row, col  # lists of site indices to connect

    @pb.hopping_energy_modifier
    def interlayer_hopping_value(energy, x1, y1, z1, x2, y2, z2, hop_id):
        """Set the value of the newly generated hoppings as a function of distance"""
        d = np.sqrt((x1-x2)**2 + (y1-y2)**2 + (z1-z2)**2)
        interlayer = (hop_id == 'interlayer')
        energy[interlayer] = 0.4 * c0 / d[interlayer]
        return energy

    return rotate, interlayer_generator, interlayer_hopping_value


model = pb.Model(
    two_graphene_monolayers(),
    pb.circle(radius=1.5),
    twist_layers(theta=21.798)
)
plt.figure(figsize=(6.5, 6.5))
model.plot(hopping=dict(width=1.6, cmap='auto'))
plt.title(r"$\theta$ = 21.798 $\degree$")
plt.show()
Twisted bilayer graphene for arbitrary angles
model = pb.Model(
    two_graphene_monolayers(),
    pb.circle(radius=1.5),
    twist_layers(theta=12.95)
)
plt.figure(figsize=(6.5, 6.5))
model.plot(hopping=dict(width=1.6, cmap='auto'))
plt.title(r"$\theta$ = 12.95 $\degree$")
plt.show()
Twisted bilayer graphene for arbitrary angles