Performance of various density-functional approximations for cohesive properties of 64 bulk solids

Our dataset comprises 64 crystalline solids with cubic structures: 24 metals and 40 semiconducting and insulating solids. The structures cover: A1 (face-centered cubic (fcc), 13 solids), A2 (body-centered cubic (bcc), 11 solids), A4 (diamond, four solids), B1 (rock-salt, 16 solids), and B3 (zincblende, 20 solids), as shown in figure 1. To explore and understand trends, we divide these 64 solids into five classes, as follows (see table 2): main-group metals (MM), transition metals (TM), semiconductors (SC), ionic crystals (IC), and transition metal carbides and nitrides (TMCN). MM includes groups IA and IIA, together with aluminum and lead; TM covers periods 4–7 of the periodic table; SC represents the solids of groups IVA, IIIA–VA, and IIB–VIA; the six ionic crystals are LiCl, LiF, NaCl, NaF, MgO, and MgS; and there are 10 compounds in the TMCN class (TMC and TMN, TM = Ti, Zr, Hf, V, and Nb). The full list of the solids studied here is shown in the supplementary material, with the calculated and experimental cohesive data tabulated in tables SII–SVII, alongside the relevant literature.

Zoom In Zoom Out Reset image size

All calculations were performed employing the FHI-aims package [63], which implements full-potential, all-electron electronic-structure theory with numeric atom-centered orbital basis functions. These basis functions are grouped together as 'tiers' (levels) in FHI-aims in addition to a minimal basis of occupied orbitals of spherically symmetric free atoms. Besides, there are different sets of numerical defaults for each atomic species, namely light, tight, and really tight settings [63]. For numerical accuracy and computational efficiency purposes, tight setting with tier 2 basis sets are recommended for obtaining cohesive properties of bulk solids [20], and are employed throughout this work. The typical tight setting convergence criteria are 10⁻⁵ electrons for the electron density, 10⁻³ eV for the eigenvalues, and 10⁻⁶ eV for the total energy of the system. An meV level of convergence is achieved for the cohesive energy.

For each solid, the equilibrium volume V₀ and bulk modulus B₀ were determined by fitting cohesive energies per unit-cell at 7–11 different points in a range of V₀ ± 20% to a third-order Birch–Murnaghan equation of state (EOS) [64]

$$\begin{eqnarray}&&E(V)={E}_{0}+\displaystyle \frac{9{V}_{0}\,{B}_{0}}{16}\left\{{\left({\left(\displaystyle \frac{{V}_{0}}{V}\right)}^{\displaystyle \frac{2}{3}}-1\right)}^{3}{B}_{0}^{\mathrm{{\prime} }}+{\left({\left(\displaystyle \frac{{V}_{0}}{V}\right)}^{\displaystyle \frac{2}{3}}-1\right)}^{2}\left(6-4{\left(\displaystyle \frac{{V}_{0}}{V}\right)}^{\displaystyle \frac{2}{3}}\right)\right\},\end{eqnarray} \tag{ 1 }$$

where V₀, E₀, B₀, and B₀' correspond to the equilibrium unit-cell volume, the cohesive energy, the bulk modulus at V₀, and the pressure derivative of the bulk modulus, respectively. The unit-cell volumes V are related to the lattice constants a by the expressions V = a³/4 for the fcc, diamond, and zincblende structures, and V = a³/2 for the bcc structure. The cohesive energy, defined as the energy per atom required to atomize the crystal, was evaluated using the LDA, PBE, PBEsol, SCAN, M06-L, and HSE06 functionals from the energies of the crystal and the constituent free atoms. Note that both SCAN and M06-L calculations were performed using the PBE orbitals, as the functionality to perform calculations with these functionals in self-consistent manner is not yet available within FHI-aims. For SCAN, the effect of non-self-consistency has been assessed by Tran et al [19], and they suggest the difference between self-consistent and non-self-consistent calculations of hard solids should be in most cases below 0.005 Å for the lattice constant and 50 meV/atom for the cohesive energy (see supplementary material).

The spin-restricted formalism was employed for bulk crystals, with two exceptions (the ferromagnetic metals: Fe and Ni), and the spin-unrestricted formalism for open-shell atoms (no fractional occupancies were allowed). The basis sets and k-point meshes in reciprocal space are fully converged. Here we used Γ-centered Monkhorst–Pack [65] k-grids, and the settings for most calculations are outlined in the supplementary material. Relativistic effects are more important for heavy elements; here for consistency a scalar-relativistic treatment using the scaled zero-order regular approximation [66] was employed for all solids. It should be noted that spin–orbit coupling, which can be important in heavy solids, is not included. However, Philipsen and Baerends [67] have compared the scalar-relativistic and fully relativistic calculations of the cohesive properties of 11 solids in four columns of the periodic table and determined that spin–orbit coupling effects barely altered lattice constants, cohesive energies, and bulk moduli, with mean absolute contributions of about 0.005 Å, 0.03 eV, and 1 GPa, respectively. They showed that the spin–orbit coupling contributions were only non-negligible for the cohesive energies and the lattice constants of Au and Bi. Therefore, the scalar-relativistic treatment is sufficient to capture the relativistic effects in the solids considered here.

vdW interactions arise from the correlated motion of electrons and constitute a large part of the long-range electron correlation energy. While most commonly used DFAs lack the ability to describe the long-range correlation energy, they provide a route to the inclusion of the vdW energy in an effective way. Numerous promising ideas have been proposed for including the vdW energy in DFT motivated by the need to treat large and complex molecular and condensed-matter systems [68–71].

One of the most widespread approaches is the pairwise-additive fragment-based one. According to second-order perturbation theory, the vdW energy of a collection of fragments (atoms or molecules) can be expressed as sum of attractive contributions, with the leading nonvanishing term being due to dipole–dipole interactions

$$\begin{eqnarray}&&{E}_{\mathrm{vdW}}=-\displaystyle \frac{1}{2}\displaystyle \sum _{i}\displaystyle \sum _{j}{f}_{\mathrm{damp}}\displaystyle \frac{{C}_{6}^{{ij}}}{{R}_{{ij}}^{6}},\end{eqnarray} \tag{ 2 }$$

where R_ij is the interatomic distance between atoms i and j of the system, ${C}_{6}^{{ij}}$ is the corresponding dipole–dipole dispersion coefficient, and f_damp is a damping function that is used to avoid singularities at short distances and ameliorate double-counting of correlation at intermediate distances [68–71]. The resulting vdW energy is then augmented to the standard KS-DFT total energy calculated by employing a local, semi-local or hybrid DFA. In principle, the accuracy of any fragment-based DFT-vdW method, referred to as 'DFT+vdW' here, depends solely on the coefficients used and the choice of the damping function.

The leading-order C₆ term can be expressed as

$$\begin{eqnarray}&&{C}_{6}^{{ij}}=\displaystyle \frac{3}{\pi }{\int }_{0}^{\infty }{\alpha }_{{i}}({\rm{i}}\omega ){\alpha }_{j}({\rm{i}}\omega ){\rm{d}}\omega ,\end{eqnarray} \tag{ 3 }$$

where the integrand contains the isotropic dipolar polarizability as a function of imaginary frequency. Thus, the polarizability is a central quantity to the description of vdW interactions. The method employed here is based upon the TS scheme [72], which computes the polarizability, the C₆ dispersion coefficient and the vdW energy by utilizing Hirshfeld partitioning [73] of the electron density to rescale free-atom values [74] obtained from self-interaction corrected time-dependent DFT calculations. The TS approach yields accurate C₆ coefficients for a set of 1225 small molecules (with an accuracy of 5.5%) due to its ability to capture local environmental effects. However, TS lacks the ability to describe the long-range electrostatic screening that extends beyond the range of the exponentially decaying atomic densities. These effects can be included microscopically by modeling the environment as a dipole field of quantum harmonic oscillators (QHOs) and solving the resulting classical electrodynamics self-consistent screening (SCS) equation [75].

In this model, a molecule (or material) is considered as an arrangement of N atoms, each of which is described by a spherical QHO [71, 75, 76]. For a given atom i in the system, one can write

$$\begin{eqnarray}&&{\alpha }_{{i}}^{{\rm{SCS}}}({\boldsymbol{r}},{\rm{i}}\omega )={\alpha }_{{i}}^{{\rm{TS}}}({\boldsymbol{r}},{\rm{i}}\omega )+{\alpha }_{{i}}^{{\rm{TS}}}({\boldsymbol{r}},{\rm{i}}\omega )\displaystyle \sum _{{i}\ne {j}}^{N}{{\bf{T}}}_{{ij}}{\alpha }_{{j}}^{{\rm{SCS}}}({\boldsymbol{r}}^{\prime} ,{\rm{i}}\omega ),\end{eqnarray} \tag{ 4 }$$

where ${\alpha }_{{i}}^{\mathrm{SCS}}({\boldsymbol{r}},{\rm{i}}\omega )$ is the fully screened atomic polarizability tensor (for a given frequency of the electric field), containing both short-range (via the input atomic polarizability tensor from the TS method) and long-range (via the SCS equation) screening effects. The dipole–dipole interaction tensor, T_ij, has the following expression

$$\begin{eqnarray}&&{{\bf{T}}}_{{ij}}={{\rm{\nabla }}}_{{{\boldsymbol{r}}}_{i}}\otimes {{\rm{\nabla }}}_{{{\boldsymbol{r}}}_{j}}V({r}_{{ij}}),\end{eqnarray} \tag{ 5 }$$

where V(r_ij) is the Coulomb potential between two spherical Gaussian distributions at the distance r_ij, with ${{\boldsymbol{r}}}_{i}$ and ${{\boldsymbol{r}}}_{j}$ representing the QHO positions. The solution of the SCS equation yields the atomic polarizability tensors, from which the isotropic polarizability can be obtained by diagonalizing and tracing the tensor. To calculate the C₆ dispersion coefficient (via equation (3)), the Gauss–Legendre approach is utilized [75].

In this work, we calculate dipolar polarizabilities and C₆ dispersion coefficients for a set of 23 semiconductor solids using the method introduced above, and use them via the DFT+vdW scheme to evaluate the role of vdW interactions in the cohesive properties of solids.

Zero-point vibrational energy effects have been evaluated using the quasi-harmonic approximation. Phonon calculations were performed using phonopy [77] interfaced with the FHI-aims package. The supercell approach [78] was combined with the quasi-harmonic approximation, i.e., for a set of different volumes, harmonic vibrational frequencies are determined using ab initio calculations, with the zero-point energies being added to the DFT ground-state energies at each volume V. In the FHI-aims package, the finite-displacement method is employed to calculate phonon spectra and the vibrational density of states (DOS), g(V, ω), and ZPE is estimated as the frequency integration over the vibrational DOS

$$\begin{eqnarray}&&\mathrm{ZPE}=\displaystyle \frac{1}{2}\int {\rm{\hslash }}\omega g({\text{}}V,\omega ){\rm{d}}\omega .\end{eqnarray} \tag{ 6 }$$

A systematic test of the supercell size was carried out and it was found that a choice of a 2×2×2 supercell with a 8×8×8 k-point grid is sufficient to ensure the convergence of zero-point vibrational energies for all the solids studied here. Within the finite-displacement approach, one may adjust the value of the small finite-displacement δ used to calculate the force response in the three Cartesian directions. This has also been systematically tested, demonstrating that a reasonable variation in δ has only a small impact on the computed ZPE. For example, the ZPE for diamond is about 0.18 eV/atom, and the change in ZPE is less than 1 meV/atom (0.2%) when varying δ from 1×10⁻³ to 1×10⁻¹ Å. In addition, it was found that the effect of the specific DFT functional on the ZPE is negligible for these systems, and as a result, ZPE was calculated at the PBE level and added on top of the ground-state DFT total energies for all six functionals.

For the analysis of our calculations, we consider statistical quantities: the mean error (ME), the mean absolute error (MAE), the mean relative error (%), and the mean absolute relative error (%), all of which are determined by comparing with reliable experimental values measured at low temperatures and/or extrapolated to 0 K.

We start by discussing the contribution of ZPE to the cohesive properties of the solids in our dataset. Figures 2 and 3 show the relative contributions of ZPE to the lattice constants a₀, the cohesive energies E₀, and the bulk moduli B₀ of selected metals and non-metallic solids, respectively, together with available theoretical values reported from previous work [12, 14, 79]. A summary is given in table 3 for all of the solids considered in the present work. The results show that including zero-point vibrations has an even stronger effect on the bulk modulus than on both the lattice constant and the cohesive energy. For 24 metals, the average contributions are 0.2%, 0.7%, and 1.4% for a₀, E₀, and B₀, respectively. In addition, the ZPE effect is, as expected, inversely proportional to nuclear mass (e.g., along the series of Li–Na–Al–K), and is found to be negligible for 'heavy' solids. When only considering the 'light' metals, located on the second to fourth rows of the periodic table, the averaged contribution is almost twice as large as the value calculated for the whole group of metals. In particular, for bulk Li, the lattice constant expands by 0.7% and the bulk modulus reduces considerably by ∼4%.

Zoom In Zoom Out Reset image size

For non-metallic solids, ZPE contributions become more pronounced. A set of 40 solids shows average contributions of 0.3% to the lattice constants, 1.1% to the cohesive energies, and 2.2% to the bulk moduli. For ionic crystals, the effect is particularly strong. The largest contribution occurs for LiF, with changes of about 1% (0.04 Å) for a₀ and 8% (5.6 GPa) for B₀, which are much larger than numerical errors in DFT calculations, and even might be comparable to the contribution due to the long-range electron correlation [80]. This suggests that for non-metallic solids, especially ionic crystals, the inclusion of ZPE becomes vital for an accurate assessment of systems where long-range vdW interactions can be of relevance.

In comparison with existing literature data, our results are in excellent agreement with the calculations reported by Schimka and co-authors [14] using the projector augmented wave method implemented in the Vienna ab initio simulation package [81, 82]. Their dataset comprises 30 solids (including seven metals, six ionic crystals, and 17 semiconductors). For some metals (like Na, Rh, and Cu) and most non-metallic solids, their corrections to both the lattice constants and cohesive energies are essentially identical to the present results, as shown in figures 2 and 3. To the best of our knowledge the evaluation of the ZPE contribution to the bulk moduli has not been performed previously for such a diverse set of bulk solids.

Zoom In Zoom Out Reset image size

The present study confirms that ZPE tends to increase lattice constants and reduce bulk moduli, with its inclusion improving the predictions of the LDA functional and worsening those of PBE (see table 1). As a result, neglecting zero-point vibrations will lead to a bias in the appraisal of these and other DFAs. Therefore, in the following sections, ZPE vibrational contributions will always be included in the discussion of functional performance.

Having assessed the impact of ZPE, we can now study the performance of local, semi-local, and hybrid DFAs for the dataset of 64 solids. The calculated a₀, E₀, and B₀ using the LDA, PBE, PBEsol, SCAN, M06-L, and HSE06 functionals are tabulated in tables SII–SVII of the supplementary material, together with their experimental values, while the MEs and MAEs are shown in table 4. To illustrate the performance of the six functionals for different types of solids, we also plot the relative errors in the bulk moduli and the cohesive energies as a function of those in the lattice constants for the specific groups of solids in figures 4 and 5, respectively. The corresponding MAREs are presented in table 5 in terms of classes of solids, properties, as well as functionals.

Zoom In Zoom Out Reset image size

Let us first discuss the performance of the LDA functional. The present work confirms the well-known trend that LDA underestimates lattice constants, overestimates bulk moduli, and considerably overestimates cohesive energies in most situations. From the MAREs of the three properties summarized in table 5, one can clearly observe that LDA shows a different performance for the different classes of solids. For the SC class, it predicts the lattice constants and the bulk moduli in good agreement with experiment, however at the expense of the cohesive energies, while for ionic crystals it yields moderately accurate results. However, for metals (groups MM and TM) the performance of LDA is unsatisfactory. Note that the performance of LDA for 4d and 5d transition metals is better than 3d and main-group metals.

Turning to the PBE functional, one can see the opposite trend to LDA: lattice constants are too large, and the bulk moduli are too small. Cohesive energies are significantly improved though, with a MARE of ∼6% versus ∼20% for LDA. Overall, PBE improves the agreement with experiment compared to the LDA functional. For the MM and 3d metals, the superiority of PBE over LDA for calculating cohesive properties becomes particularly prominent. For semiconductors and ionic crystals, PBE yields large deviations from experimental results for both the lattice constants and the bulk moduli, while the cohesive energies are reasonably well predicted. A characteristic feature of the PBE functional observed in the current study is the increase of the error with increasing lattice constant, as can be seen, e.g., along the series of C–SiC–Si–Ge–Sn, BN–BP–BAs, and AlP–AlAs–AlSb. For certain systems, the bulk moduli and the cohesive energies follow the same trend. We will return to this observation in section 3.3, when exploring the underlying causes for deviations and trends in the DFAs considered here.

The PBEsol functional was first proposed to improve the description of solids and surfaces compared to PBE [25]. For lattice constants, it does indeed yield improvements, curing the problem of underbinding. In fact, it slightly overbinds bulk crystals. The current study shows significant improvement in the lattice constants using PBEsol for semiconductors, insulators, and fcc metals (except the alkaline-earth elements, Ni, Cu, and Th). Consistently, the bulk moduli are much better predicted with PBEsol compared to PBE for the aforementioned systems. For instance, the MAREs are 0.4% and 0.5% respectively for the lattice constants of the SC and IC classes by PBEsol (versus 1.6% and 1.8% by PBE). The corresponding MAREs for the bulk moduli are very similar for the same solids, both about 6% (versus 14.5% and 16%, respectively, by PBE). For the rest of the solids, i.e., the bcc metals (except K and Rb) and TMCN, both lattice constants and bulk moduli calculated with PBEsol are only marginally better or even worse than with PBE. This agrees with the finding by Janthon and co-authors [83] in a study of 30 transition metals. The substantial change from PBE to PBEsol is a consistent shift of the bonding nature towards tightening: decreasing of the bond lengths (lattice constants) and increasing of the bulk moduli. Overall, PBEsol shows overbinding of solids. One 'expected' problem of using the PBEsol functional, which is observed in the present data, is that the calculated cohesive energies are inaccurate, only halving the LDA error and twice as large as the PBE ones. This can be traced back to inaccurate total energies of the individual atoms [25].

Turning to the meta-GGA functionals, the overall performance of M06-L for cohesive properties lies between that of LDA and PBE. However, no consistent trend can be identified. The current work suggests that M06-L significantly improves the performance on the cohesive properties of ionic crystals in comparison with LDA, GGAs, and hybrid functionals. For the six solids, it slightly overestimates the lattice constants, and reasonably predicts the cohesive energies and the bulk moduli, with resulting MAREs of 0.3%, 3.3%, and 4.9%, respectively. On the other hand, for most semiconductors, M06-L overestimates lattice constants, and underestimates cohesive energies, as well as the bulk moduli. In particular, for main-group metals, M06-L severely overbinds: the lattice constants are too small, and both the cohesive energies and the bulk moduli are too large. For transition metals, M06-L improves the lattice constants of 3d elements in better agreement with experiment, but worsens the other two properties; a significant improvement in cohesive energies is found only for the coinage metals (Pd, Pt, Ag, and Au), however at the expense of both the lattice constants and the bulk moduli. Overall, the M06-L functional performs less accurately than PBE for the TM group, as can be seen from table 5.

The recent SCAN functional is a non-empirical meta-GGA, unlike M06-L, and is the first functional to fulfill all known constraints that a meta-GGA can fulfill. The current study reveals that SCAN predicts remarkably accurate cohesive properties for a diverse set of bulk crystals and tends to strike a good balance between PBE and PBEsol. For non-metallic systems, significant and consistent improvement is found for the predictions of all three cohesive properties using the SCAN functional: the lattice constants are very accurate, with a MARE of less than 0.5% for each class of solids, while both the bulk moduli and the cohesive energies are reasonably well predicted. Our work also suggests that particularly for ionic and semiconductor solids, SCAN remains the best method of choice among the six density functionals tested here, with MAREs within 0.5%, 3%, and 5% for a₀, E₀, and B₀, respectively. For the MM group, SCAN performs equally as well as PBE does for the prediction of the lattice constants, but yields a better prediction of bulk moduli, while for the TM group, the performance of SCAN is as good as PBEsol. Overall, SCAN overestimates the bulk moduli of the bcc metals and underestimates those of the fcc metals. Especially for the 4d and 5d coinage metals (Pd, Ag, Pt, and Au), SCAN cures the underbinding problem of PBE to a large extent. To our knowledge, a theoretical assessment on the reliability of the SCAN functional was first performed by the SCAN authors themselves [34], for the lattice constants of 20 solids. To compare the present results to their data, we consider the same set of 20 solids. The mean absolute deviation and mean deviation of the lattice constants are calculated to be 0.01 Å and 0.002 Å, respectively, between the two studies. This suggests our post-SCF procedure is a good approximation for the current work. A further evaluation of SCAN for other cohesive properties should be performed in the future.

The final DFA tested in the present work is the HSE06 hybrid functional. HSE06 predicts reliable lattice constants and bulk moduli for non-metallic crystals and its performance is only slightly worse than SCAN and PBEsol. For both the ionic and semiconductor solids, the general trend that HSE06 overestimates the lattice constants, underestimates the bulk moduli, and predicts reasonably good cohesive energies is confirmed. While HSE06 improves agreement with experiment for the lattice constants of TMCN, it underestimates the cohesive energies, and considerably overestimates the bulk moduli. Our study indicates that HSE06 exhibits an overall tendency to underbind solids, and this can be pronounced for metals. The overestimation of the lattice constants is found for nearly all of the metals (with the exception of V). Especially for the main-group metals (except Al and Pb), the HSE06 lattice constants are even larger than the PBE values, with a corresponding decrease in bulk moduli. For the TM group, a large difference is seen between the fcc and bcc metals: HSE06 overestimates the bulk moduli for nearly all bcc-TM (except for Fe), whereas it typically underestimates those for fcc-TM (except for Rh). In addition, the cohesive energies show significantly increased errors compared to the PBE results. This can be explained by the fact that hybrid functionals generally increase the stability of the spin-polarized atom, thus reducing stability of bulk systems [48, 49].

Figures 4 and 5 show comparisons of the density functionals tested here for describing the cohesive properties of all five classes of solids. Among the six functionals, LDA predicts accurate lattice constants and bulk moduli for semiconductors, however, the cohesive energy remains a considerable issue at this level, with a MARE of ∼20% for the 24 solids in this group. While PBE outperforms LDA on the prediction of the cohesive energies for semiconductors, it worsens both the lattice constants and the bulk moduli. In fact, the improved performance of PBE upon LDA can be mainly seen in its better description of the cohesion in the systems having a (semi-)metallic nature. As depicted in figures 4 and 5, all the cohesive properties obtained using PBE are in much better agreement with experiment for the MM, TM and TMCN groups, compared to the LDA results. Among the six functionals, PBE remains the best for the 3d transition metals. The PBEsol functional cures the underbinding of PBE, however it somewhat overbinds, but yields reliable lattice constants and bulk moduli for ionic crystals, semiconductors, and fcc metals. Employing PBEsol, the MARE of semiconductors is reduced by a factor of four for lattice constants and two and a half for bulk moduli respectively, compared to the PBE errors. Overall, the PBEsol cohesive energies show systematic overestimation, halving the LDA deviation, with the exception of ionic crystals for which PBEsol yields accurate predictions. The HSE06 functional corrects the lattice constants, and improves agreement with experiment for ionic and semiconductor solids, however it still shows a tendency to underbind. The increased errors in cohesive energies, together with the demanding computational cost, remain significant issues for using the HSE06 functional.

Zoom In Zoom Out Reset image size

In comparison, the M06-L MGGA functional improves on LDA for TM and TMCN (though its performance is slightly worse than PBE), however it yields significant deviations from experiment for the main-group metals, as seen in table 5. This is consistent with the study by Zhao and Truhlar [61]. Compared to PBE, M06-L performs better for semiconductors with lighter elements. For C, Si, BN, BP, and AlP, the cohesive properties obtained from M06-L are in good agreement with experiment. The bulk moduli calculated with M06-L are as good as the LDA values, whereas PBE yields the worst performance among the six functionals, while the M06-L cohesive energies are greatly improved compared to the LDA results. However, when considering heavy-element systems, there is barely any improvement to be gained by using M06-L over PBE. This can be seen for the cases of GaSb, InSb, ZnTe, CdTe, etc. Similar behavior is also observed for ionic crystals. The M06-L functional significantly improves the description for the six ionic compounds that have medium mass, with a MARE of only 0.3% for the lattice constants versus those of 1.5% and 1.9% by LDA and PBE, respectively. Turning to the non-empirical MGGA functional, our study suggests that SCAN can be regarded as the method of choice for non-metallic solids and it agrees remarkably well with experiment for all the three cohesive properties. The SCAN lattice constants are considerably smaller than those of PBE, with an increase of the bulk moduli and larger cohesive energies. For metals, no significant improvement is found on the lattice constants and bulk moduli, while SCAN worsens the cohesive energies compared to PBE. For the 4d and 5d coinage metals, SCAN substantially improves on PBE.

In summary, our study shows that all of the functionals considered have varying performances in predicting the cohesive properties of solids, which strongly depend on the bonding nature of the solid. To understand these trends and behavior in more detail, we now move to discussing the relationships between pairs of cohesive properties, that is, a₀/B₀ and a₀/E₀.

Much attention has been paid in the literature to understanding the relationship between geometrical parameters and properties of various classes solids and the cohesive energy and bulk modulus. It has been found that both E₀ and B₀ are inversely related to a₀ (or the nearest-neighbor distance d) [84–89]. This inverse relation can be explained by the volume dependence of the total energy (the EOS) causing a monotonic decrease of the bulk modulus with increasing volume (see equation (1)). There are several analytical expressions reported in the literature for different families of cubic solids [84–89], e.g., semiconductors (groups IVA, IIIA–VA, and IIB–VIA), ionic crystals (groups IIA–VIA and IA–VIIA), and metals (groups IA and IIA, and noble metals), based upon empirical approaches. While such empirical methods are often not able to yield highly accurate results, they can still be very useful, particularly for illustrating trends in properties of a wide variety of materials.

Anderson and Nafe [84] first proposed an empirical relationship between bulk modulus and crystal volume of the form ${B}_{0}\approx {V}_{0}^{-x}$, where x depends on the bonding type of solid. Cohen and co-authors [85, 86] studied a number of covalent and ionic crystals, and deduced an analytical expression between the bulk modulus and the nearest-neighbor separation d for diamond and zincblende structures

$$\begin{eqnarray}&&{B}_{0}{d}^{3.5}=(1971-220\lambda ),\end{eqnarray} \tag{ 7 }$$

where λ is an empirical ionicity parameter accounting for the ionic character of the bonds. This relationship was found to be appropriate for group IVA (λ = 0), IIIA–VA (λ = 1), and IIB–VIA (λ = 2). For the IA–VIIA rock-salt compounds, Cohen proposed

$$\begin{eqnarray}&&{B}_{0}{d}^{3}=550.\end{eqnarray} \tag{ 8 }$$

A similar scaling of B₀d^3.5 for the rock-salt structure was suggested by Schlosser [87, 88], and the cohesive energy was discussed in terms of d as well. Aresti et al [89] studied the cohesive energy of zincblende solids using the form

$$\begin{eqnarray}&&{E}_{0}d\approx \mathrm{constant}.\end{eqnarray} \tag{ 9 }$$

According to the empirical expressions given by equations (7)–(9), an overestimation (underestimation) of the lattice constant corresponds to an underestimation (overestimation) of the bulk modulus or the cohesive energy. In fact, for the present dataset all of the DFAs considered broadly follow this trend.

Figures 4 and 5 show the relative errors in the bulk moduli and the cohesive energies as a function of the lattice-constant errors, using the LDA, PBE, PBEsol, SCAN, M06-L, and HSE06 functionals for the 64 solids divided into five classes. It can be seen that most data points are in the quadrants II and IV of the Cartesian plane, reflecting that the shorter the lattice constants, the larger the bulk moduli and the cohesive energies, and vice versa. This indicates that those DFAs can predict the observed experimental trends (the inverse relationships of bulk moduli and cohesive energies with lattice constants). If one further looks into the errors concerning a given type of solids, large differences can be found for specific functionals. For SC and IC, a nearly monotonic dependence is observed for the deviations of the bulk moduli from experiment upon those of the lattice constants using all the functionals tested, that is, the overestimated lattice constants is accompanied by the underestimated bulk moduli. It can also be observed that there is a better ('linear') correlation between the errors from the PBE calculations, in comparison with, e.g., LDA, applied to the SC compounds, where half of the values fall in the quadrant III, in disagreement with the empirical study based on experimental observations. Overall, the PBE functional yields the most consistent results with experiment, and it is the method that reproduces systematic trends in the cohesive properties of solids better than the other functionals.

Turning to metals (the MM and TM groups), a noticeable scatter of the errors is found, e.g., from the LDA, PBEsol, and M06-L calculations, with the inversely correlated behavior being somewhat captured by the PBE, HSE06, and SCAN functionals. This also reflects the poor performance of LDA and M06-L, particularly in the description of main-group metals. Finally, for the TMCN group, there is no clear trend shown by any method used here. We note that no clear experimental correlations have been found between any pair of cohesive properties for this class of solids (see figures 6 and 7), mainly due to the mixed nature of the metallic, covalent, and ionic bonding in these solids [90, 91].

Zoom In Zoom Out Reset image size

The above findings suggest that the errors in different cohesive properties obtained from DFT calculations are correlated and system dependent. In general, PBE shows superior performance to LDA, due to the inclusion of the density gradient term for describing nonlocality in realistic systems. The improvement of PBE over LDA is significant for main-group and 3d metals. For some semiconductors and 'heavier' 4d, 5d transition metals, the cohesive properties are better predicted by LDA than by PBE (except for the cohesive energies). When it comes to the M06-L functional, the present study shows that it can perform well for small and medium-sized solids. This can be ascribed to the fact that the parameters in the M06-L functional are obtained from fits to molecular systems, allowing it to describe mid-range interactions to some extent, which can be seen from its good performance for the ionic crystals considered here.

Furthermore, the current work suggests that investigating the relationships between errors in cohesive properties can help us to understand performance and capability of DFAs. In particular, it is found that the system-dependent behavior of the cohesive properties is best reproduced at the PBE level. Indeed, by using the PBE functional an increase of errors in bulk moduli is accompanied by that in lattice constants for many solids in the current dataset. This is consistent with the finding by Grabowski and co-authors, who discussed the dependence of deviations given by the LDA and PBE functionals for fcc metals [57]. The behavior of the increasing errors with mass, the 'sizable error' [92], can also be observed in the lattice constants. In figure 8, the deviations of the calculated lattice constants from experiment are plotted as a function of their experimental values for 24 metals and 40 non-metals using the six functionals. One can see that the monotonic relationship, indicating that the systematic error increases as the crystal unit-cell volume increases, is better reproduced by PBE and HSE06, rather than LDA, PBEsol, M06-L or SCAN. The trend is more pronounced in non-metals than metals, in good agreement with experiment (see figures 6 and 7). For 40 non-metallic solids, a nearly monotonic behavior is captured using PBE, while two separate regions are found using M06-L. Our analysis suggests the size-related error issue can be related to the lack of the long-range electron correlation in density approximations, which is discussed further below.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The accurate prediction of both structural and energetic properties of bulk solids remains a difficult task for current DFAs. The SCAN meta-GGA and HSE06 hybrid functionals, which belong to the third and the fourth rungs of Jacob's ladder, were specifically developed to improve upon conventional LDA and GGAs. As we have seen, semiconductors (24 solids) and insulators (six solids) benefit greatly from SCAN and HSE06 for predicting structural properties (a₀ and B₀): SCAN yields remarkably accurate lattice constants (MAREs: 0.5% and 0.2% versus 1.6% and 1.8% by PBE), bulk moduli (MAREs: 4.4% and 3.8% versus 16% and 14.5%), while HSE06 performs reasonably well for both properties, with the errors reduced by more than a factor of two compared to the PBE results. The SCAN functional improves on PBE for predicting all the three properties of semiconductor and insulator solids simultaneously and this is clearly a substantial improvement for a semi-local DFA.

However, for transition metals meta and hybrid functionals generally yield no improvement in accuracy over GGA functionals [28, 31, 92]. For the lattice constants, PBE yields fairly accurate results for 3d metals, but systematically overestimates them for 4d and 5d metals, consistently underestimating both bulk moduli and cohesive energies. In particular, PBE fails to describe the cohesion in coinage metals including Pd, Ag, Pt, and Au. The averaged errors are about 2% larger in the lattice constants, 17% in the bulk moduli, and 11% in the cohesive energies for those four metals. Note that the size-related issue discussed in section 3.3, the increase of the error with increasing nuclear mass, is more pronounced here. The origin of this size-related error is not fully understood, but it is likely that this can be related to the presence of self-interaction errors and/or the lack of the long-range electron correlation in semi-local DFAs. As vdW interactions arise from the correlated motion of electrons, they are larger in magnitude for larger systems (i.e., with more electrons), and constitute a significant part of the long-range electron correlation energy. The inclusion of the attractive vdW energy helps shorten bond lengths and stabilize the crystal.

The SCAN functional (as well as M06-L) captures mid-range interactions to a certain extent, which can be seen from its improved description for the interaction energies of the S22 dataset [34]. When applied to solid transition metals, our results suggest that SCAN achieves a major improvement over PBE for 4d and 5d coinage metals where the electron correlation becomes more important between the almost filled d shells. M06-L also performs better for medium-sized unit cells, but for those solids with a large unit-cell volume (heavy nuclear mass), it performs worse than PBE for cohesive properties. This is also an indication that vdW interactions are responsible for part of the deviations of modern DFAs, with the meta-GGA form of SCAN and M06-L allowing them to capture mid-range vdW interactions. This likely accounts for some of the better performance of these functionals for some systems. However, the long-range vdW interactions are still missing from these functionals [36, 35], with the present study showing noticeable deviations from experiment for solid systems including 3d metals, alkali and alkaline-earth metals.

Part of these deviations will also stem from self-interaction errors [93]. The HSE06 hybrid functional, partially counters self-interaction errors by admixing a fraction of exact exchange. As already seen in section 3.2 and other studies [48, 49], HSE06 yields improved performance on the cohesive properties of molecules, semiconductors, and insulators. While for metallic systems HSE06 performs worse than PBE, particularly for the cohesive energies of transition metals, the HSE06 results deviate significantly from experimental values, mainly due to the incorrect treatment of localized d electrons. This illustrates both the impact of self-interaction errors and the challenge of accounting for them. Therefore, given the different contributions that DFAs need to account for (i.e., adding vdW interactions and correcting self-interaction errors), it is not surprising that different semi-local and hybrid functionals underbind bulk solids and yield increased errors with increasing system size (nearest-neighbor distance, lattice constant, or nuclear mass).

As noted previously, one potential source of the size-related issues is the lack of a full description of vdW interactions in the six functionals studied, particularly, the lack of long-range vdW interactions, which form part of the correlation energy. The importance of such interactions for bulk solids has become clearer in recent years [80, 94]. In our previous study of six ionic and semiconductor solids, long-range vdW interactions were included using C₆ coefficients derived from accurate dielectric functions. The inclusion of long-range vdW interactions on top of the PBE and HSE06 functionals lead to a significant improvement for the cohesive properties of the six solids compared to the results without vdW interactions, and as a result the characteristic underbinding of PBE and HSE06 was largely remedied [80].

Here we extend this study to 23 semiconductors, comprising all of the SC group in the present study, apart from tin. The vdW interactions have been modeled using the method proposed by Tkatchenko et al [75], which accounts for both short-range hybridization and long-range screening effects, seamlessly. As briefly discussed in section 2.3, recent studies demonstrated that screening effects become important in large molecules and materials [71, 75, 76]. Figure 9 shows the MAREs in the cohesive properties (a₀, E₀, and B₀) obtained for the 23 semiconductors using the six functionals considered in the present study, alongside the results where three of the functionals are augmented with long-range vdW interactions, specifically, PBE+vdW, M06-L+vdW, and HSE06+vdW (corresponding data are shown in table SVIII). Employing PBE+vdW leads to a consistent improvement in all three cohesive properties compared to PBE alone, reducing deviations from experiment significantly, especially for the lattice constants and bulk moduli. This is unsurprising, as the trends for PBE showed a size-related dependence for the deviations from experiment, with PBE performing worse for systems with heavier elements and, hence, larger lattice constants. As such larger systems will have a larger contribution from vdW interactions, augmenting PBE with a vdW term clearly and understandably leads to improvement.

Zoom In Zoom Out Reset image size

Turning to HSE06, again, all three cohesive properties are in better agreement when the functional is augmented with a vdW term. Overall, HSE06+vdW yields the most accurate results among the DFT+vdW approaches employed, solving the underbinding issue of HSE06, and in fact leads to slight overbinding, with overestimation of both cohesive energies and bulk moduli. In contrast, augmenting M06-L with a vdW term improves lattice constants and bulk moduli but worsens the MARE in the cohesive energies, and as a result gives mixed performance. This likely stems from the fact that the functional form of M06-L contains many empirical parameters that are determined through fits to molecular data. This will partly capture some vdW interactions, but effects and contributions outside or beyond the dataset will not be captured, leading to limited transferability. Similar to M06-L, the ability to capture some mid-range vdW interactions of SCAN, which has been seen from the improved performance on e.g., the S22 set [34] and the bulk solids in this work, makes it difficult to further develop DFT+vdW methods. One of our co-authors [95] has recently illustrated this particular issue of coupling vdW interactions with semi-local functionals, pointing out that the effective range (at which distance fluctuations are correlated) of SCAN depends on system size. Capturing new or higher-level contributions will always be a challenge for empirical functionals (MGGAs or GGAs). Developing new functionals and augmenting existing ones with physically motivated contributions, such as long-range vdW interactions, may often be easier starting from a non-empirical functional that is not fitted to specific types of systems or interactions.