We report the variation of the binding energy of the Formic Acid Dimer with the size of the basis set at the Coupled Cluster with iterative Singles, Doubles and perturbatively connected Triple replacements [CCSD(T)] level of theory, estimate the Complete Basis Set (CBS) limit, and examine the validity of the Basis Set Superposition Error (BSSE)-correction for this quantity that was previously challenged by Kalescky, Kraka, and Cremer (KKC) [J. Chem. Phys. 140, 084315 (2014)]. Our results indicate that the BSSE correction, including terms that account for the substantial geometry change of the monomers due to the formation of two strong hydrogen bonds in the dimer, is indeed valid for obtaining accurate estimates for the binding energy of this system as it exhibits the expected decrease with increasing basis set size. We attribute the discrepancy between our current results and those of KKC to their use of a valence basis set in conjunction with the correlation of all electrons (i.e., including the 1s of C and O). We further show that the use of a core-valence set in conjunction with all electron correlation converges faster to the CBS limit as the BSSE correction is less than half than the valence electron/valence basis set case. The uncorrected and BSSE-corrected binding energies were found to produce the same (within 0.1 kcal/mol) CBS limits. We obtain CCSD(T)/CBS best estimates for De = − 16.1 ± 0.1 kcal/mol and for D0 = − 14.3 ± 0.1 kcal/mol, the later in excellent agreement with the experimental value of −14.22 ± 0.12 kcal/mol.
I. INTRODUCTION
The quantitative description of intermolecular interactions represents a challenging undertaking in modern electronic structure theory.1 Over the years, a path that combines advances in methods treating electron correlation combined with approaches that focus on systematically saturating the orbital basis set towards a complete limit has been proven quite successful in providing quantitative descriptions of a wide range of inter- and intra-molecular bonding interactions.2 This path amounts to the use of highly correlated methods, such as Configuration Interaction (CI)3–5 or Coupled Cluster (CC)6–11 treatments based on either single or multireference zero-order type wavefunctions with a simultaneous hierarchical expansion of the orbital basis set.
Several aspects of how the computed properties (structures, energetics, electric, and magnetic properties) vary upon increasing the basis set towards a well defined Complete Basis Set (CBS) limit have been previously reported in voluminous studies since the early 1990s12–23 using the correlation consistent basis sets of Dunning and co-workers, (aug)cc-pVnZ, n = D, T, Q, 5, 6.24–32 These basis sets are designed to systematically expand both the radial and angular parts of the respective atomic wavefunctions and are most appropriate to be used to attain a CBS limit. In this paper, we will focus on the following aspects that arise when following such a path for computing accurate intermolecular binding energies:
the convergence of the binding energy with basis set,
the quantitative account of the Basis Set Superposition Error (BSSE),
the existence of a common limit for the former two when the basis set is saturated, and
issues related to linear dependencies in the basis set as it approaches saturation and how these affect the behavior of the former three aspects.
The concept of the BSSE appeared in the literature about half a century ago.33–35 According to this concept, the basis functions on one atom “gives an extension of the basis set” of another atom.33 Therefore, a choice of balanced basis sets must always be made among all atoms of a molecule. To remedy this problem, the function counter-poise (fCP) correction is widely applied as first proposed by Boys and Bernardi,34 and extended to account for the geometrical distortions arising from the interaction in polyatomic molecular systems by Xantheas.36 Other aspects of this correction to the binding energy have been previously discussed.37–44
The general consensus from the plethora of published literature is that there exist clear limits for the binding energy evaluated with a particular theoretical method upon saturating the basis set either by applying or not the BSSE correction and that those two limits lie quite close to each other. After all, the concept of the BSSE correction is related to the oftentimes incompleteness of the basis set and as such this effect, albeit not proven via a strict mathematical approach, should decrease upon increasing the basis set and ultimately vanish at the CBS. As a consequence, a more or less “common” (depending on the distance of the largest basis set used from the CBS) limit between the uncorrected and BSSE-corrected binding energies is found when their convergence with increasing basis set is considered. The fourth aspect, namely, the one of linear dependencies especially as the basis set approaches CBS, has been rarely considered in the literature before, however, as will be shown here, it has a dramatic effect on all previous four aspects.
This traditional picture and anticipated behavior with basis set size have been recently challenged by Kalescky, Kraka, and Cremer (KKC) for the binding energy of the Formic Acid Dimer (FAD).45 These authors even confronted the validity of the BSSE correction for this system based on some earlier work by Baerends and co-workers.46 KKC reported that the BSSE correction to the binding energy of FAD with the aug-cc-pVTZ basis set is even slightly larger (1.45 kcal/mol) than the one with the smaller aug-cc-pVDZ set (1.43 kcal/mol) and, based on this finding and the earlier work of Baerends and co-workers46 albeit for a chemically different system (Be2), they claimed that “the BSSE corrections obtained at the Coupled Cluster with iterative Singles, Doubles and perturbatively connected Triple replacements [CCSD(T)] level of theory are questionable,” and that no BSSE corrections should be applied.45 They reported binding energies at the CCSD(T) from a single reference Hartree Fock (HF) wavefunction with basis sets as large as pentuple zeta (aug-cc-pV5Z), without discussing any issues related to linear dependencies that can undoubtedly arise for this very large basis set with very diffuse high angular momentum basis functions. The use of a large basis set, such as aug-cc-pV5Z, is often problematic because of the linear dependency issues that cause convergence problems during the coupled-cluster iterations. Some codes like NWChem47 can project out a specific number of linear dependencies based on a threshold for the overlap of the symmetry adapted basis functions. This, however, reduces the molecular orbital space resulting in a potential mismatch of the energies. Recently, we studied the benzene dimer using basis sets up to cc-pV5Z and aug-cc-pVQZ at the CCSD(T) level.48 In that study, we had to remove 40 linear dependencies at the CCSD(T)/aug-cc-pVQZ but none at the CCSD(T)/cc-pVQZ and only 4 at the CCSD(T)/cc-pV5Z levels, whereas we were not able to perform any CCSD(T)/aug-cc-pV5Z calculations. The large number of linear dependencies in that study resulted in a discrepancy of ∼0.15 kcal/mol (or 5.5%) in the CBS limit of the binding energy between the plain and augmented series (2.65 vs. 2.80 kcal/mol).48 Conclusively, special attention should been given to the basis set linear dependencies issue. Finally, KKC relied on usual extrapolation schemes based on a monotonic dependence of the harmonic frequencies with basis set to obtain CBS-quality harmonic frequencies and by correcting for anharmonicities at the second order vibrational perturbation theory level (VPT2) they arrived at a best estimate of De = − 15.55 kcal/mol and a D0 = − 14.32 kcal/mol, which compared favorably with the experimentally obtained value of −14.22 ± 0.22 kcal/mol from Suhm and co-workers.49
The previous study by KKC is, to the best of our knowledge, the only one that challenges the traditional behavior expected for the variation of the binding energies for hydrogen bonded systems and, if it were true, it would have a profound effect on the generally accepted consensus in the community. To resolve this issue, we carried out extensive CCSD(T) calculations employing both the plain and augmented correlation consistent basis sets with up to five-zeta quality for the formic acid monomer (FAM) and FAD. Linear dependencies needed to be removed only during calculations with the aug-cc-pV5Z basis set. The BSSE corrections were applied and the CBS limit was estimated by extrapolating both the BSSE-uncorrected and BSSE-free binding energies. Our recent results (vide infra) are in complete disagreement with the earlier conclusions of KKC, showing in contrast a normal behavior with all cases converging asymptotically to the CBS limit. For example, our BSSE corrections at the CCSD(T) level with the aug-cc-pVDZ and aug-cc-pVTZ basis sets are 2.58 kcal/mol and 1.46 kcal/mol, respectively, compared to 1.43 kcal/mol and 1.45 kcal/mol reported by KKC. A rather complete account of the literature on the FAD is given by KKC45 and will not be repeated here.
In Sec. II, we provide a detailed account of our computational approach. In Sec. III, we present our numerical results, compare them with the ones in the literature including those previously reported by KKC, and identify the source of the discrepancy between our results and those reported earlier by KKC. Finally, in Sec. IV, we summarize our findings.
II. COMPUTATIONAL DETAILS
The CCSD(T) approach based on a restricted Hartree-Fock reference wavefunction was applied to obtain optimal geometries and energies. Both the cc-pVnZ (abbreviated VnZ) and aug-cc-pVnZ (abbreviated AVnZ) families of basis sets with n = D, T, Q, 5 were employed for all atoms.24,25 We obtained optimized geometries with both families of basis sets for n = D, T, Q, while the VQZ and AVQZ geometries were used for the reported V5Z and AV5Z energies, respectively. The optimal geometries of the FAM and FAD have Cs and C2h symmetries, respectively. Due to symmetry, which results in the two monomers having identical geometries in the dimer, the BSSE uncorrected binding energy, De, is
where we follow the notation proposed earlier by Xantheas,36 in which parentheses identify molecular systems, subscripts denote optimized geometries of moieties, and superscripts are used to indicate the orbital basis sets used. The first term in the right hand side of Eq. (1) corresponds to the energy of the complex (AB) at the optimal geometry of AB (subscript AB) using the basis sets α and β on monomers A and B, respectively. The second term of the right hand side is twice the energy of the monomer (A) at its optimal geometry (subscript A) using basis functions α only on that fragment A. Using the same notation, the counterpoise-corrected binding energy, , is36
The second term in the right hand side is twice the energy of fragment (A) in the dimer AB geometry evaluated with the full basis of monomers A and B (the union of basis sets α and β). The quantity in the square brackets is the relaxation or deformation energy (Edef), namely, the difference of the energy of monomer A (or B) between its optimal geometry in isolation and its distorted geometry at the dimer’s optimal geometry. As shown earlier,36 Eqs. (1) and (2) converge to the same limit as the orbital basis approaches a complete set.
The CBS limit of the binding energy was obtained by fitting the De values for the VnZ (BSSE-corrected and uncorrected) and the AVnZ (BSSE-corrected only) series to the exponential expression12,21
The AVnZ uncorrected De values do not follow an exponential decay and thus were fitted using the expression50–52
The four different sets of values (VnZ and AVnZ, both BSSE-corrected and uncorrected) were fitted separately and four different CBS limits were obtained ranging within 0.1 kcal/mol from each other (vide infra). Alternatively, we fitted simultaneously the BSSE-corrected and uncorrected values for the VnZ series. Similarly, a “common” CBS limit for the AVnZ BSSE corrected and uncorrected values was also obtained.
Linear dependency issues emerged only for the FAD at the CCSD(T)/AV5Z level; in that case, a total of 40 linear dependencies had to be removed. This large number renders the calculated energy questionable and therefore no binding energy at this level is presently reported. However, we report the CCSD(T)/AV5Z//CCSD(T)/AVQZ energy of FAM.
Our best estimate for the anharmonic frequencies of the FAM and FAD is derived by combining the CCSD(T)/AVQZ harmonic vibrational frequencies (ωi) with the MP2/AVDZ anharmonicities (νi − ωi) obtained at the VPT2 according to
The CCSD(T)/AVQZ harmonic vibrational frequencies were calculated numerically using the Z-matrix formulation as implemented recently by Miliordos and Xantheas.53 The employed Z-matrix coordinates and the resulting normal modes are given in the supplementary material.54 The zero point energy estimated via VPT2, (ZPE)V PT2, is55,56
where νi are the fundamental (anharmonic) vibrational fre-quencies, (ZPE)harm and (ZPE)fund denote the harmonic and fundamental zero-point energies whereas χii are the diagonal anharmonicities. All CCSD(T) calculations were carried out using the NWChem47 and MOLPRO suites of codes,57 while the VPT2 anharmonicities at the MP2/AVDZ level were estimated using the Gaussian09 program.58
III. RESULTS AND DISCUSSION
Table I lists the optimal geometrical parameters defined in Figure 1 at the CCSD(T)/VnZ and CCSD(T)/AVnZ levels with n = D, T, Q for both the FAM and FAD molecular systems. The optimal bond lengths and bond angles are within 0.001 Å and 0.5° between the VQZ and AVQZ basis sets confirming that the two families (VnZ and AVnZ) lead to the same CBS limit geometries. Going from VTZ to VQZ or from AVTZ to AVQZ, the bond lengths change by a few mili-Å and the bond angles by ∼0.5°. Therefore, the VQZ and AVQZ geometries can be considered converged and are used for single point calculations with the V5Z and AV5Z basis sets. As mentioned earlier, we were able to perform CCSD(T)/AV5Z energy calculations with zero linear dependencies only for the FAM, while 40 linear dependencies had to be removed to obtain a CCSD(T)/AV5Z//CCSD(T)/AVQZ energy of −379.192 452 a.u. Interestingly, KKC did report a De value at the CCSD(T)/AV5Z level of theory;45 however, they did not provide either the absolute energies for comparison or discussed any issues related to linear dependencies with that large basis set.
Basis set . | RCH . | RCO . | ROH . | R′CO . | φHCO . | φCOH . | φ′HCO . | Rhb . | φhb . |
---|---|---|---|---|---|---|---|---|---|
HCOOH | |||||||||
VDZ | 1.1105 | 1.3526 | 0.9756 | 1.2092 | 109.32 | 105.06 | 125.52 | ||
VTZ | 1.0944 | 1.3472 | 0.9686 | 1.2023 | 109.75 | 106.14 | 125.22 | ||
VQZ | 1.0938 | 1.3437 | 0.9669 | 1.1995 | 109.92 | 106.55 | 125.15 | ||
AVDZ | 1.1063 | 1.3614 | 0.9747 | 1.2144 | 109.78 | 106.47 | 125.19 | ||
AVTZ | 1.0948 | 1.3480 | 0.9701 | 1.2040 | 109.89 | 106.60 | 125.13 | ||
AVQZ | 1.0940 | 1.3442 | 0.9675 | 1.2004 | 109.99 | 106.73 | 125.12 | ||
(HCOOH)2 | |||||||||
VDZ | 1.1092 | 1.3218 | 0.9977 | 1.2267 | 111.30 | 108.36 | 121.97 | 1.6961 | 179.42 |
VTZ | 1.0934 | 1.3143 | 0.9950 | 1.2207 | 111.67 | 109.45 | 121.90 | 1.6701 | 179.89 |
VQZ | 1.0929 | 1.3117 | 0.9925 | 1.2175 | 111.76 | 109.65 | 122.00 | 1.6799 | 179.35 |
AVDZ | 1.1048 | 1.3295 | 0.9978 | 1.2321 | 111.59 | 109.27 | 122.18 | 1.7024 | 179.99 |
AVTZ | 1.0939 | 1.3159 | 0.9955 | 1.2221 | 111.77 | 109.54 | 122.00 | 1.6787 | 179.28 |
AVQZ | 1.0932 | 1.3127 | 0.9934 | 1.2186 | 111.81 | 109.72 | 122.00 | 1.6822 | 178.75 |
Basis set . | RCH . | RCO . | ROH . | R′CO . | φHCO . | φCOH . | φ′HCO . | Rhb . | φhb . |
---|---|---|---|---|---|---|---|---|---|
HCOOH | |||||||||
VDZ | 1.1105 | 1.3526 | 0.9756 | 1.2092 | 109.32 | 105.06 | 125.52 | ||
VTZ | 1.0944 | 1.3472 | 0.9686 | 1.2023 | 109.75 | 106.14 | 125.22 | ||
VQZ | 1.0938 | 1.3437 | 0.9669 | 1.1995 | 109.92 | 106.55 | 125.15 | ||
AVDZ | 1.1063 | 1.3614 | 0.9747 | 1.2144 | 109.78 | 106.47 | 125.19 | ||
AVTZ | 1.0948 | 1.3480 | 0.9701 | 1.2040 | 109.89 | 106.60 | 125.13 | ||
AVQZ | 1.0940 | 1.3442 | 0.9675 | 1.2004 | 109.99 | 106.73 | 125.12 | ||
(HCOOH)2 | |||||||||
VDZ | 1.1092 | 1.3218 | 0.9977 | 1.2267 | 111.30 | 108.36 | 121.97 | 1.6961 | 179.42 |
VTZ | 1.0934 | 1.3143 | 0.9950 | 1.2207 | 111.67 | 109.45 | 121.90 | 1.6701 | 179.89 |
VQZ | 1.0929 | 1.3117 | 0.9925 | 1.2175 | 111.76 | 109.65 | 122.00 | 1.6799 | 179.35 |
AVDZ | 1.1048 | 1.3295 | 0.9978 | 1.2321 | 111.59 | 109.27 | 122.18 | 1.7024 | 179.99 |
AVTZ | 1.0939 | 1.3159 | 0.9955 | 1.2221 | 111.77 | 109.54 | 122.00 | 1.6787 | 179.28 |
AVQZ | 1.0932 | 1.3127 | 0.9934 | 1.2186 | 111.81 | 109.72 | 122.00 | 1.6822 | 178.75 |
The monomer geometry is largely distorted from its optimal gas phase geometry upon complexation to form the dimer, and the problems arising from its omission in estimating accurate dimerization energies have been pointed out in earlier studies.59 Except from the intact CH bond, the rest of the bonds elongates by ∼0.02-0.03 Å as a result of the two strong hydrogen bonded interactions. Specifically, the OH length increases by 0.026 Å with both the VQZ and AVQZ basis sets. Note that the three hydrogen-bonded atoms are practically aligned (φhb ≈ 180°). Our CCSD(T)/AVQZ optimal geometry is in very good agreement with that of KKC,45 except for a discrepancy of ∼0.01 Å for the OH⋯O distance, viz., Rhb = 1.669 Å45 vs. our 1.682 (AVQZ) and 1.680 Å (VQZ).
Table II collects the energies at the FAM and FAD equilibrium geometries, as well as those required for the counter-poise correction of the BSSE. An interesting observation is the large deformation energy being nearly independent of the size of the basis set with a magnitude of 1.24 ± 0.13 kcal/mol. The corresponding quantity for the water and parallel-displaced benzene dimers is less than 0.1 and 0.01 kcal/mol, respectively,36,48 the latter due to the weak π − π interaction. Even larger deformation energies of 3.98 kcal/mol (MP2/AVDZ) and 5.08 kcal/mol (MP2/AV5Z) have been previously reported for the F−(H2O) cluster due to the large binding energy that greatly distorts the water geometry in the dimer.36
. | VDZ . | VTZ . | VQZ . | V5Za . | CBS . |
---|---|---|---|---|---|
b | −378.647 257 | −379.030 858 | −379.149 593 | −379.188 176 | |
c | −189.308 960 | −189.501 934 | −189.561 659 | −189.581 158 | |
d | −189.306 930 | −189.499 745 | −189.559 629 | −189.579 181 | |
e | −189.312 507 | −189.501 959 | −189.560 505 | −189.579 501 | |
Edeff | 1.27 | 1.37 | 1.27 | 1.24 | |
Deg | −18.41 | −16.94 | −16.49 | −16.23 | −16.15h |
BSSEi | 7.00 | 2.78 | 1.10 | 0.40 | |
j | −11.41 | −14.16 | −15.39 | −15.83 | −16.21h |
(−16.18)k | |||||
AVDZ | AVTZ | AVQZ | AV5Za | CBS | |
b | −378.738 577 | −379.062 282 | −379.161 658 | −379.192 452l | |
c | −189.356 188 | −189.517 766 | −189.567 693 | −189.583 512 | |
d | −189.354 416 | −189.515 798 | −189.565 707 | −189.581 506 | |
e | −189.356 474 | −189.516 960 | −189.566 219 | … | |
Edeff | 1.11 | 1.23 | 1.25 | 1.26 | |
Deg | −16.44 | −16.79 | −16.49 | … | −16.02h |
Dem | −16.71 | −18.12 | −17.00 | ||
BSSEi | 2.58 | 1.46 | 0.64 | … | |
j | −13.86 | −15.33 | −15.84 | … | −16.11n |
(−16.07)k |
. | VDZ . | VTZ . | VQZ . | V5Za . | CBS . |
---|---|---|---|---|---|
b | −378.647 257 | −379.030 858 | −379.149 593 | −379.188 176 | |
c | −189.308 960 | −189.501 934 | −189.561 659 | −189.581 158 | |
d | −189.306 930 | −189.499 745 | −189.559 629 | −189.579 181 | |
e | −189.312 507 | −189.501 959 | −189.560 505 | −189.579 501 | |
Edeff | 1.27 | 1.37 | 1.27 | 1.24 | |
Deg | −18.41 | −16.94 | −16.49 | −16.23 | −16.15h |
BSSEi | 7.00 | 2.78 | 1.10 | 0.40 | |
j | −11.41 | −14.16 | −15.39 | −15.83 | −16.21h |
(−16.18)k | |||||
AVDZ | AVTZ | AVQZ | AV5Za | CBS | |
b | −378.738 577 | −379.062 282 | −379.161 658 | −379.192 452l | |
c | −189.356 188 | −189.517 766 | −189.567 693 | −189.583 512 | |
d | −189.354 416 | −189.515 798 | −189.565 707 | −189.581 506 | |
e | −189.356 474 | −189.516 960 | −189.566 219 | … | |
Edeff | 1.11 | 1.23 | 1.25 | 1.26 | |
Deg | −16.44 | −16.79 | −16.49 | … | −16.02h |
Dem | −16.71 | −18.12 | −17.00 | ||
BSSEi | 2.58 | 1.46 | 0.64 | … | |
j | −13.86 | −15.33 | −15.84 | … | −16.11n |
(−16.07)k |
The VQZ and AVQZ optimal geometries for V5Z and AV5Z basis sets were used.
Energy of FAD at its optimal geometry.
Energy of FAM at its optimal geometry.
Energy of FAM at its geometry in the dimer’s optimal geometry with the monomer basis.
Same as d, but with the full dimer basis.
Deformation energy, .
.
De values at the VnZ or AVnZ basis set were fitted to Eq. (3).
.
.
De and were fitted simultaneously.
Forty linear dependencies were removed, see text.
Reference 45.
De values at the AVnZ basis set were fitted to Eq. (4).
The FAD binding energies are also listed in Table II and their variation with respect to cardinal number n of the basis set is shown in Figure 2. Observe that the binding energies (both uncorrected and BSSE-corrected) with the plain VnZ basis set exhibit an exponential decay with n. On the other hand, the BSSE uncorrected values with the AVnZ sets follow a different pattern, starting from −16.44 kcal/mol (AVDZ), decreasing to −16.79 kcal/mol (AVTZ) and then increasing to −16.49 kcal/mol (AVQZ). For this reason, we used the “4-5” inverse power extrapolation formula for this case. The corresponding values reported by KKC were −16.71 (AVDZ), −18.12 (AVTZ), and −17.00 (AVQZ) kcal/mol.45 Their CCSD(T)/AVTZ binding energy differs from ours by 1.3 kcal/mol or ∼8%. Unfortunately, these authors do not provide this information or the absolute energies for a direct comparison. After correcting for BSSE, the exponential decay is recovered for the AVnZ series (see Figure 2). Additionally, unlike KKC, we found a normal decreasing trend for the BSSE with the larger AVnZ basis sets. Contrary to their reported BSSE correction of 1.43 (AVDZ) and 1.45 (AVTZ) kcal/mol—showing an unphysical increase of the BSSE-correction with increasing basis set size—our values of 2.58 (AVDZ), 1.46 (AVTZ), and 0.64 (AVQZ) kcal/mol instead support the expected decrease of this correction with basis set size. As expected, the BSSE-correction with the plain (VnZ) correlation consistent basis sets is much larger than the corresponding one with the augmented (AVnZ) ones, viz., 7.00 (VDZ), 2.78 (VTZ), 1.10 (VQZ), and finally 0.40 (V5Z) kcal/mol as can be seen from Table II.
It is of interest to identify the nature of the discrepancy between our results and those of KKC for both the uncorrected and BSSE-corrected binding energies with the different basis sets. In Table III, we list the magnitude of the BSSE correction when using various combinations of which electrons are correlated (valence only/all electrons) in conjunction with different classes (valence/core-valence) of basis sets. “Valence only” denotes correlation of just the valence electrons whereas the designation “Valence”/“Core-valence” is used for the aug-cc-pVnZ/aug-cc-pCVnZ, n = D, T, Q, basis sets. The possible combinations are Valence only electron correlation/Valence basis sets (V/V), All Electron correlation/Valence basis sets (AE/V), and All Electron correlation/Core-Valence basis sets (AE/CV). The results of KKC are also listed in Table III for comparison. The energies of the various pieces that comprise the BSSE correction, Eqs. (1) and (2), are given in the supplementary material.54 The first line of Table III shows the (V/V) results already listed in Table II that exhibit the expected monotonic decrease of the magnitude of BSSE with increasing basis set. The AE/V results, in contrast, show a non-monotonic behavior with increasing basis set size, following very closely the results reported by KKC. The fact that our results are different (albeit by < 0.2 kcal/mol) from the results of KKC is probably attributed to the fact that we used the (V/V) optimized geometries to perform the (AE/V) calculations and we speculate (as this information is not mentioned in the KKC paper) that they optimized the geometries at the (AE/V) level. Finally, the (AE/CV) results (at the (V/V) geometries) also show a monotonic (decreasing) BSSE correction pattern with basis set that is less than half of the corresponding (V/V) values. These results suggest that the nature of the discrepancy between our results and those of KKC can be attributed to them correlating all electrons using a basis set designed to recover valence correlation. We also want to point out that should KKC have estimated the BSSE with the AVQZ set (0.55 kcal/mol, see Table III), albeit via (AE/V) calculations, they would probably have changed their conclusions about the validity of this correction. Tzeli and Tsekouras have previously warned when using this combination (AE/V) to estimate the magnitude of the BSSE correction.43 Another interesting observation is that the (AE/CV) calculations converge faster to the CBS limit than the (V/V) ones, since the magnitude of the BSSE with the former is smaller (almost half) than the one for the later.
Correlation . | Basis set . | Notation . | n = 2 . | n = 3 . | n = 4 . | Reference . |
---|---|---|---|---|---|---|
Valence only | Valencea,b,c | V/V | 2.58a | 1.46b | 0.64c | This study |
All electrons | Valencea,b,c | AE/V | 1.42a | 1.43b | 0.55c | This study |
All electrons | Core-valenced,e,f | AE/CV | 1.41d | 0.67e | 0.27f | This study |
1.43a | 1.45b | n/a | KKCg |
When fitting separately the four different series of data (AVnZ or VnZ, uncorrected or BSSE-corrected) we obtained the CBS limits for the binding energy in each case. Alternatively, we fitted together the plain and counter-poise corrected values for the two families (plain or augmented) of basis sets. The CBS limits obtained by those different procedures are listed in Table II and they vary within the −16.1 ± 0.1 kcal/mol range. The CCSD(T)/CBS limit of −15.55 kcal/mol previously reported by KKC45 is slightly smaller, but their CCSD(T)-F12a/CBS or CCSD(T)-F12b/CBS values lie within our suggested range. We believe that our −16.1 ± 0.1 kcal/mol range should be considered more accurate than those earlier predictions. We also note that the extrapolation of the AE/CV binding energies produces CBS estimates that are within the error bars of the above value.
We finally attempt to separate the convergence of the HF energy (E[HF]) from the one for the correlation energy (ECorr.), the later defined as ECorr. = E[CCSD(T)] − E[HF]. These energies (a.u.) and the corresponding HF and CCSD(T) binding energies (kcal/mol) are listed in Table IV. The HF contribution to the binding energy is overbound with the smaller basis sets relatively to the HF/CBS limit whereas the correlation contribution to the binding energy is underbound with the smaller basis sets relatively to the ECorr./CBS limit. In Table V, we list the extrapolated CBS limits of the HF, ECorr., and CCSD(T) absolute energies (a.u.) of the FAM and FAD by fitting the n = 2, 3, 4 (and 5 where applicable) energies to an exponential decay function while imposing a common CBS limit for the two families of basis set series. By using the CBS limits of the various extrapolated energies, we obtain values of −11.7 ± 0.2 kcal/mol for the HF/CBS, −4.5 ± 1.0 kcal/mol for the ECorr./CBS, and −16.1 ± 0.9 kcal/mol for the CCSD(T)/CBS limits, respectively. Note that the error bars for De (±0.9 and ±1.0 kcal/mol, last column in Table V), obtained as where σFAM and σFAD are the corresponding errors in the extrapolated monomer and dimer absolute energies listed in Table V, are almost an order of magnitude larger than the one (±0.1 kcal/mol) obtained from the extrapolation of the binding energies of Table III. An alternative estimate of −16.2 ± 1.0 kcal/mol is obtained from the expression De [CCSD(T)] = De [HF] + De[ECorr.] using the absolute energies of FAM and FAD. These two estimates are within 0.1 kcal/mol, albeit with an error bar that is almost a magnitude larger, from the value obtained from extrapolating the energy differences of Table III.
Method . | n = 2 . | n = 3 . | n = 4 . | n = 5 . |
---|---|---|---|---|
(HCOOH)2 | ||||
HF/cc-pVnZ | −377.583 818 | −377.701 133 | −377.731 010 | −377.738 019 |
CCSD(T)/cc-pVnZ | −378.647 257 | −379.030 858 | −379.149 593 | −379.188 176 |
ECorr.a/cc-pVnZ | −1.063 439 | −1.329 725 | −1.418 583 | −1.450 157 |
HF/aug-cc-pVnZ | −377.609 273 | −377.706 406 | −377.732 279 | |
CCSD(T)/aug-cc-pVnZ | −378.738 577 | −379.062 282 | −379.161 658 | |
ECorr.a/aug-cc-pVnZ | −1.129 304 | −1.355 876 | −1.429 379 | |
HCOOH | ||||
HF/cc-pVnZ | −188.779 739 | −188.840 624 | −188.856 020 | −188.859 688 |
CCSD(T)/cc-pVnZ | −189.308 960 | −189.501 934 | −189.561 659 | −189.581 158 |
ECorr.a/cc-pVnZ | −0.529 221 | −0.661 310 | −0.705 639 | −0.721 470 |
HF/aug-cc-pVnZ | −188.794 689 | −188.843 641 | −188.856 789 | |
CCSD(T)/aug-cc-pVnZ | −189.356 188 | −189.517 766 | −189.567 693 | |
ECorr.a/aug-cc-pVnZ | −0.561 498 | −0.674 126 | −0.710 904 | |
De (kcal/mol) | ||||
HF/cc-pVnZ | −15.27 | −12.48 | −11.90 | |
ECorr.a/cc-pVnZ | −3.14 | −4.46 | −4.59 | |
CCSD(T)/cc-pVnZ | −18.41 | −16.94 | −16.49 | |
HF/aug-cc-pVnZ | −12.48 | −12.00 | −11.74 | |
ECorr.a/aug-cc-pVnZ | −3.96 | −4.79 | −4.75 | |
CCSD(T)/aug-cc-pVnZ | −16.44 | −16.79 | −16.49 |
Method . | n = 2 . | n = 3 . | n = 4 . | n = 5 . |
---|---|---|---|---|
(HCOOH)2 | ||||
HF/cc-pVnZ | −377.583 818 | −377.701 133 | −377.731 010 | −377.738 019 |
CCSD(T)/cc-pVnZ | −378.647 257 | −379.030 858 | −379.149 593 | −379.188 176 |
ECorr.a/cc-pVnZ | −1.063 439 | −1.329 725 | −1.418 583 | −1.450 157 |
HF/aug-cc-pVnZ | −377.609 273 | −377.706 406 | −377.732 279 | |
CCSD(T)/aug-cc-pVnZ | −378.738 577 | −379.062 282 | −379.161 658 | |
ECorr.a/aug-cc-pVnZ | −1.129 304 | −1.355 876 | −1.429 379 | |
HCOOH | ||||
HF/cc-pVnZ | −188.779 739 | −188.840 624 | −188.856 020 | −188.859 688 |
CCSD(T)/cc-pVnZ | −189.308 960 | −189.501 934 | −189.561 659 | −189.581 158 |
ECorr.a/cc-pVnZ | −0.529 221 | −0.661 310 | −0.705 639 | −0.721 470 |
HF/aug-cc-pVnZ | −188.794 689 | −188.843 641 | −188.856 789 | |
CCSD(T)/aug-cc-pVnZ | −189.356 188 | −189.517 766 | −189.567 693 | |
ECorr.a/aug-cc-pVnZ | −0.561 498 | −0.674 126 | −0.710 904 | |
De (kcal/mol) | ||||
HF/cc-pVnZ | −15.27 | −12.48 | −11.90 | |
ECorr.a/cc-pVnZ | −3.14 | −4.46 | −4.59 | |
CCSD(T)/cc-pVnZ | −18.41 | −16.94 | −16.49 | |
HF/aug-cc-pVnZ | −12.48 | −12.00 | −11.74 | |
ECorr.a/aug-cc-pVnZ | −3.96 | −4.79 | −4.75 | |
CCSD(T)/aug-cc-pVnZ | −16.44 | −16.79 | −16.49 |
ECorr. = E[CCSD(T)] − E[HF].
Method . | E[HCOOH] . | E[(HCOOH)2] . | Dea . |
---|---|---|---|
HF | −188.8612 ± 0.0001 | −377.7409 ± 0.0003 | −11.7 ± 0.2b |
ECorr.c | −0.7290 ± 0.0006 | −1.4651 ± 0.0013 | −4.5 ± 1.0b |
CCSD(T) | −189.5896 ± 0.0006 | −379.2049 ± 0.0011 | −16.1 ± 0.9b,d |
Method . | E[HCOOH] . | E[(HCOOH)2] . | Dea . |
---|---|---|---|
HF | −188.8612 ± 0.0001 | −377.7409 ± 0.0003 | −11.7 ± 0.2b |
ECorr.c | −0.7290 ± 0.0006 | −1.4651 ± 0.0013 | −4.5 ± 1.0b |
CCSD(T) | −189.5896 ± 0.0006 | −379.2049 ± 0.0011 | −16.1 ± 0.9b,d |
De = E[(HCOOH)2] − 2 ⋅ E[HCOOH].
Error bar is obtained via .
Electron correlation energy: ECorr. = E[CCSD(T)] − E[HF].
An alternative estimate is obtained via De[CCSD(T)] = De[HF] + De[ECorr.] = − 16.2 ± 1.0 kcal/mol. The two values differ due to numerical errors at the extrapolation process.
For reasons of completeness, we also estimated the D0 value by including ZPE corrections. As described in Sec. II, the ZPE was estimated using the CCSD(T)/AVQZ harmonic frequencies and correcting for anharmonicity effects at the MP2/AVDZ level. To our knowledge, the presently reported harmonic frequencies are the most accurate in the literature since KKC reported harmonic frequencies only up to CCSD(T)/AVTZ. Our harmonic frequencies and anharmonicities, together with the experimentally measured fundamentals, are listed in Table VI. Note that the CCSD(T) harmonic frequency does not vary monotonically with basis set since it increases from AVDZ to AVTZ (as reported by KKC) and then decreases with the AVQZ basis set (our results). Therefore, the extrapolation used previously by KKC to estimate the CBS limits for the harmonic frequencies based on a monotonic variation with respect to the basis set is not justified from our CCSD(T)/AVQZ results. This, together with the use of B3LYP/AVTZ anharmonicities by KKC, is the reason for the larger deviation between their best-estimated anharmonic frequencies and experiment compared to ours. Their Mean Absolute Error () amounts to 33 cm−1, whereas ours to just 12 cm−1. Their largest absolute deviation from experiment is 152 cm−1 for the asymmetric combination (ag symmetry) of the two C–H stretches, whereas ours (for the same mode with a harmonic frequency of 3099 cm−1 in Table VI) is 79 cm−1, see also the Z-matrix definition and harmonic normal mode analysis in the supplementary material.54
Mode . | ωMP2 / AVDZ . | νMP2 / AVDZ . | MP2 anh. . | ωCCSD(T) / AVQZ . | νCCSD(T) / AVQZa . | νExp.b . | Δc . |
---|---|---|---|---|---|---|---|
HCOOH | |||||||
a′ | 618 | 610 | −7 | 630 | 623 | 628 | 5 |
a″ | 674 | 633 | −41 | 671 | 630 | 639 | 9 |
a″ | 1047 | 1023 | −23 | 1055 | 1032 | 1037 | 5 |
a′ | 1116 | 1082 | −33 | 1135 | 1102 | 1101 | −1 |
a′ | 1295 | 1198 | −97 | 1314 | 1218 | 1306 | 88 |
a′ | 1396 | 1385 | −11 | 1408 | 1397 | 1384 | −13 |
a′ | 1771 | 1739 | −32 | 1810 | 1777 | 1844 | 67 |
a′ | 3138 | 2999 | −140 | 3092 | 2952 | 2956 | 4 |
a′ | 3727 | 3529 | −197 | 3761 | 3563 | 3554 | −9 |
ZPEd | 21.13 | 20.30 | 21.27 | 20.43 | 20.66 | 0.23 | |
SDAe | 0.13 | ||||||
x0 | −0.04 | ||||||
(HCOOH)2 | |||||||
au | 69 | 68 | −1 | 72 | 72 | 69 | −3 |
ag | 169 | 161 | −8 | 167 | 159 | 161 | 2 |
au | 175 | 170 | −5 | 177 | 172 | 168 | −4 |
ag | 209 | 194 | −15 | 209 | 195 | 194 | −1 |
bg | 258 | 246 | −12 | 255 | 242 | 242 | 0 |
bu | 275 | 263 | −12 | 275 | 263 | 264 | 1 |
ag | 671 | 665 | −6 | 684 | 678 | 677 | −1 |
bu | 702 | 694 | −9 | 712 | 704 | 698 | −6 |
bg | 964 | 914 | −50 | 963 | 913 | 911 | −2 |
au | 981 | 947 | −34 | 986 | 952 | 922 | −30 |
bg | 1079 | 1052 | −27 | 1079 | 1052 | 1050 | −2 |
au | 1107 | 1067 | −40 | 1099 | 1059 | 1060 | 1 |
ag | 1237 | 1209 | −28 | 1252 | 1223 | 1214 | −9 |
bu | 1244 | 1211 | −33 | 1256 | 1223 | 1218 | −5 |
bu | 1386 | 1354 | −32 | 1405 | 1374 | 1364 | −10 |
ag | 1391 | 1358 | −33 | 1409 | 1376 | 1375 | −1 |
bu | 1450 | 1400 | −50 | 1455 | 1405 | 1415 | 10 |
ag | 1475 | 1424 | −52 | 1481 | 1429 | 1454 | 25 |
ag | 1688 | 1648 | −40 | 1713 | 1673 | 1670 | −3 |
bu | 1750 | 1714 | −37 | 1777 | 1740 | 1746 | 6 |
ag | 3141 | 2863 | −278 | 3099 | 2821 | 2900 | 79 |
bu | 3148 | 2981 | −167 | 3103 | 2936 | 2939 | 3 |
ag | 3158 | 2949 | −209 | 3204 | 2995 | 2949 | −46 |
bu | 3253 | 2986 | −267 | 3306 | 3038 | 3084 | 46 |
ZPEc | 44.29 | 42.23 | 44.51 | 42.45 | 42.52 | 0.07 | |
SDAd | 0.23 | ||||||
x0 | 0.01 |
Mode . | ωMP2 / AVDZ . | νMP2 / AVDZ . | MP2 anh. . | ωCCSD(T) / AVQZ . | νCCSD(T) / AVQZa . | νExp.b . | Δc . |
---|---|---|---|---|---|---|---|
HCOOH | |||||||
a′ | 618 | 610 | −7 | 630 | 623 | 628 | 5 |
a″ | 674 | 633 | −41 | 671 | 630 | 639 | 9 |
a″ | 1047 | 1023 | −23 | 1055 | 1032 | 1037 | 5 |
a′ | 1116 | 1082 | −33 | 1135 | 1102 | 1101 | −1 |
a′ | 1295 | 1198 | −97 | 1314 | 1218 | 1306 | 88 |
a′ | 1396 | 1385 | −11 | 1408 | 1397 | 1384 | −13 |
a′ | 1771 | 1739 | −32 | 1810 | 1777 | 1844 | 67 |
a′ | 3138 | 2999 | −140 | 3092 | 2952 | 2956 | 4 |
a′ | 3727 | 3529 | −197 | 3761 | 3563 | 3554 | −9 |
ZPEd | 21.13 | 20.30 | 21.27 | 20.43 | 20.66 | 0.23 | |
SDAe | 0.13 | ||||||
x0 | −0.04 | ||||||
(HCOOH)2 | |||||||
au | 69 | 68 | −1 | 72 | 72 | 69 | −3 |
ag | 169 | 161 | −8 | 167 | 159 | 161 | 2 |
au | 175 | 170 | −5 | 177 | 172 | 168 | −4 |
ag | 209 | 194 | −15 | 209 | 195 | 194 | −1 |
bg | 258 | 246 | −12 | 255 | 242 | 242 | 0 |
bu | 275 | 263 | −12 | 275 | 263 | 264 | 1 |
ag | 671 | 665 | −6 | 684 | 678 | 677 | −1 |
bu | 702 | 694 | −9 | 712 | 704 | 698 | −6 |
bg | 964 | 914 | −50 | 963 | 913 | 911 | −2 |
au | 981 | 947 | −34 | 986 | 952 | 922 | −30 |
bg | 1079 | 1052 | −27 | 1079 | 1052 | 1050 | −2 |
au | 1107 | 1067 | −40 | 1099 | 1059 | 1060 | 1 |
ag | 1237 | 1209 | −28 | 1252 | 1223 | 1214 | −9 |
bu | 1244 | 1211 | −33 | 1256 | 1223 | 1218 | −5 |
bu | 1386 | 1354 | −32 | 1405 | 1374 | 1364 | −10 |
ag | 1391 | 1358 | −33 | 1409 | 1376 | 1375 | −1 |
bu | 1450 | 1400 | −50 | 1455 | 1405 | 1415 | 10 |
ag | 1475 | 1424 | −52 | 1481 | 1429 | 1454 | 25 |
ag | 1688 | 1648 | −40 | 1713 | 1673 | 1670 | −3 |
bu | 1750 | 1714 | −37 | 1777 | 1740 | 1746 | 6 |
ag | 3141 | 2863 | −278 | 3099 | 2821 | 2900 | 79 |
bu | 3148 | 2981 | −167 | 3103 | 2936 | 2939 | 3 |
ag | 3158 | 2949 | −209 | 3204 | 2995 | 2949 | −46 |
bu | 3253 | 2986 | −267 | 3306 | 3038 | 3084 | 46 |
ZPEc | 44.29 | 42.23 | 44.51 | 42.45 | 42.52 | 0.07 | |
SDAd | 0.23 | ||||||
x0 | 0.01 |
Calculated from Eq. (5).
Δ = νExp. − νCCSD(T) / AVQZ (cm−1).
Sum of diagonal anharmonicities xii, SDA = − 1/4∑ixii.
IV. CONCLUSIONS
We report a thorough investigation of the binding energy of the FAD at the CCSD(T) level of theory with emphasis on the effects of the size of the basis set, the geometry optimization, and the BSSE correction. The two FAMs are largely distorted in the FAD geometry due to the formation of the two strong hydrogen bonds and therefore the deformation energy needs to be included for a proper estimate of the BSSE corrections. We found the convergence patterns of both the uncorrected and the BSSE-corrected binding energies for the titled system to be very similar to those observed earlier for hydrogen-bonded systems36 or weak interactions.48 Contrary to what was reported earlier by KKC for this system, we find that the BSSE correction does decrease (as expected) upon increasing the size of the basis set. KKC based their conclusions on an example presented earlier by Baerends and co-workers46 (Be2), which has a different convergence pattern of the uncorrected and BSSE-corrected binding energies with basis set compared to the FAD. In particular, the correlated result for Be2 is underbound for all basis sets so the BSSE uncorrected result is closer to the CBS limit than the BSSE-corrected one. As has been shown in Sec. III, for the present case, the uncorrected and BSSE-corrected binding energies converge from either side of the CBS limit with basis set, i.e., the convergence pattern is quite different than the one seen for Be2. Although from a pragmatic point of view, correcting for BSSE for Be2 will always produce results that are further away from CBS than the uncorrected result, the BSSE correction for that system is indeed theoretically justified since it will produce the same CBS estimate as the uncorrected result.
The discrepancy between our results and those reported earlier by KKC is attributed to their use of a valence basis set in conjunction with all electron (i.e., including the 1s of C, O) correlation. We have shown that this can indeed yield an anomalous, non-monotonic pattern of the BSSE correction with increasing basis set. In addition, by extrapolating the uncorrected or BSSE-corrected binding energies either separately or simultaneously (when possible), we arrive at CBS limits that are within 0.1 kcal/mol from one another. We obtained best estimates for De = − 16.1 ± 0.1 kcal/mol and D0 = − 14.3 ± 0.1 kcal/mol, the later in excellent agreement with the existing experimental value of −14.22 ± 0.12 kcal/mol. The estimate for D0 was based on the most accurate to date harmonic estimates at the CCSD(T)/AVQZ level of theory. Our results further justify the validity of the BSSE correction for this system, contrary to what suggested in previous studies.45,46
Acknowledgments
We thank Professor Ernest Davidson and Dr. Edoardo Aprà of PNNL for many helpful discussions and a critical review of the manuscript. This work was supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences, and Biosciences. Pacific Northwest National Laboratory (PNNL) is a multiprogram national laboratory operated for DOE by Battelle. A portion of this research was performed using the Molecular Science Computing Facility (MSCF) in EMSL, a national scientific user facility sponsored by the Department of Energy’s Office of Biological and Environmental Research and located at PNNL. This research also used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.