Equilibrium and Rate Constants, and Reaction Mechanism of the HF

Article pubs.acs.org/JPCA

Equilibrium and Rate Constants, and Reaction Mechanism of the HF Dissociation in the HF(H2O)7 Cluster by ab Initio Rare Event Simulations Alin Marin Elena,† Simone Meloni,*,†,‡,∥ and Giovanni Ciccotti†,§ †

School of Physics, Room 302 UCD-EMSC, University College Dublin, Belfield, Dublin 4, Dublin, Ireland CINECA, Via dei Tizii 6, 00185 Roma, Italy § Dipartimento di Fisica and CNISM, Università La Sapienza, P. le A. Moro 5, 00185 Rome, Italy ‡

S Supporting Information *

ABSTRACT: We perform restrained hybrid Monte Carlo (MC) simulations to compute the equilibrium constant of the dissociation reaction of HF in HF(H2O)7. We find that the HF is a stronger acid in the cluster than in the bulk, and its acidity is higher at lower T. The latter phenomenon has a vibrational entropic origin, resulting from a counterintuitive balance of intra- and intermolecular terms. We find also a temperature dependence of the reactions mechanism. At low T (≤225 K) the dissociation reaction follows a concerted path, with the H atoms belonging to the relevant hydrogen bond chain moving synchronously. At higher T (300 K), the first two hydrogen atoms move together, forming an intermediate metastable state having the structure of an eigen ion (H9O4+), and then the third hydrogen migrates completing the reaction. We also compute the dissociation rate constant, kRP. At very low T (≤75 K) kRP depends strongly on the temperature, whereas it gets almost constant at higher T’s. With respect to the bulk, the HF dissociation in the HF(H2O)7 is about 1 order of magnitude faster. This is due to a lower free energy barrier for the dissociation in the cluster.

1. INTRODUCTION Dissociation of (weak) acids and the correlated phenomenon of proton transport are among the most important processes for (bio)chemistry in bulk water solution, atmospheric and environmental chemistry (see ref 1 for a perspective on this subject), and chemistry in confined systems (channels, interfaces, etc.; see, for example, refs 2 and 3). The study of these processes in bulk solutions has been the subject of many experimental and computational investigations (see, for example, refs 4−12). Recent progresses in experimental techniques have fostered the study of acid−water clusters, AH(H2O)n, both experimentally and computationally (see ref 13 for an exhaustive review on the subject). Concerning the acid dissociation, apart some notable exceptions (see, for example, refs 14−16), the main research objective has been to identify the minimum degree of solvation (i.e., minimum number of water molecules) allowing this process.17−28 Most of the computational investigations have been performed at zero temperature, and the identification of the dissociation path and rate, and their dependence on the temperature, have received little attention. This is surprising considering the wide range of relevant temperatures, which goes from ∼1 K, reached in the recent experiments of Gutberlet et al.,15 to up to about 200− 250 K, the typical range of atmospheric conditions. As a first step in the attempt to close this gap, we investigate the dissociation of HF in one of the forms of the cubic HF(H2O)7 cluster (Figure 1) in the range 25−300 K. In particular, we © 2013 American Chemical Society

Figure 1. Two configurations belonging to the (A) initial and (B) final states of the HF dissociation in HF(H2O)7. Ball and stick representation is used for emphasizing the molecules involved in the process. Stick representation is used for the others. Hydrogen bonds relevant to the reaction are represented by blue dashed lines (the others are in red). Fluorine is drawn in green, oxygen in red, and hydrogen in white. This color convention is used throughout the manuscript, apart in Figure 7.

focus on the most stable dissociated cubic cluster and the corresponding undissociated form.29 The reason to focus on Received: July 15, 2013 Revised: November 18, 2013 Published: November 20, 2013 13039

dx.doi.org/10.1021/jp406982h | J. Phys. Chem. A 2013, 117, 13039−13050

The Journal of Physical Chemistry A

Article

not affect our argument, but it is necessary for the interpretation of f k(z) as the derivative of Fk(z) = −β−1 ln[Ak(z)/A ]. Noting that limk→∞ exp[−βk/2(ϵ(r) − z)2]/ (2π/βk)1/2 = δ(ϵ(x) − z), in this limit we have f k(z) = ∇zFk(z) → −β−1∇z ln Pϵ(z). Here Pϵ(z) = ∫ dr δ(ϵ(r) − z) exp[−βV(r)]/A is the probability that ϵ(r) = z. Recalling that the (Landau) free energy of a CV is defined as F(z) = −β−1 ln Pϵ(z), we find that in the proper limit eq 1 is an estimate of the derivative of the free energy: ∇zF(z). The mean force in eq 1 can be estimated by RMC, and the free energy profile can be computed by its numerical integration (thermodynamic integration34). A comment is in order concerning the estimation of the statistical error on the free energy. Using the error propagation formula for an observable that is, essentially, the sum of a large number of terms (the mean force at 34 z-points) results in a large overestimate of the error. Here we follow a different approach. We split each computational sample into four sets, from which we compute four independent estimates of the mean force at each z-value, f ki (z) with i = 1−4. From them we compute, via thermodynamic integration, four free energy profiles. The error at each z-point is, then, obtained from the square root of the estimate of the variance of the free energy at each point as obtained from the four profiles. Considering that the four samples are independent, this number must be further divided by the square root of the number of samples, i.e., by 2. Quantum effects associated with the hydrogen atoms could be included in the calculations following an extension of the TAMC/RMC method recently proposed by Bonella and Vuilleumier.35,36 This method, which is based on a path integral formulation of quantum mechanics, will increase the computational cost of our simulation by 1 order of magnitude or more. Considering this fact, and that previous works concluded that quantum effects do not change the qualitative picture of proton transfer, dissociation and recombination reactions in bulk systems,4,5,10,33,37 here hydrogen atoms are treated as classical particles. The physical potential V(r) is derived from ab initio calculations. This means that at each random displacement of one atom, a Kohn−Sham self-consistent cycle must be performed, which makes the “conventional” Monte Carlo method computationally inefficient. We overcome this problem by using the hybrid Monte Carlo (hMC) method.38,39 In hMC the random displacement of an individual atom is replaced by a collective displacement of all the atoms according to a short MD trajectory started from the current configuration rO and random momenta (pO) extracted from a Maxwell distribution. The new configuration rN is accepted or rejected according to the probability α = min{1, exp[−β(H(rN,pN) − H(rO,pO)]}, where H(·, ·) is the (Born−Oppenheimer) nuclear Hamiltonian and pN is the vector of momenta after N MD steps. It is worth stressing that, at a variance with standard MC, in hMC the acceptance criterion takes into account also the initial and final kinetic energy. In the hMC the Hamiltonian governing the MD can be different from the one used in the acceptance test, the second one determining the probability density function (PDF) sampled by the simulation. We exploit this feature in our RMC simulation. We evolve rO and pO according to the physical potential V(r) and accept/reject the hMC move according to the biased potential Uk(r,z). In this way we can efficiently sample the conditional PDF mentioned above via

HF(H2O)7 is that this is the smallest hydrated HF cluster supporting dissociation.26 Studying acid dissociation (and proton transfer) in finite temperature clusters by atomistic simulations is difficult because this process is a rare event, i.e., a process occurring with a frequency that is too low to be studied by brute force approaches (“standard” molecular dynamics (MD) or Monte Carlo (MC)). In this work we use the temperature accelerated Monte Carlo (TAMC) method,30 which allows the reconstruction of the free energy of nonanalytical Collective Variables (CV; see section 2). We characterize the process according to its pKa (pKa = −log Ka, with Ka dissociation constant), reaction mechanism, and rate as a function of the temperature. In the specific case of HF, the interest for the pKa is due to its unusual dependence on T in bulk solutions: the acidity of HF increases at lower temperatures.31 This phenomenon has been attributed to the ability of the F− ion to “structure” the solution, thus reducing the entropy of the product. This seems not to be possible in the HF(H2O)7, considering that the reactant and product have a very similar structure. It is, thus, interesting to investigate the trend of pKa vs T for this system and, possibly, explain its atomistic origin. We use a novel CV based on the electronic structure of the system. The main advantage of this CV over geometrical ones used in previous works (A−−H+ distance or the coordination of A−)8,14−16,32 is that it is not based on any (crucial) a priori assumption on the dissociation path. Thus, for example, we will be able to assess whether the dissociation follows a step-like path, consisting of a series of proton migration events involving a pair of molecules at a time, or a concerted one, characterized by a collective migration of several H+ all together. The steplike path, which has been thought to be the diffusion path of H+ in bulk water, has recently been challenged by Parrinello and co-workers,33 who have proposed that self-dissociation of H2O (and possibly other weak acids) in bulk water proceeds according to a concerted mechanism. The remainder of the paper is organized as follows. In section 2 we summarize the computational method used to compute the free energy, describe our CV, give the details of the atomistic model and of the setup of the ab initio TAMC simulations. In section 3 we present and discuss our results. Finally, in section 4 we draw some conclusions.

2. THEORETICAL BACKGROUND AND COMPUTATIONAL DETAILS We perform free energy calculations using the restrained MC (RMC) approach, a specialized version of the TAMC simulation.30 In RMC, atoms are “evolved” according to the potential Uk(r,z) = V(r) + k/2(ϵ(r) − z)2, which is the sum of the physical (V(r)) and restraining (k/2(ε(r) − z)2) potentials. ϵ(r) is a suitable CV and z is its target value. Let us consider the following average: k

f (z ) = −

∫ dr k(ε(r) − z) exp[−βU k(r,z)]

A k (z ) ⎛ A k (z ) ⎞ ⎟ = −∇z β −1 ln⎜ ⎝ A ⎠

(1)

w h e r e A k (z ) = ∫ d r e x p [ − β U k ( r , z ) ] a n d A = ∫ dr exp[−βV(r)] is the canonical partition function of the real system. Because A is z-independent, its introduction does 13040



Article

apparently simpler. For example, the interpretation of the values of Δr, negative when the hydrogen is attached to one molecule, positive when it is attached to the other, and zero when it is equally shared, is simple only when the molecular species forming the hydrogen bond are the same (e.g., both water molecules), and are embedded in an equivalent environment (e.g., have the same degree of solvation). This is not the case for the problem at hand, in which we have HF and H2O molecules. Moreover, Δr is a “local” observable and does not take into account the state of the entire system, such as polarization of the “solvent”5 (the other water molecules of the cluster, in the present case) or the overcoordination of the donor/acceptor atoms, which have been identified as relevant phenomena in the dissociation of weak acids in bulk systems.33 The CV we use to study the dissociation of HF in the HF(H2O)7 cluster is derived from the observable ξ(r) introduced above. In the HF(H2O)7 structure considered in this work, the transformation from reactant to product consists in the transfer of a proton along the three H-bonds chain represented by the dashed blue segments in Figure 1. Let us denote by wix(x;r), with i = 1, 2, 3 and χ = α, β, the MLWFs associated with this chain (Figure 2). The index i identifies a specific H-bond in the chain, i = 1 denoting the F−H−O Hbond, and the index increasing along the chain. The index χ distinguishes between the two MLWFs associated with each Hbond, α being associated with the atom of the ith H-bond with the lowest counting along the H-bond chain starting from F, and β with the other (Figure 2). We can, thus, define ξi(r) = ⟨wiα|Hkswiα⟩ − ⟨wiβ|Hkswiβ⟩, which measures the “degree” of proton transfer in the ith H-bond along the relevant chain. ϵ(r) = Σ3i=1 ξi(r), which takes negative values in the reactant state and positive values in the product state, and then measures the total degree of proton transfer along the chain. Before closing this section, we must explain how to identify wiα/β(x;r). Let us consider the H-bond associated with a specific H atom. The two MLWFs associated with this H-bond are those having the centers ⟨x⟩j(r) the closest, in the Euclidean sense, to this H atom. We can assign the index α and β on the basis of the distance of the corresponding center from the O or F atoms involved in the H-bond. 2.2. Atomistic Model and Computational Setup. Calculations were performed using a customized version of the CP2K package,46,47 in which we implemented the hMC version of the TAMC method30 and the CV introduced in section 2.1. CP2K is a density functional theory (DFT) based Born−Oppenheimer molecular dynamics code. Only valence electrons are explicitly taken into account, and their interaction with the nuclei plus the associated (frozen) core electrons is described in terms of the Goedecker, Teter, and Hutter (GTH) pseudopotentials.48−50 Kohn−Sham orbitals are expanded on a localized (Gaussian) basis set, and a supplementary plane wave basis is used to expand the electron density.51−53 In this work we use the m-TZV2P (3s3p2d/s2p) Gaussian basis set54 (see also ref 55). In ref 29 it has been shown that this basis set is adequate to model the dissociation reaction studied here. The cutoff on the auxiliary planewave basis is set to ∼4000 eV. Calculations were performed using the Hamprecht et al. exchange−correlation functional HCTH120,56,57 which has been advocated to produce accurate structural results for pure and protonated water.58,59 Nevertheless, to the best of our knowledge the accuracy of this functional in modeling dissociation reactions has not been extensively investigated.

collective atomic displacements. It must be stressed that this approach is efficient in exploring the configurational space only when the CV is a smooth function of the atomic coordinate, i.e., when it does not change significantly during the short MD, so that Uk(rN,z) − Uk(r0,z) ∼ V(rN) − V(r0); otherwise, a large number of hMC move could be rejected. 2.1. Collective Variables To Study the Dehydrogenation Process. The CV we propose to use to study the HF dissociation in HF(H2O)7 is based on the expectation value of the single particle Kohn−Sham Hamiltonian (HKS) computed over a suitable set of localized orbitals, namely the maximally localized Wannier functions40,41 (MLWFs; see ref 42 for their definition), representing the covalent OH/FH bonds and the O/F− lone pairs involved in the relevant H-bond chain connecting H+ and F− (the H-bonds denoted by the dashed blue segments in Figure 1). Each H-bond has associated two MLWFs, characterized by a single particle density ρ(x;r) = |w(x;r)|2 (w(x;r) denotes a generic MLWF at the nuclear configuration r) localized around the corresponding covalent OH/FH bond or O/F− lone pair.43−45 In Figure 2 we show the

Figure 2. Same configurations of Figure 1 together with the centers of the relevant MLWFs.

“centers”, ⟨x⟩(r) = ⟨w(r)|xw(r)⟩ of the relevant MLWFs for two configurations belonging to the reactant (panel A) and product (panel B) states. How to identify these MLWFs out of the entire set of orbitals, which are in number equivalent to the Kohn−Sham orbitals of the system, will be explained at the end of this section. Here, we first concentrate on the definition of the CV. Let us denote by wα(x;r) and wβ(x;r) the two MLWF associated with a given H-bond. The expectation value of the single particle Kohn−Sham Hamiltonian over them, ⟨HKS⟩α(r) = ⟨wα(r)|HKSwα(r)⟩ and ⟨HKS⟩β(r), measures the strength of the corresponding “bond”: the more negative is ⟨HKS⟩(r), the stronger is the bond. In particular, ⟨HKS⟩(r) of the O/F−H covalent bond is more negative than the ⟨HKS⟩(r) of the lone pair involved in the given H-bond. Let us consider the observable ξ(r) = ⟨HKS⟩α(r) − ⟨HKS⟩β(r). In a configuration r′ in which the MLWF α corresponds to the covalent bond and the β to the lone pair, ξ(r′) < 0. In a configuration r″ in which the nature of wα(x) and wβ(x) is inverted, ξ(r″) > 0. When the proton is equally shared between the two molecules forming the H-bond, ⟨HKS⟩α = ⟨HKS⟩β and ξ(r) = 0. Thus, ξ(r) can be used to monitor the progress of the proton transfer within an H-bond. Simpler observables, such as Δr = rα − rβ, where rα and rβ are the distances between the hydrogen and the donor/ acceptor atoms forming the H-bond, have been used to monitor the proton transfer/acid dissociation process.8,32 These observables suffer from several limitations and are only 13041



Article

We performed a series of tests (described in detail in Supporting Information) showing that the mechanical energy profile along the dissociation of H3O+···H2O and HF···H2O clusters computed using HCTH120 exchange and correlation functional is consistent with those obtained from other largely used functionals (BLYP,60,61 B3LYP,62,63 and PBE64). We also tested HCTH120 against highly accurate quantum chemistry methods (CCSD(T)65−67) finding results in fairly good agreement. The HF(H2O)7 cluster was put in a cubic box of size 14.0 Å. The long-range interaction method for cluster systems of Martyna and Tuckerman is used to compute electron−electron and electron-ion (nuclei plus frozen core electrons) interactions.68 To study the effect of T on the equilibrium and rate constants, and reaction mechanism, the free energy profile of the CV defined above has been computed at five temperatures, 25, 75, 150, 225, and 300 K. The free energy curves have been obtained by numerical integration of the mean force computed at 34 z-points in the range ϵ ∈ [−16.6, +16.3] eV. The value of the coupling constant k of eq 1 has been set on the basis of the results of a convergence test on the dependence of the Fk(z) on this parameter at T = 75 K. Fk(z) has been calculated for k = 5, 25, 50, 100, 200 and it was found that it is already well converged for k = 50 (we consider the free energy converged when maxz|Fk′(z) − Fk″(z)| ≪ kBT, where k′ and k″ are two successive k values in the list above). Noticing that Fk(z) requires less RMC steps to converge for lower k values, we set k = 50 for all the temperatures. With this value of k, the estimation of the mean force at each z-point requires 24 000 hMC steps.

Figure 3. Free energy vs z at various temperatures. In the inset we report the same curve in kBT units to highlight the effect of the temperature on the free energy barrier. The arbitrary constant in the thermodynamic integration is set such that F(zr) = 0 (i.e., zr is the position of the low z-value free energy minimum for the reactant).

3. RESULTS AND DISCUSSION We start the presentation of our results by stressing that we did not observe any evaporation or change of conformation of our cluster even at 300 K. This may be surprising considering the relative low stability of hydrogen bonded clusters. We believe that this behavior is due to restraint forces, which limit the system to visit only states compatible with the condition ϵ(r) = z. Other cluster geometries might also be compatible with this condition; however, if the domains for which the condition is satisfied are disjoint, the system is forced to explore only the initial domain. Let us now move to our results and begin by analyzing some general features of the free energy profile at various temperatures (Figure 3). For T ≤ 225 K the system presents only two metastable states, corresponding to the undissociated state and the (dissociated) state in which the hydronium ion is on the opposite corner of the cluster with respect to F−, i.e., when the charged species are at the maximum distance (Figure 1). This picture is consistent with the T = 0 K findings of Xie et al.29 Consider the difference of internal energy between product and reactant, ΔE = ⟨H⟩zp − ⟨H⟩zr, with ⟨·⟩z denoting conditional ensemble average at z, and zp and zr are z-values of the product and reactant minimum of F(z). ΔE values calculated in the present work (ΔE ∼ 0.1−0.125 eV, Figure 4) are consistent with literature data (ΔE ∼ 0165 eV).29 At T = 300 K the system presents a shallow intermediate minimum at Z ∼ 6 eV. The transformation from a two-state to a three-state system with T is continuous. At lower T we observe a flex in the F(z) curve corresponding with the intermediate minimum. This flex becomes more evident at

Figure 4. ΔF, ΔE, and TΔS vs T.

higher T and, eventually, transforms into a minimum at T = 300 K. The atomistic origin of this change of shape of the F(z) curve with T will be explained in detail below, in the context of the discussion on the reaction mechanism. From the free energy curves of Figure 3 we compute the equilibrium constant of the dissociation reaction: Ka = Pε(zp)/ Pε(zr). In Figure 5 we report pKa vs T. As a reference, we also report experimental values obtained in bulk conditions (ref 31 and references cited therein). We notice that the acidity increases with the decrease of T, with the exception of the T = 25 K case. T = 25 K results are probably affected by artifacts related to the classical treatment of the H atoms. In fact, at this 13042



Article

controlling the dependency of the acidity of HF on T is the entropy. Because there is no difference between the configurational entropy of the reactant and product, we must conclude that the increase of HF acidity at low T is governed by the difference of vibrational entropy (Appendix A) between these two states. To confirm this conclusion, we compute the vibrational TΔS vs T in the harmonic approximation, i.e., from the normal modes of the minimum potential energy configurations corresponding to reactant and product, and compare it with the TΔS curve obtained from the free energy calculations. The harmonic vibrational TΔS is (Appendix B): 3N − 6

T ΔS hvib =

∏ νiR kBT log i3=N1− 6 P 2 ∏i = 1 νi

(2)

where νi is the vibrational frequency of the ith normal mode and R and P denote reactant and product, respectively. The set of vibrational frequencies {νi}i=1,3N−6 is determined using the “conventional” normal modes method, consisting in diagonalizing the mass-weighted potential energy Hessian matrix, Ĥ , of the reactant and product configurations. The Hessian matrix is computed within the finite difference approximations (Ĥ i,j = [Fi(r + δrj) − Fi(r + δrj)]/2|δrj|, where Fi(·) is the ith component of the force vector and δrj is a vector with nonzero value only in the jth element. The comparison between the TΔS curves obtained from free energy and normal modes calculations is shown in Figure 6. As a general remark, we

Figure 5. pKa vs T. For reference, we also report experimental data for HF dissociation in bulk water (see ref 31 and references cited therein).

temperature the de Broglie thermal wavelength of H is λ ∼ 3.5 Å, longer than the O−O/O−F distance (2.4−2.8 Å) between neighboring molecules in the cluster. At T = 75 K λ ∼ 2 Å, shorter than the O−O/O−F distances, and the classical approximation for the hydrogen atoms is less critical. In the following we discuss only the cases at T ≥ 75 K. The origin of low acidity of HF in bulk water, and of the temperature dependence of its pKa, has been longly debated. Consensus has been achieved that this is due to an entropic effect: F− is a “structure making” anion; i.e., it strengthens the H-bond network of the solution, thus reducing its configurational entropy (see Appendix A for a precise definition of configurational entropy). This results in a highly negative TΔS = T(Sp − Sr), which makes the HF dissociation a process thermodynamically unfavored.69 The entropic contribution to the free energy decreases at lower T, thus producing the observed increase of the acidity of HF. In the case of HF dissociation in HF(H2O)7, where we have only one metastable state in the reactant and product domain, the configurational ΔS is zero. Thus, we expect that HF is (i) a stronger acid than in the bulk and (ii) its strength is weakly dependent on the temperature. Our results only partly support this picture. In fact, although HF is stronger in the cluster than in bulk (Figure 5),70 the pKa still shows a marked T dependence. We analyze the origin of this dependence by calculating the various contributions to the ΔF, the internal energy and the entropic contributions. The definition of the variation of internal energy ΔE between product and reactant has already been given above; the entropic contribution to TΔF is computed via the relation TΔS = ΔE − ΔF. ΔF, ΔE, and TΔS vs T are reported in Figure 4. As a first remark, we notice that TΔS in the HF(H2O)7 cluster ([−0.01, −0.21] eV) is sizably smaller than in bulk water (∼−0.3 eV at 298 K).31 We also notice that the ΔE and TΔS trends are similar. However, the variation of TΔS with T is steeper than the variation of ΔE, resulting in a decrease of ΔF at low T (i.e., ΔF gets more negative), corresponding to an increase of pKa. In conclusion, as anticipated above, also in the cluster case the factor

Figure 6. Comparison between the TΔS obtained from free energy calculations and normal modes (see text). The dotted line represents the entropic contribution to the free energy coming from intramolecular modes, and the dashed line is relative to intermolecular modes.

notice that the curves are in qualitative agreement, confirming the hypothesis that the TΔS in our free energy calculations has a vibrational origin. It is, however, worth noticing that although the TΔS curves match well at low T, they diverge at higher temperatures. We believe that this difference is due to anharmonic effects, neglected in eq 2, which become more important at higher T’s. 13043



Article

values and, passing by green, become blue at high values of z. We notice that they move together; i.e., at intermediate z-values (green) they are all, approximately, midway along their path and complete their migration at almost the same value of z (blue). This indicates that the HF dissociation occurs via a synchronous migration of all the H atoms of the relevant Hbond chain. This mechanism is consistent with the one recently proposed by Parrinello and co-workers33 for the autodissociation of H2O in bulk water, who also proposed that this process is triggered by a shortening of the H-bond chain. To test whether the dissociation of HF in HF(H2O)7 is also triggered by a compression of the relevant H-bond chain, we measure its (Euclidean) length as a function of z (L(z)). L(z) is defined as ∑3i=1 di(z), where di(z) = |rαi (z) − rβi (z)|, i.e., the distance between the average positions of the heavy atoms (O and F) of the ith H-bond of the relevant chain (Figure 1). L(z) curves at various T are shown in Figure 8. In all cases we observe an

We further analyze the results of our normal-mode analysis to identify the origin of the negative TΔS. We are able to classify the normal modes in two sets: intramolecular modes, involving bonds stretching and bending (typically corresponding to the high frequency part of the spectrum), and intermolecular modes, corresponding to molecular librations and breathing modes of the cluster (low frequency part). The intra- and intermolecular TΔS contributions are shown as dotted and dashed lines in Figure 6, respectively. We find that the negative total TΔS is the result of two terms of opposite sign, the intramolecular positive and the intermolecular negative. The positive intramolecular entropy variation results from the presence in the product of the softer O−H covalent bonds of the H3O+ ion, replacing the stiffer O−H bonds of the former water molecule and the H−F bond of the undissociated acid of the reactant. The presence of two ions of opposite signs (H3O+ and F−) strengthen the H-bonds of the product, thus increasing the frequencies of the intermolecular modes and so reducing the entropy of the product. As discussed above, the combination of these effects results in an overall decrease of the vibrational entropy. Let us now move to the analysis of the mechanism of the HF dissociation and its dependence on the temperature. By mechanism we mean the atomistic path followed by the system along the transformation from reactant to product. This information is not readily available in our simulations as we do not generate reactive trajectories. Nevertheless, following the procedure described below, we are able to obtain relevant information about the mechanism. Let us define the average configuration at given ϵ(r) = z, r(z), as ⟨r⟩z. In the case of the reaction studied in this work, this is a meaningful quantity because the distribution of r at fixed z is monomodal, as resulted from a standard cluster analysis71 on the sampled configurations showing that they are well represented by a single configuration. r(z) is a parametric curve that can be seen as the path followed by the atoms along the reaction. In Figure 7 we show the “reactive path”, r(z), at T = 75 K. The color of the atoms involved in reaction denotes the z-value at which the corresponding r(z) is computed: atoms are in red at low z-

Figure 8. Length of the H-bond chain, L(z), vs z at the four temperatures considered in the text.

initial compression of the chain (−18 eV ≤ z ≤ −10 eV), followed by a phase in which the chain remains compressed. Finally, at z ∼ 0 eV the H-bond chain starts relaxing, and its length returns to large values (L(z) > 7.8 Å). This final phase is different for the T = 300 K case, in which we observe a sudden increase of L(z) at the z-value corresponding to the intermediate minimum of the free energy. After this, the Hbond chain shrinks again and, then, follows a trend analogous to that at the other temperatures. The above analysis confirms that also in the case of HF dissociation in HF(H2O)7 the reaction is triggered by a shortening of the H-bond chain along which the H+ defect moves. To analyze the reaction mechanism in more detail, we compute the curves ⟨ξi(r)⟩z vs z at various T (Figure 9). We recall that ⟨ξi(r)⟩z measures how much an H atom is tied up to the heavy atoms of an H-bond. As a general remark, we notice that the dissociation mechanism is not really synchronous, as was inferred by visual inspection of the r(z) “trajectory” (Figure 7). Indeed, the H atoms forming the first and second H-bond of the chain move first, followed by the H atom of the third Hbond. This is due (Figure 10) to an asynchronous compression of the H-bond chain, which takes place via an initial shortening

Figure 7. Dissociation path at 75 K. The oxygen and hydrogen atoms not involved in the dissociation process are shown at their initial configuration and drawn like in Figure 1. The other atoms are drawn in a trajectory-like plot: they are drawn at the conditional average positions rz for several z values going from the reactant (low z-values) to the product (high z-values). The color of these atoms indicates the corresponding z-value: red means low, green intermediate, and blue high. 13044



Article

Figure 9 we notice that the z value at which ⟨ξ3(r)⟩z starts to grow from its initial value increases with T. Also the lag between the z values at which ⟨ξ1/2(r)⟩z and ⟨ξ3(r)⟩z are equal to 0 increases with T. This indicates that at higher T the proton transfer of the first and second H-bond is completed before the third is started; i.e., a positive charge is accumulated on the second oxygen, O2 (see Figure 1 for the numbering convention of the O atoms). We now investigate the origin of the change of the shape of the free energy with T. We recall that at lower temperatures, T ≤ 225 K, the free energy presents two minima, whereas at T = 300 K it is characterized by three minima, the intermediate one being separated from the product by a small barrier. We characterize the reaction at various T by monitoring the distance between O2 and the oxygen of the neighboring water molecules (inset of Figure 11) vs z. We also measure the

Figure 9. ξi(z) vs z at the four temperatures considered in the text. Black arrows highlight when ε1(z) = 0, i.e., when the H atom is perfectly shared between F and O1. Red arrows highlight the same event for the H atoms between O1 and O2. The low z-value blue arrows highlight when the H atom between O2 and O3 starts changing its original state, i.e., start migrating from O2 to O3. Finally, the high zvalue blue arrows highlight when the H atom between O2 and O3 is perfectly shared between the two.

Figure 11. Length of the O−H (left) and O−O (right) bonds formed by O2 with the neighboring water molecules vs z at 225 (bottom) and 300 K (top). The image at the top-right shows the cluster with three arrows highlighting the relevant water molecules. The color of the arrows are consistent with the color of the curves in the graph. The two big arrows in the graphs show the z region of the intermediate metastable state.

distance between O2 and the hydrogen atoms forming the Hbonds with these water molecules. In Figure 11 we report these results for the T = 225 and 300 K cases. The curves at T < 225 K are similar to the 225 K case. In the region z ∼ 6 eV, corresponding to the intermediate metastable state, we observe a qualitative difference between the curves at the two temperatures. At 225 K two O−O (O−H) distances are shorter, 2.52 Å (longer 1.05 Å) than the third, 2.64 Å (0.99 Å). At 300 K the three distances are closer: two are 2.55 Å (1.05 Å), and the third is 2.60 Å (1.01 Å). The structure at 300 K looks like the well-known eigen ion (H9O4+)72 characterized by two slightly shorter (2.505 Å) and one longer (2.532 Å) Hbond.73 Our conclusion is that at high enough temperatures the larger fluctuations of the distances between the molecules forming the cluster permit the formation of the relatively stable H9O4+ ion.

Figure 10. di(z) vs z at all the temperatures considered in the text.

of the first and second H-bond followed (higher z values) by the third one. We will return to the analysis of the H-bond length vs z, and its dependence on T, later on. For the moment, we keep focusing on ⟨ξi(r)⟩z. Analyzing the results shown in 13045



Article

We now present the results of the calculation of the rate constant of HF dissociation reaction, kRP, and discuss its temperature dependence. We will use the transition state theory in the Vanden-Eijnden and Tal formulation74 (a similar formulation has been proposed also by other authors: Anderson,75 van Erp,76 and van Erp and Bolhuis77). The rate constant is the result of the product of two terms: (i) the exponential term exp[−ΔF*/kBT] (ΔF* = F(z*) − F(zr), z* is the z-value corresponding to the transition state), and (ii) the mean frequency of hopping between the reactant and product domains, κs. κs = ⟨ε̇(r) ζp/z*(r,p) ζp/z(r,−p)⟩z*, where ζp/z*(r,p) is the probability that a trajectory started at (r, p) will reach the product before crossing the transition state surface (ϵ(r) = z*), and ζr/p(r,−p) is the probability that the trajectory comes from the reactant rather than the product. The symbol ⟨·⟩z* indicates that the average is taken over the transition ensemble, i.e., over the probability density to be at r and p under the condition ϵ(r) = z*. The exponential term can be obtained from the free energy curves of Figure 3. The mean frequency of hopping is estimated using the formula κs = 1/ R/P N∑i=1,N ε̇(ri) χP(r(tp/z i *,ri,pi)) χR(r(ti ,ri,pi)) (see ref 78 and r/p * ≥ 0 and t Appendix C). Here tr/p i i < 0 are the times at which a trajectory started (t = 0) in ri and pi reaches either the product or the transition state (hyper)surface, and the reactant or product, respectively. The sum runs over N initial conditions belonging to the transition ensemble. ε̇(r) is estimated by the finite difference along the trajectory started from the transition ensemble: ε̇(r) = ϵ(r(h)) − ϵ(r(−h))/2h, where h is the MD time step. We cannot use the common formula ε̇(r) = ∇rϵ(r)ṙ because, as explained above, our CV is a nonanalytic function. Before presenting our results, we remark that at T = 300 K the system is characterized by two rate constants, one associated with the transition between the reactant and intermediate metastable state, and the second to the transition between the intermediate metastable state and the product. The second process is 3 orders of magnitude faster than the first: K1RP = 2.17 × 109 s−1 vs K1RP = 4.46 × 1012 s−1. This means that the limiting step is the first part of the process, which is the only one we discuss in the following. In Figure 12 we report the kRP as a function of the temperature. kRP grows exponentially with T in the range 25− 150 K and then remains almost constant at higher T (kRP = 2.2 × 109 s−1 at 150 K, 3.2 × 109 s−1 at 225 K, and 2.17 × 109 s−1 at 300 K). The rapid growth of kRP at low T is mainly due to the decrease of the free energy barrier (ΔF*/kBT), which passes from ∼30 to ∼7 in the range 25−150 K. Also the increase of κs with T contributes to increasing the rate constant but, being the growth of κs only a factor 5 in the 25−150 K range, its effect on kRP is much lower than the effect of the free energy barrier. The dissociation rate of the HF in HF(H2O)7 at 300 K is about 1 order of magnitude higher than experimental79 (7 × 107 s−1) and computational80 (1.5 × 108 s−1) values in bulk water. This is consistent with recent experimental results indicating that the water surface is more acid than the bulk.81 In fact, all the water molecules in our cluster, being undercoordinated, are surfacelike. The higher rate constant in the cluster is due to a lower free energy barrier (ΔF*/kBT ∼ 7.6 in the present calculations against ∼10.6 in ref 80), whereas the hopping frequency is almost the same in the cluster and bulk case (∼45 × 1012 s−1 against 6 × 1012 s−1). A higher free energy barrier in the bulk case is not unexpected. In fact, the formation of the transition state in this condition requires a significant structural change,

Figure 12. kRP vs T. Experimental79 (◆) and theoretical80 (■) data available in the literature are also reported. In the inset we also report ΔF* and κs vs T.

with an increase of the coordination number of fluorine from ∼3 to ∼5,80 whereas in the cluster it requires only a compression of the relevant H-bond chain.

4. CONCLUSION In this work we have studied the dissociation reaction of HF in the HF(H2O)7 cluster. We found that, like in the bulk, the HF is a weak acid, which increases its strength at lower temperatures. Also in cluster conditions the weakness of HF has an entropic origin, the entropy difference between the reactant and product being negative. At lower T the entropic contribution to the free energy is lower and the acid is stronger. We also analyzed the dissociation mechanism. We found that the process is “cooperative”, with all the H atoms of the relevant H-bond chain moving together from the reactant to the product configuration. This process is triggered by a compression of the H-bond chain. This finding is consistent with the observations of Parrinello and co-workers33 on the self-dissociation of water in bulk water. A closer inspection, however, reveals that the H atoms of the first and second Hbonds move first, followed by the third. We also found that at T = 300 K an intermediate state, corresponding to the formation of an Eigen ion, is formed. Concerning the rate constant of the reaction, we found that it strongly depends on the temperature at low T (≤75 K), whereas at T ≥ 150 K it is almost constant. Indeed, at T = 300 K we must consider two processes, the transformation from the reactant to the intermediate metastable state, and the transformation from the intermediate state to the product. However, the rate of the second process is 3 orders of magnitude higher than the first, thus the limiting step of the reaction results to be the first. The rate of dissociation of HF in the cluster is found to be about 1 order of magnitude higher than in the bulk. This is mostly due to a lower free energy barrier for the dissociation in the cluster. This is not surprising as the formation of the transition state in the cluster requires just a compression of the H-bond chain whereas in the bulk 13046



Article

TΔS = T(S2 − S1) and the partition function in the canonical ensemble:

case a complex rearrangement of the hydration shell of the HF is needed.80

■

T ΔS = −

APPENDIX A. DEFINITION OF CONFIGURATIONAL AND VIBRATIONAL ENTROPY The entropy is often split in several contributions. Here we consider the configurational and vibrational contributions, Sc and Sv. The explicit formulas of these two contributions in terms of the partition function can be derived according to the following argument. Let us consider a system in the microcanonical ensemble (this argument can be made general via the equivalence of ensembles) and assume that this system has N metastable states. Then, the partition function can be written as Q = ∫ dr dp δ(H(r,p) − E) = ∑i=1,N ∫ Ωidr ∫ dp δ(H(r,p) − E) where Ωi is the volume in the configuration space associated with the ith metastable state. The set {Ωi}i=1,N can be obtained, for example, from a Voronoi tessellation of the configuration space. The partition function can then be written Q = ∑i=1,N qi, where qi = ∫ Ωidr ∫ dp δ(H(r,p) − E) can be interpreted as a “local” partition function, relative to the ith metastable state. We can now multiply and divide the right hand side of this relation by N obtaining Q = Nq̃, where q̃ = 1/ N∑i=1,N qi is the average value of the “local” partition functions. Thus, the entropy can be written S = kB log Q = kB log N + kB log q̃ = Sc + Sv, where we defined Sc = kB log N and Sv = kB log q̃. The extension of the definitions of vibrational and configurational entropy to a two states system, having a reactant and a product domain, is straightforward. We first divide the entire configuration space into two parts, the reactants and products subspace. For example, according to the definition of reactant and product used in this work, we define the reactant configuration subspace as ΩR, {r|ϵ(r) ≤ z*}, and the product configuration subspace as ΩP {r|ϵ(r) > z*}. We, then, perform the same analysis described above within ΩR and ΩP. Thus, the configurational and vibrational entropies of the reactant (product) is defined as SRc = kB log NR (SPc = kB log NP) and SRc = kB log q̃R (SPv = kB log q̃P), where NR (NP) is the number of metastable states in the reactant (product) domain and q̃R = 1/NR∑i=1,NRqi (q̃R = 1/NP∑i=1,NPqi) is the average value of the corresponding “local” partition functions.

Z Z d log 2 + β −1 log 2 dβ Z1 Z1

(3)

where β is the inverse temperature and Z is the partition function. Because the number of particles and their masses are the same in both states, the only contribution to TΔS comes from the configurational part of the partition function (configuration integral), Q. Thus, in the following we will focus only on this term: Q = ∫ dr exp[−1/2βrTK̂ r]. We now perform a variable transformation q = Û Tr, where Û is a unitary matrix. We chose Û T such that the harmonic constant matrix in the q coordinates, Â = Û TK̂ Û , is diagonal: 1/2rTK̂ r = ∑jÂ jj/2qj2. Of course, 3, for an isolated system, or 6, for a condensed matter system, elements of Â are zero, corresponding to the translational and rotational degrees of freedom. Because the number of particles and their masses are the same in states 1 and 2, the translational contribution to the entropy is the same in both states and thus does not contribute to TΔS. On the contrary, the distribution of particles in these states might be different, and there can be a nonzero contribution of rotational degrees of freedom to TΔS. Here we discard this second term for two reasons: first, because we want to test the hypothesis that the TΔS computed in our simulation is of vibrational origin, and second, because in our simulations the HF(H2O)7 does not rotate and thus the rotational entropy is zero. To proceed with our derivation, we compute explicitly Q. Because the Jacobian of the r → q transformation is 1, the vibrational part of the original configurational integral can be written as 3N − 6

Q=

∏ j=1

3N − 6

∫ dqj exp( − βÂjj /2qj2 = ∏ (2πβ −1Âjj−1) j=1

(4)

Replacing eq 4 in eq 3 we get eq 2. Normally, computer codes use “mass weighted coordinates”, q = 1/2Û TM̂ r, which allow us to compare computational frequencies with experimental spectra. Because M̂ is not unitary, the Jacobian of this transformation is not 1. However, this Jacobian is the same in both states, and thus, a relation analogous to eq 2, in which the elements of Â are replaced by the elements of ω̂ = Û ′TM̂ −1K̂ Û ′, holds.

■

APPENDIX B. VIBRATIONAL ENTROPY OF A CLASSICAL SYSTEM We report here the derivation of the vibrational entropy to the free energy for a system of interacting classical harmonic oscillators (eq 2). Despite the fact that harmonic systems are widely used as examples because some exact results can be obtained for them, we were unable to find in textbooks and in the literature such a formula, apart for the trivial case of a single monodimensional case. Let us consider a system composed of N atoms with an associated Hamiltonian H(p,r) = 1/2pTM̂ −1p + 1/2rTK̂ r, where p and r are the 3N-dimensional vectors of the momenta and positions of the particles, M̂ is the (diagonal) matrix of the atomic masses, and K̂ is the (real symmetric) matrix of the harmonic constants. We assume that this system can exist in two states, 1 and 2, characterized by K̂ 1 and K̂ 2. These states may correspond, for example, to the initial and final state of the reaction discussed in this work. The starting point for the calculation of the entropic contribution to the free energy is the usual relation between

■

APPENDIX C. ALGORITHM FOR COMPUTING THE TRANSMISSION COEFFICIENT WHEN THE TRANSITION ENSEMBLE IS SAMPLED BY RESTRAINED MD/MC The calculation of the rate constant within the Vanden-Eijnden and Tal formalism74 requires the calculation of the expectation value of a suitable operator over the transition ensemble. In particular, the value of the operator is zero if a trajectory started from a point belonging to the transition state surface comes back to this surface before reaching the product domain. The estimation of the expectation value of such an observable by restrained MC/MD poses a problem. Configurations sampled by restrained MC/MD not always lie on the transition state surface. Consider, for example, the point r1 shown in Figure 13, which is in the reactant domain of the configuration space. A trajectory started from this point and ending up in the product crosses the transition state surface only once (green curve starting from r1 in Fig. 13), and give a zero contribution to the 13047



Article

No. 259 SIMBEDD−Advanced Computational Methods for Biophysics, Drug Design and Energy Research. Finally, the authors thank Ireland’s High-Performance Computing Centre (ICHEC) and the European DECIDEISA Project for the provision of computational resources.

■

(1) Vaida, V. Perspective: Water Cluster Mediated Atmospheric Chemistry. J. Chem. Phys. 2011, 135, 020901. (2) Rousseau, R.; Kleinschmidt, V.; Schmitt, U. W.; Marx, D. Assigning Protonation Patterns in Water Networks in Bacteriorhodopsin Based on Computed IR Spectra. Angew. Chem., Int. Ed. Engl. 2004, 43, 4804−4807. (3) Garczarek, F.; Gerwert, K. Functional Waters in Intraprotein Proton Transfer Monitored by FTIR Difference Spectroscopy. Nature 2005, 439, 109−112. (4) Tuckerman, M.; Laasonen, K.; Sprik, M.; Parrinello, M. Ab Initio Molecular Dynamics Simulation of the Solvation and Transport of H3O+ and OH− Ions in Water. J. Phys. Chem. 1995, 99, 5749−5752. (5) Geissler, P. L.; Dellago, C.; Chandler, D.; Hutter, J.; Parrinello, M. Autoionization in Liquid Water. Science 2001, 291, 2121−2124. (6) Mohammed, O. F.; Pines, D.; Dreyer, J.; Pines, E.; Nibbering, E. T. J. Sequential Proton Transfer Through Water Bridges in Acid-Base Reactions. Science 2005, 310, 83−86. (7) Woutersen, S.; Bakker, H. J. Ultrafast Vibrational and Structural Dynamics of the Proton in Liquid Water. Phys. Rev. Lett. 2006, 96, 138305. (8) Lee, J.-G.; Asciutto, E.; Babin, V.; Sagui, C.; Darden, T.; Roland, C. Deprotonation of Solvated Formic Acid: Car-Parrinello and Metadynamics Simulations. J. Phys. Chem. B 2006, 110, 2325−2331. (9) Iftimie, R.; Thomas, V.; Plessis, S.; Marchand, P.; Ayotte, P. Spectral Signatures and Molecular Origin of Acid Dissociation Intermediates. J. Am. Chem. Soc. 2008, 130, 5901−5907. (10) Siwick, B. J.; Cox, M. J.; Bakker, H. J. Long-Range Proton Transfer in Aqueous Acid-Base Reactions. J. Phys. Chem. B 2008, 112, 378−389. (11) Thomas, V.; Iftimie, R. Toward Understanding the Dissociation of Weak Acids in Water: 1. Using IR Spectroscopy to Identify ProtonShared Hydrogen-Bonded Ion-Pair Intermediates. J. Phys. Chem. B 2009, 113, 4152−4160. (12) Thomas, V.; Maurer, P.; Iftimie, R. On the Formation of Proton-Shared and Contact Ion Pair Forms During the Dissociation of Moderately Strong Acids: an Ab Initio Molecular Dynamics Investigation. J. Phys. Chem. B 2010, 114, 8147−8155. (13) Leopold, K. R. Hydrated Acid Clusters. Annu. Rev. Phys. Chem. 2011, 62, 327−349. (14) Kumar, P. P.; Kalinichev, A. G.; Kirkpatrick, R. J. Dissociation of Carbonic Acid: Gas Phase Energetics and Mechanism from Ab Initio Metadynamics Simulations. J. Chem. Phys. 2007, 126, 204315. (15) Gutberlet, A.; Schwaab, G.; Birer, O.; Masia, M.; Kaczmarek, A.; Forbert, H.; Havenith, M.; Marx, D. Aggregation-Induced Dissociation of HCl(H2O)4 Below 1 K: the Smallest Droplet of Acid. Science 2009, 324, 1545−1548. (16) Forbert, H.; Masia, M.; Kaczmarek-Kedziera, A.; Nair, N. N.; Marx, D. Aggregation-Induced Chemical Reactions: Acid Dissociation in Growing Water Clusters. J. Am. Chem. Soc. 2011, 133, 4062−4072. (17) Re, S.; Osamura, Y.; Morokuma, K. Coexistence of Neutral and Ion-Pair Clusters of Hydrated Sulfuric Acid H2SO4(H2O)n (n= 1−5): a Molecular Orbital Study. J. Phys. Chem. A 1999, 103, 3535−3547. (18) Re, S. Enhanced Stability of Non-Proton-Transferred Clusters of Hydrated Hydrogen Fluoride HF(H2O)n (n= 1−7): A Molecular Orbital Study. J. Phys. Chem. A 2001, 105, 9725−9735. (19) Weber, K. H.; Tao, F.-M. Ionic Dissociation of Perchloric Acid in Microsolvated Clusters. J. Phys. Chem. A 2001, 105, 1208−1213. (20) Cabaleiro-Lago, E. M.; Hermida-Ramón, J. M.; RodríguezOtero, J. Computational Study of the Dissociation of H-X acids (X=F, Cl, Br, I) in water clusters. J. Chem. Phys. 2002, 117, 3160−3168.

Figure 13. Sketch of reactive trajectories (green) that gives nonzero contribution to the estimate of the kRP and nonreactive (red) trajectories according to our revised definition of tp/z i *. The line denoted by ϵ(r) = z* represents the transition state surface.

hopping frequency. This is, of course, an artifact of the restrained MC/MD sampling as this trajectory should give a nonzero contribution to the estimation of the reactive flux. To cope with this problem, we changed the original p/z definition of tp/z i * as follows: ti * ≥ 0 is the time at which the trajectory started at ri and pi reaches the product or the transition state surface if it is started from a configuration in the product domain, or is the time at which the trajectory reaches the product or the transition state surface for a second time if it is started from a configuration in the reactant domain. According to this definition, the trajectories represented by the green lines in Figure 13 will give a nonzero contribution to the hopping frequency. On the contrary, the trajectories represented by the red lines will give a zero contribution to κs.

■

ASSOCIATED CONTENT

S Supporting Information *

Validation of the HCTH120 exchange and correlation functional. This material is available free of charge via the Internet at http://pubs.acs.org/.

■

REFERENCES

AUTHOR INFORMATION

Corresponding Author

*S. Meloni: e-mail, simone.meloni@epfl.ch. Present Address ∥

Laboratory of Computational Chemistry and Biochemistry, École Polytechnique Fédérale de Lausanne, CH-1015 Lausanne, Switzerland. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors thank Mantthias Krack for providing psuedopotentials. S.M. acknowledges financial support from the European Community under the Marie Curie Intra-European Fellowship for Career Development Grant No. 255406. The authors acknowledge financial support from SFI Grant No. 08-IN.1I1869. G.C. and S.M. acknowledge financial support from the Istituto Italiano di Tecnologia under the SEED project grant 13048



Article

(43) Romero, A. H.; Silvestrelli, P. L.; Parrinello, M. Compton Scattering and the Character of the Hydrogen Bond in Ice Ih. J. Chem. Phys. 2001, 115, 115−123. (44) Ludwig, R. Water: From Clusters to the Bulk. Angew. Chem., Int. Ed. Engl. 2001, 40, 1808−1827. (45) Fernández-Serra, M.; Artacho, E. Electrons and Hydrogen-Bond Connectivity in Liquid Water. Phys. Rev. Lett. 2006, 96, 016404. (46) -. CP2K: GPL State-Of-The-Art Methods for Ef f icient and Accurate Atomistic Simulations 2000−2012, http://www.cp2k.org/. (47) This modified version of CP2K is available in the SVN repository of the package. (48) Goedecker, S.; Teter, M.; Hutter, J. Separable Dual-Space Gaussian Pseudopotentials. Phys. Rev. B 1996, 54, 1703−1710. (49) Hartwigsen, C.; Goedecker, S.; Hutter, J. Relativistic Separable Dual-Space Gaussian Pseudopotentials from H to Rn. Phys. Rev. B 1998, 58, 3641−3662. (50) Krack, M. Pseudopotentials for H to Kr Optimized for GradientCorrected Exchange-Correlation Functionals. Theor. Chem. Acc. 2005, 114, 145−152. (51) Lippert, G.; Hutter, J.; Parrinello, M. A hybrid Gaussian and Plane Wave Density Functional Scheme. Mol. Phys. 1997, 92, 477− 488. (52) VandeVondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J. Quickstep: Fast and Accurate Density Functional Calculations Using a Mixed Gaussian and Plane Waves Approach. Comput. Phys. Commun. 2005, 167, 103−128. (53) VandeVondele, J.; Hutter, J. An Efficient Orbital Transformation Method for Electronic Structure Calculations. J. Chem. Phys. 2003, 118, 4365−4369. (54) VandeVondele, J.; Hutter, J. Gaussian Basis Sets for Accurate Calculations on Molecular Systems in Gas and Condensed Phases. J. Chem. Phys. 2007, 127, 114105. (55) The m-TZV2P basis set is composed of spherical harmonics of s, p, and d symmetry multiplied by radial functions which, in turn, are linear combination of Gaussian functions. The notation (3s3p2d/s2p) indicates that three functions of symmetry s and p, and two of symmetry d are used for non-hydrogen atoms, and one function of symmetry s and two of symmetry p are used for hydrogen. (56) Hamprecht, F. A.; Cohen, A. J.; Tozer, D. J.; Handy, N. C. Development and Assessment of New Exchange-Correlation Functionals. J. Chem. Phys. 1998, 109, 6264−6271. (57) Boese, A. D.; Doltsinis, N. L.; Handy, N. C.; Sprik, M. New Generalized Gradient Approximation Functionals. J. Chem. Phys. 2000, 112, 1670−1678. (58) Izvekov, S.; Voth, G. A. Ab Initio Molecular-Dynamics Simulation of Aqueous Proton Solvation and Transport Revisited. J. Chem. Phys. 2005, 123, 044505. (59) Izvekov, S.; Voth, G. A. Erratum: “Ab Initio MolecularDynamics Simulation of Aqueous Proton Solvation and Transport Revisited” [J. Chem. Phys. 123, 044505 (2005)]. J. Chem. Phys. 2006, 124, 039901. (60) Becke, A. D. Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior. Phys. Rev. A 1988, 38, 3098−3100. (61) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti Correlation-Energy Formula Into a Functional of the Electron Density. Phys. Rev. B 1988, 37, 785−789. (62) Kim, K.; Jordan, K. D. Comparison of Density Functional and MP2 Calculations on the Water Monomer and Dimer. J. Phys. Chem. 1994, 98, 10089−10094. (63) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab Initio Calculation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Force Fields. J. Phys. Chem. 1994, 98, 11623−11627. (64) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (65) Sinanolu, O. Many-Electron Theory of Atoms and Molecules. I. Shells, Electron Pairs vs Many-Electron Correlations. J. Chem. Phys. 1962, 36, 706.

(21) McCurdy, P. R.; Hess, W. P.; Xantheas, S. S. Nitric Acid-Water Complexes: Theoretical Calculations and Comparison to Experiment. J. Phys. Chem. A 2002, 106, 7628−7635. (22) Escribano, R.; Couceiro, M.; Gómez, P. C.; Carrasco, E.; Moreno, M. A.; Herrero, V. J. The Nitric Acid Hydrates: Ab Initio Molecular Study, and RAIR Spectra of the Solids. J. Phys. Chem. A 2003, 107, 651−661. (23) Odde, S.; Mhin, B. J.; Lee, S.; Lee, H. M.; Kim, K. S. Dissociation chemistry of Hydrogen Halides in Water. J. Chem. Phys. 2004, 120, 9524−9535. (24) Scott, J. R.; Wright, J. B. Computational Investigation of the Solvation of Nitric Acid: Formation of the NO3- and H3O+ Ion Pair. J. Phys. Chem. A 2004, 108, 10578−10585. (25) Li, S.; Tao, F.-M.; Gu, R. Theoretical Study on the Ionic Dissociation of Halosulfonic Acids in Small Water Clusters. Chem. Phys. Lett. 2006, 426, 1−7. (26) Odde, S.; Mhin, B. J.; Lee, K. H.; Lee, H. M.; Tarakeshwar, P.; Kim, K. S. Hydration and Dissociation of Hydrogen Fluoric Acid (HF). J. Phys. Chem. A 2006, 110, 7918−7924. (27) Li, S.; Qian, W.; Tao, F.-M. Ionic Dissociation of Methanesulfonic Acid in Small Water Clusters. Chem. Phys. Lett. 2007, 438, 190−195. (28) Ndongmouo, U. F. T.; Lee, M. S.; Rousseau, R.; Baletto, F.; Scandolo, S. Finite-Temperature Effects on the Stability and Infrared Spectra of HCl(H2O)6 Clusters. J. Phys. Chem. A 2007, 111, 12810− 12815. (29) Xie, Z.-Z.; Ong, Y.-S.; Kuo, J.-L. On the Effects of Basis-Set in Studying the Hydration and Dissociation of HF in Cubic HF(H2O)7 Clusters. Chem. Phys. Lett. 2008, 453, 13−17. (30) Ciccotti, G.; Meloni, S. Temperature Accelerated Monte Carlo (TAMC): a Method for Sampling the Free Energy Surface of NonAnalytical Collective Variables. Phys. Chem. Chem. Phys. 2011, 13, 5952−5959. (31) Ayotte, P.; Hébert, M.; Marchand, P. Why is Hydrofluoric Acid a Weak Acid? J. Chem. Phys. 2005, 123, 184501. (32) Park, J. M.; Laio, A.; Iannuzzi, M.; Parrinello, M. Dissociation Mechanism of Acetic Acid in Water. J. Am. Chem. Soc. 2006, 128, 11318−11319. (33) Hassanali, A.; Prakash, M. K.; Eshet, H.; Parrinello, M. On the Recombination of Hydronium and Hydroxide Ions in Water. Proc. Natl. Acad. Sci. U. S. A. 2011, 108, 20410−20415. (34) Kirkwood, J. C. Statistical Mechanics of Fluid Mixtures. J. Chem. Phys. 1935, 3, 300−313. (35) Bonella, S.; Vuilleumier, R. Path Integral Single Sweep Method for Quantum Free Energy Reconstruction. Submitted to J. Chem. Theory Comput. (36) Bonella, S.; Poma, A.; Monteferrante, M.; Ciccotti, G. Quantum Free Energy Barrier for Hydrogen Vacancy Diffusion in Na3AlH6. Phys. Chem. Chem. Phys. 2012, 14, 15458−15463. (37) Tuckerman, M. E.; Chandra, A.; Marx, D. Structure and Dynamics of OH−(aq). Acc. Chem. Res. 2006, 39, 151−158. (38) Duane, S. Hybrid Monte Carlo. Phys. Lett. B 1987, 195, 216− 222. (39) Mehlig, B.; Heermann, D.; Forrest, B. Hybrid Monte Carlo Method for Condensed-Matter Systems. Phys. Rev. B 1992, 45, 679− 685. (40) Marzari, N.; Vanderbilt, D. Maximally Localized Generalized Wannier Functions for Composite Energy Bands. Phys. Rev. B 1997, 56, 12847−12865. (41) Berghold, G.; Mundy, C.; Romero, A.; Hutter, J.; Parrinello, M. General and Efficient Algorithms for Obtaining Maximally Localized Wannier Functions. Phys. Rev. B 2000, 61, 10040−10048. (42) MLWFs are defined as wχ(x) = ∑icχ,iϕi(x), where wχ(x) denotes the χth MLWF, ϕi(x) is the ith KS orbital, and the coefficients cχ,i are those minimizing the width of the MLWFs: {c χ,i} χ,i = arg min{cχ,i}χ,i∑x⟨wχx2wχ⟩ − ⟨wχxwχ⟩2 (see ref 41 and references cited therein). 13049



Article

(66) Č ižek, J. On the Correlation Problem in Atomic and Molecular Systems. Calculation of Wavefunction Components in Ursell-Type Expansion Using Quantum-Field Theoretical Methods. J. Chem. Phys. 1966, 45, 4256. (67) Sinano lu, O.; Brueckner, K. A. Three Approaches to Electron Correlation in Atoms: A Chemistry-physics Interface; Yale University Press: New Haven, CT, 1970. (68) Martyna, G. J.; Tuckerman, M. E. A Reciprocal Space Based Method for Treating Long Range Interactions in Ab Initio and ForceField-Based Calculations in Clusters. J. Chem. Phys. 1999, 110, 2810− 2821. (69) Marcus, Y. Ion Solvation; Wiley: Chichester, U.K., 1985. (70) In our argument we assume that the comparison between experimental and computational pKa is fair. Of course, both can be affected by errors that might make the comparison difficult. In particular, theoretical results can be affected by larger systematic errors associated with the “force field”. The typical error in the computational estimate of pKa is one unit (see, for example, ref 82). Because this error is smaller than the difference between bulk and cluster values, we conclude that the observation that HF is a stronger acid in HF(H2O)7 than in bulk water is a genuine result. (71) Murtagh, F. COMPSTAT Lectures 4; Physica-Verlag: Würzburg, GE, 1985. (72) Eigen, M. Proton Transfer, Acid-Base Catalysis, and Enzymatic Hydrolysis. Part I: Elementary Processes. Angew. Chem., Int. Ed. Engl. 1964, 3, 1−19. (73) Xie, Z.; Bau, R.; Reed, C. A. A Crystalline H9O+4 Hydronium Ion Salt with a Weakly Coordinating Anion. Inorg. Chem. 1995, 34, 5403− 5405. (74) Vanden-Eijnden, E.; Tal, F. A. Transition State Theory: Variational Formulation, Dynamical Corrections, and Error Estimates. J. Chem. Phys. 2005, 123, 184103. (75) Anderson, J. B. Predicting Rare Events in Molecular Dynamics. Adv. Chem. Phys. 1995, 91, 381−431. (76) van Erp, T. S. Solvent Effects on Chemistry with Alcohols. An Ab Initio Study. Ph.D. thesis, Universiteit van Amsterdam, 2003. (77) van Erp, T. S.; Bolhuis, P. G. Elaborating Transition Interface Sampling Methods. J. Comput. Phys. 2005, 205, 157−181. (78) Monteferrante, M.; Bonella, S.; Ciccotti, G. Short Range Hydrogen Diffusion in Na3AlH6. Phys. Chem. Chem. Phys. 2011, 13, 10546−10555. (79) Eigen, M.; Kruse, W.; Maass, G.; De Maeyer, L. Rate Constants of Protolytic Reactions in Aqueous Solution. Prog. React. Kinet. 1964, 2, 285−318. (80) Joutsuka, T.; Ando, K. Hydration Structure in Dilute Hydrofluoric Acid. J. Phys. Chem. A 2011, 115, 671−677. (81) Buch, V.; Milet, A.; Vácha, R.; Jungwirth, P.; Devlin, J. P. Water Surface is Acidic. Proc. Natl. Acad. Sci. U. S. A. 2007, 104, 7342−7347. (82) Sprik, M. Computation of the pK of Liquid Water Using Coordination Constraints. Chem. Phys. 2000, 258, 139−150.

13050


Equilibrium and Rate Constants, and Reaction Mechanism of the HF

Recommend Documents