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

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

• Bit off-topic, but I personally prefer using gist.github.com for this kind of thing instead of having to keep a repo around forever. Jul 2 '21 at 16:13
• This is a standard result in HF theory, see e.g. section 3.3.1 in Szabo and Ostlund. The issue is that summing the single particle energies double counts the electron-electron interactions. Jul 2 '21 at 16:43
• For the problem at hand, the log file of the calculation would fit in its entirety into this post. However, it is completely unnecessary as in this question Gaussian is merely an example. You'd still have the same problem if you'd use Orca, Psi4, Gamess, … , Pen and Paper. In the g16 case, the input file (and version) should be enough to reproduce the calculation with miniscule differences; scientific rigor should still stand. But, since you have posted it: Do not use NBO 3.1! That is now deprecated for at least a decade, we're at 7, now. It really has issues with theory. Don't use it. Please. Jul 2 '21 at 22:27

## 1 Answer

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.

• As well as the maths you might consider a physical picture - the orbital energy of electron 1 is from electron 1 moving in the potential of all the other electrons, so it includes the interaction of electron 1 and electron 2. Similarly electron 2's orbital energy includes the interaction with all other electrons, including electron 1. Thus if you sum the orbital energies you include each e-e electron twice. Thus you need to half it to get the right energy. Jul 2 '21 at 16:56