old_contents: import numpy as np
  from itertools import product
  from functools import reduce
  
  def mp2(hamiltonian, orbs, orbe, nocc, nvirt, verbose=False):
      ints = hamiltonian.i2
      moints = ints.transform_mp2(orbs, nocc)
      Emp2 = 0
      for a,b in product(range(nocc), repeat=2):
          Eab = 0
          for r,s in product(range(nvirt), repeat=2):
              arbs, asbr = moints[a,r,b,s], moints[a,s,b,r]
              Eab += arbs*(2*arbs-asbr)/(orbe[a]+orbe[b]-orbe[nocc+r]-orbe[nocc+s])
          if verbose: print("MP2 pair energy for %d,%d: %f" % (a, b, Eab))
          Emp2 += Eab
      if verbose: print("MP2 energy: %f" % Eab)
      return Emp2
  
  if __name__ == '__main__': test_mp2()
  
new_contents: import numpy as np
  from itertools import product
  from functools import reduce
  
  
  def mp2(hamiltonian, orbs, orbe, nocc, nvirt, verbose=False):
      ints = hamiltonian.i2
      moints = ints.transform_mp2(orbs, nocc)
      Evirt, Eocc = orbe[nocc:], orbe[:nocc]
  
      denominator = 1/(Eocc.reshape(-1, 1, 1, 1) - Evirt.reshape(1, -1, 1, 1) +
                       Eocc.reshape(1, 1, -1, 1) - Evirt.reshape(1, 1, 1, -1))
  
      MP2corr_OS = np.einsum('iajb,iajb,iajb->', moints, moints, denominator)
      MP2corr_SS = np.einsum('iajb,iajb,iajb->', moints - moints.swapaxes(1,3), moints, denominator)
      return MP2corr_OS + MP2corr_SS
  
  
  if __name__ == '__main__': test_mp2()
  
subject: rewrite using numpy