vdW-DF-ahcx: a range-separated van der Waals density functional hybrid

Figure 1 shows exchange hole n_x descriptions for the semilocal PBE and PBEsol functionals, as well as for the truly nonlocal-correlation CX functional. The exchange part of the functionals are denoted PBEx, PBEsolx, and cx13. The hole shapes are here represented in the HJS or AH model framework for exchange [119], using a form ${J}_{\text{x,HJS}}^{\text{GGA}}(n(\mathbf{r}),s(\mathbf{r}))$. However, the plots for PBEx and PBEsolx exchange can be directly compared with the EP-based analysis, ${J}_{\text{x,EP}}^{\text{GGA}}(n(\mathbf{r}),s(\mathbf{r}))$ [132], shown in figure 2 of reference [132] and figure 3 of reference [161]; there are no discernible differences. There are also just small differences from the original PBEx and PBEsolx hole representations (apart from the removal of Friedel-oscillation signatures).

Zoom In Zoom Out Reset image size

We work with the HJS framework for exchange holes because it gives several advantages. We call it an AH model of GGA exchange because it permits an explicit evaluation of equation (9) once ${J}_{\text{x,HJS}}^{\text{GGA}}(n(\mathbf{r}),s(\mathbf{r}))$ is inserted in equation (9). The evaluation is stated in reference [119] and it also permits a straightforward definition of functional derivatives. This is a clear advantage when the AH analysis is used to also define and code (as summarized in the following section) an RSH hybrid HJS-PBE [119] that mirrors HSE. The new design simply involves changing from the ${J}_{\text{x,EP}}^{\text{GGA}}$ form to the ${J}_{\text{x,HJS}}^{\text{GGA}}$ form.

A key AH advantage is ease of generalization. It follows from having an analytical evaluation of equation (9). For example, HJS immediately established a plausible exchange-hole shape that reflects the PBEsol exchange enhancement factor [119]. In turn, this AH determination led to the definition of HJS-PBEsol [119, 128], that is, an RSH based on PBEsol. We shall use the AH-model for seeking generalization to the vdW-DF method in the following subsection.

To begin, we summarize the HJS AH framework for characterizing the PBE and PBEsol exchange nature. HJS revisits the mathematical framework for representing an LDA exchange hole

$$\begin{align}{J}_{\text{HJS,x}}^{\text{LDA}}(y)& =-\frac{9}{4{y}^{4}}(1-{\mathrm{e}}^{-\bar{\mathcal{A}}{y}^{2}})\\ & \quad +\left(\frac{9\bar{\mathcal{A}}}{4{y}^{2}}+\mathcal{B}+\mathcal{C}{y}^{2}+\mathcal{E}{y}^{4}\right){\mathrm{e}}^{-\mathcal{D}\enspace {y}^{2}},\end{align} \tag{ 10 }$$

in a form that is free of signatures of Friedel oscillations. The parameters $\bar{\mathcal{A}},\mathcal{B},\mathcal{C},\mathcal{D}$ and $\mathcal{E}$ are listed in the HJS reference [119]. Like in the EP model, the parameters in equation (10) are set by constraints on the exchange hole.

The next step in the HJS exchange specification involves an extension to the case of gradient corrections. In both the EP [132] and the HJS [119] hole models, this gradient-extension form retains the Gaussian form. In the HJS form, the scaled exchange hole is written

$$\begin{equation}{J}_{\text{HJS,x}}^{\text{GGA}}(s,y)=\mathcal{I}(s,y){\mathrm{e}}^{-{s}^{2}\mathcal{H}(s)\enspace {y}^{2}}\end{equation} \tag{ 11 }$$

$$\begin{align}\mathcal{I}(s,y)& =\left(\frac{9\bar{\mathcal{A}}}{4{y}^{2}}+\mathcal{B}+\mathcal{C}\bar{\mathcal{F}}(s){y}^{2}+\mathcal{E}\bar{\mathcal{G}}(s)\enspace {y}^{4}\right){\mathrm{e}}^{-\mathcal{D}\enspace {y}^{2}}\\ & \quad -\frac{4}{9{y}^{4}}(1-{\mathrm{e}}^{-\bar{\mathcal{A}}{y}^{2}}).\end{align} \tag{ 12 }$$

The $\bar{\mathcal{F}}(s)$ function (affecting the y² term of the LDA description) has the form

$$\begin{equation}\bar{\mathcal{F}}(s)=1-\frac{1}{27C}\frac{{s}^{2}}{1+{(s/2)}^{2}}-\frac{1}{2C}{s}^{2}\mathcal{H}(s).\end{equation} \tag{ 13 }$$

That is, it is given by the shape of $\mathcal{H}(s)$. Meanwhile, the final $\bar{\mathcal{G}}(s)$ modification (off of the LDA description) is set by an overall charge conservation criteria [23, 25, 145–147], i.e., a limit on the integral of the exchange hole.

The important observation is that all components of the AH exchange description, equation (12), are completely defined by the shape of $\mathcal{H}(s)$. The relation between $\bar{\mathcal{G}}(s)$ and $\mathcal{H}(s)$ can be explicitly stated in the HJS AH exchange model, see equation (40) of reference [119]. Overall, the gradient-corrections provide, effectively, an extra Gaussian damping (defined by $\mathcal{H}(s)$), compared to the LDA description. The suppression arises at large values of the scaled distance y. This damping is, however, offset by an enhancement of the exchange hole at and around y = 1 (with the enhancement growing significantly with s). Nevertheless, all exchange hole components are completely set once we pick a plausible form for $\mathcal{H}(s)$.

In practice, the HJS AH modeling, for a given exchange of the GGA family, proceeds in 3 steps. First, we assume a rational form of the Gaussian damping function

$$\begin{equation}\mathcal{H}(s)=\frac{\sum\limits _{n=2}^{7}{a}_{n}\enspace {s}^{n}}{1+\sum\limits _{m=1}^{9}{b}_{m}\enspace {s}^{m}}\end{equation} \tag{ 14 }$$

and make a formal evaluation of equation (9):

$$\begin{align}{F}_{\text{x}}(s)& =\bar{\mathcal{A}}-\frac{4}{9}\frac{\mathcal{B}}{\lambda }-\frac{4}{9}\frac{\mathcal{C}\bar{\mathcal{F}}(s)}{{\lambda }^{2}}-\frac{8}{9}\frac{\mathcal{E}\bar{\mathcal{G}}(s)}{{\lambda }^{3}}\\ & \quad +\zeta \enspace \mathrm{ln}\left(\frac{\zeta }{\lambda }\right)-\eta \enspace \mathrm{ln}\left(\frac{\eta }{\lambda }\right).\end{align} \tag{ 15 }$$

Here $\zeta ={s}^{2}\mathcal{H}(s)$, $\eta =\bar{\mathcal{A}}+\zeta$, and $\lambda =\mathcal{D}+\zeta$. Second, since the GGA exchange is fully specified by the enhancement factor F_x(s), we can invert equation (15) to get a numerical representation of $\mathcal{H}(s)$. Finally, we fit that numerical determination to the rational form, equation (14).

This HJS procedure formally works for any GGA exchange—but it is an implicit assumption that we can still view the AH form J as reflecting an implicit, constrained maximum-entropy principle [26, 132, 152, 160]. The Gaussian damping form, $\mathcal{H}(s)$, must reproduce a given GGA-type exchange enhancement factor F_x(s) while using as little structure as possible; there is no physics rationale for keeping such structure [132, 160]. In an HJS representation, equation (14), we must strive to use as few significant polynomial coefficients as possible. We must also check that we have, indeed, avoided fluctuations that could hide non-physics input.

The left and middle panel of figure 1 contrast the dependence on scaled density gradient s values of our AH models for PBEx and PBEsolx, as refitted here by us. The AH model for PBEsol exchange has a softer variation with the scaled separation y and with the scaled gradient s. In contrast, the AH for PBE exchange deepens more quickly. The GGA hole formulation is best motivated when the value of the scaled hole form J remains larger than −1 (as it reflects a suppression of the electron occupation). The PBE and PBEsol exchange hole descriptions comply with that criteria for s < 2.5, i.e., in regions that are expected to dominate the description of bonding in materials and often also among molecules [20, 26, 27, 56, 60, 160, 162].

The bottom panel of figure 2 shows the shape of the exponential suppression or damping factor $\mathcal{H}(s)$ for PBEx and PBEsolx (and for the cx13 exchange functional that is relevant for CX). For the PBEx analysis there is a close overlap with the exchange hole described in the original EP model, reference [132], as also observed in the HJS paper, reference [119]. Similarly, the PBEsolx hole matches the EP-model-based representation asserted in reference [161].

Zoom In Zoom Out Reset image size

The AH hole characterizations should be seen as trusted models of the PBEx and PBEsolx hole. The HJS AH analysis of PBEx and PBEsolx also have nearly identical forms: alone the scripted parameters denoted by an over bar in equation (12) differs slightly in their formal expression. There are difference in the precise formulation of $\mathcal{H}(s)$. However, the resulting AH descriptions of PBEx and PBEsolx can claim an important maximum entropy status because of the close similarity with the EP descriptions.

An important difference between the EP and the HJS exchange-hole descriptions lies in how these models are being used. The EP model [132] was introduced to discuss the enhancement factor PBEx in terms of a non-oscillatory exchange hole, and it led to HSE [116, 117]. The HJS model [119] allows for a simpler PBE discussion, but it is also being used for reverse engineering, i.e., being used to assert a plausible exchange-hole shape from the exchange enhancement factor. This track gave rise to both HJS-PBEsol [119, 128] and several recent long-ranged corrected hybrids [134].

This paper formulates the AHCX RSH (below) from an expectation that the HJS procedure (for reverse engineering an exchange hole form) also remains valid for both rPW86 and the LV exchange descriptions. This is plausible because cx13 is formed as a Padé interpolation of the LV exchange [54, 99, 101] (at small to medium values of the scaled density gradient s) and of rPW86 [25, 142] (at large s values).

We first note that the PW86 exchange paper, reference [25], introduced a design strategy for setting the exchange hole within the electron gas tradition. The procedure involves four steps: (1) establish the symmetry limits on how the density gradient modifies the exchange hole off of the well-understood LDA representations; (2) establish the small-s variation, for example, from MBPT results as in LV exchange; (3) impose an overall exchange-hole charge conservation criteria; and (4) extract the exchange enhancement factor from equation (9). Input beyond such physics analysis is always minimized in the electron-gas tradition that defines LDA, PW86/rPW86 exchange as well as both PBE, PBEsol and the vdW-DF method.

We focus our wider AH analysis on the exchange choices that are used in the Chalmers–Rutgers vdW-DF releases. The exchange forms are the revPBEx (used in vdW-DF1), the rPW86 (used in vdW-DF2), and cx13 (used in CX). One of these, rPW86 [142], is a refit of the PW86 that defined the strategy for making constraint-based GGA exchange design off of an model exchange hole [25–27, 116, 119, 128]. We shall give a detailed motivation for trusting the HJS-type AH exchange characterizations for rPW86 and cx13, below.

We have fitted the parameters a_2–7 and b_0–9 of equation (14) for the AH representations of the cx13, revPBEx, and rPW86. This extends our PBEx and PBEsolx discussion in the previous subsection. Additional information is available in the SI material. Table 1 summarizes the most important such parameterization results: our AH characterization for cx13 exchange and our AH refits of the HJS descriptions for PBEx and PBEsolx.

The right panel of figure 1 provides a practical motivation for trusting our AH analysis of the cx13 (LV-rPW96) exchange hole. The argument is given by noting similarities to the exchange-hole shapes of well-established exchange functionals. The cx13 shape begins (at low s) by reflecting a PBEsol nature and eventually it rolls over to the PBE-type (and rPW86-type) variation. Values greater than s = 2.5 lead to the deeper hole minima, as in PBE. The cross over is not surprising since the PBEsol builds on an MBPT analysis that is close to the LV description (with small exchange enhancements up to medium s values) while the rPW86 is known to have deep exchange holes.

Figure 2 summarizes the exchange enhancements (top panel) and our AH analysis of corresponding exchange holes (bottom panels). The top panel confirms, in terms of exchange enhancements, that cx13 starts with a PBEsolx- (and PBEx-) like behavior at small s values but transforms to the rPW86 behavior at large s values.

A more formal argument for trusting the AH analysis of cx13 exchange can be stated as follows. We note that the CX leverages the formally exact electrodynamic-response framework of the vdW-DF method, as far as it is possible [60]. The CX explicitly enforces current conservation up to the cross over in the cx13 (or LV-rPW86) Padé construction, namely at s ≈ 2.5. As such, for the lower-s LV end (up to 2–3), CX is an example of a consistent vdW-DF [60], systematically relying on a plasmon-based response description that adheres to all known constraints and sum rules [51, 58]. At the LV end, cx13 furthermore leverages the Lindhard screening logic to balance exchange and correlation [60]. There is, in fact, a strong formal connection to the PBE constraint-based design logic in that the CX emphasis on current conservation ensures an automatic compliance with the charge-conservation criteria [60].

The cx13 or LV-rPW86 is a well-motivated electron-gas construction at this low-to-mid-s LV end. Here the cx-13 also satisfies the PBE-type local implementation [26] of the Lieb–Oxford bound [163], systematically having exchange-enhancement factors F_x(s) < 1.804. For s < 2.5, the cx13 (and hence CX design) complies with the criteria on the exchange hole depth, J > −1.

At the high-s end, the cx13 exchange design rolls over into the rPW86 [142]. It is there naturally within the realm of constraint-based exchange descriptions that started with PW86, reference [25]. However, in this end, we must also discuss compliance with the actual, globally implemented, Lieb–Oxford bound [163]. This is because the PW86 was designed before the importance of the bound was understood [26, 27, 160].

For any given ground-state electron density n, we consider the ratio R of the total exchange to the total LDA exchange (in the unit cell). Rigorous bounds are R < 1.804 (as is hard wired at the local level in PBE [26]) and R < 1.174 for two-electron systems (as is hard wired at the local level of SCAN [46, 164]). We note that there are unphysical (spherical-shell) systems where the density is constructed so that the scaled density gradient is constant and where the actual Lieb–Oxford bound can only be satisfied if implemented also at the local level (as in PBE and SCAN) [165]. For such unphysical systems, the rPW86 and cx13 (and hence CX) will fail to comply with this exchange condition [163, 165], if we furthermore assume that the unphysical system has a large scaled-density value, s > 3.

The important question remains, however, whether the cx13 (and hence CX) violates the actual Lieb–Oxford bound [163] in real systems, that is, as encountered in actual DFT studies [166]. The point is, that the Lieb–Oxford bound is formulated globally and must be checked on the unit-cell level in our periodic-system calculations [163].

Answering this question is straightforward for any given problem, as long as one saves the density after completing a Quantum Espresso calculation. Our now updated ppACF post-processing code (launched in references [123, 124] and committed to the Quantum Espresso package) gives a general coupling-constant analysis, an evaluation of kinetic-correlation energy components, as well as a per-system determination of actual exchange R ratio (for nonhybrids). For the 10 fullerene calculations discussed in the introduction, we find that the cx13/CX value for the actual exchange ratio R never exceeds 1.025. In practice, the CX exchange ratio remains far below even the stringent bound (1.174) that exists for two-electron systems [164].

The bottom panel of figure 2 provides further details of the cx13 design. It does so by comparing the Gaussian suppression factors $\mathcal{H}(s)$ that arise in our AH analysis of constrained exchange functionals. The shapes of the Gaussian suppression $\mathcal{H}(s)$ are similar in all cases. The s position of the maximum for cx13 (LV-rPW86) sits essentially on the PBEsolx position, slightly above the maximum position for rPW86 and PBEx. Meanwhile, the cx13 maximum value aligns with that of the PBE description. As expected, the asymptotic cx13 behavior coincides with that of the rPW86 while for low s values, the Gaussian suppression is almost identical to that which characterizes PBEsol. This is also expected as the PBEsolx was explicitly designed to move the exchange enhancement closer to the input that also defines LV exchange. That is, the PBEsolx is closer to the MBPT analysis of exchange in the weakly perturbed electron gas [27, 54, 60, 99, 101].

Overall, figure 2 confirms that there are no wild features in the AH characterizations of exchange in the vdW-DF releases, including CX. Moreover, the bottom panel confirms our assumption, that our AH model description of cx13 can be trusted. This follows because we find that it reflects a mixture of PBEsol-like, PBE-like, and rPW86-like behaviors for the Gaussian suppression factor of the hole. In fact, like PBE, the cx13 aims to serve as a compromise of staying close to MBPT results at low s (good for solids) and a more rapid enhancement rise (as in rPW86) for larger s values (good for descriptions of molecular binding energies [142, 167]).

Table 2 reports a simple sanity check that we provide to argue robustness of our maximum-entropy description of the AH model and RSH constructions. We use the fact that reference [119] reports parameters for their AH characterization of PBE exchange, a form that is here termed HJS-PBE. We provide a modeling self-test using the ONCV-SG15 PPs at 160 Ry, noting that the details of the $\mathcal{H}(s)$ specification must not change the RSH results.

Table 2 contrasts the results of such HSE-like descriptions for 6 atomization energies. For the HJS-PBE and PBE-AH we chose values of the Fock-exchange mixing α = 0.25 and for the screening γ = 0.106 (inverse Bohr's) that exist for default HSE runs [117]. The agreement with reference data is good. The alignment of HSE and HJS-PBE (or AH-PBE) descriptions is very good.

Most importantly, table 2 illustrates that our code to extract AH descriptions for PBE fully aligns results of PBE-AH with those of HJS-PBE. We trust our AH model description of PBE exchange, as defined by parameters in table 1.

Interestingly, the details of our $\mathcal{H}$ specifications (table 1) differ from those provided by HJS [119]. This is true even if there is no difference in the resulting RSH descriptions. This finding suggests robustness in the overall AH constructions of RSH, not only for the HJS-PBE but also for the here-defined vdW-DF-based RSH, the AHCX.

For a survey of performance on molecular properties we have turned to the large parts of GMTKN55 that are immediately available for a planewave assessment using the ONCV-SG15 PP set [139, 140]. The full GMTKN55 is organized into 6 groups of benchmarks: (1) small molecule properties, (2) large molecular properties, (3) barrier properties, (4) inter-molecular noncovalent interactions (NCIs), (5) intra-molecular NCIs, and (6) total NCIs. However, we exclude the G21EA and WATER27 benchmark sets of GMTKN55 [120]. In effect, we focus on a subset, the remaining easily-accessible 'GMTKN53', in our discussion.

The G21EA benchmark is excluded because it contains negative charging of ions and small radicals (like OH⁻). Such small negatively charged systems genuinely challenge our planewave assessment of broad molecular properties exactly because we seek to maintain a high precision. The fundamental problems arise by a combination of two factors. First, the SIEs will, in these systems, push the highest-occupied molecular-orbital level up toward or above the vacuum floor, when it is done in a complete-basis set description [121]. Second, our planewave description rapidly approached that complete-basis set limit due to our efforts to also carefully control spurious inter cell electrostatics interactions down toward the 0.01 kcal mol⁻¹ limit [178]. We conclude that a direct performance assessment on G21EA is meaningless in a planewave code and it is perhaps best left for a more advanced handling [121].

We note in passing that we can approximately assert the G21EA performance by using the planewave equivalent of the so-called-moderate-basis-set approach [121], trapping the negatively charge ions in a small 6-to-8 Å-cubed boxes. We can thus obtain an approximate comparison that shows that the variation in G21EA performances has no discernible impact on the statistical performance measures that we track in GMTKN53.

We furthermore removed WATER27 for a related SIE-induced [121] problem—the path to electronic-structure convergence for the OH⁻ is exceedingly cumbersome once we seek to fully converge the WATER27 set with respect to unit-cell size. Of course, hybrids helps in the WATER27 (and in the G21EA) study—but then we can offer no fully converged comparisons. Again, using a separate (12 Å) small-box handling of the OH⁻ system makes it possible to complete an approximate WATER27 benchmarking; as will be detailed elsewhere, this assessment identifies vdW-DF2, CX, CX0p and AHCX as top performers.

Fortunately, having 53 (out of 55 sets of reference [120]) still makes for ample statistics. Accordingly, we simply work with the easily-accessible GMTKN53 and compare with literature GMTKN55 results [120]. This is fair because, if anything, we have negatively offset our AHCX performance reporting (by omitting G21EA and WATER27).

We supplement the GMTKN53 characterization both by tracking performance differences (that arise for some individual benchmarks) and by analyzing the groups of benchmarks that mainly reflect NCIs. We essentially use the same computation setup as the pilot study included in a recent focused review of consistent vdW-DFs, reference [60]. However, we have now moved off of the electron-sparse AbInit PPs and onto the electron-rich ONCV-SG15 [139, 140]; moreover, we now make systematic use of large cubic cells and a Makov–Payne electrostatics decoupling [179]; in the pilot assessment [60], the electrostatic decoupling was only done for charged systems.

We find that this upgrade of our benchmarking strategy accounts for some of the deviations that we previously incorrectly ascribed to limitations of the CX [56] and vdW-DF2-b86r [180] functionals [60]. The past limitations seem to, relatively speaking, especially impact our assessment of group 4 (intermolecular NCI) performance, as represented in references [42, 60]. This paper provides a more fair performance characterization of the performance of CX and vdW-DF2-b86r as well as of AHCX relative to other vdW-inclusive DFT versions, for example, dispersion-corrected metaGGA.

For bulk structure and cohesive energies, we used the Monkhorst–Pack scheme with an 8 × 8 × 8 k-point grid for the Brillouin-zone integration. Our demonstration survey includes 4 semiconductors, 3 ionic insulators, 1 simple metal, and 5 transition metals including Cu, Ag, and Au. Here we also used the electron-rich ONCV-SG15 [139, 140] PPs set at 160 Ry wavefunction energy cut-off. We tracked the total-energy variation with lattice constant and used a fourth-order polynomial fit [181] to identify optimal (DFT) structure, the cohesive energy, and bulk modulus.

Surface energies for Cu(111), Ag(111), and Au(111) were estimated from the energies calculated from 4 to 12 layers of the noble metals in unit-cell configurations with 15 Å vacuum. Here we sought to offset some of the high computational cost by instead using the AbInit normconserving PPs [176] for extra demonstrations of AHCX usefulness. In the surface studies, we first determine the in-surface unit cell from the calculated bulk lattice constant. Next we allow the outermost atoms to relax in a slab-geometry description. The exception is for AHCX characterizations, where we instead relied on the results of the CX structure characterizations.

To document AHCX accuracy for CO adsorption on noble-metal surfaces, we used again the electron-rich ONCV-SG15 PPs at 160 Ry. We consider the FCC and TOP sites in a 2 × 2 surface-unit-cell description and a six-layer slab geometry. We also track the adsorption-induced relaxations on the top three layers and compute the CO adsorption energy as follows,

$$\begin{equation}{E}_{\text{ads}}={E}_{\text{CO/metal}(111)}-{E}_{\mathrm{m}\mathrm{e}\mathrm{t}\mathrm{a}\mathrm{l}(111)}-{E}_{\text{CO}}.\end{equation} \tag{ 23 }$$

Here E_{CO/metal(111)} denotes the fully relaxed energy of the adsorbate system, while E_metal(111) and E_CO denote the total energy of the fully-relaxed (isolated) surface and molecule system, respectively.

Figure 3 and table S.II of the SI material detail the computation time required for completing one evaluation of the ACE Fock-exchange operator in hybrid calculations on a single 20-core node for gold and diamond, top and bottom panels, respectively. We include a comparison of HSE [117], CX0p [123], and AHCX for gold and diamond, using an 8 × 8 × 8 k point grid and the electron-rich ONCV-SG15 PPs. We track the variation with the so-called q mesh (of k point differences used in the Fock-exchange evaluation).

Zoom In Zoom Out Reset image size

There is no additional cost of doing AHCX over HSE (as in the molecular cases discussed in the introduction) and the AHCX cost is slightly lower than the cost of the corresponding unscreened hybrid CX0p. This is in part a statement of the computational efficiency of the new nonlocal-correlation RSH functional. It is also a statement that we have provided a robust implementation of the analytical hole description and that we have succeeded in leveraging the benefits of the Quantum Espresso ACE implementation [136].

An overview of scaling of AHCX (or HSE) hybrid calculations can be deduced from reference [136]. The cost per Fock exchange evaluation scales with the FFT scope (size of unit cell) and with the number of bands (that is, number of electrons per unit cell) squared [136]. Also, for hybrid studies of bulk and surfaces one needs to consider multiple k points and a q mesh of k-point differences, causing an additional scaling factor. When there is symmetry (as in simple bulk) one can cut down this factor significantly by limiting the set of k and q points that enters in an actual Fock-exchange evaluation. Importantly, however, there are communication costs when jobs are spread across multiple nodes [136]. For example, we can push more hybrid calculations with electron-sparse PPs because they require fewer bands; by limiting bands we both directly accelerate the computation speed and we reduce the memory requirements, making our jobs fit on more types of computer resources.

The AHCX is more costly than non-hybrids, including metaGGAs, for example, as illustrated in the introduction. Also, as mentioned above, hybrid calculations are costly for surfaces where relaxations enviably produce low-symmetry geometries. However, the AHCX computational costs can be handled by choosing an appropriate parallelization strategy. We try to limit the number of nodes as we also try to accommodate the memory requirements.

For a specific illustration of AHCX costs we make the following observations from our studies. We report summary of nearly 5000 AHCX and HSE-type hybrid studies of molecular properties of the GMTKN55 (having up to 37 Å-cubed unit cells to carefully control electrostatic couplings); those ONCV-SG15 jobs took us a total of 200 000 core hours (5–6 times the cost of using the corresponding CX and PBE regular functionals). The cost for an individual bulk structure study in AHCX/HSE is (due to symmetry) relatively cheap (see figure 3); here the 700+ computations of energy-versus structure variations that enters our structure optimization (of 11 simple cases and of 5 transition metals) took just over 100 000 h at the production stage. Finally, the cost for hybrid studies of noble-metal surface energies and adsorption are on a different scale. We were able to obtain the 6 frozen-geometry AHCX CO-noble-metal adsorption energies in the electron-rich ONCV-SG15 PP setup at the costs of about 100 000 core-hours each (on selected computer resources). For the additional noble-metal workfunction and surface energy demonstrators (involving many slabs) we cut the cost to about 1.5 million core hours by instead relying on the electron-sparse AbInit PPs for AHCX calculations.

Figure 4 shows a performance comparison that we provide for 10 molecular benchmarks of the full GMTKN55 suite. The figure characterizes the performance in terms of MAD values. Use of hybrid AHCX provides significant improvement over CX in the case of self-interaction problems (in the SIE4x4 benchmark set), for atomization energies (in the W4-11) set, and for several barrier-type problems, including those reflected in the PX13 and WCPT18 benchmark sets. The improvements are perhaps largest in the case of C60 isomerizations, for large deformations having stretched bonds.

Zoom In Zoom Out Reset image size

Figure 5 summarizes our survey of performance across the entire NCI groups of the GMTKN55, although excluding the WATER27 benchmark set for reasons explained above. The figure characterizes the performance in terms of a total MAD, or TMAD, value for intermolecular NCIs (abscissa position) and of a TMAD value for intramolecular NCIs (ordinate position). These effective MAD values result as we first compute the MAD values for each NCI benchmark of the GMTKN55 and then average over these MADs over the number of intermolecular and intramolecular NCI benchmarks investigated in GMTKN55 group 4 and group 5, as in reference [42]. The SI material provides details and lists the quantitative data that we provide for vdW-DF hybrids (AHCX and CX0p, circles) and for regular vdW-DFs (vdW-DF2, vdW-DF2-b86r and CX, upwards triangles), as well as for rVV10 and dispersion corrected revPBE-D3 and HSE-D3. The SI material also includes a performance characterization of vdW-DF1 [51, 52, 54] and vdW-DF-ob86 [182].

Zoom In Zoom Out Reset image size

The NCI assessments in figure 5 can be compared to figure 3 of reference [42], a study that also summarized NCI results obtained in references [120, 183]. The assessments differ in code nature (planewave versus orbital-based DFT) but we can use figure 5 for a broader discussion of performance. That is, our comparison for NCI performance also indicates how the new AHCX RSH vdW-DF hybrid fares in relation to the dispersion-corrected metaGGAs and dispersion-corrected hybrids discussed in reference [42]. We find that the AHCX performance is better than those of the dispersion-corrected metaGGAs and of ωB97X-D3 (as reported in reference [120, 183]). However, our AHCX does not fully match the performance of DSD-BLYP-D3 [120], when it comes to the TMAD value for intermolecular NCI.

For completeness, figure 5 also shows the approximate assessments (downwards triangles) of CX and vdW-DF2-b86r performance that we obtained in a pilot benchmarking [60]. The differences in CX and vdW-DF2-b86r positions of upwards and downward filled triangles reflect the improvement that we have here made in benchmarking strategy relative to reference [60]. The systematic use of decoupling of spurious dipolar inter-cell attractions is found to have a clear effect on performance for the group of intermolecular NCI benchmark sets.

We find that the use of vdW-DF hybrids leads to significant improvements in the description of problems characterized by strong NCIs. The performance of AHCX and CX0p are again nearly identical but significantly better than dispersion-corrected HSE-D3 and of the regular vdW-DFs. This is not surprising because the hybrids are expected to have a more correct description of orbitals and reshaping orbitals (and hence the density variation) will directly affect the strength of the vdW attraction [34, 184–189].

For a summary of overall progress that the hybrid vdW-DFs may bring for general molecular properties, we rely on the 'easily available GMTKN53', defined above. This is the subset of the GMTKN55 suite that omits the WATER27 and G21EA benchmarks sets, for good reasons [121]. Fortunately, with GMTKN53 we remove just one benchmark set in two groups and can still provide a balanced account. It is still meaningful to discuss functional performance in terms of per-group results as in the GMTKN55 study [120].

For comparison of performance we therefore use the reference geometry and reference data of reference [120] to determine the deviation on each of the over 2000 molecular processes that are included. We evaluate the weighted TMAD measures (WTMAD1 and WTMAD2) that Grimme and co-workers introduced to allow a comparison also among the GMTKA55 groups [120]. The adaptation of these measures to the present GMTKN53 survey simply involves adjusting the counting of participating benchmark sets for GMTKN55 groups 1 and 4 [120], see SI material for further details.

Figure 6 shows our summary functional comparison based on this 'easily accessible GMTKN53', using WTMAD1 measures. The SI material lists the data for both WTMAD1 and WTMAD2 comparisons. We find that the AHCX performance is significantly better than dispersion-corrected HSE-D3 for this assessment of broad molecular properties. The AHCX also performs systematically better than CX, which in turn performs at the level of dispersion-corrected revPBE-D3. The latter is reported to be one of the best performing dispersion corrected GGAs [120].

Zoom In Zoom Out Reset image size

Importantly, the AHCX also performs at the same level and perhaps slightly better than the unscreened CX0p hybrid. Our primary motivation for defining AHCX is to have a vdW-DF RSH that can also describe adsorption at metallic surfaces, which is not something that the unscreened CX0p is set up to do. The comparison in figure 6 shows that the AHCX is also good enough on the molecular side of interfaces. That is, it remains a candidate for also serving us for the interface problems, at least so far.

Figure 7 contrasts the performance of hybrid functionals HSE and AHCX with those of the underlying PBE and CX functionals for bulk: 6 metals (elements Al, Cu, Rh, Ag, Pt, Au), 4 semiconductors (Si, C, SiC, and GaAs) and 3 ionic insulators (MgO, LiF, and NaCl). Here we compare deviations from measurements of lattice constants a and cohesive energies ΔE (both back corrected for zero-point energy and thermal effects) as well as for bulk moduli B₀. Symbols for the statistical measures are summarized in the figure caption. The PBE (AHCX) description for cohesive energies (for bulk moduli) have Ag and Au (Rh) as statistical outliers; the SI material contains a detailed presentation of the per-bulk-material performance.

Zoom In Zoom Out Reset image size

The figure shows that AHCX gives the smallest median deviation and spread of data for the deviations for lattice parameters and bulk moduli; there are clear improvements with AHCX over both HSE and CX (which is the second-best performer overall). For cohesive energies, the CX has the smallest median deviation whereas AHCX has the smallest spread and both perform better than PBE and HSE for bulk properties.

Table 3 compares the results of PBE, CX, and AHCX descriptions of transition-metal structure and cohesion to back-corrected experiments values [22] (last column). The results are provided for the electron-rich ONCV-SG15 PPs (except where noted by an asterisk '*' or when taken from literature).

Table 3. Performance assessment on bulk properties (lattice constants a and cohesive energies ΔE) for a set of late transition metals including the noble metals. The boldface entries are back-corrected experimental values [22], the rest of the entries were obtained using, PBE, CX, and AHCX with ONCV PBE-SG15 PPs at 160 Ry (except when labeled '*' where instead AbInit PPs was used at 80 Ry). Units are Å for lattice parameters and eV for cohesive energies.

System	aPBE	aCX	aAHCX	${a}_{\text{AHCX}}^{\ast }$	aref
Cu	3.639	3.576	3.587	3.591	3.599
Ag	4.156	4.065	4.078	4.116	4.070
Au	4.165	4.101	4.098	4.113	4.067
Pt	3.970	3.929	3.910	3.939	3.917
Rh	3.832	3.786	3.760	3.797	3.786

	ΔEPBE	ΔECX	ΔEAHCX	${\Delta}{E}_{\text{AHCX}}^{\ast }$	ΔEref
Cu	3.423	3.781	3.348	3.765	3.513
Ag	2.488	2.955	2.774	2.934	2.964
Au	2.997	3.634	3.440	3.752	3.835
Pt	5.434	6.226	5.524	4.876	5.866
Rh	5.565	6.367	5.244	4.257	5.783

We first observe that use of ONCV-SG15 PP setup yields a state-of-the-art benchmarking on these transition metal systems. For PBE (and for CX) we can track deviations of our self-consistent ONCV results relative to the self-consistent PBE (non-selfconsistent CX) all-electron results, provided in reference [22], and to the fully self-consistent PAW-based assessments that one of us have previously provided [166]. For lattice constants the largest PBE (CX) deviation from self-consistent (non-selfconsistent) all-electron results [22] is 0.2% (0.1%). The ONCV characterization has but minute differences from (is spot on) from the PAW-based characterization of PBE (CX) performance [166]. Similarly, for the description of cohesive energies, we find small PBE (CX) deviations, with a Rh maximum of 3.0% (1.4%), relative to all-electron results [22]. Here the ONCV characterizations are slightly less precise overall than the previous PAW-based characterizations (but have a smaller spread). Finally, for the bulk moduli, we find that our ONCV-SG15 benchmarking differs by at most 2.2% for Rh from the all-electron transition-metal results. Overall, we find that our ONCV-SG15 bulk-structure benchmarking is highly precise, for example, for characterizations of noble-metal properties.

Comparing next the calculated results against back-corrected experimental results [22], table 3 confirms the overall impression of AHCX promise (figure 7) also holds for the transition metals. The table shows that CX performs significantly better than PBE on the transition metal lattice constants (cohesive energies); this is also expected, perhaps especially for the noble metals [22, 60, 90, 166]. A move from PBE to CX reduces the maximum absolute deviation (from experimental values) on transition-metal lattice constants (cohesive energies) from 2.4% (24.2%) to 0.8% (10.1%). Meanwhile, AHCX performance is on par with that of CX for lattice constants and a clear improvement over HSE (as detailed in the SI material). This pattern is repeated for the bulk-modulus assessments, see SI material.

Overall for AHCX we find the following deviations from experimental lattice-constant (cohesive-energy) values [22]: −0.3% (−4.7%) for Cu, 0.2% (−6.4%) for Ag, 0.8% (−10.3%) for Au, −0.2% (−5.8%) for Pt, and −0.7% (−9.3%) for Rh. The performance for cohesive energies is overall slightly worse than that for CX on transition metals, but significantly better than that of PBE and HSE, see SI material. While broader tests of bulk-transition-metal performance of AHCX, for example, extending references [22, 166] are desirable, it is also clear that AHCX is useful for descriptions of bulk properties in general and certainly for noble-metal studies.

Descriptions of (organic) molecule adsorption at metal surfaces are an important reason for defining and launching the AHCX. For metal problems, we must screen the LR Fock-exchange component. The new RSH vdW-DF is set up to reflect the perfect electrostatic screening [130]. Below we begin the discussion by looking at the AHCX ability, as an RSH vdW-DF, to describe properties of noble-metal surfaces. We stay with the robust (but expensive) electron-rich ONCV-SG15 PPs (at 160 Ry) for the regular-functional descriptions and for the description of CO adsorption, next subsection; for a demonstration of AHCX's ability to describe workfunctions and surface energies, here, we shall use and discuss calculations obtained in the more electron-sparse (but still normconserving) AbInit PPs (at 80 Ry).

We note that understanding surface energies of metallic nanoparticles is important for correctly setting up cost-effective catalysis usage in catalysis. This is because they define the Wulff shapes [194] of nanoparticles and hence the extent that they provide access to specific active sites [162, 195]. As such this additional illustration could also be relevant for DFT practitioners.

Table 4 compares our PBE, CX, and AHCX results for workfunctions and surface energies against experimental values and against literature SCAN values [174]. An asterisk '*' on a column indicates that those results are obtained using the AbInit PPs at 80 Ry, the rest are obtained by ONCV-SG15 PPs at 160 Ry. All results are based on a sequence of slab calculations and we first compute bulk lattice constants and then surface relaxations specific to the functional and PP choice. As indicated by the prime on the results in the AHCX column, however, we rely on the CX structure determination for our characterization of this CX-based RSH.

The table furthermore compares against experimental values for workfunctions (as summarized in reference [174]), and surface energies for noble-metal (111) facets [60, 174]. Experimental surface-energy values are available from observations of the liquid metal surface tension [196, 197]. As such, they are defined by an average and it is natural to focus on the major facets [174]

$$\begin{equation}{\sigma }_{\mathrm{exp}}={\sigma }_{\text{avg}},\end{equation} \tag{ 24 }$$

$$\begin{align}{\sigma }_{\text{avg}}& =\frac{1}{3}({\sigma }_{(111)}+{\sigma }_{(110)}+{\sigma }_{(100)})\\ & \equiv \frac{1}{3}{\sigma }_{(111)}+\frac{2}{3}{\sigma }_{(\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r})},\end{align} \tag{ 25 }$$

where we now introduce σ_(other) = (σ₍₁₁₀₎ + σ₍₁₀₀₎)/2. Previous works have compared measured surface energies σ_exp (for various transition metals) and compared with per-facet surface-energy results—σ₍₁₁₁₎, σ₍₁₁₀₎, and σ₍₁₀₀₎—computed in various functionals, including CX [60, 162, 174, 195]. Here, we simplify the assessment task in two steps: (a) we use our past CX results [60] to define a ratio x = σ₍₁₁₁₎/σ_(other) that reflects the relative weight of (111) in the major-facet averaging [60, 174], equation (25), and (b) we extract the experimentally-guided estimate

$$\begin{equation}{\sigma }_{(111)}^{\mathrm{exp}.}\approx x\ast {\sigma }_{(\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r})}=\frac{3x}{2+x}\ast {\sigma }_{\mathrm{exp}}.\end{equation} \tag{ 26 }$$

Table 4 shows that moving between PBE, SCAN (as asserted in literature [174]), CX, and AHCX generally increases the computed values for both workfunctions and surface energies for the noble-metal (111) facets; we note that the same trend also arise (among the PBE, CX, and AHCX descriptions) when we instead compute these with the AbInit PPs at 80 Ry. It is alone for the Cu surface energy that we find an AHCX surface energy that is smaller than the CX results (AHCX gives the same Ag surface energy as CX).

Table 4 also shows that switching between the electron sparse AbInit PPs and the electron-rich ONCV-SG15 PPs has essentially no impact on the workfunction and a 5% (10%) effect on the CX results for the surface energy of Cu and Au (of Ag). These shifts are directly correlated with the fact that the Cu and Au (Ag) lattice constant in the AbInit PP description is 1% (2%) larger than what results with the ONCV-SG15 PP choice.

Table 4 suggests that both CX and AHCX provide an accurate description of the noble-metal surface energies and workfunctions. We note that a significant uncertainty exists for surface energy reference values due to the need to interpret the input from high-temperature measurements. Unlike for PBE and SCAN, the CX and AHCX surface energies fall within the range of such observation-based estimates (for Au(111), the CX, but not the AHCX, result falls just outside the error bars). For workfunctions we find that CX and SCAN are slightly more accurate than AHCX for Cu(111), while AHCX is the best performer for the Ag(111) surface; for the Au(111) workfunction the SCAN, CX, and AHCX results all fall with the broad range of experimental values.

To argue AHCX accuracy for clean-noble-metal surface properties, we make four observations. First, our use of AbInit PPs for AHCX likely impacts the comparison with the ONCV-SG15 based PBE and CX results for the surface energy σ. Second the trend from the PBE and CX results suggests that the impact on the PP choice on σ primarily arises by the inverse dependence on the surface-unit-cell area, an area-scaling that is discussed two paragraphs above. Third, appendix A assesses the sensitivity of the CO-adsorption site-preference results to the PP choice, documenting good robustness for the noble metals; this robustness is expected by the fact that these systems have a nearly closed d band [118, 198, 199] and it makes it likely that the PP impact on the AHCX σ description will be dominated by the area-scaling effect. Fourth, adjusting our AbInit-PP based AHCX characterization upwards by a 5%–10% area effect will generally push the AHCX description closer to the back-corrected experimental surface-energy values: the AHCX result for Ag will still reside on the edge of the range characterizing the Ag experimental data. In other words, table 4 shows that AHCX is useful and that it may well be a strong performer for characterizations of clean-noble-metal surface properties.

Finally, we turn to discuss a case of actual metal-surface adsorption. We focus on CO adsorption as it allows us to use just a 2 × 2 in-surface unit-cell extension and because it may provide a hint of whether AHCX ultimately may help in understanding heterogeneous catalysis. The CO adsorption on late transition metals is a classic surface-science challenge.

Figure 8 shows the CO adsorption geometry for 111 surfaces. The CO is standing upright with a Blyholder-type binding [200] through the carbon atoms (brown ball) which sits between the metal substrate (gold atoms) and the oxygen (red ball). There is no experimental ambiguity that CO molecule adsorbs on the TOP site (as illustrated in the right panel), rather than at the FCC or hollow site (left panel) at low coverage, reference [201].

Zoom In Zoom Out Reset image size

Table 5 summarizes the CO adsorption studies that we here provide using electron-rich ONCV-SG15 PPs. The most famous CO adsorption problem arises at the Pt(111) surface, references [205–208], a problem that has repeatedly challenged DFT descriptions when used for XC energy approximations that also deliver an accurate description of the in-surface Pt lattice constant and surface energies, references [173, 209–212]. Here we focus on the noble-metal adsorption cases. The table compares our results for CO adsorption at FCC and TOP sites on Cu(111), Ag(111), and Au(111). As indicated by the prime, the AHCX characterization is again done frozen at the adsorption structure that results in the CX study. We find that the AHCX description stands out by predicting the correct site preference for CO on noble-metal (111) surfaces; none of the PBE, SCAN or CX functionals have that consistency. The AHCX prediction of the actual (TOP-site) adsorption energy is good in silver and gold; all functionals have some problems in accurately predicting the magnitude of the adsorption energy for CO/Cu(111).

The AHCX functional overestimates the binding in the case of Cu(111). Comparing with the CX description, we see that the AHCX still provides a reduction of the adsorption energy. This AHCX binding softening (compared to CX) is, in part, expected. For instance, use of either PBE0 and HSE reduce the Cu(111) adsorption energy compared with PBE, reference [209].

We also highlight an important point of the AHCX success: we find that AHCX simultaneously have accurate descriptions of molecules, bulk lattice constants, surface properties, and CO adsorption energies; there is no lucky hit for one quantity at the expense of others, compare figure 6 as well as tables 3–5. The AHCX lattice constant for Cu has a 0.3% deviation from experiment, better than CX (see SI material). The AHCX descriptions of workfunctions are accurate for Cu. Ag, and Au, and the description of Cu, Ag, and Au surface energies are good. Last but not least, the AHCX performance for predicting the CO-adsorption site preference is excellent, and (except for Cu) accurate also on the predictions of the actual CO adsorption energies.