Benchmarking Ab Initio Binding Energies of Hydrogen-Bonded

Jul 2, 2014 - Nicolai Bork†‡, Lin Du†, Heidi Reiman§, Theo Kurtén§, and Henrik G. Kjaergaard*†. † Department of Chemistry, University of ...
0 downloads 0 Views 556KB Size
Article pubs.acs.org/JPCA

Benchmarking Ab Initio Binding Energies of Hydrogen-Bonded Molecular Clusters Based on FTIR Spectroscopy Nicolai Bork,†,‡ Lin Du,† Heidi Reiman,§ Theo Kurtén,§ and Henrik G. Kjaergaard*,† †

Department of Chemistry, University of Copenhagen, Universitetsparken 5, DK-2100 Copenhagen Ø, Denmark Department of Physics, University of Helsinki, P.O. Box 64, FI-00014 Helsinki, Finland § Department of Chemistry, University of Helsinki, P.O. Box 55, FI-00014 Helsinki, Finland ‡

S Supporting Information *

ABSTRACT: Models of formation and growth of atmospheric aerosols are highly dependent on accurate cluster binding energies. These are most often calculated by ab initio electronic structure methods but remain associated with significant uncertainties. We present a computational benchmarking study of the Gibbs free binding energies in molecular complexes and clusters based on gas phase FTIR spectroscopy. The acetonitrile−HCl molecular complex is identified via its redshifted H−Cl stretching vibrational mode. We determine the Gibbs free binding energy, ΔG°295 K, to between 4.8 and 7.9 kJ mol−1 and compare this range to predictions from several widely used electronic structure methods, including five density functionals, Møller− Plesset perturbation theory, and five coupled cluster methods up to CCSDT quality, considering also the D3 dispersion correctional scheme. With some exceptions, we find that most electronic structure methods overestimate ΔG°295 K. The effects of vibrational anharmonicity is approximated using scaling factors, reducing ΔG°295 K by ca. 1.8 kJ mol−1, whereby ΔG°295 K predictions well within the experimental range can be obtained.



INTRODUCTION During the last decades it has been discovered that hydrogenbonded acid−base clusters are responsible for large fractions of the naturally occurring nanometer-scale aerosol particles in the atmosphere.1−4 Aerosol particles are known to induce significant cloud formation, and as a result of the large radiative forcing of aerosols and clouds, aerosol formation represents one of the largest uncertainties of current climate models.5−7 The growth rate of the acid−base cluster is the difference between the rate of sticking collisions and the rate of evaporation. Several parameters are somewhat uncertain (e.g., the sticking probability and the collision cross section) but the largest uncertainty is related to the Gibbs free binding energy of the clusters (ΔG°) because the evaporation rate scales exponentially with ΔG°.8,9 This problem is enhanced by the absence of reliable experimental data, leaving ab initio calculations as the main source of data. Although this is wellknown and much effort has been devoted to computational benchmarking studies, most ab initio data remain inferior to careful experimental data. At thermal equilibrium, the Gibbs free binding energy (ΔG°) and equilibrium constant (Keq) are readily determined from the concentrations of reactants (pA and pB) and product (pAB) of a given clustering reaction

° = Keq

pA /p°·pB /p°

⎛ −ΔG° ⎞ ⎟ = exp⎜ ⎝ RT ⎠

(2)

where R is the gas constant, T is the absolute temperature, and p° is the standard pressure (1 bar). Until now, mainly mass spectrometry has been used to measure the cluster and monomer concentrations of the strongly bonded systems of atmospheric relevance.10−12 By construction, mass spectrometry can only detect charged clusters, but even the most conservative charging mechanisms (e.g., proton transfer to gaseous nitrate), may induce significant cluster processing inside the spectrometer.13,14 Hereby, the detected concentrations of molecules and clusters may deviate significantly from the ambient, charge-neutral equilibrium. On the contrary, many recent studies have demonstrated that infrared spectroscopy (IR) is a viable way for obtaining accurate values for ΔG° and Keq.15−20 Moreover, the experimental and computational uncertainties are recognized, and reliable error estimates can be made. However, these studies have mainly focused on weaker bonded systems, for instance, alcohols and CHCl3 binding to amines and NH3, as well as SO2. Because atmospheric new particle formation mainly is driven by H2SO4 and nitrogen-containing bases,1−3 the results derived Received: April 16, 2014 Revised: June 20, 2014

(1)

A + B ↔ AB

pAB /p°

via the law of mass action © XXXX American Chemical Society

A

dx.doi.org/10.1021/jp5037537 | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

functional is here investigated by the D3 scheme by Grimme et al.23 In the study of DMS-HCl, we used both the 6-311+ +G(3df,3pd) and aug-cc-pV(T+d)Z [AV(T+d)Z] basis sets, but we found only minor differences in energies and configurations. Here, we use the AV(X+d)Z basis sets, X = (D, T, and Q), which are Dunning type correlation-consistent basis sets augmented with diffuse functions on all atoms with additional tight d basis functions on chlorine. The above-mentioned methods and basis sets were used to optimize geometries and calculate harmonic vibrational frequencies of HCl, MeCN, and the MeCN−HCl complex. Entropies were calculated using the harmonic oscillator and rigid rotor approximations. The Gibbs free binding energies were either used directly or corrected by coupled cluster (CC) single-point energies, yielding improved Gibbs free binding energies as

for weaker bonded systems might not be representative. We have therefore conducted an IR study of the molecular complex of acetonitrile (MeCN) and HCl, aimed at determining its Gibbs free binding energy. MeCN−HCl is a significantly more appropriate model system, since it involves a strong acid binding to a nitrogen atom, while still being weakly enough bonded that no noticeable cluster growth beyond the 1:1 cluster takes place in the gas phase. Further, the small system size enables utilization of very high quality computational methods. The MeCN−HCl complex has previously been observed in both solid argon and in solid N2 matrices,21 but has, to the best of our knowledge, not been observed in the gas phase, nor have any computational studies been published. We will use our obtained thermodynamic data to evaluate the accuracy of several of the ab initio methods most commonly applied for evaluating binding energies of clusters of atmospheric relevance. Opposed to most computational benchmarking studies regarding noncovalently bonded systems, we benchmark against an experimentally based value with known uncertainty. The tested methods include density functional theory (DFT), Møller−Plesset perturbation theory (MP2), and coupled cluster theory (CC) with basis sets up to augmented quadruple-ζ quality. The results will be compared to previous experimental and theoretical studies on weaker bonded complexes.18−20,22 In view of these results, we will draw conclusions on the ability of the tested methods in describing hydrogen-bonded systems of atmospheric relevance.

* ° − ΔE DFT + ΔECC ΔG° = ΔG DFT

(3)

where ΔE denotes the electronic energy and * denotes that the structure is not optimized at that level of theory. For this correction, we use CC2, CCSD, and CCSD(T) with the abovementioned basis sets and CCSD(T)-F12b with the VTZ-F12 basis set. Further, CCSDT calculations with the cc-pV(D+d)Z basis set were used to probe the effect of the approximate treatment of triple excitations in the CCSD(T) and CCSD(T)F12 calculations. Basis set superposition error is a well-known phenomenon artificially strengthening intermolecular interaction and is particularly important when using small basis sets on weakly interacting systems. This error may be compensated for by the counterpoise method by adding a correctional term to the binding energy of eq 3, given as



EXPERIMENTAL METHODS Anhydrous MeCN (99.8%) was purchased from Sigma-Aldrich and purified further by consecutive rounds of freezing, pumping, and thawing. Pressurized HCl (99.8%) was purchased from Air Liquide and used without further processing. IR spectra were recorded at 1.0 cm−1 resolution in the range of 1000 to 6500 cm−1 at room temperature (ca. 295 K). The spectrometer, a Vertex 70 FTIR (Bruker), was equipped with a CaF2 beam splitter, an MCT detector, and a globar type MIR light source. A 2.5 mm aperture was applied. Gaseous MeCN and HCl was loaded into 20 cm glass cells with KBr windows via a vacuum line with a base pressure lower than 10−3 Torr. To ensure complete mixing and thermal equilibrium, the samples were mixed 30 min before the spectra were recorded. All pressures were measured with Varian PCG-750 gauges, with uncertainties less than 1%. Trace amounts of HBr (less than 1% of HCl) from the reaction between HCl and the KBr windows were observed but did not interfere noticeably with the measurements.

A AB B AB ΔECP = EAB (A) − EAB (A) + EAB (B) − EAB (B)

(4)

EZY(X)

where denote the electronic energy of species X, in conformation Y, using the basis set of species Z.24,25 The applicability of this method with respect to the present system, methods, and basis set will be evaluated. All DFT and MP2 calculations were performed with the Gaussian 09 program (revision B.01). (Revision D.01 was used for the B3LYP-D3 calculations.) All coupled cluster calculations were performed using the Molpro program (version 2010.1), except the CCSDT calculations, which were performed using the NWChem program (version 6.1.1). Further descriptions and references to the methods and basis sets used in this study are available as Supporting Information. The binding energy was calculated using eq 2. The partial pressure of the MeCN−HCl complex was calculated as26



COMPUTATIONAL METHODS Benchmarking studies are fundamentally important for assessing the reliability of a given electronic structure calculation. In a previous study on the DMS−HCl molecular cluster,20 we compared the IR-based value of ΔG° to several ab initio based values. In that study, we used the B3LYP, CAMB3LYP, PW91, M06-2X, and wB97XD density functionals, including also Møller−Plesset perturbation theory calculations (MP2) which is of comparable quality and expense to the best DFT functionals. Here, we use the same methods to investigate the consistency of our previous results and compare with other theoretical and experimental studies on similar systems. Further, the effects of dispersion corrections to the B3LYP

pcomplex =

c ln(10) (l

∫ A (ν)̃ dν ̃

(5)

where c is speed of light in vacuum, l is the length of the sample cell, and A(ν̃) is the wavenumber dependent absorbance observed in our FTIR spectrum of the complex. ( is the integrated absorption coefficient given by

(=f

NAe 2 4 me c ε0

(6)

where NA, e, me, and ε0 denotes Avogadro’s constant, the elementary charge, the electron mass, and the vacuum permittivity, respectively.26 The oscillator strength, f, was B

dx.doi.org/10.1021/jp5037537 | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

calculated using the anharmonic oscillator local mode model, described in detail in the Supporting Information and in several other papers.16−20,22,27−29 For molecules, this may be less than 20% in error,30−32 but for the X−H stretching mode involved in hydrogen bonding in molecular complexes, the sparse available data suggests that f may be as much as a factor of 2 overestimated.33



RESULTS AND DISCUSSION Experimental Results. A total of seven FTIR spectra of mixtures of HCl and MeCN were recorded. Scaled reference spectra of pure HCl and pure MeCN were subtracted, yielding the spectra shown in Figure 1. This provided the actual Figure 2. (a) Integrated absorbance between 2600 and 2900 cm−1 and (b) calculated values of ΔG°295 K as functions of the product of partial pressures of MeCN and HCl. The uncertainties in (a) are from the calculated pressures of HCl and MeCN. The uncertainties in (b) are mainly from the calculated oscillator strength and the assignment of spectral features relating to this value.

mized using M06-2X/AV(T+d)Z. No stable T-shaped complex was found and the linear, hydrogen-bonded complex was found to be significantly more stable than any other configuration. The linear structure was further optimized using CCSD(T) with the AV(D+d)Z, AV(T+d)Z, and AV(Q+d)Z basis sets. The final structure is shown in Figure 3. Cartesian coordinates Figure 1. Infrared spectra of the redshifted HCl stretch frequency in the HCl−MeCN complex at varying HCl and MeCN partial pressures. Spectra of pure HCl and pure MeCN has been subtracted, yielding the pressures (in Torr) given in the inset table. Figure 3. Geometry of the MeCN−HCl molecular complex optimized at the CCSD(T)/AV(Q+d)Z level of theory, including some descriptive parameters (Cartesian coordinates are given in Table S1).

pressures of the equilibrated HCl and MeCN in the gas mixture. We assigned a 5% uncertainty to these pressures. A feature between 2600 and 2900 cm−1 was readily observable in all gas mixtures. The feature consist of a broad, weak band between ca. 2600 and 2700 cm−1 and a stronger asymmetric band between ca. 2700 and 2900 cm−1. In all gas mixtures, the relative intensities of these bands were ca. 1:10. The integrated absorbance from 2600 to 2900 cm−1 was found to scale linearly with the partial pressures of both HCl and MeCN, shown in Figure 2a, and both bands are thus concluded to originate from the MeCN−HCl complex. At the concentrations of HCl used, several of the strongest HCl modes were optically dense and hence impossible to subtract accurately, affecting especially the stronger band between 2700 and 2900 cm−1. The presence of multiple bands in this region is consistent with the matrix isolation studies by Schriver et al.21 who attributed two peaks at 2662 and 2655 cm−1 to two different conformers of the MeCN−HCl complex. One linear conformation connected by an N−H bond, and one T-shaped conformation connected by a Cl−C bond. Other studies have attributed such neighboring spectral features to coupling of vibrational modes, resulting typically in two additional bands or shoulders around a major peak.18,20,34 In the present complex, these are likely to originate from coupling of the intramolecular H−Cl stretch mode to the intermolecular N−Cl stretch mode. The latter mode could be similar to a mode of 111 cm−1 found in spectra of HCN−HCl in an argon matrix, attributed to the intermolecular N−Cl stretch.35 This was further investigated by electronic structure calculations. Initially, several starting geometries were opti-

are given in Table S1. To the best of our knowledge, no direct structural data on this complex exist, but a microwave study of the HCN−HCl complex concluded that a linear complex was vastly dominant,35 reporting a Cl−N distance of 3.40 Å comparable to the 3.29 Å found here. Because MeCN is a stronger base than HCN, the shortened bond distance is expected. Anharmonic Local Mode Model. Using the anharmonic oscillator local mode model, we calculated the frequencies of the fundamental H−Cl stretch mode in the MeCN−HCl complex. The results are summarized in Table 1. Further data are given in Tables S2 and S3. As expected, we find that the H− Cl stretch mode is redshifted in the complex compared to free gaseous HCl. The redshift is predicted to be around 220 cm−1, and we conclude that the main band with maximum intensity at 2725 cm−1 (161 cm−1 redshifted) is most likely from the fundamental H−Cl vibration in the MeCN−HCl complex. Because we found no T-shaped complex as proposed by Schriver et al.,21 we conclude that the minor band between 2600 and 2700 cm−1 is most likely from a subtractive combination mode involving the main mode at 2725 cm−1 and the intermolecular N−Cl mode. This was further verified by a harmonic frequency analysis at the CCSD(T)/AV(T+d)Z level of theory, showing an intermolecular H−N stretch mode at 116 cm−1 (Table S4). Most often, a subtractive combination mode is accompanied by additive combination mode, C

dx.doi.org/10.1021/jp5037537 | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

Table 1. Calculated Absorption Wavenumbers, v,̃ Redshifts, Δv,̃ and Oscillator Strengths, f, of the H−Cl Stretch Mode in the MeCN−HCl Complex and in Molecular HCl from the Anharmonic Oscillator Local Mode Model using CCSD(T) and the Indicated Basis Set MeCN−HCl

a

−1

HCl −4

−1

basis set

ṽ (cm )

f ( × 10 )

ṽ (cm )

f ( × 10−6)

Δṽ (cm−1)

AV(D+d)Z AV(T+d)Z AV(Q+d)Z expt

2695 2685 2678 2725

1.57 1.67 1.70

2922 2897 2895 2886

3.64 5.38 7.15 5.37a

227 212 217 161

From the HITRAN database.36

producing peaks at ν1 ± ν2 where, in this case, ν1 = 2725 cm−1 and ν2 = 116 cm−1.18,20,34 However, the presence of an additive combination band at (2725 + 116) cm−1 could neither be verified nor dismissed due to the irregular shape of the main band and due to the incomplete subtraction of the sharp HCl lines. Still, all available evidence supports the conclusion that both visible bands shown in Figure 1 originate from the linear, hydrogen-bonded MeCN−HCl molecular complex. The anharmonic oscillator local mode model further provided the oscillator strength, f (Table 1) required to determine the pressure of the MeCN−HCl complex. The small variations in calculated f suggest that the level of theory is adequate, but when the local mode model is applied to molecular clusters, fundamental uncertainties prevail, and overestimating oscillator strengths by up to a factor of 2 cannot be excluded.33 We thus conclude that, most likely, the oscillator strength of the H−Cl mode in the MeCN−HCl complex is between 0.85 × 10−4 and 1.70 × 10−4, and we use this range in the following. Spectral Deconvolution. The spectra were deconvoluted using four Lorentzian functions (Figure S1 and Table S5). One Lorentzian contributed to the small band between 2600 and 2700 cm−1 and three contributed to the main band between 2700 and 2900 cm−1. Two of these three Lorentzians contribute to the absorbance around 2725 cm−1 and are just ca. 25 cm−1 apart. These two may originate from modes split by Fermi resonance or may be a single mode poorly described by a singe Lorentzian. The third Lorentzian contributing to the main band is centered at ca. 2790 cm−1 and may represent the above-mentioned additive absorption, although very redshifted from the predicted 2840 cm−1. The area, denoted ∫ A(ν̃)dν̃ in eq 2, is obtained on the basis of this deconvolution. The calculated oscillator strength includes the fundamental H−Cl stretch mode only, and it is uncertain how it relates to the observed combination modes. Cautiously, we chose to use the sum of the two Lorentzians contributing to the absorbance around 2725 cm−1 as a lower limit to ∫ A(ν̃)dν̃ and the sum of all four Lorentzians as an upper limit to ∫ A(ν̃)dν̃. The lower limit is about 2/3 of the upper limit. We calculate the Gibbs free binding energies according to eqs 2, 4, and 5, see Figure 2b. Disregarding the sample with lowest total pressure, we find very similar predictions, limiting ΔG°295 K to a range from ca. 4.8 to ca. 7.9 kJ mol−1. A value of ΔG°295 K of 6 kJ mol−1, for example, corresponds to a value of K°eq ° of 0.087 and implies a concentration of 2.3 Torr of the MeCN−HCl complex in the gas mixture with pHCl = 280 Torr and pMeCN = 70 Torr. Ab Initio Benchmarking. Finally, we compare the obtained Gibbs free binding energy range with ab initio electronic structure calculations. First, the structures of HCl,

MeCN, and the MeCN−HCl molecular complex were optimized using the methods described in the Computational Methods section. Hereafter, coupled-cluster-based electronic energy corrections were determined, and the corrected Gibbs free binding energies were obtained according to eq 3. The counterpoise correction terms were calculated using eq 4. Counterpoise corrections less than 0.8 kJ mol−1 are found for the DFT and F12 methods, whereas values around 2.5 kJ mol−1 are found for MP2 and the remaining coupled cluster methods. For this system, however, we find that the counterpoise correctional scheme significantly worsen the agreement with the experimental result. Such an effect has been observed in several other systems,37,38 and we therefore choose to disregard basis set superposition error for the remainder of this study. The counterpoise correction terms are presented in Figure S2. The resulting values for ΔG are summarized in Figure 4. More details are given in Tables S6 and S7.

Figure 4. Summary of the calculated Gibbs free binding energies of the HCl−MeCN complex. Optimized structures and thermal correction terms are calculated using the methods indicated by the color coding. Single-point electronic energies are calculated using the methods indicated at the x-axis. All calculations utilize the AV(T+d)Z basis set, except the CCSD(T)-F12b calculations, here abbreviated F12, which utilize the VTZ-F12 basis set. The experimental range is included for comparison.

Presented in the first column of points in Figure 4, we see that most density functionals significantly overestimate ΔG°295 K. B3LYP performs the worst, predicting a value of 12.8 kJ mol−1, whereas CAM-B3LYP, M06-2X, and wB97XD predict values around 10.8 kJ mol−1. MP2, PW91, and B3LYPD3 perform very well, predicting values of 7.6, 7.4, and 6.3 kJ D

dx.doi.org/10.1021/jp5037537 | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

mol−1, respectively. This is in accordance with previous findings by Du et al.18,19,22 for several hydrogen-bonded clusters, generally finding that pure DFT functionals, and in particular B3LYP, underbind by up to 10 kJ mol−1. Values for ΔG°295 K corrected by CC2 electronic energy calculations according to eq 3 are presented in the second column of data points in Figure 4. The M06-2X- and B3LYPD3-based results predict ΔG°295 K values of 5.7 and 4.1 kJ mol−1, whereas all remaining methods predict values between 7 and 7.5 kJ mol−1. The much-improved performance for most methods is striking but is similar to the findings for the DMS− HCl complex.20 The third row of data points present CCSD electronicenergy-corrected ΔG°295 K calculations. These are seen to be much less accurate than the CC2-corrected values, despite CCSD being a computationally more costly and more advanced method. The fourth and fifth row present data corrected by CCSD(T)-F12b and CCSD(T) electronic energies, respectively. These are seen to improve the CCSD corrected results by ca. 0.6 and ca. 1.9 kJ mol−1, respectively. Also here, these results are consistent with the findings by Du et al.,18,19 whom studied the B3LYP, M06-2X, and wB97XD functionals, finding that coupled cluster corrections to the electronic energy, generally, tend to strengthen the bond by a few kJ mol−1. When corrected by CCSD(T) electronic energies, the M06-2X and B3LYP-D3 based predictions (ΔG°295 K = 8.5 and 7.1 kJ mol−1) are both bordering the experimental range, whereas the PW91based prediction (ΔG°295 K = 11.0 kJ mol−1) is the least accurate. Test calculations using CCSDT/cc-pV(D+d)Z showed that full inclusion of triple excitations into the coupled cluster scheme lowers the energy by less than 0.07 kJ mol−1 compared to CCSD(T) calculations, suggesting that higherlevel coupled-cluster calculations will have a very limited effect on the binding energy predictions. The good performance of coupled cluster corrected M06-2X results for describing noncovalent interactions and hydrogenbonded complexes has been found in several experimental18−20 and theoretical39−42 studies. Also, the much improved performance of B3LYP-D3 compared to conventional B3LYP is welldocumented for systems of relevance to atmospheric new particle formation.41,42 Further, it should be noted that pure PW91 predictions and CC2-corrected B3LYP predictions both fall within the experimental range, similar to the observations for the DMS−HCl complex.20 To the extent that the MeCN− HCl and DMS−HCl clusters are representative for the sulfuric acid−base systems responsible for atmospheric new particle formation, our results support the validity of the predictions made by the groups of Yu (using PW91)43,44 and Vehkamäki (using CC2 and B3LYP).3,45 As apparent from the large scatter of data in the fifth column in Figure 4, significant errors in the thermal correction terms are induced by the DFT- and MP2-based vibrational analysis. To assess the thermal correction terms more accurately, structural optimizations and vibrational analysis were conducted using CCSD and CCSD(T) with the AV(T+d)Z basis set and CCSD(T)-F12b with the VTZ-F12 basis set. The resulting Gibbs free binding energies are shown as black points in Figure 4 and are all seen to be within 1.5 kJ mol−1 of the experimental range. At each level of theory, these results are on average improving the corrected DFT- and MP2-based results (i.e., the colored points) by about 3 kJ mol−1, although it should be noted that the B3LYP-D3-based results are in slightly better accordance with the experimental range. Also, a

consistent improvement from CCSD over CCSD(T)-F12b to CCSD(T) of about 1.2 kJ mol−1 in total is seen with the ΔG°295 K prediction from CCSD(T)/AV(T+d)Z falling within the experimental range. Despite the good agreement between CCSD(T) calculations and experimental results, the harmonic oscillator approximation used for determining vibrational entropy is a well-known source of error.46 This approximation may in principle be avoided using molecular dynamics simulations,47 but this approach is very computationally demanding and is in practice unavailable for systems of atmospheric interest. Vibrational perturbation theory (VPT2)48 has recently been implemented in the Gaussian software package and is capable of providing the anharmonic frequencies required to calculate the partition functions. Although both HCl and MeCN could be treated without problems, we found the MeCN−HCl cluster could not be treated. This may be due to the many linear atoms or the wide and shallow potentials of the lowest lying vibrational modes, but the precise reason remains unclear. A simpler, empirical-based approach to account for vibrational anharmonicity is manual scaling of the vibrational frequencies,46,49,50 requiring no additional electronic structure calculations. Temelso et al.50 used MP2/aug-cc-pVDZ and compared harmonic and VPT2 calculations of water clusters. They suggested to scale the zero point vibrational energy (ZPVE) by 0.976 and the vibrational frequencies by 0.844 and 0.868 when entering the vibrational entropy and enthalpy terms, respectively. For the MeCN−HCl system, this approach lowers ΔZPVE by ca. 0.06 kJ mol−1 and lowers the Gibbs free binding energy by ca. 1.8 kJ mol−1, regardless of ab initio method. We find that this energy lowering originates almost entirely from the scaling of the two lowest frequency modes (between 28 and 35 cm−1, depending on method), emphasizing the importance of calculating the low frequency modes accurately. By this approach, all of the results from coupled-cluster calculations only (black points in Figure 4) are well within the experimental range. This agreement adds to the reliability of both the experimental result and a vibrational scaling approach.



CONCLUSIONS Motivated by improving models of formation and growth of atmospheric aerosols, we have studied the MeCN−HCl molecular complex using infrared spectroscopy and a range of electronic structure methods. The complex was formed in mixtures of gaseous MeCN and HCl and identified via the redshifted H−Cl vibrational mode at 2725 cm−1. The oscillator strength of this mode was calculated using the anharmonic oscillator local mode model. We determined values for the Gibbs free binding energy (ΔG°295 K) of the complex to between 4.8 and 7.9 kJ mol−1 at 295 K. The uncertainties are mainly from the calculated value of the oscillator strength and the assignment of spectral features relating to this value. To test the accuracy of some of the most widely used electronic structure methods, we calculated ΔG°295 K using five density functionals, perturbation theory, and five coupled cluster methods up to CCSDT quality, considering also the D3 dispersion correctional scheme. This study is thus the most comprehensive benchmarking study of the performance of electronic structure methods on hydrogen-bonded clusters to date. The existing data for similar complexes suggest that these findings are consistent, at the very least, for molecular complexes involving HCl, but probably to hydrogen-bonded clusters in general. E

dx.doi.org/10.1021/jp5037537 | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

Schobesberger, S.; Rantala, P.; et al. Direct Observations of Atmospheric Aerosol Nucleation. Science 2013, 339, 943−946. (3) Almeida, J.; Schobesberger, S.; Kürten, A.; Ortega, I. K.; Kupiainen-Mäaẗ tä, O.; Praplan, A. P.; Adamov, A.; Amorim, A.; Bianchi, F.; Breitenlechner, M.; et al. Molecular Understanding of Sulphuric Acid-Amine Particle Nucleation in the Atmosphere. Nature 2013, 502, 359−363. (4) Zhang, R.; Khalizov, A.; Wang, L.; Hu, M.; Xu, W. Nucleation and Growth of Nanoparticles in the Atmosphere. Chem. Rev. 2012, 112, 1957−2011. (5) Kazil, J.; Stier, P.; Zhang, K.; Quaas, J.; Kinne, S.; O’Donnell, D.; Rast, S.; Esch, M.; Ferrachat, S.; Lohmann, U.; Feichter, J. Aerosol Aucleation and its Role for Clouds and Earth’s Radiative Forcing in the Aerosol-Climate Model ECHAM5-HAM. Atmos. Chem. Phys. 2010, 10, 10733−10752. (6) Kaufman, Y. J.; Tanré, D.; Boucher, O. A Satellite View of Aerosols in the Climate System. Nature 2002, 419, 215−223. (7) Solomon, S.; Qin, D.; Manning, M.; Chen, Z.; Marquis, M.; Averyt, K. B.; Tignor, M.; Miller, H. L. IPCC, 2007: Climate Change 2007: The Physical Science Basis. Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change. 2007. (8) Vehkamäki, H. Classical Nucleation Theory in Multicomponent Systems; Springer: Berlin, 2006. (9) Laaksonen, A.; Talanquer, V.; Oxtoby, D. W. Nucleation: Measurements, Theory, and Atmospheric Applications. Annu. Rev. Phys. Chem. 1995, 46, 489−524. (10) Petäjä, T.; Sipilä, M.; Paasonen, P.; Nieminen, T.; Kurtén, T.; Ortega, I. K.; Stratmann, F.; Vehkamäki, H.; Berndt, T.; Kulmala, M. Experimental Observation of Strongly Bound Dimers of Sulfuric Acid: Application to Nucleation in the Atmosphere. Phys. Rev. Lett. 2011, 106, 228302. (11) Froyd, K. D.; Lovejoy, E. R. Experimental Thermodynamics of Cluster Ions Composed of H2SO4 and H2O. 1. Positive Ions. J. Phys. Chem. A 2003, 107, 9800−9811. (12) Curtius, J.; Froyd, K. D.; Lovejoy, E. R. Cluster Ion Thermal Decomposition (I): Experimental Kinetics Study and Ab Initio Calculations for HSO4−(H2SO4)x(HNO3)y. J. Phys. Chem. A 2001, 105, 10867−10873. (13) Jokinen, T.; Sipilä, M.; Junninen, H.; Ehn, M.; Lönn, G.; Hakala, J.; Petäjä, T.; Mauldin, R., III; Kulmala, M.; Worsnop, D. Atmospheric Sulphuric Acid and Neutral Cluster Measurements Using CI-APiTOF. Atmos. Chem. Phys. 2012, 12, 4117−4125. (14) Olenius, T.; Schobesberger, S.; Kupiainen-Maatta, O.; Franchin, A.; Junninen, H.; Ortega, I. K.; Kurtén, T.; Loukonen, V.; Worsnop, D. R.; Kulmala, M.; et al. Comparing Simulated and Experimental Molecular Cluster Distributions. Faraday Discuss. 2013, 165, 75−89. (15) Dom, J. J.; van der Veken, B. J.; Michielsen, B.; Jacobs, S.; Xue, Z.; Hesse, S.; Loritz, H.-M.; Suhm, M. A.; Herrebout, W. A. On the Weakly C−H···π Hydrogen Bonded Complexes of Sevoflurane and Benzene. Phys. Chem. Chem. Phys. 2011, 13, 14142−14152. (16) Hippler, M.; Hesse, S.; Suhm, M. A. Quantum-Chemical Study and FTIR Jet Spectroscopy of CHCl3−NH3 Association in the Gas Phase. Phys. Chem. Chem. Phys. 2010, 12, 13555−13565. (17) Hippler, M. Quantum Chemical Study and Infrared Spectroscopy of Hydrogen-Bonded CHCl3−NH3 in the Gas Phase. J. Chem. Phys. 2007, 127, 084306. (18) Du, L.; Mackeprang, K.; Kjaergaard, H. G. Fundamental and Overtone Vibrational Spectroscopy, Enthalpy of Hydrogen Bond Formation and Equilibrium Constant Determination of the Methanol−Dimethylamine Complex. Phys. Chem. Chem. Phys. 2013, 15, 10194−10206. (19) Du, L.; Lane, J. R.; Kjaergaard, H. G. Identification of the Dimethylamine−Trimethylamine Complex in the Gas Phase. J. Chem. Phys. 2012, 136, 184305. (20) Bork, N.; Du, L.; Kjaergaard, H. G. Identification and Characterization of the HCl−DMS Gas Phase Molecular Complex via Infrared Spectroscopy and Electronic Structure Calculations. J. Phys. Chem. A 2014, 118, 1384−1389.

We have determined the following points of interest: • Most electronic structure methods overestimate ΔG°295 K. • Of the non coupled cluster corrected predictions, B3LYP was the furthest from the experiment and only PW91, MP2, and B3LYP-D3 were within the experimental range. • Most DFT and MP2 results corrected by CC2 electronic energies are within the experimental range. • Among the tested MP2 and DFT methods, M06-2X and B3LYP-D3 provided the most accurate thermal correction terms. • Vibrational scaling factors significantly improve most results, decreasing the Gibbs free binding energy by ca. 2 kJ mol−1. Further, we find that using coupled-cluster calculations only, either CCSD, CCSD(T)-F12b, or CCSD(T), is a viable way for providing high-quality binding energies, albeit only applicable to small systems. After applying vibrational scaling to account for anharmonicity, these result are all well within the experimental range.



ASSOCIATED CONTENT

S Supporting Information *

The following is available as Supporting Information: descriptions and references to all ab initio methods used in this study; details on the anharmonic oscillator local mode model; potential energies and dipole moments used in the local mode model calculated by CCSD(T)/AV(Q+d)Z; Morse parameters for the local mode calculations calculated by CCSD(T)/AV(T+d)Z; XYZ coordinates of the DMS−HCl molecular cluster optimized calculated by CCSD(T)/AV(Q +d)Z; harmonic frequencies of HCl, MeCN, and MeCN−HCl calculated by CCSD(T)/AV(T+d)Z; a representative deconvolution of the FTIR spectra; a figure of the counterpoise correction energies; tables of electronic energies; and Gibbs free energies and entropies of the clustering reaction. This material is available free of charge via the Internet at http:// pubs.acs.org/.



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS We acknowledge the Villum foundation and the Danish Council for Independent Research−Natural Sciences for funding and the CSC-IT Center for Science Ltd. and the Danish Center for Scientific Computing for allocation of computational resources. T.K. thanks the Academy of Finland for funding (grant no. 266388).



REFERENCES

(1) Kirkby, J.; Curtius, J.; Almeida, J.; Dunne, E.; Duplissy, J.; Ehrhart, S.; Franchin, A.; Gagne, S.; Ickes, L.; Kurten, A.; et al. Role of Sulphuric Acid, Ammonia and Galactic Cosmic Rays in Atmospheric Aerosol Nucleation. Nature 2011, 476, 429−435. (2) Kulmala, M.; Kontkanen, J.; Junninen, H.; Lehtipalo, K.; Manninen, H. E.; Nieminen, T.; Petä j ä , T.; Sipilä , M.; F

dx.doi.org/10.1021/jp5037537 | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

(21) Schriver, L.; Schriver, A.; Perchard, J.-P. Infrared MatrixIsolation Study of the Molecular Complexes Between Acetonitrile and Hydrogen Halides. Evidence for Reactivity Differences Within the Halide Series. J. Chem. Soc., Faraday Trans. 2 1985, 81, 1407−1425. (22) Du, L.; Kjaergaard, H. G. Fourier Transform Infrared Spectroscopy and Theoretical Study of Dimethylamine Dimer in the Gas Phase. J. Phys. Chem. A 2011, 115, 12097−12104. (23) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (24) Boys, S. F.; Bernardi, F. d. The Calculation of Small Molecular Interactions by the Differences of Separate Total Energies. Some Procedures With Reduced Errors. Mol. Phys. 1970, 19, 553−566. (25) van Duijneveldt, F. B.; van Duijneveldt-van de Rijdt, J. G.; van Lenthe, J. H. State of the Art in Counterpoise Theory. Chem. Rev. 1994, 94, 1873−1885. (26) Atkins, P. W.; de Paula, J. Physical Chemistry; Oxford University Press, 2006. (27) Chung, S.; Hippler, M. Infrared Spectroscopy of HydrogenBonded CHCl3-SO2 in the Gas Phase. J. Chem. Phys. 2006, 124, 214316. (28) Howard, D. L.; Kjaergaard, H. G. Vapor Phase Near Infrared Spectroscopy of the Hydrogen Bonded Methanol−Trimethylamine Complex. J. Phys. Chem. A 2006, 110, 9597−9601. (29) Howard, D. L.; Kjaergaard, H. G. Hydrogen Bonding to Divalent Sulfur. Phys. Chem. Chem. Phys. 2008, 10, 4113−4118. (30) Lane, J. R.; Kjaergaard, H. G.; Plath, K. L.; Vaida, V. Overtone Spectroscopy of Sulfonic Acid Derivatives. J. Phys. Chem. A 2007, 111, 5434−5440. (31) Lane, J. R.; Kjaergaard, H. G. XH-Stretching Overtone Transitions Calculated Using Explicitly Correlated Coupled Cluster Methods. J. Chem. Phys. 2010, 132, 174304. (32) Miller, B. J.; Du, L.; Steel, T. J.; Paul, A. J.; Södergren, A. H.; Lane, J. R.; Henry, B. R.; Kjaergaard, H. G. Absolute Intensities of NHStretching Transitions in Dimethylamine and Pyrrole. J. Phys. Chem. A 2011, 116, 290−296. (33) Kjaergaard, H. G.; Garden, A. L.; Chaban, G. M.; Gerber, R. B.; Matthews, D. A.; Stanton, J. F. Calculation of Vibrational Transition Frequencies and Intensities in Water Dimer: Comparison of Different Vibrational Approaches. J. Phys. Chem. A 2008, 112, 4324−4335. (34) Millen, D.; Zabicky, J. 565. Hydrogen Bonding in Gaseous Mixtures. Part V. Infrared Spectra of Amine−Alcohol Systems. J. Chem. Soc. (Resumed) 1965, 3080−3085. (35) Legon, A.; Campbell, E.; Flygare, W. The Rotational Spectrum and Molecular Properties of a Hydrogen-Bonded Complex Formed Between Hydrogen Cyanide and Hydrogen Chloride. J. Chem. Phys. 1982, 76, 2267. (36) Rothman, L.; Gordon, I.; Barbe, A.; Benner, D.; Bernath, P.; Birk, M.; Boudon, V.; Brown, L.; Campargue, A.; Champion, J.-P.; et al. The HITRAN 2008 Molecular Spectroscopic Database. J. Quant. Spectrosc. Radiat. Transfer 2009, 110, 533−572. (37) Galano, A.; Alvarez-Idaboy, J. R. A New Approach to Counterpoise Correction to BSSE. J. Comput. Chem. 2006, 27, 1203−1210. (38) de Lange, K. M.; Lane, J. R. Explicit Correlation and Intermolecular Interactions: Investigating Carbon Dioxide Complexes with the CCSD(T)-F12 Method. J. Chem. Phys. 2011, 134, 034301. (39) Zhao, Y.; Truhlar, D. G. The M06 Suite of Density Functionals for Main Group Thermochemistry, Thermochemical Kinetics, Noncovalent Interactions, Excited States, and Transition Elements: Two New Functionals and Systematic Testing of Four M06-Class Functionals and 12 Other Functionals. Theor. Chem. Acc. 2008, 120, 215−241. (40) Elm, J.; Bilde, M.; Mikkelsen, K. V. Assessment of Density Functional Theory in Predicting Structures and Free Energies of Reaction of Atmospheric Prenucleation Clusters. J. Chem. Theory Comput. 2012, 8, 2071−2077.

(41) Elm, J.; Bilde, M.; Mikkelsen, K. V. Assessment of Binding Energies of Atmospherically Relevant Clusters. Phys. Chem. Chem. Phys. 2013, 15, 16442−16445. (42) Leverentz, H. R.; Siepmann, J. I.; Truhlar, D. G.; Loukonen, V.; Vehkamäki, H. Energetics of Atmospherically Implicated Clusters Made of Sulfuric Acid, Ammonia, and Dimethyl Amine. J. Phys. Chem. A 2013, 117, 3819−3825. (43) Herb, J.; Xu, Y.; Yu, F.; Nadykto, A. Large Hydrogen-Bonded Pre-nucleation (HSO4−)(H2SO4)m(H2O)k and (HSO4−)(NH3)(H2SO4)m(H2O)k Clusters in the Earth’s Atmosphere. J. Phys. Chem. A 2012, 117, 133−152. (44) Nadykto, A. B.; Herb, J.; Yu, F.; Xu, Y. Enhancement in the Production of Nucleating Clusters Due to Dimethylamine and Large Uncertainties in the Thermochemistry of Amine-Enhanced Nucleation. Chem. Phys. Lett. 2014, DOI: 10.1016/j.cplett.2014.03.036. (45) Ortega, I.; Kupiainen, O.; Kurtén, T.; Olenius, T.; Wilkman, O.; McGrath, M.; Loukonen, V.; Vehkamäki, H. From Quantum Chemical Formation Free Energies to Evaporation Rates. Atmos. Chem. Phys. 2012, 12, 225−235. (46) Temelso, B.; Shields, G. C. The Role of Anharmonicity in Hydrogen-Bonded Systems: The Case of Water Clusters. J. Chem. Theory Comput. 2011, 7, 2804−2817. (47) Lin, S.-T.; Maiti, P. K.; Goddard, W. A., III Two-Phase Thermodynamic Model for Efficient and Accurate Absolute Entropy of Water From Molecular Dynamics Simulations. J. Phys. Chem. B 2010, 114, 8191−8198. (48) Bégué, D.; Baraille, I.; Garrain, P.; Dargelos, A.; Tassaing, T. Calculation of IR Frequencies and Intensities in Electrical and Mechanical Anharmonicity Approximations: Application to Small Water Clusters. J. Chem. Phys. 2010, 133, 034102. (49) Dunn, M. E.; Evans, T. M.; Kirschner, K. N.; Shields, G. C. Prediction of Accurate Anharmonic Experimental Vibrational Frequencies for Water Clusters,(H2O)n, n = 2−5. J. Phys. Chem. A 2006, 110, 303−309. (50) Temelso, B.; Archer, K. A.; Shields, G. C. Benchmark Structures and Binding Energies of Small Water Clusters with Anharmonicity Corrections. J. Phys. Chem. A 2011, 115, 12034−12046.

G

dx.doi.org/10.1021/jp5037537 | J. Phys. Chem. A XXXX, XXX, XXX−XXX