docs: def unperturbed_hamiltonian(states):
      r"""Return the unperturbed atomic hamiltonian for given states.
  
      We calcualte the atomic hamiltonian in the basis of the ground states of \
      rubidium 87 (in GHz).
      >>> g = State("Rb", 87, 5, 0, 1/Integer(2))
      >>> magnetic_states = make_list_of_states([g], "magnetic")
      >>> print(np.diag(unperturbed_hamiltonian(magnetic_states))/hbar/2/pi*1e-9)
      [-4.2717+0.j -4.2717+0.j -4.2717+0.j  2.563 +0.j  2.563 +0.j  2.563 +0.j
        2.563 +0.j  2.563 +0.j]
  
      """
      Ne = len(states)
      H0 = np.zeros((Ne, Ne), complex)
      for i in range(Ne):
          H0[i, i] = hbar*states[i].omega
      return H0