9
$\begingroup$

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

$\endgroup$
3
  • 3
    $\begingroup$ 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. $\endgroup$
    – orthocresol
    Jul 2 at 16:13
  • 3
    $\begingroup$ 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. $\endgroup$
    – Ian Bush
    Jul 2 at 16:43
  • 1
    $\begingroup$ 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. $\endgroup$ Jul 2 at 22:27
15
$\begingroup$

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.

$\endgroup$
1
  • 5
    $\begingroup$ 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. $\endgroup$
    – Ian Bush
    Jul 2 at 16:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.