3.5 Integration, grids, and partitioning

See grid_partitioning_method and Ref. [128] for details regarding integration batches.

This option need not be changed (or invoked) in any normal runs, since there is no quantitative difference between integrals with Delley’s and Lebedev’s tabulated grids to our knowledge.

Lebedev’s grids may be explicitly invoked when denser angular grids than 1202 points (already very dense!) per radial integration shell around each species are required. In detail, grids with the following numbers of grid points are provided:

• •

  Delley : 6, 14, 26, 50, 110, 194, 302, 434, 590, 770, 974, 1202

• •

  Lebedev : 6, 14, 26, 38, 50, 86, 110, 146, 170, 194, 302, 350, 434, 590, 770, 974, 1202, 1454, 1730, 2030, 2354, 2702, 3074, 3470, 3890, 4334, 4802, 5294, 5810

These numbers of grid points can be invoked in the subtags of the angular_grids specified description for fixed angular grids (the default in the preconstructed species_defaults files), and in further tags such as angular or angular_min.

Partitioning the integration grid properly can be performance-critical for the expensive grid-based Hamiltonian integration and charge density update steps. Details on these methods are given in Ref. [128]. In particular, we support:

• •

  maxmin : The default for serial computations: the “grid-adapted cut-plane” method described in Ref. [128]

• •

  parallel_hash+maxmin : The default for all parallel runs. We first hash the grid points to tasks by the geometric location and then run a maxmin algorithm locally in each task.

• •

  parallel_maxmin : Memory-parallel implementation of the maxmin method. However, the exact implementations requires rather much communication, and has been superseded by parallel_hash+maxmin.

• •

  octree : The “octree” method described in Ref. [128]

• •

  parallel_octree Parallel version of the “octree” method described in Ref. [128]. Only useful for research purposes—superseded by parallel_hash+maxmin for all practical applications.

• •

  octant Simple partitioning of the grid into “octants” of each radial integration shell arround each atom.

Note that additional parameters may be invoked to specify details for these methods, most importantly batch_size_limit and points_in_batch.

We mention for completeness that FHI-aims supports further, experimental grid batching methods, including the possibility to link to external libraries. The associated method strings are tetgen+metis, qhull+metis, nearest+metis, and group. As discussed in Ref. [128], the conceptually simpler maxmin method performs as well or even better than these “bottom-up” type approaches, and should be preferred.

As documented in issue 614 of the FHI-aims Gitlab, https://aims-git.rz-berlin.mpg.de/aims/FHIaims/-/issues/614, it was noted that the Graphite exfoliation PES had a discontinuous energy jump at certain interlayer separations, with a lightdense or lightdenser basis set and xdm. It was determined that this also affected many_body_dispersion, many_body_dispersion_nl, vdw_correction_hirshfeld, and vdw_ts, and was due to an incomplete list of centers used in the calculation of the Hirshfeld partition weights. This incomplete list occurs when the outer extent of the radial grid (see the radial_base keyword) is larger than the outer extent of the cutoff potential (see the cut_pot keyword) for periodic systems with sparse unit cells.

This keyword enables the bugfix to this problem and was introduced in FHI-aims version 250702. It permits use of an enlarged centers list to calculate the Hirshfeld partition weights. Setting hirshfeld_partition_centers_2025 to .false. recovers the behaviour before this change, in case it is needed for backwards compatibility or checking of old behaviour.

As this is a bugfix, providing there are no problems with this change of centers list, it should become default behaviour in FHI-aims. Therefore, this keyword will be deprecated in approximately 1 year, and it will be no longer possible to use the old centers list.

Please also see the keyword partition_centers_2025, and issue 737 of the FHI-aims Gitlab, that this addresses, https://aims-git.rz-berlin.mpg.de/aims/FHIaims/-/issues/737, as these problems have essentially the same origin.

No need to tweak for standard production calculations. See grid_partitioning_method for details regarding integration batches.

See Ref. [36] for details regarding “partitioning of unity” of the charge density in integrations and the Hartree potential. The partition functions $p_{\text{at}}(𝐫)$ are only calculated if their denominator (the norm; e.g. ${\sum }_{\text{at}}{\rho }_{\text{at}}/{|𝐫-{𝐑}_{\text{at}}|}^2$ is greater than value, else that integration point is ignored.

Notice that this type of partitioning is strictly rigorous for integrands that extend no further than the free-atom like densities used to define our partition functions. This is always true for DFT-LDA/GGA, with NAO’s but if you suspect (e.g., with very diffuse Gaussian basis functions) some kind of integration noise, reducing threshold may be a good first test.

As documented in issue 737 of the FHI-aims Gitlab, https://aims-git.rz-berlin.mpg.de/aims/FHIaims/-/issues/737, it was noted that the atomic centers list used by the different flavours of the Stratmann partition function was insufficent when the outer extent of the radial grid (see the radial_base keyword) is larger than the outer extent of the cutoff potential (see the cut_pot keyword). This is because this list of centers is populated based on distances between centers being less than some criteria that depends on the outer extent of the cutoff potential. This meant that some centers in cells adjacent to the 0-cell may not be picked up, despite having grid points nearby. In many cases, this does not happen, however, it is particularly problematic for sparse unit cells and can lead to energy jumps at certain unit cell sizes, as was demonstrated in issue 737 for He in different sized unit cells.

The solution to this problem is to expand the list of centers used by the different flavours of the Stratmann partition function, by using the outer extent of the radial grid as a distance criteria, rather than the outer edge of the cutoff potential, in cases where the radial grid is larger. This keyword partition_centers_2025 automatically applies this fix, and was introduced in FHI-aims version 250626, via merge request 1857 (https://aims-git.rz-berlin.mpg.de/aims/FHIaims/-/merge_requests/1857). Setting partition_centers_2025 to .false. recovers the behaviour before this change, in case it is needed for backwards compatibility or checking of old behaviour.

As this is a bugfix, providing there are no problems with this change of centers list, it should become default behaviour in FHI-aims. Therefore, this keyword will be deprecated in approximately 1 year, and it will be no longer possible to use the old centers list.

Usually, this tag need not be modified from the default. See the Computer Physics Communications description of FHI-aims for a description of the numerical integration technique used in FHI-aims.

In brief, each extended three-dimensional integrand is broken down into atom-centered pieces, using a set of localized, atom-centered partition functions:

$$p_{\text{at}}(𝐫)=\frac{g_{\text{at}}(𝐫)}{{\sum }_{{\text{at}}^{\prime }}{g}_{{\text{at}}^{\prime }}(𝐫)}.$$

(3.26)

where $g_{\text{at}}(𝐫)$ is an atom-centered weight function. The following options for type are available:

• •

  rho_r2 : $g_{\text{at}}(𝐫)=n_{\text{at}}^{\text{free}}(r)/r^2$
  (first suggested by Delley [73]).

• •

  rho_r : $g_{\text{at}}(𝐫)=n_{\text{at}}^{\text{free}}(r)/r$

• •

  rho : $g_{\text{at}}(𝐫)=n_{\text{at}}^{\text{free}}(r)$

• •

  fermi : Deprecated—do not use. A Fermi-function like approach, requires two additional parameters.

• •

  stratmann: The shape suggested by Stratmann et al., Ref. [287]. This saves $\approx$10-20 % of the numerical effort compared to rho_r2. More importantly, however, our recent testing shows that stratmann is also significantly accurate in some corner cases where the effects of integration accuracy even make a difference. Note the properly bounded “stratmann_smoother” default function below. Straight “stratmann” should not be used.

• •

  stratmann_smooth: Partial update to guarantee a smooth edge at the “outer radius” or atoms.

• •

  stratmann_smoother: Corrected version of the stratmann partition table. The following explanation refers to the prescription given in Eqs. (8), (11), and (14) of Ref. [287]. The actual (normalized) partition function is given by Eq. (10). At each grid point, it depends on a product of cell functions Eq. (9) over potentially all atoms in the system—unless its cell function is equal to one, a faraway atom may contribute to the partition function at a given grid point. Through the definition of ${\mu }_{ik}$ in Eq. (4) of Ref. [287] and the limitation of the cell function to unity for through Eqs. (11) / (14), the distance from which an atom can contribute is restricted, but potentially to a very large radius indeed. This becomes a problem for periodic systems (in our case, a theoretical radius of 25 Å would have resulted even for light settings and the farthest grid points from each atom, set to 5 Å for light settings).
  To avoid an overly large volume of contributing atoms, we restrict the list of contributing atoms to only those whose free-atom charge density would not be zero at the integration point in question. To that end, Eq. (8) of Ref. [287] is additionally multiplied with a function that smoothly interpolates between the original $s$ of Stratmann and coworkers and unity. The interpolation is done only for atom distances between 0.8 and 1.0 times their free-atom radius, and uses a $[1-\cos (2x)]$-like interpolating function. The bottom line is that we get the benefits of both the Stratmann table and a restricted atom list without any discontinuities or wiggles as a function of atomic positions or unit cell vectors — which is as it should be.

• •

  stratmann_sparse: This version of the Stratmann partition table is the same as stratmann_smoother, but it stores the relevant interatomic distances in a memory saving form.

Note that the free atom electron density $n_{\text{at}}^{\text{free}}(r)$ still determines the extent of many partition function types. This is controlled by the cut_free_atom keyword. See also the hartree_partition_type keyword, which presently must have the same setting as the partition_type keyword.

Please also see the partition_centers_2025 for a change in how the stratmann, stratmann_smooth, stratmann_smoother, and stratmann_sparse partition functions are computed, as of version 250626.

See grid_partitioning_method and Ref. [128] for details regarding integration batches.

By default, the code stops when the charge integration error at the end of a SCF iteration is larger than ${10}^{-5}$. This is not always an indication for an erroneous calculation, but has indicated faulty calculations in the past, where the MPI communication was wrongly set up on HPC systems. When set, the code assumes that the user must know what they are doing.

This option is only meaningful for self-adapting angular grids, which are not the recommended default for production calculations with FHI-aims – (i) because these grids are often rather dense, and (ii) because they are meaningful only for cluster-type geometries. In order to specify self-adapting angular grids anyway, you must also set the keywords angular_min and angular_acc.

The available values of integration points in given angular grids are listed with the keyword force_lebedev.

If threshold is not zero, this option invokes the self-adaptation of all angular integration grids, within the limits given by angular_min and angular. The adaption criteria are the initial Hamiltonian / overlap matrix integrals.

In all preconstructed species_default files, we specify reliable angular integration grids for all elements for DFT. No adaptation is required. For the curious, our own grids are adapted for symmetric dimers at a tight bond distance, using threshold = ${10}^{-8}$.

The standard species_default files provided with FHI-aims provide specified angular grids (on the safe side, i.e., rather dense) for each radial_base integration shell around an atom. The line:
angular_grids specified
must be immediately followed by a series of lines with
division […]
outer_grid […]
tag(s). These contain the actual grid specification.

If method auto is given, appropriate specifications for self-adapting grids should be included in control.in (keywords angular, angular_min, angular_acc).

For specified angular_grids, acts as a lower bound for the number of points per radial shell (specified grids will be increased accordingly).

For self-adapting angular grids, use together with the angular and angular_acc keywords.

In practice, value will be reduced to the next-highest available Lebedev integration grid (see force_lebedev tag for possible values).

Although this is a technical parameter (ideally, no influence on self-consistent, converged results), it has important implications for a variety of numerical tasks in the code:

• •

  It influences (slightly) the basis-defining potential for the minimal basis, and for confined basis functions.

• •

  It limits the radius of the free-atom density, which in turn limits the extent of the default integration partition table. For DFT-LDA/GGA, this extent need must not be smaller than the radius of the most extended basis function, but it also need not be larger, since all integrands are zero outside anyway. This is not the case for the two-electron Coulomb operator, which is needed for Hartree-Fock, hybrid functionals, $GW$, etc, in which case the default is currently infinite (no cutoff potential applied).

• •

  It also limits the extent of the partition table used for the Hartree potential.Especially in periodic calculations, it is vital that the real-space part of the Hartree potential is kept small. In that case, it is thus critical to keep radius as small as possible.

Usually, the default specified in the code should be accurate for all requirements. If, however, you suspect some kind of integration noise which is not related to the grid, increasing the cut_free_atom value may be a good test.

Use the outer_grid tag to specify the number of angular grid points used outside the outermost division radius.

If, after on-site orthonormalization, a radial function’s innermost extremum is inside the radial grid shell number, counting from the nucleus, that radial function is rejected in order to prevent inaccurate integrations.

The number of logarithmic grid shells, $N$, is uniquely determined by r_min, r_max, and increment. Specifying a dense logarithmic grid is not performance-critical.

Use the division tag to specify the number of angular grid points used for radial shells within specified radii.

The location of the number radial shells is given by

$$r(i)=r_{\text{outer}}\cdot \frac{\log \{1-{[i/(N+1)]}^2\}}{\log \{1-{[N/(N+1)]}^2\}}$$

(3.27)

With this prescription, shell $i$=0 would be located exactly at $r(i)=0$, and shell $i$=$N+1$ would be located exactly at $r(i)=\mathrm{\infty }$, i.e., this provides an exact mapping of the interval [0,$\mathrm{\infty }$].

The FHI-aims species_default files provide values for number according to the formula $N$=1.2$\times$14(nucleus+2)^1/3, as determined empirically in Ref. [21]. These “basic” grids are can then be augmented by adding uniform subdivisions, using the radial_multiplier keyword described below.

The basic grid of $N$ radial shells (see radial_base definition) is evenly subdivided number-1 times, to yield number$\cdot (N+1)-1$ shells in the actually used integration grid. Thus, the radial_multiplier tag allows to systematically increase the number of radial shells (by factors). For example, number=2 (the default) would yield $2N+1$ shells total.

Note that some all-electron Gaussian basis sets contain either very high or very low exponents. If such basis sets are used for test purposes, it may be necessary to test the convergence of the radial integration grid directly by increasing the radial_multiplier.

The effect of theradial_multiplier is explained in Ref. [327] (open access at http://iopscience.iop.org/1367-2630/15/12/123033/article) Look at Figure A.1 and the accompanying explanation in that reference.