Main-group test set for materials science and engineering with user-friendly graphical tools for error analysis: systematic benchmark of the numerical and intrinsic errors in state-of-the-art electronic-structure approximations

The importance of reliable test sets for the success of quantum-chemistry methods was first realized by Pople and co-workers [1, 26, 41–44]. Along with the development of the Gaussian-n theories [1, 26, 41–44], a hierarchy of extrapolating levels of correlation and basis sets have been developed to obtain increasingly accurate thermochemistry. A sequence of pioneering test sets for quantum chemistry were also developed and used for methdological validation. These test sets are now widely recognized as the Gn test sets, addressing atomization energies and other energetic properties of molecules with increasing numbers of accurate reference values obtained from experiments [1, 26, 41–44]. One of these test sets, G3/99 [26] generated in 1999, has been widely used in the development of density functional approximations (DFA) to describe covalent bonding in the main group molecules [1, 25–38]. Since then, many well-established test sets have been proposed for different, mainly ground-state properties of molecules and molecular processes [1–15, 26, 42]. For example, Hobza and co-workers designed the S22 set, comprising 22 non-covalent binding complexes of biologic relevance, which can be further divided into 7 systems for hydrogen bonds, 8 for dispersion bonds, and 7 for mixed bonds [7]. Moreover as one of the most popular standard benchmark sets, the BH76 set proposed by the Truhlar group consists of forward and reverse barrier heights of 19 hydrogen transfer reactions, 6 heavy atom transfer reactions, 8 nucleophilic substitution reactions, and 5 unimolecular and association reactions [4]. Most recently, the Grimme group compiled a comprehensive benchmark test set, GMTKN30, which includes 30 subsets collected from the literature, covering a large section of chemically relevant properties of the main-group molecules [14, 15].

It may appear plausible that reference data, i.e. the accurate values of the relevant chemical and physical properties of atoms and molecules, can be acquired from experiments [1–3, 6, 11–13, 26, 42], but the required accuracy is not always achievable. The influence of electron-vibrational coupling is often unclear, which can hamper the comparison with the directly computed values. Thus, reference data would ideally be created by accurate quantum chemistry methods such as the coupled-cluster (CC) method with single, double, and perturbative triple excitations, CCSD(T) [45–47], or single, double, triple and perturbative quadruple excitations, CCSDT(Q), or full configurational interaction (full-CI) theory extrapolated to the complete basis set (CBS) limit [48–50]. For the test sets aforementioned, the corresponding reference data are either carefully chosen from accurate experiments [1–3, 6, 11–13, 26, 42] or obtained from high-level first-principles calculations [4, 5, 7–10].

As the most popular choice in quantum chemistry, the atom-centered Gaussian-type orbital (GTO) basis sets provide well-converged total energies of atoms and molecules with a reasonable basis set size for the mean-field methods, including the local-density approximation (LDA), generalized gradient approximations (GGAs) and hybrid functionals in DFT, and the Hartree–Fock method in WFT. For an all-electron approach and for heavier atoms (e.g. Z > 18), i.e. when wave-functions oscillate more strongly at the core region, it is more convenient, efficient, and accurate to use numerical atom-centered orbital (NAO) basis sets [51, 52].

For advanced correlation methods which require unoccupied single-particle states, e.g. the second-order Møller–Plesset perturbation theory (MP2) [53], the random-phase approximation (RPA) [54–57], or CCSD(T) [45–47, 58], the slow basis set convergence is a more serious problem [23, 24, 59–62]. It can be ascribed to the inaccurate description of the electron-electron Coulomb cusps using smooth orbital product expansions [63, 64]. The so-called correlation consistent basis sets have been proposed [24, 59], which allow for an analytic extrapolation to the CBS limit. Alternatively, by introducing an explicit dependence on the inter-electronic distance into the wave-function (F12 strategies [64–68]), it is possible to treat the cusp explicitely and suppress the basis set incompleteness error at a finite basis set size. Both techniques were proposed to address the basis set incompleteness errors in advanced correlation methods for both atoms and molecules, and they have been demonstrated to be very successful for light elements (e.g. Z < 18) and for ground-state properties [24, 59, 64, 68]. Therefore, accurate chemical or physical properties, e.g. reaction energies, reaction barrier heights and isomerization energies, can be obtained without concern for the numerical error cancellation originating from a given finite basis set. However, the quantum-chemistry test sets are mainly available only for the ground-state properties of light elements and small molecules, as it remains a challenge to obtain accurate reference data with advanced correlation methods for heavy elements (e.g. Z > 18), excited states, and large systems.

As mentioned earlier, in computational materials science, the situation is more complex and less developed than in the quantum chemistry of molecules, and there is an urgent need for accurate test sets to support the developement of advanced DFAs for solids. Staroverov and co-workers considered lattice constants, bulk moduli, and cohesive energies of 18 solids, as well as jellium surface energies to benchmark the TPSS functional, a meta generalized gradient approximation (meta-GGA) [69]. They also presented results for LDA, GGA PBE, and meta-GGA PKZB [70]. Heyd and Scuseria [71] used lattice constants and bulk moduli for 21 solids and band gaps for 8 semiconductors out of this set to investigate the screened Hartree–Fock-exchange properties in the HSE hybrid functional. Comparison was also made with LDA, PBE, and TPSS. This test set (or part of it) has been used to benchmark the HSE and HSEsol as implemented in VASP [72, 73], to understand the failure of B3LYP for solids [74], and to develop improved DFAs for solids [75–77]. A set of 60 cubic solids was used by Haas, Tran, and Blaha to compare the performance of different GGA functionals, LDA and TPSS in describing lattice constants and bulk moduli [78–80]. To compare the accuracy of different van der Waals (vdW) functionals applied to solids, Klimes and co-workers determined cohesive properties for a set of 23 bulk solids [81], based on a set of materials and properties used by Csonka and co-workers when testing the accuracy of GGA and meta-GGA functionals [82]. Recently, Tran et al performed an extensive test on the lattice constant, bulk modulus, and cohesive energy for a large set of XC approximations belonging to rungs 1–4 of 'Jacob's ladder of DFT' [83], performed in a non-self-consistent way using PBE orbitals. The test set used contain 44 strongly and 5 weakly bound solids, most of which were cubic, except for hexagonal graphite and h-BN. These test sets of materials with different properties have already been playing an important role in developing and benchmarking the DFAs for solids. These test sets exclusively rely on reference values from experiment; no quasi-exact calculated results exist, which can be used as reference data for solids.

For reasons of comparability and consistency, reference values from theory provide the unique opportunity to compare calculations based on exactly the same atomic structure and exclude finite-temperature, zero-point vibration, relativistic and electron-vibration coupling effects. Recently, the 'gold standard' method in quantum chemistry, CCSD(T), has gained attention in materials science [84–89]. A significant step towards an exact description of solids was demonstrated by Booth et al by performing a full-CI quality calculation with the aid of the Quantum Monte Carlo (QMC) stochastic strategy [85]. However, the numerical difficulty and computational complexity for describing solids has strongly limited the applications of these highly accurate quantum chemistry methods to only model systems with very few electrons, basis functions, and ${\boldsymbol{k}}$ points [84, 85, 88, 90]. Other QMC stochastic strategies, such as Variational Monte Carlo (VMC) [91] and Diffusion Monte Carto [92] feature a low computational scaling with respect to the system size and a weak dependence on basis set. QMC has been demonstrated to be as accurate as aforementioned quantum chemistry methods, e.g. CCSD(T), if certain numerical issues could be considered properly, which includes the fixed-node error and the form and optimizaiton of the trial function [93–95].

Unlike the molecular systems investigated by quantum chemistry, extended periodic materials are often simulated in reciprocal space with periodic boundary conditions (PBC) so as to fully leverage the periodic symmetry. However, it also introduces extra numerical difficulties: the first Brillouin zone in reciprocal space must be sampled by a finite number of ${\boldsymbol{k}}$ points, for which different kinds of ${\boldsymbol{k}}$-mesh have been proposed to provide an efficient sampling according to different periodic symmetries [96–100]. Basis set incompleteness error is another big issue: the atom-centered GTO-type basis sets that are dominant in quantum chemistry become cumbersome for condensed matter systems [101, 102]. In particular, large and especially diffuse GTO basis sets might cause severe ill-conditioning problems due to the lack of rigorous orthogonality for GTO basis functions. The plane-wave basis sets are a popular alternative choice for solids as they provide an intrinsically improvable description of the electronic orbitals with only one parameter: the plane-wave cutoff energy. However, the delocalized nature of plane waves makes them inefficient to describe the localized core electrons surrounding atomic nuclei. In practice, the pseudo-potential approximation, which removes the core electrons from explicit electronic-structure calculations, are often used together with the plane-wave basis sets in the projector augmented wave (PAW) framework [103, 104]. More recently, it has been demonstrated that the use of compact NAO basis sets [52, 105] or the linearized augmented plane-wave (LAPW) strategies [106] are able to converge the mean-field approximations in all-electron and full-potential resolution with relativistic effects explicitly included as opposed to implicitly as is the case for pseudo-potential methods.

Such diversity of methods leads to a need for a standard data set that is comparable among different numerical implementations. This is indeed the purpose of the Δ-value concept introduced by Cottenier et al [107, 108], which is currently focusing on the implementation of the PBE method for solids. The numerically well-converged PBE results from the all-electron, full-potential LAPW-based program WIEN2k [108, 109] are taken as reference. To quantify the numerical errors in the description of equations of state, the Δ-value is defined as the root-mean-squared (RMS) energy difference between the equations of state of two codes, averaged over an exhaustive test set of crystalline solids, containing all ground-state elemental crystals.

Despite the enormous success achieved by the (semi)-local DFA, e.g. LDA and GGA PBE, there are several notorious failures of these methods. For instance, there exist serious self-interaction errors and significant underestimation of vdW interactions for these functionals, which demand more sophisticated approximations. The meta-GGAs, in which the orbital kinetic energy density is added to the density functional evaluation, belong to a higher rung of the 'Jacob's Ladder in DFT'. Representative are the nonempirical TPSS [69] and the multiparameter empirical M06-L [110]. A recently proposed nonempirical meta-GGA, SCAN, which satisfies 17 known exact constraints appropriate to a semilocal functional [111], has been found to show a big step ahead in accuracy for manifold chemical and physical properties and has attracted increasing interest in computational materials science. The so-called hybrid functionals, e.g. HSE06 [112], incorporate the information of the occupied orbitals in Hartree–Fock-like 'exact exchange'. The next level of complexity is to derive DFAs using the information of virtual orbitals. Two methods at such level, MP2 [53] and RPA [18, 34, 35, 54–56, 113], are state-of-the-art in computational materials science [34, 76, 77, 85, 87, 114–117]. The numerical errors in these methods can either be inherited from the aforementioned algorithms to solve the one-electron Kohn–Sham (or Hartree–Fock) equations, or arise from extra algorithms, such as the choice of the self-consistent Kohn–Sham orbitals for the post-processing evaluations [118], the resolution-of-identity technique to handle the two-electron four-center integrals [18, 119–121], and the localization approximations [114, 122, 123] to reduce the computational scaling in these advanced correlation methods.

In this context, it becomes imperative for computational materials science to have a representative test set with numerically well-converged reference values at various levels of theory. In spirit of the existing quantum-chemistry test sets and the Δ-value concept, we introduce in the following our main group test set for MSE, coined MSE, which is based on results acquired using density functional methods from LDA, PBE, PBEsol, SCAN, TPSS, M06-L, HSE06; and state-of-the-art MP2 and RPA. The numerical convergence of these methods is investigated in terms of the total energy. Cohesive energies, lattice constants and bulk moduli are then reported. A comprehensive understanding of the numerical errors, particularly in MP2 and RPA, is discussed in order to aid the community's pursuit of a numerically stable implementation of CCSD(T) and full-CI for the exact description of solids.

At the MSE home page, an overview of the MSE test set project is given along with a table that summarizes all materials available for a given DFA. The materials are sorted by their crystalline structures with a layout similar to figure 1. The default DFA is PBE; the table for other DFAs is obtained via a drop-down select box. A search engine allows to access a group of results for a given material, structure, or DFA, and/or combination of the above. Figure 2 shows a snapshot of the filtered table produced by the search engine, listing the MP2 reference data currently available. This multi-mode search framework, together with a well-organized data structure, guarantees easy access to the more than 10 000 calculations in the current MSE test set.

The MSE web site provides a quick and easy-to-use visual display of the test set data. For any individual reference value calculated for a given material and using a specific DFA, the crystal structure and the well-converged equation of state are displayed. If available, one can also visualize the numerical convergence towards this reference data; this includes convergence tests for basis set, ${\boldsymbol{k}}$-mesh, and internal numerical parameters. Furthermore, a statistic analysis of the basis set convergence for the whole test set can be visualized, which is derived using the mean absolute deviation (MAD).

As an example of this functionality, figure 3 presents the basis set incompleteness error in MP2 cohesive energies per atom for (A) Si diamond and (B) 14 materials representing covalent, ionic, vdW, and mixed bonding characters. Clearly, NAO-VCC-nZ basis sets provide an improved description for advanced correlation methods with the increase of the basis set size, which allows for the extrapolation to the CBS limit from NAO-VCC-3Z and 4Z (i.e. CBS[3,4]). As a recommendation for practical use, we note that the basis set extrapolation from NAO-VCC-2Z and 3Z (CBS[23]) guarantees a near NAO-VCC-4Z quality, but with much less computational cost.

Zoom In Zoom Out Reset image size

For the benchmark studies of the DFAs, the MAD and the root-mean-square deviation (RMSD) are often used to quantify the numerical errors in calculations. For a given group of materials with the number of materials of n and the targeted observables of {Y_i}, if the computational data are {y_i} with the errors distributed as ${x}_{i}={Y}_{i}-{y}_{i}$, we then define the relevant MAD and RMSD as

$$\begin{eqnarray}&&{\rm{MAD}}=\displaystyle \frac{1}{n}\underset{i=1}{\overset{n}{\displaystyle \sum }}\left|{x}_{i}\right|{\rm{;}}\,{\rm{RMSD}}=\sqrt{\displaystyle \frac{{\displaystyle \sum }_{i}^{n}{x}_{i}^{2}}{n}}.\end{eqnarray} \tag{ 1 }$$

Despite both MAD and RMSD being measures of accuracy, we note that larger errors have a disproportionately large effect on RMSD, making it more sensitive to outliers.

A linear regression by means of a least-squares fit allows us to separate the error into predictable (or systematic) and material-specific (or residual) parts [125, 126]. The resulting linear model

$$\begin{eqnarray}&&\mathop{L}\limits^{\unicode{x002C6}}({y}_{i})=\beta {y}_{i}+\alpha \end{eqnarray} \tag{ 2 }$$

allows the analysis of the material-specific deviations {_i} and the corresponding root-mean-square deviation ($\overline{{\rm{RMSD}}}$) as

$$\begin{eqnarray}&&{Y}_{i}-\mathop{L}\limits^{\unicode{x002C6}}({y}_{i})={\epsilon }_{i}{\rm{;}}\overline{{\rm{RMSD}}}=\sqrt{\displaystyle \frac{{\displaystyle \sum }_{i}^{n}{\epsilon }_{i}^{2}}{n}}.\end{eqnarray} \tag{ 3 }$$

The systematic error can be determined by the difference between ${\rm{RMSD}}$ and $\overline{{\rm{RMSD}}}$. It was argued that the material-specific deviations represent a true measure of the inadequacy of the method or the numerical incompleteness of the basis set and ${\boldsymbol{k}}$-mesh [125, 126]. To make the most of this analysis method, we have equipped the MSE web site with a linear-regression analysis tool. In three examples, we will demonstrate how to gain insight into the numerical and intrinsic limitations of a state-of-the-art DFAs.

3.3.1. Basis set convergence in the MP2 lattice constants

Slow basis set convergence is a well-documented problem for advanced correlation methods like MP2 and RPA. Figure 3 presents the basis set incompleteness errors in MP2 cohesive energies using the valence correlation consistent basis set in the NAO framework. Taking the MP2 results in the CBS limit as the reference, table 1 lists the RMSDs and $\overline{{\rm{RMSD}}}{\rm{s}}$ of the incompleteness errors for NAO-VCC-nZ with n = 2, 3, and 4. Despite the basis set errors systematically decreasing with the increase of the basis set size, RMSD and $\overline{{\rm{RMSD}}}$ with NAO-VCC-3Z remain about 0.1 Å, which is unacceptable for real applications.

Table 1.  RMSDs and material-specific $\overline{{\rm{RMSD}}}$ of the basis set incompleteness errors for the MP2 lattice constants. The values in parentheses exclude Ne and Ar. (unit: Å)

Method	RMSD	$\overline{\mathrm{RMSD}}$
N2Z	0.148 (0.043)	0.126 (0.015)
N3Z	0.089 (0.017)	0.071 (0.009)
N4Z	0.032 (0.016)	0.021 (0.007)

Figure 4 shows the linear regression analysis. It clearly suggests that there are two outliers, which are Ne and Ar with fcc structure. We note that Ne and Ar crystals are bonded by vdW forces, which are overestimated by the MP2 method; however it is an intrinsic error of the MP2 method and thus not addressed here. For isolated molecules, it has been discussed in previous literature [18, 24] that an accurate description of weak interactions using advanced DFAs needs either a counterpoise correction scheme or very large basis sets with diffuse basis functions to address the notorious basis set superposition error. As suggested by our results in this work, this conclusion is also true for solids governed largely by vdW interactions. Excluding these two cases from the linear regression analysis results in a much better quantitative accuracy (values in parentheses in table 1). The RMSD for NAO-VCC-3Z and $\overline{{\rm{RMSD}}}$ for NAO-VCC-2Z are below 0.02 Å, which clearly suggests that NAO-VCC-3Z or even 2Z with a systematic correction would be good enough to converge the MP2 lattice constants for strongly bonded systems, while the weak interactions must be treated with special care.

Zoom In Zoom Out Reset image size

3.3.2. Starting point (SP) influence on meta-GGAs and RPA

The common practice of evaluating meta-GGAs, like TPSS, M06-L and SCAN, and advanced many-body perturbation methods, like RPA, is by performing an energy evaluation a posteriori using LDA, GGA or hybrid GGA orbitals [35, 83]. Meanwhile, the Hartree–Fock orbitals are used for the MP2 method generally in quantum chemistry. The self-consistent implementation of meta-GGAs has been realized in many computational packages, including FHI-aims. In the MSE test set, we provide accurate reference data of TPSS, M06-L and SCAN for both (a) self-consistent and (b) non-self-consistent using PBE orbitals, calculations. We also examine the starting point influence on RPA cohesive energies by using PBE and PBE0 orbitals. For simplicity and consistency, we use 'DFA@SP' in this discussion to denote which method is based on which starting point (SP) orbitals.

Table 2 shows the influence of the starting point on the cohesive energies per atom in terms of the RMSD. The linear regression was performed to extract the material-specific part of the deviations. Clearly, the starting point influence is mild for meta-GGAs: SCAN shows the smallest influence, leading to a RMSD of only 15 meV. Meanwhile, the linear regression analysis suggests that these errors are almost material-specific. The corresponding systematic errors (${\rm{RMSD}}-\overline{{\rm{RMSD}}}$) are less than 2 meV. In contrast, RPA is much more sensitive to the choice of the starting point orbitals, though the influence is quite systematic: the material-specific error $\overline{{\rm{RMSD}}}$ is only about 30 meV. This small error indicates that a careful choice of the starting point orbitals can improve the RPA performance because of its systematic underestimation of cohesive energies [118] and weak interactions [34].

Table 2.  RMSDs and material-specific $\overline{{\rm{RMSD}}}$ of the starting-point influence on cohesive energies per atom. 'DFA@SP' denotes which method is based on which starting point orbitals. 'sc-DFA' denotes the self-consistent study. (unit: eV)

Methods to compare	RMSD	$\overline{{\rm{RMSD}}}$
sc-SCAN versus SCAN@PBE	0.015	0.015
sc-M06-L versus M06-L@PBE	0.049	0.043
sc-TPSS versus TPSS@PBE	0.035	0.033
RPA@PBE versus RPA@PBE0	0.119	0.030

For meta-GGAs and RPA, the starting point influences on lattice constant and bulk modulus are similar: the material-specific deviations are mild, though the systematic deviation is quite large for RPA results. The readers interested in this topic can easily access the data and perform the linear regression online by themselves. It is also easy to investigate the influence of relativistic effects on different methods and properties in the same manner.

3.3.3. Cross-over comparison between different methods

Figure 5 summarizes the cross-over comparison of cohesive energies between different methods. RMSDs and the material-specific $\overline{{\rm{RMSD}}}{\rm{s}}$ are shown in different subtables. In this comparison, RPA@PBE reference data were used and meta-GGAs were calculated self-consistently. The nine DFAs investigated cover all five rungs of the 'Jacob's ladder of DFT'. From the many-body perturbation theory point of view, MP2 and RPA consider the electronic correlations explicitly in the many-body interaction picture, while others are all mean-field approximations.

Zoom In Zoom Out Reset image size

LDA and RPA are the density functional approximations on the lowest and highest rungs of the 'Jacob's ladder', respectively, and they show the largest deviation in the cohesive energy calculations, leading to a RMSD of over 1 eV. In fact, LDA shows a large disagreement with all other methods, with the average RMSDs about 0.7 eV. This observation agrees with a widely accepted argument that the local density approximation derived from homogeneous electron gas completely misses the high order density derivative or response information that is important for a proper description of chemical bonding. However, our linear regression analysis suggests that this error is quite systematic: the material-specific $\overline{{\rm{RMSD}}}$ is only 0.16 eV between LDA and RPA.

A similar analysis can be applied to study the intrinsic limitations in MP2 and RPA. It is well-known that the RPA method performs unsatisfactorily in the calculations of cohesive energy for solids [35, 37, 127, 128]. From the many-body perturbation point of view, this error is due to the lack of an infinite summation of the second-order exchange diagrams, which is necessary to eliminate the notorious one-electron self-interaction error in the RPA method [35, 37]. On the other hand, the MP2 method itself is exactly free from such self-interaction error, but completely ignores high-order perturbative diagrams. In consequence, a large RMSD of about 0.68eV between MP2 and RPA cohesive energies is observed, which can be traced back to the dissimilarity of underlying physical models used in the two methods. Recently, much effort has been devoted to improve RPA and MP2 from different perspectives, but these approaches often exacerbate the computational complexity dramatically. Our linear regression analysis suggests that the difference between MP2 and RPA cohesive energies is also quite systematic, leading to a material-specific deviation $\overline{{\rm{RMSD}}}$ of only about 100 meV.

We also compared the results of different methods with experimental data that were corrected for thermal and zero-point vibrational effects [83, 107]. In-line with the observation from the above cross-over comparision, most of the errors in different theretical approximations are quite systematic, leading to the material-specific errors $\overline{{\rm{RMSD}}}{\rm{s}}$ in 9 methods are lower than 210 meV, with M06-L presenting the largest difference. The visualization of the linear regression analysis suggests that abnormally large errors in M06-L occur in metallic systems (Na and Al), which can be ascribed to an incorrect oscillation of the exchange enhancement factor of M06-L when approaching to the uniform electron gas limit [129]. Taking these two points out, the RMSD and $\overline{{\rm{RMSD}}}$ of M06-L are reduced to 131 meV and 95 meV, respecitvely.

To summarize this section, we have presented three examples that demonstrate the usage of our MSE web site. In the context of on-going innovation in computational materials science, driven by data technology, there is a growing awareness of the importance of effective data sharing and recycling. Here, we argue that a dedicated test-set web site should be about more than an easy, static access of the reference data. In order to liberate the power of test sets, of key importance is providing a friendly analysis interface that facilitate the users to play with the online data, repeat the observation in the original papers, and even gain new insights from the data by themselves.

In the framework of the mean-field approximations, the evaluation of the XC energy needs only the electron density and/or occupied orbitals. However, the advanced correlation methods, i.e. MP2 and RPA, explicitly consider the excitations that are between occupied and virtual orbitals and that cross-over different ${\boldsymbol{k}}$ points. Therefore, it is not surprising to see that MP2 and RPA calculations for solids have a much higher demand on the sampling of the first Brillouin zone. The ${\boldsymbol{k}}$-mesh convergence of the periodic MP2 method has been investigated with the Γ-center sampling strategy by Grüneis et al [76, 130]. In FHI-aims, we adopt a hybrid strategy to generate the Coulomb matrices in MP2 and RPA calculations, i.e. using the cut-Coulomb operator only for the Γ point and the full Coulomb operator for the rest of the ${\boldsymbol{k}}$ points. For MP2 and RPA, this hybrid choice shows a faster ${\boldsymbol{k}}$-mesh convergence compared to the cut-Coulomb strategy used in HSE06 calculations [105].

Taking C diamond and MgO rocksalt as examples, figure 6 shows the ${\boldsymbol{k}}$-mesh incompleteness error of total energies per atom for five methods:

$$\begin{eqnarray}&&{\rm{\Delta }}{E}_{{\rm{total}}}[{n}_{{\boldsymbol{k}}}]={E}_{{\rm{total}}}[{n}_{{\boldsymbol{k}}}]-{E}_{{\rm{total}}}[\infty ],\end{eqnarray} \tag{ 4 }$$

where ${n}_{{\boldsymbol{k}}}$ is the number of ${\boldsymbol{k}}$ points in each direction (${n}_{{\bf{k}}}\times {n}_{{\bf{k}}}\times {n}_{{\bf{k}}}$). Our results reveals that the total energies per atom calculated at MP2 and RPA levels converge dramatically slower than those of LDA, PBE and HSE06. The ${\boldsymbol{k}}$-mesh incompleteness error in C diamond remains to be 15 meV for both MP2 and RPA total energies with a dense 8 × 8 × 8 ${\boldsymbol{k}}$-mesh. We notice that other off-set evenly sampled ${\boldsymbol{k}}$ meshes [96, 97], the so-called mean-value point strategy [98], and smearing sampling methods, such as the method of Methfessel and Paxton [99], the improved linear tetrahedron method [100], etc., have all been proposed to provide a better sampling with fewer (or even one) ${\boldsymbol{k}}$ points, but at the semi-local DFT level of theory only. It is interesting to investigate the influence of these sampling strategies on advanced correlation methods.

Zoom In Zoom Out Reset image size

The computational expense of the canonical MP2 (and RPA) scales as ${N}_{{\boldsymbol{k}}}^{3}$ (and ${N}_{{\boldsymbol{k}}}^{2}$), where ${N}_{{\boldsymbol{k}}}={n}_{{\boldsymbol{k}}}^{3}$ is the total grid number in a given ${\boldsymbol{k}}$ mesh. For C diamond and NAO-VCC-2Z basis set (with 28 basis functions per unit cell), the MP2 calculation with the 8 × 8 × 8 ${\boldsymbol{k}}$-mesh requires 10 hours to complete using 320 CPU cores of an Infiniband-connected Intel cluster with Intel Xeon E5-2680 v2 cores (2.8 GHz, 20 cores per node). The periodic MP2 implementation in FHI-aims is a ${\boldsymbol{k}}$-mesh-oriented parallelization in the framework of the message passing interface, guaranteeing efficient MP2 calculations with thousands of cores. However, with the optimistic assumption of a perfect parallel scalability, it would need over 2000 CPU cores to finish the MP2 calculation with a denser ${\boldsymbol{k}}$-mesh (10 × 10 × 10) in a comparable amount of time. While the RPA calculations with such a dense ${\boldsymbol{k}}$-mesh and the NAO-VCC-2Z basis set can be carried out at a reasonable cost, the remaining ${\boldsymbol{k}}$-mesh error is about 7 meV for C diamond. The RPA calculations with denser ${\boldsymbol{k}}$-meshes and larger basis sets become infeasible as well.

A practical way to approach the complete ${\boldsymbol{k}}$-mesh limit (CKM) for advanced correlation methods is the two-point extrapolation in terms of an inverse relation between ${\rm{\Delta }}{E}_{{\rm{total}}}[{n}_{{\boldsymbol{k}}}]$ and ${n}_{{\boldsymbol{k}}}$:

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}{E}_{{\rm{total}}}[{n}_{{\boldsymbol{k}}}] & = & A/{n}_{{\boldsymbol{k}}}^{\alpha }\\ \mathrm{log}\left({\rm{\Delta }}{E}_{{\rm{total}}}[{n}_{{\boldsymbol{k}}}]\right) & = & \mathrm{log}\left(A\right)-\alpha \mathrm{log}\left({n}_{{\boldsymbol{k}}}\right)\end{array}\end{eqnarray} \tag{ 5 }$$

The exponential 'α' determines the speed of the ${\boldsymbol{k}}$-mesh convergence of different methods. On a logarithmic scale (see equation (5) and figure 6), 'α' is the negative value of the slope of ${\rm{\Delta }}{E}_{{\rm{total}}}[{n}_{{\boldsymbol{k}}}]$ towards the CKM limit. Previously, the exponential α = 3.00 has been successfully used to extrapolated MP2 and CCSD(T) cohesive energies for several materials [85].

Table 3 shows the performance of such ${\boldsymbol{k}}$-mesh extrapolation for five selected materials (two elemental and three binary solids) from the MSE test set. The complete ${\boldsymbol{k}}$-mesh extrapolation from 8 × 8 × 8 and 10 × 10 × 10, CKM(8, 10), are taken as the reference for the periodic RPA calculations. The deviations of the CKM(6, 8) values are smaller than 1.5 meV for all five materials and the MAD is only 1.1 meV. Our results confirm the validity of the practical extrapolation (equation (5)) with α = 3.00 for advanced methods provided that reasonable dense ${\boldsymbol{k}}$-meshes are utilized. However, a notable decrease of the performance can be observed in the combination of CKM(4, 6) and α = 3.00. The MAD remains at about 8.0 meV (with the maximum error of 14 meV for C diamond) for the RPA total energies per atom. In this project, we optimize 'α' for the CKM(4, 6) extrapolation, and find that α = 3.95 improves the CKM(4, 6) extrapolation consistently, resulting in the MAD of only 1.6 meV.

The CKM(6, 8) values with α = 3.00 are chosen as the reference for MP2. Table 3 indicates that the CKM(4, 6) extrapolation is not good enough when using the default choice of α = 3.00, and it can be effectively improved with an exponential value α = 3.95 as optimized for RPA. Our results reveal that advanced correlation methods share a similar ${\boldsymbol{k}}$-mesh convergence behavior, and thus the CKM(6, 8) values with α = 3.00 are good enough to converge the error within 2.0 meV for insulators and semiconductors, which is choosen for both MP2 and RPA calculations in this work (see table 1 in supporting information).

We notice that the periodic MP2 and CCSD(T) results reported in the literature were often extrapolated using CKM(5, 6) or even sparser CKM(3, 4) with α = 3.00 [85]. It is an understandable compromise considering the computational expense with denser ${\boldsymbol{k}}$-meshes and larger basis sets for these methods; however, in order to achieve ${\boldsymbol{k}}$-mesh-converged values with the numerical uncertainty within 10 meV, an adapted exponential factor for CKM extrapolation appears necessary from our benchmarking.

The NAO basis sets of NAO-VCC-nZ, used in FHI-aims, hold the promise to provide numerically well-converged total energies for advanced correlation methods in frozen-core approximations. (see [24] for details). For MP2 and RPA calculations, the complete basis set (CBS) results are extrapolated using a two-point extrapolation formula [24, 60] based on NAO-VCC-3Z and NAO-VCC-4Z, namely CBS(3, 4):

$$\begin{eqnarray}&&E[\infty ]=\displaystyle \frac{E[3]{3}^{3}-E[4]{4}^{3}}{{3}^{3}-{4}^{3}}\end{eqnarray} \tag{ 6 }$$

Here E[n] is the RPA or MP2 total energy per atom using NAO-VCC-nZ basis sets (with n = 3 and 4). While the extrapolation from NAO-VCC-4Z and NAO-VCC-5Z shall further reduce the numerical uncertainty of the reference data [24], NAO-VCC-5Z is too expensive to be used for periodic MP2 and RPA calculations at present.

The ${\boldsymbol{k}}$-mesh convergence of MP2 and RPA methods tested in the previous section is carried out together with the smallest NAO-VCC-2Z basis sets (see figure 6). However, it would be impossible to perform MP2 calculations with an 8 × 8 × 8 k-mesh for NAO-VCC-4Z, and sometimes even a NAO-VCC-3Z calculation is too expensive, because of the unfavorable scaling of the MP2 method with respect to the number of ${\boldsymbol{k}}$ points.

Taking C diamond as an example, table 4 shows the ${\boldsymbol{k}}$-mesh convergence of the MP2 total energies per atom for different basis sets. The most time consuming result in this table is the combination of NAO-VCC-3Z and the 8 × 8 × 8 ${\boldsymbol{k}}$-mesh, which requires about 4 days using 320 CPU cores of an Infiniband-connected Intel cluster with Intel Xeon E5-2680 v2 cores (2.8 GHz, 20 cores per node). As illustrated in table 4, while the error of a 4 × 4 × 4 k-mesh varies with the basis set employed by several meV, it becomes almost independent of the choice of basis sets for 6 × 6 × 6 and 8 × 8 × 8 meshes. Ohnishi et al [131] and Grüneis et al [132] had a similar observation in MP2 and CCSD calculations using the GTO-type and plane-wave basis sets, respectively. In consequence, they concluded that the long-range behavior of the correlation energy depends mostly on the low-lying excitations, and proposed their progressive down-sampling technique to approach the complete ${\boldsymbol{k}}$-mesh and basis set limit for advanced correlation methods.

Table 4.  Errors in the MP2 total energies per atom using different basis sets for C diamond (in meV). The reference is the CKM(6, 8) extrapolated values with α = 3.00 (see equation (5)).

${\boldsymbol{k}}$-mesh	NAO-VCC-2Z	NAO-VCC-3Z	NAO-VCC-4Z
(4 × 4 × 4)	117	114	115
(6 × 6 × 6)	22	22	22
(8 × 8 × 8)	9	9	9a

^aFor NAO-VCC-4Z, the MP2 total energy per atom with the 8 × 8 × 8 ${\boldsymbol{k}}$-mesh is extrapolated using the energy difference between ${\boldsymbol{k}}$-grids 6 × 6 × 6 and 8 × 8 × 8 with the NAO-VCC-3Z basis set.

In this work, we adopt a similar down-sampling concept to estimate the MP2 total energies per atom with the NAO-VCC-4Z basis set and the 8 × 8 × 8 ${\boldsymbol{k}}$-grid: ${E}_{\mathrm{total}}^{{\rm{4Z}}}[{n}_{{\boldsymbol{k}}}=8]$

$$\begin{eqnarray}&&{E}_{\mathrm{total}}^{{\rm{4Z}}}[{n}_{{\boldsymbol{k}}}=8]={E}_{\mathrm{total}}^{{\rm{4Z}}}[{n}_{{\boldsymbol{k}}}=6]+{\rm{\Delta }}{E}_{\mathrm{total}}^{n{\rm{Z}}}[{n}_{{\boldsymbol{k}}}=8]\end{eqnarray} \tag{ 7 }$$

using the energy change between 6 × 6 × 6 and 8 × 8 × 8 k-meshes with the NAO-VCC-3Z (or 2Z) basis set ${\rm{\Delta }}{E}_{\mathrm{total}}^{n{\rm{Z}}}[{n}_{{\boldsymbol{k}}}=8]$ (n = 3 or 2) of:

$$\begin{eqnarray}&&{\rm{\Delta }}{E}_{\mathrm{total}}^{n{\rm{Z}}}[{n}_{{\boldsymbol{k}}}=8]={E}_{\mathrm{total}}^{n{\rm{Z}}}[{n}_{{\boldsymbol{k}}}=8]-{E}_{\mathrm{total}}^{n{\rm{Z}}}[{n}_{{\boldsymbol{k}}}=6]\end{eqnarray} \tag{ 8 }$$

In summary, the reference data for the MP2 and RPA total energies (and also valence correlation energies) per atom are obtained by the combination strategy of CKM(6, 8), CBS(3, 4) and ${\rm{\Delta }}{E}_{\mathrm{total}}^{{\rm{nZ}}}[{n}_{{\boldsymbol{k}}}=8]$ with (n = 3 or 2). The error bar is estimated to be 20 meV for these MP2 and RPA reference data. The numerical uncertainty in the combination strategy is dominated by the CBS(3, 4) extrapolation, i.e. 15 meV for CBS(3, 4), 2 meV for CKM(6, 8), and less than 1 meV for ${\rm{\Delta }}{E}_{\mathrm{total}}^{{\rm{3Z}}}[{n}_{{\boldsymbol{k}}}=8]$ or ${\rm{\Delta }}{E}_{\mathrm{total}}^{{\rm{2Z}}}[{n}_{{\boldsymbol{k}}}=8]$, respectively.

Based on the efforts discussed above to achieve the complete ${\boldsymbol{k}}$-mesh and basis set convergence in terms of absolute total energy for various levels of theory, we next focus on the numerical convergence of MP2 and RPA calculations on lattice constants, bulk moduli and cohesive energies.

The post-processing MP2 and RPA correlations are frozen-core, but evaluated based on the Hartree–Fock, PBE and/or PBE0 orbitals in all-electron description. The basis set convergence of these properties are well-documented for (semi-) local and hybrid DFAs. However, the slow basis set convergence together with an unfavorable computational scaling makes it a big challenge to perform a systematic investigation of the basis set convergence of these materials' properties at MP2 and RPA levels. We present here the MP2 results of five selected materials from the MSE test sets in table 5. To the best of our knowledge, it is the first report of the basis set convergence of MP2 for condensed matter systems using NAO basis sets. The readers can easily access the MP2 and RPA data on convergence test for all 13 materials in our MSE web site.

Table 5 reveals that the basis set convergence is relatively fast in the MP2 calculations of lattice constants. NAO-VCC-3Z is good enough to provide a well-converged lattice constant (the convergence criterion is 0.01 Å) for all five materials (for details see also the first example in Section 3). In agreement with the previous investigations using plane-wave basis sets [76, 85] or a Gaussian and plane waves hybrid approach [117], the slow basis set convergence of MP2 cohesive energies is observed with NAO basis sets as well. Compared with CBS(3, 4) values, the basis set incompleteness errors at NAO-VCC-4Z still remain about 60-130 meV. Such slow convergence arises from the slower basis set convergence of the MP2 total energies in bulks than that in free atoms. It is well-documented in quantum chemistry that the accurate geometry information for the MP2 and other advanced correlation methods can be obtained with the basis sets of triple-zeta quality [133], but the converged atomization energy or other energy differences cannot be achieved with finite basis sets [41, 134, 135]. Our results confirm that this conclusion is also valid for solids. While the bulk moduli obtained at the NAO-VCC-4Z basis set are not well-converged, the basis set error remains about 3–7 GPa.

Inspecting table 5 also reveals the capability of NAO-VCC-nZ to provide a consistently improvable description of cohesive properties. In general, the calculated MP2 cohesive energies increase with the cardinal index (n) of NAO-VCC-nZ. On the one hand, it allows for an accurate extrapolation to the CBS limit; and, on the other hand, the NAO-VCC-4Z values set up a quite rigorous lower bound for the converged values of MP2. In other words, any results smaller than the NAO-VCC-4Z values might be unreliable.

The canonical MP2 method was implemented in VASP for periodic systems in 2009 [130]. Later, Grüneis et al calculated the MP2 cohesive properties of 13 materials in 2010 [76]. To approach the CBS limit in the framework of the projector-augmented-wave (PAW) method using plane-wave basis sets, they proposed an inverse relation between the MP2 correlation energy and the energy cutoff E_χ that represents the overlap charge densities for the CBS extrapolation. Meanwhile, an efficient composite scheme was proposed to reduce the computational cost due to the unfavorable scaling of the ${\boldsymbol{k}}$-point number. To be specific, the Hartree–Fock total energy is calculated with a dense 10 × 10 × 10 Γ-centered ${\boldsymbol{k}}$-mesh, but the MP2 second-order direct and exchange terms are evaluated with 6 × 6 × 6 and 3 × 3 × 3 ${\boldsymbol{k}}$-meshes, respectively. In 2013, Booth et al updated the MP2 cohesive energies of three materials [85], including C diamond, BN and AlP with zincblende structure. The main difference is to replace the composite scheme by an extrapolation scheme from 5 × 5 × 5 and 6 × 6 × 6, i.e. CKM(5, 6) with α = 3, to approach the converged ${\boldsymbol{k}}$-mesh sampling. Table 5 also lists the canonical MP2 cohesive properties calculated with VASP. FHI-aims predicts very similar cohesive energies to the VASP values with the ${\boldsymbol{k}}$-grid extrapolation scheme. The difference is about 40 meV for C diamond and AlP zincblende, and 100 meV for BN in zincblende. For MgO and Si, the discrepancy of over 170 meV between FHI-aims and VASP values shall be ascribed to the ${\boldsymbol{k}}$-grid incompleteness of the composite scheme used in the previous VASP calculations [76].

The ${\boldsymbol{k}}$-mesh convergence for LDA and PBE methods is investigated together with Γ-center (GC) and Monkhorst-Pack (MP) sampling strategies [136–138], both of which span the first Brillouin zone in an evenly spaced manner. The MP ${\boldsymbol{k}}$-mesh coincides with the GC one for an uneven number of grid points [137, 138]. In this project, we choose a set of ${\boldsymbol{k}}$-meshes with even numbers of grid points. Our observation is consistent with previous investigations where the MP ${\boldsymbol{k}}$-mesh shows faster convergence behavior than the GC ${\boldsymbol{k}}$-mesh for insulators and semiconductors at the semilocal density functional level of theory [137, 138]. In order to achieve 1 meV accuracy, the 4 × 4 × 4 MP grid is sufficient for rare-gas crystals (Ne and Ar with fcc structure), and the 6 × 6 × 6 MP grid for other ionic- and covalent-bond systems. To achieve the same accuracy, an 8 × 8 × 8 (sometimes even denser) grid is often necessary with the GC ${\boldsymbol{k}}$-mesh. However, the slow ${\boldsymbol{k}}$-mesh convergence in the calculation of metallic systems is more serious, and cannot be improved using the MP ${\boldsymbol{k}}$-mesh. Based on LDA, PBE, and three meta-GGAs convergence benchmarks, we summarize the ${\boldsymbol{k}}$-mesh setting in table A1, recommended for all semi-local DFAs investigated in this paper: LDA, PBE, PBEsol, TPSS, M06-L and SCAN.

Table A1.  The ${\boldsymbol{k}}$-meshes that are used for different methods. 'GC' denotes an evenly spaced Γ-center ${\boldsymbol{k}}$-mesh and 'MP' for the Monkhorst-Pack one. For mean-field approximations, the convergence criteria is 1 meV (or tighter) in the calculations of the total energy per atom. For MP2 and RPA, the complete ${\boldsymbol{k}}$-mesh extrapolation from 6 × 6 × 6 and 8 × 8 × 8, denoted as CKM(6, 8) is adopted with α = 3 (see equation (5) in the main text).

		(semi-)local DFAs		Hybrid DFAs		MP2 and RPA
	Structure	Type	${\boldsymbol{k}}$-points	Type	${\boldsymbol{k}}$-points	Type	${\boldsymbol{k}}$-points
Li	bcc	GC	16 × 16 × 16	GC	16 × 16 × 16	—	—
Na	bcc	GC	16 × 16 × 16	GC	16 × 16 × 16	—	—
Al	fcc	GC	20 × 20 × 20	GC	20 × 20 × 20	—	—
Ne	fcc	MP	4 × 4 × 4	GC	6 × 6 × 6	GC	CKM(6, 8)
Ar	fcc	MP	4 × 4 × 4	GC	6 × 6 × 6	GC	CKM(6, 8)
C	diamond	MP	6 × 6 × 6	GC	10 × 10 × 10	GC	CKM(6, 8)
Si	diamond	MP	6 × 6 × 6	GC	10 × 10 × 10	GC	CKM(6, 8)
LiF	rocksalt	MP	6 × 6 × 6	GC	6 × 6 × 6	GC	CKM(6, 8)
LiCl	rocksalt	MP	6 × 6 × 6	GC	6 × 6 × 6	GC	CKM(6, 8)
NaF	rocksalt	MP	6 × 6 × 6	GC	6 × 6 × 6	—	—
NaCl	rocksalt	MP	6 × 6 × 6	GC	6 × 6 × 6	—	—
MgO	rocksalt	MP	6 × 6 × 6	GC	6 × 6 × 6	GC	CKM(6, 8)
MgS	rocksalt	MP	6 × 6 × 6	GC	8 × 8 × 8	GC	CKM(6, 8)
BeS	zincblende	MP	6 × 6 × 6	GC	8 × 8 × 8	GC	CKM(6, 8)
BN	zincblende	MP	6 × 6 × 6	GC	8 × 8 × 8	GC	CKM(6, 8)
BP	zincblende	MP	6 × 6 × 6	GC	10 × 10 × 10	GC	CKM(6, 8)
SiC	zincblende	MP	6 × 6 × 6	GC	10 × 10 × 10	GC	CKM(6, 8)
AlP	zincblende	MP	6 × 6 × 6	GC	8 × 8 × 8	GC	CKM(6, 8)
LiH	rocksalt	MP	6 × 6 × 6	GC	8 × 8 × 8	GC	CKM(6, 8)

A linear-scaling HSE06 method for condensed matter systems has been implemented in FHI-aims [105]. In order to fully utilize the sparsity of the density matrix in real space, this HSE06 implementation needs a projection of density matrix in a Γ-center ${\boldsymbol{k}}$-mesh to the corresponding Born-von-Karman (BvK) supercell by means of Fourier transformation. Furthermore, FHI-aims features a massively parallel, in-memory implementation of canonical MP2 and RPA methods for periodic systems. At present, HSE06, MP2 and RPA calculations can be performed only with the Γ-center ${\boldsymbol{k}}$-mesh in FHI-aims. MP2 and RPA calculations have been reported using ${\boldsymbol{k}}$-meshes centered at Γ-point for the PAW method and plane-wave basis sets [76, 77, 85, 87, 88, 132].

Table A1 lists the ${\boldsymbol{k}}$-meshes for converged HSE06 calculations: generally speaking, HSE06 shows a similar convergence behavior as LDA and PBE using the Γ-center ${\boldsymbol{k}}$-mesh. Notable differences happen only with very rough grids, e.g. in the rare-gas crystals (Ne and Ar) with 2 × 2 × 2 and ionic crystals with 4 × 4 × 4. Such abnormality can be ascribed to an improper description of the integrable singularity of the Coulomb potential in reciprocal space [72, 139], which is mathematically equivalent to a slow decay of Coulomb interaction in condensed matter systems. In FHI-aims, a cut-Coulomb operator is used for HSE06 to generate the Coulomb matrices at every ${\boldsymbol{k}}$-point [105, 112]. However, a reasonably dense ${\boldsymbol{k}}$-mesh is still required. Taking C diamond and MgO rocksalt as examples, figure 6 in the main text shows the ${\boldsymbol{k}}$-mesh incompleteness error of total energies per atom for five methods:

$$\begin{eqnarray}&&{\rm{\Delta }}{E}_{{\rm{total}}}[{n}_{{\boldsymbol{k}}}]={E}_{{\rm{total}}}[{n}_{{\boldsymbol{k}}}]-{E}_{{\rm{total}}}[\infty ],\end{eqnarray} \tag{ A.1 }$$

where ${n}_{{\boldsymbol{k}}}$ is the number of ${\boldsymbol{k}}$ points in each direction (${n}_{{\bf{k}}}\times {n}_{{\bf{k}}}\times {n}_{{\bf{k}}}$). Our results reveal that a Γ-center mesh with 6 × 6 × 6 grid points is enough for HSE06 to address this singularity issue, thus providing a similar ${\boldsymbol{k}}$-mesh sampling quality as LDA and PBE.

The difficulty of using the GTO basis sets, the de facto choice in quantum chemistry, to extended systems has been studied extensively in different kinds of crystals mainly using the CRYSTAL program [101, 140, 141]. To quote the statement in their User's Manual (Chapter 10) [102], 'as a rule, extended atomic basis sets, or 'triple-zeta' type basis sets should be avoided (...), because the outer functions are too diffuse.' On the other hand, despite the predominate use in calculations on periodic systems, the plane-wave basis set is ineffective to describe the localized core-electron states and usually needs to be used in conjunction with pseudo-potentials.

The NAO basis sets of FHI-aims [24, 52], including tier-n and NAO-VCC-nZ, hold the promise to provide numerically well-converged total energies for mean-field approximations in all-electron description, as the minimal basis of both tier-n and NAO-VCC-nZ is composed of the exact core and valence orbitals of spherically symmetric free atoms. In the meantime, a confining potential can be used to generate the primitive NAOs that are exactly localized in a certain region, so that any numerical instability caused by an unnecessary overlap between the diffuse atom-centered basis functions can be significantly suppressed (see [24, 52] for details).

Table A2 lists the basis set incompleteness errors in HSE06 total energies per atom using the tier-n series, which are the default choice for mean-field approximations in FHI-aims. The tier-2 basis set is recommended, and tier-3 is considered to be good enough to reach the CBS limit [52]. While the tier-n basis sets are optimized with respect to small molecules, i.e. symmetric dimers for all elements [52], our results demonstrate their transferability to the extended systems with diverse chemical environments. The average (maximum) error in tier-2 and tier-3are about 8 meV (16 meV) and 2 meV (9 meV) in HSE06 total energies per atom of light elements' materials. Note that the basis size of tier-2 is slightly larger than Dunning's GTO triple-zeta basis set cc-pVTZ [59], which has been recommended to be avoided for calculations on solids [102]. In our work, we do not encounter any self-consistent field (SCF) convergence problem with tier-2 in all our calculations with GGAs and hybrid GGAs, including LDA, PBE, PBEsol, and HSE06. For meta-GGAs, the SCF convergence can be achieved with tier-2 as well, but some code-specific numerical parameters must be set carefully. For instance, we should use a finer numerical integration grid for the vdW systems [142], and it is necessary to introduce a threshold for the kinetic-energy related variable τ, which removes values that do not effect the total energy but do cause sigularities when calculating the potential in the NAO framework [143]. Both parameters have been tested extensively, with values chosen such that the errors these parameters introduce are in the sub-meV regime.

FHI-aims provides a larger tier basis set (tier-4) for some light elements. We notice that for the binaries with mainly ionic bonding characters (LiF—MgS), the SCF procedure with tier-4 can fail. In this work, we report HSE06 total energies per atom using tier-4 for C and Si with diamond structure, SiC and AlP with zincblende structure. Compared with tier-3 results, a noticeable change can be observed only for AlP (about 3 meV). It suggests that tier-3 is competent to provide numerically well-converged HSE06 total energies for most of materials. Concerning the two worst cases of tier-3, i.e. Na bcc and Ar fcc, the basis set incompleteness errors are about 5 meV and 9 meV, respectively. For Na and Ar, the tier-4 basis sets are unavailable. We then introduce 2 s-type and 2 p-type hydrogen-like functions from NAO-VCC-3Z as an additional s, p group to the tier-3, denoted as tier-3+. Table A2 and figure 7 reveal that such s, p group effectively compensates for the basis set incompleteness in tier-3, resulting in the HSE06 total energies per atom with CBS quality for Na bcc and Ar fcc.

Zoom In Zoom Out Reset image size

The NAO basis functions used in the tier-n series are apt to saturate the bonding region in the middle of two atoms, as they are optimized to minimize the LDA total energies of symmetric dimers [52]. In contrast, the sequence of NAO-VCC-nZ basis sets is determined by minimizing the RPA total energies of atoms [24], thus introducing more compact NAO basis functions than the tier-n series. With similar basis sizes, NAO-VCC-3Z delivers a larger basis set incompleteness error (18 meV on average) than tier-2. It indicates the necessity of introducing diffuse atom-centered basis functions in small basis sets to balance the aforementioned incompleteness in core and bonding regions, which renders a better performance for mean-field approximations. However, such discrepancy can be reduced with the increase of the index 'n' in both sequences, the average error being only 2 meV and 3 meV for tier-3 and NAO-VCC-4Z with similar basis sizes.

NAO-VCC-5Z is the largest NAO basis set in FHI-aims, and yields the lowest HSE06 total energies per atom for most of the materials in the test set. However, SCF convergence problem occurs at Al fcc, C diamond, LiF and MgO with rocksalt structure, BN, BP and SiC with zincblende structure. For binary crystals 'AB', we tried a hybrid basis set strategy N(n₁/n₂)Z, which is a short-hand notation of using NAO-VCC-n₁Z for 'A' and NAO-VCC-n₂Z for 'B'. Compared with NAO-VCC-4Z results, better HSE06 total energies per atom can be obtained without SCF convergence problem for LiF, MgO, BP and LiH with N(4/5)Z and for SiC with N(5/4)Z.

Figure 7 presents the differences of the HSE06 total energies per atom (Δ E_total) between the best NAO-VCC-nZ and tier-n basis sets. The discrepancy Δ E_total is very small, about 1 meV on average. Considering that NAO-VCC-nZ and tier-n are designed for completely different purposes, such an agreement with each other clearly indicates that the CBS limit has been approached for the calculations of HSE06 total energies per atom in all-electron description.

Other (semi-)local DFAs share a similar basis set convergence as HSE06. The observation and discussion here are then transferable. In addition to HSE06, one can find a comprehensive basis set convergence test with LDA, PBE, TPSS, M06-L, and SCAN in the MSE web site. For PBEsol, the best NAO basis sets determined by PBE are used to calculate the reference data directly.

To find the optimized lattice parameters for each material and each method, we calculate seven points within a range of ±5% around the initial value of the lattice constant, and fit the respective (volume, energy) data points to the Birch-Murnaghan equation of state to obtain the optimized lattice constant. If all 7 lattice constants used to generate the equation of state lie within a range of ±7% around the optimized value, the value is taken as the final, optimized lattice constant, a₀. Otherwise, this optimized lattice constant is used as a new starting point and more (volume, energy) data points are calculated, so that finally all materials' optimized lattice constants are obtained from an equation of state fitted to seven points with lattice constant values in a range of ±5% to ±7% around the final, optimized lattice constant a₀. We have carefully tested the stability of this procedure using the PBE functional. The obtained optimized lattice constants and bulk moduli are accurate within 0.02 Å and 0.1 GPa, respectively.

The cohesive energies E^M_coh of the materials are then calculated at the optimized lattice constant a₀:

$$\begin{eqnarray}&&{E}_{{\rm{coh}}}^{M}={E}^{M}/{N}_{{\rm{atom}}}-\mathop{\displaystyle \sum }\limits_{{\rm{atom}}}{E}^{{\rm{atom}}}/{N}_{{\rm{atom}}},\end{eqnarray} \tag{ A.2 }$$

where N_atom is the number of atoms in the unit cell, E^M the total energy of the unit cell for material M, and the sum is taken over the total energies E^atom of the constituent atoms in their spin-polarized symmetry-broken ground-state, i.e. no fractional occupancies.