Dear Kwant developers, I am a new user of Kwant, and I want to build a polyatomic model by directly loading the coordinate and Hamiltonian file. Since each atom site could have more than one orbital (like s, px, py, pz ..), and different kind of atom sites could have a different number of orbitals. Can I just assign the Hamiltonian for each atom by just giving the exact coordinates? For example, for a model system: H 0.0 0.0 0.0 C 2.0 0.0 0.0 C 8.0 0.0 0.0 with the basis function H(H-1s) C(C-2s,C-2px,C-2py, C-2pz), and I set the "family" of the atomic site first, basis[0] -> H-1s basis[1] -> C-2s basis[2] -> C-2px basis[3] -> C-2py basis[4] -> C-2pz and then I set the onsite energy and offsite coupling for each atom by directly assign the atomic coordinate, i.e. for the 2nd C atom, I use the "family" of basis[1] basis[2] basis[3] basis[4] with the coordinate tag to be '2.0 , 0.0 , 0.0'. Enclosed please find my simple kwant script. Thanks in advance for your help. Hang Zang ================================================================================== import kwant import numpy as np import scipy primitive_vectors = [ (1.0,0.0,0.0), (0.0,1.0,0.0), (0.0,0.0,1.0) ] #coordinates of all atoms # H-1s, C-2s, C-2px, C-3py, C-4pz # index 0 1 2 3 4 basis_sites = [(0.0,0.0,0.0), (0.0,0.0,0.0), (0.0,0.0,0.0), (0.0,0.0,0.0), (0.0,0.0,0.0)] lat = kwant.lattice.Polyatomic(primitive_vectors, basis_sites) syst = kwant.Builder() basis = lat.sublattices #onsite energies, or the diagonal terms of the Hamiltonian matrix syst[basis[0](0.0 , 0.0 , 0.0)] = 1.0 # for H1 at 0.0 0.0 0.0 syst[basis[1](2.0 , 0.0 , 0.0)] = 1.0 # for C1 at 2.0 0.0 0.0 syst[basis[2](2.0 , 0.0 , 0.0)] = 1.0 # for C1 at 2.0 0.0 0.0 syst[basis[3](2.0 , 0.0 , 0.0)] = 1.0 # for C1 at 2.0 0.0 0.0 syst[basis[4](2.0 , 0.0 , 0.0)] = 1.0 # for C1 at 2.0 0.0 0.0 syst[basis[1](8.0 , 0.0 , 0.0)] = 1.0 # for C2 at 8.0 0.0 0.0 syst[basis[2](8.0 , 0.0 , 0.0)] = 1.0 # for C2 at 8.0 0.0 0.0 syst[basis[3](8.0 , 0.0 , 0.0)] = 1.0 # for C2 at 8.0 0.0 0.0 syst[basis[4](8.0 , 0.0 , 0.0)] = 1.0 # for C2 at 8.0 0.0 0.0 #off-site copling, or the off-diagonal terms of the Hamiltonian matrix syst[basis[0](0.0 , 0.0 , 0.0) , basis[1](2.0 , 0.0 , 0.0)] = -2.0 syst[basis[0](0.0 , 0.0 , 0.0) , basis[2](2.0 , 0.0 , 0.0)] = -2.0 syst[basis[0](0.0 , 0.0 , 0.0) , basis[3](2.0 , 0.0 , 0.0)] = -2.0 syst[basis[0](0.0 , 0.0 , 0.0) , basis[4](2.0 , 0.0 , 0.0)] = -2.0 syst[basis[0](0.0 , 0.0 , 0.0) , basis[1](8.0 , 0.0 , 0.0)] = -4.0 syst[basis[0](0.0 , 0.0 , 0.0) , basis[2](8.0 , 0.0 , 0.0)] = -4.0 syst[basis[0](0.0 , 0.0 , 0.0) , basis[3](8.0 , 0.0 , 0.0)] = -4.0 syst[basis[0](0.0 , 0.0 , 0.0) , basis[4](8.0 , 0.0 , 0.0)] = -4.0 kwant.plot(syst) syst = builder.finalized() sparse_mat = syst.hamiltonian_submatrix(sparse=False) print(np.matrix(sparse_mat)) evs = scipy.sparse.linalg.eigs(sparse_mat,3)[0] print(evs.real) ===============================================================================