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.
"""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()
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()