# Why is my hydrogen energy not equal to -0.5 hartrees?

I'm trying to calculate the HOMO of various molecules in water solution. I first tried to see how PySCF calculates the HOMO of hydrogen. I ran the following code:

mol_h = pyscf.gto.M(atom = 'H 0 0 0', spin = 1, basis = 'ccpvdz')
rks_h = mol_h.ROKS(xc = 'b3lyp')
rks_h.kernel()
rks_h.analyze(verbose=5)

And I got the following output:

**** SCF Summaries ****
Total Energy =                          -0.501257936931332
Nuclear Repulsion Energy =               0.000000000000000
One-electron Energy =                   -0.499271486999867
Two-electron Coulomb Energy =            0.311475159768718
DFT Exchange-Correlation Energy =       -0.313461609700183
**** MO energy ****
Roothaan           | alpha              | beta
MO #1   energy= -0.148854429623937 | -0.319481818297586 | 0.02178298730584   occ= 1
MO #2   energy= 0.520550055619181  | 0.4180250336286    | 0.623065049353633  occ= 0
MO #3   energy= 1.3343779636294    | 1.2099070476478    | 1.458848879611     occ= 0
MO #4   energy= 1.33437796363348   | 1.2099070476478    | 1.45884887961915   occ= 0
MO #5   energy= 1.33437796363348   | 1.2099070476478    | 1.45884887961915   occ= 0
** MO coefficients (expansion on meta-Lowdin AOs) **
#1        #2        #3        #4        #5
0 H 1s           0.99722  -0.07457  -0.00000  -0.00000   0.00000
0 H 2s           0.07457   0.99722  -0.00000  -0.00000   0.00000
0 H 2px          0.00000   0.00000   0.00005   0.82979  -0.55808
0 H 2py          0.00000   0.00000  -0.00001   0.55808   0.82979
0 H 2pz          0.00000   0.00000   1.00000  -0.00003   0.00004
** Mulliken pop on meta-lowdin orthogonal AOs  **
** Mulliken pop  **
pop of  0 H 1s            0.99444
pop of  0 H 2s            0.00556
pop of  0 H 2px           0.00000
pop of  0 H 2py           0.00000
pop of  0 H 2pz           0.00000
** Mulliken atomic charges  **
charge of    0H =      0.00000
Dipole moment(X, Y, Z, Debye): -0.00000, -0.00000, -0.00000

((array([9.94438767e-01, 5.56123333e-03, 2.85005342e-34, 2.87485689e-35,
1.36070362e-34]),
array([3.33066907e-16])),
array([-5.71059869e-17, -1.81369083e-17, -3.94581675e-17]))

According to PySCF calculation, it seems the HOMO should be -0.14885. However, isn't the correct HOMO value of hydrogen -0.5 (according to this answer)?

Yes, the exact orbital energy of the hydrogen atom is $$-1/2\ E_h$$. The problems in the calculation are as follows.

First, you are running ROKS. ROKS has many definitions in the literature; PySCF appears to use the code for ROHF (which likewise has many, many definitions in the literature; see Kreb's 141-page review in Comput. Phys. Commun. 116, 137 (1999), for example for ways to handle the off-diagonal Lagrange multipliers as well as Table I in J. Chem. Phys. 140, 014102 (2014) for various definitions of the ROHF Hamiltonian) in combination with the Kohn-Sham Fock matrix instead of the Hartree-Fock Fock matrix. As far as I remember, it is generally a bad idea to do this in density functional theory, since in general it will break one of the fundamental theorems of density functional theory, which is the density cusp on the nuclei.

ROKS with the B3LYP functional you used yields an orbital energy $$-0.148854 E_h$$. If you switch to UHF, you get a much better value $$-0.319475 E_h$$. Finally, the reason why this is still far from the exact value is self-interaction error: the Coulomb and exchange terms don't cancel out perfectly in DFT, which becomes especially clear in the case of one-electron systems whose exact solution is easy to find, such as in the case of hydrogen; we recently published a benchmark on this topic in J. Chem. Phys. 157, 174113 (2022).

Hartree-Fock is exact for one-electron systems, and either RHF or UHF works fine and you will find the orbital energy $$-0.499278 E_h$$ in the cc-pVDZ basis you used. The remaining error is due to basis set truncation error: the $$1s$$ orbital is not exactly reproduced by the cc-pVDZ basis set.

Note that Hartree-Fock still gets all the other orbitals completely wrong. We know that the exact energy for primary quantum number $$n$$ is $$E_n = - \frac 1 {2n^2}$$ which means that one should find the Rydberg spectrum with energies $$-0.125000E_h$$, $$-0.055556E_h$$, $$-0.031250E_h$$, etc. Instead, the Hartree-Fock calculation gives the first unoccupied orbital energy as $$0.181933 E_h$$. This is because the potential for the unoccupied orbitals given by Hartree-Fock (and Kohn-Sham DFT, too!) is wrong: it decays exponentially, instead of $$-1/r$$ like the exact potential is known to behave.

The exact non-relativistic ground state electronic energy of a hydrogen-like atom, assuming that the nucleus is "clamped", point-sized, and with an infinitely large mass compared to the mass of an electron, and ignoring any fine, hyperfine, or finer than hyperfine effects, and in a vacuum containing nothing else, is -0.5 Eh according to Schroedinger's 1926 equation.

You have approximately found the lowest eigenvalue of the Schroedinger equation for the above-mentioned hydrogen atom, using the cc-pVDZ basis set and the B3LYP functional (and the ROKS mean-field method).

When I do HF calculations for larger basis sets for an H atom, I get the following energies, and this shows that if your basis set were big enough, and you did an HF calculation, you would get -0.5 Eh:

aug-cc-pV6Z:    -0.499999 276396 663  (PySCF, 127 functions), default settings
aug-cc-pV7Z:    -0.499999 752574 58   (PySCF, 189 functions), default settings
aug-cc-pV8Z:    -0.499999 966893 035  (PySCF, 268 functions), default settings
aug-cc-pV9Z:    -0.499999 977357 467  (PySCF, 367 functions), default settings
aug-cc-pV10Z:   -0.499999 980721 463  (PySCF, 485 functions), default settings

You can see this data and more here. You can also see the original input/output files here.