Dear Adel, Thank you so much for your help. I have some questions about it. Is the graphene "Finit" periodic in y direction? Can I attach the leads to it? I also have other questions related to current-voltage curve and applying magnetic field. I wrote the following code for this aspect but I need to be sure about it. Could you please help me? Thank you in advance ###current-voltage#### def fermi_dirac(energy, mu, T): return 1 / (np.exp((energy - mu) / (k_B * T)) + 1) for voltage in voltages: mu_left = voltage / 2 mu_right = -voltage / 2 # Compute the current using the Landauer-Büttiker formula I = 0 # Loop over energy intervals and compute contributions for e1, e2, g1, g2 in zip(energies[:-1], energies[1:], conductance[:-1], conductance[1:]): # Compute Fermi-Dirac distributions at e1 and e2 for left and right reservoirs f_left_e1 = fermi_dirac(e1, mu_left, temperature) f_right_e1 = fermi_dirac(e1, mu_right, temperature) f_left_e2 = fermi_dirac(e2, mu_left, temperature) f_right_e2 = fermi_dirac(e2, mu_right, temperature) # Current calculation using trapezoidal integration delta_E = e2 - e1 I += ((g1 + g2) / 2) * ((f_left_e1 - f_right_e1) + (f_left_e2 - f_right_e2)) / 2 * delta_E # Convert from eV to Amperes current.append(I * e/h * 1e6) # Multiply by e to convert the current to Amperes return current ###magnetic field### def hopping_B(site1, site2, phi): x1,y1=site1.pos x2,y2=site2.pos y = 0.5 * (y1 + y2) y *= 1e-10 A_x = phi * y peierls = A_x * (x1 - x2) peierls *= 1e-10 return -t * np.exp(1j * 2*np.pi/phi0 * peierls) On Thu, 14 Nov 2024, 10:43 Adel Belayadi, <adelphys@gmail.com> wrote:

Dear Masoumeh, I share an example of graphene with 04 atoms in the unit cell. You can also have the zigzag and armchair ribbon from that. I hope this helps Best Adel

import kwant import matplotlib.pyplot as plt from math import sqrt

def graphene_04_atoms (a=1.42, t=-1, Ep=0, Lx=5, Ly=5, boundary ="zigzag"): a1 = a*sqrt(3) a2 = 3*a

### here premitive vectors are perpondicular: primitive_vector_1.primitive_vector_1 = 0 primitive_vector_1 = (a1, 0) primitive_vector_2 = (0, a2)

## Now let us define the atoms positions of the unit cell Pos_A1 = ( 0, -a/2) Pos_B1 = ( 0, a/2) Pos_A2 = ( a1/2, a) Pos_B2 = ( a1/2, 2*a)

lat = kwant.lattice.general([primitive_vector_1, primitive_vector_2], [Pos_A1, Pos_B1, Pos_A2, Pos_B2]) Sub_A1, Sub_B1, Sub_A2, Sub_B2 = lat.sublattices

#................................................................................... if boundary == "Finit": syst= kwant.Builder() def square(pos): x, y = pos return abs(x)<Lx and abs(y)<Ly syst[lat.shape(square, (0,0))]= Ep #nearest neighbors............................................................. ## Be carfull if you consider second nearest neighbors ## ========================================================= ## Hopping within unit cell ================================ syst[[kwant.builder.HoppingKind((0,0),Sub_A1,Sub_B1)]] = t syst[[kwant.builder.HoppingKind((0,0),Sub_B1,Sub_A2)]] = t syst[[kwant.builder.HoppingKind((0,0),Sub_A2,Sub_B2)]] = t ## Hopping between neighbouring unit cell======================= syst[[kwant.builder.HoppingKind((+1, +1),Sub_A1,Sub_B2)]] = t syst[[kwant.builder.HoppingKind(( 0, +1),Sub_A1,Sub_B2)]] = t syst[[kwant.builder.HoppingKind((+1, 0),Sub_B1,Sub_A2)]] = t

syst.eradicate_dangling() # Plot the closed system without leads. kwant.plot(syst);

if boundary == "zigzag": def square(pos): x, y = pos return abs(y)<Ly sym = kwant.TranslationalSymmetry(primitive_vector_1) syst = kwant.Builder(sym) syst[lat.shape(square, (0,0))]= Ep syst[[kwant.builder.HoppingKind((0,0),Sub_A1,Sub_B1)]] = t syst[[kwant.builder.HoppingKind((0,0),Sub_B1,Sub_A2)]] = t syst[[kwant.builder.HoppingKind((0,0),Sub_A2,Sub_B2)]] = t ## Hopping between neighbouring unit cell======================= syst[[kwant.builder.HoppingKind((+1, +1),Sub_A1,Sub_B2)]] = t syst[[kwant.builder.HoppingKind(( 0, +1),Sub_A1,Sub_B2)]] = t syst[[kwant.builder.HoppingKind((+1, 0),Sub_B1,Sub_A2)]] = t

syst.eradicate_dangling() # Plot the closed system without leads. kwant.plot(syst);

if boundary == "armchair": def square(pos): x, y = pos return abs(x)<Lx sym = kwant.TranslationalSymmetry(primitive_vector_2) syst = kwant.Builder(sym) syst[lat.shape(square, (0,0))]= Ep syst[[kwant.builder.HoppingKind((0,0),Sub_A1,Sub_B1)]] = t syst[[kwant.builder.HoppingKind((0,0),Sub_B1,Sub_A2)]] = t syst[[kwant.builder.HoppingKind((0,0),Sub_A2,Sub_B2)]] = t ## Hopping between neighbouring unit cell======================= syst[[kwant.builder.HoppingKind((+1, +1),Sub_A1,Sub_B2)]] = t syst[[kwant.builder.HoppingKind(( 0, +1),Sub_A1,Sub_B2)]] = t syst[[kwant.builder.HoppingKind((+1, 0),Sub_B1,Sub_A2)]] = t

syst.eradicate_dangling() # Plot the closed system without leads. kwant.plot(syst);

graphene_04_atoms (a=1.42, t=-1, Ep=0, Lx=5, Ly=15, boundary ="Finit"); graphene_04_atoms (a=1.42, t=-1, Ep=0, Lx=5, Ly=15, boundary ="zigzag"); graphene_04_atoms (a=1.42, t=-1, Ep=0, Lx=15, Ly=5, boundary ="armchair");

Le lun. 11 nov. 2024 à 14:20, Masoumeh Alihosseini via Kwant-discuss < kwant-discuss@python.org> a écrit :

...

Hi dear all, I have a question, I want to work with rectangular graphene cell with 4 atoms in the unitcell. Like zigzag nanoribbon but I want to apply periodicity in y direction for both scattering region and leads. Could you please help me? Thank you, Masoumeh