Home
Calculating the total energy of free atoms sounds like a simple thing to do, but in density-functional theory, it is not as simple as it may seem.
Read on on the present page for the underlying physical problem and also for an example of a control.in file that defines the problem more clearly and thus solves problems of numerical instability.
On free atoms
Most unsuspecting users assume that free atoms are inherently spherical in their ground state and therefore, there is only one unique ground state.
THIS IS NOT THE CASE.
You may easily check this for the case of an O atom, using the PBE functional and apply what is written down below. You will find that the energy difference between a spherically symmetric state and the actual, symmetry breaking ground state amounts to 0.7 eV (written from memory). That's a lot. How come?
In DFT and/or Hartree-Fock theory, we only have a single Slater determinant. Unless all orbitals of an angular momentum shell are fully occupied, we would have to use fractional occupation numbers to define the density:
where f_orbital between 0 and 1 in the spin-polarized case. In particular, in the case of O, one of the spin channels would have to be occupied with 1/3, 1/3, 1/3 occupation in the m=-1,0,1 orbitals to create a spherically symmetric density.
This may sound ok from a formal point of view (in the basic foundations of DFT, fractional orbitals are justified as statistical mixtures of ensembles of atoms in different pure states). However, actual density functional approximation do NOT necessarily produce the average of pure states when fractional occupations are used. If you want something approaching a pure state, you would need (in the simplest case) a single determinant with only completely filled or completely empty orbitals. For something like O, this would inherently break the spherical symmetry, as only one of the p orbitals will be occupied.
This is the formal aspect. From a practical DFT point of view, all that matters is which occupation has the lowest total energy. And in practice, the lowest total energy is often given by the "pure state" occupation, i.e., f=0 or 1. So we need to break the spherical symmetry of an atom.
But this is not straightforward. You will note that since there are three p orbitals for O, there must be more than one lowest-energy state - i.e., more than one lowest-energy solution. Any electronic structure code, during its s.c.f. optimization, needs to pick one of those solutions to fall into as a stationary point of the s.c.f. problem. In other words, the s.c.f. problem, which is a nonlinear optimization problem, has more than one stationary point.
AND THIS CAN PRODUCE PROBLEMS OF NUMERICAL UNIQUENESS AND/OR STABILITY, ESPECIALLY IN MORE COMPLICATED d OR f SHELL systems.
The basic problem is that we need to start somewhere. Just like in geometry optimization, if we are searching for a local minimum, which minimum we find depends on where we start. There might be more than one local minimum in the total energy surface, one or more of which are the total energy minimum.
FHI-aims, by definition, starts with a superposition of free atom densities created from spherical free atoms with fractional occupation numbers. You will note that this may, in fact, not be a minimum of the total energy surface, but because of its symmetry, it would have to be a saddle point or a maximum of the total energy.
As an important consequence, starting from such a point WILL mean that we can fall into different local minima depending on which direction of the slope our calculation decides to go. Or, worse, since we use an optimization scheme to find a single local minimum, the density mixing scheme might be confused by the fact that more than one minimum is present and not converge at all. Even more interestingly, we use fractional occupation numbers (occupation_type) for numerical reasons in our calculation. If we introduce a large broadening, we might find fractional occupation numbers and a non-pure states. Only if we mandate a very low broadening will we safely find integer occupation numbers, but then, we face the problem of multiple allowed solutions and potential numerical instabilities.
Because, by default, we start from a saddle point or maximum, we may even fall into different local minima starting from tiny numerical changes, such as the number of MPI tasks, which changes the order of summation in the integration of matrix elements. At double precision, this leads to initial discrepancies of 1e-12 or so at some point in the s.c.f. mixing cycle, but this difference is large enough for the Pulay mixer to pick it up and determine which specific solution it will end up in after the s.c.f. cycle. But the solution(s) found are all valid from an s.c.f. point of view. It's just that there is more than one of them.
Why do we not just start from a symmetry broken state?
Because this would mean we would have to pick FOR YOU which symmetry broken state you want.
In a problem with multiple minima, we would force you to use one of them without telling you, and (a) we might not get it right or (b) you might not want the specific state we picked for you because you were interested in something else.
Either way, this is a user choice - if we were to legislate which electronic structure result you get without telling you, we would be doing something wrong.
The real way out, like in all problems of multiple minima, is that you would need to pick which minimum you want - or search through them all, systematically. This is generally difficult for any global optimization problem, but for free atoms, there is a way.
In FHI-aims, you can constrain certain occupation numbers (really, Mulliken occupations) for certain valence orbitals by hand, in control.in , using the force_occupation_basis keyword. Here is a specific example enforcing a specific occupation of 4p valence orbitals for the Se atom
force_occupation_basis 1 2 atomic 4 1 1 0. 20
force_occupation_basis 1 2 atomic 4 1 0 0. 20
Please read the manual regarding this keyword. Don't just copy these lines into every other free atom calculation - again - they enforce ONE specific shell occupation for one specific type of atom, Se. The chosen syntax for other atoms WILL BE different.
In short, we legislate Mulliken occupations of 0. in the first 20 Kohn Sham states for atom 1, spin channel 2, atomic orbital n=4 l=1 and m=1 (first line) and m=0 (second line).
This means that the electron in the 4 p spin channel 2 will have to be found completely in the m=-1 orbital.
The choice of the "20" states is also system specific. From the point of view of the mathematical constraint, we need to project the 4 p 0 and 4 p 1 orbitals out of the occupied Kohn-Sham states but not out of the unoccupied Kohn Sham states. (The basis functions define a vector space and if you calculate eigenvectors and eigenstates, each part of the underlying basis needs to go somewhere or else the dimension of your eigenspace is not the same as that of the underlying basis space.)
In longer notes, the control.in file for this particular instance of Se was this:
xc pw-lda
spin collinear
default_initial_moment hund
relativistic atomic_zora scalar
#
mixer pulay
charge_mix_param 0.05
sc_accuracy_rho 1e-5
#
occupation_type gaussian 1e-6
#
# Enforce occupation of specific atomic orbital basis functions
#
# Usage: force_occupation_basis i_atom spin basis_type basis_n basis_l basis_m occ_number max_KS_state
#
# Purpose: Flag originally programmed to compute core-hole spectroscopy simulations.
# In practice, it constrains the occupation of a specific energy level of an specific atom,
# being also able to break the symmetry of an atom.
#
# See manual for more information
#
force_occupation_basis 1 2 atomic 4 1 1 0. 20
force_occupation_basis 1 2 atomic 4 1 0 0. 20
#
################################################################################
#
# FHI-aims code project
# Volker Blum, Fritz Haber Institute Berlin, 2009
#
# Suggested "tight" defaults for Se atom (to be pasted into control.in file)
#
################################################################################
species Se
# global species definitions
nucleus 34
mass 78.96
#
l_hartree 6
#
cut_pot 8.0 2.0 1.0
basis_dep_cutoff 0.0
#
radial_base 55 7.0
radial_multiplier 2
angular_grids specified
division 0.0830 110
division 0.1357 194
division 0.7377 302
division 0.8610 434
# division 0.9640 590
# division 1.0773 770
# division 2.5539 974
# outer_grid 974
outer_grid 434
################################################################################
#
# Definition of "minimal" basis
#
################################################################################
# valence basis states
valence 4 s 2.
valence 4 p 4.
valence 3 d 10.
# ion occupancy
ion_occ 4 s 1.
ion_occ 4 p 3.
ion_occ 3 d 10.
################################################################################
#
# Suggested additional basis functions. For production calculations,
# uncomment them one after another (the most important basis functions are
# listed first).
#
# Constructed for dimers: 1.85 A, 2.15 A, 2.50 A, 3.00 A, 4.00 A
#
################################################################################
# "First tier" - improvements: -336.21 meV to -36.85 meV
hydro 3 d 4.3
hydro 2 p 1.6
hydro 4 f 7.2
ionic 4 s auto
# "Second tier" - improvements: -14.39 meV to -1.49 meV
hydro 5 g 10.4
hydro 4 p 4.5
hydro 4 d 6.2
hydro 4 f 11.2
hydro 1 s 0.65
hydro 6 h 15.2
# Third tier - improvements: -0.46 meV and below
hydro 5 g 14.4
ionic 4 d auto
hydro 4 f 22.4
hydro 5 f 7.4
hydro 5 p 8
hydro 5 s 9.4
# Fourth tier - improvements: -0.12 meV and below
# hydro 5 d 11.6
# hydro 6 h 18
# hydro 4 p 4
# hydro 5 f 16
# hydro 4 s 3.9