Dear Paweł, The issue lies within the `lead.substituted` part. Until that bit is fixed, you can go ahead and define the leads' functions and Peierl's phase functions with different names for each lead. I have checked that the below example works, it will be easy to adapt it to your case. Best regards, Pablo ``` import numpy as np import kwant def hopping(a, b, t, peierls): return -t * peierls(a, b) def hopping_lead_1(a, b, t, peierls_lead_1): return -t * peierls_lead_1(a, b) def hopping_lead_2(a, b, t, peierls_lead_2): return -t * peierls_lead_2(a, b) size = 30 def make_system(hopping): lat = kwant.lattice.square(norbs=1) syst = kwant.Builder() syst[(lat(i, j) for i in range(size) for j in range(size))] = 0 syst[lat.neighbors()] = hopping return syst def make_lead(hopping_lead, lattice_direction=(-1, 0)): symmetry = kwant.TranslationalSymmetry(lattice_direction) lat = kwant.lattice.square(norbs=1) lead = kwant.Builder(symmetry=symmetry) lead[(lat(i, j) for j in range(size) for i in range(abs(lattice_direction[0])))] = 0 lead[lat.neighbors()] = hopping_lead return lead syst = make_system(hopping) lead_1 = make_lead(hopping_lead_1, lattice_direction=(1, 0)) lead_2 = make_lead(hopping_lead_2, lattice_direction=(-1, 0)) syst.attach_lead(lead_1) syst.attach_lead(lead_2) syst = syst.finalized() gauge = kwant.physics.magnetic_gauge(syst) def B_syst(pos): return np.pi / 10 peierls_syst, peierls_lead_1, peierls_lead_2 = gauge(B_syst, B_syst, B_syst) params = dict(t=1, peierls=peierls_syst, peierls_lead_1=peierls_lead_1, peierls_lead_2=peierls_lead_2) smatrix = kwant.smatrix(syst, params=params, energy=1) print(smatrix.transmission(0, 1)) ```