Skip to content

Oxygen atom reference energy

FHI-aims reference energy calculation for a tricky case

The O atom's valence occupation is 2s2 2p4 ... which means ( Hund's rules for the spin-polarized free atom) that:

  • the s-derived orbitals will be filled
  • the p-derived spin-up orbitals will be filled
  • and there will be one electron left over for the three spin-down p orbitals.

When starting from a spherical free-atom initialization, the only choice we have is to distribute the one electron evenly between the three p orbitals. If the p orbital eigenvalues are exactly equal, this will inevitably lead to exactly 1/3 electron in each of the orbitals. The resulting free atom is still spherical, but alas, the chemists will teach you, not quite natural. There are no split atoms in nature.

This simple problem would not matter much in Kohn-Sham density functional theory, where the desired outcome is the density (and lowest possible total energy), rather than some pesky single-particle eigenvalues. As long as the fractional occupation and spherical atom give the minimum total energy density, why worry?

Well. Alas, it turns out that the spherically symmetric solution for free atoms is usually not the lowest-energy solution. Breaking the symmetry gains energy. While we just dismissed those Kohn_Sham orbitals, that may have been a little quick - one rigorous definition is the minimum total energy after all, however that was obtained.

Within a given density functional (here, the Perdew Wang LDA), we must thus break the symmetry of the free atom. One way to achieve this end is to prevent the existence of fractional occupation numbers altogether. Thus, let's use a very small "smearing" at the Fermi level:

occupation_type gaussian 1-6

... and hope for very small differences between the formally equal 2p eigenvalues, simply due to some very small numerical symmetry breaking.

In practice, this strategy usually works, but since we now have three different, degenerate (non-symmetric) Kohn-Sham eigenvalue solutions, a general mixing strategy may search for a while before settling for one of them.

For the O atom, Perdew Wang LDA, atomic_zora scalar, and the radial grids used in tight and really tight species defaults, here is an example of input and output files that do converge to the proper, non-symmetric solution:

control.in

geometry.in

O-atom.pw-lda.atomic_zora.rm_2.gauss_1e-6.out

Indeed, by looking at the last(!), atomic_zora scalar corrected total energy, we find:

  Spin-up eigenvalues:

  State    Occupation    Eigenvalue [Ha]    Eigenvalue [eV]
      1       1.00000         -18.794530         -511.42517
      2       1.00000          -0.917511          -24.96674
      3       1.00000          -0.402694          -10.95785
      4       1.00000          -0.402693          -10.95784
      5       1.00000          -0.336333           -9.15210
      6       0.00000          -0.000143           -0.00388
      7       0.00000           0.079015            2.15011

  Spin-down eigenvalues:

  State    Occupation    Eigenvalue [Ha]    Eigenvalue [eV]
      1       1.00000         -18.740306         -509.94968
      2       1.00000          -0.796765          -21.68107
      3       1.00000          -0.273406           -7.43976
      4       0.00000          -0.262342           -7.13868
      5       0.00000          -0.262341           -7.13867
      6       0.00000           0.004767            0.12972
      7       0.00000           0.085184            2.31798

The spin-up channel shows the three filled, p-derived levels. In the spin-down channel, only one of the p-derived levels is filled, but with a single full occupation number (1.00000).

The total energy is found right below:

  | Total energy uncorrected      :         -0.203033254607510E+04 eV
  | Total energy corrected        :         -0.203033254607510E+04 eV  <-- do not rely on this value for anything but (periodic) metals
  | Electronic free energy        :         -0.203033254607510E+04 eV

In case of doubt, it really is the first line that matters. Only fractional occupation numbers would lead to a difference between the three values, and even then, the free energy expressions or T->0 extrapolations for the free electron gas (the lower two lines) really do not apply for a single, isolated atom.

Contrast this with the spherical solution, which is not the ground state

This is the spherically symmetric solution for the free O atom, Perdew Wang LDA, atomic_zora scalar, which is '''NOT''' the lowest energy solution - '''BUT''' this is a free-atom solution that is also self-consistent, and that you will find with almost all standard settings if you are not really careful with breaking the symmetry.

The following values were obtained with the exact same input files as for the lowest-energy, symmetry broken solution, except we set:

occupation_type gaussian 1e-5

The resulting occupations of the spin channels are:

  Spin-up eigenvalues:

  State    Occupation    Eigenvalue [Ha]    Eigenvalue [eV]
      1       1.00000         -18.795266         -511.44521
      2       1.00000          -0.917001          -24.95288
      3       1.00000          -0.380566          -10.35574
      4       1.00000          -0.380566          -10.35574
      5       1.00000          -0.380566          -10.35574
      6       0.00000          -0.000620           -0.01687
      7       0.00000           0.079176            2.15449

  Spin-down eigenvalues:

  State    Occupation    Eigenvalue [Ha]    Eigenvalue [eV]
      1       1.00000         -18.741865         -509.99209
      2       1.00000          -0.803353          -21.86036
      3       0.33350          -0.272197           -7.40687
      4       0.33346          -0.272197           -7.40687
      5       0.33304          -0.272197           -7.40687
      6       0.00000           0.003743            0.10185
      7       0.00000           0.086887            2.36431

and the total energy is:

  | Total energy uncorrected      :         -0.203022743299952E+04 eV
  | Total energy corrected        :         -0.203022743685608E+04 eV  <-- do not rely on this value for anything but (periodic) metals
  | Electronic free energy        :         -0.203022744071263E+04 eV

This is a perfectly innocent-looking, spherically symmetric O atom, except the total energy is about 0.11 eV higher than the lowest-energy, symmetry broken one.

Using the PBE GGA, the difference becomes much greater, leading to serious discrepancies in published cohesive energies for oxygen-containing materials.