Communication: xDH double hybrid functionals can be qualitatively incorrect for non-equilibrium geometries: Dipole moment inversion and barriers to radical-radical association using XYG3 and XYGJ-OS

Density functional theory is widely used for electronic structure calculations as it tends to yield sufficiently accurate results for a significantly lower computational cost relative to correlated wave function approaches.^1,2 However, approximate local functionals are often known to fail in describing systems with charge delocalization^3–5 on account of missing important nonlocal information. Orbital dependent nonlocal exchange and correlation components are thus often hybridized with local exchange and correlation to produce more sophisticated functionals in the hope of attaining higher accuracy. Such functionals are typically classified as hybrids if they only contain nonlocal Hartree-Fock like exchange and as double hybrids (DHs) if both nonlocal orbital dependent exchange and correlation [typically from second order perturbation theory (PT2) or RPA] are included. DH functionals consequently occupy the fifth rung of Jacob’s ladder⁶ and are amongst the most accurate density functionals developed so far, yielding highly accurate predictions of energies,^7,8 dipole moments,⁹ and other properties.^10,11 This makes DH functionals very attractive for electronic structure calculations, despite a somewhat higher computational complexity due to the dependence on unoccupied orbitals for nonlocal correlation.

The energy of a DH functional constructed from a hybrid part, E_hyb, and the scaled second order PT2 correlation energy, E_PT2, may be written in two equivalent ways

The first generic form does not specify how the orbitals are optimized, while the second identifies the energy minimized via orbital optimization (OO) as $EDHOO$ plus a non-SCF correction, ΔE_DH. The first form is standard, while the second form is appropriate for considering the evaluation of molecular properties, as will be our purpose later.

With regard to choice of orbitals, three distinct OO approaches have been developed so far, which we will label tDH, xDH, and OO-DH. The truncated DH (tDH) examples like B2PLYP¹² optimize orbitals with a truncated Kohn-Sham (KS)¹³ functional that is simply $EtDHOO=Ehyb$ so that the non-SCF correction is solely the PT2 term (i.e., ΔE_tDH = E_PT2). A later approach proposed by Zhang, Xu, and Goddard¹⁴ (xDH) uses orbitals optimized by a reference (ref) lower rung functional, like B3LYP;¹⁵ so that $ExDHOO=Eref$⁠, and the non-SCF correction becomes $ΔExDH=(Ehyb−Eref)+EPT2$⁠. XYG3¹⁴ and XYGJ-OS¹⁶ are two well-known examples of this category. Finally, it is also possible¹⁷ to use the total DH energy for orbital optimization (OO-DH) in a manner similar to orbital optimized MP2^18,19 so that ΔE_OO−DH = 0. OO-DH functionals are not yet widely used.

Both the tDH and xDH approaches have certain drawbacks. The truncated KS functional used for orbital optimization in tDH, $EtDHOO=Ehyb$⁠, does not have a complete description of correlation, which will cause effects that include early onset of spin polarization. On the other hand, the xDH approach relies on the adiabatic connection,¹⁴ which may not hold if the lower rung orbital generating functional like B3LYP yields insufficiently accurate densities, as has been suggested by a number of recent studies.^9,20,21 Nonetheless, both approaches have yielded functionals highly accurate for equilibrium properties^8,9,14,16 like B2GPPLYP²² (tDH) and XYG3 (xDH).

Our recent study on dipole moments,⁹ however, suggested that the picture might not be so rosy for xDH functionals away from equilibrium. Both XYG3 and XYGJ-OS yielded dipoles on the order of 1-2 D for the HF molecule at internuclear separations of 2.75–4.25 Å, versus a coupled-cluster singles and doubles, with second order perturbation corrections [CCSD(2)]^23,24 reference value of 0.1 D or less. The polarity furthermore was reversed from H⁺F⁻ to the unphysical H⁻F⁺. Similar unphysical behavior was seen for stretched HCl, ClF, LiH, and CH₃F (along the C–F bond breaking coordinate), suggesting that such density artifacts for stretched bonds are a rather general issue with these xDH functionals. Practical implications aside, this also raises several other questions that we address here: (i) Are only densities affected or can other xDH properties also behave poorly at stretched geometries? (ii) What is the origin of this peculiar behavior? (iii) Are tDH functionals affected as well?

An examination of potential energy surfaces generated by XYG3 and XYGJ-OS for these species revealed that several have an unphysical local maximum at stretched geometries, resulting in a spurious barrier for homolytic bond formation. This is illustrated in Fig. 1 for the HF molecule. Similar features were seen even for some non-polar species like F₂. The energy barrier is not very large (≈3–15 kJ/mol) relative to the bond strength but is chemically significant. Furthermore, the local energy maxima and the unphysical dipole moments occur in roughly the same region of the dissociation curve (as seen in Fig. 1), suggesting that they may have similar origins. tDH functionals employing self-consistently optimized truncated KS orbitals like B2GPPLYP do not appear to have such features, indicating that the origin may lie in the xDH recipe itself.

We investigate the origins of this unphysical behavior through a detailed study on the HF dissociation curve. Figure 2 plots XYG3 dipoles and potential energies for this system, along with results from the reference (B3LYP), a tDH functional (B2GPPLYP), and a reliable wave function method [CCSD(2)]. While not shown, XYGJ-OS exhibits similar behavior. The dipole curves predicted by the two DH functionals have discontinuities for small bond-stretches, which arise from known N-representability violations in unrestricted MP2 and DH functionals^9,25 around the Coulson-Fischer (CF)²⁶ point of the SCF method. The CF point lies between 1.35 Å (Hartree-Fock) and 1.6 Å (B3LYP) for HF, where the sign of the dipole moment is still correct. The unphysical local extrema of XYG3 are found beyond 2 Å internuclear separation, indicating that the dipole inversion (and the energy maximum) is not directly connected to this N-representability issue, but arises from something else. The innocence of E_PT2 itself is made evident by the observation from Fig. 2 that the unphysical local extremum is more pronounced when E_PT2 is deleted. E_PT2 was therefore partially correcting these issues in the full functional.

There are two possible sources for this unphysical behavior of the non-PT2 portion, E_hyb, of the xDH functionals: either E_hyb is intrinsically unphysical or there is a major mismatch between the reference B3LYP orbitals and those which minimize E_hyb. The former can be ruled out because turning XYG3 into a tDH functional (i.e., replacing B3LYP orbitals by E_hyb optimized orbitals) yields energies and dipoles that are qualitatively acceptable, as shown in Fig. 3. Lack of self-consistency between B3LYP orbitals and E_hyb is thus the origin of these unphysical local extrema in the dissociation curves. A strong dependence of the dipole curves on the choice of reference orbitals is established in Fig. 3, where BLYP^27,28 (0% exact exchange) orbitals yield far worse behavior than XYG3’s default use of B3LYP (20% exact exchange) orbitals. On the other hand, using orbitals with higher exact exchange content (i.e., closer to XYG3’s E_hyb with ≈80% exact exchange) such as B5050LYP²⁹ (50% exact exchange) and HFLYP²⁸ (100% exact exchange) entirely removes the issue. This, in combination with an earlier study that found optimal xDH performance using orbitals optimized with 50%–70% exact exchange,³⁰ appears to indicate that reference orbitals optimized with a large fraction of exact exchange might be a better fit for xDH type functionals in general.

Having numerically established the origin of the density and property artifact, let us next analyze it more closely. Returning to the second form of (1), a first order property, $EDHx$⁠, like the dipole moment or force includes direct (Hellman-Feynman) derivative contributions that we can denote as $EDHx$⁠, and non-Hellman-Feynman terms that reflect how the orbital degrees of freedom, θ, change as a function of the perturbation, x. These terms arise only from the non-SCF part, ΔE_DH, of E_DH. In detail, the first order derivative (property) is given by

In the second form, z is the so-called z-vector³¹ which corresponds to the occupied-virtual or response block of the DH density matrix and is the solution to the response equation in terms of the orbital-optimized energy $EDHOOθθ′z+ΔEDHθ=0$ (where $EDHOOθθ′$ is the Hessian of the orbital optimized energy $EDHOO$⁠, with respect to orbital rotation and $ΔEDHθ$ is the rate of change of the non-orbital optimized energy with orbital rotation).

The response density, z, is strictly zero for OO-DH functionals (one of their main benefits). It is evidently small for the tDH functionals tested here, but large for the xDH functionals, XYG3 and XYGJ-OS, that use B3LYP reference orbitals. The B3LYP orbitals give qualitatively correct behavior on their own [though Fig. 2 suggests that they exhibit significant density delocalization and energy lowering at stretched bond lengths relative to the CCSD(2) reference values]. The critical thing, however, is that they are not optimized for the non-PT2 truncated KS part of xDH functionals, leading to a non-Hellman-Feynman contribution [second term of (2)] to both dipoles and forces from the non-PT2 part as well. This term can be large if there is a major mismatch between the orbitals and the non-PT2 part of the functional. This is completely consistent with the numerically identified origin of the unphysical behavior, as presented in Fig. 3.

$z=−EDHOOθθ′−1ΔEDHθ$ indicates that z might be unphysically large either due to a poorly conditioned $EDHOOθθ′$ or a large $ΔEDHθ$⁠. The smallest eigenvalues of $EDHOOθθ′$ however do not appear to be particularly different between methods in Fig. 3, indicating that the unphysicalities in the z stem at least in part from $ΔEDHθ$⁠, enabling a chemical interpretation of the mathematics. At the very stretched geometries where XYG3 shows problems, the B3LYP orbitals are too delocalizing relative to the truncated hybrid functional (which has ≈80% exact exchange vs 20% in B3LYP). The orbital gradient of the XC difference, $ΔEDHθ$⁠, therefore corrects the reference orbitals (and density) to be more localizing. Self-consistent iterations that treat XYG3 as an OO-DH rigorously zero this quantity. $ΔEDHθ$ is also made smaller by choices of reference orbitals that are closer to OO-DH such as by treating XYG3 as a tDH functional. However, in XYG3 used as an xDH with B3LYP orbitals, the critical problem is that this linear response term leads to a substantial overcorrection in the density for R > 2.5 Å, which in turn leads to the density and force artifacts already discussed. Overall, the artifacts identified in this work originate from a breakdown in the linear response approximation to OO-DH under conditions such as highly stretched bonds, where there can be a substantial orbital mismatch between the reference orbitals and the optimal ones. Asymptotically incorrect reference densities could in particular lead to a wrong dissociation limit, as is observed for the case of C–F bond breaking in CH₃F, where partial charges persist even at large separations for XYG3 applied on a B3LYP reference (which itself has residual charge of the opposite polarity). A large non-Hellman-Feynman term may even cause the one-particle density matrix to become non-N-representable in extreme cases.

While the breakdown of XYG3 and XYGJ-OS is interesting, and potentially valuable as a guide to developing DH functionals in the future, it is important to clarify two points. First, XYG3 and XYGJ-OS are generally very good functionals for the study of chemistry around equilibrium molecular configurations, as is consistent with existing studies.^10,11,14,16 That is consistent with the way in which the few parameters in these xDH functionals were trained, and is consistent with the expectation that there will be only a small orbital mismatch under such conditions. Thus our present results therefore are not in themselves a reason to abandon routine use of XYG3 and XYGJ-OS in the equilibrium regime, or, when there is no reason to expect a significant orbital mismatch. Our results are certainly a reason to recommend that future attempts to develop widely applicable DH type functionals (especially xDH) should carefully take into account behavior at non-equilibrium configurations.

A second point to clarify is that our results do not imply that unphysical behavior will always occur at non-equilibrium configurations as it is quite possible to have a small non-Hellman-Feynman term. Indeed, dissociation curves for H₂ and Li₂ appear to be perfectly well described by both XYG3 and XYGJ-OS. Furthermore, the PT2 contribution also appears to act in the opposite direction from the non-Hellman-Feynman KS term and may in fact be successful in ameliorating unphysical behavior in some cases. It is beyond the scope of this letter to develop a quantitative metric that can be employed to determine a priori whether XYG functionals (or indeed any xDH functional) will give physical results or not. However, it is clear from our present results that there is reason to be cautious wherever a large difference between reference densities and xDH truncated KS self-consistent densities might be suspected.

All calculations performed with a development version of the Q-Chem 4.0 package,³² and unrestricted orbitals were employed in each case. Density functional calculations were done with the aug-cc-pCVQZ^33–36 basis (except for CH₃F, for which aug-cc-pCVTZ was used), and all electron CCSD(2) calculations were done with the aug-cc-pCVTZ basis.

See supplementary material for energies and dipole moments for all problematic species.

This research was supported by the Director, Office of Science, Office of Basic Energy Sciences, of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. D.H. was also funded by a Berkeley Fellowship.