Accuracy and limitations of second-order many-body perturbation theory for predicting vertical detachment energies of solvated-electron clusters

John M. Herbert * and Martin Head-Gordon *
Department of Chemistry, University of California, Berkeley, CA 94720, USA. E-mail: jherbert@berkeley.edu; mhg@bastille.cchem.berkeley.edu

Received 15th September 2005 , Accepted 26th October 2005

First published on 14th November 2005


Abstract

Vertical electron detachment energies (VDEs) are calculated for a variety of (H2O)n and (HF)n isomers, using different electronic structure methodologies but focusing in particular on a comparison between second-order Møller–Plesset perturbation theory (MP2) and coupled-cluster theory with noniterative triples, CCSD(T). For the surface-bound electrons that characterize small (H2O)n clusters (n ≤ 7), the correlation energy associated with the unpaired electron grows linearly as a function of the VDE but is unrelated to the number of monomers, n. In every example considered here, including strongly-bound “cavity” isomers of (H2O)24, the correlation energy associated with the unpaired electron is significantly smaller than that associated with typical valence electrons. As a result, the error in the MP2 detachment energy, as a fraction of the CCSD(T) value, approaches a limit of about −7% for (H2O)n clusters with VDEs larger than about 0.4 eV. CCSD(T) detachment energies are bounded from below by MP2 values and from above by VDEs calculated using second-order many-body perturbation theory with molecular orbitals obtained from density functional theory. For a variety of both strongly- and weakly-bound isomers of (H2O)20 and (H2O)24, including both surface states and cavity states, these bounds afford typical error bars of ±0.1 eV. We have found only one case where the Hartree–Fock and density functional orbitals differ qualitatively; in this case the aforementioned bounds lie 0.4 eV apart, and second-order perturbation theory may not be reliable.


1. Introduction

Clusters comprised of closed-shell, polar molecules such as HF or H2O can form stable (or at least metastable) supramolecular radical anions that are often termed dipole-bound anions,1,2 a reference to the fact that the dipole moment μ0 of the neutral cluster frequently provides an important stabilizing interaction for the excess electron. In many such systems, however, dispersion interactions contribute significantly to electron binding,2–4 and even clusters with μ0 ≈ 0 can bind an extra electron quite strongly within a cavity.5 As such, accurate calculation of vertical electron detachment energies (VDEs) for this type of cluster anion is a challenging problem for electronic structure theory, because electron correlation is responsible for a significant fraction of the electron binding energy.1–3,6–9 Water cluster anions, in particular, have been intensely scrutinized by electronic structure theorists10 in the two decades since they were first observed in a molecular beam.11,12 Recent experiments13–16 have detected two new classes of (H2O)n isomers with VDEs that are significantly smaller, for a given n, than those of the dominant isomers observed in earlier experiments.17 These new data have rekindled an old controversy15–24 regarding whether the excess electron is solvated within the cluster or at its surface, and how this solvation mechanism evolves as a function of n.

To this debate, we have recently contributed a benchmark study5 of electronic structure predictions of VDEs for small water cluster anions (n ≤ 6), demonstrating that second-order Møller–Plesset perturbation theory (MP2), in conjunction with a modest one-electron basis set, underestimates VDEs by 0.03–0.05 eV, an accuracy that is sufficient to assign the experimental photoelectron spectra.13,25 A preliminary examination of selected isomers of (H2O)20 and (H2O)24, using the same methodology, cast considerable doubt on all previous electronic structure calculations for large water cluster anions, all of which utilize density functional theory (DFT) and probably overestimate VDEs by 0.2–0.4 eV.5 At the time, we speculated that these were probably the most accurate VDEs to date for cluster anions of this size, but it is difficult to quantify this statement absolutely since electronic structure calculations beyond the MP2 level are tremendously difficult (and rapidly become impossible) beyond six or seven water monomers. Moreover, for n > 6 there is no agreement as to which isomers of (H2O)n are actually interrogated in the experimental photoelectron spectra.

In the present study, we follow up on our previous work, with the goal of quantifying the accuracy of MP2 detachment energies for large water cluster anions. This can only be accomplished indirectly, for the reasons mentioned above; nevertheless, by extensively benchmarking MP2 detachment energies against coupled-cluster results, both for (H2O)n (with n ≤ 7) and also for (HF)n (with n ≤ 8), we can establish whether, and to what extent, the accuracy of MP2 theory degrades as a function of cluster size. In fact, we find that it does not; rather, the error in MP2 detachment energy predictions is an approximately constant fraction of the VDEs themselves. Scaled MP2 results provide VDEs for large cluster anions with typical estimated errors of ±0.1 eV.

2. Technical details

Drawing upon what we learned in ref. 5, all calculations reported in this work employ the 6-31(1+,3+)G* basis set unless otherwise stated.26 This basis augments the standard 6-31++G* basis with two additional diffuse s functions on each hydrogen atom, whose exponents are decremented by successive factors of 3.32. Comparison to larger and more diffuse basis sets indicates that this basis affords VDEs at both the MP2 and the CCSD(T) levels that are 0.01–0.02 eV below their respective complete-basis limits.5 For example, the cis-(H2O)2, cyclic (H2O)3, and linear (H2O)3 clusters introduced below afford MP2/6-31(1+,3+)G* VDEs of 0.026, 0.007, and 0.141 eV, respectively, whereas the corresponding MP2/TZVPP+3s2p values27 are 0.031, 0.021, and 0.139 eV.

Again consistent with our previous work,5 the calculations discussed here use drop tolerances of 10−12–10−14 au. Eigenvectors of the overlap matrix that correspond to eigenvalues smaller than the square root of the drop tolerance are projected out of the atomic orbital basis. Numerical integrations required for DFT exploit the SG-1 grid,29 and all self-consistent field (SCF) calculations are converged to a maximum occupied-virtual Fock matrix element of 10−8 au. Except where noted, all geometries are optimized using the B3LYP30,31 density functional with no constraints on point-group symmetry. Electronic structure calculations were performed using Q-Chem,32 and the Molden33 and VMD34 programs were used for visualization.

We have confirmed through numerous examples that sundry standard procedures for generating an SCF guess all lead to the same SCF solution, except in one high-symmetry case that is discussed in Section 3.2. In contrast, when converged molecular orbitals for (H2O)n are used as an initial guess for (H2O)n at the same geometry, we have found cases where the SCF calculation for the anion converges to a higher-energy (and qualitatively different) solution than that obtained with a standard SCF guess. For this reason, we never employ neutral molecular orbitals as an initial guess.

Fig. 1 depicts a catalogue of B3LYP/6-31(1+,3+)G* optimized geometries for (H2O)n isomers with n ≤ 7 and possessing CCSD(T)/6-31(1+,3+)G* detachment energies that range from nearly zero up to about 0.75 eV. This will serve as our test set for the calculations described in Section 3.1. In order to facilitate future benchmarking, the geometries and absolute energies (at various levels of theory) for all (H2O)n isomers discussed in this work are available as electronic supplementary information (ESI).


Isomers of (H2O)n− considered in this work, optimized at the B3LYP/6-31(1+,3+)G* level. Also indicated for each isomer is the dipole moment μ0 of the neutral cluster and the changes in natural atomic charges that accompany electron detachment, both calculated at the Hartree–Fock/6-31(1+,3+)G* level. The δqi values in black correspond to hydrogen atoms while those in gray are for oxygen atoms.
Fig. 1 Isomers of (H2O)n considered in this work, optimized at the B3LYP/6-31(1+,3+)G* level. Also indicated for each isomer is the dipole moment μ0 of the neutral cluster and the changes in natural atomic charges that accompany electron detachment, both calculated at the Hartree–Fock/6-31(1+,3+)G* level. The δqi values in black correspond to hydrogen atoms while those in gray are for oxygen atoms.

3. Results and discussion

3.1. Benchmark results for small clusters

Natural population analysis35,36 provides a convenient atomic partition of the singly-occupied molecular orbital (SOMO) coefficients from an SCF calculation, and is one way to assess and quantify which atoms are most directly responsible for solvating the extra electron. (Unlike Mulliken populations, natural populations are reasonably stable with respect to basis-set expansion.35) In our catalogue of (H2O)n isomers, Fig. 1, we indicate the changes δqi in natural atomic charges that accompany electron detachment. For clarity, δqi is given explicitly only for those atoms i where |δqi| ≥ 0.050, whereas atoms with 0.025 ≤ |δqi| < 0.050 are indicated by asterisks, and cases where |δqi| < 0.025 are omitted.

Except for certain weak-binding cyclic isomers, in which the excess electron is delocalized over many dangling hydrogen atoms, natural population analysis reveals that the excess electron is typically localized on two or three hydrogen atoms. In cases where the electron is bound to both hydrogen atoms on a single water molecule, this is accompanied by a small change in charge on the oxygen atom (δq ≈ 0.05–0.10). This is the so-called double hydrogen-bond acceptor or “AA” binding motif, and isomers of this sort represent the dominant feature in the dimer through hexamer photoelectron spectra.14,25,37,38 This inherent localization of the SOMO is consistent with previous electronic structure results indicating that diffuse basis functions on the hydrogen atoms are much more important in VDE calculations than are diffuse functions on the oxygen atoms, at least in (H2O)n clusters.5

In Fig. 2, we plot both the absolute and the fractional errors in the MP2 detachment energies for each of the isomers in Fig. 1. Here and elsewhere in this work, we define “error” to mean the difference with respect to the CCSD(T) value in the same basis set. As the predicted VDE increases, so does the difference between the MP2 and the CCSD(T) detachment energies, but this error is apparently uncorrelated with the number of monomers, n. On the other hand, the fraction of the CCSD(T) detachment energy recovered at the MP2 level appears to approach an asymptotic limit as the VDE increases, and for the ten (H2O)n isomers in our data set that have VDEs greater than 0.4 eV, MP2 theory recovers an average of 93.3 ± 1.0% of the CCSD(T) detachment energy. (Uncertainties reported in this work represent one standard deviation.) Moreover, this asymptotic behavior appears to be rather general: as shown in Fig. 3, VDEs obtained from full fourth-order Møller–Plesset perturbation theory (MP4) as well as from coupled-cluster theory at the singles-doubles level (CCSD) also approach limiting fractions of the CCSD(T) detachment energy as the VDE increases. These data indicate that MP4 theory recovers a larger fraction of the VDE than does MP2 theory, with CCSD recovering a larger fraction still. As demonstrated below, the electron correlation contribution to the VDE increases as the VDE itself increases, hence we expect this trend to continue in larger, higher-binding isomers.


Absolute errors (upper panel) and fractional errors (lower panel) in MP2 detachment energies, with respect to CCSD(T) values, for various isomers of (H2O)n−. All calculations use the 6-31(1+,3+)G* basis set.
Fig. 2 Absolute errors (upper panel) and fractional errors (lower panel) in MP2 detachment energies, with respect to CCSD(T) values, for various isomers of (H2O)n. All calculations use the 6-31(1+,3+)G* basis set.

Trends in the MP2 error noted above can best be understood by examining trends in the correlation energy Δ associated with the unpaired electron,

Δ = Ecorr(anion) − Ecorr(neutral).
We will use the notation ΔX to indicate a particular level of theory, X. Note that ΔX is also the correction to the Hartree–Fock VDE. The quantities ΔMP2 and ΔCCSD(T) are plotted in Fig. 4, and again we find no correlation whatsoever with cluster size but a clear trend with increasing VDE, namely, Δ is an approximately linear function of the VDE. Hartree–Fock VDEs and ΔCCSD(T) corrections for our catalogue of (H2O)n isomers are listed in Table 1.


Fractional errors, with respect to CCSD(T) values, in VDEs calculated for the (H2O)n− isomers in Fig. 1.
Fig. 3 Fractional errors, with respect to CCSD(T) values, in VDEs calculated for the (H2O)n isomers in Fig. 1.
Table 1 Hartree–Fock VDEs and CCSD(T) corrections to the VDE, ΔCCSD(T), for isomers of (H2O)n
Isomer Hartree–Fock VDE/eV Δ CCSD(T)/eV
n = 2 0.000 0.037
n = 3(A) −0.018 0.039
n = 3(B) 0.092 0.069
n = 4(A) 0.248 0.131
n = 4(B) 0.128 0.101
n = 4(C) 0.010 0.051
n = 4(D) 0.170 0.108
n = 4(E) 0.127 0.091
n = 4(F) 0.189 0.100
n = 5(A) 0.298 0.141
n = 5(B) 0.038 0.063
n = 5(C) 0.325 0.154
n = 5(D) 0.189 0.130
n = 5(E) 0.263 0.141
n = 5(F) 0.264 0.124
n = 6(A) −0.019 0.052
n = 6(B) 0.513 0.227
n = 6(C) −0.011 0.058
n = 6(D) 0.056 0.071
n = 6(E) 0.343 0.162
n = 6(F) 0.371 0.170
n = 6(G) 0.322 0.154
n = 7(A) 0.477 0.198
n = 7(B) 0.434 0.173
n = 7(C) 0.348 0.136


Averaging over all (H2O)n cluster sizes and isomers, we find that MP2 theory recovers 94.90 ± 0.07% of the CCSD(T) correlation energy, and 95.00 ± 0.07% in the case of the corresponding neutral clusters. (Uncertainties represent one standard deviation.) That these fractions are nearly geometry- and size-independent merely reflects the fact that the total correlation energy is dominated by the sum of monomer correlation energies, but together with the linear relationship between Δ and the VDE, this explains the asymptotic fractional error observed in Fig. 2. The increasing importance of electron correlation as a function of the VDE eventually overshadows any differences between isomers at the Hartree–Fock level, such that in strong-binding isomers, the error in the MP2 calculation is dominated solely by the difference between ΔMP2 and ΔCCSD(T).

In addition to wave function-based electronic structure methods, we also calculate VDEs using the BHLYP5,39 density functional (EBHLYPXC = 0.5ELDAX + 0.5EHFX + ELYPC), which tends to be the least overbinding of the common density functionals. Fig. 3 suggests that, for (H2O)n isomers with VDEs in excess of about 0.4 eV, BHLYP overestimates the VDE by about 15%, though the scatter among data points in the asymptotic region is greater than in the case of wave function-based methods. Somewhat more consistent behavior is obtained when BHLYP molecular orbitals are used in second-order many-body perturbation theory (MBPT2), including single excitations. We refer to this procedure as the MBPT2(BHLYP) method.40 As shown in Fig. 5, the CCSD(T) detachment energy is bounded from below by the MP2 value, which averages 93.2 ± 1.1% of the CCSD(T) detachment energy, and from above by MBPT2(BHLYP), at an average of 111.1 ± 1.4% of CCSD(T). In calculating these averages and standard deviations, we have included only those isomers whose CCSD(T) detachment energies exceed 0.4 eV, since we demonstrated above that it is only in these cases that MP2 affords an approximately constant fractional error. Our conclusions regarding the accuracy of second-order perturbation theory thus apply mainly to strongly-bound isomers (VDE ≳ 0.4 eV).


Correlation energy ΔX [X = MP2 or CCSD(T)] associated with the unpaired electron, for various isomers of (H2O)n−. Lines indicate best fits to the data.
Fig. 4 Correlation energy ΔX [X = MP2 or CCSD(T)] associated with the unpaired electron, for various isomers of (H2O)n. Lines indicate best fits to the data.

Having only a small number of data points above 0.4 eV, we cannot be certain that the apparently asymptotic behavior observed in the 0.4–0.8 eV range persists to higher detachment energies. For additional data, we turn to (HF)n clusters, since a relatively small number of HF molecules can bind an electron much more strongly than can a small water cluster. Optimized geometries for (HF)n, n = 3–8, are depicted in Fig. 6, and clearly show the evolution of a nascent cavity within the cluster, in which protons are oriented toward the excess electron. None of these isomers has a neutral dipole moment larger than 2.7 D and in most cases μ0 is considerably smaller than the critical value (≈2.4 D; see ref. 2) necessary to obtain anionic states that are rigorously dipole-bound. Nevertheless, these clusters exhibit VDEs as large as 4.0 eV. Along with their cavity-like structures, this fact implies that these clusters are true “solvated electron” systems, in contrast to small water cluster anions, which are sometimes termed “wet electron” systems41 to emphasize that the electron is not solvated in the bulk sense. The existence of solvated electron or “cavity states” (as we shall call them) at such a small number of monomers is important because it allows us to perform coupled-cluster calculations for such species, which are anticipated to be important in large water cluster anions but cannot be located in clusters small enough to be treated with coupled-cluster theory.


Predicted VDEs for the (H2O)n− isomers shown in Fig. 1. CCSD(T) VDEs lie along the diagonal (solid line); the upper and lower broken lines represent 111.1% and 93.2% of the CCSD(T) VDE, which represent the average fractions recovered by the MBPT2(BHLYP) and MP2 methods, respectively.
Fig. 5 Predicted VDEs for the (H2O)n isomers shown in Fig. 1. CCSD(T) VDEs lie along the diagonal (solid line); the upper and lower broken lines represent 111.1% and 93.2% of the CCSD(T) VDE, which represent the average fractions recovered by the MBPT2(BHLYP) and MP2 methods, respectively.

We digress briefly here to discuss some technical details related to electronic structure calculations for (HF)n. First of all, we observe that these clusters are not quite as insensitive to basis sets effects as are (H2O)n clusters. In (HF)3, for example, expansion of the basis from 6-31(m+,n+)G* to 6-311(m+,n+)G** decreases the MP2 detachment energy by about 0.06 eV, while additional diffuse shells added to the fluorine atoms increase the VDE by about 0.02 eV. The fact that basis set effects of this magnitude are not observed in (H2O)n is consistent with our earlier explanation5 that the bent geometry of a water molecule allows the SOMO (essentially an sp hybrid) to acquire p character by mixing together two s functions on different hydrogen atoms; as such, diffuse p functions are not absolutely required in the case of (H2O)n. The same argument does not apply to (HF)n, and for these systems it appears that p functions on the hydrogen atoms have an important effect. Furthermore, as demonstrated by the natural charge differences given in Fig. 6, the excess electron density is shared by a far greater number of atoms in the (HF)n cavity states than in most of the (H2O)n isomers shown in Fig. 1. Despite these considerations, we will stick to the 6-31(1+,3+)G* basis set for consistency and because our primary interest lies in (H2O)n, keeping in mind that our calculated VDEs for (HF)n isomers likely lie further below the complete-basis limit than they do in the case of (H2O)n.

With this in mind, it is interesting to note that our calculated VDEs, especially for the larger (HF)n cavity isomers, are substantially greater than those reported in older calculations.2,3 One possible explanation for the discrepancy is the fact that previous calculations focused on high-symmetry geometries, e.g., Oh symmetry for the hexamer anion.2 Although we do obtain fully-relaxed D∞h, D3h, and Td stationary points for the dimer, trimer, and tetramer anions, respectively, we failed in repeated attempts to locate a D3h stationary point for (HF)5 or an Oh stationary point for (HF)6, both of which have been reported elsewhere.2,3 At the B3LYP/6-31(1+,3+)G* level, both the pentamer and the hexamer anions adopt structures in which the excess electron is tetrahedrally coordinated, with the additional HF molecule(s) forming an incipient second solvation shell. A stationary point of (HF)6 with Oh symmetry can be located if point-group symmetry is constrained, and the resulting Oh structure has a VDE of 2.76 eV, compared to 3.07 eV for the fully-relaxed structure in Fig. 6. [Both VDEs are CCSD(T)/6-31(1+,3+)G* values.] In comparison, Jordan and Wang2 report a VDE of 2.69 eV for their Oh structure, also at the CCSD(T) level. On the other hand, our D3h trimer anion has a predicted VDE of 0.93 eV, much larger than the values of ≈0.6 eV reported by Gutowski and co-workers3,42 for comparable D3h structures. Further calculations—in particular, geometry optimizations at higher levels of theory, and a more thorough exploration of basis set effects—are necessary in order to resolve the true structures of these clusters, but for the purpose of benchmarking VDE calculations, the structures presented in Fig. 6 will suffice.

In addition to these cavity isomers of (HF)n, we also consider HF⋯HF⋯HF⋯ chains, artificially constrained to linear geometries but otherwise optimized at the B3LYP/6-31(1+,3+)G* level. These isomers bind an electron very weakly at the H-terminus and, as compared to the cavity isomers, are more analogous to the surface-bound anions obtained from small water clusters. Fig. 7 bolsters this analogy, as it shows that the fractional error in the MP2 detachment energy behaves in much the same way for these linear (HF)n isomers as it does for small (H2O)n isomers. (Numerical data corresponding to this figure are given in Table 2.) In contrast, for cavity isomers of (HF)n, MP2 predicts a larger VDE than does CCSD(T). Given their very different electron binding motifs, there is no reason why the MP2 fractional errors should behave the same way for linear versus globular (HF)n; the important observation is that, in both cases, the MP2 fractional error appears to approach an asymptotic limit as the VDE increases.


Optimized B3LYP/6-31(1+,3+)G* geometries of (HF)n−
						“cavity” isomers, along with VDEs [CCSD(T)/6-31(1+,3+)G* level] and selected δqi values [Hartree–Fock/6-31(1+,3+)G* level].
Fig. 6 Optimized B3LYP/6-31(1+,3+)G* geometries of (HF)n “cavity” isomers, along with VDEs [CCSD(T)/6-31(1+,3+)G* level] and selected δqi values [Hartree–Fock/6-31(1+,3+)G* level].
Table 2 Hartree–Fock VDEs and CCSD(T) corrections to the VDE, ΔCCSD(T), for isomers of (HF)n
Isomer Hartree–Fock VDE/eV Δ CCSD(T)/eV
n = 2 (linear) 0.004 0.023
n = 3 (linear) 0.106 0.041
n = 4 (linear) 0.183 0.053
n = 5 (linear) 0.240 0.061
n = 6 (linear) 0.280 0.065
n = 7 (linear) 0.310 0.068
n = 8 (linear) 0.332 0.069
n = 2 (cavity) −0.009 0.123
n = 3 (cavity) 0.721 0.213
n = 4 (cavity) 1.595 0.196
n = 5 (cavity) 2.314 0.168
n = 6 (cavity) 2.917 0.150
n = 7 (cavity) 3.406 0.144
n = 8 (cavity) 3.868 0.137


Tables 1 and 2 list ΔCCSD(T) for each (H2O)n and (HF)n isomer examined so far. In no case does this quantity exceed 0.23 eV, which is substantially smaller than the standard a priori estimate of ∼1 eV of correlation energy per valence electron pair. Not surprisingly, ΔCCSD(T) is somewhat larger for the cavity isomers of (HF)n than it is for the linear isomers, since the excess electron is bound more tightly and interacts more strongly with the cluster in the former case, but while ΔCCSD(T) (and, similarly, ΔMP2) demonstrates a clear surface state/internal state correlation, this quantity does not correlate with the number of H2O or HF monomers.

As a final set of small-cluster benchmarks, in Fig. 8 we compare absolute errors obtained from intermediate-level treatments of electron correlation, including CCSD, full MP4, and MP4 sans triple excitations, MP4(SDQ). Notably, the CCSD, MP4, and MP4(SDQ) methods exhibit essentially constant absolute errors for the (HF)n cavity isomers, where the VDEs are large. The same does not hold for the more weakly-binding (H2O)n isomers, although these methods remain more accurate than MP2. In future work, we shall attempt to quantify the importance of various individual correlations beyond second order in perturbation theory.


Fractional errors in MP2 detachment energies, with respect to CCSD(T) values, for isomers of (HF)n−.
Fig. 7 Fractional errors in MP2 detachment energies, with respect to CCSD(T) values, for isomers of (HF)n.

3.2. Comparison to large (H2O)n clusters

The results of the previous section provide strong evidence that MP2/6-31(1+,3+)G* detachment energies for (H2O)n lie within about 7% of CCDS(T) values for isomers with VDEs greater than 0.4 eV. For the dominant features in the photoelectron spectra, this régime is broached starting with n = 6, while larger clusters have even greater VDEs: 0.7 eV for (H2O)11, and 1.0–1.2 eV for (H2O)20 and (H2O)24, for example.15,23 In this section we attempt to address the accuracy of MP2 detachment energies for larger clusters, taking (H2O)20 and (H2O)24 as representative examples.

Fig. 9 presents the results of natural population analysis for several isomers of (H2O)20 and (H2O)24. These isomers were first studied in ref. 5, and the isomeric labels introduced in that work, which indicate the number of four-, five-, and six-member rings that comprise the cluster, are retained here. In order to make a fair comparison with the calculations in the previous section, we have re-optimized the geometries obtained in ref. 5 at the B3LYP/6-31(1+,3+)G* level, though this has no qualitative impact on the structures.

The three structures shown in the upper part of Fig. 9 are pentagonal dodecahedral or “512” isomers of (H2O)20 that bind the extra electron on the surface of the cluster, but via quite different morphologies, depending upon how the ten dangling hydrogen atoms are distributed. In 512 isomers A and B, the dangling hydrogens are asymmetrically distributed and the corresponding neutral clusters have enormous dipole moments (28.9 D and 21.9 D, respectively, at the Hartree–Fock level). As a result, the SOMO in each of these isomers is localized at one end of the cluster, which is clearly reflected in the δqi values. Significant polarization is found only on two surface hydrogen atoms, and none of the hydrogen atoms involved in hydrogen bonding is significantly polarized by electron attachment. Charge penetration into these (H2O)20 isomers does not appear to be any more significant than in the smaller water cluster anions depicted in Fig. 1, implying similar electronic structure in both cases.

Isomer 512F of (H2O)20 is markedly different from A and B. In this isomer, the ten dangling hydrogen atoms are located symmetrically around the equator of the cluster, so that μ0 ≈ 0 and the SOMO is a torus that is completely delocalized over these hydrogens [see Fig. 13(c) of ref. 5]. None of these hydrogens has a |δqi| value exceeding 0.050 but a large number of atomic charges lie in the range 0.025 ≤ |δqi| < 0.050. For this particular isomer, the electron binding mechanism appears to differ significantly from any that we observed in small cluster anions. On the other hand, the predicted VDE for this isomer (see Table 3) is far too small to be consistent with the experimental estimates of 1.0–1.2 eV.15,23

Table 3 Estimated VDEs for various isomers of (H2O)20 (upper data set) and (H2O)24 (lower data set), at X3LYP/6-31++G* optimized geometries (geometry I, from ref. 5); B3LYP/6-31(1+,3+)G* optimized geometries (geometry II); and B3LYP/6-31G* optimized geometries (geometry III)
    VDE/eV
Isomer   Hartree–Fock MP2 Corrected MP2a BHLYP MBPT2(BHLYP)
              Raw Scaledb
a Corrections to the MP2 VDE obtained by scaling the MP2 correlation energies by 1.053 (left column), or by scaling the MP2 VDE by 1.070 (right column). b VDE scaled by 0.8715.
445462 I 0.026 0.074 0.077 0.079 0.181 0.219 0.191
512A I 0.883 1.083 1.094 1.159 1.286 1.245 1.085
  II 0.903 1.101 1.111 1.178 1.291 1.254 1.093
512B I 0.723 0.908 0.918 0.972 1.084 1.062 0.925
  II 0.704 0.883 0.892 0.945 1.050 1.030 0.898
512C I 0.497 0.657 0.666 0.703 0.818 0.810 0.706
512D I 0.379 0.515 0.522 0.551 0.668 0.672 0.586
512E I −0.051 −0.022 −0.021 0.111 0.135 0.118
512F II 0.053 0.150 0.155 0.160 0.438 0.375 0.327
  III 0.093 0.227 0.234 0.243 0.559 0.483 0.421
                 
4668A I −0.086 −0.066 −0.065 0.010 0.084 0.073
4668B I 0.080 0.575 0.601 0.615 0.863 0.726 0.632
  II 0.177 0.583 0.605 0.624 0.838 0.716 0.624
41464A I −0.035 0.003 0.004 0.004 0.171 0.201 0.175
41464B I −0.035 0.004 0.004 0.004 0.421 0.346 0.302
51262A I −0.076 −0.051 −0.050 0.079 0.118 0.103
51262B I 0.254 0.793 0.821 0.849 1.100 0.927 0.808
  II 0.337 0.782 0.807 0.837 1.068 0.901 0.785
51262C I 0.049 0.136 0.192 0.146 0.430 0.362 0.316
414 I −0.077 −0.055 −0.055 0.048 0.141 0.123


Cavity isomers of large water cluster anions have different characteristics altogether, and Fig. 9 depicts two such isomers of (H2O)24. The electron binding mechanism is again apparent from the natural population analysis: δqi is quite large for each of the four hydrogen atoms that directly solvates the extra electron, but charge rearrangement is not limited to these atoms and is observed throughout the cluster. The SOMOs for these isomers (see Fig. 13 of ref. 5) corroborate this picture: the unpaired electron density has spherical symmetry and originates at the center of the cluster, but a non-negligible part of the SOMO extends beyond the cavity and into the water network. This same sort of through-cluster polarization was observed in cavity states of (HF)n, however, and did not affect the conclusion that MP2 calculations recover a constant fraction of the CCSD(T) correlation energy.

Hartree–Fock, MP2, and BHLYP detachment energies for these and several other (H2O)20 and (H2O)24 isomers studied in our previous work5 are listed in Table 3. For certain of these clusters, we have experimented with different density functionals (B3LYP and X3LYP43) and different basis sets for performing the geometry optimizations, in order to evaluate how much this aspect affects the VDE. Note also that ΔMP2 is given by the difference between the Hartree–Fock and MP2 detachment energies. In most cases ΔMP2 ≲ 0.02 eV, just as we found in the small cluster examples of Section 3.1. The exceptions, (H2O)24 isomers 4868B and 51262B, are cases where the extra electron is bound inside of the cluster. For these isomers, ΔMP2 ≈ 0.5 eV, still considerably smaller than the correlation energy associated with valence electron pairs.

Several of the isomers listed in Table 3 possess VDEs that are slightly negative at the Hartree–Fock level. This represents a basis set artifact; a complete basis would correctly describe continuum states and therefore any unbound anion would have a VDE of exactly zero. In the present cases, a higher-quality basis set might bind some of the anions at the Hartree–Fock level, but even without further calculations we can say with confidence that the resulting VDEs would be quite small, since our previous work5 suggests that basis set effects beyond 6-31(1+,3+)G* are typically on the order of 0.05 eV or less. As an example, we performed a Hartree–Fock VDE calculation for (H2O)24 isomer 41464B using the much larger TZVPP + 3s2p basis introduced in Section 2, which ought to be close to basis-set saturation at the SCF level. The resulting VDE, −0.027 eV, is nearly equal to the value of −0.035 eV obtained at the Hartree–Fock/6-31(1+,3+)G* level. On the one hand, this demonstrates that it is essentially impossible, in practice, to obtain a correct description of continuum states using atom-centered Gaussian basis functions. More importantly, however, this calculation offers additional evidence that improvements upon the 6-31(1+,3+)G* basis set are not likely to alter the conclusions reached in this study.

MP2 corrections to the negative Hartree–Fock VDEs in question are also quite small, affording VDEs that mostly remain negative, and in no case are greater than +0.004 eV. Clearly, the anions in question are either very weakly bound or else not bound at all, though a treatment of electron correlation beyond the MP2 level is required in order to make a definitive case for one of these possibilities over the other.

Table 3 also presents two corrected values of each MP2 detachment energy. The first is obtained by scaling the correlation energies of both the neutral and the anionic cluster by a factor of 1.053, the average ratio of the CCSD(T) to MP2 correlation energies for all of the data listed in Table 1. (ΔMP2 is an extensive quantity, so in order to preserve the size-consistency of the MP2 method, the neutral and anionic correlation energies must be scaled by the same factor.) This correction amounts to less than 0.03 eV in every case except isomer 51262C of (H2O)24, where it is almost 0.06 eV. In light of the consistent results obtained for small clusters in Section 3.1, and given that several other empirical scaling procedures for MP2 correlation energies have been shown to improve calculated properties significantly,44–48 we consider these scaled MP2 values to be our most accurate estimates of the VDEs for these clusters. An alternative correction provided in Table 3 is obtained by scaling the VDE itself by a factor of 1.070, obtained by a fit of the fractional error data (Fig. 2) for the nine isomers of (H2O)n whose MP2 detachment energies exceed 0.4 eV. This procedure results in a somewhat larger correction to the MP2 detachment energy, increasing the largest VDEs by 0.05–0.06 eV.

BHLYP detachment energies for our (H2O)20 and (H2O)24 isomers are also listed in Table 3, and, consistent with the results for smaller clusters, are significantly larger than the corresponding (raw) MP2 values. The MBPT2(BHLYP) detachment energy can be larger or smaller than the BHLYP value, with the two results generally falling within 0.1–0.2 eV of one another. Upon scaling the MBPT2(BHLYP) detachment energies by a factor of 0.8715, equal to the average ratio of MBPT2(BHLYP) and CCSD(T) detachment energies for the small (H2O)n isomers in our test set, all but one of the resulting VDEs are within 0.18 eV of the scaled MP2 values. We regard the electronic structure of these isomers to be a solved problem, and consider the adjusted MP2 results (using scaled correlation energies) to be the most accurate existing estimates of their VDEs.

The lone outlier is (H2O)24 isomer 41464B, which has a VDE of essentially zero at the MP2 level (corrected or uncorrected), but has an MBPT2(BHLYP) detachment energy of 0.302 eV after downward scaling. While even the BHLYP detachment energy for this isomer is far too small to account for the experimental photoelectron spectrum, the discrepancy between MP2 and MBPT2(BHLYP) is nevertheless worrisome insofar as it may portend a breakdown in the accuracy of second-order perturbation theory.

Examination of the Hartree–Fock and BHLYP SOMOs (Fig. 10) reveals that each is qualitatively consistent with the VDE obtained by subsequent second-order perturbation theory. In particular, the Hartree–Fock SOMO exhibits two diffuse lobes on the surface of the cluster, indicative of a weakly-bound anion. By changing the SCF initial guess (from a superposition of atomic densities to diagonalization of the core Hamiltonian), one can locate a second Hartree–Fock solution in which the SOMO consists of only one of these lobes, but which is only 2 × 10−5 au above the two-lobe solution shown in Fig. 10(a), and only 5 × 10−5 au higher in energy at the MP2 level. These energy differences are small but clearly distinct, given the thresholds discussed in Section 2, and indicate that this cluster exhibits a form of “electronic tautomerism”, at least at the Hartree–Fock and post-Hartree–Fock levels of theory. This tautomerism may be similar to the “bi-dipole-bound anion” effect, a symmetry artifact first discovered by Gutowski and co-workers49 in linear (HCN)n⋯HCCH⋯(NCH)n chains. In those systems, the SOMOs have σg symmetry with lobes located on either end of the chain, but in each case the SOMO is nearly degenerate with a corresponding σu virtual orbital. As such, there exist low-lying Hartree–Fock solutions, corresponding to excitation of the unpaired electron into σg ± σu, that localize the unpaired electron on one end of the molecule.


Errors in VDEs at different levels of theory, for (H2O)n− isomers (upper panel) and (HF)n−cavity isomers (lower panel). For (H2O)n− isomers, full MP4 calculations were performed only for n
						≤ 6.
Fig. 8 Errors in VDEs at different levels of theory, for (H2O)n isomers (upper panel) and (HF)ncavity isomers (lower panel). For (H2O)n isomers, full MP4 calculations were performed only for n ≤ 6.

Isomers of (H2O)20− (top row) and of (H2O)24− (bottom row), along with Hartree–Fock μ0 and δqi values. The δqi values in black correspond to hydrogen atoms while those in gray are for oxygen atoms.
Fig. 9 Isomers of (H2O)20 (top row) and of (H2O)24 (bottom row), along with Hartree–Fock μ0 and δqi values. The δqi values in black correspond to hydrogen atoms while those in gray are for oxygen atoms.

SOMOs for (H2O)24− isomer 41464B, calculated at (a) the Hartree–Fock/6-31(1+,3+)G* level and (b) the BHLYP/6-31(1+,3+)G* level. Both plots use a contour value of 0.005 au.
Fig. 10 SOMOs for (H2O)24 isomer 41464B, calculated at (a) the Hartree–Fock/6-31(1+,3+)G* level and (b) the BHLYP/6-31(1+,3+)G* level. Both plots use a contour value of 0.005 au.

In contrast to these Hartree–Fock solutions, the BHLYP SOMO contains significant density in the center of the cluster and penetrates throughout the water network, consistent with a nontrivial VDE. Starting with BHLYP orbitals and using an SCF algorithm based on direct minimization,50 we attempted to locate a Hartree–Fock solution whose SOMO resembles the BHLYP one, but this calculation merely relaxed to the previously-obtained Hartree–Fock solution. It remains unclear whether the Hartree–Fock or the BHLYP orbitals are a better model of physical reality in this cluster isomer. In such cases of electronic tautomerism, where different SCF procedures (e.g., Hartree–Fock versus DFT) result in qualitatively different molecular orbitals, finite-order perturbation theory predictions are probably not sufficient to determine the electronic structure. Such examples appear to be rare, as this is the only one we have found despite extensive calculations for this and previous work.5

4. Summary

We have examined whether the accuracy of MP2 detachment energy predictions for (H2O)n and (HF)n clusters degrades as a function of cluster size. In fact, the accuracy correlates not with n but with the VDE itself, a behavior that we attribute to the fact that, for (H2O)n isomers, the electron correlation energy Δ associated with the unpaired electron grows linearly as a function of the VDE. Even in cavity states of (H2O)24, which ought to be relatively strongly correlated, Δ ≲ 0.5 eV, only half as large as the correlation energy associated with typical valence electron pairs. This weak electron correlation explains the success of MP2 theory for VDE predictions.

Because the total correlation energy of a cluster is dominated by the sum of monomer contributions, MP2 theory recovers an essentially constant fraction of the CCSD(T) correlation energy in (H2O)n, (HF)n, and their anions, and as such the MP2 results can be easily adjusted in order to achieve better accuracy. An independent estimate of the VDE can be obtained by applying second-order many-body perturbation theory to DFT orbitals (obtained with the BHLYP functional), and in the small cluster isomers studied here, this procedure always affords an upper bound to the CCSD(T) detachment energy. For a wide variety of (H2O)20 and (H2O)24 isomers, these two estimates agree to within 0.15 eV, thus providing an error bar for MP2 detachment energies. The combination of MP2 and MBPT2(BHLYP) provides reasonably accurate VDE predictions, at low enough cost to be applied to fairly large clusters.

We have thus far discovered only one example in which MP2 and MBPT2(BHLYP) detachment energies differ by significantly more than 0.15 eV. In this case, Hartree–Fock theory exhibits a form of electronic tautomerism (which may be a symmetry artifact) and, more importantly, the SOMO is qualitatively different at the Hartree–Fock and BHLYP levels of theory. In this case, the error bar established by MP2 and MBPT2(BHLYP) is about 0.4 eV.

5. Note added in proof

We have recently calculated the VDE of isomer 41464B from section 3.2 at the RI-MP2/TZVPP+3s2p level, where RI indicates the resolution of identity approximation to MP2. (Our own tests and those in ref. 28 indicate that the RI approximation introduces errors of < 0.001 eV in VDEs.) The RI-MP2/TZVPP+3s2p result is 0.002 eV, as compared to 0.004 eV at the MP2/6-31(1+,3+)G* level.

Acknowledgements

This work was primarily supported by NSF grant CHE-9981997 (to M.H.-G.) and by an NSF postdoctoral fellowship (to J.M.H.). Additional support came from the Director, Office of Energy Research, Office of Basic Energy Sciences, Chemical Sciences Division of the US Department of Energy under Contract DE-AC03-76SF00098, and from supercomputer time from NERSC.

References

  1. M. Gutowski, K. D. Jordan and P. Skurski, J. Phys. Chem. A, 1998, 102, 2624 CrossRef CAS.
  2. K. D. Jordan and F. Wang, Annu. Rev. Phys. Chem., 2003, 54, 367 CrossRef CAS.
  3. M. Gutowski and P. Skurski, J. Phys. Chem. B, 1997, 101, 9143 CrossRef CAS.
  4. F. Wang and K. D. Jordan, J. Chem. Phys., 2001, 114, 10717 CrossRef CAS.
  5. J. M. Herbert and M. Head-Gordon, J. Phys. Chem. A, 2005, 109, 5217 CrossRef CAS.
  6. M. Gutowski, P. Skurski, A. I. Boldyrev, J. Simons and K. D. Jordan, Phys. Rev. A, 1996, 54, 1906 CrossRef CAS.
  7. M. Gutowski, P. Skurski, K. D. Jordan and J. Simons, Int. J. Quantum Chem., 1997, 107, 2968 CAS.
  8. P. Skurski and M. Gutowski, J. Chem. Phys., 1999, 111, 3004 CrossRef CAS.
  9. P. Skurski, M. Gutowski and J. Simons, J. Chem. Phys., 2001, 114, 7443 CrossRef CAS.
  10. Consult ref. 5 for a bibliography of electronic structure calculations for (H2O)n clusters.
  11. H. Haberland, C. Ludewigt, H.-G. Schindler and D. R. Worsnop, J. Chem. Phys., 1984, 81, 3742 CrossRef.
  12. H. Haberland, H. Langosch, H.-G. Schindler and D. R. Worsnop, J. Phys. Chem., 1984, 88, 3903 CrossRef CAS.
  13. J. Kim, I. Becker, O. Cheshnovsky and M. A. Johnson, Chem. Phys. Lett., 1998, 297, 90 CrossRef CAS.
  14. N. I. Hammer, J. R. Roscioli and M. A. Johnson, J. Phys. Chem. A, 2005, 109, 7896 CrossRef CAS.
  15. J. R. R. Verlet, A. E. Bragg, A. Kammrath, O. Cheshnovsky and D. M. Neumark, Science, 2005, 307, 93 CrossRef CAS.
  16. A. E. Bragg, J. R. R. Verlet, A. Kammrath, O. Cheshnovsky and D. M. Neumark, J. Am. Chem. Soc., 2005, 127, 15283 CrossRef CAS.
  17. J. V. Coe, G. H. Lee, J. G. Eaton, S. T. Arnold, H. W. Sarkas, K. H. Bowen, C. Ludewigt, H. Haberland and D. R. Worsnop, J. Chem. Phys., 1990, 92, 3980 CrossRef CAS.
  18. R. N. Barnett, U. Landman, C. L. Cleveland and J. Jortner, J. Chem. Phys., 1988, 88, 4429 CrossRef CAS.
  19. R. N. Barnett, U. Landman, C. L. Cleveland and J. Jortner, Chem. Phys. Lett., 1988, 145, 382 CrossRef CAS.
  20. R. N. Barnett, U. Landman, G. Rajagopal and A. Nitzan, Israel J. Chem., 1990, 30, 85 Search PubMed.
  21. P. Ayotte and M. A. Johnson, J. Chem. Phys., 1997, 106, 811 CrossRef CAS.
  22. D. M. Bartels, J. Chem. Phys., 2001, 115, 4404 CrossRef CAS.
  23. J. V. Coe, Int. Rev. Phys. Chem., 2001, 20, 33 CrossRef CAS.
  24. L. Turi, W.-S. Sheu and P. J. Rossky, Science, 2005, 309, 914 CrossRef CAS.
  25. J.-W. Shin, N. I. Hammer, J. M. Headrick and M. A. Johnson, Chem. Phys. Lett., 2004, 399, 349 CrossRef CAS.
  26. The exponents on these additional diffuse functions differ slightly from the 6-31(1+,3+)G* exponents defined in our previous work (ref. 5), which were decremented by factors of 3.00 rather than 3.32.
  27. The TZVPP + 3s2p basis consists of Ahlrichs’ doubly-polarized, valence triple zeta (TZVPP) basis for the oxygen atoms, the singly-polarized TZVP basis for the hydrogen atoms and a 3s2p set of diffuse basis functions centered on each oxygen and each hydrogen atom, with exponents taken from ref. 28.
  28. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 1999, 1, 4537 RSC.
  29. P. M. W. Gill, B. G. Johnson and J. A. Pople, Chem. Phys. Lett., 1993, 209, 506 CrossRef CAS.
  30. A. D. Becke, J. Chem. Phys., 1993, 98, 5648 CrossRef CAS.
  31. P. J. Stephens, J. F. Devlin, C. F. Chabalowski and M. J. Frisch, J. Phys. Chem., 1994, 98, 11623 CrossRef CAS.
  32. J. Kong, C. A. White, A. I. Krylov, C. D. Sherrill, R. D. Adamson, T. R. Furlani, M. S. Lee, A. M. Lee, S. R. Gwaltney, T. R. Adams, C. Ochsenfeld, A. T. B. Gilbert, G. S. Kedziora, V. A. Rassolov, D. R. Maurice, N. Nair, Y. Shao, N. A. Besley, P. E. Maslen, J. P. Dombroski, H. Daschel, W. Zhang, P. P. Korambath, J. Baker, E. F. C. Byrd, T. Van Voorhis, M. Oumi, S. Hirata, C. Hsu, N. Ishikawa, J. Florian, A. Warshel, B. G. Johnson, P. M. W. Gill, M. Head-Gordon and J. A. Pople, J. Comput. Chem., 2000, 21, 1532 CrossRef CAS.
  33. G. Schaftenaar and J. H. Noordik, J. Comput. Aided Mol. Des., 2000, 14, 123 CrossRef CAS.
  34. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33 CrossRef.
  35. A. E. Reed, R. B. Weinstock and F. Weinhold, J. Chem. Phys., 1985, 83, 735 CrossRef CAS.
  36. F. Weinhold, in Encyclopedia of Computational Chemistry, eds. P. v. R. Schleyer, N. L. Allinger, T. Clark, J. Gasteiger, P. A. Kollman, H. F. Schaefer III and P. R. Schreiner, Wiley, Chichester, United Kingdom, 1999, Vol. 3, p. 1792 Search PubMed.
  37. N. I. Hammer, J.-W. Shin, J. M. Headrick, E. G. Diken, J. R. Roscioli, G. H. Weddle and M. A. Johnson, Science, 2004, 306, 675 CrossRef CAS.
  38. N. I. Hammer, J. R. Rosciolli, M. A. Johnson, E. M. Myshakin and K. D. Jordan, J. Phys. Chem. A DOI:10.1021/jp053769c.
  39. J. C. Rienstra-Kiracofe, G. S. Tschumper, H. F. Schaefer III, N. Sreela and G. B. Ellison, Chem. Rev., 2002, 102, 231 CrossRef CAS.
  40. All of the MP2 calculations reported here employ frozen core orbitals, but in our MBPT2(BHLYP) calculations all orbitals are correlated.
  41. K. S. Kim, I. Park, S. Lee, K. Cho, J. Y. Lee, J. Kim and J. D. Joannopoulos, Phys. Rev. Lett., 1996, 76, 956 CrossRef CAS.
  42. M. Gutowski, C. S. Hall, L. Adamowicz, J. H. Hendricks, H. L. de Clercq, S. A. Lyapustina, J. M. Nilles, S.-J. Xu and K. H. Bowen Jr, Phys. Rev. Lett., 2002, 88, 143001 CrossRef CAS.
  43. X. Xu and W. A. Goddard III, Proc. Natl. Acad. Sci. USA, 2004, 101, 2673 CrossRef CAS.
  44. M. Gerenkamp and S. Grimme, Chem. Phys. Lett., 2004, 392, 229 CrossRef CAS.
  45. I. Hyla-Kryspin and S. Grimme, Organometallics, 2004, 23, 5581 CrossRef CAS.
  46. S. Grimme, J. Phys. Chem. A, 2005, 109, 3067 CrossRef CAS.
  47. Y. Jung, R. C. Lochan, A. D. Dutoi and M. Head-Gordon, J. Chem. Phys., 2004, 121, 9793 CrossRef CAS.
  48. R. C. Lochan, Y. Jung and M. Head-Gordon, J. Phys. Chem. A, 2005, 109, 7598 CrossRef CAS.
  49. M. Gutowski, P. Skurski and J. Simons, Int. J. Mass Spectrom., 2000, 201, 245 CrossRef CAS.
  50. T. Van Voorhis and M. Head-Gordon, Mol. Phys., 2002, 100, 1713 CrossRef.

Footnote

Electronic supplementary information (ESI) available: Geometries and energies. See DOI: 10.1039/b513098k

This journal is © the Owner Societies 2006
Click here to see how this site uses Cookies. View our privacy policy here.