We show that the AM05 functional [Armiento and Mattsson, Phys. Rev. B72, 085108 (2005)] has the same excellent performance for solids as the hybrid density functionals tested in Paier et al [J. Chem. Phys.124, 154709 (2006); 125, 249901 (2006)]. This confirms the original finding that AM05 performs exceptionally well for solids and surfaces. Hartree–Fock hybrid calculations are typically an order of magnitude slower than local or semilocal density functionals such as AM05, which is of a regular semilocal generalized gradient approximation form. The performance of AM05 is on average found to be superior to selecting the best of local density approximation and PBE for each solid. By comparing data from several different electronic-structure codes, we have determined that the numerical errors in this study are equal to or smaller than the corresponding experimental uncertainties.

Density functional theory1,2 (DFT) has become the foundation of most large-scale quantum mechanical simulations. The success stems from the theory’s good quantitative results for a broad range of systems in combination with its relatively low computational cost. At the core of every Kohn–Sham DFT calculation lies the exchange-correlation (XC) functional. It is, in principle, the only limiting factor for the accuracy of the calculations.63 The development of new functionals is, therefore, of utmost importance to the progress of computational materials science, -physics, -chemistry, and -biology. Despite the importance, significant improvements of the exchange and correlation treatment have been few. Since DFT is increasingly being employed for large systems (several hundred atoms) and for long (tens of picoseconds) molecular dynamics (MD) simulations, the trade-off between speed and accuracy is arising as an additional major concern in functional development.

In this article, we assess the performance of the Armiento–Mattsson 2005 functional3 (AM05) for a large set of crystalline solids. We show that AM05 systematically improves upon earlier functionals of the same class [the density and gradient based functionals local density approximation (LDA),2 PBE,4 BLYP,5,6 and RPBE7]. The AM05 functional, in fact, performs as well as the hybrid functionals PBE0 (Refs. 8 and 9) and HSE06,10 see Table I. In addition to two functionals commonly used for solid-state applications, LDA and PBE, we chose to include BLYP and RPBE because of the large and growing interest in water-solid interactions.14 PBE and BLYP are both used extensively for water14–19 and RPBE has been suggested to give a water structure in good agreement with experimental results.15 The performance of AM05 for water will be presented elsewhere,20 but preliminary results show that the behavior is similar to that of PBE. For water-solid interactions to be modeled correctly, however, it is crucial to model the bulk solid well to begin with. Crystallization of ice on a substrate is one example. Formation of ice depends sensitively on the matching of lattices, hence, the lattice constant of the solid has to be described with high accuracy.

Table I.

Mean error (ME), mean absolute error (MAE), and root mean square error (RMSE) for lattice constants a0 (Å) and bulk moduli B0 (GPa), for the functionals tested in this work [AM05, PBE, LDA, RPBE, BLYP, and the best of LDA or PBE (LoP) with respect to lattice constant or bulk modulus] and in Ref. 10 [PBE0, HSE06, and PBE, here denoted as PBE ()]. We have used VASP 5.1 (Refs. 10–13) and used the same PAW core potentials as were employed in Ref. 10. The excellent agreement between the present results and those of and Ref. 10 can be seen by comparing PBE and PBE (). Of the nonhybrid functionals, only AM05 performs as well as the hybrids. However, the computational cost using AM05 is only a fraction of that using hybrids.

 a0 (Å)B0 (GPa)
MEMAERMSEMEMAERMSE
AM05 0.001 0.025 0.033 4.48 8.10 11.2 
PBE0 0.007 0.022 0.029 0.1 7.9 11.3 
HSE06 0.010 0.023 0.030 1.6 8.6 12.8 
LoP (a0) 0.006 0.040 0.048 6.67 11.6 16.6 
LoP (B0) 0.003 0.049 0.053 2.45 7.26 9.79 
PBE () 0.039 0.045 0.054 12.3 12.4 16.4 
PBE 0.039 0.046 0.056 14.1 14.2 18.3 
LDA 0.070 0.070 0.082 7.48 10.7 15.2 
RPBE 0.090 0.091 0.113 17.9 20.5 24.7 
BLYP 0.093 0.100 0.114 26.0 26.1 32.2 
 a0 (Å)B0 (GPa)
MEMAERMSEMEMAERMSE
AM05 0.001 0.025 0.033 4.48 8.10 11.2 
PBE0 0.007 0.022 0.029 0.1 7.9 11.3 
HSE06 0.010 0.023 0.030 1.6 8.6 12.8 
LoP (a0) 0.006 0.040 0.048 6.67 11.6 16.6 
LoP (B0) 0.003 0.049 0.053 2.45 7.26 9.79 
PBE () 0.039 0.045 0.054 12.3 12.4 16.4 
PBE 0.039 0.046 0.056 14.1 14.2 18.3 
LDA 0.070 0.070 0.082 7.48 10.7 15.2 
RPBE 0.090 0.091 0.113 17.9 20.5 24.7 
BLYP 0.093 0.100 0.114 26.0 26.1 32.2 

While the XC functional determines the fundamental accuracy of the calculation, there is a second source of errors, the numerical precision in solving the Kohn–Sham equations. Implementation-related approximations, such as basis sets, pseudopotentials, approximate matrix diagonalization methods, plane-wave cutoff energies, etc. (see, e.g., Ref. 21) can all be successively improved by increasing the computational expense (i.e., by “converging” the calculations). The concept of precision is particularly important in the area of functional development. When the differences between functionals are small,22,23 the precision of the calculations has to be high enough to resolve them.

In the following, we will give a brief theoretical background that covers the different functionals used in our study. Following this, we present a number of calculations of lattice constants and bulk moduli. Finally, the results are analyzed and discussed.

In the Kohn–Sham DFT computational scheme,2 the ground state electron energy is obtained via the solution of the Kohn–Sham (KS) equations. These equations resemble the Hartree–Fock (HF) equations. However, whereas the Hartree–Fock equations have a fully nonlocal exchange potential, the KS equations instead have an XC potential vxc(r) that is diagonal and local in real space. This locality means that the equations are solved with significantly less computational expense. Also, in difference to the Hartree–Fock equations, the transformation of the many-body electron problem into the KS equations introduces no approximations in itself. All correlation effects can formally be accounted for within the XC potential vxc(r). In practice, the crucial question is how good the approximation of this quantity is.

The XC potential vxc(r) is obtained through functional differentiation of the DFT XC energy Exc[n] as a functional of the ground state electron density n(r). The XC energy is itself obtained from an integration of the exchange-correlation energy per particle ϵxc(r;[n]),

Exc[n]=n(r)ϵxc(r;[n])dr.
(1)

A DFT XC functional is usually understood to mean an approximation to ϵxc([n];r). This quantity is further split in separate exchange and correlation parts ϵxc=ϵx+ϵc. This separation originates from the choice that the exchange part should give the energy obtained from the Hartree–Fock exchange expression when the one-particle orbitals ϕi(r) obtained from solving the KS equations are inserted (in hartree atomic units, for a spin-unpolarized system, i.e., spin up and down orbitals are always filled equally),

ϵx(r;[n])=1n(r)rrioccϕi(r)ϕi*(r)2dr.
(2)

However, any transformation, e.g., by integration by parts, of Eq. (1), gives alternative definitions of the exchange energy per particle24 that are equally valid.

The LDA XC functional was presented already in the early works on DFT.2 LDA obtains the exchange energy from that of a uniform electron gas having a density equal to the local density n(r) at each spatial point r. Using this model system, the exchange energy can be derived exactly, giving the following exchange energy per particle,

ϵxLDA(n(r))=34π(3π2n(r))13.
(3)

The correlation energy per particle ϵcLDA is obtained as a parametrization of quantum Monte Carlo data25 for the uniform electron gas at different densities. In this work we use the parameterization of Perdew and Wang26 for the all-electron full-potential calculations performed with the computer code named RSPT and that of Perdew and Zunger27 for the VASP calculations, but these and other parameterizations are largely equivalent.

Despite its simple construction, LDA has proven successful for many applications, in particular, for solids. To improve upon the LDA form, more degrees of freedom than only the local value of the electron density n(r) must be introduced in the approximation for the XC energy per particle. The usual way to extend the LDA form is to introduce the electron density gradient n(r), giving a generalized gradient approximation (GGA) form for the XC energy approximation,

Exc[n]=n(r)ϵxcGGA(n(r),n(r))dr.
(4)

Historically, attempts were made to use a straightforward expansion of the XC energy in the density to produce a systematic improvement of LDA of this form. However, the outcome was disappointing and generally did not improve upon LDA results. Instead, several other approaches have been pursued. The ones important for this work are briefly discussed in the following sections.

The method of employing model systems builds upon the strength of LDA: The exchange and correlation energy expressions stem from a single model system, the uniform electron gas, for which the exact results are known. This leads to an internal consistency between the exchange and correlation approximations, which makes the combined XC quantity more widely applicable than the individual approximations. We refer to this property as “compatible” exchange and correlation functionals. The compatibility manifests itself as a strong cancellation of errors between the exchange and correlation parts of LDA.

Kohn and Mattsson discussed the creation of a XC functional from a surface-oriented model system and its possible combination with another treatment where this model was unsuitable.28,29 The approach was formalized and generalized in the subsystem functional scheme by Armiento and Mattsson.24 The idea is to create separate functionals from different model systems and merge them using a density functional index30 that locally determines the nature of the system. These ideas were made concrete in the AM05 functional.3 It involves the following two model systems: For regions that are locally bulklike, the uniform electron gas is used; for regions that are locally surfacelike, a surface functional is derived from the Airy gas28 in combination with the jellium surfaces.31 

The AM05 surface exchange functional is a parameterization of the Airy electron gas data.28 Such a parameterization was made by Vitos et al.,32 but the AM05 exchange functional improves on it by imposing the correct limiting behavior for large scaled density gradients. The local Airy approximation (LAA) parameterization used in AM05 is given by

ϵxLAA(n(r),n(r))=ϵxLDA(n(r))FxLAA(s),
(5)

where the scaled density gradient s=n(r)[2(3π2)13n43(r)], and the refinement function FxLAA(s) is defined as

FxLAA(s)=(cs2+1)(cs2Fxb+1),
(6)

where c=0.7168 is a fitted constant. The form of FxLAA(s) is chosen to impose a correct uniform limit onto Fxb, which is constructed as an analytical interpolation between two known limits of the Airy refinement function as discussed below.

An effective scaled z coordinate in the Airy gas model,3,28 exact in the limit of high values of s, can be defined as

ζ̃(s)=[32W(s3226)]23,
(7)

where W is the Lambert W function.33 (This function can be calculated to machine accuracy by a few iteration steps implemented in a short piece of code; a routine is available from the authors.)

The properties of ζ̃(s) makes it possible to create an Airy refinement function that simultaneously satisfy both the true high and low s limits,

Fxb=1[ϵxLDA(ñ0(s))2ζ͌(s)],
(8)

where ζ̃̃(s) is a suitable interpolation between the two limits,

ζ͌(s)=([(43)132π3]4ζ̃(s)2+ζ̃(s)4)14,
(9)

and where n0(s) is an effective density defined from the effective z coordinate ζ̃(s),

ñ0(s)=ζ̃(s)323π2s3.
(10)

For a fully compatible setup, the ϵxLAA exchange functional should be combined with a correlation functional obtained as a parameterization of, e.g., quantum Monte Carlo data for the Airy gas system. However, with no such data available, Armiento and Mattsson derived a semicompatible surface correlation from the XC data available for jellium surface models, based on the idea that both the Airy gas and jellium surface models are related to similar surface physics. The basic idea is to correct for the difference in surface correlation, as compared to bulk correlation, by a simple scaling factor to the LDA expression. The scaling factor was obtained from a fit to the combined exchange and correlation energies for surface jellium data,34 giving

ϵc(r;[n])=γϵcLDA(n(r)),γ=0.8098.
(11)

This fit used the Perdew–Wang parametrization of the LDA correlation,26 which should thus always be used within AM05.

The surface functional is combined via the subsystem functional scheme with LDA for bulklike regions using a density index,

X=1αs2(1+αs2),
(12)

where α=2.804 was obtained simultaneously with γ in the fit to surface jellium energy data [cf. comment to Eq. (11)]. The final composed expression for the AM05 exchange and correlation functionals are

(13)
ϵxLAA(r;[n])=ϵxLDA(n(r))[X+(1X)FxLAA(s)],
ϵcLAA(r;[n])=ϵcLDA(n(r))[X+(1X)γ].
Note that, although AM05 is of the same form as a GGA, and can be easily implemented35 in a code using the same input quantities as, e.g., PBE does, it stems from a different theoretical background than functionals originally referred to as GGAs (PW91, BLYP, and PBE).

The AM05 functional is the first functional constructed according to the subsystem functional scheme. It can be seen as a consistent improvement over LDA in the sense that it reproduces the exact XC energy for two types of model systems: The uniform electron gas and the jellium surfaces, describing two situations with fundamentally different physics. The subsystem functional approach outlines a series of further improved functionals, achieved by keeping the current exact XC model systems and adding others. In the application of AM05 to the solid-state systems of the current work (see Table I), we see how this formal improvement over LDA also translates into significantly improved numerical results.

When LDA was analyzed to better understand its success for systems dominated by electron densities far from uniform, it was observed that LDA fulfills a number of exact constraints that one can show also the exact XC functional fulfills. Several works, most prominently the ones of Perdew and co-workers, focus on this observation to argue that improved functionals can be constructed by retaining the constraints LDA fulfills and adding new ones.

One of the most prominent examples of a GGA XC functional is the popular and successful one of Perdew, Burke, and Ernzerhof (PBE).4 For systems such as atoms and molecules, PBE constitutes a significant improvement over LDA. The derivation explicitly focuses on seven constraints that are fulfilled by the construction. Some of these are argued to be the energetically significant constraints of LDA, and others are new ones fulfilled in addition (see, e.g., the discussion in Ref. 36). However, despite PBE being constructed to fulfill the important constraints of LDA, it is not a uniform improvement on LDA for solids. Rather, for combined XC, PBE is, for example, less accurate than LDA for surface jellium,36 a traditional solid-state model system. Although PBE often gives lattice constants in better agreement with experiments than LDA does, the same is not true for the bulk moduli (see Table I).

Since the publication of AM05 in 2005, other authors have also attempted to create functionals with improved performance for solids. Wu and Cohen’s approach37 results in the exchange term having a weaker dependence on the scaled gradient s than the PBE exchange has; a relatively weak s dependence is a property of both LAG32 and AM05.3,32 A very recent functional by Perdew et al.38 also has an exchange term with a weak s dependence. For correlation, they follow closely the construction of AM05, and fit to the combined exchange and correlation energy for surface jellium. The resulting functional gives values close to AM05 over a range of s typically encountered in real solids, and, thus, we predict that this new functional will show a performance similar to that of AM05.

By allowing additional degrees of freedom in the XC energy per particle approximation than the gradient of the density, as used in the GGA form Eq. (4), further constraints can be satisfied. This leads to the concept of meta-GGAs which include, e.g., the kinetic energy density of the KS orbitals. The TPSS39 meta-GGA functional is reported to give improved results over PBE, but since the testing was done on a slightly different set of solids and with a different code, it is not possible to make direct comparisons in this work (see Sec. IV B).

An alternative to the constraint-based functional construction is the more pragmatic approach of empirical functionals. The governing principle is that good exchange and correlation functionals can be obtained from suitable generic expressions fitted to known energies. Traditionally, this approach has been more prevalent for functionals aimed at chemical systems than solid-state systems; perhaps due to the more readily available high-quality total energy data of, e.g., atoms that can be used for such fitting.

One of the most widely used functionals created from this idea is the BLYP XC functional. It is composed by the B88 exchange functional5 and the LYP correlation.6 Both these functionals rely on fits of their free parameters to atomic data. The BLYP functional has proven very successful for various applications in chemistry. However, from a theoretical standpoint, the LYP functional has been strongly criticized for a number of shortcomings: (i) It does not reproduce LDA in the limit of slowly varying densities; (ii) it gives zero correlation energy for any fully spin-polarized system; and (iii) its derivation involves theoretical problems related to non-normalized wavefunctions.40 BLYP was implemented in VASP as a part of a recent study of the properties of B3LYP.41 

There have also been attempts to turn PBE into a semiempirical functional. Zhang and Yang created revPBE by giving up one of the PBE constraints, and instead refit one of its parameters to data on atoms ranging from helium to argon.42 Hammer et al. observed that the refit not only gave improved chemisorption energies,7 but also that this improvement could still be achieved without giving up any of the original PBE constraints by changing the form of the expression for the exchange functional. The result, the RPBE functional,7 is not explicitly empirical in the sense that it does not directly involve any fits. However, the new exchange functional form was chosen because it reproduces relevant behavior of the exchange of revPBE, which was fitted to atomic data. As opposed to the LAG and AM05 functionals, the dependence on the scaled gradient s is stronger for the RPBE functional than for the PBE functional. Because of this, RPBE is not expected to improve on PBE for solid-state systems. On the contrary, because of the relatively strong s dependence the volume should increase compared to the PBE functional. This is indeed what we find in our tests (see Table I). It is, however, still relevant to include RPBE in the comparisons, since the choice of functional for mixed systems involves a trade-off with the goal of performing well enough for all parts. For such applications, our results may help assess how large the trade-off is in using RPBE for the solids.

Hybrid functionals, which are characterized by the admixture of a certain amount of nonlocal Fock exchange to a part of local or semilocal DFT exchange, are extensively applied in the field of quantum chemistry. Results obtained using hybrid functionals are usually in significantly better agreement with experiments than those obtained using local or semilocal DFT exchange-correlation functionals (see, e.g., Ref. 43).

Becke was the first to successfully formulate a true Hartree–Fock/DFT hybrid scheme, known as “half-and-half functional.”44 In a subsequent publication, the semiempirical three-parameter hybrid functionals were introduced.45 Although Becke’s concept is motivated by the adiabatic connection formula for the exchange-correlation energy,46–49 the free parameters are determined by a least-squares fit to experimental atomization energies, electron and proton affinities, and ionization potentials of atomic and molecular species in the G2 test set. The B3LYP hybrid functional has become one of the most popular semiempirical hybrid functionals in computational chemistry.

In Ref. 41, the performance of the B3LYP hybrid functional applied to crystalline solids is scrutinized. There, Paier et al. show that B3LYP performs significantly worse than even simpler and computationally less expensive gradient corrected functionals (such as PBE) for the prediction of almost any relevant property (lattice constants, bulk moduli, or cohesive energies) of systems containing elements beyond the second row. Moreover, cohesive energies for metals are wrong by up to typically 50%.41 The bad performance for metals is not surprising since the B3LYP functional does not describe the homogeneous electron gas exactly.

A concept to rationalize the amount of Fock exchange admixed to the standard DFT exchange energy was devised by Perdew and co-workers.8,50,51 Hybrid functionals motivated by this work are, often termed nonempirical or parameter-free and can be based on GGAs (PBE0, Refs. 9 and 52) or meta-GGAs (TPSSh, Ref. 53). The accuracy of nonempirical hybrid functionals has been shown to be not far from that of semiempirical ones for molecular systems. The way the functionals are derived and the lack of empirical parameters fitted to specific properties make these functionals applicable to both quantum chemistry and condensed matter physics. Comprehensive studies of the performance of PBE08,9 and HSE0655 for crystalline solids were published recently.10,41,54,55

A major limitation of hybrid schemes, however, is a steep increase in computational cost for periodic systems. Exact exchange HF calculations are typically an order of magnitude slower than local or semilocal density functional calculations. The actual difference in computational cost depends on program package and the electronic structure of the system (metallic, semiconducting, or insulating). At present, the applicability of hybrids to large systems and/or long MD simulations is substantially limited by computational resources.

Our VASP results are summarized in Tables II and III. The pseudopotential, plane-wave calculations were done with VASP 5.1,10,11,13 using projector augmented wave (PAW) core potentials.56 We applied the same PAW cores as those used by Paier et al. in their assessment of hybrid functionals for solids.10 The PAW implementation in VASP 5.1 allows use of multiple XC functionals on the same set of core potentials13 while retaining high precision. Although the core wave functions are frozen in the configuration determined as the PAW core is constructed (say, an LDA atom), the core-valence interaction is consistently recalculated with the selected functional. Transferability errors are hence reduced.13 By performing the AM05 and BLYP calculations on both the LDA and PBE cores, we conclude that the errors thereby introduced are insignificant. In addition, we did full-potential all-electron calculations for all solids, the results of which confirm that the PAW cores developed for VASP are of very high fidelity.

Table II.

Lattice constants a0 (Å), ME, MAE, RMSE, and mean absolute relative error (MARE), obtained with the AM05, LDA, PBE, BLYP, and RPBE functionals, using VASP. The experimental (Exp) results are the same as used in Refs. 10 and 61.

SolidExpAM05LDAPBEBLYPRPBE
Li 3.477 3.455 3.359 3.433 3.421 3.476 
Na 4.225 4.212 4.052 4.201 4.210 4.295 
Al 4.032 4.004 3.984 4.041 4.116 4.064 
BN 3.616 3.605 3.583 3.627 3.647 3.646 
BP 4.538 4.516 4.491 4.548 4.592 4.573 
3.567 3.551 3.534 3.573 3.598 3.590 
Si 5.430 5.431 5.403 5.467 5.532 5.499 
SiC 4.358 4.350 4.330 4.377 4.411 4.398 
β-GaN 4.520 4.492 4.460 4.548 4.611 4.511 
GaP 5.451 5.441 5.394 5.509 5.607 5.556 
GaAs 5.648 5.672 5.611 5.755 5.871 5.812 
LiF 4.010 4.039 3.908 4.065 4.084 4.146 
LiCl 5.106 5.119 4.962 5.150 5.232 5.254 
NaF 4.609 4.686 4.508 4.708 4.716 4.824 
NaCl 5.595 5.686 5.466 5.702 5.763 5.847 
MgO 4.207 4.232 4.168 4.259 4.281 4.302 
Cu 3.603 3.565 3.523 3.637 3.711 3.682 
Rh 3.798 3.773 3.757 3.833 3.905 3.857 
Pd 3.881 3.872 3.844 3.946 4.034 3.984 
Ag 4.069 4.054 4.002 4.150 4.262 4.215 
ME ⋯ 0.001 0.070 0.039 0.093 0.090 
MAE ⋯ 0.025 0.070 0.046 0.100 0.091 
RMSE ⋯ 0.033 0.082 0.056 0.114 0.113 
MARE ⋯ 0.6% 1.6% 1.0% 2.2% 2.0% 
SolidExpAM05LDAPBEBLYPRPBE
Li 3.477 3.455 3.359 3.433 3.421 3.476 
Na 4.225 4.212 4.052 4.201 4.210 4.295 
Al 4.032 4.004 3.984 4.041 4.116 4.064 
BN 3.616 3.605 3.583 3.627 3.647 3.646 
BP 4.538 4.516 4.491 4.548 4.592 4.573 
3.567 3.551 3.534 3.573 3.598 3.590 
Si 5.430 5.431 5.403 5.467 5.532 5.499 
SiC 4.358 4.350 4.330 4.377 4.411 4.398 
β-GaN 4.520 4.492 4.460 4.548 4.611 4.511 
GaP 5.451 5.441 5.394 5.509 5.607 5.556 
GaAs 5.648 5.672 5.611 5.755 5.871 5.812 
LiF 4.010 4.039 3.908 4.065 4.084 4.146 
LiCl 5.106 5.119 4.962 5.150 5.232 5.254 
NaF 4.609 4.686 4.508 4.708 4.716 4.824 
NaCl 5.595 5.686 5.466 5.702 5.763 5.847 
MgO 4.207 4.232 4.168 4.259 4.281 4.302 
Cu 3.603 3.565 3.523 3.637 3.711 3.682 
Rh 3.798 3.773 3.757 3.833 3.905 3.857 
Pd 3.881 3.872 3.844 3.946 4.034 3.984 
Ag 4.069 4.054 4.002 4.150 4.262 4.215 
ME ⋯ 0.001 0.070 0.039 0.093 0.090 
MAE ⋯ 0.025 0.070 0.046 0.100 0.091 
RMSE ⋯ 0.033 0.082 0.056 0.114 0.113 
MARE ⋯ 0.6% 1.6% 1.0% 2.2% 2.0% 
Table III.

Bulk moduli B0 (GPa), obtained with the AM05, LDA, PBE, BLYP, and RPBE functionals, using VASP. The experimental (Exp) results are the same as used in Refs. 10 and 61. For the materials marked with a star “,” different experimental data are frequently quoted in the III–V literature, see Sec. IV C.

SolidExpAM05LDAPBEBLYPRPBE
Li 13.0 13.0 15.1 13.7 13.7 13.1 
Na 7.5 7.36 9.22 7.62 7.08 6.94 
Al 79.4 83.9 81.4 75.2 54.9 73.7 
BN 400* 378 394 365 350 353 
BP 165 165 171 158 146 152 
443 442 456 424 399 410 
Si 99.2 90.2 93.6 86.4 77.0 83.1 
SiC 225 217 224 208 194 201 
β-GaN 210* 181 196 166 152 237 
GaP 88.7 80.2 87.0 74.3 64.3 69.4 
GaAs 75.6 65.1 71.8 59.7 50.6 55.2 
LiF 69.8 65.8 85.7 66.7 65.5 59.3 
LiCl 35.4 30.3 40.4 31.2 28.9 27.3 
NaF 51.4 43.2 60.1 44.5 44.3 38.3 
NaCl 26.6 22.0 31.4 23.4 22.0 20.1 
MgO 165 152 169 148 145 139 
Cu 142 157 180 134 112 120 
Rh 269 285 304 249 214 232 
Pd 195 194 216 165 137 148 
Ag 109 109 132 88.9 71.2 74.4 
ME ⋯ 4.48 7.48 14.1 26.0 17.9 
MAE ⋯ 8.10 10.7 14.2 26.1 20.5 
RMSE ⋯ 11.2 15.2 18.3 32.2 24.7 
MARE ⋯ 7.1% 10.8% 10.4% 18.7% 16.0% 
SolidExpAM05LDAPBEBLYPRPBE
Li 13.0 13.0 15.1 13.7 13.7 13.1 
Na 7.5 7.36 9.22 7.62 7.08 6.94 
Al 79.4 83.9 81.4 75.2 54.9 73.7 
BN 400* 378 394 365 350 353 
BP 165 165 171 158 146 152 
443 442 456 424 399 410 
Si 99.2 90.2 93.6 86.4 77.0 83.1 
SiC 225 217 224 208 194 201 
β-GaN 210* 181 196 166 152 237 
GaP 88.7 80.2 87.0 74.3 64.3 69.4 
GaAs 75.6 65.1 71.8 59.7 50.6 55.2 
LiF 69.8 65.8 85.7 66.7 65.5 59.3 
LiCl 35.4 30.3 40.4 31.2 28.9 27.3 
NaF 51.4 43.2 60.1 44.5 44.3 38.3 
NaCl 26.6 22.0 31.4 23.4 22.0 20.1 
MgO 165 152 169 148 145 139 
Cu 142 157 180 134 112 120 
Rh 269 285 304 249 214 232 
Pd 195 194 216 165 137 148 
Ag 109 109 132 88.9 71.2 74.4 
ME ⋯ 4.48 7.48 14.1 26.0 17.9 
MAE ⋯ 8.10 10.7 14.2 26.1 20.5 
RMSE ⋯ 11.2 15.2 18.3 32.2 24.7 
MARE ⋯ 7.1% 10.8% 10.4% 18.7% 16.0% 

The k-point integration was performed using the tetrahedron method with Blöchl corrections57 on an 18×18×18 uniform grid58 centered at the gamma point. All real-space cells are cubic with two atoms for bcc, four atoms for fcc, and eight atoms for the zinc blende and diamond cells. We applied energy cutoffs ranging between 600 and 1000eV: 600eV (Al, Ag, Pd, Rh, Cu, GaAs, GaP), 800eV (Na, NaF, NaCl, MgO, SiC, Si, C, GaN, BN, BP), and 1000eV (Li, LiF, LiCl) with a convergence criterion of 1.0×105eV for the self-consistent loop. The precision in the calculations is thus overall very high. The PBE results can be directly compared to those of Paier et al.;10 the discrepancies are minute and caused by slight differences in the volume range governing the Murnaghan fits. None of the differences affect any conclusions.

To confirm the validity of using PAW core potentials created with LDA or PBE together with AM05 in VASP 5.1, we have performed all-electron calculations for the same set of solids using the RSPT (Ref. 59) code. To further compare the results given by the two codes, we have also performed LDA and PBE calculations with RSPT. The results are compared in Tables IV and V.

Table IV.

Lattice constants a0 (Å), obtained with AM05, LDA, and PBE, using VASP and RSPT. The AM05 VASP results are calculated using LDA and PBE PAW core potentials. As shown, the results are nearly identical. The AM05 values given in Tables II and III are the mean of the AM05 values obtained with LDA and PBE PAW core potentials. The two different codes also give similar results, showing that using PBE or LDA PAW core potentials for AM05 calculations in VASP is a valid approach.

 AM05LDAPBE
  VASPRSPTVASPRSPTVASPRSPT
SolidLDA PAWPBE PAW LDA PAW PBE PAW 
Li 3.4539 3.4559 3.456 3.359 3.362 3.433 3.434 
Na 4.2124 4.2125 4.222 4.052 4.053 4.201 4.196 
Al 4.0003 4.0076 4.008 3.984 3.986 4.041 4.043 
BN 3.6026 3.6071 3.604 3.583 3.583 3.627 3.625 
BP 4.5118 4.5203 4.520 4.491 4.495 4.548 4.553 
3.5497 3.5529 3.551 3.534 3.534 3.573 3.573 
Si 5.4306 5.4317 5.436 5.403 5.405 5.467 5.474 
SiC 4.3491 4.3514 4.361 4.330 4.337 4.377 4.386 
β-GaN 4.4914 4.4921 4.506 4.460 4.465 4.548 4.553 
GaP 5.4385 5.4435 5.457 5.394 5.405 5.509 5.518 
GaAs 5.6689 5.6747 5.686 5.611 5.620 5.755 5.761 
LiF 4.0364 4.0420 4.041 3.908 3.912 4.065 4.065 
LiCl 5.1163 5.1223 5.114 4.962 4.966 5.150 5.149 
NaF 4.6860 4.6866 4.685 4.508 4.507 4.708 4.692 
NaCl 5.6844 5.6877 5.693 5.466 5.467 5.702 5.692 
MgO 4.2352 4.2291 4.221 4.168 4.164 4.259 4.253 
Cu 3.5641 3.5668 3.564 3.523 3.522 3.637 3.633 
Rh 3.7729 3.7729 3.786 3.757 3.769 3.833 3.845 
Pd 3.8713 3.8727 3.880 3.844 3.852 3.946 3.953 
Ag 4.0538 4.0549 4.062 4.002 4.010 4.150 4.155 
ME 0.001 0.002 0.006 0.070 0.066 0.039 0.041 
MAE 0.026 0.025 0.022 0.070 0.066 0.046 0.048 
RMSE 0.033 0.033 0.033 0.082 0.079 0.056 0.056 
MARE 0.6% 0.6% 0.5% 1.6% 1.5% 1.0% 1.1% 
 AM05LDAPBE
  VASPRSPTVASPRSPTVASPRSPT
SolidLDA PAWPBE PAW LDA PAW PBE PAW 
Li 3.4539 3.4559 3.456 3.359 3.362 3.433 3.434 
Na 4.2124 4.2125 4.222 4.052 4.053 4.201 4.196 
Al 4.0003 4.0076 4.008 3.984 3.986 4.041 4.043 
BN 3.6026 3.6071 3.604 3.583 3.583 3.627 3.625 
BP 4.5118 4.5203 4.520 4.491 4.495 4.548 4.553 
3.5497 3.5529 3.551 3.534 3.534 3.573 3.573 
Si 5.4306 5.4317 5.436 5.403 5.405 5.467 5.474 
SiC 4.3491 4.3514 4.361 4.330 4.337 4.377 4.386 
β-GaN 4.4914 4.4921 4.506 4.460 4.465 4.548 4.553 
GaP 5.4385 5.4435 5.457 5.394 5.405 5.509 5.518 
GaAs 5.6689 5.6747 5.686 5.611 5.620 5.755 5.761 
LiF 4.0364 4.0420 4.041 3.908 3.912 4.065 4.065 
LiCl 5.1163 5.1223 5.114 4.962 4.966 5.150 5.149 
NaF 4.6860 4.6866 4.685 4.508 4.507 4.708 4.692 
NaCl 5.6844 5.6877 5.693 5.466 5.467 5.702 5.692 
MgO 4.2352 4.2291 4.221 4.168 4.164 4.259 4.253 
Cu 3.5641 3.5668 3.564 3.523 3.522 3.637 3.633 
Rh 3.7729 3.7729 3.786 3.757 3.769 3.833 3.845 
Pd 3.8713 3.8727 3.880 3.844 3.852 3.946 3.953 
Ag 4.0538 4.0549 4.062 4.002 4.010 4.150 4.155 
ME 0.001 0.002 0.006 0.070 0.066 0.039 0.041 
MAE 0.026 0.025 0.022 0.070 0.066 0.046 0.048 
RMSE 0.033 0.033 0.033 0.082 0.079 0.056 0.056 
MARE 0.6% 0.6% 0.5% 1.6% 1.5% 1.0% 1.1% 
Table V.

Bulk moduli B0 (GPa), obtained with AM05, LDA, and PBE, using VASP and RSPT. The AM05 VASP results are calculated using LDA and PBE PAW core potentials. As shown, the results are nearly identical. The two different codes also give similar results, showing that using PBE or LDA PAW core potentials for AM05 calculations in VASP is a valid approach.

SolidAM05LDAPBE
VASPRSPTVASPRSPTVASPRSPT
LDA PAWPBE PAW LDA PAW PBE PAW 
Li 13.01 12.99 13.2 15.1 15.0 13.7 13.9 
Na 7.363 7.361 7.65 9.22 9.16 7.62 7.74 
Al 84.08 83.63 86.2 81.4 82.5 75.2 77.1 
BN 378.5 377.5 384 394 400 365 370 
BP 165.1 164.3 168 171 174 158 160 
442.5 441.4 450 456 465 424 431 
Si 90.30 90.11 92.0 93.6 95.4 86.4 87.5 
SiC 216.9 216.3 217 224 226 208 208 
β-GaN 180.6 180.5 183 196 199 166 170 
GaP 80.31 80.13 81.1 87.0 88.2 74.3 75.1 
GaAs 65.08 65.07 65.4 71.8 72.4 59.7 59.4 
LiF 65.85 65.82 65.8 85.7 86.2 66.7 67.5 
LiCl 30.31 30.25 30.7 40.4 41.0 31.2 31.9 
NaF 43.27 43.08 42.4 60.1 60.4 44.5 45.6 
NaCl 22.04 21.99 21.0 31.4 31.5 23.4 23.7 
MgO 151.3 151.9 154 169 171 148 149 
Cu 157.4 157.3 165 180 187 134 140 
Rh 285.3 285.5 293 304 312 249 253 
Pd 194.2 193.9 202 216 224 165 167 
Ag 108.6 108.9 114 132 137 88.9 90.2 
ME 4.38 4.59 1.80 7.48 10.4 14.1 12.1 
MAE 8.03 8.19 9.27 10.7 12.3 14.2 12.2 
RMSE 11.2 11.3 11.9 15.2 18.4 18.3 16.2 
MARE 7.1% 7.2% 8.1% 10.8% 11.8% 10.4% 9.3% 
SolidAM05LDAPBE
VASPRSPTVASPRSPTVASPRSPT
LDA PAWPBE PAW LDA PAW PBE PAW 
Li 13.01 12.99 13.2 15.1 15.0 13.7 13.9 
Na 7.363 7.361 7.65 9.22 9.16 7.62 7.74 
Al 84.08 83.63 86.2 81.4 82.5 75.2 77.1 
BN 378.5 377.5 384 394 400 365 370 
BP 165.1 164.3 168 171 174 158 160 
442.5 441.4 450 456 465 424 431 
Si 90.30 90.11 92.0 93.6 95.4 86.4 87.5 
SiC 216.9 216.3 217 224 226 208 208 
β-GaN 180.6 180.5 183 196 199 166 170 
GaP 80.31 80.13 81.1 87.0 88.2 74.3 75.1 
GaAs 65.08 65.07 65.4 71.8 72.4 59.7 59.4 
LiF 65.85 65.82 65.8 85.7 86.2 66.7 67.5 
LiCl 30.31 30.25 30.7 40.4 41.0 31.2 31.9 
NaF 43.27 43.08 42.4 60.1 60.4 44.5 45.6 
NaCl 22.04 21.99 21.0 31.4 31.5 23.4 23.7 
MgO 151.3 151.9 154 169 171 148 149 
Cu 157.4 157.3 165 180 187 134 140 
Rh 285.3 285.5 293 304 312 249 253 
Pd 194.2 193.9 202 216 224 165 167 
Ag 108.6 108.9 114 132 137 88.9 90.2 
ME 4.38 4.59 1.80 7.48 10.4 14.1 12.1 
MAE 8.03 8.19 9.27 10.7 12.3 14.2 12.2 
RMSE 11.2 11.3 11.9 15.2 18.4 18.3 16.2 
MARE 7.1% 7.2% 8.1% 10.8% 11.8% 10.4% 9.3% 

The RSPT code59 is a full-potential linear muffin tin orbital (LMTO) code. Since it uses an efficient smaller basis set, it is fast compared to other all-electron codes, while the flexible basis, that has been specially built for every single solid, permits very well converged results. Our primary concern has been to construct a basis with small “leakage,” that is, to only keep as core orbitals the orbitals that do not significantly contribute to the density outside of the muffin tin spheres. We have thus, in most cases, used a larger number of valence orbitals than is needed for production runs. We have used the same basis for a specific solid for all different functionals. For each of the seven volume points considered in the fits to the Murnaghan60 equation of state, the muffin tin radius was varied so as to give a muffin tin sphere with a specific fraction of the cell volume. The volumes have been centered around the equilibrium volume given by the VASP calculations and spaced between ±3% of the equilibrium lattice constants. The cells used have all been primitive. We used a 24×24×24k-point grid shifted by (12,12,12) for all solids and we used the tetrahedron method for integrating the k space. The Fourier grid was 30×30×30.

Tables IV and V show that the differences using LDA or PBE core potentials in the VASP AM05 calculations are minute and that either set of results compares well with the all-electron RSPT calculations. Thus, if it is not possible to do both sets of AM05 calculations and use the mean, as we have done in this work, the choice between using LDA or PBE core potentials together with AM05, could be left as a practical consideration. For example, if LDA will also be used, LDA core potentials could be employed also for AM05.64 

The comparison between RSPT and VASP LDA and PBE values also confirms that the VASP LDA and PBE core potentials almost always reproduce the all-electron lattice constants to within 0.1%–0.2%.

With RSPt, we have also confirmed that PBE and PW91 indeed give nearly identical results for lattice constants and bulk moduli of solids, something that led users to believe that PBE and PW91 would perform equally on all types of systems. However, this has been shown not to be the case (see Refs. 22 and 23).

In addition to the set of 20 solids, we next turn to a few additional metals: A heavy bcc metal (W), two heavy fcc elements (Pt and Au), and a light hcp metal (Be). Be is an hcp metal with experimental lattice constants a=2.29Å and c=1.567Å. Our results for (a,c) are the following: (2.232,1.58) for LDA, (2.264,1.58) for PBE, and (2.255,1.58) for AM05. The results for W, Pt, and Au are presented in Table VI. The PAW core potentials include scalar relativistic corrections and treat the semicore p states as valence. The number of valence electrons are thus 12 for W, 11 for Au, and 10 for Pt. It is well known that LDA is superior to PBE for these heavy elements (BLYP and RPBE both give results worse than PBE). AM05 yields lattice constants as well as bulk moduli in very good agreement with experimental data. The addition in AM05 of a surface model system to the LDA bulk region thus did not negatively affect the performance for this class of elements.

Table VI.

VASP results for the heavy metals W, Pt, and Au: AM05 performs as well, or better, than LDA.

 a0 (Å)B0 (GPa)
ExpAM05LDAPBEExpAM05LDAPBE
3.16 3.153 3.142 3.190 310 333 335 310 
Pt 3.92 3.915 3.906 3.977 283 292 307 251 
Au 4.065 4.076 4.053 4.161 166–171 164 183 137 
 a0 (Å)B0 (GPa)
ExpAM05LDAPBEExpAM05LDAPBE
3.16 3.153 3.142 3.190 310 333 335 310 
Pt 3.92 3.915 3.906 3.977 283 292 307 251 
Au 4.065 4.076 4.053 4.161 166–171 164 183 137 

For most realistic solid-state applications, such as size-converged calculations of defect formation energies and migration energies, the use of a hybrid functional is not practical because of the prohibitive computational cost. This is true, in particular, for investigating mixed systems (alloys or solid/molecular systems) at nonzero temperatures. Such calculations require both DFT based molecular dynamics (DFT-MD) and large supercells. For these demanding applications, a functional based only on quantities that are easily calculated, such as density and density derivatives, is currently paramount.

It could be argued that there is no need for new functionals for solid-state systems since the existing LDA and PBE do yield good results. In particular, by selecting the best of LDA or PBE, for a specific system, the result can be further improved. However, such a choice between LDA and PBE is not possible if reliable experimental data is not available. Resorting to using either LDA or PBE depending on application is, therefore, not an approach with predictive power. The approach is particularly uncertain in mixed systems where LDA would be preferred for one component of the system and PBE for another. However, even when disregarding these shortcomings, Table VII shows that relying only on AM05, in general, is a better alternative than choosing the result of LDA or PBE that is closest to experiment. AM05 is significantly better for lattice constants, and, within the relevant precision of bulk moduli calculations, AM05 matches the results obtained when picking LDA or PBE depending on the system.

Table VII.

VASP results for lattice constants a0 (Å) and bulk moduli B0 (GPa), obtained with the best of LDA or PBE. The experimental (Exp) results are the same as used in Refs. 10 and 61. This choice between LDA and PBE is a common practice but the result is less accurate than using AM05.

Solida0 (Å)B0 (GPa)
ExpBest a0a0 for best B0ExpB0 for best a0Best B0
Li 3.477 3.433 3.433 13.0 13.7 13.7 
Na 4.225 4.201 4.201 7.5 7.62 7.62 
Al 4.032 4.041 3.984 79.4 75.2 81.4 
BN 3.616 3.627 3.583 400 365 394 
BP 4.538 4.548 4.491 165 158 171 
3.567 3.573 3.534 443 424 456 
Si 5.430 5.403 5.403 99.2 93.6 93.6 
SiC 4.358 4.377 4.330 225 208 224 
β-GaN 4.520 4.548 4.460 210 166 196 
GaP 5.451 5.394 5.394 88.7 87.0 87.0 
GaAs 5.648 5.611 5.611 75.6 71.8 71.8 
LiF 4.010 4.065 4.065 69.8 66.7 66.7 
LiCl 5.106 5.150 5.150 35.4 31.2 31.2 
NaF 4.609 4.708 4.708 51.4 44.5 44.5 
NaCl 5.595 5.702 5.702 26.6 23.4 23.4 
MgO 4.207 4.168 4.168 165 169 169 
Cu 3.603 3.637 3.637 142 134 134 
Rh 3.798 3.833 3.833 269 249 249 
Pd 3.881 3.844 3.844 195 216 216 
Ag 4.069 4.002 4.150 109 132 88.9 
ME ⋯ 0.006 0.003 ⋯ 6.67 2.45 
MAE ⋯ 0.040 0.049 ⋯ 11.6 7.26 
RMSE ⋯ 0.048 0.053 ⋯ 16.6 9.79 
MARE ⋯ 0.9% 1.1% ⋯ 8.0% 6.2% 
Solida0 (Å)B0 (GPa)
ExpBest a0a0 for best B0ExpB0 for best a0Best B0
Li 3.477 3.433 3.433 13.0 13.7 13.7 
Na 4.225 4.201 4.201 7.5 7.62 7.62 
Al 4.032 4.041 3.984 79.4 75.2 81.4 
BN 3.616 3.627 3.583 400 365 394 
BP 4.538 4.548 4.491 165 158 171 
3.567 3.573 3.534 443 424 456 
Si 5.430 5.403 5.403 99.2 93.6 93.6 
SiC 4.358 4.377 4.330 225 208 224 
β-GaN 4.520 4.548 4.460 210 166 196 
GaP 5.451 5.394 5.394 88.7 87.0 87.0 
GaAs 5.648 5.611 5.611 75.6 71.8 71.8 
LiF 4.010 4.065 4.065 69.8 66.7 66.7 
LiCl 5.106 5.150 5.150 35.4 31.2 31.2 
NaF 4.609 4.708 4.708 51.4 44.5 44.5 
NaCl 5.595 5.702 5.702 26.6 23.4 23.4 
MgO 4.207 4.168 4.168 165 169 169 
Cu 3.603 3.637 3.637 142 134 134 
Rh 3.798 3.833 3.833 269 249 249 
Pd 3.881 3.844 3.844 195 216 216 
Ag 4.069 4.002 4.150 109 132 88.9 
ME ⋯ 0.006 0.003 ⋯ 6.67 2.45 
MAE ⋯ 0.040 0.049 ⋯ 11.6 7.26 
RMSE ⋯ 0.048 0.053 ⋯ 16.6 9.79 
MARE ⋯ 0.9% 1.1% ⋯ 8.0% 6.2% 

Solid-state codes solve the KS equations using approximations that lead to less than perfect precision. When dealing with and comparing such small errors as those of AM05 and hybrids (see Table I), it is relevant to ask whether it is possible to determine which one is more accurate, and whether it is at all possible to conclude that another functional performs better or worse than AM05 (and hybrids) for the present solid-state systems.

Table V of Ref. 10 gives PBE lattice constants and bulk moduli for a subset of the 20 solids studied in this work, calculated using VASP, a full-potential LAPW code, and a Gaussian-type-orbital (GTO) code.

Our RSPt values compare extremely well with the PBE APW+lo values and the VASP values in Table V of Ref. 10, although it seems that the error bars of the RSPt calculations are slightly larger. Specifically, the RSPT lattice constants deviate by up to 0.3% from the VASP values and from the APW+lo values, whereas the difference between the VASP and APW+lo lattice constants is not exceeding 0.1%. It is still clear that all three codes are able to give highly converged results, and the results mutually support each other. In particular, for LDA and PBE, the VASP and RSPT statistical errors are practically identical. For the AM05 case, the use of LDA and PBE cores might slightly increase the VASP error bars, as indicated by the comparison of AM05 for RSPT and VASP. However, the error bars are still smaller than approximately 0.005Å in the lattice constants and 3GPa in the bulk moduli. The same error bars are likely to apply to the HSE03 case, covered originally in Ref. 10. Hence, the difference in accuracy between the hybrids and AM05 cannot be resolved.

The small numerical error bars for the VASP 5.1 PAW calculations in this work and the work of Paier et al.10 are not common to pseudopotential codes.22 Comparing the data in Table I to data obtained with different codes and different types of pseudopotentials can thus not be done with the same small error bars, but the error bars should be expected to be significantly larger, if conventional pseudopotentials are applied.

The same conclusion holds for results obtained with all-electron GTO codes. The results in Table 5 of Ref. 10 (collected from Ref. 61) obtained with a GTO code do not compare as well with the LAPW results as the RSPT ones do.

The precision of the best codes and accuracy of the best functionals for solid-state calculations are now at a level where further improvement in accuracy likely cannot be distinguished without a better understanding of experimental error bars. Although not conclusive in itself, we note that the three best functionals in this study (AM05, HSE06, and PBE0) give the same accuracy within numerical uncertainties, indicating that the remaining disagreement might be a sign of experimental errors and not a consequence of functional accuracy.

An analysis of the VASP results in Table III reveals that two solids alone account for a large part of the mean bulk moduli errors between AM05 and experimental data. The solids, with their respective errors in parenthesis, are GaN (29.5GPa) and BN (22.0GPa). Both errors are surprisingly large, at first suggesting that AM05 has a systematic problem with III-V nitrides. However, theoretical studies more specifically targeting these kinds of system quote other experimental values for the bulk modulus, namely, 190GPa for GaN and 369GPa for BN.62 Using these values instead of the ones provided in Ref. 61 changes the AM05 errors for these two solids to the more expected 9GPa for GaN and +9GPa for BN. These two changes in experimental values yield a mean error (ME) of 1.93GPa, a mean absolute error (MAE) of 6.45GPa, a root mean square error (RMSE) of 8.19GPa, and a mean absolute relative error (MARE) of 6.5% for AM05. Adapting this change in experimental values, the PBE0 and HSE06 results for these quantities would be slightly increased. This observation does not imply that AM05 has better accuracy for the bulk moduli than do the hybrids, but it does illustrate the difficulty in distinguishing between the performance of different functionals which have an accuracy at the level of AM05 and the hybrids. Although a comprehensive review of available experimental results is of interest, not only for GaN and BN, but for all solids, such an investigation is outside the scope of the present work.

A similar examination of the two solids with unusually large AM05 errors in lattice constants, NaF and NaCl, does not provide such a simple explanation. We should note that the AM05 lattice constants for these materials compare better to the experimental values than LDA and PBE, and that both hybrids’ results, though closer, still have substantial deviations. Since RSPT confirms the VASP values, the PAW core potentials cannot be at fault. Considering these caveats, one can draw only few definite conclusions on the relative performance of the HSE06 and AM05 functionals. Generally, AM05 predicts slightly smaller lattice constants than the two hybrid functionals do. Interestingly, it also predicts smaller bulk moduli than HSE06, although the energy volume curvature is usually expected to increase at smaller volumes. Concomitantly, the average bulk moduli are slightly underestimated by AM05 compared to experiment. The underestimation is modest for the metals and semiconductors, but clearly increases towards more ionic compounds and becomes as large as 20% for the two most ionic compounds, NaF and NaCl. As already mentioned, these are the compounds with the largest errors in the AM05 lattice constants. The hybrid functionals performed very well for these two systems yielding a much larger—albeit still too small—curvature at the equilibrium volume. We note that these systems, as well as other ionic compounds, are largely exchange dominated in the sense of the adiabatic coupling theorem, i.e., the adiabatic coupling theorem and the GW method suggest that one should include a large fraction of the nonlocal exchange to account for the physics in these systems. In summary, ionic compounds are the only bulk systems where hybrid functionals might offer an advantage over the AM05 functional. The opposite applies to the metallic systems, where the AM05 functional seems to give an overall slightly better description than the hybrid functionals, in particular, for the bulk moduli. More reliable error bars on the experimental values are required before further definite statements can be made.

The subsystem functional AM05 is based on two exact reference systems: The uniform electron gas and the surface jellium. AM05 hence constitutes a systematic improvement upon LDA by adding terms depending on the gradient of the density while maintaining the exact exchange-correlation limit of LDA. The systematic improvement is confirmed numerically by our careful solid-state benchmarks (cf. Table I). Not only is AM05 found to be a significant improvement over LDA, but also for the studied quantities over other often-used functionals, e.g., PBE, BLYP, and RPBE. Within the high numerical precision of this study, AM05 is found to be as accurate as the most advanced hybrid functionals proposed to date, PBE0 and HSE06. The only exception are ionic systems, where AM05 clearly improves upon PBE, but still yields much too large volumes and much too small bulk moduli.

In further analysis, we find the performance of AM05, on average, to be even better than what a computational user would reach by choosing between LDA and PBE with guidance from experimental knowledge. Hence, AM05 provides a predictive power neither of these two functionals possess. Our comparison of results from different codes also shows that the level of precision available in modern codes, and the size of experimental errors for solid-state systems, will make it difficult to register a further improved functional for these applications without going far beyond the benchmark systems and properties studied in this work. Note that the AM05 construction was made solely on a theoretical and nonempirical basis. The excellent results thus strongly points to it having a sound theoretical basis. Furthermore, we would like to emphasize that the theoretical foundation of AM05 is fundamentally different from those of other available functionals (cf. Sec. II B). In studies where the calculated value is to be relied upon, in particular, when experimental data are either unavailable or have large uncertainties, more than one functional is often applied as a means to assess the accuracy. It is in such cases beneficial to employ functionals that are based on different fundamental principles.

The high speed and precision of the computer code VASP 5.1 makes it an excellent choice for functional assessment, expanding the possibilities of initial testing from simple properties for few atoms to include also large systems and DFT-MD. For the applications of this work, we found its PAW core potentials to be general enough to be interchangeable between different functionals, which saves the otherwise tedious work of generating pseudopotentials. Note, however, that this is a property pertaining to the PAW scheme and the implementation used by this code. Pseudopotentials generated according to other schemes are generally not interchangeable and doing so may severely affect the computational results.

We thank Peter Feibelman for useful conversations. We have had valuable discussions regarding bulk moduli and equations of state with Michael Desjarlais and John Carpenter. Per Andersson provided advice in how to work with the RSPT code. R.A. gratefully acknowledges support from the Alexander von Humboldt Foundation. A.E.M. and T.R.M. acknowledge support from the LDRD office at Sandia. Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the United States Department of Energy’s National Nuclear Security Administration under Contract No. DE-AC04-94AL85000. Los Alamos National Laboratory is operated by Los Alamos National Security, LLC, for the National Nuclear Security Administration of the U.S. Department of Energy under Contract No. DE-AC52-06NA25396. Support by the Austrian Science Fund (FWF) is gratefully acknowledged.

1.
P.
Hohenberg
and
W.
Kohn
,
Phys. Rev.
136
,
B864
(
1964
).
2.
W.
Kohn
and
L. J.
Sham
,
Phys. Rev.
140
,
A1133
(
1964
).
3.
R.
Armiento
and
A. E.
Mattsson
,
Phys. Rev. B
72
,
085108
(
2005
).
4.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
77
,
3865
(
1996
).
5.
A. D.
Becke
,
Phys. Rev. A
38
,
3098
(
1988
).
6.
C.
Lee
,
W.
Yang
, and
R. G.
Parr
,
Phys. Rev. B
37
,
785
(
1988
).
7.
B.
Hammer
,
L. B.
Hansen
, and
J. K.
Nørskov
,
Phys. Rev. B
59
,
7413
(
1999
).
8.
J. P.
Perdew
,
M.
Ernzerhof
, and
K.
Burke
,
J. Chem. Phys.
105
,
9982
(
1996
).
9.
C.
Adamo
and
V.
Barone
,
J. Chem. Phys.
110
,
6158
(
1999
).
10.
J.
Paier
,
M.
Marsman
,
K.
Hummer
,
G.
Kresse
,
I. C.
Gerber
, and
J. G.
Ángyán
,
J. Chem. Phys.
124
,
154709
(
2006
);
[PubMed]
J.
Paier
,
M.
Marsman
,
K.
Hummer
,
G.
Kresse
,
I. C.
Gerber
, and
J. G.
Ángyán
,
J. Chem. Phys.
125
,
249901
(
2006
).
11.
G.
Kresse
and
J.
Hafner
,
Phys. Rev. B
47
,
R558
(
1993
);
G.
Kresse
and
J.
Hafner
,
Phys. Rev. B
49
,
14251
(
1994
);
G.
Kresse
and
J.
Furthmüller
,
Phys. Rev. B
54
,
11169
(
1996
).
12.
The calculations were made using version 5.1.38 of the VASP code.
13.
J.
Paier
,
R.
Hirschl
,
M.
Marsman
, and
G.
Kresse
,
J. Chem. Phys.
122
,
234102
(
2005
).
14.
P. J.
Feibelman
,
Science
295
,
5552
(
2002
).
15.
D.
Asthagiri
,
L. R.
Pratt
, and
J. D.
Kress
,
Phys. Rev. E
68
,
041505
(
2003
).
16.
E.
Schwegler
,
J. C.
Grossman
,
F.
Gygi
, and
G.
Galli
,
J. Chem. Phys.
121
,
5400
(
2004
).
17.
J.
VandeVondele
,
F.
Mohamed
,
M.
Krack
,
J.
Hutter
, and
M.
Parrinello
,
J. Chem. Phys.
122
,
014515
(
2005
).
18.
H.-L.
Sit
and
N.
Marzari
,
J. Chem. Phys.
122
,
204510
(
2005
).
19.
H. S.
Lee
and
M.
Tuckerman
,
J. Chem. Phys.
125
,
154507
(
2006
).
20.
A. E.
Mattsson
and
T. R.
Mattsson
(unpublished).
21.
A. E.
Mattsson
,
P. A.
Schultz
,
M. P.
Desjarlais
,
T. R.
Mattsson
, and
K.
Leung
,
Modell. Simul. Mater. Sci. Eng.
13
,
R1
(
2005
).
22.
A. E.
Mattsson
,
R.
Armiento
,
P. A.
Schultz
, and
T. R.
Mattsson
,
Phys. Rev. B
73
,
195123
(
2006
).
23.
B.
Santra
,
A.
Michaelides
, and
M.
Scheffler
,
J. Chem. Phys.
127
,
184104
(
2007
).
24.
R.
Armiento
and
A. E.
Mattsson
,
Phys. Rev. B
66
,
165117
(
2002
).
25.
D. M.
Ceperley
and
B. J.
Alder
,
Phys. Rev. Lett.
45
,
566
(
1980
).
26.
J. P.
Perdew
and
Y.
Wang
,
Phys. Rev. B
45
,
13244
(
1992
).
27.
J. P.
Perdew
and
A.
Zunger
,
Phys. Rev. B
23
,
5048
(
1981
).
28.
W.
Kohn
and
A. E.
Mattsson
,
Phys. Rev. Lett.
81
,
3487
(
1998
).
29.
A. E.
Mattsson
and
W.
Kohn
,
J. Chem. Phys.
115
,
3441
(
2001
).
30.
J. P.
Perdew
,
J.
Tao
, and
R.
Armiento
,
Acta Phys. Chim. Debrecina
36
,
25
(
2003
).
31.
N. D.
Lang
and
W.
Kohn
,
Phys. Rev. B
1
,
4555
(
1970
).
32.
L.
Vitos
,
B.
Johansson
,
J.
Kollár
, and
H. L.
Skriver
,
Phys. Rev. B
62
,
10046
(
2000
).
33.
G. H.
Gonnet
,
D. E. G.
Hare
,
D. J.
Jeffrey
, and
D. E.
Knuth
,
Adv. Comput. Math.
5
,
329
(
1996
).
34.
Z.
Yan
,
J. P.
Perdew
, and
S.
Kurth
,
Phys. Rev. B
61
,
16430
(
2000
).
35.

Subroutines for AM05 are available from two of the authors (A.E.M. and R.A.).

36.
S.
Kurth
,
J. P.
Perdew
, and
P.
Blaha
,
Int. J. Quantum Chem.
75
,
889
(
1999
).
37.
Z.
Wu
and
R. E.
Cohen
,
Phys. Rev. B
73
,
235116
(
2006
);
F.
Tran
,
R.
Laskowski
,
P.
Blaha
, and
K.
Schwarz
,
Phys. Rev. B
75
,
115131
(
2007
).
38.
J. P.
Perdew
,
A.
Ruzsinszky
,
G. I.
Csonka
,
O. A.
Vydrov
,
G. E.
Scuseria
,
L. A.
Constantin
,
X.
Zhou
, and
K.
Burke
(unpublished).
39.
J.
Tao
,
J. P.
Perdew
,
V. N.
Staroverov
, and
G. E.
Scuseria
,
Phys. Rev. Lett.
91
,
146401
(
2003
);
[PubMed]
V. N.
Staroverov
,
G. E.
Scuseria
,
J. T.
Tao
, and
J. P.
Perdew
,
Phys. Rev. B
69
,
075102
(
2004
).
40.
X.-Y.
Pan
,
V.
Sahni
, and
L.
Massa
,
Int. J. Quantum Chem.
107
,
816
(
2006
).
41.
J.
Paier
,
M.
Marsman
, and
G.
Kresse
,
J. Chem. Phys.
127
,
024103
(
2007
).
42.
Y.
Zhang
and
W.
Yang
,
Phys. Rev. Lett.
80
,
890
(
1998
).
43.
L. A.
Curtiss
,
K.
Raghavachari
,
P. C.
Redfern
, and
J. A.
Pople
,
J. Chem. Phys.
106
,
1063
(
1997
).
44.
A. D.
Becke
,
J. Chem. Phys.
98
,
1372
(
1993
).
45.
A. D.
Becke
,
J. Chem. Phys.
98
,
5648
(
1993
).
46.
J.
Harris
and
R. O.
Jones
,
J. Phys. F: Met. Phys.
4
,
1170
(
1974
).
47.
O.
Gunnarsson
and
B. I.
Lundqvist
,
Phys. Rev. B
13
,
4274
(
1976
).
48.
D. C.
Langreth
and
J. P.
Perdew
,
J. Phys. F: Met. Phys.
15
,
2884
(
1977
).
49.
50.
M.
Ernzerhof
,
Chem. Phys. Lett.
263
,
499
(
1996
).
51.
M.
Ernzerhof
,
J. P.
Perdew
, and
K.
Burke
,
Int. J. Quantum Chem.
64
,
285
(
1996
).
52.
M.
Ernzerhof
and
G. E.
Scuseria
,
J. Chem. Phys.
110
,
5029
(
1999
).
53.
V. N.
Staroverov
,
G. E.
Scuseria
,
J.
Tao
, and
J. P.
Perdew
,
J. Chem. Phys.
119
,
12129
(
2003
).
54.
J.
Heyd
,
J. E.
Peralta
,
G. E.
Scuseria
, and
R. L.
Martin
,
J. Chem. Phys.
123
,
174101
(
2005
).
55.
A. V.
Krukau
and
O. A.
Vydrov
,
A. F.
Izmaylov
, and
G. E.
Scuseria
,
J. Chem. Phys.
125
,
224106
(
2006
).
56.
P. E.
Blöchl
,
Phys. Rev. B
50
,
17953
(
1994
);
G.
Kresse
and
D.
Joubert
,
Phys. Rev. B
59
,
1758
(
1999
).
57.
P. E.
Blöchl
,
O.
Jepsen
, and
O. K.
Andersen
,
Phys. Rev. B
49
,
16223
(
1994
).
58.
H. J.
Monkhorst
and
J. D.
Pack
,
Phys. Rev. B
13
,
5188
(
1976
).
59.
J. M.
Wills
,
O.
Eriksson
,
M.
Alouani
, and
D. L.
Price
,
Lect. Notes Phys.
535
,
148
(
2000
);
AM05 is implemented in the version available at http://www.rspt.net
60.
F. D.
Murnaghan
,
Proc. Natl. Acad. Sci. U.S.A.
30
,
244
(
1944
).
61.
J.
Heyd
and
G. E.
Scuseria
,
J. Chem. Phys.
121
,
1187
(
2004
).
62.
H. M.
Tütüncü
,
S.
Bağci
,
G. P.
Srivastava
,
A. T.
Albudak
, and
G.
Uğur
,
Phys. Rev. B
71
,
195309
(
2005
).
63.

In practice, basis sets, k points, pseudopotentials, and other numerical approximations come into play. There is, however, a fundamental difference between the two sources of errors: Shortcomings in the XC functional cannot be mitigated by improvements in the numerical precision.

64.
AM05 is implemented in the main version of VASP 5.1 and its use is invoked by setting GGA=AM in the INCAR file.