Hydrogen atom: orbital energy is different from the total electronic energy
Vsevolod Dergachev
Dear DIRAC community,
I am doing atomic calculations in DIRAC. As a test for the one-electron system, I have calculated the ground state energy of the hydrogen atom using the SCF module of DIRAC. I have found that the 1s1/2 orbital energy is different from the total electronic energy, ~ -0.5 Hartree. I am a bit confused since there are no other contributions to the energy, like nuclear-nuclear repulsion or electron correlation. For comparison, I ran SCF for H atom in GAMESS and got the total electronic energy equal to the energy of the 1s orbital, as I would expect. Could you please elucidate this situation in DIRAC to me?
My only guess is that this situation is specific only to one-electron calculations. To test it, I have run SCF for the Helium atom in DIRAC and GAMESS and got the match up to four digits in respective 1s orbital energies.
I am attaching the output files for both H and He atoms.
Thank you very much,
Seva Dergachev
Trond SAUE
Dear Seva,
in your DIRAC calculation the total energy is reported as -0.50000640133863117 au, which is close to the analytic value. There are (at least) two sources of difference:
1) You are using a finite basis set, but it is high-quality (4Z), so the error from this is very small.
2) By default, DIRAC uses a finite nucleus model, not a point charge as used in the analytically solvable case. A bit down in the output you find
Nuclear model: Gaussian charge distribution.
The orbital energy is
* Open shell #1, f = 0.5000 -0.343753242239 ( 2)
and indeed surprising. For open-shell cases, as the hydrogen
atom, DIRAC uses average-of-configuration HF, as explained here
http://www.diracprogram.org/doc/release-19/tutorials/open_shell_scf/aoc.html
We usually think of HF calculations as cycles of Fock matrix construction and diagonalization, and this is also how things are done in DIRAC. However, if the energy is parametrized in terms of orbital rotation
http://www.diracprogram.org/doc/release-19/tutorials/open_shell_scf/aoc.html#orbital-rotations
which has the advantage of avoiding orthogonality constraints, then one can show that only rotations between different orbital classes are redundant and that the gradient with respect to such orbital rotations defines the OFF-diagonal blocks (between orbital classes) of the Fock matrix. For DIAGONAL blocks there is some freedom, hence allowing for convergence tricks such as level shifts. One way to fix the diagonal blocks is to insist that the orbital energies should obey Koopman's theorem, also for open-shell cases, as discussed here
However, the resulting definition is not necessary the one that provides the best convergence. A bit down in your output you find
* Scaling of active-active block correction to open shell Fock operator 0.500000 to improve convergence (default value). The final open-shell orbital energies are recalculated with 1.0 scaling, such that all occupied orbital energies correspond to Koopmans' theorem ionization energies.
The final statement is in fact wrong; orbital energies are presently not corrected at the end of calculations, so this shall be fixed. If you really want to have Koopmans' theorem obeyed, you should set under *SCF
.OPENFAC
1.0
Present default, as indicated in the output, is 0.5.
All the best,
Trond
Trond Saue
Laboratoire de Chimie et Physique Quantiques
UMR 5626 CNRS --- Université Toulouse III-Paul Sabatier
118 route de Narbonne, F-31062 Toulouse, France
Phone : +33/561556361 Fax: +33/561556065
Mail : trond...@irsamc.ups-tlse.fr
Web : http://dirac.ups-tlse.fr/saue
DIRAC : http://www.diracprogram.org/
ESQC : http://www.esqc.org/
Book: Principles and Practices of Molecular Properties: Theory, Modeling, and Simulations