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, <adelp...@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 >> >