Why don't the MO energies sum up to Hartree-Fock SCF value?

Question

I would like to calculate Hartree-Fock energy of $\ce{H3^+}$ ion from orbital energies. So I performed the following calculation in Gaussian 16:

# hf/sto-3g
1 1
H                     0.        0.        0.
H                     0.        0.        0.85261
H                     0.        0.       -0.85261

Then I took these two energy values from the log:

 nuclear repulsion energy         1.5516354770 Hartrees
...
Alpha  occ. eigenvalues --   -1.07569

Unexpectedly, the sum 1.5516354770 - 2 * 1.07569 = -0.599765 Hartree is not even close to "SCF Done:" value computed by Gaussian:

SCF Done:  E(RHF) =  -1.20555771776 

Since I didn't explicitly require any corrections or solvation model, I wonder what contribution to SCF energy did I miss?

Input file and log are uploaded on https://github.com/nkrivoshchapov/h3plus_calculation

Answer 1

The key misconception you have is that the sum of the orbital energies is equal to the electronic energy. Compare the formulas for the two (taken from Chapter 3 of Szabo and Ostlund's Modern Quantum Chemistry)

$$E_0=\sum_i^Nh_{ii}+\frac{1}{2}\sum_{ij}^N\langle ij||ij\rangle$$ $$\sum_i^N\epsilon_i=\sum_i^Nh_{ii}+\sum_{ij}^N\langle ij||ij\rangle$$

If you simply sum up the orbital energies, you are double counting the electron-electron repulsion interactions.

For a closed-shell system, we can write the electronic energy in terms of the orbital energies, but we need additional terms: $$E_0=\sum_i^{N/2}(h_{ii}+\epsilon_i)$$ To get Gaussian to print the core Hamiltonian (along with its components), you need to include the keyword iop(3/33=1) (or iop(3/33=5) if you also wanted to see the two electron integrals for some reason). You will also need pop=full to print out the MO coefficients so you can transform the core Hamiltonian to MO basis.

For convenience, I'll just copy the results here:

MO coefficients

0.58017   0.00000   1.42841
0.31436   0.78557  -0.94873
0.31436  -0.78557  -0.94873

Core Hamiltonian

-0.163989E+01 -0.107667E+01 -0.107667E+01    
-0.107667E+01 -0.136306E+01 -0.377849E+00    
-0.107667E+01 -0.377849E+00 -0.136306E+01   

When this is done, we can now evaluate the Hartree-Fock Energy, which I did using Python:

>>> import numpy as np
>>> C = np.loadtxt("C.txt")
>>> H = np.loadtxt("H.txt")
>>> H_MO = C.T @ H @ C
>>> H_MO[0,0]
-1.6815253284170897
>>> e_MO = -1.07569
>>> E_NR = 1.5516354770
>>> E_HF = H_MO[0,0] + e_MO + E_NR
>>> E_HF
-1.2055798514170895

This is in much better agreement and the only reason for the deviation is that the MO energies/coefficients and core Hamiltonian aren't printed to sufficient precision.

You can also get the Hartree-Fock energy by summing the individual components of the energy printed after the SCF, along with nuclear repulsion energy (requires #p in the input line)

 SCF Done:  E(RHF) =  -1.20555771776     A.U. after    4 cycles
            NFock=  4  Conv=0.11D-09     -V/T= 2.1301
 KE= 1.066737227463D+00 PE=-4.429751734382D+00 EE= 6.058213121791D-01

Where KE/PE/EE is the kinetic/potential/electronic energy.