Hello everyone, i would appreciate any suggestion you could give me to solve this. Well, im trying to create a structure with 3 basis atoms in the unit cell, each atom has a different number of orbitals a=3, b=6 and c=3. I have tried with polyatomic module and creating one lattice for each atom with no success The thing is that i cant figure it out how to create a lattice with different orbitals in each site. I have the on-site and hopping energies matrix so i need to be able to introduce this values. Im pretty new using kwant so any advice is appreciated. Here is the code i've working on (ignore the matrix elements for the moment): import kwant import tinyarray import numpy from matplotlib import pyplot VecPrim = [(3.176, 0, 0), (1.588, 2.7504, 0), (0, 0, 17.49)] base = [(0, 0, 0), (1.588, 0.9168, 0.8745), (0, 0, 1.749)] ### This is the part that i cant get it right, i have tried kwant.lattice.polyatomic ### and defining 3 lattices one for each atom, but no success a = kwant.lattice.general(prim_vecs = VecPrim, basis = (0, 0, 0)) Cpx, Cpy, Cpz, = a.sublattices b = kwant.lattice.general(prim_vecs = VecPrim, basis = (1.588, 0.9168, 0.8745)) xy, xz, yz, x2y2, z2, s = b.sublattices c = kwant.lattice.general(prim_vecs = VecPrim, basis = (0, 0, 1.749)) Npx, Npy, Npz = c.sublattices ### Ignore from here ### def make_cuboid(t=1.0, a=15, b=10, c=5): def cuboid_shape(pos): x, y, z = pos return 0 <= x < a and 0 <= y < b and 0 <= z < c syst = kwant.Builder() syst[lat.shape(cuboid_shape, (0, 0, 0))] = 4 * t syst[lat(0, 0, 0)] = numpy.array([[-1.888842, -0.014064-0.000409j, 0.026908+0.003422j, 0, 0, 0, 0, 0, 0, 0, 0, 0], [-0.014064+0.000409j, -1.784936, -0.031384+0.021565j, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0.026908-0.003422j, -0.031384-0.021565j, 0.710443, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]) syst[lat(0, 1, 0)] = numpy.array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 7.168134, -0.005189-0.047376j, 0.001352+0.000103j, -0.001189-0.00095j, -0.048166-0.007663j, 0.429882+0.00052j, 0, 0, 0], [0, 0, 0, -0.005189+0.047376j, 0.245358, 0.113162-0.072937j, 0.660849-0.06348j, -0.017571+0.035587j, 0.020961-0.017679j, 0, 0, 0], [0, 0, 0, 0.001352-0.000103j, 0.113162+0.072937j, 1.338059, -0.003126-0.016285j, -0.624702-0.065652j, 0.002353+0.001639j, 0, 0, 0], [0, 0, 0, -0.001189+0.00095j, 0.660849+0.06348j, -0.003126+0.016285j, 1.325337, 0.110694+0.01548j, 0.009041+0.004483j, 0, 0, 0], [0, 0, 0, -0.048166+0.007663j, -0.017571-0.035587j, -0.624702+0.065652j, 0.110694-0.01548j, 0.362102, -0.030934+0.002401j, 0, 0, 0], [0, 0, 0, 0.429882-0.00052j, 0.020961+0.017679j, 0.002353-0.001639j, 0.009041-0.004483j, -0.030934-0.002401j, -0.650839, 0, 0, 0] [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]) syst[lat(0, 0, 1)] = numpy.array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, -3.398186, 0.001839+0.001631j, -0.117366+0.00498j], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0.001839-0.001631j, -3.401815, 0.023339+0.009067j], [0, 0, 0, 0, 0, 0, 0, 0, 0, -0.117366-0.00498j, 0.023339-0.009067j, -0.637207]]) syst[lat(0, 0, 1), lat (0, 0, 0)] = ([[0, 0, 0, 0, 0, 0, 0, 0, 0, -1.852931+0.147524j, 0.274586-0.034451j, -0.015051-0.001677j], [0, 0, 0, 0, 0, 0, 0, 0, 0, -0.274338-0.015064j, -1.796449-0.055331j, -0.034396+0.019407j], [0, 0, 0, 0, 0, 0, 0, 0, 0, -0.103127+0.001438j, 0.004816-0.00237j, 3.410475+0.001248j], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [-1.852931-0.147524j, -0.274338+0.015064j, -0.103127-0.001438j, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0.274586+0.034451j, -1.796449+0.055331j, 0.004816+0.00237j, 0, 0, 0, 0, 0, 0, 0, 0, 0], [-0.015051+0.001677j, -0.034396-0.019407j, 3.410475-0.001248j, 0, 0, 0, 0, 0, 0, 0, 0, 0]]) syst[lat(0, 1, 0), lat(0, 0, 0)] = ([[0, 0, 0, 0.003644+0.00088j, -0.305706-0.050514j, 0.638425-0.083597j, -1.105645+0.086829j, 1.463462-0.012994j, -0.01338-0.00264j, 0, 0, 0], [0, 0, 0, -0.046131+0.04229j, -1.963829-0.280977j, -1.171301-0.088146j, -1.112788-0.056971j, 0.174769+0.106872j, -0.058159-0.007987j, 0, 0, 0], [0, 0, 0, 0.23188+0.001924j, -0.597505-0.078772j, -0.146976-0.012943j, -0.210584+0.002248j, 0.177022+0.031485j, 0.953363-0.001855j, 0, 0, 0], [0.003644-0.00088j, -0.046131-0.04229j, 0.23188-0.001924j, 0, 0, 0, 0, 0, 0, 0, 0, 0], [-0.305706+0.050514j, -1.963829+0.280977j, -0.597505+0.078772j, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0.638425+0.083597j, -1.171301+0.088146j, -0.146976+0.012943j, 0, 0, 0, 0, 0, 0, 0, 0, 0], [-1.105645-0.086829j, -1.112788+0.056971j, -0.210584-0.002248j, 0, 0, 0, 0, 0, 0, 0, 0, 0], [1.463462+0.012994j, 0.174769-0.106872j, 0.177022-0.031485j, 0, 0, 0, 0, 0, 0, 0, 0, 0], [-0.01338+0.00264j, -0.058159+0.007987j, 0.953363+0.001855j, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]) syst[lat(0, 0, 1), lat(0, 1, 0)] = ([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0.184152+0.002778j, 0.321743-0.010571j, 0.183538-0.00067j], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0.275145+0.007691j, -0.54144+0.012815j, 0.045847+0.015379j], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0.47922+0.015775j, 1.256012+0.007326j, -0.523174+0.00585j], [0, 0, 0, 0, 0, 0, 0, 0, 0, 1.255561-0.002255j, 1.903063+0.010647j, -0.876392+0.010246j], [0, 0, 0, 0, 0, 0, 0, 0, 0, 1.282837-0.091022j, -0.54898+0.001161j, 0.057341-0.004095j], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0.02765-0.001793j, -0.002902-0.001866j, -0.746426-0.000748j], [0, 0, 0, 0.184152-0.002778j, 0.275145-0.007691j, 0.47922-0.015775j, 1.255561+0.002255j, 1.282837+0.091022j, 0.02765+0.001793j, 0, 0, 0], [0, 0, 0, 0.321743+0.010571j, -0.54144-0.012815j, 1.256012-0.007326j, 1.903063-0.010647j, -0.54898-0.001161j,-0.002902+0.001866j, 0, 0, 0], [0, 0, 0, 0.183538+0.00067j, 0.045847-0.015379j, -0.523174-0.00585j, -0.876392-0.010246j, 0.057341+0.004095j,-0.746426+0.000748j, 0, 0, 0]]) ### 'Till here ##### syst[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = -6 lead = kwant.Builder(kwant.TranslationalSymmetry((-3.176, 0, 0))) def lead_shape(pos): return 0 <= pos[1] < b and 0 <= pos[2] < c lead[lat.shape(lead_shape, (0, 0, 0))] = 4 * t lead[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = -6 syst.attach_lead(lead) syst.attach_lead(lead.reversed()) return syst def plot_conductance(syst, energies): data = [] for energy in energies: smatrix = kwant.smatrix(syst, energy) data.append(smatrix.transmission(1, 0)) pyplot.figure() pyplot.plot(energies, data) pyplot.xlabel("energy [t]") pyplot.ylabel("conductance [e^2/h]") pyplot.show() def main(): syst = make_cuboid() kwant.plot(syst) syst = make_cuboid(a=100, b=28, c=4) def family_colors(site): if site.family == d: return 'yellow' elif site.family == e: return 'gray' else: return 'blue' kwant.plot(syst, site_size=0.25, site_lw=0.025, hop_lw=0.05, site_color=family_colors) syst = syst.finalized() plot_conductance(syst, energies=[0.01 * i - 0.3 for i in range(100)]) if __name__ == '__main__': main()