Dear Sir, I have two questions. (1) I want to study the intralayer intrinsic spin-orbit coupling in bilayer graphene. For the monolayer case, it is easy to put the nnn hopping by, sys[graphene.neighbors(2)] = spin_orbit But for the bilayer case, as there are two planes and the intralayer coupling which I want to calculate should be in those two planes, Can I use the above code to get the effect. If not, then how can I put the effect in bilayer system. (2) When I am calculating the Hamiltonian submatrix between two sites to see the full matrix without spin-orbit ( only for tight-binding in bilayer where t=1 and interlayer coupling ,t_prep= 0.2) it is showing [[7.+0.j 0.+0.j 0.+0.j 0.+0.j] [0.+0.j 7.+0.j 0.+0.j 0.+0.j] [0.+0.j 0.+0.j 7.+0.j 0.+0.j] [0.+0.j 0.+0.j 0.+0.j 7.+0.j]] where 7 is for onsite element (sys[graphene.shape(rect, (0, 0, 0))] = 7 * sigma_0). Thanks in advance. sincerely, Priyanka Sinha