Ab Initio Calculations of Thermochemical ... - ACS Publications

Article pubs.acs.org/JPCA

Ab Initio Calculations of Thermochemical Properties of Methanol Clusters Muhammad Umer and Kai Leonhard* Lehrstuhl für Technische Thermodynamik, RWTH Aachen University, Schinkelstrasse 8, 52062 Aachen, Germany S Supporting Information *

ABSTRACT: Highly accurate ab initio calculations of binding enthalpies and entropies of gas phase clusters of methanol have been performed, yielding uncertainties smaller than 1 kJ/ mol per hydrogen bond in the Gibbs free energy of reaction. This requires quantum chemical RIMP2 and CCSD(T) postHartree−Fock methods with basis sets up to aug-cc-pV5Z for energy calculations. An analysis of topological symmetry and hindered rotor effects proves necessary for reliable entropies. This approach goes beyond the rigid rotor plus harmonic oscillator method implemented in standard quantum mechanics software tools. The results demonstrate that (1) thermochemical methanol cluster properties can nowadays be obtained by ab initio methods with an accuracy comparable to or even better than that of the experimental data available, especially for larger species that cannot be studied directly by experiments and (2) cooperativity effects and state-dependent cluster distributions cause a strongly varying average enthalpy and entropy per bond as a function of temperature and density for methanol.

■

INTRODUCTION The formation of hydrogen-bonded alcohol clusters is crucial for solvation effects, potential biochemical and fuel-cell applications, phase equilibria, nucleation and condensation phenomena, etc. Therefore, hydrogen bonds have been studied using several theoretical and experimental approaches. For an overview see, for example, Laksmono et al.1 and references therein. Although the main interest has been in the structure of these clusters and their hydrogen bonds, much less work has been done on the thermodynamic properties of methanol clusters, despite the fact that these properties are relevant for all above-mentioned applications and phenomena. Reasons may be the difficulty to obtain accurate experimental data and a lack of theoretical methods for accurate entropy calculations that were available until recently. We are presently interested in accurate enthalpies and entropies for association reactions to develop a predictive equation of state for the computation of phase equilibria that requires such data. Therefore, we briefly review the available experimental thermochemical data, the relevant theoretical studies, and new methods for accurate entropy calculation in this introduction before we evaluate these methods in the main part of the present paper. Curtiss and Blander2 reviewed the experimentally determined thermodynamic data on methanol association available until 1988. These data include enthalpy, entropy, and equilibrium constants for the dimerization from seven sources determined by different methods (pVT data, infrared and NMR spectroscopy, heat capacity, and thermal conductivity measurements) that show quite some scatter. Often the data do not agree within the specified experimental errors or no error estimates are given. The review contains also results for the trimerization, the tetramerization and the octamerization, most © XXXX American Chemical Society

of which are determined from pVT data. In this case, one should keep in mind that for an analysis of these data, certain species have to be postulated and that missing species lead to systematic errors in the resulting equilibrium constants. Additionally, enthalpy and entropy are derived from these constants and therefore have larger uncertainties. There is a more recent dimerization constant,3 but we have not found any more recent direct experimental thermochemical results regarding larger methanol clusters. Several papers appeared, however, on spectroscopic studies since then, e.g., refs 12−18 in Laksmono et al.1 The difficulties to obtain accurate experimental thermochemical data for methanol cluster formation raises the question whether theoretical methods may be better suited to tackle this problem. Though it may be hard to judge their accuracy by the existing experimental thermochemical data, the situation is different for the density of methanol vapor, because an accurate IUPAC reference equation exisits.4 As indicated above, such density data do not allow for an unambigous determination of thermochemical properties of all existing clusters. The opposite approach, however, is possible. If an ab initio cluster model is set up, the agreement of the predicted densities with experimental ones is a very good indicator for the quality of the model. Many researchers have been interested in understanding the electronic structure of methanol clusters, in analyses of the hydrogen bonds, and in shifts in the frequencies of the OH and CO vibration. Huang and Day,5 e.g., studied the structure and energy of methanol clusters, performed NBO (natural bond Received: September 7, 2012 Revised: January 16, 2013

A

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

The Journal of Physical Chemistry A

Article

orbital) analyses and QTAIM (quantum theory of atoms in molecules) to characterize the H bonds and discussed OH frequency shifts at the harmonic level. Similarly, Gonzalez et al.6 carried out calculations of geometries, frequencies and energies for ethanol dimers and trimers with a medium sized basis set. They found the all-trans conformation to be the most stable one in the dimer and trimer and compared calculated and experimentally determined OH frequencies of these species. Although most studies discuss the OH and CO frequencies at the harmonic level, Buck et al.7 developed a semiempirical potential model to investigate structure and energy of methanol clusters. They studied the OH and CO vibrational modes including third- and fourth-order anharmonicities. The authors stress that “The necessity of an anharmonic treatment cannot be overemphasized”. This is even more true for the large amplitude motions in alcohol clusters, as we will show in this work. Muñoz-Caro and Niño8 considered the anharmonicities in the water dimer by calculating variational solutions of the Schrödinger equation for several motions of the atoms in a water dimer. All studied motions are 1-dimensional. The focus was on anharmonic vibrations and their comparison to spectroscopic results but not on low-frequency modes particularly relevant for thermochemical properties. Nevertheless, they calculated the equilibrium constant K for dimerization and found it to be about a factor 3 too small compared to experimental results, indicating either a too weak energetic attraction or too small entropy contributions from the large amplitude motions in the dimer. However, in most studies on alcohols no attempt was made to compute thermochemical properties like enthalpy and entropy even though they determine the dominating cluster species. An exception is the work by Sum and Sandler who calculated Hatree−Fock (HF) geometries with MP2/aug-ccpVDZ single point energies (SPE) for methanol, ethanol, 1propanol, and methanethiol clusters.9 Because of computational limitations, they assumed the entropy change for the addition of a molecule to a cluster be cluster-size independent and thus obtained equilibrium constant ratios for clusters of different sizes. However, they identified this assumption as a severe simplification. Moreover, absolute equilibrium constants are not obtainable by this method. The focus of the present work is mainly on thermochemical properties of methanol clusters. Such properties are interesting not only to describe the cluster distribution but also to determine association parameters for the PCP-SAFT (perturbed-chain polar statistical associating fluid theory) equation of state (EOS).10 This could lead to a fully predictive EOS because we recently developed an approach to obtain size, shape, and dispersion energy parameters of PCP-SAFT, i.e., all parameters not related to association, from quantum mechanically obtained molecular properties.11 We identified an approach by Wolbach and Sandler12 as a promising starting point to obtain the additional association parameters required for hydrogen-bonding compounds. In its original form, the dimerization approach shows four main weaknesses that are mainly due to the limitation of available methods and computational capacity when the approach was developed by Wolbach and Sandler in 1997.12 First, the quantum mechanical calculations were conducted applying the HF method and density functional theory (DFT) that used Becke’s threeparameter functional combined with the gradient corrected correlation functional of Lee, Yang, and Parr (B3LYP).13−16

This method has recently been shown to perform worse than average for hydrogen-bonded complexes compared to other DFT and post-HF methods available today.17−19 Second, a rigid rotor, harmonic oscillator (RRHO) approximation was applied to model the effect of the motion of the nuclei on thermodynamic functions, in which rotations of parts of a molecule or cluster, also called hindered rotations, were neglected. Third, only one conformation of each species (monomers and dimers only) was considered and, fourth, cooperativity was thus neglected. Especially the neglect of cooperativity effects shows a major impact on the calculation of the enthalpy change ΔH resulting in a difference of more than 10 kJ/mol per hydrogen bond for methanol clusters.20 In comparison, a sensitivity analysis revealed that reaching the same 10 K standard error in the normal boiling temperature computed by an EOS for hydrogen-bonding compounds as was found for non hydrogen-bonding ones in ref.,11 requires an accuracy of 1 kJ/mol in ΔG for short-chain 1-alkanols. The discussion of the methods employed by Wolbach and Sandler showed that these methods were not sufficient to compute thermochemical data of dimerization with the accuracy required to judge the validity of their basic assumption, namely that the transfer of thermochemical gas phase hydrogen bond data into SAFT parameters is possible. Although much progress has been made regarding the accuracy of the electronic energies,9,20,21 cooperativity,6,7,9,20−22 and the consideration of multiple conformations,21 the studies of thermochemical properties of alcohol clusters we are aware of8,9,12 do not provide results that are accurate enough for our purpose. Therefore, we investigate in this paper whether the required accuracy can be achieved with modern quantum mechanical and statistical thermodynamic methods. Methanol clusters are a good starting point due to their computational simplicity, because of the availability of experimental data for this molecule and also to develop the methodology to study larger alcohols. A comparison to experimental gas densities shows that the required accuracy can indeed be reached.

■

METHODS Quantum Mechanics. To compute thermochemical properties of methanol monomers and clusters, the following procedure was applied. All geometry optimizations and frequency calculations were performed with the B3LYP functional13−16 and TZVP basis set23,24 with “tight” convergence criteria and an ultrafine grid in Gaussian 0925 because this method is known to give accurate geometries and frequencies at moderate computational costs.26 Starting from the initially found (local) minimum of the potential energy surface (PES), rotations of dihedral angles were performed to find all possible conformers including the global minimum. Because the average error of B3LYP/TZVP for intermolecular interactions is about 3 kJ/mol,17,27 it does not achieve the required accuracy of 1 kJ/mol. Therefore, we employed RIMP228,29 SPE calculations using the counterpoise correction (CPC) for the basis set superposition error.30 Encpc = En(R n) −

∑ E MB

Mi

i

−

∑ E MB (R n) n i

i

B

i

(R M i) +

∑ E MB

Mi i

(R n)

i

(1)



Article

Here Ecpc n is the counterpoise corrected hydrogen bond energy of an n-mer at a given basis set and En is the uncorrected energy of that cluster where Rn indicates the geometry of the n-mer and RMi indicates the lowest-energy geometry of monomer i. EBMni is a monomer energy computed with the n-mer basis set and, finally, EBMMii is a monomer energy in the appropriate monomer basis. In eq 1, the first two terms represent the uncorrected interaction energy and the last two terms represent the counterpoise correction, ΔEcpc n . Initial calculations used the aug-cc-pVTZ basis set31 (X1 = 3). For the most-relevant species, additional calculations using the aug-cc-pVQZ basis sets31 (X2 = 4) have been performed and used to extrapolate the results to the complete basis set limit 32 (Ecbs n ). Encbs ≈ E aug‐cc‐pVX2Z + − E aug‐cc‐pVX1Z)

X13 X 2 3 − X13

bonded clusters and to sort out irrelevant clusters. For species for which we found several conformations, we took the effect of their Boltzmann distribution on enthalpy and entropy into account when using the RRHO model. Large Amplitude Motions. Even though the RRHO model is used in most calculations as a standard method, it is nowadays consensus in the literature that the anharmonicity of large amplitude motions has to be considered for accurate thermochemical and kinetic data.36−38 Large amplitude motions occur typically in molecular systems with weak bonds, i.e., in transition states of reactions and in intermolecular (hydrogen bonded) clusters. Standard software, like Gaussian 09, can handle anharmonicity by Taylor expanding the PES up to fourth order around the energy minimum. This leads to improved vibrational frequencies for relatively rigid molecules,39 but turned out to be inadequate for the large amplitude motions of intermolecular complexes in some preliminary calculations in the present work. The state-of-the-art way is to replace the harmonic oscillator contributions of such motions by those of uncoupled 1D hindered rotations (HR).36−38 The reason for the uncoupled 1D treatment is 2-fold: (1) a fulldimensional treatment is not practical for all but the smallest molecules40,41 and (2) the 1D approach has proven accurate for some special cases (n-alkanes, alcohols, thiols, ethers, and thioethers) by comparison to the coupled 2D case.42,43 The largest rigorously tested system is the hydrogen peroxide monomer,41 which is a demanding test case and probably gives an upper bound for the error of decoupling one HR from all other anharmonic modes. The error in the partition function when using an exact 1D treatment decoupled from all other motions compared to the exact full dimensional case is 24% (averaged from 200 to 5000 K), corresponding to about 0.5 kJ/ mol at ambient temperature. To study anharmonic degrees of freedom, we employ the TAMkin package.38 TAMkin is an open source python package that allows the full QM treatment of 1D HR. For this purpose, it contains functions to read in the results of Gaussian 09 calculations of rotational scans, geometry and Hesse matrix information. It fits an analytical function to the scanned data, sets up the Schrödinger equation for hindered rotation, expands the wave functions in sine and cosine functions, and solves the equation by a variational approach.42 From the resulting energy spectrum, TAMkin computes the partition function and the thermodynamic properties. Consideration of the lowestfrequency mode, i.e., the rotation around the hydrogen bond only, did not achieve the desired accuracy, as a comparison with experimental virial coefficient data showed. Therefore, we treat all cluster modes that originate from the lost translational and rotational monomer degrees of freedom (3 + 3) and the methyl rotations (2) anharmonically. This gives eight modes in the case of a methanol dimer. To our knowledge, such a sophisticated treatment of anharmonic degrees of freedom has not yet been described in the literature before. Because TAMkin was not able to handle all these modes correctly, we implemented the following extensions. In TAMkin, bond bends are always treated harmonically. However, in the normal-mode analysis of a methanol dimer, e.g., mode 4 is a large-amplitude bond bend (see Results and Discussion section). If this mode is scanned with all other degrees of freedom relaxed, the bond distance varies such that a full rotation of one molecule around the other is possible and the resulting potential has the same symmetry as that of a dihedral scan. Therefore, we extended the TAMkin input

(E aug‐cc‐pVX2Z (2)

Such extrapolation results are labeled “aug-cc-pV(X1X2)Z”. From the estimated energy, Ecbs n , of a cluster consisting of n molecules having HBNn intra cluster hydrogen bonds, the hydrogen bond energy per bond, ΔEHB is ΔEHB =

Encbs − nE1cbs HB

Nn

(3)

Similar equations hold for the thermodynamic contributions to the enthalpy, ΔH, and entropy, ΔS. For the systems studied, RIMP2 is about 30 and 100 times faster than conventional MP2 with the aug-cc-pVTZ and augcc-pVQZ basis sets, respectively. The deviations between RIMP2 and MP2 for small clusters were on average below 0.2 kJ/mol with a maximum deviation of 0.5 kJ/mol and are therefore not critical compared to the required accuracy. As a rule of thumb, using a basis set for RIMP2 with a cardinal number one larger than that of a chosen basis set for an MP2 calculation, the RIMP2 calculation is both faster and more accurate. For the methanol dimer, we additionally investigated the effect of truncating electron correlation and basis set by a CCSD(T)/aug-cc-pVTZ and an RIMP2/aug-cc-pV5Z calculation that yielded negligible corrections of 0.03 and 0.1 kJ/mol, respectively. Due to their high computational demand, a factor 100 compared to that of RIMP2/aug-cc-pVQZ, they were applied to the dimer only. Thermodynamics. Quantum chemistry yields the electronic energy states of molecular systems at zero Kelvin but does not include the concept of temperature. Therefore, we have to define models for the kinetic energy of the nuclei in molecular translation, rotation, and vibration as well as for their potential energy in vibrations and use these models to obtain thermodynamic functions. A standard model for this purpose is the RRHO model, which is described in textbooks,33,34 implemented in common software,25,35 and therefore most widely used for this purpose. In this model, the PES is Taylor expanded to second order at its minimum. Hence, all intramolecular motions become independent harmonic oscillations. It is simple to use, is computationally cheap, and has proven to yield reasonable results for typical, stable molecules that are relatively rigid.35 Therefore, it was used by Wolbach and Sandler12 and most other researchers. We apply it here as low-level model to compute the thermodynamics of hydrogen C



Article

thermochemical data obtained from theoretical calculations by experimentally determined pVT data at relatively low densities. (2) For future application in PC-SAFT, average or effective reaction enthalpies and entropies are required because PC-SAFT association parameters are always cluster size independent. The cluster model allows for the calculation of such average quantities. Because the literature20,48 suggests that linear (l), cyclic (c), and branched (b) clusters of up to six molecules may be relevant for methanol thermochemistry, reactions to form all these species were included into the model:

routine to read also bond angle scans and store the scan type information as additional information. Originally the moment of inertia is computed around the axis denoted by vector a connecting the second and the third atom of the dihedral scan with positions r2 and r3, respectively (a = r3 − r2). For bond bends, we compute it around an axis going through the second atom and being perpendicular to the plane containing the three atoms that define the scan (a = (r1 − r2) × (r3 − r2)). The calculation of the moment of inertia itself and the solution of the Schrödinger equation is then done in exactly the same way as before. To model the stretching motion of the hydrogen bond, we scanned the O−O distance, relaxed all other modes and fitted a Morse potential to the so-obtained data points. The focus was on a good fit near the energy minimum. Larger deviations at larger energies and separations do not affect the final results much because the poorly described high-energy states are almost not populated at the relevant temperatures. We then analytically calculated the possible energy states,44 summed them up to get the partition function and replaced the respective harmonic oscillator contributions. For molecular and intramolecular rotations, a careful consideration of symmetry is essential. Consider, for example, the hindered rotation of the OH group with respect to the CH3 group in a methanol monomer. Here, one finds a symmetry number of three because a rotation by 120° creates an identical conformation because all three H atoms are equivalent. The partition function obtained without taking symmetry into account then has to be divided by the symmetry number to correct for the overcounting of identical quantum states. Although these facts are textbook knowledge in principle,34 the determination of correct symmetry numbers is by far not trivial in the general case, especially for flexible molecules and complexes, where so-called topological symmetry in contrast to geometrical symmetry is needed. This is described very well in two recent articles.45,46 We will describe in the Results and Discussion section how we applied the procedures described therein to determine the symmetry of the motions of the methanol dimer. Because the topological symmetry determined this way does not always correspond to the symmetry of the potential employed in TAMkin, we did not use any symmetry there but divided the resulting partition function by the topological symmetry number afterward. We note that the use of (all) HRs automatically includes all relevant conformers of the clusters in the partition function. Therefore, an explicit averaging over conformers is not necessary in this approach and all results are given with respect to the lowest-energy conformer. Gas Phase Cluster Model. Because experimental thermochemical data are most often indirectly derived from experimental data, their accuracy is sometimes questionable and results from different sources often do not agree within the reported error bars.2 For a direct comparison of our results to experimental results, we employ a model of reactive clusters. Such cluster models have been combined with EOS models already in the 1970s, see, for example, Heidemann and Prausnitz.47 Due to the large number of required thermochemical data inaccessible at that time, crude simplifications were required to obtain simple analytical models with few adjustable parameters. Hence, the models were later replaced by physical association models like SAFT. In this work, we use the clustermodel for two reasons. (1) In combination with the ideal gas EOS it provides an unambiguous means to validate

nCH3OH ⇌ (CH3OH)(nα)

n = 2, ..., 6 and α = l, c, b (4)

Note that there are no branched and cyclic dimers. For each valid reaction, the law of mass action is given by eq 5. K n(α)

1−n ⎛ ΔH (α) − T ΔS(α) ⎞ xn(α) ⎛ p ⎞ n n ⎟ ⎜ = exp⎜ ⎟ = xn ⎜ ⊖⎟ RT ⎝ ⎠ 1 ⎝p ⎠

(5)

Here ΔH and ΔS are the enthalpy and entropy of the reaction, respectively, K is the equilibrium constant, x is the mole fraction, p is the pressure and p⊖ = 1 atm is the reference pressure used in the definition of the entropies, R is the universal gas constant, and T is the temperature. The subscript n denotes the number of monomers in a cluster and superscript (α) its type. Using the normalization condition for the mole fractions,

∑ ∑ xn(α) = 1 n

α

(6)

we obtain a nonlinear system of as many equations as mole fractions to be determined. The desired cluster mole fractions can be obtained by solving the system of equations iteratively. By coupling the cluster model with the ideal gas law, we assume that the change of particles due to association has an effect on the pressure, but the effect of molecular volume, dispersion, and electrostatics is negligible. By calculating the respective contributions from PC-SAFT, we have shown that these effects are at least two orders of magnitude smaller than the effect of association at all conditions where the cluster model was applied. The expression for the density ρ for a given pressure by the cluster + ideal gas model reads49 ρ=

■

p RT

6

∑ ∑ nxn(α) n=1 α

(7)

RESULTS AND DISCUSSION RRHO Results. The RRHO+RIMP2/aug-cc-pVTZ// B3LYP/TZVP method was used as low level model for initial computations of methanol clusters to determine the relevant species for further ab initio calculations. For this purpose, linear (l), cyclic (c), and branched (b) clusters were considered for methanol consisting of up to six molecules, as suggested in the literature.48 All species investigated at this level are shown in Scheme 1. Two stable cyclic tetramers were found (with methyl groups up−down−up−down (udud), and uudd with respect to the plane containing all oxygen atoms) as well as two hexamers (ududud and uududd). For some species additional corrections were added: for the dimer and all other linear clusters, R ln 2 was added to the entropy because two dimers that are mirror images to each other exist. For the longer linear clusters, this D



Article

Scheme 1. All Methanol Clusters Studied in This Worka

The entropy change to form a cyclic cluster is generally more negative than that for a linear one, but the entropy change per bond is less negative, again due to the higher number of hydrogen bonds. The decrease in enthalpy and Gibbs free energy change with increasing cluster size hence shows the importance of cooperativity for this system. The branched clusters are generally less stable than linear and cyclic structures (see Ohno et al.48 for the energies of larger branched clusters). This is due to unfavorable electrostatic interactions of the two donor OH groups pointing toward the same acceptor. The just described thermochemical data have been used in the cluster model described in section Gas Phase Cluster Model to assess the accuracy of the data and to decide which species are the most-relevant ones in case more time-consuming higher-level model calculations are needed. The results show that dimer, trimers, c-tetramer, and c-pentamer are the most relevant clusters of methanol at densities where the model is applicable; cf. Table 2. The dimer is dominant at saturation pressure but c-tetramer and c-pentamer show potentially relevant mole fractions at increased pressures. A comparison of the density obtained from the cluster model eq 7 with experimental densities4 showed that the increase in densities over those of the ideal gas is roughly a factor 10 too small compared to the experimental data. Therefore, additional higher level calculations have been performed for the most important species at the RRHO+RIMP2/aug-cc-pVQZ// B3LYP/TZVP level and the RIMP2 energies have been extrapolated to the basis set limit; see Table 3. At this level, hydrogen bonds are about 1−2 kJ/mol more attractive than at the lower level but the resulting densities still do not agree well with experimental ones (now roughly a factor 5 too small). As described in the Methods, we checked with an even larger basis set and coupled-cluster calculations that the quantum chemical calculations are very accurate and that the remaining error cannot be an explanation for the observed differences to the experimental data. Because an HR treatment is required for accurate entropies as discussed in the Methods, such corrections have been performed for the most relevant clusters to check for possible errors of the RRHO approximation in this specific application. HR Results for Monomers. Because a considerable improvement of the accuracy of the standard entropy of ethanol was reported in the literature when changing from the RRHO to HR model,43 we used this model for methanol

a

The B3LYP/TZVP geometries of the clusters are provided in the Supporting Information.

may seem incorrect because each hydrogen acceptor oxygen is chiral. The diastereomers thus created, however, converge to other species when being optimized and are therefore not counted again. The linear trimer, for example, changes to the cyclic trimer when the mirror image of only one acceptor oxygen is created. For the cyclic structures consisting of n monomers we applied a rotational symmetry number of n due to topological symmetry46 because vibrations can couple to rotations and create n identical arrangements of equivalent nuclei during one rotation. This is called “topological symmetry” because this is not possible by simple geometric symmetry operations. For the same reason, a topological symmetry number of two was applied to the branched trimer. The results in Table 1 show that cyclic structures are more stable than the linear ones of the same size due to a higher number of hydrogen bonds and, except for the trimer, which has a high ring tension, also due to stronger hydrogen bonds.

Table 1. Comparison of Methanol Clusters at the RRHO+RIMP2/aug-cc-pVTZ//B3LYP/TZVP level at 298 Ka methanol cluster

rotational symmetry no.

(total) topological symmetry no.

ΔEHB, kJ/mol

ΔHHB, kJ/mol

ΔSHB, J/(mol K)

ΔGHB, kJ/mol

dimer l-trimer c-trimer b-trimer l-tetramer c-tetramer (udud) c-tetramer (uudd) l-pentamer c-pentamer l-hexamer c-hexamer (ududud) c-hexamer (uududd)

1 1 1 1 1 2 1 1 1 1 3 1

1 1 3 2 1 4 4 1 5 1 6 6

−22.9 −28.7 −23.1 −21.0 −28.6 −30.6 −29.5 −30.7 −32.3 −32.1 −33.4 −32.9

−16.1 −21.8 −18.2 −14.2 −21.8 −25.2 −24.2 −23.9 −26.7 −25.3 −27.7 −27.2

−93.6 −114.4 −87.2 −102.8 −110.1 −101.1 −98.6 −114.1 −104.8 −121.0 −108.8 −111.0

11.8 12.3 7.8 16.4 11.1 5.0 5.2 10.1 4.5 10.8 4.7 5.9

a

All thermodynamic data are reported per hydrogen bond. Table S1 in the Supporting Information provides geometry information and additional results for all species. E



Article

Table 2. Methanol Cluster Distribution at RRHO+RIMP2/aug-cc-pVTZ//B3LYP/TZVP at 298 Ka pressure, bar

x2

x(3,l)

x(3,c)

x(4,c)

x(5,c)

0.169 1.00

1.4 × 10−3 8.3 × 10−3

1.4 × 10−6 4.7 × 10−5

2.2 × 10−6 7.5 × 10−5

2.8 × 10−6 5.6 × 10−4

9.0 × 10−8 1.1 × 10−4

a

The remaining species are mainly monomers. Mole fractions of all other species are at least 2 orders of magnitude smaller than the smallest displayed one. 0.169 bar is the saturation pressure.

1.51 J/(mol K) and enthalpy change of −0.095 kJ/mol at 298 K. The results presented in Table 4 show the improvement achieved by including the HR treatment.

Table 3. Relevant Methanol Clusters at the RRHO+RIMP2/ aug-cc-pV(TQ)Z//B3LYP/TZVP Level (Extrapolations) at 298 Ka methanol cluster

ΔEHB, kJ/mol

ΔHHB, kJ/mol

ΔSHB, J/(mol K)

ΔGHB, kJ/mol

dimer l-trimer c-trimer c-tetramer c-pentamer

−24.1 −29.8 −24.1 −31.5 −34.1

−17.2 −22.9 −19.1 −26.1 −28.4

−93.6 −114.4 −87.2 −98.8 −104.8

10.7 11.2 6.9 3.4 2.9

Table 4. Comparison of Methanol Monomer Entropy from RRHO and RRHO + HR Models and from Experimental Results50

a

Table S2 in the Supporting Information provides additional results for all species.

methanol

S(RRHO), J/(mol K)

S(RRHO +HR), J/(mol K)

S(Experiment), J/(mol K)

ΔS(internal rotation), J/(mol K)

237.84

239.35

239.87

1.51

HR Results for the Dimer. A normal-mode analysis with Gaussian 09 for the methanol dimer yielded the eight lowestfrequency modes that we decided to treat anharmonically. Two of the eight modes originate from the CH3 rotations already treated as HR in the monomer and six modes originate from three translational and rotational degrees of motion of one of the formerly unbonded molecules. The scans to obtain the 1D potentials of these modes were sometimes not trivial and are therefore described here. To follow a mode, we visually inspected a normal mode vibration in gaussview and tried to identify a suitable internal rotation to replace the normal mode, as suggested in the literature.36,38,51 When the identified dihedral angle is then scanned and all other internal coordinates are relaxed, the scan often does not follow the desired rotation but changes to some other movement due to the movement coupling with others at large displacements from the minimum energy structure. Therefore, we often had to freeze one or two internal coordinates to make sure that only the scanned internal coordinate changed substantially. As suggested in the literature,38 we subtracted the HO contribution corresponding

monomers to check our procedures because accurate experimental entropies are available for these species. In the methanol monomer, the OH group and the CH3 group can rotate with respect to each other, as indicated by the arrows in Figure 1. The scan of the rotation, which has a symmetry number of 3, is straightforward, and the results are also shown in Figure 1. All figures of large-amplitude motions indicate the motion on the left-hand side and show the results on the right-hand side. The arrows indicate the corresponding normal mode, the axis of rotation is always (almost) perpendicular to the paper plane, except for the tetramer. In the monomer case, the dihedral angle defined by atoms 1, 2, 3, and 4 is scanned. Symbols represent individual points of the scan (all performed at the B3LYP/TZVP level), the solid line is the fitted potential function, the thin horizontal lines represent the energy states of the motion, and the length of the thicker horizontal lines indicates the relative population of the respective state at 300 K. The replacement of the corresponding HO by a HR results in an entropy change of

Figure 1. Mode 1 of the methanol monomer: OH rotation. F



Article

are longer than some on the donor. The potential of this motion was determined by scanning D(7,3,2,1). When the donor rotates counterclockwise, i.e., opposite to the movement shown in Figure 3, the hydrogen bond gets completely broken and unfavorable dipole−dipole interactions push the molecules apart when the original acceptor’s orientation is fixed. When we allowed relaxation of the acceptor orientation, a favorable dipole−dipole interaction occurred, but atoms 7, 3, and 2 arranged linearly and a further scan of this mode became impossible. A clockwise rotation of the donor causes the OH group of the acceptor to rotate toward the donor oxygen atom. After passing a 5 kJ/mol energy barrier, this leads to a second energy minimum, where the former donor is now the acceptor and vice versa. Closer inspection shows that all angles and bonds are the same, except for the sign of some dihedral angles. This means that each species with D(7,3,2,1) is a mirror image of another one with D′(7,3,2,1) = 305° − D(7,3,2,1). Therefore, in total only one-half of the rotation was scanned. TAMkin was extended to apply a translation of −55° to D(7,3,2,1) so that an even symmetry could be imposed on the potential and the wave functions, making both minima exactly identical regarding energy and curvature. The plot in Figure 3 shows the correct angles after a back-translation was performed to plot this figure. Because all species along the rotation are distinctly different, the contribution of this mode is considerable due to the second minimum; cf. Table 5. Mode 3, Acceptor Rotation. For the scan of D(2,9,7,8) two other dihedrals, D(3,2,7,9) and D(3,2,1,7), have to be frozen to avoid coupling with other modes, especially mode 1. The scan has to be performed from the stable minimum in both directions because it terminates at zero or 360° due to a collinearity of the atoms 2, 7, and 8. The possibly poor representation of the potentials at this angle has, however, no relevant consequences due to its high energy in this region (>20 kJ/mol, cf. Figure 4). The effect of the anharmonicity of this mode on enthalpy, and entropy is comparably small. Mode 4, COO Bond Bend. The scan of the angle defined by atoms 3, 2, and 7 (A(3,2,7)) has to be performed in both directions from the equilibrium angle of 111° again. Decreasing the angle leads to a movement opposed to the arrows in Figure 5 and to a breakage of the hydrogen bond corresponding to the barrier around 0°. The barrier is “only” about 22 kJ/mol

to the scanned internal rotation, not the one corresponding to the most similar normal mode. The resulting thermodynamic properties of all modes can be found in Table 5. Table 5. Properties of Large Amplitude Motions in the Methanol Dimera mode

ΔH, kJ/mol

ΔS, J/(mol K)

1: OHO rotation 2: donor rotation 3: acceptor rotation 4: COO bond bend 5: donor methyl rotation 6: hydrogen bond stretch 7: acceptor H rotation 8: donor H rotation topological symmetry

0.04 0.12 −0.21 −0.21 −0.11 0.16 0.90 0.43 0

4.69 8.11 0.49 0.84 0.99 0.45 8.00 9.88 −R ln 2 = −5.76

ΔH and ΔS are the differences between the HR and the HO models in enthalpy and entropy at 298 K, respectively. a

Mode 1 is the rotation around the OHO bond; cf. Figure 2. It is described by the dihedral angle defined through atoms 9, 7, 2, and 3 (this will be abbreviated by D(9,7,2,3)). To prevent an acceptor rotation such that the hydrogen atom 8 flips to the left when the methanol monomer in the back rotates clockwise, D(2,9,7,8) was frozen. This latter movement corresponds to mode 3, acceptor rotation, and will be described below. The local minimum at about 40° and 2.5 kJ/mol appears only in this projection, it disappears when all degrees of freedom are relaxed (this is true for any minimum discussed in this section, HR Results for the Dimer, except for the global one). We can see that the potential widens considerably at about 1.6 and 2.5 kJ/mol compared to the harmonic one. Therefore, the states get considerably closer in this region. At even higher energies corresponding roughly to a free rotor, the spacing increases again. Altogether, the effects lead to a substantial increase in entropy and a minimal increase in enthalpy compared to the harmonic case at intermediate temperatures (200−400 K); see also Table 5. Mode 2, Donor Rotation. Figure 3 shows the rotation of the hydrogen bond donor as a whole around its OC bond in a methanol dimer. The acceptor just performs some compensating movement, even though the arrows on some of its atoms

Figure 2. Mode 1 of the methanol dimer: rotation around OHO bond. G



Article

Figure 3. Mode 2 of the methanol dimer: The donor rotation.

Figure 4. Mode 3 of the methanol dimer: The acceptor rotation.

Figure 5. Mode 4 of the methanol dimer: COO bond bend.

Figure 6. Mode 5 of the methanol dimer: Donor methyl rotation.

H



Article

Figure 7. Mode 6 of the methanol dimer: Hydrogen bond stretch.

Figure 8. Mode 7 of the methanol dimer: Acceptor H rotation.

is sufficient for the partition function around 300 K. At larger distances, long-range Coulomb interactions are relevant but neglected in the Morse potential. The effect of the anharmonicity of this mode is small; cf. Table 5. Mode 7, Acceptor H Rotation. We scanned D(2,9,7,8) (in both directions) and froze D(2,7,9,12) to obtain the potential for this mode (Figure 8). The scan along arrows (decreasing dihedral angle) leads to bond breakage and roughly antiparallel dipoles next to each other, hence the relatively low barrier (because a coupling with the rotation of H1 is inevitable). The minimum near 0° corresponds to an exchange of donor and acceptor that is topologically identical to the initial structure. A further rotation from this position corresponds to the donor H rotation. For an increasing angle H8 rotates alone. At the first maximum, at 180°, H1 is in between the 2 lone pairs of O7. To go beyond 200°, D(2,3,7,9) has to be frozen again (up to 280°, then 3 of these atoms are collinear). At 280° we have again a broken H bond with antiparallel dipoles next to each other, but slightly different from the orientation at 60°. At 360°, H8 is the donor (as at zero). Because this rotation cannot be decoupled from mode 8, a 1D HR is problematic. The minima at 0 and 140° are topologically identical; still a symmetry number of 2 cannot be assigned because the 240° minimum is different from the others. In the future, we will investigate these modes as 2D HRs to test the present approximations. Mode 8, Donor H Rotation. A scan of D(7,3,2,1) with D(7,2,3,4) frozen yields this motion. Starting from 0° and increasing the dihedral leads to bond breakage and termination of the scan when atoms 2, 3, and 7 are collinear. Scanning in the other direction leads to coupling with the acceptor rotation;

because the distance of the two monomers increases during this motion, an overlap of atoms 3 and 7 is thus avoided. For a movement in the direction of the arrows, i.e., an increase in A(3,2,7), we fixed the dihedrals D(1,2,7,9) and D(4,2,7,9) to avoid (coupling with) a hydrogen rotation that could make donor and acceptor exchange their roles. This motion leads also to hydrogen bond breakage. Because an angle cannot be larger than 180° by definition, we had to follow the motion beyond zero and 180° to continue the scan and subtract the actual bond angle from 360° to describe the full motion. The rotation could also be called “donor rotation 2” because it is the second rotation of the donor. It has mainly technical reasons that we treat it as bond bend instead of as a torsion because there is no atom to define the required dihedral angle properly. For close-to-collinear arrangements, the geometry optimizations did not converge. Due to these missing points and due to its flatness, the potential fit contains spurious minima at high energies. However, the high-energy states resulting from these minima do not affect the partition functions at the temperatures we are interested in. Mode 5, Donor Methyl Rotation. This scan (D(1,2,3,4)) can easily be performed in one pass. We like to note the symmetry number of three and the fact that the barrier of the rotation in the dimer is lowered compared to the monomer (Figure 1 and Figure 6). At ambient temperature, this leads to a negative contribution to the entropy of dimerization, even though the correction from HO to HR is positive. Mode 6, Hydrogen Bond Stretch. The scan of the OO distance is again straightforward. A Morse potential was fitted and evaluated as described in the Methods. As can be seen from Figure 7, the fit is accurate at energies up to 15 kJ/mol and becomes inaccurate at higher energies and larger distances. This I



Article

Figure 9. Mode 8 of the methanol dimer: Donor H rotation.

but the wrong symmetry of the motion. As noted above, a coupled treatment especially including modes 7 and 8 is highly interesting but presently not possible. Thermodynamic Results for the Dimer and Comparison with Experimental Data. From the discussion of the individual modes we have learned that those modes that have several minima not related by symmetry contribute most to the entropy change from the RRHO model, cf. Table 5. This means that modes 2, 7, and 8 are the most relevant ones, but the others also contribute more than 7 J/(mol K). Therefore, all these modes have to be considered. Table 6 shows that

hence the maximum around 300° is again one with antiparallel dipoles next to each other (like the 60° maximum in mode 7). The shapes of both peaks are not identical because different degrees of freedom had to be frozen in both scans. At 240°, H8 is the donor and after exchanging topologically identical nuclei the minimum is a mirror image of the 0° minimum. If we now continue rotation of H1 by scanning D(7,3,2,1), the maximum at 180° is identical to that of mode 7 at the same angle. A further rotation leads to another exchange of donor and acceptor to recover the original arrangement; cf. Figure 9. As indicated above, the 1D treatments of mode 7 and 8 may be a crude approximation. Because the 2 minima in the 2D plot are distinct (mirror images of each other after exchange of topological identical nuclei), we do not apply a symmetry number to these motions. However, the minima found in mode 2 appear also in modes 7 or 8. Therefore, we assign an overall topological symmetry number of 245,46 to the system. In case it may seem counterintuitive to treat these two mathematically distinct geometrical arrangements as identical, we remind the reader that a partition function for two unbound monomers contains a factor of 2! in the denominator because we integrate over all mathematically possible arrangements and have to account for indistinguishability. In the dimer case, we do exactly the same. One of the reviewers raised the valid question whether the freezing of various dihedral angles leads to too-high barriers and therefore to wrong thermochemistry. Additional studies of the effect of freezing modes showed, however, that this is not largely the case. Without freezing any dihedrals in mode 1, the barrier in this mode remains the same, but a new transition state connecting the minima of modes 1 and 3 appears that is roughly 1.7 kJ/mol high. For modes 2, 4, 5, and 6 no degrees of freedom were frozen. When D(3,2,7,9) and D(3,2,1,7) are not frozen during the scan of mode 3, the motion changes into that of mode 1 for a large part of the rotation. Obviously, the barrier will be much lower, but mode 1 would effectively be counted twice in that case. Therefore, we consider the present treatment more realistic than the unfrozen one. For mode 7, the highest barrier is substantially reduced from 12 to 5 kJ/mol. However, the reduced barrier height comes from the methyl group orienting itself always in the lowest-energy position with respect to the rotating hydrogen atom. Because the moment of inertia of the methyl group is much larger than that of an H atom, we doubt that this will always happen and expect the effective barrier to lie between the two extremes, probably closer to the frozen one. Mode 8 would become identical to mode 2 if D(7,2,3,4) was not frozen leading to a similar barrier

Table 6. Thermochemistry for Methanol Dimerizationa ΔH, kJ/mol dimerization in the RRHO approximation dimerization with 8 HRs experimental dimerization (range) experimental dimerization (median)

ΔS, J/(mol K)

K

−17.2

−93.6

0.013

−15.9

−74.7

0.078

−13.5 to −18.0

−69.0 to −104.1

−14.6

−72.7

0.072

ΔH and ΔS are the differences between the enthalpy and entropy at 298 K, respectively, of the dimer and two monomers. Energy calculations have been performed at the RIMP2/aug-cc-pV(TQ)Z level. Exp. enthalpy and entropy data are from ref 2 and the equilibrium constant is from density data from ref 4. a

experimental results suggest a large possible range for the enthalpy and entropy of dimerization compatible with the RRHO and the HR results. However, the HR results agree much better with most of the experimental data. In addition to that, the constant of dimerization has much lower experimental uncertainties than enthalpy and entropy have. This property shows clearly that the HR results are much more accurate than the RRHO results and that they fulfill the accuracy requirements to calculate PC-SAFT association parameters. Thermodynamics for the Trimers, the c-Tetramer and c-Pentamer. In the trimers one would expect 15 (=3 + 2·6) modes that require anharmonic treatment each, in the cyclic tetramer 22 modes and in the pentamer even 29. To avoid a huge computational effort, we assumed that the thermodynamic HR corrections for modes of the dimer whose visualized normal modes look similar to modes of larger clusters can simply be added to the RRHO results of these clusters. The J



Article

Table 7. Assignment of Modes from the Dimer or Cyclic Tetramer to Similar Modes in Other Clustersa methanol cluster

linear trimer

1: OHO rotation 2: donor rotation 3: acceptor rotation 4: COO bond bend 5: donor methyl rotation 6: OO stretch 7: acceptor H rotation 8: donor H rotation uuud-to-udud 4-mer rotation CH3 4-mer rotation

2 1 2 2 3 2 1 2

× × × × × × × ×

cyclic trimer

(1,2) (3) (6,7) (4,5) (8−10) (11,12) (13) (14,15)

cyclic tetramer

cyclic pentamer

3 × (4,5,9)

6 × (6,7,11−14)

8 × (6,7)

3 × (10−12)

4 × (15−18)

6 × ()b

3 × (13−15) 3 × (1−3) 3 × (6−8)

4 × (19−22) 4 × (1−4) 4 × (5,8−10)

5 × ()b 5 × (1−5) 5 × (8−12)

a The numbers in parentheses denote the normal modes that are replaced by the dimer or tetramer modes, respectively. bIn the cyclic pentamer, modes 6−18 are linear combinations of CH3 rotations and OOC bond bends with varying contribution of each mode. We assume that there are five independent methyl rotations and attribute all other modes of this kind to the COO bond bend.

Figure 10. Left side: mode 1 of the cyclic methanol tetramer. This mode is a linear combination of rotations around the OHO bonds of all four monomers in the tetramer. Right side: scan of the rotation of one monomer in the tetramer around the OHO bond to its neighbor (up−down−up− down (udud) to up−up-up−down (uuud) transition).

and c-pentamer based on the RRHO results and HR corrections discussed before. Compared to the RRHO model, the entropy is considerably larger (less negative) for all species, as expected. Though the enthalpy is also less negative for most species due to the consideration of higher-energy local energy minima or almost-minima for most species, the c-tetramer is the only species with a more negative enthalpy. The reason for this behavior is that the uudd conformer was already considered in the RRHO treatment, and this treatment seems to overestimate its influence compared to the more accurate HR treatment. Table 9 contains also experimental data available for the trimer and tetramer. In the experimental studies, a differentiation between linear and cyclic species was not made. A comparison suggests that our calculations underestimate the association to trimers and tetramers considerably. However, most of the experimental data were derived from heat capacity and vapor density data using a cluster model neglecting pentamers. As we show in the next section, pentamers are very relevant and our model, which is fully ab initio and does not contain any assumptions on the existing species, agrees also very well with measured densities. Therefore, we expect that our model is more reliable and accurate than the existing data derived from experimental data by using additional assumptions for their analysis. Cluster Model for Methanol. The data in Table 10 show that the dimer is the overwhelmingly dominating cluster at low

modes that we actually transferred in this way from one species to another are compiled in Table 7. Additionally, we scanned one of the four equivalent methyl rotations in the cyclic tetramer (not shown) and one of the 4 similar vibrations that lead to an up−down transition of a methyl group (up−down−up−down, udud, to up−up-up− down, uuud, transition, Figure 10) to account for modes not present in the dimer. The contributions of individual modes for these clusters can be found in Table 8. By counting the udud to Table 8. Thermodynamic Corrections for the Explicitly Treated Hindered Rotations in the Cyclic Methanol Tetramer at 298 K rotation

ΔH, kJ/mol

ΔS, J/(mol K)

methyl rotation udud-to-uuud transition

−1.1 1.8

1.8 17.6

uuud transition four times, we generate 16 entropy-increasing plateaus for the tetramer of which only 6 are distinct after consideration of topological symmetry. However, exactly those conformations that have (three additional) topologically equivalent conformations have a topological rotational symmetry number of unity instead of four like the most stable tetramer. Therefore, these effects cancel exactly and it is not necessary to consider these topological equivalences. Table 9 shows the final thermochemical results for trimers, c-tetramer K



Article

Table 9. Thermochemical Data (Enthalpy and Entropy per Bond) of Relevant Methanol Clusters at the HR+RIMP2/aug-ccpV(TQ)Z//B3LYP/TZVP Level (Extrapolations) at 298 K and Equilibrium Constants at 373 Ka

a

methanol cluster

ΔH, kJ/mol

ΔS, J/(mol K)

ΔG, kJ/mol

dimer l-trimer c-trimer trimer, experimental result 1 and 2 c-tetramer c-tetramer, experimental range c-tetramer, experimental median c-pentamer

−16.0 −22.2 −18.7 −6.4, −17.4 −26.9 −18.8 to −25.9 −23.6 −28.2

−74.4 −96.1 −72.4 −24.8, −61.6 −91.6 −81.7 to −85.0 −83.4 −95.1

6.2 6.5 2.9 0.1, 1.0 0.4 −0.6 to +0.3 −0.1 0.1

K(373 K)

0.00008 0.00024−0.00026 0.00025

Experimental data are from ref 2.

Table 10. Pressure Dependent Mole Fraction of Relevant Methanol Clusters at the HR+RIMP2/aug-cc-pV(TQ)Z//B3LYP/ TZVP Level (Extrapolations) at 298 K and Ratio of Predicted and Experimentally Observed4 Density Change by Association Compared to the Ideal Gas (ρmodel − ρideal)/(ρexperimental − ρideal) p, bar 0.005 0.050 0.100 0.169

x2 4.1 4.1 8.0 1.3

× × × ×

10−4 10−3 10−3 10−2

x3,l 1.4 1.3 5.3 1.4

× × × ×

10−7 10−5 10−5 10−4

x3,c 7.9 7.8 3.1 8.1

× × × ×

x4,c

10−7 10−5 10−4 10−4

6.4 6.3 4.9 2.3

x5,c

10−8 10−5 10−4 10−3

5.0 4.9 7.7 5.3

× × × ×

10−10 10−6 10−5 10−4

density change ratio 1.08 1.15 1.07 0.89

obtain accurate entropies with ab initio methods. An interesting observation is that the intermolecular rotations do couple very strongly, but nevertheless our 1D HR results agree with the experimental ones within 1 kJ/mol for the clusters for which accurate experimental data are available. Whether this is by chance or generally valid needs to be investigated in the future. Also, the possibility to utilize the present results to compute PC-SAFT association parameters will be investigated in the future.

densities. At saturation pressure, 0.169 bar, the dimer still has the highest concentration, but especially the cyclic tetramer becomes also relevant. At this pressure, formation of the dimer captures only 50% of the density effect, demonstrating the relevance of the larger clusters. The predicted density change is always within 20% of the experimental one, corresponding to an error of 0.5 kJ/mol in ΔG. This error is well within the accuracy required for use in PC-SAFT. The composition of data in Table 10 together with the thermochemical data of Table 9 demonstrate that the average enthalpy and entropy of a hydrogen bond in methanol will be strongly temperature and density dependent. This is due to cooperativity and opposed to the usual assumption of state independent and hence effective association parameters in SAFT EOSs. Therefore, an extension of Wolbach and Sandler’s method of coupling both models at zero density seems necessary but is beyond the scope of the present contribution. In contrast to our results for the gas phase, chains seem to dominate in liquid methanol according to Huang et al.5 This may be due to unfavorable packing effects of the cyclic tetramers. In the gas phase, cyclic tetramers are more relevant than linear ones according to the same work, in agreement with our findings. Laksmono et al.1 tried to model the cluster formation observed in jet condensation experiments using various cluster sizes but failed to reproduce their experimental results. Although the cyclic pentamer was absent in all of their models taken from other works, they suggested that it may play an important role in the nucleation process. Our results corroborate their conclusion, because the pentamer can become the most important cluster for supersaturated gas densities at low temperature according to our findings.

■

× × × ×

■

ASSOCIATED CONTENT

S Supporting Information *

Energies and thermodynamic contributions of methanol clusters. Geometries of all used stable species. This material is available free of charge via the Internet at http://pubs.acs. org/.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: +49 (0)241 8098174. Fax: +49 (0)241 8092255. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors thank the Cluster of Excellence Tailor Made Fuels f rom Biomass, which is funded by the Excellence Initiative by the German federal and state governments to promote science and research at German universities for financial support and G. Sadowski and K. Albers for inspiring discussions regarding this work.

■

CONCLUSION

We have calculated accurate thermochemical properties of the most relevant methanol clusters in the gas phase. Besides large basis sets, counterpoise corrections and reasonable methods for electron correlation, the HR treatment was applied to methanol clusters for the first time and turned out to be essential to

REFERENCES

(1) Laksmono, H.; Tanimura, S.; Allen, H. C.; Wilemski, G.; Zahniser, M. S.; Shorter, J. H.; Nelson, D. D.; McManus, J. B.; Wyslouzil, B. E. Monomer, Clusters, Liquid: an Integrated Spectroscopic Study of Methanol Condensation. Phys. Chem. Chem. Phys. 2011, 13, 5855−5871.

L



Article

(25) Frisch, M. J.; et al. Gaussian 09 Revision A.1; Gaussian Inc.: Wallingford, CT, 2009. (26) Irikura, K. K.; Johnson, R. D., III; Kacker, R. N. Uncertainties in Scaling Factors for ab initio Vibrational Frequencies. J. Phys. Chem. A 2005, 109, 8430−8437. (27) Tsuzuki, S.; Uchimaru, T.; Matsumura, K.; Mikami, M.; Tanabe, K. Effects of Basis Set and Electron Correlation on the Calculated Interaction Energies of Hydrogen Bonding Complexes: MP2/cc-pV5Z Calculations of H2O−MeOH, H2O−Me2O, H2O−H2CO, MeOH− MeOH, and HCOOH−HCOOH complexes. J. Chem. Phys. 1999, 110, 11906−11910. (28) Feyereisen, M.; Fitzgerald, G.; Komornicki, A. Use of Approximate Integrals in ab initio Theory. An Application in MP2 Energy Calculations. Chem. Phys. Lett. 1993, 208, 359−363. (29) Weigend, F.; Haser, M.; Patzelt, H.; Ahlrichs, R. RI-MP2: Optimized Auxiliary Basis Sets and Demonstration of Efficiency. Chem. Phys. Lett. 1998, 294, 143−152. (30) Boys, S. F.; Bernardi, F. The Calculation of Small Molecular Interactions by the Differences of Separate Total Energies. Some Procedures with Reduced Errors. Mol. Phys. 1970, 19, 553−565. (31) Dunning, T. H., Jr. Gaussian Basis Sets for Use in Correlated Molecular Calculations. I. The Atoms Boron through Neon and Hydrogen. J. Chem. Phys. 1989, 90, 1007−1023. (32) Bak, K. L.; Halkier, A.; Jørgensen, P.; Olsen, J.; Helgaker, T.; Klopper, W. Chemical Accuracy from ‘Coulomb Hole’ Extrapolated Molecular Quantum-Mechanical Calculations. J. Mol. Struct. 2001, 567, 375−384. (33) McQuarrie, D. A. Statistical Thermodynamics; Harper & Row: New York, 1999. (34) Lucas, K. Molecular Models for Fluids; Cambridge University Press: Cambridge, U.K., 2007. (35) Foresman, J. B.; Frisch, A. Exploring Chemistry with Electronic Structure Methods; Gaussian Inc.: Wallingford, CT, 1993; pp 141−143. (36) Pfaendtner, J.; Yu, X.; Broadbelt, L. J. The 1-D Hindered Rotor Approximation. Theor. Chem. Acc. 2007, 118, 881−898. (37) Szakacs, P.; Csontos, J.; Das, S.; Kallay, M. High-Accuracy Theoretical Thermochemistry of Atmospherically Important Nitrogen Oxide Derivatives. J. Phys. Chem. A 2011, 115, 3144−3153. (38) Ghysels, A.; Verstraelen, T.; Hemelsoet, K.; Waroquier, M.; Van Speybroeck, V. TAMkin: A Versatile Package for Vibrational Analysis and Chemical Kinetics. J. Chem. Inf. Model. 2010, 50, 1736−1750. (39) Barone, V. Vibrational Zero-Point Energies and Thermodynamical Funktions beyond the Harmonic Approximation. J. Chem. Phys. 2004, 120, 3059−3065. (40) Lynch, V.; Mielke, S.; Truhlar, D. High-Precision Quantum Thermochemistry on Nonquasiharmonic Potentials: Converged PathIntegral Free Energies and a Systematically Convergent Family of Generalized Pitzer−Gwinn Approximations. J. Phys. Chem. A 2005, 109, 10092−10099. (41) Ellingson, B. A.; Lynch, V. A.; Mielke, S. L.; Truhlar, D. G. Statistical Thermodynamics of Bond Torsional Modes: Tests of Separable, Almost-Separable, and Improved Pitzer−Gwinn Approximations. J. Chem. Phys. 2006, 125, 084305. (42) Van Speybroeck, V.; Vansteenkiste, P.; Van Neck, D.; Waroquier, M. Why Does the Uncoupled Hindered Rotor Model Work Well for the Thermodynamics of n-Alkanes? Chem. Phys. Lett. 2005, 402, 479−484. (43) Vansteenkiste, P.; Verstraelen, T.; Van Speybroeck, V.; Waroquier, M. Ab initio Calculation of Entropy and Heat Capacity of Gas-Phase n-Alkanes with Hetero-Elements O and S: Ethers/ Alcohols and Sulfides/Thiols. Chem. Phys. 2006, 328, 251−258. (44) Atkins, P. W. Physical Chemistry; Oxford University: Oxford, U.K., 1994. (45) Fernández-Ramos, A.; Ellingson, B. A.; Meana-Pañeda, R.; Marques, J. M. C.; Truhlar, D. G. Symmetry Numbers and Chemical Reaction Rates. Theor. Chem. Acc. 2007, 118, 813−826. (46) Gilson, M. K.; Irikura, K. K. Symmetry Numbers for Rigid, Flexible, and Fluxional Molecules: Theory and Applications. J. Phys. Chem. B 2010, 114, 16304−16317.

(2) Curtiss, L. A.; Blander, M. Thermodynamic Properties of GasPhase Hydrogen-Bonded Complexes. Chem. Rev. 1988, 88, 827−841. (3) Haase, R.; Tillmann, W. Mixing Properties of the Liquid Systems Methanol + 2-Propanol and 1-Propanol + 2-Propanol. Z. Phys. Chem. 1995, 192, 121−131. (4) de Reuck, K.; Craven, R. Methanol, International Thermodynamic Tables of the Fluid State  12; IUPAC, Blackwell Scientific Publications: London, 1993. (5) Huang, Z.; Yu, L.; Dai, Y. Combined DFT with NBO and QTAIM Studies on the Hydrogen Bonds in (CH3OH)n (n = 2 − 8) Clusters. Struct. Chem. 2010, 21, 565−572. (6) Gonzalez, L.; Mo, O.; Yanez, M. Density Functional Theory Study on Ethanol Dimers and Cyclic Ethanol Trimers. J. Chem. Phys. 1999, 111, 3855−3861. (7) Buck, U.; Siebers, J.; Wheatley, R. Structure and Vibrational Spectra of Methanol Clusters from a new Potential Model. J. Chem. Phys. 1998, 108, 20−32. (8) Muñoz-Caro, C.; Nino, A. Effect of Anharmonicities on the Thermodynamic Properties of the Water Dimer. J. Phys. Chem. A 1997, 101, 4128−4135. (9) Sum, A.; Sandler, S. Ab initio Calculations of Cooperativity Effects on Clusters of Methanol, Ethanol, 1-Propanol and Methanethiol. J. Phys. Chem. A 2000, 104, 1121−1129. (10) Gross, J.; Sadowski, G. Perturbed-Chain SAFT: An Equation of State Based on a Perturbation Theory for Chain Molecules. Ind. Eng. Chem. Res. 2001, 40, 1244−1260. (11) Van Nhu, N.; Singh, M.; Leonhard, K. Quantum Mechanically Based Estimation of Perturbed-Chain Polar Statistical Associating Fluid Theory Parameters for Analyzing Their Physical Significance and Predicting Properties. J. Phys. Chem. B 2008, 112, 5693−5701. (12) Wolbach, J. P.; Sandler, S. I. Using Molecular Orbital Calculations to Describe the Phase Behavior of Hydrogen-Bonding Fluids. Ind. Eng. Chem. Res. 1997, 36, 4041−4051. (13) Becke, A. D.; Density-Functional Thermochemistry, I. The Effect of the Exchange-only Gra-dient Correction. J. Chem. Phys. 1992, 96, 2155−2160. (14) Becke, A. D. Density-Functional Thermochemistry. II. The Effect of the Perdew-Wang Generalized-Gradient Correlation Correction. J. Chem. Phys. 1992, 97, 9173−9177. (15) Becke, A. D. Density-Functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648−5652. (16) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle− Salvetti Correlation-Energy Formula into a Functional of Electron Density. Phys. Rev. B. 1988, 37, 785−789. (17) Zhao, Y.; Truhlar, D. G. Benchmark Databases for Nonbonded Interactions and their Use to Test Density Functional Theory. J. Chem. Theory Comput. 2005, 1, 415−432. (18) Santra, B.; Michaelides, A.; Scheffler, M. On the Accuracy of Density-Functional Theory Exchange−Correlation Functionals for H Bonds in Small Water Clusters: Benchmarks Approaching the Complete Basis Set Limit. J. Chem. Phys. 2007, 127, 184104−1−−9. (19) Sousa, S. F.; Fernandes, P. A.; Ramos, M. J. General Performance of Density Functionals. J. Phys. Chem. A 2007, 111, 10439−10452. (20) Vasiltsova, T.; Heintz, A. New Statistical Mechanical Model for Calculating Kirkwood Factors in Self-Associating Liquid Systems and its Application to Alkanol plus Cyclohexane Mixtures. J. Chem. Phys. 2007, 127, 114501. (21) Dyczmons, V. Dimers of Ethanol. J. Phys. Chem. A 2004, 108, 2080−2086. (22) Boyd, S. L.; Boyd, R. J. A Density Functional Study of Methanol Clusters. J. Chem. Theory Comput. 2007, 3, 54−61. (23) Schaefer, A.; Huber, C.; Ahlrichs, R. Fully Optimized Contracted Gaussian-Basis Sets of Triple Zeta Valence Quality for Atoms Li to Kr. J. Chem. Phys. 1994, 100, 5829−5835. (24) Weigend, F.; Ahlrichs, R. Balanced Basis Sets of Split Valence, Triple Zeta Valence and Quadruple Zeta Valence Quality for H to Rn: Design and Assessment of Accuracy. Phys. Chem. Chem. Phys. 2005, 7, 3297−3305. M



Article

(47) Heidemann, R. A.; Prausnitz, J. M. A van der Waals-Type Equation of State for Fluids with Associating Molecules. Proc. Natl. Acad. Sci. 1976, 73, 1973−1976. (48) Ohno, K.; Shimoaka, T.; Akai, N.; Katsumoto, Y. Relationship between the Broad OH Stretching Band of Methanol and HydrogenBondind Patterns in the Liquid Phase. J. Phys. Chem. A 2008, 112, 7342−7348. (49) Economou, I. G.; Donohue, M. D. Chemical, Quasi-Chemical and Perturbation Theories for Associating Fluids. AIChE J. 1991, 37, 1875−1894. (50) Gurvich, L. V. Reference Books and Data Banks on the Thermodynamic Properties of Individual Substances. Pure Appl. Chem. 1989, 61, 1027−1031. (51) Zheng, J.; Yu, T.; Papajak, E.; Alecu, I. M.; Mielke, S. L.; Truhlar, D. G. Practical Methods for Including Torsional Anharmonicity in Thermochemical Calculations on Complex Molecules: The Internal-Coordinate Multi-Structural Approximation. Phys. Chem. Chem. Phys. 2011, 13, 10885−10907.

N


Ab Initio Calculations of Thermochemical ... - ACS Publications

Recommend Documents