Quantum Chemistry in Solution by Combining 3D Integral Equation

J. Phys. Chem. B 2008, 112, 4337-4343

4337

Quantum Chemistry in Solution by Combining 3D Integral Equation Theory with a Cluster Embedding Approach Thomas Kloss, Jochen Heil, and Stefan M. Kast* Eduard-Zintl-Institut fu¨r Anorganische und Physikalische Chemie, Technische UniVersita¨t Darmstadt, Petersenstraβe 20, 64287 Darmstadt, Germany ReceiVed: NoVember 7, 2007; In Final Form: January 24, 2008

Free energy changes associated with chemical reactions in solution are treated by integral equation theory in the form of the 3D reference interaction site model (RISM) in combination with quantum-chemical calculations via an embedded cluster approach (EC-RISM). The electronic structure of the solute is computed self-consistently with the solvent structure by mapping the charge distribution of the solvent onto a set of discrete background point charges that are added to the molecular Hamiltonian. The EC-RISM procedure yields chemical accuracy in free energy predictions for several benchmark systems without adjusting empirical parameters. We apply the method to the standard reaction free energy for the gauche-trans equilibrium of 1,2-dichloroethane in water and to pKa shift calculations for trifluoroacetic acid/acetic acid and 4-nitroaniline/ aniline in water.

Introduction The adequate description of solvation phenomena is the key to successful theoretical and computational modeling of molecular interactions and chemical reactions in the liquid or solution phase.1 A proper account for solvent effects on the electronic structure of a solute is particularly important. The available methods range from a more or less first-principle, fully quantum-mechanical treatment such as the Car-Parrinello molecular dynamics2,3 (CPMD) simulation technique over mixed quantum-mechanical/molecular-mechanical models (QM/MM, see, e.g., refs 4 and 5) to continuum electrostatics approaches in various flavors. Representative for the latter so-called selfconsistent reaction field (SCRF) techniques are for instance the polarizable continuum model6,7 (PCM) and the conductor-like screening model8 (COSMO). The choice of a particular model, i.e., level of detail of the solvent description, hinges upon consideration of the balance between computational speed and accuracy. In situations where the solvent granularity becomes important, for instance when H-bonds between solute and solvent or other strongly directional properties are relevant, statistical theories of solvation such as the reference interaction site model (RISM) integral equations are the models of choice. These can be combined with selfconsistent field (SCF) quantum-chemical calculations based on a joint variational principle for the quantum-mechanical energy and the solvation free energy. Both 1D site-site theories9-11 as well as 3D grid-based theories are available.12-15 As an alternative to RISM-SCF, we demonstrate here that a simple cluster-embedding strategy applied to quantum-chemical calculations can be combined with the 3D RISM theory in such a way that predictions for solution-phase reaction thermodynamics with chemical accuracy is possible. Cluster embedding strategies are an instance of the general family of QM/MM methods that are based on coupling the electronic structure of * Corresponding author. Telephone: +49 (6151) 165397. Fax: +49 (6151) 164298. E-mail: [email protected].

a species, the QM part, to the environment that is described by a classical-mechanical (MM) model. We adapt the cluster embedding strategy that we used earlier for optimizing site charges in a solid-state environment16 in order to make it applicable to a liquid-state medium. The approach will be termed the embedded cluster reference interaction site model (ECRISM) and is characterized by a self-consistent procedure for the determination of the electronic structure of the solute and the embedding solvent. The solvent distribution functions are computed via 3D RISM integral equation theory, the resulting electrostatic potential distribution is mapped onto background charges that act on the solute electronic structure which in turn influences the solvent response. Such an approach has many advantages. The EC-RISM framework is easily coupled to an existing quantum chemistry code that allows the specification of discrete background charges. The iterative cycles can be controlled entirely by scripting languages, do not require modifications to the quantumchemical program, and converge quite quickly. The methodology works seamlessly with variational and non-variational theories, both for the integral equation and for the quantum chemistry part. As in the RISM-SCF approach, the solvent granularity and H-bonding are adequately accounted for. The framework can directly be applied to heterogeneous systems such as a reacting, solvent-exposed species in the active site of an enzyme. The EC-RISM approach represents a substantial enhancement of earlier attempts to couple electronic structure theory with integral equation methods via external fields. Gao and coworkers17 used 1D RISM theory in conjunction with semiempirical quantum chemistry and adjustable parameters in order to compute hydration free energies and the reaction free energy of the gauche-trans equilibrium of 1,2-dichloroethane in various solvents. Du and Wei18 described a cluster embedding approach in combination with 3D RISM theory but did not develop it into a tool for free energy predictions. In the next sections, we will describe the theoretical and computational background of the EC-RISM procedure and apply

10.1021/jp710680m CCC: $40.75 © 2008 American Chemical Society Published on Web 03/15/2008

4338 J. Phys. Chem. B, Vol. 112, No. 14, 2008

Kloss et al.

Figure 2. Thermodynamic cycle for a solution-phase reaction involving reactant A and product B.

thermodynamic cycle shown in Figure 2. ∆G can consequently be expressed as

∆G ) -∆Gsolv(A) + ∆Gvac + ∆Gsolv(B) Figure 1. Flowchart for self-consistent EC-RISM computations. Starting with the vacuum wave function, Ψvac, of the solute molecule, the set of solute partial charges qu is obtained by fitting with respect to the electrostatic potential, and the solvent site distribution functions g(r) around the solute with the new charges is calculated using 3D RISM integral equation theory. The embedding charge background is constructed based on g(r) and the solute wave function is recalculated in the presence of the background field. The cycle is repeated until convergence of the wave function Ψtot of the molecule in the presence of the background charges is reached, yielding an approximation to the wave function for the solvated system, Ψsolv.

it to a number of benchmark systems that have also been the subject of RISM-SCF studies,19-24 namely a conformational equilibrium and some important acid-base pairs. The first example is the gauche-trans equilibrium of 1,2-dichloroethane in water for which a wealth of experimental and theoretical reference data is available.17,19,25-30 We will test several levels of theory with respect to their performance for predicting the standard reaction free energy in water. The second example concerns thermodynamic calculations for acid-base equilibria that are a challenge for any quantum-mechanical solvation model due to the change of charge over the course of the reaction. Besides a pKa shift calculation for the aniline/4-nitroaniline pair, we intentionally chose the pathological example of trifluoroacetic acid/acetic acid for which even a high-level PCM treatment fails to give correct predictions.31 Theory Cluster Embedding Strategy. Similar to typical solid state applications, the solvent effect on the electronic structure of a solute (“cluster”) is approximated by the field originating from a discrete set of embedding point charges. Since the background charge distribution is perturbed by the solute wave function and vice versa, both must be solved in a self-consistent way, and Figure 1 illustrates the iteration cycles. The key component of the strategy is to map the continuous solvent site distribution onto the discrete background charge set as will be discussed below. While we could have used for instance molecular dynamics simulations for computing the solvent distribution functions, the use of approximate 3D RISM theory has many advantages such as a much less computational effort and, for several levels of 3D RISM theory, the availability of closedform expressions for the chemical potential of the solvated species. The latter are components of the total free energy change for chemical reactions in solution which are discussed next. Reaction Free Energy. For a general reaction Asolv f Bsolv involving two solvated molecules, each in a single conformation, we determine the Gibbs energy of reaction, ∆G, from the

(1)

where ∆Gvac ) ∆Hvac - T∆Svac is the vacuum free energy difference that is related to the vacuum enthalpy difference, ∆Hvac and entropy difference, ∆Svac, at temperature T, and ∆Gsolv(i) is the Gibbs energy of solvation for species i. Under the assumptions that (i) the nuclear degrees of freedom can be described classically, (ii) the number of rotational degrees of freedom remains constant over the course of the reaction, (iii) the vibrational entropy contributions cancel (the differences can indeed be shown to be negligibly small for our examples, which corresponds with other studies19), and (iv) pressure-volume work changes can be neglected, we obtain

∆G ≈ -∆Ereorg(A) - µex(A) + Evac(B) - Evac(A) + ∆Ereorg(B) + µex(B) (2) where

∆Ereorg(i) ) Esolv(i) - Evac(i)

(3)

represents the quantum-mechanical reorganization energy upon solvating a molecule of species i. Inserting eq 3 into eq 2 shows that the vacuum energy terms cancel and we end up with a relation (similar to ref 19) that can be evaluated by quantumchemical and statistical-mechanical tools

∆G ≈ ∆Esolv + ∆µex

(4)

Here, ∆Esolv ) Esolv(B) - Esolv(A) represents the difference of the quantum-mechanical energies of A and B in solution (see below), and ∆µex is the respective excess chemical potential difference. µex corresponds to the free energy change associated with transferring the solute (kept rigid and with solution-phase wave function) from the gas phase into solution. This quantity can be evaluated by integral equation theory (see below). If the solvent used in such a computation corresponds to a pressure of 1 bar, then ∆G is equivalent to the standard reaction free energy ∆G0. Equation 4 is immediately applicable to the gauche-trans equilibrium of dichloroethane in water, t-C2H4Cl2(aq) f g-C2H4Cl2(aq), that is characterized by ∆G0(g,t), i.e., our first numerical example. EC-RISM computations can also be used for the thermodynamic characterization of acid-base equilibria by constructing a thermodynamic cycle for relative pKa values of different acid-base pairs, so-called pKa shifts, ∆pKa. Specializing to the water case, we have from HA(aq) f H+(aq) + A-(aq) with ∆G0(A) and HB(aq) f H+(aq) + B-(aq) with ∆G0(B) for each reaction

pKa(i) )

1 ∆G0(i) RT log(10)

(5)

Quantum Chemistry in Solution

J. Phys. Chem. B, Vol. 112, No. 14, 2008 4339

with the molar gas constant R. The shift ∆pKa ) pKa(A) pKa(B) corresponds to the reaction HA(aq) + B-(aq) f HB(aq) + A-(aq) with ∆∆G0(A,B) ) ∆G0(A) - ∆G0(B). In this formulation, the terms concerning the free energy of the proton cancel out. With the assumptions given above we have

∆∆G(A,B) ) Esolv(HB) + Esolv(A-) - Esolv(HA) Esolv(B-) + µex(HB) + µex(A-) - µex(HA) - µex(B-) (6) and finally

∆pKa(A,B) )

1 ∆∆G0(A,B) RT log(10)

(7)

Quantum-Mechanical Energy of the Solute. Within the ECRISM scheme, Esolv is accessible from a quantum-chemical calculation of the molecule in the background charge field. The total Hamiltonian, H ˆ tot, accounts for the interactions between the nuclei (n), the electrons (e), and the point charges (q), respectively. We write

H ˆ tot ) H ˆ1 + H ˆ2

(8)

ˆ ne + H ˆ ee + H ˆ nn H ˆ1 ) H

(9)

with

representing pure solute properties, and

H ˆ2 ) H ˆ nq + H ˆ eq + H ˆ qq

(10)

for the molecule-charge interactions. Esolv is related to H ˆ 1 by

Esolv ) 〈Ψtot|H ˆ 1|Ψtot〉

(11)

ˆ tot. where Ψ tot is the wave function associated with H In practice we have to subtract the contribution of H ˆ 2 from the total energy Etot of the Hamiltonian H ˆ tot according to

Esolv ) Etot - E2

(12)

E2 ) Eself + Eq

(13)

with

Here, Eself is the self-energy of the background charges arising from H ˆ qq, the computations of which is computationally demanding and scales as O(N2) with the number of background charge centers. It is therefore advantageous to suppress the evaluation of this term in quantum chemistry codes such as Gaussian 0332 that we have used throughout in this study. The remaining contribution Eq associated with H ˆ nq + H ˆ eq describes the electrostatic interaction of the solute molecule with the background charges and can be written as

Eq )

∫ dr Fq(r)φ(r)

(14)

with the quantum-mechanical electrostatic potential of the solute φ at spatial points r and the charge density Fq. This density can be written as

Fq(r) )

∑γ qγFγgγ(r)

(15)

where the summation runs over solvent sites with index γ, qγ is the solvent site charge, Fg the solvent site density, and gγ the pair distribution function between solute and solvent sites that

has to be consistent with the electronic structure of the solute after convergence of the iterative cycle. To this end, the continuous charge density is mapped onto the discrete background charges at points ri via

q(ri) ) Fq(ri)∆V

(16)

where ∆V ) ∆x∆y∆z represents the volume occupied by a background charge center with grid spacing ∆x, ∆y, and ∆z, taken to be uniform in each spatial direction throughout. This procedure differs substantially from that of ref 18. As long as the computation of the self-energy is avoided, the number of background charge sites can be taken as large as the dimension of the computational grid used in 3D RISM theory without significant performance loss. 3D RISM Integral Equation Theory. We proceed by making the approximation that the solvent site density can be described by 3D RISM theory which was developed by several authors.33-36 Total (h) and direct (c) correlation functions between solute molecule and solvent sites are connected via

Fγhγ(r) )

cγ′ χγγ′(r) ∑ γ′

(17)

where χ represents the solvent susceptibility that is obtained from 1D RISM37 or the dielectrically corrected DRISM theory,38,39 respectively, and the star means convolution. The pair distribution function g is related to the total correlation function by gγ ) hγ + 1. The integral equations must be combined with another so-called closure relation to get a closed set of equations. Several approximate closures exist that have been tested in this work. The most important one is the 3D generalization of the hypernetted chain (HNC) closure40,41

hγ(r) ) exp[-βuγ(r) + hγ(r) - cγ(r)] - 1

(18)

where uγ is the interaction potential between solvent site γ and the solute molecule, β ) (kBT)-1 where kB is the Boltzmann constant. The electrostatic part of uγ is determined selfconsistently during EC-RISM cycles. Although it is in principle possible to use the exact potential that is derived from the wave function directly, we opted for practical reasons to approximate the long-range electrostatic solute-solvent interactions by discrete partial solute site charges. These are obtained by fitting to the quantum-mechanical electrostatic potential (ESP) via standard procedures implemented into Gaussian 03.32 Notwithstanding the quality of distribution functions derived from HNC theory, it is frequently not possible to obtain a numerical solution to the set of equations due to the strongly nonlinear character. Kovalenko and Hirata42,43 (KH) therefore formulated a partially linearized closure

hγ(r) )

{

exp[dγ(r)] - 1 T dγ(r) e0 dγ(r) T dγ(r) > 0

dγ(r) ) - βuγ(r) + hγ(r) - cγ(r)

(19)

which is numerically much more stable. Furthermore, we tested the 3D repulsive bridge correction44 (RBC) in a generalized closure which is defined as

hγ(r) ) exp[ - βuγ(r) + hγ(r) - cγ(r) + bγ(r)] - 1

(20)

with the exponential of the bridge function given by

exp[bγ(r)] )

ωRγ exp[-βurep ∏ R (r)] R*γ

(21)

4340 J. Phys. Chem. B, Vol. 112, No. 14, 2008

Kloss et al.

Here, urep is the repulsive r-12 part of the Lennard-Jones potential and ωaγ are elements of the intramolecular correlation functions matrix of the rigid solvent molecules. For several closures, closed-form expressions for the excess chemical potentials exist, in the case of HNC41,45

µex HNC ) β-1

[

1

]

1

∑γ Fγ ∫ dr 2 hγ2(r) - cγ(r) - 2 hγ(r)cγ(r)

(22)

and for KH41

µex KH ) β-1

[

1

1

∑γ Fγ ∫ dr 2 hγ2(r)Θ(-hγ(r)) - cγ(r) - 2 hγ(r)cγ(r)

]

(23) where Θ is the Heaviside step function. These are derived by coupling parameter integration and exhibit well-known deficiencies with respect to the prediction of absolute numbers. The h2 term in HNC and KH expressions increases the values as compared to exact results from free energy simulations; results from 1D theory particularly suffer from such behavior due to an unphysical dependence on auxiliary solute sites.46,47 In 3D theory the solute is always represented as a single (anisotropic) site while the number of solvent sites is small and typically constant. Even though 3D RISM free energies are therefore less seriously affected than 1D data, this means that it is for instance not useful to construct such a model based on the absolute hydration free energy of the proton48 for the purpose of pKa predictions. In practice however, frequently only relative differences are important. We have shown in the past how HNC theory can be used for quantitative predictions of relative free energy data, see, e.g., refs 49-51. For RBC no unique closed-form expression is possible due to the dependence of the value for the chemical potential on the thermodynamic integration path.52,53 It is, however, useful to apply a perturbation expansion around the HNC solution as an approximation44 ex ≈ µex µRBC HNC +

β-1

(r) + 1] ∑γ Fγ ∫ dr[ exp(bγ(r)) - 1][hHNC γ

(24)

In practice, the RBC closure can be used during EC-RISM cycles, while the computation of the chemical potential requires an additional HNC solution with the solute charge set obtained after EC-RISM convergence. With any level of 3D theory, the correlation functions for charged solutes (as required in the pKa case) need to be properly renormalized in order to remove periodicity artifacts. For the aqueous solutions treated in this study the compensating background charge method developed by Kovalenko and Hirata41 is sufficient. Computational Details The solvent susceptibility χ was calculated with the dielectrically corrected DRISM/HNC38,39 equations for several water models at T ) 298.15 K, a bulk density of F ) 0.03334 Å-3, and with the dielectric constant set to  ) 78. The equations were solved iteratively with the modified direct inversion of iterative subspace (MDIIS) technique54 on a logarithmically spaced grid55,56 ranging from 0.0059 to 164.02 Å until the maximum residual of the direct correlation functions fell below

Figure 3. Typical convergence behavior of an EC-RISM computation. The dipole moment (solid line) and the total energy change ∆Etot (dashed, on a logarithmic scale) between two successive iterations are shown as a function of the number of steps for the g-C2H4Cl2 isomer on the HNC/MP2 level of theory with the 6-311+G(d,p) basis set. The first step depicted corresponds to the vacuum situation. In solution the dipole moment of the solute increases; stable results are obtained after 6-7 iteration cycles.

10-8. The data reported in the tables in the next section were computed on the basis of a variant of the SPC/E water model57 with the Lennard-Jones parameter σ of hydrogen set to 1.0 Å.58,59 We additionally experimented with similar variants of the TIP3P model.60 Solute-solvent interactions were described by the sum of Lennard-Jones (LJ) and Coulomb potential using standard Lorentz-Berthelot mixing rules for the former. The LJ parameters for dichloroethane were taken from ref 19, whereas for pKa shift predictions, the Amber force field61 parameters were used. All force field data are summarized as Supporting Information. The solute structures were obtained by geometry optimization in vacuum using Gaussian 0332 with different basis sets within Hartree-Fock (HF), second-order Møller-Plesset (MP2) perturbation, and B3LYP density functional theory (DFT) approximations. The geometries were subsequently kept rigid during the EC-RISM iteration cycles. Dichloroethane was optimized at the same level of theory (see next section) that has been used in EC-RISM cycles, whereas EC-RISM computations for all acid-base pairs were performed on the basis of structures resulting from MP2/6-311+G(d,p) optimizations. Two acid-base pairs were examined, 4-nitroaniline/aniline and trifluoroacetic acid/acetic acid. Only the global optima for the acid forms of the acetic acid derivatives were investigated. PCM calculations with the UA0 (united atom topological) model were furthermore performed for dichloroethane, using gas-phase structures optimized at the same level of theory that has been applied in solution. An additional set of PCM calculations has been done taking into account structures obtained from geometry optimization in PCM solution. For 3D RISM calculations, long-range electrostatic interactions were treated by Ewald summation taking into account a uniform background charge where necessary.41 The integral equations were solved on rectangular grids of 643 points with a spacing of ∆x ) ∆y ) ∆z ) 0.4 Å in the case of dichloroethane, and of 1403 points and 0.18 Å spacing in all other cases. We again used the MDIIS method54 up to a convergence threshold of 10-6 in the direct correlation functions. For some of the acid base pairs, a low-pass filter was applied in order to remove oscillations beyond a certain absolute value for reciprocal space vectors, kmax, such that convergent solutions are found reliably: such a treatment was necessary for CH3COOH (HNC, kmax ) 9 Å-1), CF3COOH (HNC, 7 Å-1), and all protonated anilines (HNC/KH, 10 Å-1). The energetic influence of this modification is negligible: for g-dichloroethane (optimized at HF/6-311+G(d,p)) µex changes from 16.9740 to 16.9746 kcal mol-1 upon application of the 8 Å-1 cutoff and to 16.9705 for 7 Å-1. We find for a Na+ ion (LJ parameters 

Quantum Chemistry in Solution

J. Phys. Chem. B, Vol. 112, No. 14, 2008 4341

TABLE 1: Free Energy Differences ∆G0(g,t) (298.15 K, in kcal mol-1)a for the Gauche-Trans Equilibrium of 1,2-C2H4Cl2 in Water from EC-RISM, RISM-SCF, and PCM Calculations theory level/basis set

HNCb

RBCb

KHb

RISM-SCFc

PCMd

HF/6-311G HF/6-311+G(d,p) MP2/6-311G MP2(FC)/6-311+G(d,p) MP2/6-311+G(d,p) MP2(FC)/aug-cc-pVTZ B3LYP/6-311+G(d,p)

0.18 0.27 0.29 (0.24) 0.13 (0.11) 0.12 (0.10) -0.05 0.09

0.39 0.48 0.42 (0.46) 0.29 (0.34) 0.28 (0.33) 0.19 0.30

0.20 0.28 0.30 (0.26) 0.14 (0.12) 0.13 (0.11) -0.04 0.11

-0.58 -0.74 -0.33

0.40 (0.16) 0.53 (0.33) 0.51 (0.32) 0.37 (0.22) 0.36 (0.21) 0.20 0.32 (0.14)

-0.79

a Experimental consensus results:28-30 -0.05 ( 0.05 kcal mol-1. b For MP2 EC-RISM calculations: data from using either the HF wave function for solute charges and potential, or from corresponding MP2 properties (in parentheses). c From ref 19. d PCM results: including electrostatic and nonelectrostatic contributions for gas-phase geometry, and with only electrostatic contributions from the structure optimized under PCM conditions at the respective level (in parentheses).

TABLE 2: Contributions to ∆G0(g,t) (298.15 K, in kcal mol-1)a for the Gauche-Trans Equilibrium of 1,2-C2H4Cl2 in Water from EC-RISM-HNC Calculations theory level/basis set

∆Etot - ∆Eself

∆Eq

∆µex

HF/6-311G HF/6-311+G(d,p) MP2/6-311G MP2(FC)/6-311+G(d,p) MP2/6-311+G(d,p) MP2(FC)/aug-cc-pVTZ B3LYP/6-311+G(d,p)

-2.77 -2.03 -2.93 (-2.13) -2.17 (-1.71) -2.16 (-1.70) -2.19 -2.08

-5.89 -4.55 -6.40 (-4.76) -4.60 (-3.65) -4.57 (-3.62) -4.29 -4.30

-2.94 -2.26 -3.18 (-2.29) -2.29 (-1.83) -2.29 (-1.82) -2.15 -2.13

a For MP2 EC-RISM calculations: data from using either the HF wave function for solute charges and potential, or from corresponding MP2 properties (in parentheses).

TABLE 3: Dipole Moment (in Debye) of gauche-1,2-C2H4Cl2 in Vacuuma and Water from EC-RISM and PCM Calculations theory level/basis set

HNCb

RBCb

KHb

PCMc

vacuumb

HF/6-311G HF/6-311+G(d,p) MP2/6-311G MP2(FC)/6-311+G(d,p) MP2/6-311+G(d,p) MP2(FC)/aug-cc-pVTZ B3LYP/6-311+G(d,p)

5.14 4.12 5.33 (4.41) 4.09 (3.54) 4.08 (3.53) 4.04 3.89

4.83 3.86 5.00 (4.14) 3.83 (3.31) 3.82 (3.29) 3.75 3.63

5.12 4.11 5.31 (4.39) 4.08 (3.53) 4.07 (3.51) 4.02 3.87

5.16 4.12 5.36 4.10 4.09 4.03 3.88

3.99 3.15 4.11 (3.39) 3.12 (2.69) 3.12 (2.68) 2.99 (2.67) 2.92

a Experimental result from refs 62 and 63: 3.12 D. b For MP2 calculations: data from using either the HF wave function for solute charges and potential, or from corresponding MP2 properties (in parentheses). c With structures optimized in the gas phase at the respective level.

) 0.0095 kcal mol-1, σ ) 2.85 Å) -85.3482 (no cutoff), -85.3476 (kmax ) 8 Å-1), and -85.3879 kcal mol-1 (7 Å-1), respectively. All 3D integrals occurring in energy evaluations were approximated by discrete sums over the grid. For the ECRISM iterations, we used the change in the total energy as convergence criterion, which should fall below 6.3 × 10-4 kcal mol-1 between two subsequent cycle steps. Typical trends for the convergence behavior are shown in Figure 3. Results and Discussion 1,2-Dichloroethane. Free energy results from using a number of different levels of theory and methods are summarized in Table 1, and Table 2 shows the contribution of various terms from converged EC-RISM-HNC calculations to the total value of ∆G0 in order to elucidate possible sources of error compensation. This number can be reconstructed according to ∆G0 ) ∆Etot - ∆Eself - ∆Eq + ∆µex and should be compared with the experimental consensus value of 0 ( 0.05 kcal mol-1.28-30 Basically all closure approximations and quantum-chemical levels of theory yield similar EC-RISM results within a range of ca. 0.4 kcal mol-1. In contrast to RISM-SCF results,19 there is a clear trend toward the experimental value with increasing level of theory. PCM theory also approaches this number, but the values are comparable to the worst-performing EC-RISM closure, the RBC. HNC and KH give consistently better results, with a slight advantage for HNC. Purely electrostatic PCM free

energies with geometries optimized in solution are somewhat better than full PCM free energies with gas-phase geometry but still do not reach the quality of EC-RISM-HNC predictions. It is important, however, not to put to much weight on the absolute numbers and not to overinterpret the results, given the small range of free energy values obtained by the various theoretical approaches. The influence of error compensation concerning the quantummechanical level is clearly visible by inspection of free energy contributions in Table 2 and dipole moments (experimental data taken from refs 62 and 63) in Table 3. As expected, dipole moments increase in aqueous solution relative to the gas phase. Adding diffuse functions is important for reliable MP2 calculations, whereas the effect of the frozen core approximation is negligible. HF theory with small basis set clearly benefits from error compensation if taken together with the solvation effect. DFT with the B3LYP model appears to be a useful alternative to MP2 calculations. The fact that chemical potential and electronic structure are mutually dependent adds robustness to the EC-RISM procedure, as can be seen by comparing MP2 results obtained by either taking the HF or the MP2 potential for describing the electrostatic interactions. At least for HNC and KH, the free energy predictions are slightly improved (more favorable for the gauche form) in the latter case, although the gauche dipole moments from MP2 wave functions are systematically smaller. This leads to less gauche-favorable chemical

4342 J. Phys. Chem. B, Vol. 112, No. 14, 2008

Kloss et al.

TABLE 4: pKa Shifts (298.15 K) for Different Acid-Base Pairs in Aqueous Solution ∆pKa

a

system

theory level/basis set

HNC

KH

exp.a

CF3COOH-CH3COOH CF3COOH-CH3COOH aniline-nitroaniline aniline-nitroaniline

B3LYP/6-311+G(d,p) MP2/6-311+G(d,p) B3LYP/6-311+G(d,p) MP2/6-311+G(d,p)

-4.01 -4.25 5.19 2.94

-5.54 -5.61 5.66 3.33

-4.53 -4.53 3.6 3.6

From refs 31 and 64.

TABLE 5: Contributions to ∆∆G0(A,B) (298.15 K, in kcal mol-1) for Different Acid-Base Pairs in Aqueous Solution from EC-RISM-HNC Calculations system

theory level/basis set

∆∆Etot - ∆∆Eself

∆∆Eq

∆∆µex

CF3COOH-CH3COOH CF3COOH-CH3COOH aniline-nitroaniline aniline-nitroaniline

B3LYP/6-311+G(d,p) MP2/6-311+G(d,p) B3LYP/6-311+G(d,p) MP2/6-311+G(d,p)

-23.04 -24.67 5.13 11.21

-55.68 -58.83 24.35 30.50

-27.17 -28.36 12.14 15.28

potentials that are counter-balanced by more gauche-favorable quantum-mechanical energies. Dipole moments from RBC calculations are, however, always smaller than from HNC or KH closures which leads to stronger deviation from the experimental free energy value. This indicates that solvent polarization is not as effective since the solvent distribution functions exhibit much lower peaks as compared to HNC results. Although we did not fit any empirical parameters in this study, it is worthwhile to examine the influence of the water model. The SPC/E or the TIP3P models with different Lennard-Jones parameters σH assigned to hydrogen yield basically identical results in the dichloroethane case. A value of σH ) 0.4 Å, which is frequently used in RISM-based studies, gives for MP2/6311+G(d,p)/EC-RISM-HNC with HF density for electrostatics the result ∆G0(298.15 K) ) 0.12 kcal mol-1 with the SPC/E and 0.13 kcal mol-1 with the TIP3P model; TIP3P with σH ) 1.0 Å yields 0.13 kcal mol-1. These data compare well with the SPC/E result (1.0 Å) of 0.12 kcal mol-1. Geometrical and force field differences among water models play apparently a negligible role here. pKa Shift Calculations. With the information gathered on the dichloroethane system, we applied the methodology separately to a neutral base (aniline/nitroaniline) and a neutral acid system (acetic/trifluoroacetic acid) for the prediction of pKa shifts. Particularly in the latter difficult case, PCM models even with a high level of quantum-chemical theory fail, deviating by 3 or more pKa units.31 Since EC-RISM-HNC and KH theories with the MP2 (and HF wave function for electrostatic interactions) or B3LYP/6-311+G(d,p) levels represents a good compromise between computational effort and accuracy for the dichloroethane gauche-trans equilibrium, we applied these approximations to the acid-base pairs. pKa results are summarized in Table 4 (experimental data taken from refs 31 and 64), associated free energy components are shown in Table 5. We find that the MP2 level used with the HNC closure gives very accurate results, while the prediction of the pKa shift for the acids works better than for the bases. KH is in general less reliable, which is not unexpected since the peak size in distribution functions for H-bonding site pairs is substantially influenced by the partial linearization of the closure. B3LYP results show a less clear tendency. Reliable DFT calculations would require a broader dataset of test cases in order to examine useful combinations of basis sets and exchange-correlation functionals. The fact that B3LYP yields a good prediction for the acid pair is apparently related to the consistency of free energy components shown in Table 5 among different levels of theory. This is not the case for the bases and might be a sign

for the general difficulty of pKa calculations for amines which has also been reported elsewhere.65 The choice of the water model plays a bigger role for pKa shift calculations as compared to the dichloroethane situation which is related to the stronger polar character of the species involved and the resulting more pronounced effect of H-bonding. For the MP2/HNC level of theory in the CF3COOH-CH3COOH case, we find with the TIP3P model and σH ) 1.0 Å ∆pKa ) -4.21, i.e., basically the same result as with SPC/E and the identical hydrogen parameter. In contrast, the analogous calculation with SPC/E and σH ) 0.4 Å yields ∆pKa ) -1.21 and with TIP3P (σH ) 0.4 Å) -1.20. This means that only hydrogen’s Lennard-Jones parameters are critically important and σH should be set to 1.0 Å for a quantitative description whereas other structural and force field details play an insignificant role. Concluding Remarks The intention of this study was to demonstrate the efficiency of the combination of 3D RISM integral equation theory with cluster-embedding electronic structure calculation in order to treat reaction thermodynamics in solution. We have confined ourselves to the presentation of the methodology and to the application to a few important benchmark cases. For these test scenarios, we find basically quantitative chemical accuracy without the need to adjust some empirical parameters. Instead, we have used well-established water models and force fields for the nonelectrostatic part of the Hamiltonian that describes the solute-solvent interaction. Depending on the polarity of the solute species, only the choice of Lennard-Jones parameters for hydrogen of the water model might be important, which should be sufficiently large. The best description of solutionphase thermodynamics is consistently possible by the HNC closure used together with post-HF quantum chemistry and a sufficiently diffuse basis set. Since increased levels of theory in both solvation and electronic structure terms lead to increased quality of the prediction, we can be confident that the EC-RISM methodology has a sound physical basis and can be systematically improved. Only the application to a much larger dataset of test compounds will enable us to assess and optimize the general usability of the EC-RISM approach. We are therefore presently studying such large sets in order to establish guidelines for suitable levels of theory, particularly for pKa shift calculations. Error cancellation will most likely play a role for the RISMbased chemical potential prediction which means that for instance pKa shifts between very dissimilar species will lead to

Quantum Chemistry in Solution increased error. Better free energy functionals as well as correction by fitting a linear correlation equation to experimental data along the lines of, e.g., refs 64-67, might allow the calculation of absolute quantities. Furthermore, we are investigating the limitations of point charge models for describing solute-solvent electrostatics. Although polarizability is in principle accounted for by the EC-RISM procedure, the choice of only a small number of atom-based charge centers might be too restrictive. All these aspects will be important for optimizing the EC-RISM methodology by balancing accuracy and computational effort. Acknowledgment. We thank the Deutsche Forschungsgemeinschaft, the Fonds der Chemischen Industrie, and the AdolfMesser-Stiftung for financial support. Supporting Information Available: Force field parameters for the species examined in this study. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Hirata, F., Ed.; Molecular Theory of SolVation; Kluwer Academic Publishers: Boston, MA, 2003. (2) Car, R.; Parrinello, M. Phys. ReV. Lett. 1985, 55, 2471. (3) Marx, D.; Hutter, J. Modern Methods and Algorithms of Quantum Chemistry; Grotendorst, J., Ed.; NIC, FZ Ju¨lich, 2000; p 301. (4) Dungsrikaew, V.; Limtrakul, J.; Hermansson, K.; Probst, M. J. Quant. Chem. 2004, 96, 17. (5) Bakowies, D.; Thiel, W. J. Phys. Chem. 1996, 100, 10580. (6) Cramer, C. J.; Truhlar, D. G. Chem. ReV. 1999, 99, 2161. (7) Tomasi, J.; Mennucci, B.; Cammi, R. J. Chem. ReV. 2005, 105, 2999. (8) Klamt, A. COSMO-RS: From Quantum Chemistry to Fluid Phase Thermodynamics and Drug Design; Elsevier Science Ltd.: Amsterdam, 2005. (9) Ten-no, S.; Hirata, F.; Kato, S. Chem. Phys. Lett. 1993, 214, 391. (10) Ten-no, S.; Hirata, F.; Kato, S. J. Chem. Phys. 1994, 100, 7443. (11) Sato, H.; Hirata, F.; Kato, S. J. Chem. Phys. 1996, 105, 1546. (12) Sato, H.; Kovalenko, A.; Hirata, F. J. Chem. Phys. 2000, 112, 9463. (13) Minezawa, N.; Kato, S. J. Chem. Phys. 2007, 126, 054511. (14) Yoshida, N.; Hirata, F. J. Comp. Chem. 2006, 27, 453. (15) Casanova, D.; Gusarov, S.; Kovalenko, A.; Ziegler, T. J. Chem. Theo. Comput. 2007, 3, 458. (16) Hauptmann, S.; Dufner, H.; Brickmann, J.; Kast, S. M.; Berry, R. S. Phys. Chem. Chem. Phys. 2003, 5, 635. (17) Shao, L.; Yu. H. A.; Gao, J. J. Phys. Chem. A 1998, 102, 10366. (18) Du, Q.; Wei, D. J. Phys. Chem. B 2003, 107, 13463. (19) Lee, J. Y.; Yoshida, N.; Hirata, F. J. Phys. Chem. B 2006, 110, 16018. (20) Kawata, M.; Ten-no, S.; Kato, S.; Hirata, F. Chem. Phys. Lett. 1995, 240, 199. (21) Kawata, M.; Ten-no, S.; Kato, S.; Hirata, F. J. Am. Chem. Soc. 1995, 117, 1638. (22) Sato, H.; Hirata, F. J. Am. Chem. Soc. 1999, 121, 3460. (23) Kawata, M.; Ten-no, S.; Kato, S.; Hirata, F. J. Phys. Chem. 1996, 100, 1111.

J. Phys. Chem. B, Vol. 112, No. 14, 2008 4343 (24) Kawata, M.; Ten-no, S.; Kato, S.; Hirata, F. Chem. Phys. 1996, 203, 53. (25) Mele´ndez-Paga´n, Y.; Taylor, B. E.; Ben-Amotz, D. J. Phys. Chem. B 2001, 105, 520. (26) Madurga, S.; Vilaseca, E. J. Phys. Chem. A 2004, 108, 8439. (27) Madurga, S.; Vilaseca, E. Chem. Phys. Lett. 2005, 406, 409. (28) Cappelli, C.; Corni, S.; Tomasi, J. J. Phys. Chem. A 2001, 105, 10807. (29) Wiberg, K. B.; Keith, T. A.; Frisch, M. J.; Murcko, M. J. Phys. Chem. 1995, 99, 9072. (30) Kato, M.; Abe, I.; Taniguchi, Y. J. Chem. Phys. 1999, 110, 11982. (31) Namazian, M.; Kalantary-Fotooh, F.; Noorbala, M. R.; Searles, D. J.; Coote, M. L. J. Mol. Struct. THEOCHEM 2006, 758, 275. (32) Frisch, M. J.; et al. Gaussian 03, revision C.02; Gaussian, Inc.: Wallingford, CT, 2004. (33) Ikeguchi, M.; Doi, J. J. Chem. Phys. 1995, 103, 5011. (34) Beglov, D.; Roux, B. J. Phys. Chem. B 1997, 101, 7821. (35) Du, Q.; Beglov, D.; Roux, B. J. Phys. Chem. B 2000, 104, 796. (36) Kovalenko, A.; Hirata, F. Chem. Phys. Lett. 1998, 290, 237. (37) Chandler, D.; Andersen, H. C. J. Chem. Phys. 1972, 57, 1930. (38) Perkyns, J.; Pettitt, B. M. Chem. Phys. Lett. 1992, 190, 626. (39) Perkyns, J.; Pettitt, B. M. J. Chem. Phys. 1992, 97, 7656. (40) Morita, T.; Hiroike, K. Prog. Theor. Phys. 1961, 25, 537. (41) Kovalenko, A.; Hirata, F. J. Chem. Phys. 2000, 112, 10391. (42) Kovalenko, A.; Hirata, F. J. Phys. Chem. B 1999, 103, 7942. (43) Kovalenko, A.; Hirata, F. J. Chem. Phys. 1999, 110, 10095. (44) Kovalenko, A.; Hirata, F. J. Chem. Phys. 2000, 113, 2793. (45) Singer, S., J.; Chandler, D. Mol. Phys. 1985, 55, 621. (46) Freedman, H.; Truong, T. N. J. Chem. Phys. 2004, 121, 2187. (47) Sato, K.; Chuman, H.; Ten-no, S. J. Phys. Chem. B 2005, 109, 17290. (48) Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem. B 2006, 110, 16066. (49) Schmidt, K. F.; Kast, S. M. J. Phys. Chem. B 2002, 106, 6289. (50) Kast, S. M. Chem. Phys. Chem. 2004, 5, 449. (51) Schilling, B.; Brickmann, J.; Kast, S. M. Phys. Chem. Chem. Phys. 2006, 8, 1086. (52) Kast, S. M. Phys. ReV. E 2003, 67, 041203. (53) Kast, S. M. Phys. ReV. E 2006, 73, 012201. (54) Kovalenko, A.; Ten-no, S.; Hirata, F. J. Comput. Chem. 1999, 20, 928. (55) Talman, J. D. J. Comput. Phys. 1978, 29, 35. (56) Rossky, P. J.; Friedman, H. L. J. Chem. Phys. 1980, 72, 5694. (57) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91 6269. (58) Maw, S.; Sato, H.; Ten-no, S.; Hirata, F. Chem. Phys. Lett. 1997, 276, 20. (59) Sato, H.; Hirata, F. J. Chem. Phys. 1999, 111, 8545. (60) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R., W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (61) Fox, T.; Kollman, P. A. J. Phys. Chem. B 1998, 102, 8070. (62) Taniguchi, Y.; Takaya, H.; Wong, P. T. T.; Whalley, E. J. Chem. Phys. 1981, 75, 4815. (63) Takaya, H.; Taniguchi, Y.; Wong, P. T. T.; Whalley, E. J. Chem. Phys. 1981, 75, 4823. (64) Eckert, F.; Klamt, A. J. Comput. Chem. 2006, 27, 11. (65) Klamt, A.; Eckert, F.; Diedenhofen, M.; Beck, M. E. J. Phys. Chem. A 2003, 107, 9380. (66) Klicic´, J. J.; Friesner, R. A.; Liu, S.-Y.; Guida, W. C. J. Phys. Chem. A 2002, 106, 1327. (67) Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem. A 2006, 110, 2493.