Can Popular DFT Approximations and Truncated Coupled Cluster Theory Describe the Potential Energy Surface of the Beryllium Dimer?*

The beryllium dimer (in the X ¹∑⁺ ground state) is a pathologically multireference weakly bound molecule that has eluded electronic structure methods since the 1960s (see ref. [1] for a comprehensive review of the previous literature).^[1–15] The multireference character of this elusive diatomic system stems from the near-degeneracy of the 2s and 2p orbitals of the beryllium atom, which is a highly multireference system on its own right.^[10,15] The weak bond in Be₂ is due to a mixture of dynamic and static correlation effects.^[5] The combination of a pathologically multireference system and a weak binding energy make Be₂ a notoriously challenging problem for electronic structure methods. The small bond dissociation energy (BDE) of ~900 cm⁻¹ means that an error of 100 cm⁻¹ ≈ 1 kJ mol⁻¹ (which is the target accuracy of highly accurate composite methods such as W4 theory)^[16–20] translates to an error ~10 % in the BDE. In order to achieve a more respectable error of ~1 % in the BDE, one has to calculate the BDE to within ~10 cm⁻¹. In the present work, we define benchmark accuracy as errors of ~1 % from the reference value. The significant multireference character means that high-level ab initio methods have to be employed in conjunction with large basis sets in order to achieve this level of accuracy.

Several recent computational studies investigated the PES of Be₂ around the equilibrium bond distance using high levels of theory.^[2,3,6] These studies have shown that very high levels of theory as well as secondary energetic contributions (e.g. core-valence, scalar relativistic, and diagonal Born–Oppenheimer corrections) are needed in order to reproduce the experimental energetic, structural, and spectroscopic parameters.

Benchmarking approximate ab initio and density functional theory (DFT) procedures against high-level ab initio data has become an important general field over the past two decades (see for example refs [21–25] for an overview). In this context, it is of interest to examine the performance of truncated coupled cluster (CC) and DFT methods for the challenging potential energy surface (PES) of the Be₂ dimer. To this end, we calculate the entire PES of Be₂ at the valence full-configuration interaction (FCI) complete basis set (CBS) limit. As there are only four valence electrons, we are able to extrapolate the FCI energy to the CBS limit from the cc-pV5Z and cc-pV6Z basis sets. We use the valence FCI/cc-pV{5,6}Z reference data to evaluate the performance of the following levels of theory for the entire PES of the Be₂ dimer:

1. FCI/cc-pVnZ (n = D, T, Q, 5, and 6)

2. HF/cc-pV{5,6}Z, CCSD/cc-pV{5,6}Z, CCSD(T)/cc-pV{5,6}Z, CCSDT/cc-pV{5,6}Z, and CCSDT(Q)/cc-pV{5,6}Z

3. Conventional DFT functionals from each rung of Jacob’s Ladder in conjunction with the cc-pV5Z basis set

4. Double-hybrid DFT (DHDFT) methods as well as standard and modified Møller–Plesset perturbation theory (MPn) methods (n = 2–5) in conjunction with the cc-pV5Z basis set.

All the ab initio calculations were carried out using the MRCC program suite on the Linux cluster of the Karton group at the University of Western Australia.^[26,27] In all cases, the correlation-consistent basis sets of Dunning and coworkers were used.^[28–30] Several DFT methods were evaluated for their performance in reproducing the FCI/cc-pV{5,6}Z PES. The cc-pV5Z basis set was used in all the DFT and DHDFT calculations. The considered functionals include: (i) the pure generalized gradient approximation (GGA) functionals BLYP,^[31,32] PBE,^[33,34] revPBE,^[35] B97-D3,^[36] and N12^[37]; (ii) the hybrid GGAs B3LYP,^[31,38,39] B3PW91,^[38,40] PBE0,^[41] and BHandHLYP^[42]; (iii) the hybrid meta-GGAs M06,^[43,44] M06-2X,^[43,44] M06-HF,^[45] PW6B95,^[46] and MN15,^[47] (iv) the range-separated hybrids LC-BLYP,^[48] LC-PBE,^[48] LC-wPBE,^[48] CAM-B3LYP,^[49] and wB97X-D,^[50] and (v) the double hybrids^[51] B2-PLYP,^[52] mPW2-PLYP,^[53] B2GP-PLYP,^[54] the spin-component-scaled double hybrids DSD-BLYP,^[55] DSD-PBEP86,^[56] and the parameter-free PBE0-DH functional.^[57] In some cases, empirical D3 dispersion corrections^[58,59] were included using the Becke–Johnson^[60] damping potential as recommended in ref. 36 (denoted by the suffix -D3BJ).

In addition, the performance of Møller–Plesset perturbation theory is evaluated. We consider the MPn methods (n = 2, 3, 4, 5), as well as the MPn.5 methods (n = 2, 3, and 4). The latter are defined as the average of MPn and MPn₊₁.^[61,62] All of the DFT and MPn calculations were performed using the Gaussian 16 program suite.^[63]

For all levels of theory, a 95-point potential energy curve was calculated. The single-point energy calculations were carried out at bond distances r = r_eq ± x, where r_eq is the equilibrium bond distance at the valence FCI/cc-pV{5,6}Z level of theory and x is varied at 0.001 Å intervals between x = −0.01 and +0.01 Å, at 0.01 Å intervals between x = −1.00 and +2.00 Å, at 0.02 Å intervals between x = −1.60 and +4.00 Å, and at 0.5 Å intervals between x = +4.0 and +14.0 Å. The absolute energies at all the levels of theory considered are given in Table S1 in the Supplementary Material.

Vibrational transition frequencies were obtained using the Duo computer program.^[64] In brief, Duo uses the sinc discrete variable representation (DVR) method to find the direct variational solution to the rotational-less diatomic nuclear motion Schrödinger equation.^[65] Natural quintic spline interpolation was used between the data points to provide the full potential energy curve. These calculations used 5001 grid points in the range 4.0–19.0 a.u.; the results provided are converged to within 0.01 cm⁻¹.

Let us begin by examining the basis set convergence of the valence FCI energy. The PESs calculated at the FCI/cc-pVnZ levels of theory (n = D, T, Q, 5, and 6) are presented in Fig. 1. Errors in the equilibrium bond distance (Δr_e) and BDE (ΔD_e) as well as squared correlation coefficients (R²) with the FCI/cc-pV{5,6}Z PES are given in Table 1. Inspection of Fig. 1 reveals that the shape of the FCI/cc-pVDZ PES is fundamentally flawed. It exhibits two local minima, a deeper one at r_e = 4.855 a.u. and a shallow one at approximately r_e = 8 a.u. Previous studies using more approximate correlated ab initio methods have found that the shallow minimum at large distance is an artefact of the small basis set.^[5,12,13] The cc-pVTZ basis set results in the correct shape of the PES; however, the PES is too narrow and the binding energy is underestimated by as much as 231 cm⁻¹. The shape of the FCI/cc-pVQZ PES is closer to that of the FCI/cc-pV{5,6}Z PES; however, it is still visibly narrower and the binding energy is still underestimated by 105 cm⁻¹. It should be pointed out that extrapolating the PES from the cc-pV{T,Q}Z basis set pair, at the same computational cost, results in significantly better performance and reduces the error in the binding energy by an order of magnitude. The cc-pV5Z basis set still results in an appreciable error of 49 cm⁻¹ in the binding energy. Extrapolating the PES from the cc-pV{Q,5}Z basis set pair slightly overcorrects for the deficiencies of the cc-pV5Z basis set and results in a binding energy that is too large by 10 cm⁻¹. The cc-pV6Z basis set still results in an appreciable error in the binding energy of 28 cm⁻¹. In summary, only the cc-pV{T,Q}Z and cc-pV{Q,5}Z extrapolations result in benchmark accuracy.

Fig. 1.  Basis set convergence of the Be₂ FCI potential energy surface. The PESs are calculated at the FCI/cc-pVnZ (n = D, T, Q, 5, and 6) and FCI/cc-pV{n,n+1}Z (n = T, Q, and 5) levels of theory (in a.u. and cm⁻¹).

Table 1.  Basis set convergence of the Be₂ FCI potential energy surface, equilibrium bond distances (r_e), bond dissociation energies (D_e), and vibrational transition frequencies
The reference values are calculated at the FCI/cc-pV{5,6}Z level of theory. Bond distances are in a.u. and bond energies and frequencies are in cm^−1A

In contrast to the binding energy, which converges exceedingly slowly to the basis set limit, r_e converges much faster. Considering errors smaller than 1 % in the equilibrium bond distance as benchmark accuracy (i.e. errors smaller than 0.0467 a.u.), all the basis sets apart from cc-pVDZ and cc-pVTZ achieve sub-benchmark accuracy. The cc-pV5Z and cc-pV6Z basis set achieve errors that are one order of magnitude smaller than the above target error and the cc-pV{Q,5}Z extrapolation reproduces the cc-pV{5,6}Z r_e exactly.

It is also of interest to examine the basis set convergence of the vibrational transition frequencies (v = 1–5) obtained from the FCI/cc-pVnZ potential energy curves (n = D, T, Q, 5, and 6). These results are presented in Table 1. For all the basis sets considered, the errors with respect to the FCI/cc-pV{5,6}Z vibrational transition frequencies increase in the order Δv₁ < Δv₂ < Δv₃ < Δv₄ < Δv₅. The largest errors (obtained for v₅) are −399.1 (cc-pVDZ), −74.7 (cc-pVTZ), −23.4 (cc-pVQZ), −9.6 (cc-pV5Z), and −5.5 (cc-pV6Z) cm⁻¹. Table S2 in the Supplementary Material lists the relative errors in the vibrational transition frequencies. Inspection of these results reveals that the relative errors are fairly constant across all the vibrational transition frequencies (v₁–v₅). Namely, they are ~81 % (cc-pVDZ), ~12 % (cc-pVTZ), ~3 % (cc-pVQZ), ~1.5 % (cc-pV5Z), and ~1.0 % (cc-pV6Z) cm⁻¹. Overall, the results presented in Fig. 1 and Table 1 demonstrate that the FCI PES converges smoothly, albeit slowly, with the basis set size. The cc-pV{T,Q}Z and cc-pV{Q,5}Z basis set extrapolations provide significantly better performance than the cc-pVQZ and cc-pV5Z basis sets, respectively, at the same computational cost. The FCI/cc-pV{T,Q}Z PES provides a good compromise between accuracy and computational cost and results in errors in the BDE and vibrational transition frequencies that are smaller than or equal to ~10 cm⁻¹.

It is of interest to examine the performance of truncated CC theory in reproducing the entire PES for the beryllium dimer. Fig. 2 compares the HF, CCSD, CCSD(T), CCSDT, CCSDT(Q), and FCI PESs. Errors in the bond distance, BDE, and vibrational transition frequencies with respect to the FCI/cc-pV{5,6}Z values, as well as squared correlation coefficients (R²) with the FCI/cc-pV{5,6}Z PES are given in Table 2. In all cases, the energies are extrapolated to the complete basis set limit from the cc-pV5Z and cc-pV6Z basis sets. As expected,^[14] the Hartree–Fock PES is purely repulsive over the entire PES. However, we note in passing that adding the original D3 dispersion correction^[58] to the HF/cc-pV{5,6}Z energies results in a minimum of 312 cm⁻¹ at 7.2 a.u. (Fig. S1 in the Supplementary Material). Replacing the zero-damping function in the original D3 procedure with the finite Becke–Johnson damping function leads to a more attractive PES, with a minimum of 2055 cm⁻¹ at 5.1 a.u. (see Supplementary Material).^[60] The valence CCSD/cc-pV{5,6}Z PES is strongly repulsive in the bonding region;^[66] however, it exhibits a shallow minimum of 60.4 cm⁻¹ at approximately 8.3 a.u. This is consistent with previous results obtained with smaller basis sets (see for example refs [8] and [14]). The CCSD(T)/cc-pV{5,6}Z level of theory reproduces the correct shape of the PES. The predicted bond distance is longer by 0.039 a.u. relative to the FCI/cc-pV{5,6}Z result; however, the PES is much too narrow and the binding energy is underestimated by as much as 219.5 cm⁻¹. Consideration of higher-order connected triple excitations in the CCSDT/cc-pV{5,6}Z PES leads to an improvement; however, there is still a noticeable difference between the CCSDT and FCI PESs (with R² = 0.9932) and the binding energy is still underestimated by 82.4 cm⁻¹. The non-iterative connected quadruple excitations have a significant effect on the PES and the CCSDT(Q)/cc-pV{5,6}Z PES reproduces the FCI/cc-pV{5,6}Z curve almost perfectly (R² = 0.9998), except in the vicinity of the equilibrium distance where there is still a visible difference between the curves (Fig. 2). In particular, the binding energy is underestimated by 16.6 cm⁻¹.

Fig. 2.  Potential energy surfaces for the Be₂ dimer calculated at the HF/cc-pV{5,6}Z, CCSD/cc-pV{5,6}Z, CCSD(T)/cc-pV{5,6}Z, CCSDT/cc-pV{5,6}Z, CCSDT(Q)/cc-pV{5,6}Z, and FCI/cc-pV{5,6}Z levels of theory (in a.u. and cm⁻¹).

Table 2.  Evaluation of truncated coupled cluster methods for the shape of the potential energy surface of Be₂
The tabulated values are errors in the equilibrium bond distances (Δr_e), bond dissociation energies (ΔD_e), and vibrational transition frequencies (Δv_n). The reference values are calculated at the FCI/cc-pV{5,6}Z level of theory. Bond distances are in a.u. and bond energies and frequencies are in cm^−1A

For all the truncated CC methods, the errors with respect to the FCI/cc-pV{5,6}Z vibrational transition frequencies increase in the order Δv₁ < Δv₂ < Δv₃ < Δv₄ < Δv₅. Thus, it is informative to look at the largest errors obtained for v₅: they are −459.5 (HF/CBS), −426.3 (CCSD/CBS), −92.4 (CCSD(T)/CBS), −35.6 (CCSDT/CBS), and −6.8 (CCSDT(Q)/CBS) cm⁻¹. However, the relative errors in the vibrational transition frequencies (given in Table S3 in the Supplementary Material) are fairly constant across all the vibrational transition frequencies (v₁–v₅). Namely, they are ~98 % (HF/CBS), ~89 % (CCSD/CBS), ~14 % (CCSD(T)/CBS), ~6 % (CCSDT/CBS), and ~1 % (CCSDT(Q)/CBS) cm⁻¹. Thus, it is evident, both from the perspective of the absolute errors and the relative errors, that CCSDT(Q) is the only method that achieves benchmark accuracy. We note, however, that the equilibrium distance has already converged to within 1 % with the CCSD(T) and CCSDT methods (namely, these methods attain errors of 0.039 and 0.027 a.u., respectively).

Fig. 3 shows the PESs of a few GGA functionals (Table S4 in the Supplementary Material lists errors in r_e, D_e, and v_n). All the GGA functionals result in a PES that is much too wide and severely overbind the beryllium dimer. In particular, BLYP and PBE result in BDEs of 2163.0 and 3449.2 cm⁻¹, respectively. revPBE results in a PES that is almost identical to that of PBE. As expected, addition of dispersion corrections increases the binding energies further and results in significantly larger BDEs (Fig. 3). In the remainder of this subsection, we will not consider dispersion corrections for DFT functionals that overbind the beryllium dimer. We have also considered the local non-separable gradient approximation N12 functional, which results in a BDE sandwiched between those of the BLYP and PBE functionals. Interestingly, both the BLYP and N12 functionals exhibit a shallow maximum at long distances. For BLYP, the height of this maximum is 47.5 cm⁻¹ around r_e = 9.6 a.u. and for N12 it is 136.7 cm⁻¹ around r_e = 8.9 a.u.

Fig. 3.  Potential energy surfaces for the Be₂ dimer calculated with several GGA functionals in conjunction with the cc-pV5Z basis set (in a.u. and cm⁻¹).

Despite the fact that the GGA functionals in Fig. 3 give a very poor description of the FCI/cc-pV{5,6}Z PES, their errors in the equilibrium bond distances are not very large. In particular, BLYP and PBE underestimate r_e by ~0.08 a.u. and N12 overestimates it by ~0.03 a.u.

Hybrid GGA functionals are expected to perform better than GGAs because the purely repulsive Hartree–Fock potential compensates for the severe overbinding of the GGAs. Fig. 4 gives the PESs for several HGGAs. Inspection of this figure reveals that HGGAs with low percentages of HF exchange still severely overbind the Be₂ dimer, albeit to a lesser extent than the GGAs. The popular B3LYP functional (20 % HF exchange) predicts a binding energy of 1492.1 cm⁻¹. This represents a significant improvement over BLYP, which predicts a BDE of 2163.0 cm⁻¹. Replacing the LYP correlation functionals with PW91 results in a significantly wider PES and a binding energy of 2085.9 cm⁻¹. In contrast, BH&HLYP with 50 % exact exchange severely underbinds the dimer, with a binding energy of 374.8 cm⁻¹. Inclusion of the D3BJ dispersion correction overcorrects for this deficiency and leads to a binding energy of 1288.2 cm⁻¹. The PBE0 functional (25 % HF exchange) results in a very wide PES and a binding energy of 2307.7 cm⁻¹. This represents a significant improvement over PBE, which predicts a binding energy of 3449.2 cm⁻¹.

Fig. 4.  Potential energy surfaces for the Be₂ dimer calculated with several hybrid GGA functionals in conjunction with the cc-pV5Z basis set (in a.u. and cm⁻¹).

Fig. 5 depicts the PESs obtained with several hybrid-meta GGA from the Truhlar group. It is instructive to compare the PESs obtained with M06 (27 % exact exchange), M06–2X (54 % exact exchange), and M06-HF (100 % exact exchange). M06 severally overbinds the Be₂ dimer and predicts a binding energy of 1822.8 cm⁻¹. M06-HF underbinds the dimer with a binding energy of 760.7 cm⁻¹. However, M06-2X predicts a binding energy of 928.8 cm⁻¹, which deviates from the FCI values by only 55.9 cm⁻¹. In terms of predicting the binding energy, M06-2X shows the best performance of all the conventional DFT functionals (i.e. excluding the DHDFT methods). Nevertheless, it should be pointed out that M06-2X does not give a good representation of the shape of the FCI PES as demonstrated, for example, from a squared correlation coefficient of R² = 0.8783 with the FCI/cc-pV{5,6}Z PES.

Fig. 5.  Potential energy surfaces for the Be₂ dimer calculated with several hybrid-meta-GGA functionals in conjunction with the cc-pV5Z basis set (in a.u. and cm⁻¹).

Nearly all of the GGA, HGGA, and HMGGA functionals considered so far tend to overbind the Be₂ dimer (the only exceptions being BH&HLYP and M06-HF). In contrast, the considered range-separated functionals tend to systematically underbind the dimer. These results are presented in Fig. 6. LC-BLYP gives a very shallow PES with a minimum of 95.8 cm⁻¹ at ~7 a.u. LC-PBE better represents the FCI PES, but still shows poor performance with a wide PES and a binding energy of 677.5 cm⁻¹. These results should be compared with the results for BLYP and PBE, which severely overbind the Be₂ dimer (Fig. 3). Similarly, CAM-B3LYP results in a BDE of 418.8 cm⁻¹, in contrast to B3LYP, which predicts a binding energy of 1492.1 cm⁻¹ (Fig. 4). Finally, we note that wB97X-D gives the best performance, with a BDE of 735.9 cm⁻¹, which is lower than the FCI BDE by 137.0 cm⁻¹.

Fig. 6.  Potential energy surfaces for the Be₂ dimer calculated with several range-separated DFT functionals in conjunction with the cc-pV5Z basis set (in a.u. and cm⁻¹).

Fig. 7 gives the PESs obtained with several double-hybrid DFT functionals. With the exception of three functionals (PBE0-DH, B2GP-PLYP, and DSD-BLYP), all the DHDFT functionals considered show surprisingly good performance in both reproducing the shape and binding energy of the FCI/CBS PES. Let us begin with the poor performers. Both DSD-BLYP and B2GP-PLYP underbind the Be₂ dimer; for example, B2GP-PLYP underbinds it by 201.9 cm⁻¹ with a BDE of 670.9 cm⁻¹. Addition of the D3BJ dispersion correction overcorrects for this deficiency and results in a very wide PES and a BDE of 1207.4 cm⁻¹. PBE0-DH gives an even wider PES with a binding energy of 1626.3 cm⁻¹. All of the other DHDFT functionals reproduce the FCI/CBS PES very well. The older-generation B2-PLYP and mPW2-PLYP predict binding energies that are too low by 46.7 and 40.3 cm⁻¹, respectively. It is also worth pointing out that mPW2-PLYP predicts the FCI/CBS vibrational transition frequencies to within 10 cm⁻¹ (Table 3). Quite remarkably, the more recently developed spin-component-scaled DHDFT DSD-PBEP86 functional reproduces the FCI/CBS PES spot on with a binding energy of 872.3 cm⁻¹, less than 1 cm⁻¹ from the FCI/CBS binding energy.

Fig. 7.  Potential energy surfaces for the Be₂ dimer calculated with several double-hybrid DFT functionals in conjunction with the cc-pV5Z basis set (in a.u. and cm⁻¹).

Table 3.  Evaluation of DHDFT methods for the shape of the potential energy surface of Be₂
The tabulated values are errors in the equilibrium bond distances (Δr_e), bond dissociation energies (ΔD_e), and vibrational transition frequencies (Δv_n). The reference values are calculated at the FCI/cc-pV{5,6}Z level of theory. Bond distances are in a.u. and bond energies and frequencies are in cm^−1A

Finally, it is of interest to compare the performance of DHDFT with MPn theory. These results are presented in Fig. 8. Inspection of Fig. 8 reveals that the MPn series up to MP4 converges monotonically towards the FCI solution in the order MP2 → MP2.5 → MP3 → MP3.5 → MP4. For example, the following binding energies are obtained: 440.3 (MP2), 524.4 (MP2.5), 614.4 (MP3), 693.5 (MP3.5), and 783.3 (MP4) cm⁻¹. However, the MP4.5 and MP5 PESs provide no further improvement and are practically identical to the MP4 PES (e.g. the predict binding energies of 786.2 (MP4.5) and 790.7 (MP5) cm⁻¹). In light of these results, it is clear that DHDFT provides a much better description of the PES compared with the MPn methods.

Fig. 8.  Potential energy surfaces for the Be₂ dimer calculated with the MPn methods in conjunction with the cc-pV5Z basis set (in a.u. and cm⁻¹).

We have obtained the entire PES of the beryllium dimer at the FCI/cc-pV{5,6}Z level of theory. We use this reference data to evaluate (i) the basis set convergence of the valence FCI method; (ii) the performance of truncated CC methods at the infinite basis-set limit; (iii) the performance of DFT methods across the rungs of Jacob’s Ladder; and (iv) the performance of standard and modified MPn methods. Our main findings can be summarized as follows:

• The FCI/cc-pV{T,Q}Z basis set extrapolation reproduces the shape of the FCI/cc-pV{5,6}Z PES as well as the binding energy and vibrational transition frequencies to within ~10 cm⁻¹.

• Of the truncated CC methods, CCSDT(Q)/cc-pV{5,6}Z provides a very good representation of the FCI/cc-pV{5,6}Z PES.

• GGA functionals, as well as hybrid-GGA and hybrid-meta-GGA functionals with low percentages of exact exchange tend to severely overbind the Be₂ dimer over the entire PES.

• Range-separated functionals tend to underbind the Be₂ dimer.

• DHDFT functionals show exceptionally good performance relative to their computational cost.

• Møller–Plesset perturbation theory converges smoothly up to fourth order; however, fifth-order corrections have a minor effect on the PES.

Absolute energies at all the considered levels of theory (Table S1); relative errors in the vibrational transition frequencies for the ab initio methods (Tables S2 and S3); errors in the bond distance, BDE, and vibrational transition frequencies for the DFT and MPn methods with respect to the FCI/cc-pV{5,6}Z values, as well as squared correlation coefficients (R²) with the FCI/cc-pV{5,6}Z PES (Table S4); and PES calculated at the HF-D3/cc-pV{5,6}Z and HF-D3BJ/cc-pV{5,6}Z levels of theory (Fig. S1) are available on the Journal’s website.

The authors declare no conflicts of interest.