Ab Initio Structural and Vibrational Investigation of Sulfuric Acid

Article pubs.acs.org/JPCA

Ab Initio Structural and Vibrational Investigation of Sulfuric Acid Monohydrate Lauri Partanen, Vesa Han̈ ninen,* and Lauri Halonen Laboratory of Physical Chemistry, Department of Chemistry, University of Helsinki, P.O. Box 55 (A.I. Virtasen aukio 1), FIN-00014 Finland S Supporting Information *

ABSTRACT: We employ ab initio methods to find stable geometries and to calculate potential energy surfaces and vibrational wavenumbers for sulfuric acid monohydrate. Geometry optimizations are carried out with the explicitly correlated coupled-cluster approach that includes single, double, and perturbative triple excitations (CCSD(T)-F12a) with a valence double-ζ basis set (VDZ-F12). Four different stable geometries are found, and the two lowest are within 0.41 kJ mol−1 (or 34 cm−1) of each other. Vibrational harmonic wavenumbers are calculated at both the density-fitted local spin component scaled second-order Møller−Plesset perturbation theory (DF-SCSLMP2) with the aug-cc-pV(T+d)Z basis set and the CCSD-F12/VDZ-F12 level. Water O−H stretching vibrations and two highly anharmonic largeamplitude motions connecting the three lowest potential energy minima are considered by limiting the dimensionality of the corresponding potential energy surfaces to small two- or three-dimensional subspaces that contain only strongly coupled vibrational degrees of freedom. In these anharmonic domains, the vibrational problem is solved variationally using potential energy surfaces calculated at the CCSD(T)-F12a/VDZ-F12 level.

■

INTRODUCTION In the atmosphere, sulfuric acid is one of the key components at the first stages of cloud formation.1,2 A significant part of new particle formation in the atmosphere is thought to begin with the complexation of sulfuric acid with water, ammonia, and other species.3−5 Most of the sulfuric acid molecules in the atmosphere are contained in hydrates, which has a large impact on the nucleation rates. Therefore, in order to devise reliable nucleation models, it is crucial that the equilibrium constants for sulfuric acid hydration reactions be known accurately.6 Both ab initio and density functional theory calculations have been previously used to obtain structural and thermodynamical parameters for sulfuric acid hydrates.7−11 Unfortunately, the variance of the results between different methods is high, and because the size of the clusters can be large, the level of the calculations has not generally been good. Most importantly, there is an additional need for improvement because the equilibrium constants are sensitive to the accuracy of the calculations. Recently, major developments have been achieved in the different explicitly correlated methods that take into account the singularities in the potential energy at points where two electrons coincide.12−16 Specifically, the explicitly correlated coupled-cluster singles and doubles method with perturbative triples CCSD(T)-F12, has made it possible to obtain accurate electronic energies and equilibrium structures with smaller basis sets than before. This has substantially reduced © 2012 American Chemical Society

computational time when compared with standard CCSD(T) methods.16,17 In the Møller−Plesset second-order perturbation theory (MP2) calculations, a careful use of localized orbitals can be used to reduce the computational scaling with respect to the molecular size down to linear, and by construction, these local methods are virtually free of the basis set superposition error.18−21 Additional reductions in computational time of about an order of magnitude can be achieved with densityfitting approximations, which reduce the steep fourth-order scaling with respect to the basis set size per atom.22 Furthermore, to account for the biased description of singlet and triplet electron pairs already present at the Hartree−Fock (HF) level, the spin component scaling (SCS) methods separately scale the singlet and triplet MP2, reducing the errors usually present in this method.23−25 For atmospheric species like the sulfuric acid monohydrate, the electronic and rotational parts of the molecular partition function 26 can be evaluated satisfactorily with existing methods. The hardest problems lie in the computation of the vibrational energies. Previously, the effects of using different computational methods for the harmonic wavenumbers have been studied, 27 and even the effect of anharmonic frequencies has been considered Received: November 1, 2011 Revised: January 19, 2012 Published: January 19, 2012 2867

dx.doi.org/10.1021/jp210489f | J. Phys. Chem. A 2012, 116, 2867−2879

The Journal of Physical Chemistry A

Article

with methods such as the correlation-corrected vibrational self-consitent field approach (CC-VSCF).9,28,29 However, all of these treatments fail to asses the nature of floppy large-amplitude motions that extend over several local potential energy minima present in H2SO4·H 2O. In this study, anharmonic effects were considered by limiting the dimensionality of the potential energy surface (PES) to small two- or three-dimensional subspaces that contained a few strongly coupled vibrational degrees of freedom. In the next section, we introduce the vibrational model used for the anharmonic domain calculations and give a brief overview of the computational details. The results of our quantum chemical and variational vibrational calculations are given in the third section. We use the CCSD(T)-F12a method with the valence double-ζ VDZ-F12 basis set30 for the calculation of electronic energies and equilibrium structures. For vibrational harmonic wavenumbers, we use either the same VDZ-F12 basis combined with a simplified CCSD-F12a method, where the perturbative triples are missing, or a density-fitted local spin component scaled MP2 (DF-SCS-LMP2) calculation with the Dunning augmented correlation-consistent polarized valence triple-ζ basis set with an additional high exponent function added for the third-period atoms (aug-cc-pV(T+d)Z).31−33 To ascertain the validity of our results, we compare these with earlier computational studies and experimental results for equilibrium geometries7,34−38 and harmonic frequencies.27,36,37,39,40 The anharmonic domain PES calculations are done at the CCSD(T)-F12a/VDZ-F12 level. For these anharmonic domains, the vibrational problem is solved variationally. Because the domains are completely uncoupled to the rest of the vibrational degrees of freedom, great care is taken to confirm the reliability of results with experimental data37,40−43 and other anharmonic calculations when experimental data are not available.28 Conclusions are given in the final section.

■

configuration and the high-frequency vibrations that have the largest effect on the zero-point energies, specifically the two water stretches, and the water bend to which these are coupled. In the high-frequency case, the three-dimensional PES for water in the monohydrate complex was calculated by keeping the coordinates of sulfuric acid and intermolecular modes constant and varying the internal coordinates of water. Following the approach adopted in ref 44 after the PES had been computed, the vibrational problem of the water molecule was solved variationally. For the variational calculation, the nuclear Hamiltonian was expressed in terms of the curvilinear internal coordinates. With the omission of angular momentum components in the Hamiltonian, it takes the form44−46 Ĥ = T̂ + V̂ =−

ℏ2 2

(q , q ) ⎞ 2 ⎤ ∂g i j ⎟ ∂ (q , q ) ∂ ⎥ +g i j + V̂ ′(q) + V̂ (q) ∂qi ⎟⎠ ∂qj ∂qi∂qj ⎥ ij ⎣⎝ ⎦

⎡⎛

∑ ⎢⎢⎜⎜

(1)

where the q vector consists of the curvilinear internal coordinates and g(qi,qj) terms are the elements of the contravariant metric tensor and have the form45,46 (q , q ) g i j =

Nn

∑ α

1 (∇α qi) ·(∇α qj) mα

(2)

where mα quantities are the masses of the nuclei, α is a summation index over the nuclei, and Nn is the number of the nuclei. The pseudopotential term V̂ ′(q) is independent of the momentum operators47,48 and is of quantum mechanical origin. Its effect on the energies is usually small. Finally, V̂ (q) is the potential energy operator given in the curvilinear internal coordinates. The weight function in the volume element of integration is 1. Curvilinear internal coordinates give an accurate representation of the PES and make the surface parameters independent of the different isotopic species within the Born−Oppenheimer approximation.47 The drawback of the curvilinear internal coordinates compared to rectilinear normal coordinates is such that the kinetic energy operator becomes complicated in curvilinear coordinates, as can be seen from eq 1. The potential energy operator was expressed in terms of the Morse variables z b and z f for the stretches and the curvilinear internal displacement coordinate θ for the bend, as in ref 47. The Morse variable z b is defined as z b = 1 − e −abrb, where r b = R b−Rb,e and ab is the Morse parameter.49 The quantity R b is the instantaneous bond length of the hydrogen-bonded OH bond in the water molecule, and Rb,e is the equilibrium bond length. The variable z f is defined analogously for the bond between the oxygen and free hydrogen atom. The bond angle displacement coordinate is θ = ϕ − ϕe, where ϕ is an instantaneous bond angle and ϕ e is the equilibrium value. In the isolated water molecule, Rb,e = Rf,e and ab = af for the Morse parameters. However, due to the asymmetry of the sulfuric acid environment, the OH bonds of H2O in the sulfuric acid monohydrate possess different equilibrium bond lengths

METHODS

Vibrational Model. For vibrational motions with a single deep minimum, the harmonic approximation is sometimes used for the few lowest levels. Problems arise when the vibrational quantum number increases because the equally spaced harmonic energy level ladder fails to describe the ever increasing density of vibrational states as the dissociation limit becomes closer. Also, in systems with multiple potential energy minima, the splitting of states caused by tunneling effects causes the harmonic approximation to fail. As the size of the system increases, the dimensionality of the PES increases by three with each added atom. Unfortunately, the complete variational anharmonic treatment requires calculations of vast areas of the PES, and this already makes it impossible for systems such as the sulfuric acid molecule and sulfuric acid monohydrate complex, which have 15 and 24 vibrational degrees of freedom, respectively. For small systems such as a single water molecule with its three vibrational degrees of freedom, complete anharmonic treatment is doable by a full variational calculation. In this work, we divide the vibrational problem of the sulfuric acid monohydrate into smaller uncoupled subsystems that can be treated with the variational approach. From now on, this will be known as the anharmonic domain approximation. The anharmonic domains were used here for those large-amplitude modes that could result in a change of 2868



Article

and Morse parameters. To account for this, the potential energy operator adopted takes the more general form 7

V (z b , z f , θ ) =

7

∑ k=2

Tk(b)z bk +

∑

Tk(f)z fk +

k=2

fr r bf z bz f a ba f

⎛ fr r ⎞ 1 fr r r + ⎜⎜ b2 b f + b f ⎟⎟z b2z f 2 ⎝ ab af a ba f ⎠ +

+ + + + + +

⎛ fr r ⎞ 1 ⎜ fr brf rf 1 bf⎟ z bz f2 + fθθ θ 2 + ⎜ ⎟ 2 2 ⎝ a ba f 2 a ba f ⎠ fr θ 1 1 fθθθ θ 3 + fθθθθ θ 4 + b z bθ 6 24 ab fr θθ fr θθ fr θ f z f θ + b z bθ 2 + f z f θ 2 2a b 2a f af 1 ( fr r θ + a b fr θ)z b2θ b 2a b2 b b 1 ( fr r θ + a f fr θ)z f2θ 2 f f f 2a f 1 ( fr r θθ + ab fr θθ)z b2θ 2 b 4a b2 b b 1 ( fr r θθ + a f fr θθ )z f2θ 2 f 4a f2 f f (3)

Figure 1. The four sulfuric acid monohydrate conformers starting from the upper left corner labeled as g1, g2, g3, and g4 in the order of increasing energy. The angles and bond lengths specified refer to the values in Table 1.

angles in eq 5 were constrained to the equilibrium values. The potential energy term V(τ,ϑ) is of the form 6

ℏ2 ∂ 2 ℏ2 ∂ 2 Ĥ = − − + V (τ , ϑ) 2Iτ ∂τ 2 2Iϑ ∂ϑ2

(A m cos mτ + Bm sin mτ )

m=1 8

+

∑

Ckϑk + fτϑ (τ − τe)(ϑ − ϑe)

k=2

+ fττϑ (τ − τe)2 (ϑ − ϑe) + fτϑϑ (τ − τe)(ϑ − ϑe)2 + fττϑϑ (τ − τe)2 (ϑ − ϑe)2 + fτϑϑϑ (τ − τe)(ϑ − ϑe)3 + fτττϑ (τ − τe)3 (ϑ − ϑe)

(6)

where τe and ϑe are equilibrium values for the torsional and wagging angles, respectively. Clearly, a Fourier series50 up to the sixth order was used to describe the one-dimensional part in the torsion angle τ, while an ordinary eighth-order polynomial was used for ϑ. Coupling terms were included up to fourth order. The wave function of the system was approximated by a harmonic oscillator basis set for the wagging motion and a free rotor basis set for the torsional motion. Computational Details. Most of the geometry optimizations were done in three stages. First, DF-SCS-LMP2 calculations with the Dunning aug-cc-pV(D+d)Z basis set31−33 were performed starting from a structural guess. Structures were then refined with the same method using a higher-quality aug-cc-pV(T+d)Z basis set. All DF-SCS-LMP2 geometry optimizations made use of the spin component scaled MP2 energy functional. The localization was done using the Pipek−Mezey method.51 The problems associated with the localization tails arising from the diffuse basis functions52 were avoided by ignoring the most diffuse basis function of each

(4)

where τ is the torsional angle and ϑ is the wagging angle. The moments of inertia Iτ and Iϑ are defined by 2 Ii = mHr H, i sin θ H, i

∑

V (τ , ϑ)/hc = A 0 +

This equation is of similar form as the potential energy operator in ref 44. As can be seen, a polynomial up to the seventh order was used in the Morse variables, while a fourth-order polynomial was needed for θ. Coupling terms between the stretches and bend were included up to the fourth order and up to the third order between the two stretches. The eigenvalues were obtained variationally using the Hamiltonian in eq 1. A basis set of harmonic oscillator wave functions was used for the bend. Morse oscillator eigenfunctions49 were employed for the stretches. The analytical matrix elements resulting from the Hamiltonian adopted in the calculation are found elsewhere.44,47 Two of the large-amplitude motions of the sulfuric acid monohydrate were included in the anharmonic treatment. These were chosen so that the transformation from one conformer to another was described naturally. The motions were the torsional rotation of the non-hydrogen-bonded hydrogen in the sulfuric acid molecule (transformation between g1 and g3) and the wagging motion of the free hydrogen in the water molecule (transformation between g1 and g2). See Figure 1 for details. The vibrational energy levels of the problem described above were computed variationally with an approximate Hamiltonian given by

(5)

where mH is the mass of the hydrogen atom, rH,i is the bond length of the OH bond ri, θH,i is the angle between the bond and the rotational axis, and i = τ or ϑ. The bond lengths and 2869



Article

For the sulfuric acid monohydrate, a total of four different stable geometries were found and are given in Figure 1. These structures have already been reported in earlier studies.7−11 The optimized geometries for the different conformers are given in Table 1. In the first three conformers, water and sulfuric acid form a six-membered ring with water acting both as a donor and an acceptor. As can be seen from Figure 1, the most stable g1 geometry can be converted into g3 by a change of the dH7O6S3O2 dihedral angle. It can be thought of as resulting from the complexation of a g2 sulfuric acid conformer (see Figure 2) with water. The geometries g1 and g2 are connected by a rotation of the water plane with respect to the sulfuric acid species. The fourth conformer of the complex differs radically from the first three because while the six-membered ring is still present, the water molecule acts solely as a donor and the sulfuric acid molecule as an acceptor in this geometry. Additional geometries were also investigated, for example, a structure with water as in g2 in Figure 1 and sulfuric acid in the g2 conformation of Figure 2, but none of these structures were stable. Looking at Table 1, compared to the experimental geometry by Fiacco et al.,34 both the CCSD(T)-F12a/VDZ-F12 and DFSCS-LMP2/AV(T+d)Z methods yield structures similar to those from the B3LYP calculation by Bandy and Ianni.8 As is seen in Figure 1, the g1, g2, and g3 conformers of the sulfuric acid monohydrate are held together by a strong hydrogen bond, where the sulfuric acid molecule acts as a donor and the water molecule as an acceptor. In the g1 conformer, this bond length is 1.677 Å at the CCSD(T)-F12a level of theory (1 Å = 10−10 m). From the length of the hydrogen bond, it is possible to estimate the bond enthalpy,62 which for this bond is 31 kJ mol−1. In the second hydrogen bond (r10), water acts as a donor to the doubly bonded oxygen atom of sulfuric acid. The bond has a CCSD(T)-F12a bond length of 2.138 Å in the g1 conformer, corresponding to a bond enthalpy of around 15 kJ mol−1. While it is not as strong as the r7 bond, this bond is nevertheless important in providing ridigity to the monohydrate molecule. Curiously, r10 is shortest in the g3 geometry with the CCSD(T)-F12a value of 2.120 Å, where r7, on the other hand, possesses the largest value. For all three geometries, largest differences between the DFSCS-LMP2/AV(T+d)Z bond lengths and the CCSD(T)-F12a/ VDZ-F12 ones are found in the two hydrogen bonds. On average, when compared with the CCSD(T)-F12a results, the DF-SCS-LMP2 method overestimates the hydrogen bonds by 0.05 Å for r7 and 0.282 Å for r10, while the average difference for the rest of the bonds is only 0.005 Å. We think that the reason for this is the relatively large decrease in total energy when moving from DF-SCS-LMP2 to CCSD(T)-F12a, which translates to tighter, that is, more energy containing, bonds. This effect is strongest in loose bonds, which means that the hydrogen bonds will be most affected. In the fourth observed geometry with the water molecule located on the top of the sulfuric acid molecule, the two hydrogen bonds are significantly weaker than those in all the other geometries, which can be seen from the respective longer bond lengths of 2.396 and 2.386 Å at the CCSD(T)-F12a level of theory. The energies of the geometry-optimized structures for the two sulfuric acid conformers and the four monohydrate ones are shown in Table 2. The absolute energies are listed in hartrees (1 hartree = 4.36 × 10−19 J), and the energies relative to the g1 sulfuric acid geometry are given in kJ mol−1. Due to

angular momentum type for the aug-cc-pV(D+d)Z basis set and two of the most diffuse basis functions of each angular momentum type for the aug-cc-pV(T+d)Z basis set in the localization. The local domains were determined automatically using the procedure of Boughton and Pulay53 with the completeness criterion set to 0.98 for the aug-cc-pV(D+d)Z basis set and 0.985 for the aug-cc-pV(T+d)Z basis set. Finally, all intermolecular pairs were treated as strong in the sulfuric acid monohydrate complex, and the convergence threshold for the HF calculation was set to 10−14 because of the sensitivity of the LMP2 gradients to the accuracy of the self-consistent field (SCF) convergence. After the DF-SCS-LMP2 calculations, the final geometries and energy differences were obtained from a CCSD(T)-F12a calculation with the VDZ-F12 basis of Peterson et al.30 Of the two simplifying approximations to the standard CCSD(T)-F12 method present in the MOLPRO package,54 the aapproximation was chosen because, even though it slightly overestimates the correlation energies, it has been shown to give better results for the aug-cc-pVDZ and aug-cc-pVTZ basis sets than its counterpart CCSD(T)-F12b.16 The VDZF12 basis is optimized for use in F12 calculations and has a similar size as the aug-cc-pVDZ basis set. The CCSD-F12a/ VDZ-F12 geometries were a byproduct of vibrational harmonic frequency calculations done at this level. From this point onward, we shall denote the aug-cc-pV(D+d)Z basis by AV(D+d)Z and the aug-cc-pV(T+d)Z by AV(T+d)Z. The counterpoise correction was not made in this study because both the local and F12 methods greatly decrease basis set superposition errors.16,55 The aug-cc-pVDZ/MP2FIT auxiliary basis sets of Weigend et al.56,57 were used in the density fitting of the CCSD-F12a and CCSD(T)-F12a calculations, except for the calculation of the exchange and Fock operators, where a VDZ/ JKFIT basis set was used. The VDZ-F12/optri auxiliary basis set of Yusef and Peterson58 was employed for the resolution of identity parts. This choice of the basis has been found to work well for other hydrogen-bonded complexes.59,60 All geometry optimizations, single-point calculations, and frequency calculations were carried out with the MOLPRO suit of programs.54

■

RESULTS AND DISCUSSION Geometry Optimizations and Energies. The geometries for both the sulfuric acid and water molecules are given in the Supporting Information, where a more detailed interpretation of the geometry optimizations for these molecules is found. The two stable conformers of sulfuric acid reported by several other authors7,11,61 were detected in this study as well (see Figure 2). The CCSD(T)-F12a/VDZ-F12 and

Figure 2. The two sulfuric acid conformers labeled as g1 and g2 in the order of increasing energy. The bond lengths and angles specified refer to Table 1 in the Supporting Information.

DF-SCS-LMP2/AV(T+d)Z methods are shown to give excellent results for both of these geometries and the water molecule. 2870



Article

Table 1. Equilibrium Geometries of Different Sulfuric Acid Monohydrate Conformersa DF-SCS-LMP2/AV(T+d)Z r1 r2 r3 r4 r5 r6 r7 r8 r9 r10 a1 a2 a3 a4 a5 a6 a7 a8 dO6S3O2H1 dH7O6S3O2 dO8H7O6S3 a

g1

g2

g3

0.969 1.594 1.423 1.429 1.566 0.992 1.730 0.967 0.962 2.293 108.3 108.6 122.8 109.3 102.8 108.7 165.2 105.8 82.4 85.0 12.2

0.969 1.600 1.423 1.427 1.567 0.992 1.715 0.965 0.963 2.680 108.1 108.8 122.9 109.8 102.4 108.6 171.0 105.0 72.8 78.8 29.6

0.968 1.593 1.416 1.438 1.567 0.989 1.742 0.967 0.962 2.276 108.5 105.7 122.7 108.5 102.9 109.0 165.7 106.1 −87.2 96.0 6.2

CCSD(T)-F12a/VDZ-F12 g4 0.970 1.587 1.423 2.529 2.526 0.963

108.7 102.4 109.1 123.4 103.2

g1

g2

g3

0.966 1.582 1.416 1.424 1.554 0.994 1.677 0.966 0.959 2.138 108.4 108.4 122.4 109.1 103.0 108.9 164.3 106.0 84.1 89.4 14.5

0.966 1.585 1.416 1.424 1.553 0.994 1.677 0.966 0.959 2.145 108.2 108.2 122.4 109.1 102.9 108.9 164.2 105.9 87.3 89.0 19.9

0.965 1.580 1.409 1.433 1.556 0.992 1.681 0.966 0.959 2.120 109.0 105.7 122.1 108.3 103.0 109.3 165.6 106.4 −85.8 102.8 12.9

g4 0.967 1.575 1.417 2.396 2.386 0.960

109.0 102.6 108.9 123.1 103.2

exptl.

other calc.

g134

g1b

0.95 1.578 1.410 1.464 1.567 1.04 1.645 0.98 0.98 2.050 108.5

1.611 1.430 1.439 0.999 1.682 0.968 0.962 2.207

123.3 106.7 101.8 108.6 165.2 107

108.7 106.5

13.8

The angles are in degrees, and the bond lengths are in Å. bB3LYP/6-311++G(2d,2p) results are by Bandy and Ianni.8

Table 2. Energies and Relative Energies for the Geometry-Optimized Structures for the Different Sulfuric Acid Monohydrate Conformersa DF-SCS-LMP2/AV(D+d)Z H2SO4 g1 g2 H2SO4·H2O g1 g2 g3 g4 a

DF-SCS-LMP2/AV(T+d)Z

CCSD(T)-F12a/VDZ

E

ΔE

E

ΔE

E

ΔE

−699.09952 −699.09750

5.30

−699.43739 −699.43538

5.28

−699.59581 −699.59381

5.25

−775.37640 −775.37571 −775.37447 −775.36395

1.79 5.05 32.67

−775.78408 −775.78317 −775.78215 −775.77003

2.40 5.09 36.90

−775.97921 −775.97905 −775.97732 −775.96493

0.41 4.97 37.50

CCSD-F12a/VDZ E −699.55526

−775.93093 −775.92901

The E values are reported in hartrees, and the ΔE values are in kJ mol−1.

be verified by comparing relative Boltzmann populations between the geometries g3 and g4 at the average tropospheric temperature of 268 K.63 Vibrational Harmonic Wavenumbers. Harmonic wavenumbers were calculated for water and all different conformers of sulfuric acid and the sulfuric acid monohydrate at the DFSCS-LMP2/AV(T+d)Z level of theory. For the g1 geometry, the wavenumbers were also calculated at the CCSD-F12a/ VDZ-F12 level of theory. Unfortunately, the harmonic wavenumber calculations with the CCSD(T)-F12a/VDZ-F12 method proved to be computationally too demanding for the H2SO4·H2O system. Water and Sulfuric Acid Molecules. As a first application, the harmonic wavenumbers calculated were used to distinguish between true local minima and saddle points. This was possible because at saddle points, the harmonic wavenumber corresponding to the reaction coordinate is imaginary. Furthermore, the T1 diagnostic was used to estimate whether the single determinantal wave function of HF methods gave an accurate enough representation of the system so that the application of correlation methods using just a single determinant reference

the simplicity of the sulfuric acid system, the local methods give almost the same energy difference as the CCSD(T)-F12a/ VDZ-F12 treatment, but our results are somewhat larger than the energy difference of 3.28 kJ mol−1 predicted by Demaison et al.61 at the MP2/AV(Q+d)Z level of theory. Of the monohydrate conformers, the g1 geometry was found to be most stable, but it is only 0.41 kJ mol−1 (34 cm−1) below the g2 conformer in energy at the CCSD(T)-F12a/VDZ-F12 level of theory. The difference between g1 and g3 is 4.96 kJ mol−1. Both of the different DF-SCS-LMP2 methods give the g3 − g1 energy difference correctly, while larger discrepancies exist for the g2 − g1 difference. The results are not discouraging as the difference between the DF-SCS-LMP2/ AV(T+d)Z and the CCSD(T)-F12a/VDZ-F12 calculations is only 2 kJ mol−1. Although not shown in Table 2, the CCSD-F12a energy difference between these two states is 5.04 kJ mol−1. All of the different correlation methods show that the g4 geometry is much higher in energy than any of the other conformations. This means that even though the state exists, its population at atmospheric temperatures is so small that it has little effect on thermodynamic properties. This observation can 2871



Article

Table 3. Harmonic Wavenumbers for Sulfuric Acid and Water in cm−l this study H2SO4

SO2 rock OSO bend SO2 wag S−OH symmetric stretch S−OH asymmetric stretch symmetric SOH bend asymmetric SOH bend SO2 symmetric stretch SO2 asymmetric stretch asymmetric OH stretch symmetric OH stretch H2O HOH bend symmetric OH stretch asymmetric OH stretch

exptl.

gla

g2a

glb

glc

258 338 381 445 502 548 558 834 881 1173 1188 1241 1492 3773 3777

108 348 371 428 505 547 559 834 884 1159 1184 1238 1493 3784 3789

246 321 391 453 518 569 581 883 933 1180 1204 1275 1525 3832 3836

237 311 381 443 506 553 565 853 906 1167 1190 1247 1498 3785 3789

1663 3870 3976

1652 3834 3942

1645 3811 3934

gl65

551 834 891 1157 1220 1465 3609 ref 37 1648 3832 3943

other calc. gl39

568 550 834 883 1138 1159 1223 1450 3610

gld

g1e

380 379 459 489 525 519 480 779 821 1118 1128 1183 1434 3500 3590 ref 36

221 321 341 413 463 499 513 745 801 1138 1153 1160 1426 3627 3632 ref 66

1648 3832 3942

1649 3833 3944

a

DF-SCS-LMP2/AV(T+d)Z. bCCSD-F12a/VDZ-F12. cCCSD(T)-F12a/VDZ-F12. dCC-VSCF/MP2-TZP anharmonic calculation byMiller et al.28 PW91/TZP calculation by Al Natsheh et al.7

e

was justified. The T1 diagnostic is defined with the help of the coupled-cluster singles amplitude vector t1 by64 T1 =

|t1| Nelec

The OH stretches are overestimated with both the CCSD(T)-F12a and DF-SCS-LMP2 methods by more than 150 cm−1. As can be seen from the results by Natsheh et al.,7 DFT methods suffer no such drawback but tend to show worse perfomance for the middle-range frequencies. Miller et al. suggested that this success results from the cancellation of errors between the overly soft DFT frequencies and the increases in frequency caused by the harmonic approximation.28 A method that gives good predictions in the midinfrared range is preferred because the middle-frequency modes have a higher impact on the vibrational partition function than the high-frequency modes and because the high-frequency OH stretches are often less coupled to the other motions67 and thus easier to treat anharmonically. Because of anharmonicity, the harmonic wavenumbers are not fully comparable to the experimental ones. However, as long as the anharmonic effects can be assumed to be small, experimental data can be used for the selection of the superior method. To investigate the applicability of the harmonic approximation, Table 3 also lists anharmonic fundamental wavenumbers calculated at the CC-VSCF/MP2-TZP level by Miller et al.28 For the high-frequency OH stretches, these are much closer to the experimental ones than any of the ab initio results. Surprisingly, for some of the middle-range wavenumbers, these show large deviations from the experimental data. The problems with the anharmonic method are associated with the accurate prediction of the S−(OH)2 stretches and the SO2 symmetric stretch. In the case of the water molecule, the equilibrium geometry has a T1 value of 0.008. We see from Table 3 that the CCSD(T)-F12a/VDZ-F12 method produces results close to the experimental harmonic wavenumbers.37 With this method, the difference between the computed and experimental harmonic wavenumber is at most 4 cm−1. The DF-SCSLMP2 method performs in a similar manner as it did in the sulfuric acid case. As before, the DF-SCS-LMP2 values are in

(7)

where Nelec is the number of electrons. A single determinantal reference is adequate if T1 < 0.020. The harmonic wavenumbers for the stable sulfuric acid geometries g1 and g2 (see Figure 2) are given in Table 3. The T1 diagnostics for these geometries were 0.015 for each conformer. We see from the DF-SCS-LMP2/AV(T+d)Z calculations that the harmonic wavenumbers for the different conformers are usually within 10 cm−1 of each other for all but the lowest torsional mode, where the difference is 150 cm−1. Comparing the DF-SCS-LMP2/AV(T+d)Z calculations with the CCSD(T)-F12a/VDZ-F12 ones and experimental values,39,65 one finds that the DF-SCS-LMP2 values underestimate the experimental ones below 1173 cm−1 and overestimate the rest. Because the harmonic approximation tends to overestimate the observed anharmonic wavenumbers, this behavior implies that below 1173 cm−1, the DF-SCS-LMP2/AV(T+d)Z method tends to underestimate the true wavenumbers. The CCSD(T)-F12a values show the expected behavior as these tend to overestimate the experimental values everywhere. In addition to the harmonic approximation, the reason for the overestimation might be that the CCSD(T) method has a habit of overestimating the vibrational harmonic wavenumbers. Normally, this is compensated for by the basis set incompleteness error, which is greatly diminished in the CCSD(T)-F12 method.59 The rms deviations from the experimental values are similar for both the DF-SCS-LMP2 and CCSD(T)-F12a methods. The CCSD-F12a calculations overestimate the CCSD(T)-F12a values by more than 9 cm−1 everywhere and are thus of lower quality. 2872



Article

Table 4. Harmonic wavenumbers for sulfuric acid monohydrate in cm−1 this work

SO2 rock OSO bend SO2 wag S−OH symmetric stretch S−OH asymmetric stretch free SOH bend SO2 symmetric stretch bound SOH bend SO2 asymmetric stretch HOH bend on H2O bound OH stretch on H2SO4 symmetric OH stretch on H2O free OH stretch on H2SO4 asymmetric OH stretch on H2O a

b

exptl.

g1a

g2a

g3a

g1b

48 120 191 211 271 314 388 428 465 530 557 565 766 844 921 1179 1230 1359 1484 1645 3329 3768 3781 3903

14 53 114 215 278 315 406 442 373 523 553 562 867 830 916 1184 1237 1332 1482 1656 3324 3788 3776 3904

30 124 195 223 235 318 381 425 459 535 553 566 745 845 920 1175 1218 1351 1486 1637 3368 3771 3787 3905

1 128 178 212 246 318 395 442 484 547 577 588 797 893 969 1190 1261 1382 1512 1658 3332 3810 3839 3944

c

g143

other calc. g140

554

839 929

834 887 1205 1449 1600

2722 3582

3640

3697

3745

g1c

g1d

g1e

104 180 529 331 397 529 600 486 540 536 514 476 787 834 868 1153 1172 1336 1436 1573 2973 3593 3628 3738

49 119 209 229 244 486 353 389 462 327 497 528 709 792 808 1090 1124 1281 1401 1595 2970 3580 3609 3720

48 167 212 258 284 339 389 411 491 509 527 635 761 856 917 1131 1140 1343 1467 1607 2771 3508 3638 3730

DF-SCS-LMP2/AV(T+d)Z. CCSD-F12/VDZ-F12. CC-VSCF/MP2-TZP anharmonic calculation by Miller et al.28 calculation by Nadykto et al.27 ePW91/TZP calculation by Al Natsheh et al.7

d

BLYP/aug-cc-pVTZ

these directly to estimate the bond enthalpies of the hydrogen bonds.62,68 In the case of the strong hydrogen bond, the estimated bond enthalpy is ΔHr7 = 27 kJ mol−1, which is a little lower than the value predicted from the bond length. The largest changes are observed in the vibrational modes such as the SOH bends that have changed their character as a result of the complexations. Comparison of harmonic wavenumbers with the experimental data is problematic in the sulfuric acid monohydrate case because the large-amplitude motions have significant anharmonicity and all experimental data have been obtained from matrix isolation studies using argon matrixes.40,43 The system is difficult to study in a precise manner both experimentally and computationally. There is a lot more variation in the experimental values than in the monomers. For the high-frequency motions, the DF-SCS-LMP2/ AV(T+d)Z method again overshoots the experimental values by over 100 cm−1, while the different DFT methods of Nadykto et al.27 and Natsheh et al.7 give values at most only some tens of reciprocal centimeters away from the experimental ones. Because the harmonic approximation overestimates the vibrational wavenumbers, the results show that the DFT wavenumbers in themselves are not a good representation of the true values. The anharmonicity of the low-wavenumber large-amplitude intermolecular modes is seen in the CC-VSCF/MP2-TZP results of Miller et al.,28 where differences of at least 50 cm−1 with the harmonic values are common. As before, the anharmonic calculation yields far better results than any of the other harmonic ab initio methods for the OH stretches, but in the middle range of wavenumbers, the agreement is about as good as that with the DF-SCS-LMP2 method.

general smaller than the CCSD(T)-F12a ones, and the CCSDF12a values are always larger than the CCSD(T)-F12a ones. Compared with the water monomer results by Kim et al.,36 we see that their CCSD(T) calculation using the nonstandard basis set with the (13,8,4,2/8,4,2) contraction gives results that agree even better with experimental values than the CCSD(T)F12a/VDZ-F12 calculation. The same is true for the harmonic wavenumbers obtained from a CVRQD PES by Barletta et al.66 Sulfuric Acid Monohydrate. Harmonic wavenumbers for the different stable sulfuric acid monohydrate conformers are given in Table 4. Because of its low impact on the thermodynamic properties, the structurally different g4 conformer has been omitted from the table. For the three presented structures, the T1 diagnostic values were all 0.014, well under the 0.02 limit. By comparing the different geometries, we see that the g1 and g3 conformers have vibrational wavenumbers mostly within 10 cm−1 of each other with some exceptions, such as the stretching motion of the bound sulfuric acid OH group. For the g1 and g2 conformers, the differences are in general larger as already many of the low-frequency modes differ by tens of reciprocal centimeters. It should be noted that for large systems such as the sulfuric acid monohydrate, it becomes difficult to assign exact desciptions to some of the low-frequency modes; therefore, some caution is recommended when comparing our results for the unnamed modes with those of previous studies. The greatest difference between the harmonic spectrum of the bound complex and the spectra of its free constituents is the red shift of 445 cm−1 in the hydrogen-bonded OH group. These kinds of large red shifts have been known from experiments to accompany strong hydrogen bonds. Therefore, it is possible to use 2873



Article

Anharmonic Vibrational Wavenumbers. The two factors that have the largest effect on the vibrational partition function of sulfuric acid monohydrate are the large-amplitude and low-frequency vibrations, for which the contributions of the excited vibrational states to the partition function are the largest, and the high-frequency vibrations, which have the largest effect on the zero-point energy. We have reserved anharmonic vibrational treatment only to two subsystems for the sulfuric acid monohydrate. The first one is a threedimensional space spanned by the internal coordinates of the two stretches and the bend on the doubly hydrogen-bonded water molecule. The second one is a two-dimensional space consisting of the wagging motion of the free water hydrogen, which connects the conformers g1 and g2, and the free HOSO torsion, which connects the conformers g1 and g3. Three-Dimensional Vibrational H 2 O Problem in H2SO4·H2O and Water. To obtain the parameters present in eq 3, single-point calculations were performed around the g1 equilibrium structure by varying one or two of the three water coordinates. The three two-dimensional surfaces were obtained by varying the bond lengths from −0.375 to 0.3375 Å and the angle from −20 to 20° around the equilibrium geometry. Each two-dimensional surface consisted of 169 points. Additional points were calculated for the one-dimensional cuts of the PESs. In the end, these consisted of 25 points for the stretches ranging from −0.375 to 0.4875 Å and 20 points for the bend with the same −20 to 20° range as in the two-dimensional surfaces. After the surface calculations were complete, the parameters of the complex were obtained by the least-squares optimization of the potential energy operator parameters of eq 3 to the surface point data. In order to obtain the best possible fit near the equilibrium geometry, this was done by first finding the parameters for the one-dimensional parts of the surface in eq 3. In the case of the hydrogen-bonded hydrogen atom, the best fit was obtained with a function containing the terms Tn(b)z bn with n = 2, 3, ..., 6, while a polynomial containing the terms Tn(f)z fn with n = 2, 3, ..., 7 gave the best fit for the freely vibrating hydrogen. For both fits, the Morse parameter a was first estimated by least-squares fitting a fourth-order polynomial to the one-dimensional data. The bending PES cut proved to be well-described by a fourth-order polynomial in the displacement coordinate θ. Because the Taylor series of eq 3 is evaluated at the equilibrium, the first-order terms make no contribution to it. With all one-dimensional parameters known, the two-dimensional parameters were obtained by fitting eq 3 to the two-dimensional PESs with the one-dimensional parameters held fixed. All leastsquares fittings were done with Mathematica.69 The results for the conformer g1 are in Table 5. In the isolated water monomer case, the PESs were calculated as described above, and a sixth-order polynomial with all of the different terms was used for the OH stretches, while a fourth-order polynomial was used for the bend. Because of the symmetry of the monomer system, some of the parameters, such as f rbrfrf and f rbrbrf, possess the same values. The final parameters are shown in Table 5. In the variational energy level calculation, for both the Morse oscillator and the harmonic oscillator wave functions, the maximum value of the quantum numbers was set to 8. To further limit the basis in the calculation, it was required that the sum of the hydrogen-bonded OH stretch quantum number nb and the free OH stretch quantum number nf was less than or equal to 8. After the results were obtained, basis set convergence was checked by performing the same variational calculation

Table 5. PES Parameters in Equation 3 for Water in the g1 Conformer of the Sulfuric Acid Monohydrate Complex H2O(c) and for the Isolated Water Molecule H2O(i) H2O(c) −1

H2O(i) −2

H2O(c)

H2O(i)

ab (Å )

1.8099

1.9950

f rbrf (aJ Å )

−0.0896

−0.1244

T2(b)(aJ)

1.2215

1.0634

f rbrfrf (aJ Å−3)

0.0896

0.1357

T3(b)(aJ)

−0.3532

−0.1687

f rbrbrf (aJ Å−3)

0.1365

T4(b)(aJ) T5(b)(aJ) T6(b)(aJ)

0.1175 0.0145 0.0188

0.0836 0.0121 0.0054

af (Å−1)

1.8852

f rbrbθ (aJ Å−2)

−0.2324

−0.1612

T2(f)(aJ) T3(f)(aJ) T4(f)(aJ) T5(f)(aJ) T6(f)(aJ) T7(f)(aJ)

1.1872

f rbθθ (aJ Å−1)

−0.4472

−0.3884

f rbrbθθ (aJ Å−2)

−0.2573

−0.0269

−0.2652

fθθ (aJ) fθθθ (aJ) fθθθθ (aJ) f rbθ (aJ Å−1)

0.6999 −0.6302 −0.6505 0.1374

0.1171

f rfθ (aJ Å−1)

0.0125

f rfrfθ (aJ Å−2)

−0.1281

−0.0124

f rfθθ (aJ Å−1)

−0.3509

−0.0150

f rfrfθθ (aJ Å−2)

−0.0339

0.7061 −0.7323 −0.6385 0.2664

0.2429

with the maximum values set to 9. The states were considered converged if the wavenumber change upon the increase of the basis set size was less than 1 cm−1. It should be noted that the changes were below 0.1 cm−1 for the seven lowest states. Some of the states calculated from these parameters are shown in Table 6 for the monohydrate. In Table 6, rf denotes the free r9 bond, and rb refers to the bond between the hydrogen-bonded H9 and O8 atoms (see Figure 1). The state of the system is given as a linear combination of different |i⟩f|j⟩b|k⟩ basis functions, where i is the quantum number for the free OH stretching vibrational Morse oscillator, j is the quantum number for the bound OH stretching vibrational Morse oscsillator, and k is the quantum number for the harmonic oscillator describing the bending motion. The percentage behind the term describes its weight in the state. The effect of the hydrogen bond can be seen in the percentage shift from the 50/50 distribution between the |0⟩f|i⟩b|0⟩ and |i⟩f|0⟩b|0⟩ states to one where one state always dominates the other. Therefore, instead of a symmetric or an asymmetric stretch, it is better to talk about the stretches of certain bonds, that is, to use the local description of molecular vibrations. As the quantum numbers increase, the states become more and more mixed, and, for example, the third bending overtone already possesses 11% of the |0⟩f|0⟩b|5⟩ character and 81% of the |0⟩f|0⟩b|4⟩ character. The fundamental wavenumbers of water are again given in Table 7. To estimate the effect of the argon matrix on the experimental sulfuric acid monohydrate values, both the most recent gas-phase and argon matrix values for water are shown. In the isolated case, the gas-phase value for the bend is about 5 cm−1 larger than the matrix value, while for both OH stretches, the gas-phase values are larger by 20 cm−1. Thus, the gas-phase values for the bound water can be expected to be somewhat larger than the experimental matrix values given in Table 7. Even when taking this into account, the difference between the experimental40,43 and anharmonic values is at most around 20 cm−1, which has little effect on the vibrational partition function. As can be seen from Table 4, the experimental value reported for the bonded OH stretch by Givan et al.40 is almost 60 cm−1 larger than the value by Rozenberg and Loewenschuss43 shown in Table 7. Thus, especially if the effect of the argon 2874



Article

The two final columns in Table 7 report the shifts in wavenumbers upon complexation for both the anharmonic and experimental values. As expected, the largest red shift occurs in the OH bond where the hydrogen is hydrogen-bonded to the O5 oxygen in the sulfuric acid. This red shift is a little smaller than the red shift encountered in the hydrogen-bonded OH stretching of the water dimer, which, from the results of Salmi et al.,45,46 can be predicted to be 97 cm−1 at the CCSD(T)/ AVQZ level of theory. Accordingly, this second hydrogen bond (including water hydrogen) formed in sulfuric acid monohydrate is somewhat weaker than the one formed in the water dimer. This is supported by the fact that the bond in the water dimer is also a little shorter with a length of 1.95 Å, compared with 2.14 Å of the second hydrogen bond in the monohydrate. If the effect of the argon matrix on the wavenumbers is approximately constant, the shifts between the anharmonic and experimental values should be comparable. With regard to the shifts in wavenumbers, the best description is achieved for the free OH. This is logical as this degree of freedom is the least coupled to the intermolecular modes and the vibrational modes in sulfuric acid. Because of the hydrogen bond, the strongest couplings are most likely found with the hydrogen-bonded OH stretch. This is supported by experiment and our computational data shown in Table 7, where the shifts are larger for the bonded OH stretch compared to those of the free one. Because the HOH bend has the effect of bending the hydrogen-bonded hydrogen toward the sulfuric acid molecule, rather strong couplings can be expected for this degree of freedom as well, a notion that is supported by the 16 cm−1 difference between the anharmonic and experimental shifts. Two-Dimensional Large-Amplitude Problem in H2SO4·H2O. The two-dimensional PES consists of the torsional rotation of the H1 hydrogen around the S3−O2 axis and the wagging motion of the H10 hydrogen atom. These two motions connect the conformers g1, g2, and g3, and therefore, their anharmonic treatment is of fundamental importance if accurate values for the thermodynamic properties are to be obtained. The two-dimensional PES was determined with the CCSD(T)F12a/VDZ-F12 level of theory. It consisted of 117 points. The torsional angle τ varied from −150° to 180° and the wagging angle θ from −90° to 90°, with the average increment being 30° for both angles. As in the three-dimensional H2O case, the parameters in eq 6 were obtained first by performing two separate least-squares optimizations of the coefficients of the eighth-order polynomial and Fourier series using onedimensional cuts of the surface with the other variable held at its equilibrium value. For the torsional motion, the cut contained 36 points, and for the wagging angle, the number of points was 21. Because the geometry optimizations at the CCSD(T)-F12a level of theory are demanding computationally, the structure could not

Table 6. Vibrational states and wavenumbers of water in the g1 conformer of sulfuric acid monohydrate. All calculations have been done at the CCSD(T)-F12a/VDZ-F12 level wavenumber/ cm−1

state |0⟩f|0⟩b|1⟩ |0⟩f|0⟩b|2⟩ |0⟩f|1⟩b|0⟩ |1⟩f|0⟩b|0⟩ |0⟩f|0⟩b|3⟩ |0⟩f|1⟩b|1⟩ |1⟩f|0⟩b|1⟩ |0⟩f|0⟩b|4⟩ |0⟩f|1⟩b|2⟩ |1⟩f|0⟩b|2⟩ |0⟩f|2⟩b|0⟩ |2⟩f|0⟩b|0⟩ |1⟩f|1⟩b|0⟩ |0⟩f|2⟩b|1⟩ |2⟩f|0⟩b|1⟩ |1⟩f|1⟩b|1⟩ |0⟩f|3⟩b|0⟩ |3⟩f|0⟩b|0⟩ |1⟩f|2⟩b|0⟩ |2⟩f|1⟩b|0⟩ |0⟩f|3⟩b|1⟩ |3⟩f|0⟩b|1⟩ |1⟩f|2⟩b|1⟩ |2⟩f|1⟩b|1⟩ |0⟩f|4⟩b|0⟩ |4⟩f|0⟩b|0⟩ |1⟩f|3⟩b|0⟩ |3⟩f|1⟩b|0⟩ |2⟩f|2⟩b|0⟩ |0⟩f|4⟩b|1⟩

1589 3146 3579 3727 4669 5139 5294 6157 6665 6829 7018 7214 7358 8547 8758 8899 10300 10606 10713 10913 11800 12130 12225 12429 13430 13854 13976 14183 14378 14901

(85%) + |1⟩f|0⟩b|0⟩ (13%) (85%) + |0⟩f|1⟩b|0⟩ (13%) (84%) (84%) (81%) (80%) (82%)

+ + + + +

|1⟩f|0⟩b|1⟩ |0⟩f|1⟩b|1⟩ |0⟩f|0⟩b|5⟩ |1⟩f|0⟩b|2⟩ |0⟩f|1⟩b|2⟩

(11%) (12%) (11%) (9%) (11%)

(67%) + |1⟩f|1⟩b|0⟩ (24%) (67%) + |2⟩f|0⟩b|0⟩ (28%) (61%) + |1⟩f|1⟩b|1⟩ (26%) (62%) + |2⟩f|0⟩b|1⟩ (30%) (70%) + |2⟩f|1⟩b|0⟩ (14%) (64%) + |3⟩f|0⟩b|0⟩ (16%) + |2⟩f|1⟩b|0⟩ (10%) (70%) + |1⟩f|2⟩b|0⟩ (20%) (59%) (59%) (67%) (79%) (77%) (70%) (51%) (56%) (71%)

+ + + +

|2⟩f|1⟩b|1⟩ |3⟩f|0⟩b|1⟩ |1⟩f|2⟩b|1⟩ |0⟩f|5⟩b|0⟩

(14%) + |1⟩f|2⟩b|1⟩ (9%) (18%) (18%) (9%)

+ + + +

|2⟩f|2⟩b|0⟩ |2⟩f|2⟩b|0⟩ |3⟩f|1⟩b|0⟩ |0⟩f|5⟩b|1⟩

(12%) (22%) + |1⟩f|3⟩b|0⟩ (11%) (32%) (12%)

matrix is to lower the wavenumbers, our computational results support the more recent data in ref 43. On the other hand, the quality of our computed results justifies the approximation that the water molecule can be treated independently while holding the rest of the system rigid. On the basis of the fundamental wavenumbers and the first vibrational overtone of the water monomer, the variational calculation tends to produce more accurate values for the bending motions than those for the stretches, as was the case in all calculations for the water dimer by Salmi et al.45,46 This propably has to do with the form of the polynomial chosen for the stretches. A more balanced treatment could possibly be achieved by introducing higher-order coupling terms between the stretches and the bend in eq 3.

Table 7. Fundamental Wavenumbers for Free Water and Water Bound in the g1 Conformer of the Sulfuric Acid Monohydrate Complex Compared with Experimental Valuesa isolated H2O anharm. HOH bend sym. stretch (bonded) asym. stretch (free)

1594.5 3652.4 3761.5

exptl.b 1589.1 3638.0 3734.3

bound H2O

exptl.c

other calc.d

1594.7 3657.1 3755.9

1595.5 3651.5 3756.4

anharm. 1589.2 3578.8 3726.9

exptl. e

1599.7 3581.9f 3696.6f

Δν̃a

Δν̃e

5.3 73.6 34.6

−10.6 56.1 37.7

All values are in cm−1. The shifts in wavenumber for the anharmonic Δν̃a and experimental Δν̃e values are defined as the free water fundamental wavenumber minus the fundamental wavenumber in the hydrate. bEngdahl and Nelander42 in an Ar matrix. cTennyson et al.70 in the gas phase. d CCSD(T)/AVQZ calculation by Salmi et al.45,46 eGivan et al.40 in an Ar matrix. fRozenberg and Loewenschuss43 in an Ar matrix. a

2875



Article

for the vibrational quantum number was 12. The torsional quantum number possessed integer values ranging from −11 to 11. In Figure 4, the states below 1000 cm−1 are shown along with the two one-dimensional cuts from Figure 3 obtained by varying only one of the variables and keeping the other fixed at its equilibrium value. All of the states up to 1500 cm−1 are also given in Table 9, along with the one-dimensional results where all of the couplings have been ignored for both the torsional and the wagging motion calculated with the CCSD(T)-F12a method. Figure 4 is interesting in several respects. We see that the barrier height between the g1 and g2 geometries is only 384 cm−1, which corresponds to 4.6 kJ mol−1. On the other hand, the barrier height between the g1 and g3 geometries is much larger, 996 cm−1 or 11.9 kJ mol−1. We also see that the zero-point energy is 297 cm−1, which is almost double the value of 160 cm−1, predicted by the harmonic approximation. As shown in Table 9, the transition lowest in wavenumber is 34 cm−1, which deviates from the harmonic prediction of 48 cm−1 for the vibrational degree of freedom that best describes this motion (see Table 4). From Figure 4, we see that this transition arises from the splitting of levels caused by quantum mechanical tunneling. Looking at the next vibrational overtone, a harmonic transition would occur at 96 cm−1, whereas the anharmonic picture predicts this value to be 237 cm−1. This total failure of the harmonic approximation can be attributed to the multistationary nature of the PES. According to our results, the anharmonic treatment leads to a radical change in the distribution of states in the low -wavenumber region. Moving from the harmonic approximation to the anharmonic treatment, there is a flush of states from the low frequencies to higher ones, which, again, is most likely the result of the multistationary nature of the PES. This is in stark contrast with the commonly held belief, based on semirigid systems and stretching vibrations, that adding anharmonic corrections leads to a lowering of the wavenumbers. The change in the density of states at low wavenumbers observed in this work will have a large impact on the vibrational partition function and thus on the entropy term of the Gibbs free energy. It is interesting to note that the fundamental anharmonic transition for the water wagging motion at a wavenumber of 34 cm−1 is also different from the fundamental anharmonic transition of 104.3 cm−1 predicted by Miller et al.28 with the CC-VSCF/ MP2-TZP level of theory. Similar kinds of large discrepancies are obtained for the torsion as the values are 259.4 and 396.5 cm−1, respectively. The underlying reason for this disagreement probably arises because the CC-VSCF/MP2-TZP calculation had to be performed for the g1 conformation only and, thus, had not accounted for the presence of the other conformers. This would, for example, explain the large difference between the wagging fundamentals as the first transition arises from the tunneling effect that is ultimately the result of the multistationary nature of the wagging PES. Thus, our results would indicate that the ab initio VSCF method, while showing good performance for the simple OH stretches, is not useful for the prediction of intermolecular modes, at least when there are several minima present. The tunneling effect is also visible in the splittings of the higher excited states. As expected, the density of states increases with increasing wavenumber. In Figure 4, it can be seen that in the g1 structure, the global minima on the one-dimensional wagging and torsional curves are not exactly on the same level. This happens because the sixth-order Fourier series has trouble predicting the sharp curve

be allowed to relax at each point. Instead, for the one-dimensional wagging calculations, the geometry was frozen to an average between the two minimum conformers g1 and g2. According to test calculations, while this approximation raised energies of both the g1 and the g2 states somewhat, it introduced only a minimal error in the energy difference between the two states. For the torsional motion, a linear mapping was introduced to stepwise change the geometry from the average one used in the wagging calculation to the g3 one. Several different forms for the PES were tried, but the one in eq 6 was chosen because it gave the best results for the onedimensional fits. Table 8 shows the parameter values obtained from the potential surface fitting. The optimized PES is given in Figure 3. The three identified minima can be distinguished with Table 8. PES Parameters in cm−1 for the Two-Dimensional Large-Amplitude Problem in H2SO4·H2O parameter

value

parameter

value

A0 A1 A2 A3 A4 A5 A6 B1 B2 B3 B4 B5 B6 ϑe

1056.67 −90.70 458.43 −71.39 5.59 4.27 0.14 −218.65 20.60 −15.26 4.09 −0.48 0.07 0.977

C2 C3 C4 C5 C6 C7 C8 fτϑ fττϑ fτϑϑ fττϑϑ fτττϑ fτϑϑϑ τe

−834.58 28.35 533.63 −11.12 −71.04 2.21 7.57 −6.16 12.10 0.52 10.66 4.08 6.44 −1.674

Figure 3. The fitted PES for the torsional rotation of H1 and the wagging motion of H10. The τ and ϑ axes are in radians.

the lowest g1 conformation in the front of the picture and the g2 and g3 conformers corresponding to the minima on the right- and the left-hand sides, respectively. The matrix elements were calculated via numerical integration. The energy eigenvalues could be obtained by diagonalizing the Hamiltonian matrix. The maximum number 2876



Article

Figure 4. The 16 lowest states of the variational calculation. The blue curve is the wagging angle PES cut, and the purple one refers to the torsional PES cut with the ϑ angle held fixed. The horizontal axis is in radians.

The inclusion of the coupling terms has large impact only on some of the states, such as the one at 791.5 cm−1 in the twodimensional picture.

Table 9. Vibrational States for the One- and TwoDimensional Treatment of the Large-Amplitude Motionsa CCSD(T)-F12/ VDZ-F12

a

2-dim

2-dim continued

1-dim

1-dim continued

wagging

torsion

297 332 535 557 598 627 707 792 796 800 841 890 939 959 982 999 1006 1035 1060 1062 1122 1130

1168 1178 1184 1189 1229 1245 1245 1248 1272 1301 1331 1340 1346 1372 1391 1403 1410 1429 1436 1442 1452 1484

301 333 539 560 592 629 702 734 794 798 801 826 888 934 940 966 983 1001 1030 1032 1033 1060

1118 1122 1150 1172 1178 1190 1202 1210 1239 1242 1262 1270 1294 1302 1329 1357 1371 1384 1403 1411 1416 1434

156 188 395 485 656 839 1046 1267

144 403 545 638 778 844 962 1022 1114 1214 1278 1432 1451

■

CONCLUSION In this study, the CCSD(T)-F12a method is used to optimize the geometries of the sulfuric acid monohydrate complex and its individual monomer components and to calculate the equilibrium energies as well as PESs. This method is more accurate than those used in previous works for this molecule and has been shown to give results comparable with calculations done at the CCSD(T)/ AVQZ level of theory. The geometry optimizations confirmed the previously found ground-state structures, but we also found that there exists a low-lying stable configuration only 0.41 kJ mol−1 above the most stable equilibrium structure for the sulfuric acid monohydrate given by g2 in Figure 1. While this geometry has been reported by other studies,7,28 its presence is yet to be accounted for in the calculation of the thermodynamic properties, despite the low-energy difference between the two geometries. The harmonic frequencies used in this study were calculated with the DF-SCS-LMP2/AV(T+d)Z and CCSD-F12a/VDZF12 methods. Due to the neglect of anharmonicity, the ab initio wavenumbers usually show large discrepancies from experimental values with the high-frequency OH stretches, but the DF-SCS-LMP2/AV(T+d)Z method showed good agreement in the middle range of frequencies. For the large-amplitude motions, large discrepancies between the methods used and the anharmonic calculations were again detected. For the most important sources of anharmonicity, we made use of the anharmonic domain approximation, where we isolated a few of the degrees of freedom and solved the complete vibrational problem variationally for these while holding the rest of the molecule fixed to some reference geometry. The domains chosen were the two-dimensional space spanned by the wagging angle of the free hydrogen in the hydrogen-bound water molecule and the HOSO torsional angle of the free OH group in the sulfuric acid molecule, together with the three-dimensional space spanned by the internal coordinates of the bound water molecule. Results for

All values are in cm−1.

in the immidiate vicinity of the minimum. The 8.1 cm−1 error arising from this effect was corrected in the final results by subtracting it from the eighth-order polynomial of the torsional motion before the calculation of the matrix elements. From Table 9, we see that the one-dimensional CCSD(T)F12a calculations predict a zero-point energy that is only 3.5 cm−1 from the two-dimensional one. In addition, up to around 750 cm−1, the states are all less than 10 cm−1 apart from the corresponding ones in the two-dimensional calculation. 2877



■

the two-dimensional surface showed radical deviation from the harmonic case as the zero-point energy almost doubled and in general the number of low-lying states was greatly decreased, which would cause large changes in the vibrational partition function. This is in stark contrast with the commonly held belief that adding anharmonic corrections leads to a lowering of the wavenumbers. Our results indicate that for a system with many energetically close minima, it is essential to employ an anharmonic treatment that takes all of the minima into account. For large systems, this might pose a problem, but according to our results, the zero-point energies and a few of the lowest-energy states can be obtained just from one-dimensional calculations with enough accuracy to be used for the calculation of thermodynamic properties. This claim is supported by the fact that the introduction of additional couplings between the different degrees of freedom does not seem to add new low-lying states to the system. The vibrational wavenumbers obtained from the three-dimensional calculation on complexed and free water molecules showed good agreement with experimental results, though the lack of couplings to the other vibrational degrees of freedom was evident in the HOH bend and the OH stretch of the bonded hydrogen in the complexed water molecule. To account for the use of matrix isolation techniques in the experimental values,40,43 the effects of the Ar matrix were estimated by comparing the wavenumber shifts upon complexation for the different fundamental wavenumbers, which were all found to be below 20 cm−1. Because the vibrational partition function is less dependent on the high-wavenumber vibrations, this accuracy is fine. According to our results, approaches akin to the anharmonic domain approximation provide a simple and efficient way to obtain large numbers of fundamental and overtone anharmonic states, with enough accuracy to obtain the thermodynamic properties. It should be noted, though, that to obtain more than just the first few vibrational states accurately, the domains must include all motions with strong couplings to each other. The vibrational approach used in this study together with our electronic energy calculations provides a systematic way through which more accurate thermodynamic properties can be obtained for systems such as the sulfuric acid monohydrate complex. Our future research aims at extending the anharmonic domain approximation to include the two stable conformers in sulfuric acid and finding a suitably simple anharmonic treatment for the rest of the vibrational degrees of freedom. As shown by the results of Miller et al., large anharmonic effects are also present in the remaining harmonically treated large-amplitude motions.28 This problem could be solved by either using some kind of anharmonic corrections to the calculated DF-SCS-LMP2/AV(T+d)Z values or by methods such as the CC-VSCF. At least with the latter, the problem is such that it might be difficult to calculate sufficiently many vibrational overtones28 for the large-amplitude motions, which are essential to obtain an accurate presentation of the vibrational partition function. After a suitable anharmonic treatment is found, we intend to calculate the thermodynamic properties for the sulfuric acid monohydrate. We did not wish to report any preliminary calculations for the thermodynamic properties in this study, but our test calculations show that the inclusion of the other geometries will have a significant effect on the equilibrium constant.

■

Article

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors would like to thank E. Sälli, T. Salmi, T. Kúrten, H. Kjaergaard, and A. Császár for instructive discussions. This project was financially supported by the Academy of Finland (the Centre of Excellence in CMS). We also thank CSC Ltd for computer time.

■

REFERENCES

(1) Finlayson-Pitts, B. J.; Pitts Jr., J. N. Chemistry of the Upper and Lower Atmosphere; Academic Press: New York, 2000. (2) vanLoon, G. W.; Duffy, S. J. Environmental Chemistry, 2nd ed.; Oxford University Press: New York, 2005. (3) Yu, F.; Turco, R. P. Geophys. Res. Lett. 2000, 27, 883−886. (4) Korhonen, P.; Kulmala, M.; Laaksonen, A.; Viisanen, Y.; McGraw, R.; Seinfeld, J. J. Geophys Res. 1999, 104, 26349−26353. (5) Zhang, R.; Suh, I.; Zhao, J.; Zhang, D.; Fortner, E. C.; Tie, X.; Molina, L. T. Science 2004, 304, 1487−1490. (6) Noppel, M.; Vehkamäki, H.; Kulmala, M. J. Chem. Phys. 2002, 116, 218−228. (7) Al Natsheh, A.; Nadykto, A. B.; Mikkelsen, K. V.; Yu, F.; Ruuskanen, J. J. Phys. Chem. A 2004, 108, 8914−8929. (8) Bandy, A. R.; Ianni, J. C. J. Phys. Chem. A 1998, 102, 6533−6539. (9) Kúrten, T.; Sundberg, M. R.; Vehkamäki, H.; Noppel, M.; Blomqvist, J.; Kulmala, M. J. Phys. Chem. A 2006, 110, 7178−7188. (10) Arstila, H.; Laasonen, K.; Laaksonen, A. J. Chem. Phys. 1998, 108, 1031−1039. (11) Re, S.; Osamura, Y.; Morokuma, K. J. Phys. Chem. A 1999, 103, 3535−3547. (12) Kutzelnigg, W.; Klopper, W. J. Chem. Phys. 1991, 94, 1985− 2001. (13) Tew, D. P.; Klopper, W. J. Chem. Phys. 2005, 123, 074101. (14) May, A. J.; Valeev, E.; Polly, R.; Manby, F. R. Phys. Chem. Chem. Phys. 2005, 7, 2710−2713. (15) Ten-no, S.; Manby, F. R. J. Chem. Phys. 2003, 119, 5358−5363. (16) Knizia, G.; Adler, T. B.; Werner, H.-J. J. Chem. Phys. 2009, 130, 054104. (17) Adler, T. B.; Knizia, G.; Werner, H.-J. J. Chem. Phys. 2007, 127, 221106. (18) Saebø, S.; Pulay, P. Annu. Rev. Phys. Chem. 1993, 44, 213−236. (19) Schütz, M.; Hetzer, G.; Werner, H.-J. J. Chem. Phys. 1999, 111, 5691−5705. (20) Hetzer, G.; Schütz, M.; Stoll, H.; Werner, H.-J. J. Chem. Phys. 2000, 113, 9443−9455. (21) Hetzer, G.; Pulay, P.; Werner, H.-J. Chem. Phys. Lett. 1998, 290, 143−149. (22) Werner, H.-J.; Manby, F. R.; Knowles, P. J. J. Chem. Phys. 2003, 118, 8149−8160. (23) Gerenkamp, M.; Grimme, S. Chem. Phys. Lett. 2004, 392, 229− 235. (24) Grimme, S. J. Chem. Phys. 2003, 118, 9095−9102. (25) Distasio, R. A. Jr.; Head-Gordon, M. Mol. Phys. 2007, 105, 1073−1083. (26) McQuarrie, D. A. Statistical Mechanics; Harper and Row: New York, 1976. (27) Nadykto, A. B.; Du, H.; Yu, F. Vib. Spectrosc. 2007, 44, 286−296. (28) Miller, Y.; Chaban, G. M.; Gerber, R. B. J. Phys. Chem. A 2005, 109, 6565−6574. (29) Kúrten, T.; Noppel, M.; Vehkamäki, H.; Salonen, M.; Kulmala, M. Boreal Environ. Res. 2007, 12, 431−453. (30) Peterson, K. A.; Adler, T. B.; Werner, H.-J. J. Chem. Phys. 2008, 128, 084102.

ASSOCIATED CONTENT

S Supporting Information *

The geometries for both the sulfuric acid and water molecules. This material is available free of charge via the Internet at http://pubs.acs.org. 2878



Article

(31) Dunning, T. H. Jr. J. Chem. Phys. 1989, 90, 1007−1023. (32) Dunning, T. H. Jr.; Peterson, K. A.; Wilson, A. K. J. Chem. Phys. 2001, 114, 9244−9253. (33) Kendall, R. A.; Dunning, T. H. Jr.; Harrison, R. J. J. Chem. Phys. 1992, 96, 6796−6806. (34) Fiacco, D. L.; Hunt, S. W.; Leopold, K. R. J. Am. Chem. Soc. 2002, 124, 4504−4511. (35) Kuczkowski, R. L.; Suenram, R. D.; Lovas, F. J. J. Am. Chem. Soc. 1981, 103, 2561−2566. (36) Kim, J.; Lee, J. Y.; Lee, S.; Mhin, B. J.; Kim, K. S. J. Chem. Phys. 1995, 102, 310−317. (37) Benedict, W. S.; Gailar, N.; Plyler, E. K. J. Chem. Phys. 1956, 24, 1139−1165. (38) Beichert, P.; Schrems, O. J. Phys. Chem. A 1998, 102, 10540− 10544. (39) Chackalackal, S. M.; Stafford, F. E. J. Am. Chem. Soc. 1966, 88, 723−728. (40) Givan, A.; Larsen, L. A.; Loewenschuss, A.; Nielsen, C. J. J. Chem. Soc., Faraday Trans. 1998, 94, 827−835. (41) Bykov, A. D.; Makushkin, Y. S.; Ulenikov, O. J. Mol. Spectrosc. 1983, 99, 221−227. (42) Engdahl, A.; Nelander, B. J. Chem. Phys. 1989, 91, 6604−6612. (43) Rozenberg, M.; Loewenschuss, A. J. Phys. Chem. A 2009, 113, 4963−4971. (44) Kauppi, E.; Halonen, L. J. Phys. Chem. 1990, 94, 5779−5785. (45) Salmi, T.; Hänninen, V.; Garden, A. L.; Kjaergaard, H. G.; Tennyson, J.; Halonen, L. J. Phys. Chem. A 2008, 112, 6305−6312. (46) Salmi, T.; Hänninen, V.; Garden, A. L.; Kjaergaard, H. G.; Tennyson, J.; Halonen, L. J. Phys. Chem. A 2012, 116, 796−797. (47) Halonen, L.; Carrington, T. Jr. J. Phys. Chem. 1988, 88, 4171− 4185. (48) Kauppi, E.-J. Studies of Overtone Spectra and Potential Energy Surfaces of Some Polyatomic Molecules. Ph.D. Thesis, University of Helsinki, Finland, 1992. (49) Morse, P. M. Phys. Rev. 1929, 34, 57−64. (50) Honkonen, J. Fysiikan Matemaattiset Menetelmät I, 2nd ed.; Limes ry: Finland, 2005. (51) Pipek, J.; Mezey, P. G. J. Chem. Phys. 1989, 90, 4916−4926. (52) Hill, J. G.; Platts, J. A.; Werner, H.-J. Phys. Chem. Chem. Phys. 2006, 8, 4072−4078. (53) Boughton, J. W.; Pulay, P. J. Comput. Chem. 1993, 14, 736−740. (54) Werner, H.-J. et al. MOLPRO, version 2010.1, a package of ab initio programs; http://www.molpro.net (2010). (55) Saebø, S.; Tong, W.; Pulay, P. J. Chem. Phys. 1993, 98, 2170− 2175. (56) Weigend, F.; Köhn, A.; Hättig, C. J. Chem. Phys. 2002, 116, 3175−3183. (57) Weigend, F. Phys. Chem. Chem. Phys. 2002, 4, 4285−4291. (58) Yousaf, K. E.; Peterson, K. A. J. Chem. Phys. 2008, 129, 184108. (59) Lane, J. R.; Kjaergaard, H. G. J. Chem. Phys. 2010, 132, 174304. (60) Lane, J. R.; Kjaergaard, H. G. J. Chem. Phys. 2009, 131, 034307. (61) Demaison, J.; Herman, M.; Lievin, J.; Rudolph, H. D. J. Phys. Chem. A 2007, 111, 2602−2609. (62) Rozenberg, M.; Loewenschuss, A.; Marcus, Y. Phys. Chem. Chem. Phys. 2000, 2, 2699−2702. (63) Atkins, P.; De Paula, J. Atkins Physical Chemistry, 8th ed.; Oxford University Press: New York, 2006. (64) Jensen, F. Introduction to Computational Chemistry, 2nd ed.; John Wiley and Sons, Ltd.: New York, 2007. (65) Hintze, P. E.; Kjaergaard, H. G.; Vaida, V.; Burkholder, J. B. J. Phys. Chem. A 2003, 107, 1112−1118. (66) Barletta, P.; Shirin, S. V.; Zobov, N. F.; Polyansky, O. L.; Tennyson, J.; Valeev, E. F.; Császár, A. G. J. Chem. Phys. 2006, 125, 204307. (67) Kjaergaard, H. G.; Garden, A. L.; Chaban, G. M.; Gerber, R. B.; Matthews, D. A.; Stanton, J. F. J. Phys. Chem. A 2008, 112, 4324−4335. (68) Iogansen, A. V. Spectrochim. Acta 1999, 55, 1585−1612. (69) Mathematica, version 8.0; Wolfram Research Inc.: Champaign, IL, 2011.

(70) Tennyson, J.; Zobov, N. F.; Williamson, R.; Polyansky, O. L.; Bernath, P. J. Phys. Chem. Ref. Data 2001, 30, 735−831.

2879


Ab Initio Structural and Vibrational Investigation of Sulfuric Acid

Recommend Documents