3.4 Specifying the basis (functions, empty sites, k-points, …)

This allows to place extra basis functions at specified locations outside the actual atoms, e.g., allowing for a counterpoise correction of basis set superposition errors. Although empty sites can be used with periodic boundary conditions, a special care should be taken when mixing periodic and non-periodic reference systems for counterpoise correction.

The definition of free atom radial functions used as the “minimal basis” of FHI-aims affects the absolute convergence of total energies calculated by FHI-aims. The radial shape of the free-atom core and valence functions towards the nucleus is nearly exact also for bonded structures if the same exchange-correlation functional is used and thus, e.g., using DFT-PBE generated radial functions when using xc pbe will improve the convergence of total energies. (Other quantities, such as atomization energies or other energy differences, will often exhibit better convergence than the total energy, since the effect of the exact shape of the minimal basis functions near the nucleus often cancels to a reasonable extent. See, e.g., Ref. [157] for a study of these effects for Gaussian-type and NAO basis sets compared to accurate reference values.)

Options for string:

sratom: FHI-aims’ default solver for the electronic structure of spherical free atoms on a dense logarithmic grid, called “sratom” (for “scalar relativistic atom”) is the same solver as used in the Fritz Haber Institute 1998 pseudopotential generation code by Martin Fuchs and coworkers, reference [101]. It was modified by Timo Jacob to incorporate ZORA scalar relativity when needed. This solver can handle semilocal density functionals but not exact exchange. Thus, for hybrid density functionals in FHI-aims, the “minimal basis” radial functions produced by sratom are semilocal DFT, not hybrid functional basis functions. This leads to a slower convergence of absolute total energies with hybrid functionals in FHI-aims (the error in atomization energies cancels to a good extent). See Figure 7 in Ref. [257] for the magnitude of this effect for the example of the convergence of Hartree-Fock calculations for Au₂ with an LDA-generated minimal basis compared to the same calculation, but with a minimal basis derived using the Krieger-Li-Iafrate (KLI) approximation [177, 178, 179] to the exact-exchange optimized effective potential. “sratom” is the current default solver for radial functions in FHI-aims.

atom_sphere: This solver for spherical free atoms was developed in Stefan Goedecker’s group for many years (beginning with Ref. [106]) and includes support for semilocal and hybrid density functionals. In FHI-aims, “atom_sphere” is used as a library under the Lesser General Public License (LGPL). As of this writing (August 2017), only non-relativistic calculations are supported by “atom_sphere”. In this case, the resulting total energies for hybrid density functionals converge precisely as well as their semilocal DFT equivalents, as shown, e.g., in Ref. [157]. “atom_sphere” will only work if support for libxc was compiled into the FHI-aims binary used.

One specific option for this keyword is:

KLI

for the Krieger-Li-Iafrate (KLI) approximation [177, 178, 179] to the exact-exchange optimized effective potential. Figure 7 in Ref. [257] is an example of the effect of the choice of different minimal basis functions on the convergence of the total energy of a Hartree-Fock calculation for a heavy element. In that case, the KLI approximation is closer to the converged result since the radial behavior of the minimal basis functions towards the nucleus is closer to the Hartree-Fock result, and because the radial function behavior of the free atom near the nucleus is nearly identical to the near-nuclear behavior of the Kohn-Sham orbitals in the bonded structures. In other words, for Hartree-Fock, KLI represents the nuclear cusp better than the local-density approximation (which is the present default choice of minimal basis radial functions for Hartree-Fock).

In addition, the following functionals have been implemented for the atomic solver and may be specified even if a different functional is chosen for the xc keyword:

pz-lda, pw-lda, vwn, vwn-gauss, pw91_gga, pbe, rpbe, revpbe, pbeint, pbesol, blyp

In all cases where atomic_solver_xc is not specified, the choice of the exchange-correlation functional for free atoms varies with the chosen setting for the xc keyword, but is not necessarily identical. For local and semilocal density functionals, the choice of functional for the minimal basis is generally identical to keyword xc, but for hybrid functionals, the choice may vary and is currently “best” documented in the source code (in subroutine get_free_atoms.f90).

An analagous keyword for the free-ion-like radial functions in the ionic basis part of the NAO basis sets is given by ionic_solver_xc.

This keyword automates a specialized version of the counterpoise correction and should only be used with great care. A general counterpoise correction for molecules can be implemented manually and in separate steps using the empty keyword. The atomization BSSE correction implemented above, if used, must be checked very carefully to ensure that all single-atom reference calculations carried out in the process reached the exact same atomic reference state. Many important atoms have more than one self-consistent solution, and mixing different self-consistent solutions may produce completely erratic results. Our recommendation therefore to carry out any counterpoise correction manually instead, by doing separate single-point FHI-aims calculations with different molecular fragments and different basis definitions.

The atomization BSSE correction for a molecular structure is defined as:

$${\mathrm{\Delta }}_{ac}=\sum _x[E^x(x)-E^x(sys)]$$

(3.14)

where $E^x(x)$ is the energy of the atom x calculated using only its basis set and $E^x(sys)$ is its energy calculated with the basis set of the whole structure. The BSSE corrected total energy is then

$$E^{BSSE}=E^{sys}(sys)+{\mathrm{\Delta }}_{ac}$$

(3.15)

The atomization BSSE correction is usually small in the case of the LDA/GGA functionals, but can become significant for methods with explicit correlation like RPA and MP2.

The BSSE corrected atomization energy implemented here, on the other hand, is

$$E_{at}^{BSSE}=E^{sys}(sys)-\sum _x{E}^x(sys)$$

(3.16)

In the case where relative energies between different conformations of the same molecule are needed, ${(E_{at}^{BSSE})}_{rel}\equiv {(E^{BSSE})}_{rel}$.

If the calculation of the full system is performed without spin polarization, the total energy of each atom will also be calculated without spin polarization (and vice versa). In this case, specially when performing the scf cycle with HF, symmetry breaking of the electronic configuration of certain atoms may occur, which might lead to wrong conclusions.

The keyword is currently implemented for use with cluster geometries only. Both RPA and MP2 as the total_energy_method are supported.

This tag should be kept at the default value unless for testing purposes (e.g., comparing to other codes).

Note that the order of n1, n2, n3 must directly correspond to the order in which the lattice_vectors are listed in geometry.in, through the definition and order of the reciprocal lattice. By default, the resulting k_grid is centered around the $\mathrm{\Gamma }$-point, but can be shifted using the k_offset keyword below.

The keyword symmetry_reduced_k_grid now allows to make use of time-reversal symmetry to reduce the number of $k$ points by a factor of nearly two. This option is the default.

This keyword generates a k_grid based on:

$n_i=\lceil \text{k\_grid\_density}\cdot {\displaystyle \frac{2\pi }{|{𝐚}_i|}}\rceil$

(3.17)

, where $n_i$ is the number of k-points along the $i$-th lattice vector, ${𝐚}_i$ is the $i$-th lattice vector, and $\lceil \cdot \rceil$ is the ceiling function.

Consider the following example. We specify a k_grid_density of 5.0 for the cubic Si diamond unit cell with a lattice constant of $a=5.43$ Å :

	${\displaystyle \frac{2\pi }{|{𝐚}_i|}}$	$={\displaystyle \frac{2\pi }{a}}={\displaystyle \frac{2\pi }{5.43\text{ Å}}}=1.157{\text{ Å}}^{-1}$		(3.18)
	$n_i$	$=\lceil \text{k\_grid\_density}\cdot {\displaystyle \frac{2\pi }{|{𝐚}_i|}}\rceil$		(3.19)
		$=\lceil 5.0\cdot 1.157\rceil$		(3.20)
		$=\lceil 5.785\rceil$		(3.21)
		$=6.$		(3.22)

This results in a k_grid of 6x6x6.

Often, the sampling of the Brillouin zone based on the input geometry is not optimal (especially when you have very skewed systems). The autoGR[224] tries to find a symmetry conserving supercell of the input geometry (=Born-von Karman cell) that results in a better sampling of the Brillouin zone.

Can be used to shift the grid off-$\mathrm{\Gamma }$ for better $k$-space sampling. For example, (0.5, 0.5, 0.5) together with even n$i$ for k_grid defines a Monkhorst-Pack [223] type grid. See Sec. 4.3 for details.

The k_offset cannot be used with Hartree-Fock or hybrid density functionals as the current implementation assumes a $\mathrm{\Gamma }$-centered grid.

This option is useful to specify an uneven special $k$-point set [59], etc.

The k_points_external keyword cannot presently be used with Hartree-Fock or hybrid density functionals as the current implementation the presence of the k_grid keyword and information in control.in.

If specified, k_points_external expects a separate input file k_list.in to be provided in the same working directory as all other input files. The format of k_list.in is as follows:
line 1: n1 n2 n3 line 2: Nk lines 3 - Nk+2: k1 k2 k3 weight

All lines are read as free format.
n1, n2, n3 are integers, specifying the descriptors of a 3D $k$-point grid. During s.c.f., these values are for reference only. However, they may be used for postprocessing after the s.c.f. cycle is complete, e.g., for keyword output postscf_eigenvalues, where the external k-point list k_list.in is not supported.
Nk is the (integer) number of $k$-points following in the file.
For each $k$-point, a separate line is expected, including via k1, k2, k3 the coordinates of this point in units of the reciprocal lattice vectors, and via weight the integration weight of this $k$-point in all Brillouin zone integrals.

See keyword output onsite_integrands for the integrals that are actually evaluated.

For “light” settings, you can most likely ignore the associated warning. Possibly, treat it as a reminder to check final results with “tight” settings or similar, just in case.

FHI-aims writes the respective onsite integral values for all its radial functions side by side after the setup of all radial functions is complete. Too large deviations between the calculated values on the radial and logarithmic grids can indicate accuracy problems, which is a particular concern for high-accuracy benchmark calculations. In particular, Gaussian-type orbital basis sets of the Dunning type that are used for high-level reference calculations need grids that are far more accurate than normal NAO basis sets, due to the unphysical wiggles that contracted Gaussian functions with high exponents introduce near the nucleus.

If the associated warning strikes, the grid accuracy can be improved either using the radial keyword or (simpler) the radial_multiplier keyword.

However, if onsite_accuracy_threshold triggers a warning for “light” settings, most likely you can safely ignore the warning. With light settings, the point may be to do a reduced-accuracy calculation, which should still be safe for its original purposes. Just check the radial grid when in doubt.

If the flag override_integration_accuracy is toggled to .false., however, FHI-aims does stop whenever it encounters a radial function whose onsite integrals are not deemed accurate enough by the onsite_accuracy_threshold criterion.

A radial function is considered zero (not evaluated) beyond the radius where $u(r)$ and its first and second derivatives become smaller than threshold. The default is chosen such as to not affect any results at all.

For electron densities or orbitals plotted for visualization using cube files (output cube functionality), a too high value of wave_threshold can sometimes lead to small but visible discontinuities. Thus, the default threshold is lowered to 10^-8 if outputcube is requested.

Before any calculation, all radial functions for a single species are Gram-Schmidt orthonormalized. If the norm of the function after orthonormalization is smaller than threshold, that function is omitted.

If not .false., the onset of the basis confining potential (see cut_pot tag below) is adjusted separately for each basis function, such that the norm of this basis function outside $r_{\text{onset}}$ is smaller that threshold. The maximum possible onset radius is still given by the value explicitly specified by the cut_pot tag.

The defining potential for this basis function type consists of the non-spinpolarized, self-consistent spherical free-atom potential (possibly itself confined, using the cut_free_atom tag), and a confining potential. The shape of the confining potential is the same for all basis functions of a given species, and set using the cutoff_type and cut_pot subtags.

Currently not needed for production calculations, but listed here because the “core” infrastructure is currently being reworked and may see useful additions in the near future. This flag defines which electrons of the species are considered “core” electrons, and which enter as explicit valence electrons.

This keyword can be used to introduce a non-spinpolarized, radially symmetric core hole into a species definition. The minimal basis functions for this species will then be generated in the presence of the core hole.

Example:
core_shell_occ 2 p 5.0

specifies the occupation of the 2p core shell with 5 electrons instead of the usual 6. Whenever core_shell_occ is used, the occupation of the valence levels must be adjusted accordingly using the keyword valence, to keep the free atom neutral.

Note that this keyword only affects the behaviour of the atomic solver. It does not modify the occupation of the core levels in the actual FHI-aims calculation - for that, one can use deltascf_basis or deltascf_projector.

Experimental at present, not needed for any production purposes. See also the core keyword. The core_states keyword should interact with the core keyword, but does not yet, since any associated functionality is still under development. Both keywords are listed here for future reference only.

This keyword applies only to the setting specified by keyword basis_dep_cutoff. Do not enable it routinely without thorough testing.

By default, the minimal basis functions in FHI-aims are subject to the cutoff potential with the fixed onset specified by the cut_pot keyword. However, the more restrictive basis_dep_cutoff keyword does not apply to the minimal basis by default.

This can be changed by setting cut_atomic_basis to .true., but the associated total energy changes are significantly larger than for other basis functions. By default, we therefore do not recommend adding a tighter cutoff to the minimal basis functions at this time. It is, however, possible that this effect is mainly a systematic error between the core states, and after further testing, we may yet choose to enable this keyword if we can guaranteed that its use is safe.

This tag is mandatory, since it specifies onset, a critical parameter that allows to tune the efficiency of a calculation for a given target accuracy. Unless reduced by the basis_dep_cutoff tag, onset is the default onset radius used to construct all valence (minimal) and hydrogen-like basis functions of this species. In addition, any confined free-atom or free-ion like radial functions use this onset radius if auto is used in their specification.

Notes: The functional form of $v_c(r)$ can be selected using the cutoff_type keyword, and width and scale apply to this shape. Modifying these latter parameters is usually not necessary for a production calculation, but the onset value should be verified at least as a quick numerical check.

All confinement potentials in FHI-aims are characterized by the rigorous boundaries $v_c(r)$=0 for and $v_c(r)$=$\mathrm{\infty }$ for $r>r_{\text{cut}}=r_{\text{onset}}+w$, where $r_{\text{onset}}$ may depend on the basis function, and $w$ is the width specified by the cut_pot tag. In addition, each shape contains a scaling parameter $s$, also specified via the cut_pot tag.

Available confinement potential shapes (identifier) for are:

• •

  exp(1_x)_(1-x)2 :

  $$v_c(r)=\exp (-\frac{w}{r-r_{\text{onset}}})\cdot \frac{1}{{(r-r_{\text{cut}})}^2}$$

  (the default in FHI-aims)

• •

  junquera :

  $$v_c(r)=\exp (-\frac{w}{r-r_{\text{onset}}})\cdot \frac{1}{(r-r_{\text{cut}})}$$

  (the form originally suggested by Junquera et al. [161])

• •

  x2_(1-x2)

  $$v_c(r)={(r-r_{\text{onset}})}^2\cdot \frac{1}{{(r-r_{\text{cut}})}^2}$$

Note that the default shape is Equation (9) in reference [36], except that the “$-$” sign in the exponential is missing in the reference. The source code, of course, has the “$-$” sign as it should be, ensuring that the confinement potential and all its derivatives go to zero at $r=r_{\text{onset}}$. Thanks go to Susi Lehtola for spotting this.

FHI-aims allows to use Gaussian-based radial functions to compare to existing popular Gaussian-based implementations of quantum chemistry. These functions can either be primitive Gaussians,

$$u(r)=B(L,\alpha )r^{L+1}\exp (-\alpha r^2),$$

(3.23)

where $B(L,\alpha )$ is the normalization constant for a primitive gaussian, or non-primitive linear combinations.

$$u(r)=\frac{1}{\mathrm{𝑁𝑜𝑟𝑚}}{r}^{L+1}\cdot \sum _{i=1}^{i=N}{g}_i{B}_i(L,{\alpha }_i)\exp (-{\alpha }_i{r}^2).$$

(3.24)

In existing quantum chemistry codes, Gaussian basis functions can be defined either as spherical Gaussians [Eq. (3.23) above], or as cartesian Gaussians,

$$\varphi (𝐫)=x^k{y}^m{z}^n\exp (-\alpha r^2),\mathrm{where}k+m+n=L.$$

(3.25)

This behavior can be mimicked using the pure_gauss tag (see below). Finally, note that in order to use an exclusively Gaussian-based basis set, you must prevent the use of the minimal free-atom like NAO basis functions using the include_min_basis tag.

This flag is normally only useful to compare explicitly with basis sets from other methods, usually Gaussian basis sets.

With Gaussian basis sets and any other basis sets that should not include any radial functions of the numeric atom-centered orbital spherical free atom, include_min_basis must be set to .false.

If include_min_basis is set to .true., the basis set will additionally include the numerical spherical free-atom radial functions. For our standard FHI-aims numeric atom-centered orbital basis sets, this behavior is desired and part of their definition. However, if these basis functions are added to a standard Gaussian-type basis set, then total energies and other computed quantities will not be the same as, say, with a pure Gaussian-type orbital code.

Note that, due to the presence of the free-atom radial functions, our numeric atom-centered basis sets usually give lower total energies for DFT than standard Gaussian basis sets, because the description of the region near the nucleus is more precise. This is demonstrated quantitatively, including figures, in the appendix of Ref. [327].

This tag defines an ionic configuration, but does not actually add any ionic functions to the basis used in the present calculation. Actual basis functions are added by the ionic keyword.

Only the topmost valence shell of each angular momentum channel of the ion is specified. All lower-lying shells are assumed to be completely filled for the self-consistent spherical free-ion calculation.

See keyword gaussian for the distinction between spherical and cartesian Gaussian functions. In short, cartesian and spherical Gaussian functions are equivalent except that for a given $L$, cartesian Gaussians add a so-called angular momentum contamination. If pure_gauss is specified, this angular momentum contamination is mimicked by FHI-aims.

Consider the simple (textbook!) three-dimensional harmonic oscillator in quantum mechanics. This can be solved either in cartesian coordinates, or in spherical coordinates. If solves in cartesian coordinates, you will find that there are six degenerate solutions for the principal quantum number 2, five of which correspond to $l$=2 ($d$ channel), but one of which corresponds to $l$=0 ($s$ channel). This is the exact angular momentum contamination exhibited by the cartesian definition of a Gaussian basis function.

Only the topmost valence shell of each angular momentum channel of the atom is specified. All lower-lying shells are assumed to be completely filled for the self-consistent spherical free-atom calculation. The valence occupation must be defined explicitly for each species.

The self-consistent free-atom potential generated by this calculation is used to generate all minimal and confined basis functions used for this species, after the confining potential is added.

The self-consistent free-atom calculation can itself be confined by a different confining potential, the onset of which is specified by the cut_free_atom keyword.

For DFT-LDA/GGA, the same xc functional that is used in the full three-dimensional calculation is also used to define the self-consistent free atom. For any methods involving Hartree-Fock exchange (e.g., hybrid functionals), the free atom is generated using the pw-lda LDA functional.

The self-consistent free atom density generated here is also used in the construction of partition functions for the Hamiltonian integrals and the Hartree potential, as well as to build the initial charge density (unless otherwise requested!) and the reference charge density subtracted before constructing the Hartree potential.

This keyword is only currently implemented forspecies_default_type minimal$+$s. This line in the species defaults enables the force, energy, and stress corrections for use with minimal$+$s type species defaults – see species default section for more information on these basis sets. To change the $c$, and $r_{cut}$ parameters of these corrections, see the entry of the small_basis_param keyword.