Calculations of chemical processes in solutions - ACS Publications

1640

The Journal of Physical Chemistry, Vol. 83, No. 12, 1979

quencies is comparable in the two cases and the force constants show that the general tendency of the added electrons is to strengthen single bonds (C-C) and weaken the multiple bonds (C=N, C=C). The spectroscopic evidence also strongly suggests that the TCNE and TCNQ dianion and trianion salts are insulators or, possibly, weak semiconductors, with one dramatic exception. No resistance studies have yet been made on the metal-rich Na3TCNE salts that show a low concentration of TCNE4-, but the remarkable low temperature enhancement of the a,-mode infrared intensities is accompanied by an intense broad far-infrared absorption that may signal onset of metallic behavior.16 The absence of a similar behavior for Na3TCNQ suggests some significant structural difference between the two trianion salts. The only possibility suggested by the available data is that the TCNQ trianion, judging from the large range in C r N stretching frequencies reported for Na3TCNQ (262 cm-l),8 may be much more distorted from the original parent molecule’s planar D 2 h structure. By contrast, the range of C=N stretching frequencies for TCNE3- is only 82 cm-l. Acknowledgment. This research has been sponsored by the National Science Foundation under Grant CHE-7709653.

Arieh Warshel

References and Notes W. D. Phillips, J. C. Rowell, and S. I. Weissman, J . Chem. fhys., 33, 626 (1960). (a) M. Itoh, J . Am. Chem. Soc., 92, 886 (1970); (b) Bull. Chem. SOC.Jpn., 45, 1947 (1972). J. Stanley, D. Smith, B. Latimer, and J. P. Devlin, J . fhys. Chem., 70, 2011 (1966). J. C.Moore, D. Smith, Y. Youhne, and J. P. Devlin, J . Phys. Chem., 75, 325 (1971). J. J. Hinkel and J. P. Devlin, J . Chem. fhys., 58, 4750 (1973). G. R. Anderson and R. L. McNeeb, Paper R5, Symposium on Molecular Structure and Spectroscopy, Ohio State University, June, 1963. F. A. Miller, 0. Sala, P. Devlin, J. Overend, E. Lippert, W. Luder, H. Moser, and J. Varchim, Spectrochim. Acta, 20, 1233 (1964). M. Khatkale and J. P. Devlin, J . Chem. fhys., 70, 1851 (1979). M. R. Suchanski and R. P. Van Duyne, J . Am. Chem. Soc., 98, 250 (1976). A. Bieber and J. J. Andre, Chem. Phys., 5, 166 (1974). G. R. Anderson and J. P. Devlin, J . fhys. Chem., 79, 1100 (1975). R. Bozio, A. Girlando, and C. Pecile, Chem. fhys., 21, 257 (1977). (a) R. Bozio and C. Pecile, J. Chem. fhys., 67, 3864 (1977); (b) M. J. Rice, N. 0. Lipari, and S. Strassler, fhys. Rev. Lett., 39, 1359 (1977). B. R. Penfold and W. N. Lipscomb, Acta. Crystalkgr., 14, 589 (1961). The electron enters the seventh MO of the conjugated T system which is either of b or b symmetry (J. J. Hinkle, Ph.D. Thesis, Oklahoma State Univeri&, 19?4). In either case, the orbital is antibonding CF3-4 and nonbonding otherwise. See, for example, (a) A. Brau, P. Bruesch, J. P. Farges, W. Hinz, and D. Kuse, fhys. Status Solidi, 62, 615 (1974); (b) A. A. Bright, A. F. Gariio, and A. J. Heeger, Sold State Commun., 13, 943 (1973).

Calculations of Chemical Processes in Solutions Arieh Warshel Department of Chemistry, University of Southern California, Los Angeles, California 90007 (Received September 29, 1978; Revised Manuscript Received February 27, 1979) Publicailon costs assisted by the National Institutes of Health

A new model for the calculation of solvation energies of polar solutes is presented. The model, referred to here as the “surface constrained soft sphere dipoles” (SCSSD) model, overcomes some of the problems in the theoretical treatment of solvent effects. On one hand, it avoids the problems of dielectric continuum approaches by explicitly including solvent molecules. On the other hand, it overcomes the extreme inefficiency of the supermolecule approaches by drastically simplifyingthe solvent molecules while retaining their main physical features. This is done by representing the solvent molecules as soft sphere point dipoles and minimizing the solute-solvent energy with respect to the orientations and positions of these dipoles while constraining the surface dipoles to have the positions of the bulk molecules. The model allows evaluation of the microscopic energy balance of charge separation in polar liquids. It is shown that the effective dielectric is more than 20, even when the charges are separated by 3.5 A. This microscopic dielectric effect is due to the fact that the solvent system has a large number of degrees of freedom and any change in the solute-solute electrostatic interactions can be balanced by a rearrangement of the solvent molecules. The SCSSD model makes it possible to perform quantum mechanical and empirical calculations for many chemical phenomena in solutions. This is demonstrated here by evaluating solvation enthalpies of ions with different charges and radii, dissociation enthalpies of acids in aqueous solution, parallel base stacking of the cytosine dimer, and electronic excitation energies. The good agreement between the calculated and experimental results indicates that the model can be used for quantitative studies of solvent effects in chemistry.

Introduction Since most quantum mechanical and empirical calculation schemes are designed to treat isolated molecules, they cannot be applied to many chemical and biochemical processes in solution. For example, the gas phase energy of ionic reactions can amount to 100-200 kcal/mol, while the corresponding energies in solution are only 10-30 kcal/mol. The difference between these energies is due

exclusively to the medium. Omitting the medium from the calculations amounts, in some cases, to dealing with only half of the energy contributions of the complete system. One possible way to overcome part of this problem is to treat the solvent as a dielectric continuum.’ Such an approach was used by Sinanoglu and co-workers,2Claverie et and more recently by Huron and C l a ~ e r i eBev,~

0022-3654/79/2083-1640$01.00/00 1979 American Chemical Society

Calculations of Chemical Processes in Solutions

eridge et al.,5aand Hylton et a1.6 Unfortunately, the dielectric models present several basic problems: (1)the use of a bulk dielectric constant is unjustified due to saturation effects in the first solvation shell; (2) the solvent cavity is treated as a sphere with no attempt to treat its realistic shape; (3) the solute charge distribution is represented (in most cases) by the dipole approximation; (4) the inclusion of the field of the solvent dipoles in the solute quantum mechanical calculations is complicated. Probably the most crucial problem of the continuum models is the determination of the size of the solute cavity. Mathematically the cavity starts where the dielectric constant changes from t = 1 to the bulk dielectric value of t = tO. Actually no such point exists, and the cavity radius is an arbitrary parameter. The continuum models might be useful in estimating conformational energies of polar molecules. However, they are of limited validity for the quantitative treatment of ionic reactions and for studies of the balance of forces between the solute-solute and solute-olvent interactions.’ An alternative approach is to include a limited number of solvent molecules in the calculations. In this case the system of the solute and a small number of solvent molecules is treated as a “supermolecule”. Such approaches have been used extensively by Pullman and Pullman and their co-workers.* Recent related studies have been reported by J ~ r g e n s e n .Combinations ~ of supermolecule and continuum approaches have been reported by Newton,loKollman and Kuntz,ll and Schnuelle and Beveridgeasb The supermolecule approach may be useful in determining specific solvation sites. However, for most processes in solution the overall solvation energy rather than the optimal solvation sites is of interest. In fact, as will be shown here, the solvation energy can be obtained only by a treatment which includes many solvent molecules. Furthermore, the calculated solvation energy might be meaningless if its evaluation did not allow solvent rearrangement. Optimization of the solvent coordinates is extremely expensive in terms of computation time in the supermolecule approach. An attempt to simplify the supermolecule approach by treating the solvent molecules as fixed point charges on atomic centers was reported by Noel1 and Morokuma.12 At the present stage this point charges model is not useful for practical calculations since it keeps the solvent molecules in arbitrarily fixed positions. The energy of the model is expected to collapse to negative infinity upon minimization since it does not include the proper repulsive terms. A promising approach has emerged from formal treatments of polar solvents13which provide asymptotic expressions for the mean force between ions. However, such approaches cannot provide the solvation energies of ion pairs which are crucial for analysis of ionic reactions. Probably the ultimate promising approach for studying solvent effects will be the brute force method of computer simulation of the given solution. This can be based either on Monte Carlo14 br molecular dynamics calculations including explicitly the solute and many solvent molecules. The potential energy of the system can be represented by analytical potential functions obtained either by quantum mechanical15or by empirical fitting procedures.16 While such approaches have been applied successfully to pure polar solvents16they are still far too expensive for accurate evaluation of solvation energies.l’ The main problem is that the multidimensional potential surface of the solution possesses too many local minima with quite different energies. Such a surface cannot be scanned efficiently by

The Journal of Physical Chemistty, Vol. 83,No. 12, 1979

1041

the available simulation approaches. Here we present a new method18 which can be considered as a refinement of our recent microscopic dielectric m e t h ~ d . This ~ ~ ?method ~~ retains the basic philosophy of the simulation approaches by explicitly including in the calculations many solvent molecules. However, in order to make the calculations practical it introduces the following drastic simplifications. (i) The dimensionality of the problem, and therefore the number of local minima of the solution system, is reduced by representing the solvent molecules as soft sphere point dipoles. This simplification is somewhat similar to the one introduced in a recent computer simulation of protein folding.21 (ii) The solvation enthalpy is approximated by the difference between the minimum potential energy of the solventsolute system and that of the bulk solvents. (iii) As a boundary condition the centers of the solvent molecules of the outer solvation shells are fixed in the positions they have in the bulk solvent. The interaction between this surface of fixed-center dipoles and the rest of the system is included in the calculations. These simplifications provide a model which includes explicitly the solvent molecules (thus avoiding the problems of the continuum models) and yet represents them in a sufficiently simple manner to allow practical studies of solvent effects. The applicability and reliability of the model are demonstrated here in studies of various chemical processes in solutions. Section I1 describes in detail our theoretical model. In the first part of section I11 we evaluate the energy balance of charge separation in polar liquids and use the observed dissociation enthalpies of dicarboxylic acids and amino acids as an experimental verification of the calculations. The rest of section I11 presents calculations of different types of chemical processes in solution. This includes evaluation of potential surfaces of ionic reactions in solution, determination of acid dissociation enthalpies, evaluation of solvent effects on conformational energies, and calculations of excitation energies. Section IV discusses the validity of our model and considers various possible applications. 11. Theoretical Model

A. General Considerations. Most previous treatments of the dielectric effect1a,22 considered the relation between the macroscopic polarizability of the bulk solvent and the magnitude of the dipole moment of the solvent molecules. From a chemical viewpoint the microscopic interaction between the solute and solvent is of much direct interest but it is not obvious how to evaluate this interaction using the available dielectric theories.22 In fact, a consideration of the average interaction between nearest neighbor solvent molecules22bmight be impossible when the solvent molecules are strongly polarized by the solute molecule. It still seems that the solute-solvent interactions can be treated (in some cases) qualitatively by the reaction field approaches.’“ However, these approaches imply that the solvent around the solute is an effective continuum, an assumption which ignores the fact that the size of the solute and solvent molecules is comparable. The value of this approximation is therefore highly questionable.2s In this work we decided to abandon the macroscopic dielectric approaches and to evaluate the microscopic energy contributions in the solute-solvent system. Since the dependence of the solvation energy on the solute coordinates (R) is of interest, we start formally by considering R as adiabatic coordinates and expressing the classical partition function for a given R as Z(R) = A ~ ~ / 2 S e x p ( - V ( R , r ) / k T d rJ

( 1)

1642

7%

Journal of Phwlcar Ctmmisity, Vol. 83,No. 12, 1979

Arleh

Warshe1

where N is the number of solvent molecules, A is a constant, Vis the potential energy of the system, and r and T designate, respectively, the position vector and volume element in the space of the Cartesian coordinates of the solvent molecules. The corresponding enthalpy and free energy are given hy H ( R ) = k'ln d In (Z(R))/aT =

[ JV(R,r)

(V) +

exp{-V(R,r)/kq dr]/Z(R) + 3NkT/2 = 3NkT/2 = V(R,r,,) (AV) + 3NkT/2 (2)

+

G(R) = -kT In (Z(R)) where ro is the minimum energy configuration for the solvent molecules (at the given R), AV is the difference between V(R,r) and V(R,ro),and ( ) designates a Boltzmann average. In order to evaluate the solvation energy it is useful to consider a reference system consisting of the bulk solvent and the isolated solute molecule (at the given R). For this system the solvation enthalpy and free energy are given by AH@) = V(R,r,) - P"L(robU1k) [(Av(R)) - (AP"")l - VwiUdR)(3)

+

AG(r) = V(R,ro) - Vbdk(rohlk)- kT In [2(R)/Zba] VwiudR) where is the potential energy of the bulk solvent, and rObU1li is the minimum energy configuration of the bulk solvent. For a given model of the solvent and its corresponding potential V, H(R)can be evaluated using a Monte Carlo approach (e.g., the Metropolis approach) or by molecular dynamics calculations. However, these approaches are at present far too expensive to be used as general and practical calculation schemes. It might be reasonable to approximate AH by its leading terms (V(R,ro)- Vbdk(robdk) - VwlUk(R)).The validity of such an approximation was examined, for solvated H30*at 300 K by preliminary Monte Carlo calculations, by using the SCSSD potential which will be described below. These calculations involve finding ten different local minima and evaluating AH by sampling the neighborhood of each minima by lo5 steps in the Metropolis approach. The Boltzmann average of these ten AH values was found to be -105 kcal/mol where only -0.5 kcal/mol Contribution comes from ((AV) - ( A P k ) ) (the average contributions of (AV) and (AP"lk) were 6.5 and 7 kcal/mol, respectively). Additional support for the above approximation was obtained by performing many different minimizations of the potential of solvated H30+ with different initial conditions. It was found that for 300 K the Boltzmann averaging of the potential, over the configuration spanned, gave an average energy which is only 0.4 kcal/mol higher than the value of the potential at the lowest minimum. It appears that there are many minima in the potential surface of polar solutions with energy values quite close to V(R,r,,). This means the presence of some type of degeneracy which increases the statistical weight of the contributions of V(R,r0).%It seems to us that in cases of polar solutes (V(R,r,,) - VbU1k(robdk)) is the leading term even for AG(R). In fact, for aqueous ionic solutions at room temperature the solvation enthalpy contributes several hundred kilocalories per mole to the solvation free energy while the entropy contributes less than 10% of the total value of the solvation enthalpy. We therefore argue that in cases of polar solutes the entropic contribution to solvation free energy are merely second-order effects that should be studied after (rather than before) obtaining a

Figure 1. The definition of the dilferent regions in the SCSSD model. Regions I, 11. and Ill are, respectively. the solute region, the strong field solvation shell region, and the bulk surface. Region I V is the surrounding continuum. The solute in the ligure Is an H,O+ ion.

better understanding of the huge first-order effects of solvation enthalpies. In view of the above arguments we took a semiempirical approach, chwsing a functional form for V and calibrating (V(R,ro)- V'"lk(robdk)) to reproduce solvation enthalpies at 300 K. The inductive approach to selecting a functional form for V involved examination of several microscopic models of increasing sophistication, beginning with the semiempirical model of ref 19, referred to here as the "fixed center Langevin dipoles" (FCLD) model. It became apparent (as will be discussed in detail in section 111) that in order to reproduce solvation energies of ions with different charges it was essential to include both orientational and translational degrees of freedom of the solvent molecules. Furthermore, it was found that in order to reproduce the energetics of charge separation it was essential to introduce some type of surface boundary condition to mimic the effect of the bulk phase surrounding the inner solvation shells (see section Ma). Such a boundary condition is introduced in the "surface constraint soft sphere dipoles" (SCSSD) model which is presented in this work. All the above discussion concerned mainly polar solvents. Treating ionic processes in nonpolar solvents appears to be quite simplelg and is described in section IIC. B. Treatment of Polar Solutions by the SCSSD Model. ( I ) The Main Features of the Model. The SCSSD model is described schematically in Figure 1. The model is designed as follows: the solute molecules (region I) are surrounded by a cluster of solvent molecules (region 11). Each solvent molecule is represented by a point dipole attached to the center of a compressible soft sphere. The solvent molecules in region I1 are surrounded by a surface of solvent molecules (region 111) constrained to fixed centers which correspond to the structure of the bulk s0lvent.2~The complete system of regions I, 11, and I11 is embedded in a continuum with the dielectric constant of the hulk solvent. The interaction between the pairs of solvent molecules is represented by dipole-dipole and van der Waals potential functions. The interaction between the solvent

The Journal of Physical Chemistty, Vol. 83, No. 12, 1979


TABLE I: Potential Parameters of the SCSSD Model parameter value parameter value e*(H,O), kcal/mol r*(H,O), a

0.16 3.66 1.88 6.0 1.6 0.8

PO

ra, A

AP’(H,O-H,O), D Afi’(H,O-solute), D e*( H), kcal/mol

C), kcal/mol e*( N), kcal/mol e*(

I

0.018 0.067 0.136

TABLE 11: Propertiesu Included in the Refinement of the Parameters of the Solvent-Solvent Potential obsd value

calod value

-8.6 -8.0 -9.1 -9.4

-7.7 -8.2 -8.8 -9.1

property

e*(O), kcal/mol e*(O-), kcal/mol

0.052 0.145

1643

cluster potential ( V N / N ) ~ N = 20 N = 40 N=80 N = 100 N = 120 no. density,

1/30

1/33

molecules and the solute atoms is represented by dipole-charge and van der Waals potential functions. Thus the total potential energy of the SCSSD model is given by

Vj is the potential of a cluster a Energy in kcal/mol. of N water molecules. The corresponding experimental estimate is given in eq 16.

V(R,r) = VsoiudR)- Z Q t r i l ~ j / ~+z ~

the external surface. This procedure is described in the Appendix. A general understanding of the results of this work does not require reading the detailed description in Appendix I. However, the notation in some of the tables is based on the notation presented in the Appendix. (C) Nonpolar Solvents. The treatment of nonpolar solvents with known structure (e.g., proteins) is described in detail in ref 19. For general organic solvents the treatment involves the following steps: (1)The solvent atoms are placed on a cubic grid with a spacing which reproduces the solvent density. Any solvent atom which is within van der Waals distance from the solute atoms is deleted. (2) Point-induced dipoles are assigned to each solvent atom by the relation Pl = ffE, (7)

41

cP][v(PfrJ]’/r]]’3)]

Pi’

+ cvVdW(rlj) + 11

vVdW(r]Y)

+I’

+

Vcontinuurn

(4)

where V,, is the intramolecular solute energy (see section II.D), the Qi are the solute atomic charges, and the 1.1, are the solvent dipoles in regions I1 and 111. Vvdwis the van der Waals potential function and Vconthumis the solvation energy of the charges in region I by a continuum (region IV) with the bulk dielectric constant. (For details see next section.) As discussed before, the solvation enthalpy is approximated by the minimum values of the corresponding potentials. That is, we use the relation = V(R,ro) - Plk(rgbulk) - VsoiutJR)

(5) where Vbdk is composed of the third and fifth terms of eq 4 and roand robdkare the minimum energy configurations of the solution and the bulk solvent, respectively. Equation 5 expresses the overall solvation energy. This energy includes the charge-dipole and dipole-dipole energies as well as the solute-solvent van der Waals energy needed to create the solvent cavity around the solute. For studies of ionic reactions it is convenient to define the total energy of the solute-solvent systems relative to their reference state by AVsolvation(R)

AVtot

=

AVsoluta(R)

+ AVsolvation(R)

where a is the atomic polarizability, and E,is the field on the solvent atom resulting from the solute charges and the other solvent dipoles. Equation 7 is solved iteratively starting with the field resulting from the solute charges and including the field due to the other dipoles in the second iteration. This treatment converges very fast (approximately seven iterations) and does not lead to computational problems. The resulting potential is

v = vsolute

-

+

~ & i P ~ ~ z j / ~ i j ~

11

(6)

where AV,olut,(R) is the energy needed to change the geometry of the isolated solute molecules from a reference geometry Ro to the given R. (2) Computational Details of the SCSSD Model. The practical computational aspects of the SCSSD model are presented in detail in Appendix I. The corresponding computer program is available upon request from the author. As described in the Appendix the model uses van der Waals, dipole-dipole, and dipole-charge potential functions. The dipole-dipole and dipole-charge potential functions include an expression for the deviation of the point dipole approximation from the behavior of a finite dipole. As discussed in the Appendix, the potential parameters used (Table I) were refined by fitting various calculated and observed properties. The properties used to refine the parameters of the solventsolvent interaction potential are listed in Table 11. The properties used to refine the single undetermined parameter in the solute-solvent interaction potential are the hydration enthalpies of monoatomic ions given in Table V. The practical implementation of the SCSSD model requires a special computational procedure for constructing

where the fi1 are the final induced dipoles of the solution atoms and the Qi are the charges of the solute atoms. Notice that classically the atoms can be considered as conducting spheres whose volume is equal to the corresponding atomic polarizability. If two spheres, representing two atoms, overlap (as in cases of bonded atoms) then the dipole-dipole interaction should be eliminated from the calculation. D. T h e Solute Molecules. The energy (Vsolute)and atomic charges (Q,) of the solute molecule can be evaluated by any LCAO-MO-SCF quantum mechanical treatment. The crucial modification of the quantum mechanical calculations is the incorporation of the effect of the solvent dipoles on the solute energy (the reduction energy in Sinanoglu’s theory). Instead of using Onsager reaction field consideration^,^ it is more straightforward to modify the core potential (the effective ionization potential) of each solute atom by the potential from the solvent ,dip o l e ~ That . ~ ~ is, the diagonal SCF matrix element p;”is modified by =

+

solution

C pjril/riF I

where v designates orbitals, i designates atoms,

(9)

is the

1644 The Journal of Physical Chemistry, Vol. 83, No. 12, 1979

Arleh Warshel

TABLE 111: Calculated Energy Balance of Charge Separation in Watera

R

VQQ

AVQ,

2.5 3.0 3.5 4.0 5.0 6.0 7.0

- 133

-111 -95 -83 - 66 - 55 - 48

-112 - 137 - 161 - 183 - 232 - 266 -311

26 27 32 27 76 100 137

4 8 13 29 15 17 22

-2 -3 -4 -5 -7 -10 -14

-215 -214 -214 -213 -213 -213

- 184 - 197 - 206 - 208

- 194 - 204 - 210 - 216

0

- 368

188

4

- 40

-213

- 210

- 206

m

Av&v

AVpp

Vcontinuum

AVtot SCSSD

AvFCLD tot

AVf$TD

- 217

- 189 - 187

- 199 - 198 - 194

- 186

a Distance in A , energies in kcal/mol. R is the charge separation distance and VQQ is the charge-charge interaction. A V Q ~Ab’, , , AVvdw, and Vcontinuum are, respectively, the charge-dipole, dipole-dipole, van der Waals, and continuum contributions t o V - Vo in eq 19. AVfZsSD includes all the energy contributions of the previous columns and a small contribution of 1-2 kcal of ( + $Vc). AVFCLD tot and AVf:fD are the total energies obtained by using the FCLD and SSPD models. The calculations were done for two ions with the size of water molecules (the same r* and e * ) and with opposite charges. The solvent includes 64 dipoles in region I1 and 1 6 dipoles in region 111. The values for = are twice the corresponding contribution for an isolated ion.

Vzption

diagonal SCF matrix element in vacuo, and the pj are the solvent dipoles. Combining this modification of the SCF calculation with the SCSSD model makes it possible to perform quantum mechanical calculations for chemical processes in solution. In the next sections we will demonstrate the applicability of our approach to different phenomena in which the solvation effects are of major importance. 111. Results A. Charge Separation in Polar Liquids. One of the most fundamental questions concerning the microscopic nature of the dielectric effect is how the electrostatic energy is so uniquely balanced in polar liquids. The puzzling fact is that in polar liquids separation of charges from 4 A to infinity costs practically no energy. Clearly this effect will be obtained if one uses Coulomb’s law and divides the electrostatic energy by 80. However, the basic question is, how the loss of -83 kcal of charge-charge electrostatic energy is being almost completely compensated for by the solvation energy. The compensation effect is not due to the reduction of the field between the solute charges by the “reaction field” from the surrounding dipoles (as might be understood from macroscopic dielectric approaches). Here we try to elucidate the mechanism of the compensation effect by detailed examination of the microscopic energy balance of charge separation. The process of examination involves comparing different microscopic approximations with increasing degree of sophistication. Table I11 summarizes these calculations and presents the different energy contributions as a function of the charge separation distance (R). The simplest model used is the “fixed center Langevin dipoles” (FCLD) model of ref 19. This model gives an artificial “activation energy” of -24 kcal/mol. The next approximation involves replacing the FCLD effective Langevin dipoles by permanent point dipoles. Although this model (referred to here as a “fixed center point dipoles” (FCPD) model) includes explicitly the dipole-dipole interaction term. the results obtained are similar to that of the FCLD model. The use of a “soft sphere point dipoles” (SSPD) model (the SCSSD model wkhout the surface constraint) lowers the activation barrier and makes it much smoother, However, this model gives a minimum in the energy at R 7 A, 10 kcal/mol less than the energy at R a. Apparently upon bringing two solvated ions together the chargecharge attraction increases faster than the solvation energy decrease, indicating that the two processes are

-

-

unrealistically unbalanced. When the distance between the two ions becomes shorter than 6 8,the solvent dipoles in the strong field region between the ions are expelled and the energy of the system increases leading to an activation barrier at R 4 A. The same type of behavior is obtained even when 200 solvent molecules are included. This effect is due to the fact that, for a limited number of solution dipoles, even dipoles in a weak field region are polarized toward the solute charges. Although the SSPD model is a simplified model it seems that this model simulates the real interaction in a system of an ion pair solvated by a limited number of water molecules. In this respect it will be very interesting to examine experimentally the geometry of such gas phase solvated ion pairs. The SSPD model does not appear to simulate the energy balance of an ion pair in a solution of many solvent molecules. That is, relaxation rnea~urements~~ of ion recombination reactions have shown that the rate of these reactions is diffusion controlled so that the corresponding activation barrier must be small. In order to simulate the situation in real solution we introduce the SCSSD model. This model constrains the surface dipoles to be in the positions of the bulk solvent thereby resolving the problem of excess stabilization of an ion pair at medium distances. The results obtained by the SCSSD model are given in Table I11 and Figure 2. Figure 2 also presents calculated equilibrium geometries of the solute-solvent system along the charge separation path. As seen from the figure the solvent dipoles in region I1 cannot be completely oriented toward the solute charges under the constraint of the bulk orientation in region 111. With the surface constraint modification it appears that the microscopic dielectric effect is due mainly to the availability of the many degrees of freedom which can balance any electrostatic force and maintain a nearly constant total energy. That is, the forces created by a charge separation from R to R + hfi can be compensated by a small reorientation and compression of the solution molecules without significant change in the total energy of the system. It should be noted that when the solvent molecules were represented by spheres with a radius smaller than that of the solute ions the results of the SSPD and the SCSSD model were similar, In this case the model simulates the continuum models. Introducing very small solvent molecules is equivalent to numerical integration over volume elements of a continuum solvent. However, it seems that the reason that the E in Coulomb’s law ( V = QiQj/(q,)) works on a microscopic scale is not that the solvent

-



1

I

I

I

I

I

I

2

3

4

5

6

7

R(%)

1645

I

a

Figure 2. Energy balance of charge separation in water. VaO and A VsDlvatbn are respectively the charge-charge and solvation energies. A Vtot is the calculated total energy as obtained by the SCSSD model. The upper part of the figure shows the calculated orientation of the solvent point dipoles for charge separations of 3, 5, and 7 A, respectively.

TABLE IV: Calculated and Observed Energies of Charge Separation in Dicarboxylic Acids and Amino Acidsa malonic acid succinic acid adipic acid glycine y-aminobutyric acid e-aminocaproic acid

3.7 4.5 6.0 3.5 4.7 5.6

87 74 55 - 92 -70 - 59

- 140 - 140 - 139 - 144 -141 - 140

- 70 - 69 - 68 -70 - 68 - 68

2 1 1 -2 -1 0

1.2 0.7 0.3 - 2.9 - 1.1

(3.12) (1.09) (0.52) (-2.15) (-0.98) (-0.52)

a Units as in Table 111, i? is the distance between the charge centers. VQQis the charge-charge energy, A V ( ' ) is the difference between AVtot (the sum of the solvation energy and VQQ)of the doubly ionized molecule and the solvation energy of the un-ionized molecule. A V(') is the difference between the solvation energy of the singly ionized molecule and the solvation energy of the un-ionized molecule. A A Vc&d is the calculated difference between the dissociation energy of the singly ionized molecule and the dissociation energy of the un-ionized molecule. For dicarboxylic acids A H o b d is the difference between the dissociation enthalpies (at 300 K ) of the singly ionized and un-ionized molecule.26 For amino acids A H o b d is the difference between the dissociation enthalpies of the amino acid (-0OC. . .NH:) and the corresponding ester (ROOC. . . The values in brackets are the corresponding A G o b a at 300 K evaluated from the corresponding ApK,. NH:).

molecules can be considered a continuum but rather that the solvent has sufficient degrees of freedom to compensate for any change in charge-charge interactions by rearrangement. The large compensation of charge-charge interactions by solvation energy can be verified by considering the effect of a charged group on the dissociation enthalpy of a neighboring group. Table IV shows this effect for several dicarboxylic acids and amino acids. As seen from Table IV both the experimental and theoretical results indicate that in aqueous solution the solute electrostatic interac-

tions are almost completely compensated for by the solvation energy even at short distances. Apparently the effective dielectric constant in water is greater than 20 even for a 3.5-A separation. It should be pointed out here that the macroscopic dielectric theory of Kirkwood and Westheimerlbgives quite successful estimates of the ApK,'salB of dicarboxylic acids. This might mean that continuum models which involve empirical parametrization of the cavity size can be used to evaluate effective dielectric constants at short distances. However, this success may be fortuitous. Any theory which

1646


Arieh Warshel

TABLE V : Calculated and Observed Hydration Enthalpies (kcal/mol) of Monoatomic Ionsa charge

rion

(r*/2)

AV D t::

A Gobsd

AHobsd

vFCLD calcd

138 126 +1 0.60 (0.8) 133 (123) t1 1.4 (1.8) 86-122 (81-114) 110 105 -1 1.4 (1.8) 113 (104) 110 105 t 2 0.7 (0.8) 520-477 (400-455) 540 452 t 2 1.4 (1.9) 320 (306) 373 330 t 3 0.6 (0.8) 1147-1073 (1106-1035) 1160 758 t 3 0.9 (1.8) 1028-920 (976-890) 1004 660 The calculated ionic radius is given (in A ) by rion = req(I) - 1.4 where rea(’) is the calculated average radius of the first solvation shell. The values of the parameters r* used are given in parentheses, the values of the corresponding E * are 0.16 kcal/mol in all cases. AVSCSSD and AVFCLD are the calculated solvation energies using the SCSSD and FCLD models, respectively. The temperature used in the FCLD model is 300 K. The values of AHobsd and A G o b d are taken from ref 30 for representative ions with the corresponding radii.

evaluates the sum of the solvation and charge-charge contributions rather than the isolated contributions will work. Having a large effective dielectric will yield a small electrostatic energy which is the desired result. Furthermore, the continuum model assumes, rather than derives, reduction of the electrostatic energy with increasing distance. The continuum model does not provide much insight into the origin of the energy balance of the solute-solute system. Such insight is important both for a better understanding of the microscopic process and for confidence in the predictive power of the calculations. B. Solvation Enthalpies of Ions in Aqueous Solutions. Solvation enthalpies of charged species are clearly dominated by electrostatic forces. Thus the SCSSD model is expected to give good results when applied to the solvation energies of ions in aqueous solutions. Table V compares the calculated and observed30solvation entalpies of ions with different charges and radii. The solute-solvent van der Waals parameters are the e* of the water-water interaction and the r* which gave the indicated ionic radii. Table V also includes the results of the FCLD model. As seen from the table the FCLD model underestimates the solvation energy of polyvalent ions. This underestimation is due to the fact that polyvalent ions exert an enormous electrostatic force on the surrounding water molecule^.^^ This electrostatic force strongly compresses the first few solvation shells and produces a very large solvation energy. The fact that the SCSSD model allows for compressibility of the solvent dipoles is the reason for the drastic improvement in the results. The SCSSD results are encouraging considering the fact that only one parameter (Ap’ of the dipole-charge interaction term) was used to fit the solvation enthalpies of the polyvalent ions. C. Calculations of Ionic Reactions in Water. The dissociation of water into OH- and H30+is one of the most fundamental reactions in aqueous solution. Yet, to the best of our knowledge, there is no quantum mechanical calculation of this reaction. We chose the water dissociation reaction as a test case for our SCSSD calculations. Two water molecules are considered as the solute and are surrounded by the SCSSD solvent molecules. The energy of the solute molecules can be evaluated by any LCAOMO-SCF quantum mechanical treatment, modifying the diagonal SCF matrix elements using eq 9. Here we are more concerned with the solvent effect than with the reliability of the quantum mechanical calculation scheme. Thus, we used the MINDO/2 approach32to evaluate the energy and charge distribution of the solute molecules, making the following correction (10) Vsoiute = V:%$e + A v ( 1 - p ) where Vi,$$ is the calculated MIND0/2 solute energy, AV is the difference between the calculated and observed energies of the gas phase proton transfer reaction (HzO +

-

HzO OH- + H30+),and P is the bond order of the OH bond which is broken in the reaction. This type of correction ensures proper behavior in the asymptotic region and leaves the exact shape of the reaction potential surface to more accurate calculations of the solute energy. In addition, a useful empirical estimate of the energy of the ionic states of the solute is given by pBoiute

= CQi&,/r,I + A +

pvdw

(11)

41

where the Q’s are the atomic partial charges of the solute fragments in the given ionic state at infinite separation, A is the observed gas phase proton transfer energy, and pv:dw is an empirical function that represents the van der Waals nonbonded interaction between the solute fragments. This approach is expected to give an accurate potential surface for ionic states at intermediate and large separations of the solute fragments (the small overlap regions). The results of the calculations are summarized in Figure 3 and Table VI. The figure and the table present the total reaction energy and the separate contributions of the solute and solvent. As seen from the figure the total energy is almost constant when the charge separation distance ( R ) is larger than 4.7 A. For shorter distances the MINDOIB potential gives a significant activation barrier while the empirical potential (F‘)gives a more realistic result. We believe that the empirical potential gives the correct representation for the ionic charge transfer state as long as R > 3.5 A since in this region the van der Waals interaction between H’ and O2 is small. A more accurate evaluation of the potential surface for R between 3.5 and 2.7 A might require further work. It is apparent that the MINDOIB potential overestimates the energy of the charge transfer state in the region 3.0 5 R 5 4.5 A. This is probably due to an improper balance between the electron-electron repulsive integral and the core potential (the MIND0 refinement did not include experimental information about this weak overlap region). In view of the success of the empirical potential P it might be worthwhile to calculate the potential surfaces for ionic reactions in solutions using an empirical valence bond approach in which the diagonal matrix elements of the covalent states are given by Morse type potential functions and the matrix elements of the ionic states are given by (Vsolute+ AVso~v~,~~,). A related treatment of the hydrogen bond in the gas phase was introduced by Coulson and Daniel~son.~~ It should be pointed out that we did not address ourselves to the question of whether the recombination of OHand H30+is accomplished via “proton hopping” between the different oxygen sites2Ior by a direct recombingtion of the solvated ions. However, the configurations R = 3 A and R = 5 A are identical in both cases (the configu-



1647

TABLE VI: Calculated Energy Balance in Water Dissociationu

R

b

QH,o~

Vsolute

0.20 0.70 1.0 1.0 1.0 1.0 1.0 1.0

10 0 15 63 127 128 148 159 166 220

2.75

1.10 1.10 1.40 1.65 1.65 2.6 3.6 4.6 5.6

m

2.95 2.75 2.75 3.75 4.75 5.75 6.75 cc

m

VQQ+ A

AVsolvation

AVtot

- 121 - 146

18

25 25 22

- 166

15 16 16

18 18 18

- 46 95 128 150 160 166 220

AqOt

8 0 9 33 45 23

- 18 - 16 - 22 - 96

- 160 -220

0

fz are the 0,-H'

and the 0,-0, distances, respectively (where H' is being transferred from 0, t o 0 2 )QH$o+ , is the net positive charge on the H 3 0 t fragment. AVsolute is the MINDOI2 solute energy (corrected using eq 1 0 ) relative t o the energy of two separated H,O molecules. When Q H ~ Q +is equal to unity the MINDO/2 calculations are constrained t o given the energy of a pure charge transfer state (by fixing the charge density and not doing SCF iterations). VQQis the calculated charge-charge energy using the charge distribution of the system for = -, and A is the observed enthalpy of gas phase proton transfer. AVsolvationis the calculated solvation energy. The MINDO/2 total energy, AVtot, is given by Avtot = Vsolute + (AVsolvagon + 1 6 ) where the 16 kcal/mol contribution gives the results relative to the solvation energy of two water molecules at R = _. The total empirical energy, A F o t r is given by A v b t = VQQ + A + v v d w + (Avsolvation + 18) (see eq ll), where the 18 kcal/mol contribution gives the results relative to the solvation energy of a water dimer. a b and

TABLE VII: Calculated and Observed Dissociation Enthalpies in Solution and in Gas Phasea acid

(AvET)calcd

H,O CH,OH HCOOH

16 10 5

)obsd

13(19) ll(22) O(5)

(AH8T)obxi

(Av$T)eq 12

(AV8T)MINDO/2

220 211 177

217 212 172

218 190 147

(AvET)c&d and (Aff$T)obd are, respectively, the calculated and observedz6enthalpies of proton transfer (PT) from the indicated acid to water molecule in aqueous solution (S). (Ai?$T)obd, ( A V$T),g12, and ( A vfpT)MINDO/2 are, respectively, the observed gas phase proton transfer enthalpies33+% and the corresponding calculated values using eq 12 and using MINDOI2. In parentheses are the observed free energy for the proton transfer reaction in solution.

rations R = 3 A and R = 5 A in the direct recombination model correspond in the proton hopping model to the proton being on the first and second neighbor sites, respectively). Thus it seems that only kinetic effects make the proton jump model more probable than the direct recombination model. D. Calculations of Dissociation Enthalpies in Solution. The SCSSD model offers the possibility of quantitative evaluation of the factors which determine the stability of ionized groups in solution. To demonstrate this point we calculated the dissociation enthalpies of several typical acids using eq 10, where the proton transfer reaction involved is RH + H 2 0 R- H30+. A comparison of the calculated and observed acidities of HzO, CH30H, and HCOOH are given in Table VII. The table shows good agreement between the calculated and observed dissociation enthalpies considering that the results involve a balance between very large contributions. In view of the effort needed to determine the gas phase acidity difference between H20 and HCOOH33s34it seems that our calculations can be used as a preliminary estimate when gas phase measurement are not available. Knowing that the solvation calculations give the correct order of magnitude one could use the relation A V ~ T A@T - AVsolvation (12) where A@, is the gas phase proton transfer energy, A@, is the observed proton transfer enthalpy in solution, and AVsolvationis the calculated difference in the solvation energy of the products and reactants. The accuracy of this approach is demonstrated in Table VII. E. Conformational Energies in Solution. The use of the SCSSD model to evaluate the conformational energies of polar molecules in solution is, in principle, straightforward. The energy of the complete solute-solvent system is simply calculated for different conformations of the solute. This point can best be demonstrated by consid-

-

+

Figure 3. Dissociation of water In aqueous solution. The 0 line ( Vme) is the MINDO/2 solute energy (corrected using eq IO). I n the region 2.75 < R < 4.5 A the MINDO/2 solute energy is described by two 0 lines: the lowest one Is the ground state energy obtalned using eq 9 and a SCF procedure, the upper one is the MIND0/2 energy of the pure charge transfer state (see Table VI). The iline (P ),, Is the empirical estimate of the energy of the solute charge transfer state. A V,,~a,,tlon - 200 is the solvation energy relative to the solvation of a water dimer (-18 kcal/mol). The potential surface AV,,, Is the sum of the MINDO/2 SCF energy and the corresponding solvation energy. A PI,, Is the total empirical energy of the solvated charge transfer state. is the distance between 0,and O2and b' is the H'---O, distance.

ering the stacking of the cytosine dimer in solution. The role of the electrostatic interaction in parallel base stacking

1048

The Journal of Physical Chemistry, Vol. 83,No. 12, 1979

Arieh Warshel

TABLE VIII: Parallel Stacking of Cytosine Monomers in WateP

-

Z,

AVQ, - 79

AVsolute

3.5 4.0 5.0 6.0

9 6 3 1

m

0

- 77

- 76 - 80 - 52

AViW

Avvdw

Vcontinuum

qNc

45 44 44 49 24

-3 -3 -1 1 -2

-5 -5 -4 -4 -2

6.6 6.8 7.3 7.6 9.9

V;;l?tion

-7.7 -7.8 - 8.3 -8.9 -9.8

AVtot

2 0

1 1 4

a Energy in kcal/mol, distances in A , notation as in Table I11 and eq 19. z i s the distance between the monomer planes where the relative orientation of the monomers is given in Figure 4. AV is given relative to AVtot for the 2 = 4 A case. AVsolute and the solute charges were evaluated by the QCFF/PI method.&’

‘4

+‘

5.0 6.0 7 .O Z(%) Figure 4. Energy balance of parallel base stacking in solution. The figure presents the intramolecular solute energy ( VsOste)and the total solute-solvent energy (A Vbt) for the cytosine dimer, which is oriented as shown in the upper part of the figure. 3.0

4.0

(especiallyin the case of the cytqsine dimer) has never been explained s a t i s f a ~ t o r i l y .The ~ ~ ~cytosine ~ dimer consists of two very large dipoles ( - 6 D each) which cannot be stacked in a parallel configuration according to most dielectric theories.36 Recently we studied the stacking of cytosine monophosphate ( C P C ) ~using ~ the QCFF/P13s method to evaluate the solute potential surface and charge densities, and using the SCSSD model to evaluate the solvation enthalpies. This work shows that parallel stacking can be accomplished due to compensation of the solute charge-charge interaction by the solute-solvent interaction. This point is demonstrated in Figure 4 and Table VIII. The figure and table present the energy of the cytosine dimer as a function of the distance between the monomer planes. As seen from the figure and table the calculated dimer energy surface has a minimum at about 4 A. The calculations are obviously not yet reliable enough to allow accurate prediction of both the minimum energy position and its depth but they clearly demonstrate that stacking of the parallel dipoles is a realistic possibility. It should be noted that in the real CpC system the ribose phosphate bridge holds the bases in a much less parallel orientation than in the present calculation. In this case the interaction between the charges of the two bases is not as repulsive as in the case of the isolated bases. Thus, the role of the solvent in allowing for parallel stacking is less crucial in the CpC case than in the dicytosine case. F. Solvent Effect on Electronic Transitions. The customary approach for evaluation of solvent effects on electronic transitions relates the “solvent shift” to the solute transition moment and the radius of the solvent ~ a v i t y However, . ~ ~ ~ ~the dipole approximation might not

be valid in many cases and the effective cavity radius is not well defined. When the system involves the solvation of a chromophore and its counterion by a polar solvent41 the calculations become even more questionable. Here we offer a straightforward and much more reliable approach for evaluation of solvent shifts. This approach involves calculation of the energy of the solute-solvent system in both the ground and excited states. That is, the energy of the transition between the ground state, G, and the Mth excited state is evaluated by A ~ G - M= VM= E l u t e - E o l u t e + Avkhation A E l v a t i o n (13) The calculation for nonpolar solvents42involves the following steps: (1) A cubic grid of solvent atoms is built around the solute (see section 1I.C) and induced dipoles with the solvent atomic polarizability are assigned to each grid point. (2) The charge density and energies of the solute ground state are evaluated by the QCFF/PI met h ~ d (3) . ~ The ~ orientations of the induced dipoles and their interactions with the solvent charges are evaluated as described in section 1I.C. (4) The solute charges and energies are reevaluated using eq 9 and step 3 is repeated. (5) Steps 1, 3, and 4 are carried out for the excited electronic state and eq 13 is evaluated. In a nonpolar solvent the induced dipoles of the solvent respond instantly to the corresponding excited state charge redistribution, and eq 13 gives directly the excitation energy of the solute-solvent system. In polar liquids, on the other hand, the solvent dipoles cannot reorient themselves instantly. Therefore, while AV$lv,,i,, is obtained by the SCSSD model the solvation of the excited state is evaluated by keeping the permanent dipoles in their ground state SCSSD configurations and allowing only the induced dipoles (which are otherwise omitted) to be reoriented. In this case one has to attach induced dipoles with the average polarizabilities of the solvent molecules to the centers of the SCSSD dipoles (e.g., a = 1.4 A-3 for water). The orientations of these induced dipoles are determined by the method of section 1I.C taking into account the field from the solute charges as well as the solvent permanent and induced dipoles. The application of our approach is demonstrated in Figure 5 and Table 1X for the protonated Schiff base of retinal in a nonpolar solvent. A part of this calculation is discussed in ref 42 and 43 and more details will be given elsewhere. As shown in Table IX, the solvent contributes in two ways to the change in the excitation energy. First, the solvation energy changes on excitation (AAVs,~vation). Second, the “reaction field” from the solvent dipoles affects the ground and excited state energies of the solute in different ways (AAVs,,ut,). The evaluation of the second contribution is quite difficult when one uses a continuum model for the solvent. Here it is evaluated as the change in the ground and excited state intramolecular solute energies as a result of including the solvent dipoles in the SCF matrix (eq 9). As shown in Table IX, the change in



TABLE IX: Solvent Effectsa on the n

-

R 3.0 4.0 m

-f

n* Electronic Excitation Energies of Protonated Schiff Base of Retinal* n+n* AVn*n*

AVES AAVsolute 2 1 167 19 250 1 8 890

1849

36 30 20

AVSvation

AVsnoTvation

-1877 -4690 - 4660

- 2197 -4935 -4680

AAVsolvation

-320 - 245 - 20

tot

20 847 19 005 18 870

a The calculations were performed for a nonpolar solvent with an average atomic polarizability of 1 A - 3 and an average Distance in A and energies in cm”. R is the distance between the Schiff base nitrogen and the density of (1/8)atom/A3. are, respectively, the intramolecular excitation energy in solution and counterion (see Figure 5). AVZoT$Z and AAV;$ A Vare, ~respectively, ~ ~ the change in this energy between the solution and the gas phase. AV$l,ationr AV:Gvation, and A the solvation energies of the ground and excited states and the difference between these solvation energies. AVtot is the total excitation energy in solution. The charge of the counterion and the degree of protonation were determined by all valence electron calculations using OH- as a c ~ u n t e r i o n . ~ ~ GROUND STATE (G)

EXCITED STATE (r”)

-

Flgure 5. Solvent effect on the a T * electronic excitation energies of a protonated Schiff base of 1l-cis-retinal and a counterion. The Figure shows the charge densities of the ground and the a xstates and the correspondlng solvation energies (A Vsolvatlon), The change in the solute energy due to the back field from the solvent dipoles is included in the solute excitation energy. The counterion (OH-) and the Schiff base nitrogen are designated by 0 and 0,respectively. The Schiff base proton is placed between the nitrogen and the counterion.

solvation energy upon excitation is much larger than the reaction field contribution. It is important to notice that the general reduction of electrostatic interactions by the solvent is also valid when one deals with excitation energies. Thus, the shift in the absorption spectra of chromophores upon interaction with external charges is strongly overestimated by calculations which do not include the solvent.44

IV. Discussion This work presents a reliable and efficient method for the calculation of solvation enthalpies and related properties of polar solutes. The method is based on representing the solvent molecules by soft sphere point dipoles, minimizing the energy of the solute-solvent system with respect to the orientations and positions of the dipoles in the strong field region, and constraining the centers of the surface dipoles to stay in the positions they have in the bulk solvent. This explicit method avoids the longstanding problems of the continuum models and allows calculation of the microscopic energy balance of electro-

static interactions in polar solvents. The present work approximates the solvation enthalpy a t 300 K by taking the difference between the minimum potential energy of the solute-solvent system and the minimum potential energy of the bulk solvent. This approximation was examined in a preliminary way by performing Monte-Carlo calculations for a solvated H30+ ion.49 The difference between the Boltzmann averages of the potentials of the solvated ion system and the bulk solvent was found to be -105 kcal/mol. This value is only 0.5 kcal/mol lower than the result obtained by taking the difference between the minimum potential energies of the solvated ion system and the bulk solvent. The present work did not evaluate solvation free energies. However, in most cases of ionic systems the difference between solvation enthalpy and solvation free energy is less than 10% of the enthalpy involved. Thus, the calculations can be used as rough estimates of solvation free energies. Furthermore, it is possible to calibrate our ( V - PIk) to reproduce AG a t 300 K, as was done successfully in ref 48. The above discussion does not imply in any way that one should abolish fundamental studies of entropic effects, but only that solvation entropy is a second-order contribution to the solvation free energy of polar solutes and therefore can be absorbed into the leading term, the solvation enthalpy, by empirical parametrization. In this respect it is instructive to note that macroscopic dielectric theories use the bulk dielectric constant to evaluate free energies rather than enthalpies. Nonpolar solutes are quite different since entropic effects might contribute significantly to the solvation free energy. However, even such “hydrophobic” effects can be incorporated into the empirical calculation, using eq 19 and replacing the term qN, by aA where u is an empirical estimate of the surface free energy and A is the solute surface area. A value of CT = 0.025 kcal mol-’ (based on the results of ref 50 and 51) gives reasonable agreement between calculated and observed solvation free energies of nonpolar solutes in water. It appears that the SCSSD model possesses many minima of similar energy. This is the reason why it is sufficient to locate only a few local energy minima in estimating the value of the lowest one. The availability of many similar minima is not an artifact of the model. It is most probably an important feature of polar solvent systems which involve motion of the solvent molecules between different configurations along many shallow valleys of similar height. An important property of the SCSSD model is the ability of the calculations to converge rapidly to low energy minima without getting trapped at higher local minima. Such convergence is somewhat surprising for a system with so many variables. It is possible that the fact that we use a relatively smooth intermolecular potential for the solvent

~

~

1650


molecules is responsible for this desirable behavior. A similar simplification appears to be very successful in studies of the protein folding process.21 Molecular dynamics and Monte Carlo simulations of pure water give energies which are very sensitive to changes in the boundary conditions. In the present model, however, this instability does not exist. That is, changing the radius of region I11 from 6 to 8 A changes the solvation energy of a monovalent ion by less than 5 kcal/mol. Changing the shape of region I11 from spherical to ellipsoidal changes the solvation energy of an ion pair by less than 3 kcal/mol. The present work indicates that a very important feature of polar solvents is that any changes in the solute electrostatic energy is almost exactly compensated for by the solvation energy. The reason for this balance is due in part to the existence of sufficient degrees of freedom in the solvent system. Including the main degrees of freedom in the calculations, we find that in aqueous solutions the effective dielectric is larger than 20 when the distance between the solute charges is larger than 3.5 A. This means there is an almost complete compensation of the electrostatic interactions in most practical cases. The present model, and any other dipole model, cannot reproduce the hydrogen bond network of aqueous solutions. In fact, the solvent structures obtained by the SCSSD model more closely resemble methanol solutions than aqueous solutions. However, by adjusting the density and dipole length it is possible to constrain the dipole model to reproduce the solvation energies observed in aqueous solutions. This is not surprising since the energetics of charge separation in water and other polar solvents is similar, despite the significant differences in the corresponding structures. For example, although the dielectric constants of water and methanol are quite different (80 and 36, respectively) both dielectric constants are sufficiently high to give almost complete energy compensation for charge separation processes. Most available quantum mechanical schemes are designed mainly to treat isolated molecules. Our method makes it possible to extend quantum mechanical calculations to chemical phenomena in solutions. This point is demonstrated in this work by calculations of ionic reactions, electronic transitions, and base-stacking in polar solutions. The calculations show that the evaluation of the energy balance in the complete solute-solvent system is essential for quantitative results. This energy balance is only obtained when the energy of the system is allowed to relax with respect to the coordinates of the solvent dipoles. This indicates that the supermolecule quantum mechanical calculations of solvent effects might not be adequate (regardless of whether they are semiempiricalor ab initio) as long as they do not involve optimization of the positions and orientations of the solvent molecules. The calculations presented in this work are of a preliminary nature. However, it seems that a whole class of problems in chemistry and biology can now be studied in a quantitative manner. We list here some possible problems: (1)Many chemical reactions in solution are characterized by ionic transition states. The calculation of the corresponding solvation energy and the empirical solute potential (P in eq 11) can provide a detailed analysis of such transition states. (2) The problem of equilibrium conformations of polar molecules in solution was considered here only for the cytosine dimer. There are still many other unanswered questions concerning the conformational preferences of molecules of biological interest. In this respect it would be very interesting to apply the

Arieh Warshel

SCSSD model to the study of the equilibrium conformations of different peptides in solution. (3) The SCSSD model can be applied to problems which are related to the redox potential of transition metals in aqueous solution and in proteins. (4) The treatments of nonpolar systems1g can be extended to studies of ionic charge transfer cryst a l ~ This . ~ ~problem, which is of great practical interest, can be treated by incorporating the self-consistent polarization treatmentlg with crystal calculation approaches.4e ( 5 ) Applications to studies of enzymatic reactions were already described in a previous and further studies are now in progress. (6) Physical organic chemistry is one of the fields in which a microscopic dielectric approach can lead to a deeper quantitative understanding. The customary approaches in this field employ linear free energy relationships4' which correlate the reactivity of different systems with substitution effects. Such approaches essentially involve extrapolation of the sum of the solute and solvent energy without considering the nature of the balance between the different opposing energy contributions. Calculation of the energy balance of reaction in solution might clearly increase the extrapolative power of physical organic chemistry. This point is especially important when one deals with effects of charge stabilization in low dielectric systems. For example, recent microscopic dielectric calculations of enzymatic reactions showed that charge stabilization of transition states is the main catalytic factor in enzymatic reactions.48 This conclusion seems contradictory to the fact that no significant charge effect is observed in reactions of model compounds in solutions. This discrepancy can be understood in view of the results of the present work which demonstrate that in polar solutions one cannot see strong electrostatic interactions even at short distances. This is due to the ability of the solution dipoles to reorient themselves and minimize any electrostatic interaction. In enzymes on the other hand, the dipoles and charges of the active site are fixed in the orientation which gives maximum stabilization to the transition state.

Acknowledgment. This work was supported in part by grants from the National Eye Institute (EY01766)and the National Institutes of Health (GM 24492). Appendix I Computational Details of the SCSSD Treatment of Aqueous Solutions. ( a ) Representation and Potential Functions. The interactions between solvent molecules are described by dipoledipole and van der Waals potential functions (see the corresponding contributions in eq 4). The van der Waals potential functions were chosen as

Vvdw(rk])= 2t"((r*/rk])' - (3,/2)(r*/~kj)~)(14) where rkl is the distance between the indicated point dipoles, and where the refinement of the parameters e* and r* is described below. In order to represent the finite dipole properties of the water molecules by point dipoles, we replace the length of the dipole moment vector in the dipole-dipole interaction term of eq 4 by

]PI

= PO +

@'/(I

+ (rjy/rJ2)

(15)

The same type of function is also used for the dipolecharge interaction term (using different Ap'). For reasons of economy we chose to evaluate dipole-dipole interactions only for dipoles which are less than 5 A apart. The compensation for this cutoff is absorbed into the parameter Ap' in eq 15. In addition we introduced a term linear in rkJinto the expression for Vvdw of the solvent molecules


Calculations of Chemical Processes In Solutions

when ri’ < 2.7 A to prevent excessive “dimerization”. ( b ) dalibration of the Potential Parameters. The potential of the bulk water includes the five parameters r*, e*, po, Ah’, and ro (see eq 14 and 15). The parameter po was taken as the dipole moment of an isolated water molecule, Le., po = 1.88 D, and ro was fixed at 6 A. The parameters r*, e*, and Ap‘ were obtained by fitting the calculated and observed properties listed in Table 11. This parametrization include the following: (i) Fitting the calculated potentials of clusters with different number of water molecules to an estimate of the corresponding experimental enthalpies. The experimental estimate used is (in kcal/mol) V, = 10.4N - 0.1A (16) where N is the number of water molecules and A is the surface area (in A2) of the given cluster. This very rough estimate (which is expected to be reliable only for large N ) was obtained using the enthalpy of vaporization of the bulk water (10.4 kcal/mol) and the surface enthalpy of the bulk water (0.1 kcal/mol/A2). (ii) Fitting the calculated overage number density of large clusters of water molecules to the experimental number density of bulk water. The calibration of the solutesolvent interaction function includes only a refinement of the corresponding Ap’. This parameter was obtained by fitting the calculated and observed solvation enthalpies of ions of different charges and radii (Table V). The van der Waals parameters of the solute-solvent interaction were determined using the arithmetic mean law for r* and the geometric mean law for e*. The solute e* and r* are taken without change from ref 19. For each system studied the solute charges are calculated quantum mechanically. Different quantum mechanical methods are used for different systems. (The MIND0/2 method is used in sections 1II.C and 1II.D and the QCFF/PI method is used in sections 1II.E and 1II.F.) The parameters used in the present calculations are summarized in Table I. (c) T h e Operaton of the Model. The model involves formation of a surface of fixed dipoles whose centers are fixed in the configuration of the bulk solvent around the region of a strong solute field. The surface starts at a radius Rn which is determined, somewhat arbitraril~?~ by the relation (PCl[~IO(RII,T)l)= 0.2Po

(17)

where ( p ) is the average projection of the solvent dipoles (at a radius RII) on the direction of the field from region I, p, is the solvent dipole moment, and E t is the maximum field from region I at a distance RII. ( p ) is determined by the Langevin type expression:lg ( I * ) = po(coth ( 3 ~ ) - 1/x) (18)

x = E0l*o/(hT(1+ The surface is constructed by the following three stage procedure: (i) A cubic grid of the solvent dipoles is built around the solute molecules, deleting every solvent dipole within van der Waals distance of a solute atom. (ii) The solute charges (Q, in eq 4) are set to zero and the van der Waals and dipole-dipole energies of the solute-solvent system are minimized with respect to the positions and orientations of the solvent dipoles. This minimization (which does not include in this stage any solvent-solute dipole-charge interactions) forms a solvent ensemble with a structure close to that of the bulk solvent. (iii) The solvent system is divided into strong field solvation shell dipoles (region 11) and surface dipoles (region 111). The

1651

solute charges are reintroduced so that the total potential energy now includes the solute-solvent dipole-charge interactions. The complete potential is minimized with respect to the positions and orientations of the dipoles in region I1 and the orientations of the dipoles in region 111, leaving their positions fixed. Freezing the positions of the dipoles in region I11 is sufficient to constrain the orientations of these dipoles close to their bulk orientations, while the small reorientation of these dipoles is sifficient to provide the contribution of region I11 to the total solvation energy. The above three-stage procedure appears to provide stable minimum energies for the solute and the first few solvation shells under the bulk surface boundary condition. The solvation energy of the system of solute and solvent dipoles (regions I, 11, and 111) in the bulk solvent (region IV) is taken as the electrostatic energy of a sphere with radius RIII + 1.5 A and the charge distribution of region I imbedded in a continuum with the bulk dielectric constant. This calculation is done by the approach of ref 5. The solvation energy can be evaluated now by the difference between the minimum energy obtained by the above approach and the minimum energy of the reference bulk system. ( d ) Practical Simplification. The above procedure provides a powerful tool for the evaluation of solvation energies of charged species and for studying the energy balance in processes which involve charge separation. However, when the calculations involve comparison of different solute conformations with very different solvent cavity sizes (e.g., the case of base stacking which is presented in section 111) they become expensive due to convergence problems. In order to reliably obtain the dependence of the calculated solvation energy on the size of the solvent cavity it is often necessary to include a large number of solvent molecules. It was found that fewer solvent molecules can be included if the effect of the cavity size is cancelled out in the calculations and incorporated only in the final result by the relation AVsolvation(R) = V(R,ro) - VO(R,r:) + V9gation+ vNc (19) where VO is the potential obtained from eq 4 when the dipole-charge interaction is not included, V is the total potential of eq 4, and r: and roare the minimum energy configurations of VO and V, respectively. V$Ftion is the contribution to VO from the calculated van der Waals interaction between the solute and the solvent. N, is the number of solvent molecules in the first solvation shell. v is evaluated by calculating the solvation energies (AV,,o~v,tion) of nonpolar solutes of different radii and taking vsolvation )/Ncas 9. This the average value of (AVsolvation - vdw gives a value of q = 0.35 kcal/mol. ( e ) Minimization Procedure. The potential energy of the SCSSD model is minimized by the steepest descent method. The minimization uses the analytical derivatives of the potential with respect to the solvent coordinates (i.e., X , Y, 2, pz, py, and pz for each solvent dipole). For 60 solvent molecules the minimization converges after about 230 steps for both stages of the SCSSD procedure (see section c). Since the calculations may involve up to 600 degrees of freedom it is important to ensure that the energy minimum obtained is the lowest one possible. To determine this, the calculation is repeated about 10 times varying the initial coordinates of the solvent molecules. Also as a standard check for convergence we compare the minimum energy obtained in stage ii of the SCSSD procedure (VO) to an empirical estimate of the energy of

1652


a cluster of the same number of solvent molecules (eq 16). For minimization with different initial conditions the values of V(R,r,) change much less than the corresponding values of VO(R,r$). This is due to a compensation of changes of the solvent-solvent interactions by changes of the solvent-solute interactions. Because of this behavior it is necessary to repeat the minimization more times for evaluation of the lowest minimum of VO than for evaluation of the lowest minimum of V.

Arieh Warshel

(24)

(25)

References and Notes (1) (a) L. Onsager, J. Am. Chem. SOC.,58, 1486 (1936); (b) J. G. Kirkwood and F. H. Westheimer, J . Cbem. Pbys., 6, 506 (1938); (c) 0. Sinanogiu in “Molecular Association in Biology”, B. Pullman, Ed., Academic Press, New York, 1968, p 427. (d) For review see 8. Linder, Adv. Cbem. Pbys., 12, 225 (1967). (2) (a) 0. Sinanoglu and S. Abduinur, Pbotochem. Pbotoblol., 3, 333 (1964); (b) T. Halicioglu and 0. Sinanogiu, Ann. N . Y . Acad. Sci., 158, 308 (1969). (3) P. Claverie, 6. Pullman, and J. Cailiet, J. Theor. Biol., 12, 419 (1966). (4) M. J. Huron and P. Claverie, J . Phys. Chem., 76, 2123 (1972). (5) (a) D. L. Beveridge, M. M. Kelly, and R. J. Radna, J. Am. Chem. SOC.,96, 3769 (1974); (b) G. W. Schnueie and D. L. Beveridge, J. Pbys. Chem., 79, 2566 (1975). (6) J. Hylton, R. E. Christoffersen, and G. Hail, Chem. Phys. Lett., 28, 501 (1974). (7) For example, the solvation of a H30+-OH- ion pair separated by 3 A is about 100 kcailmol (see section 111) while the solvation as evaluated by the continuum models is about 55 kcailmol for a reasonable cavity radlus of 3 A. (8) A. Pullman and B. Pullman, 0 . Rev. Blophys., 7, 506 (1975). (9) W. L. Joraensen. J. Am. Cbem. Soc.. 78. 1057 (1978). (a) M. D. iewton, J. Chem. Phys., 58,’5833 (1973); (bj J. Phys. Chem., 79, 2795 (1975). (1 1) P. Koliman and I.D. Kuntz, J. Am. Cbem. SOC.,98, 6820 (1976). (12) J. 0. Noeil and K. Morokuma, Chem. Phys. Lett., 36, 465 (1975). (13) . . D. Chandler. J . Chem. Phvs., 63. 1113 (1977): . . S. A. Adelman and J. M. Deutch, ibid., 59, 3971 (1973). (14) J. C. Dwlcki and H. A. Scheraga, J . Am. Chem. Soc., 99, 7413 (1977). (15) E. Clementi, F. Ciavalione, and R. Scordamgiia, J . Am. Chem. Soc., 09, 5531 (1977). (16) G. H. Stillinger and A. Rahaman, J. Chem. Pbys., 60, 1545 (1974). (17) Recent extensive Monte Carlo studiesi4 of the solvation energy of CH, gave the discouraging result of ( A H ) = -1 1 A 13 kcai/mol (the observed solvation enthalpy is -2 kcal/mol). I t seems that trying to obtaln equal total energies for solvated ion pairs at 6- and 7-A separation is at present beyond the scope of the Monte Carlo approaches. (18) A preliminary account of this work was given in A. Warshei, Chem. Pbys. Lett., 55, 454 (1978). (19) A. Warshei and M. Levitt, J. Mol. Biol., 103, 227 (1976). (20) A. Warshel, Phil. Trans. R. SOC.London, Ser. B , 278, 111 (1977). (21) M. Levitt and A. Warshei, Nature (London), 253, 694 (1975). (22) (a) P. Debye, “Polar Molecules”, Dover Publications, New York, 1929. (b) J. G. Kirkwood, J . Chem. Phys., 7, 911 (1939); (c) A. D. Buckingham, Proc. R. SOC. London, Ser. A , 238, 235 (1956). (23) The basis of the continuum dielectric theories’ is the solution of Laplace’s equation for the electrostatic potential $ inside and outslde a cavity. The boundary conditions are (a $/ar)r=a-o= e(alE//ar),=,+,

(iOi

-

(26) (27)

(28) (29) (30) (31) (32) (33) (34) (35)

(36) (37) (38)

(39) (40) (41) (42) (43) (44)

(45) (46) (47) (48) (49) (50) (51) (52)

(see ref la) where a is the cavity radius. However, even this fundamental equation is only a rough approximation for a microscopic environment since t is not a constant. The high entropy contribution associated with the availability of many minima of similar energy does not appear In the observed solvatlon free energy since it is cancelled out in the difference between the bulk and solution contributions. The selection of the surface boundary condltlon represents the fact that in the weak field region the solvent molecules are randomly polarized. When the dipoles in region I1 are polarized toward the charges in region I their polarizationleads to an unfavorable interaction between regions I1 and 111. This “penalty” for extra polarization represents the energy involved in breaking the strucure of the bulk solvent upon solvation of polar systems. “Handbook of Biochemistry”, H. A. Sober, Ed., The Chemical Rubber Co., Cleveland, Ohio, 1968. M. Eiegen and L. De Maeyer in ”Technique of Organic Chemistry”, Voi. 111, ”Rate and Mechanisms of Reactions”, Part 11, S. L. Frless, E.S. Lewis, and A. Weissberger, Ed., Wiley-Interscience, New York, 1963, pp 1031. F. H. Westhelmer and M. W. Shookhoff, J. Am. Cbem. Soc., 61, 555 (1939). C. Tanford, J. Am. Cbem. Soc., 79, 5348 (1957). H. L. Friedman and C. V. Krishnan in “Water-A Comprehensive Treatise”, Vol. 3, F. Franks, Ed., Plenum Press, New York. L. G. Hepler, J. Phys. Chem., 61, 1426 (1957). M. J. S. Dewar and E. Haseibach, J. Am. Chem. Soc.,92,590 (1970). J. E. Bartmess and R. T. McIver, Jr., in “Gas-Phase Ion Chemlstry”, M. T. Bowers, Ed., Academic Press, New York, 1979. P. Kebarle, Annu. Rev. Pbys. Chem., 28, 495 (1977). (a) V. Renugopalakrishnan, A. V. Lokshminarayanan, and V. Sasisekheran, Biopolymers, 110, 1159 (1969); (b) w. D. S.Motherweli and N. W. Isaacs, J. Mol. Biol., 71, 231 (1972); (c) S. D. Steiiman, B. Hingerty, S. B. Broyde, E. Subromanlen, T. Sato, and R. Langridge, Biopolymers, 12, 2731 (1973). (a) G. N. Ramachandran and V. Sasisekharan, Adv. Protein Chem., 23, 283 (1968); (b) V. Sasisekharan and V. Lokshminarayanan, Biopolymers, 8, 505 (1969). A. Warshei and N. Kaiienbach, to be published. (a) A. Warshel and M. Karplus, J. Am. C k m . Soc., 94, 5612 (1972); (b) A. Warshei in “Modern Theoretical Chemistry”, Voi. 7, G. Segai, Ed., Plenum Press, New York, 1977. (c) A Warshei and M. Levltt, Quantum Chemistry Program Exchange No. 247, Indlana University. E. G. McRae, J. Chem. Pbys., 61, 562 (1957). A. T. Amos and B. L. Burrows, Adv. Quantum Cbem., 1, 289 (1973). F. J. Kamps, Cbem. Phys. Left., 26, 334 (1974). A. Warshel, Proc. Natl. Acad. Sci. U.S.A., 75, 2558 (1978). A. Warshei, Nature (London), 289, 679 (1976). (a) A. Walch and L. L. Ingraham, Arch. Biocbem. Biopbys., 156, 261 (1973); (b) S. J. Milder and D.S. Kliger, Pbotochem. Photobbl., 25, 287 (1977); (c) B. Honig, A. D. Greenberg, U. Dlnur, and T. 0. Ebrey, Biochemistry, in press. A. G. Soos, Annu. Rev. Pbys. Cbem., 25, 121 (1974). (a) A. Warshei and S. Lifson, J. Chem. Pbys., 53, 8582 (1970); (b) E. Huler arid A. Warshel, Acta Crystallogr., Sect. B, 30, 1822 (1974). See, for example, J. E. Leffeler and E. Grunwald, “Rate and Equilibria of Organic Reactions”, Wiiey, New York, 1963. A. Warshel, Proc. Natl. Acad. Sci. U.S.A., 75, 5250 (1978). D. Greenberg and A. Warshel, to be published. R. B. Hermann, J . Cbem. Phys., 70, 2754 (1972). J. A. Reymolds, D. B. Gilbert, and C. Tanford, Proc. Natl. Acad. Sci. U.S.A., 71, 2925 (1974). C. A. Coulson and U. Danielsson, Ark. Fys., 8, 245 (1954).

Calculations of chemical processes in solutions - ACS Publications

Recommend Documents