I was doing some simple calculations and notice something was off with energy levels.
The systems were just single atoms (H,C,N,O,F) computed with a custom basis set (decontracted) and ROHF, symmetry C1 (because I still ignore group theory). I wanted to see the atomic energy levels!
For carbon, 2s2 2p2 (multiplicity 3) I got these (in Hartree):
1s = -11.324757
2s = -0.703903
2p = -0.126023 for the singly occupied and 0.061888 for the unoccupied one.
The energies were printed from the wavefunction.epsilon_a() vector.
Then I compared to an older database (https://doi.org/10.1006/adnd.1993.1003 or http://www.tcm.phy.cam.ac.uk/~mdt26/RHF.TABLES.txt and also http://www.few.vu.nl/~visscher/FiniteNuclei/Table3.html)
where the 2p is listed as -0.433341. That one reference was done with Slaters and not GTOs, and the other was from Dirac-Fock… but i guess they should not be too different.
So I fired up Orca and computed C with the same basis as in Psi4 and got -0.43151 for the 2p, similar to that old reference.
The 1s and 2s are quite similar… nothing that would make me suspicious. But the 2p energies are off quite a lot, and for all the organics I tested.
Using the coefficients I calculated the expectation value or r, r^2 and 1/r for all the obtained orbitals and they are in quite good agreement with the reference tables.
Did I misunderstand the meaning of wavefunction.epsilon_a() or there is some other issue?
My best guess here is differing definitions of the diagonal blocks of the ROHF Fock matrix. (Only the off-diagonal blocks are uniquely defined. Did you use ROHF or UHF in your Orca calculation? if you choose UHF in PSI4, your 2p orbital energies will be closer to what you’ve expected.
I did ROHF in both Orca and PSI4. Interestingly, Orca gave three orbitals obtained from the basis p ones, with identical energies and a weird fractional occupation number.
I now tried UHF with PSI4 and got the right answer… -0.43 hartree for the 2p of carbon, but I am puzzled. Could you explain in simple words what is going on?
Should I never trust MO energies from ROHF?
Thanks a lot!
ROHF MO energies aren’t uniquely defined quantities, unlike purely canonical RHF or UHF orbital energies. In the latter, you have one set of Hartree-Fock conditions (that first-order mixings between the occupied and unoccupied orbitals must be zero), and the canonicalization condition that the Fock matrix be diagonal is simply a more restrictive version thereof. In ROHF, you have three conditions that must be simultaneously satisfied: mixings must be zero between doubly occupied and singly occupied orbitals, between doubly occupied and unoccupied orbitals, and between singly occupied and unoccupied orbitals. Different quantum chemistry programs typically choose a definition of the diagonal blocks of the Fock matrix that is computationally convenient. (PSI, for example, at one point chose diagonal blocks that were an average between alpha- and beta-spin Fock operators because it led to simpler forms of the orbital Z-vector equations in configuration interaction gradient theory, but that’s a long story.)
Should you “never trust MO energies from ROHF”? That depends on your purpose. Given that MO energies aren’t observable quantities anyway, I wouldn’t tend to “trust” them for pretty much anything. However, if you want to compare them to ionization energies or electron affinities, for example, just be aware that not all definitions for ROHF MOs yield expressions that are comparable to the corresponding Koopmans-like equations. This depends on the quantum chemistry code you’re using.
This article is probably more technical than you’re looking for, but may help you with your analysis and understanding of ROHF orbital energies.
Thanks a lot to everyone for the explanation.