docs: def spin_hamiltonian(self, h, l):
          """Constructs the single site spin Hamiltonian"""
          h_spin  = np.einsum('i,ijk', h[1], self.oper['O'])
          h_spin += np.einsum('i,ijk', h[0], self.oper['O_d'])
          h_spin += np.einsum('i,ijk', l, self.oper['Sz+1/2'])
  
          h_spin += self.oper['Hint']
          return h_spin