Problem with attaching leads to bilayer graphene

[joyphysics]

Hello everyone. I was trying to attach leads to bilayer graphene but it gives me error.

import pybinding as pb import matplotlib.pyplot as plt from pybinding.repository import graphene from math import pi, sqrt pb.pltutils.use_style()

def bilayer_graphene(): lat = pb.Lattice(a1=[graphene.a, 0], a2=[0.5graphene.a, 0.5sqrt(3)graphene.a]) c0 = 0.335 # [nm] interlayer spacing lat.add_sublattices(('A1', [0, -graphene.a_cc/2, 0]), ('B1', [0, graphene.a_cc/2, 0]), ('A2', [0, graphene.a_cc/2, -c0]), ('B2', [0, 3graphene.a_cc/2, -c0])) lat.register_hopping_energies({'t': graphene.t, 't_layer': -0.4}) lat.add_hoppings(

layer 1

    ([ 0,  0], 'A1', 'B1', 't'),
    ([ 1, -1], 'A1', 'B1', 't'),
    ([ 0, -1], 'A1', 'B1', 't'),
    # layer 2
    ([ 0,  0], 'A2', 'B2', 't'),
    ([ 1, -1], 'A2', 'B2', 't'),
    ([ 0, -1], 'A2', 'B2', 't'),
    # interlayer
    ([ 0,  0], 'B1', 'A2', 't_layer')
)
lat.min_neighbors = 2
return lat

model = pb.Model( bilayer_graphene(), pb.rectangle(2,1) # nm )

model.attach_lead(direction=-1, contact=pb.line([-.5, -0.5 ], [-.5, 0.6]))

model.attach_lead(direction=1, contact=pb.Polygon([[1,0.4, 0], [1,-0.4, 0], [1,-0.4, -0.33], [1,0.4, -0.33]])) plt.figure(figsize=(10, 6)) # make the figure wider model.plot() plt.show()

Error: RuntimeError: Can't attach lead: no sites in lead junction

can anyone suggest me how to solve this problem ?

Also can one attach leads to only a single layer ?

[BertJorissen]

The Polygon-shape only works in 2D, in the x-y-plane. For 2D lattices like yours, you have to use the Line.

The Leads use the Lattice specified in the model. It is thus not straight forward to add a second layer. What you could do, is fold a nanoribbon and add the additional hoppings yourself. I took the liberty to code this for you:

import pybinding as pb
import numpy as np
from pybinding.repository.graphene import monolayer, a, a_cc

c0 = 0.335
eps = 10 ** -2
lead_length = 5
overlap_length = 5

def separate(overlap):
    @pb.site_position_modifier
    def double_layer(x, y, z, sub_id):
        bottom_layer = x > 0
        x[bottom_layer] = x[bottom_layer] - a / 2 - a * np.ceil(overlap / a)
        y[bottom_layer] = y[bottom_layer] + a_cc / 2
        z[bottom_layer] = z[bottom_layer] - c0
        return x, y, z

    @pb.hopping_generator("t_layer", -0.4)
    def interlayer_generator(x, y, z):
        overlap_region = np.logical_and(x > - overlap, x < 0)
        layer1 = np.logical_and(z == 0, overlap_region)
        layer2 = np.logical_and(z != 0, overlap_region)
        idx_1 = np.flatnonzero(layer1)
        idx_2 = np.flatnonzero(layer2)
        from_list = []
        to_list = []
        to_x = x[layer2]
        to_y = y[layer2]
        for from_x, from_y, idx_f in zip(x[layer1], y[layer1], idx_1):
            idx_x = idx_2[np.abs(from_x - to_x) < eps]
            idx_y = idx_2[np.abs(from_y - to_y) < eps]
            if len(idx_x) > 0 and len(idx_y) > 0:
                idx = np.intersect1d(idx_x, idx_y)
                if len(idx) == 1:
                    from_list.append(idx_f)
                    to_list.append(idx[0])
        print(len(from_list))
        return from_list, to_list

    @pb.site_state_modifier
    def remove_bond(state, x, y, z):
        state[np.logical_and(z != 0, x < - overlap)] = False
        return state
    return double_layer, remove_bond, interlayer_generator

def make_model(overlap, single_lead, width):
    model_out = pb.Model(
        monolayer(),
        pb.rectangle(overlap * 2 + single_lead * 2, width),
        separate(overlap)
    )
    model_out.attach_lead(direction=-1, contact=pb.line((-overlap - single_lead, -width/2),
                                                        (-overlap - single_lead, width/2)))
    model_out.attach_lead(direction=+1, contact=pb.line((single_lead, -width/2+a_cc / 2),
                                                        (single_lead, width/2+a_cc / 2)))
    return model_out

model = make_model(a * 5, a * 5, a_cc*2*5*1.1)
model.plot()
model.leads[0].plot_contact()
model.leads[1].plot_contact()

Best Bert

[joyphysics]

Thank you so much. This is exactly what I am looking for. You are awesome... Also the last thing I am looking for is that, can we make the leads hopping different from the system? Or can we take a different lattice to construct the leads like what we can do in KWANT ?

[BertJorissen]

Currently, the lead always is the same material as in the 'lattice' used by pb.Model(). I close here this issue, please refer to the Gitter for further questions about the problems you encounter.