Dear Colleagues, Imagine, I have two different 2D lattices in a heterostructure and want to set up interlayer hoppings between the close-by atoms. I solve this problem in the way presented below. This works, but to my surprise it works so slowly that it became a bottle-neck of my whole computation. Is there a way to improve the code below? (I suspect, it may be the tag comparison that slows down the thing, but I am not sure) def connect(sys,crit): for site1, site2 in it.product(sys.sites(), repeat=2): if crit(site1, site2): yield (site1, site2) def critInterlayer(site1,site2): ## gamma1 hopping integral diff = site1.pos-site2.pos [x,y] = diff if (abs(x)+abs(y)<0.1) and (site1.tag!=site2.tag): return True else: return False sys[connect(sys,critInterlayer)]=hopping Best wishes, Sergey