2.3 A very quick guide to ensuring numerical convergence with FHI-aims

The key physical choice to ensure converged results in FHI-aims is the list of radial functions (and their angular momenta) that are used. Beyond the minimal basis of free-atom like radial functions, we always recommend adding at least a single set of further radial functions that are optimized to describe a chemical bond efficiently. These basis functions can be found as a list (line by line) at the end of each species defaults file. For Oxygen / light, the list reads like this:

#  "First tier" - improvements: -699.05 meV to -159.38 meV
     hydro 2 p 1.8
     hydro 3 d 7.6
     hydro 3 s 6.4
#  "Second tier" - improvements: -49.91 meV to -5.39 meV
#     hydro 4 f 11.6
#     hydro 3 p 6.2
#     hydro 3 d 5.6
#     hydro 5 g 17.6
#     hydro 1 s 0.75
#  "Third tier" - improvements: -2.83 meV to -0.50 meV
#     ionic 2 p auto
#     hydro 4 f 10.8
#     hydro 4 d 4.7
#     hydro 2 s 6.8
 [...]

Obviously, only a single set of radial functions (one for each angular momentum $s$, $p$, $d$) is active (not commented!) beyond the minimal basis. Since the minimal basis already contains one additional valence $s$ and $p$ function, this choice is often called “double numeric plus polarization” basis set in the literature (where $d$ is a so-called polarization function as it does not appear as a valence angular momentum of the free atom). We call this level “tier 1”.

In order to increase the accuracy of the basis, further radial functions may be added, simply by uncommenting more lines in order! We recommend to normally proceed in order of full “tiers”, not function by function, but adding specific individual functions on their own can sometimes capture the essence of a problem at lower cost. For example, tier 2 may be added by uncommenting:

#  "First tier" - improvements: -699.05 meV to -159.38 meV
     hydro 2 p 1.8
     hydro 3 d 7.6
     hydro 3 s 6.4
#  "Second tier" - improvements: -49.91 meV to -5.39 meV
     hydro 4 f 11.6
     hydro 3 p 6.2
     hydro 3 d 5.6
     hydro 5 g 17.6
     hydro 1 s 0.75
#  "Third tier" - improvements: -2.83 meV to -0.50 meV
#     ionic 2 p auto
#     hydro 4 f 10.8
#     hydro 4 d 4.7
#     hydro 2 s 6.8
 [...]

tier 2 is the default choice of our tight settings for Oxygen, but may be very expensive for hybrid functionals. In such cases, consider the intermediate settings for a more economical alternative.

Beyond selecting appropriate radial functions, determining the correct confinement radius that all basis functions will experience is a critical step. Ensuring that each radial function goes to zero in a controlled way beyond a certain, given value is critical for numerical efficiency, but on the other hand, you do not want to reduce this confinement radius too much in order to preserve the accuracy of your basis set.

By default, the confinement radius of each potential is specified by the following line:

    cut_pot             3.5  1.5  1.0

This means (see also the CPC publication on FHI-aims, Ref. [36]) that each radial function is constructed with a confinement potential that begins 3.5 Å away from the nucleus, and smoothly pushes the radial function to zero over a width of 1.5 Å. The full extent of each radial function is thus 5 Å.

Of course, this setting is chosen to give good total energy accuracy at the light level, but the convergence of the confinement potential must still be tested, especially in situations where a strong confinement might be unphysical. Such questions include:

• •

  Accurate free atom calculations for reference purposes: choose 8 Å or higher for the onset of the confinement, or something similarly high—for a single free atom, the CPU time will not matter, and you will get all the tails of your radial functions right without much thinking.

• •

  Surfaces— e.g., low electron densities above the surface for STM simulations must not be abbreviated by the onset of the confinement potential—even if this no longer affects the total energy.

• •

  Neutral alkali atoms, or any negatively charged ions. Those are tricky—the outermost electron shell may decay very slowly to zero with distance, requiring explicit convergence tests.

As the corresponding tight setting, we use:

    cut_pot             4.0  2.0  1.0

Although the modification does not seem large, CPU times for periodic systems are significantly affected by this change of the full extent of each radial function from 5 Å to 6 Å. For example, in a densely packed solid, the density of basis functions per volume increases as $R^3$ with the full extent of each radial function, and thus the time to set up the Hamiltonian matrix should increase as $R^6$. Very often, the effect on the total energy is completely negligible, but again, explicit convergence tests are always possible to make sure.

Finally, there is the line

    basis_dep_cutoff    1e-4

If this criterion is set above zero (10^-4 in our light settings), all radial functions are checked individually, and their tails are cut off at a point where they are already close to zero.

You should note that the basis_dep_cutoff criterion usually does not matter at all, but for very large systematic basis set convergence studies (going to tier 3, tier 4, etc., and/or testing the cutoff potential explicitly), this value should be set to zero—as is done in the really_tight settings, for example.

A note on calculations beyond DFT (e.g., RPA, MP2, $GW$, BSE) or calculations of particularly sensitive observables (NMR). Here, specialized basis sets are available and tested with FHI-aims, with more information available in the literature:

• •

  Any standard Gaussian-type basis sets: See https://www.basissetexchange.org/. These basis sets work fine with FHI-aims. Our own variants of the associated species defaults are usually provided with numerically overconverged grid settings (they were intended for benchmarking) but efficient species defaults for Gaussian basis sets can also be created using the same general concepts described in this chapter.

• •

  RPA and MP2 with NAOs for light elements: NAO-VCC-nZ by Igor Ying-Zhang et al. 2013 [327]

• •

  Bethe-Salpeter Equation based on $GW$, for molecules (valence/conduction states): Chi Liu et al. 2020 [201]

• •

  Core level excitations by $GW$ or BSE based on $GW$: Yi Yao et al. [319]

• •

  Species defaults that appear to work with acceptable precision for periodic $GW$ are included in defaults_2020 as described in Sec. 2.2. See, however, the 2021 paper by Xinguo Ren et al. [260], which describes further modifications to our standard basis sets by adding tight Slater-type orbital basis functions (if needed).

The Hartree potential in FHI-aims is determined by a multipole decomposition of the electron density. The critical parameter here is the order (highest angular momentum) used in the expansion (all higher components are neglected). This value is chosen by:

    l_hartree           4

With this choice, energy differences are usually sub-meV converged also for large systems, but total energy differences, vibrational frequencies at the cm^-1 level, etc., may require more. Our tight settings,

    l_hartree           6

provide sub-meV/atom converged total energies in all our tests, but you may simply wish to test for yourself …

FHI-aims integrates its Hamiltonian matrix elements numerically, on a grid. However, this is an all-electron code: Performing integrations on an even-spaced grid (as is done in many pseudopotential codes) would provide terrible integration accuracy near the nucleus (sharply peaked, deep Coulomb potential and strongly oscillating basis functions).

Instead, we use what is a standard choice also in other codes: Each atom is given a series of radial spheres (shells) around it, and we distribute a certain number of actual grid points on each shell. Obviously, increasing the number of grid points (“angular” points) on each shell will improve the integration accuracy, but at the price of a linear increase in computational cost.

The fact that the integration spheres will overlap does not matter—we remedy this automatically by choosing appropriate integration weights (partitioning of unity, see CPC paper).

The number and basic location of radial shells is chosen by

    radial_base         36 5.0
    radial_multiplier   1

which means that we here choose 36 grid shells, and the outermost shell is located at 5 Å (this happens to be the outermost radius of each basis function, as dictated by the confinement potential).

The radial_base tag allows to increase the radial grid density systematically by adding shells inbetween those specified in the radial_base line. For example, we choose

    radial_multiplier   2

in our tight species defaults (for all practical purposes, this is converged), which means that we add one shell between the zero and the (former) first shell, one between the first and second, etc., and finally one between the (former) outermost shell and infinity … two times 36 plus one shells in total.

For an illustration of the effect of the radial_multiplier keyword on the density of the radial grid shells, go to Ref. [327] (http://iopscience.iop.org/1367-2630/15/12/123033/article, open access) and look at Figure A.1 and the accompanying explanation.

The distribution of actual grid points on these shells is done using so-called Lebedev grids, which are designed to integrate all angular momenta up to a certain order $l$ exactly. They come with fixed numbers of grid points (50, 110, 194, etc.). As a rule, fewer grid points will be needed in the inner grid shells, and more will be needed at the (more extended) faraway grid shells. We specify the increase the number of grid points per radial shell in steps, by writing:

     angular_grids specified
      division   0.2659   50
      division   0.4451  110
      division   0.6052  194
#      division   0.7543  302
#      division   0.8014  434
#      division   0.8507  590
#      division   0.8762  770
#      division   0.9023  974
#      division   1.2339 1202
#      outer_grid 974
      outer_grid 194

This example pertains to the light_194 settings (and very light elements), and means that only 50 points will be used on all grid shells inside a radius of 0.2659 Å, 110 grid points are used on all shells within 0.4451 Å, 194 grid points will be used on all shells inside 0.6052 Å—and that’s it! No more division tags are uncommented, and all shells outside 0.6052 Å also get 194 grid shells, as given by the uncommented outer_grid tag.

We note that the form of the increase of the number of points per radial shell near the nucleus, as well as the maximum number of angular grid points used outside a given radius are critical for the numerical accuracy in FHI-aims. When suspecting numerical noise anywhere in the calculations, the specification of the angular grid points should be checked first. This can be done by uncommenting further division tags with larger numbers of grid points, as well as a suitably increased outer_grid value. In particular, the choice of only 194 grid points max. per radial shell (only for the lightest elements!) is a rather aggressive choice, but in our experience still enables very reasonable geometry relaxations, structure searches or molecular dynamics for most purposes. However, the first thing to check in order to provide better convergence would be to set outer_grid to 302 (regular light settings). If this produces a noticeable change of the quantity you are calculating, be careful.

Of course, one can always introduce the denser grids provided in the tight settings, which (for reference) are

     angular_grids specified
      division   0.1817   50
      division   0.3417  110
      division   0.4949  194
      division   0.6251  302
      division   0.8014  434
#      division   0.8507  590
#      division   0.8762  770
#      division   0.9023  974
#      division   1.2339 1202
#      outer_grid 974
      outer_grid  434

These grids alone are roughly twice as expensive as the light_194 ones above, and should provide reasonable accuracy for pretty much any purpose. Nonetheless, of course one can still go and check explicitly, simply by increasing the number of grid points per shell by hand.