J. Phys. Chem. 1993,97, 11901-1 1907
11901
Multicavity Reaction Field Method for the Solvent Effect Description in Flexible Molecular Systems Mati Karelson,*-tToomas Tamm,t and Michael C. Zemer*** Department of Chemistry, University of Tartu. 2 Jakobi, St., EE2400, Tartu, Estonia, and Department of Chemistry and Quantum Theory Project, University of Florida, Gainesville, Florida 3261 1 Received: June 10, 1993"
A multicavity reaction field method is proposed for the quantum-chemical calculation of medium dielectric effects on the electronic structure of molecules in solution. This method is based on the gauge-independent partitioning of the molecule total reaction field in the polarizable medium into the partial reaction fields which belong to the rotationally or inversionally flexible groups in the molecule. The results of the calculations of the solvation free energies in aqueous solution of molecules containing such groups are presented and compared with the results obtained using the single spherical cavity self-consistent reaction field (SCRF)method.
1. Introduction
The general application of modern quantum-chemical methods to many applied problems of chemical technology, to molecular design of new biologically and pharmacologically active compounds, and to a prediction of physical and spectroscopic properties of chemical substances in solutions depends much upon the availability of a simple and efficient model of solvation applicable within the framework of modem quantum theory. Whereas classical and semiclassical molecular dynamics and Monte Carlo simulationsgive valuable information about the solvationprocess, especially in free energy terms, the need in computational time is generally prohibitively high even for the medium-size chemical systems using these methods (particularly if one needs to account for the solute-solvent electronic relaxation) and still resides in the hands of specialists. Therefore, the search for other less expensive theoretical methods that can easily be used and that account for the description of solvent effects remains very important. The interaction between a solute molecule and a solvent is usually broken down into two parts, the longer term interactions with energies comparable to or greater than kT for the solution and the transient interactions where the energies are less than kT. The stronger interactions are assumed to involve formation of a complex with a rather fixed geometrical relationship between the solute molecule and one or more solvent molecules (specific or microscopic solvation). Hydrogen bonding is a good example of this type of interaction. More complex cases relate to the formation of charge-transfer complexes, molecular ion pairs, and aggregates of polyelectrolytes. The intermolecular relationships are usually well understood in these systems, and the geometry of these complexes may be optimized by conventional quantumchemical gradient methods in the supermolecule comprising the solute molecule and electrostatically or hydrogen bonded addenda. The other type of interaction is due to the local orientational and electronic polarization of the solvent in the electric field of the solute molecule (nonspecific or macroscopic solvation). The statistical distribution over the relative orientations of the molecules controls the strength of the interaction, but for energies less than kT, the geometrical relationship is transient. A theoretical treatment of such interactions may proceed by using the macroscopic dielectric properties of the solvent to represent the observable polarizability of a randomly oriented assembly of solvent molecules in the vicinity of solute molecule as discussed To whom correspondence should be addressed. University of Tartu. t University of Florida. Abstract published in Aduance ACS Abstracrs, October 15, 1993. t
0022-365419312097- 11901$04.00/0
below.' A proper quantum-chemical model should be able to account for both interactions; thus, the macroscopic solvent polarity and polarizability (dielectric) effects must be introduced into the molecular Hamiltonian in the calculation of the solutesolvent supermolecule complex. A simple quantum-chemical reaction field model for the introduction of solvent dielectric effects into molecular Hamiltonians has been based on the classical Kirkwood4nsager theory of the electrical polarization in liquids.*-3According to this theory, the solvent is described as a polarizable dielectric continuum, characterized by its static dielectric constant e. The solute molecule is represented by Mpoint charges (ei), situated at fixed points (ri) inside of a sphere of radius a,. By solving the appropriate electrostatic equations inside and outside the sphere and applying the boundary conditions that are needed, one can obtain the following equation for the total electrostatic energy due to the molecular charge distribution p(r) interacting with the polarized medium:
In this equation, the summation extends over all charged particles in the solute molecule, and P/(cos0,) is the Legendre polynomial of 1-th order. In a quantum chemical model, this interaction energy can be represented by the respective additional terms in the molecular Hamiltonian which can be divided into nuclear-nuclear, nuclearelectronic, and electron-lectron terms.3 For most applications, (ionic and dipolar molecules) the two first terms in the infinite expansion of eq 1 are most important, corresponding to the Born ionic solvation (Vo) and the dipolar solvation terms ( V I ) . Therefore, the following modified Hamiltonian H'(1) can be used in the Schriidinger equation for a polar solute molecule in a polarizable medium:
H' = H,
+ v, + v,= H , + roQZ+ rr2
(2) where Ho is the Hamiltonian for the isolated molecule, Q the net charge, and I.( the dipole moment operator. The multipliers, P and r, in the perturbational term (the reaction field tensors) are functions of both the dielectric properties of the solvent and the size and shape of the solute molecule. In general, every atomic nucleus and electron in the molecule may have its own P and r value (cf. general Born solvation theory4-10). In the framework of the classical Kirkwood-Onsager theory, corresponding to a spherical cavity with radius a, into which the solute molecule is 0 1993 American Chemical Society
11902 The Journal of Physical Chemistry, Vol. 97, No. 46, 1993
embedded, common values of these variables are assumed for every charged particle in the solute molecule, and from eq 1
Karelson et al.
SCHEME I (dpded A-dd d A)
(charge dcrt d A)
r=
- 1 1 / ( 2 ~+ i ) ~ , 3
The use of the macroscopic value of the dielectric constant of the solvent in this formula is justified only in the case of the timeaveraged orientational polarization of the solvent molecules in the field of the solute molecule, i.e., systems in equilibrium with the solvent. If the physical or chemical phenomenon investigated by this method is characterized by a lifetime shorter than the orientational relaxation time of the solvent, an alternate value of e (and I') should be used. In the limit of the purely electronic polarization of the solvent, the dielectric constant e is often related to the square of the refraction index nD2 of the solvent in eqs 3a and 3b. The electronic energy of the solute molecule in the dielectric medium, Eel, may be found then by the solution of the respective one-electron Fock equations using the self-consistent reaction field (SCRF) procedure. Two principally different versions, based on different definitions of the total energy variational functional, are applicable for proceeding (described as model A and model B in ref 11, respectively) and lead to the corresponding sets of Hartree-Fock reaction field (HFRF) equations:
wherefo is the Fock operator for the unperturbed solute molecule; J, denotes the total molecular wave function (i.e., represents the expectation value of the total dipole moment of molecule); and ei and r#Ji are the one-electron eigenvalues and eigenvectors (orbitals), respectively. In model A,I1-l5the resulting quantum-chemical energy does not include the energy cost of solvent polarization, and thus the respective correction term E', = I'2 has to be added to obtained total energy. In model B,11JG19 the SCRF total energy includes this term explicitly. The two models often give close values of total energies and molecular geometries, but the oneelectron properties may vary significantly.20 Use of the reaction field in quantum-chemical theory requires that the shape and volume of the solute molecule be defined uniquely for any set of compounds. Various approaches to calculate these are known in which presumably more realistic shapes of molecules are introduced into the reaction field approach.21 These include the ellipsoidal shape of the cavity for the solute m o l e c ~ l e ~or ~ -the * ~use of the "polarizable continuum model" (PCM).25-33 The latter proceeds from the division of a molecule into atomic overlapping spheres, each of which has the dielectric polarization on the surface of sphere due to the electrostatic potential on the surface created by the charge distribution on that particular atom and other atoms in the solute molecule. The respective electrostatic reaction field is found by a direct numerical integration over the surface finite elements and is added as a perturbation to the molecular Hamiltonian. Calculation of the electrostatic potential over the surface is complex, and approximations that relate this only toatomiccharge terms do seem satisfa~tory.~5>3~ Nevertheless, in spite of the more *realistic" shape of the molecular cavity in the polarizable medium, this model has been shown to be very sensitive to the choice of the basis set in quantum-chemical calculation and has failed in several cases, using a judiciously chosen cavity redius, where the single spherical cavity SCRF still gives correct solvational e n e r g i e ~ . ~It~has - ~ ~been shown by us3743 and others44-49 that this latter simple SCRF is satisfactory in the case of comparatively
(dpde d 8-dd d A) (charge d B-cd d A)
(charge d B-drf d A) (dpde d A-dd d B)
rigid molecules with large dipole moments (aromatic heterocycles, e.g.), whereas it fails in the case of flexible aliphatic-chain molecules with freely or weakly restricted rotating groups along the bonds of molecules with a net dipole moment of zero or close to zero. Of course, the dipole model described above will fail in those cases in which the dipole moment is small or zero and higher moments dominate the total interaction. These higher moments can easily be included, a t least in principle, in eq 1 and subsequent equations, but this has not been done in the present work. 2. Theoretical Framework
In the present work, we suggest an extension of the selfconsistent reaction field approach with the partitioning of a molecule into fragments having independent but interacting reaction fields. A solute molecule might be partitioned as a single fragment, regenerating the simple SCRF model itself, with two or more fragments based on physical grounds, or each individual atom might be the center of a fragment generating a method that resembles the PCM or various virtual charge models.50151 However, the basic idea of the present partitioning scheme is the coupling of the rotation and inversion in the solute molecule with the solvent orientation and electronic relaxation. Therefore, only these groups in the molecule which are rotationally free or have small rotation or inversion barriers (flexible groups) are considered to interact independently with the dielectric medium. Each of these fragments can be characterized by its partial charge q A and partial dipole moment MA (and higher electrical moments) and is embedded into a spherical or ellipsoidal cavity surrounded by a dielectric continuum characterized by its macroscopic dielectric permittivity e ( w ) . The latter is a function of the electric field frequency w , and, strictly speaking, this dependence should also be accounted for in the coupling of the rotation-inversion flexibility found in the solute with the orientational and librational m o v e m e n t ~ ~ ~of- 5the ~ solvent molecule(s). For processes that have time scales sufficiently separated from the time scale of the solvent dielectric (orientational) relaxation, one may restrict the calculation with two values oft, corresponding to the full solvent relaxation (e, = e(w=O)) and the solvent electronic relaxation (e= e ( w = m ) ) . In the case of ellipsoidal cavities, a single value of the reaction field multiplier (3) is changed by the respective (diagonal) tensor. We proceed by partitioning the solute molecule into M fragments, indexed as A, B, ..., which are embedded into M spherical cavities with radii aOi( i = A, B, ...). Assuming the
Solvent Effect in Flexible Molecular Systems
The Journal of Physical Chemistry, Vol. 97, No. 46, 1993 11903
Kirkwood-Onsager expansion (1) for each of these fragments and using the terms up to the dipole56 interaction, one obtains the representation of the total energy of solute in dielectric medium with the macroscopic dielectric permittivity toshown in Scheme I. In Scheme I,
(7) and
are the reaction field tensor elements for the charge and dipole interactions, respectively; E,denotes the energy of the electrostatic interaction between the charges in cavities, and RABis the distance between the centers of fragments A and B. For these centers, one might use the center of mass of each fragment or the center of charge of each fragment. In the case of off-cavity-center electrical moments on fragments, the reaction field image moment is still localized in the center of the cavity.'' Therefore, if the distance between the mass or charge centers of two fragments is smaller than the sum of the respective cavity radii, the RABhas to be calculated as the distance between the charge or mass center of one fragment and the cavity center of the other. We notice, however, that the vector connecting neighboring fragments A and B will usually not pass through the dielectric medium. Therefore, the multicavity representation of the reaction field, whereas rigorous mathematically, will contain physical simplifications. Nevertheless, due to the presence of the interaction terms between thechargedistribution of one fragment and reaction field of the other, the dielectric representation of the solvent is not homogeneous and contains perturbations by the different fragments surrounding the charge distribution of a given fragment. A somewhat similar problem exists in the PCM model:25?26the electrostatic potential on the surface of atomic-centered spheres is generated from fields that sometimes pass through thedielectric medium and sometimes do not. The potential is, however, calculated only from fields that are assumed not to pass through the medium. In this work, the field calculated on one fragment is generated from the fields of other fragments which pass through the dielectric medium. The terms in eq 6 (Scheme I) correspond to the energies arising from the interaction of the partial charge and the dipole moment of a molecular fragment with the dipole reaction field (drf) and the charge reaction field (crf) of its own and other fragments as well from the interaction between the reaction fields of different fragments. The interactions between the charge and drf, dipole and crf, and crf and drf of a given fragment are zero because of the form of expansion (1). The first two terms in eq 6 can be referred to as first-order reaction field energies, the next four terms as second-order energies, and the last three terms as thirdorder energies. We also notice that the use of the dielectric permittivity e, in eq 6 assumes full orientational and electronic relaxation of the solvent in the field of solute fragment and therefore corresponds to the full coupling between rotationalinversional movements of the solute fragment and the solvent molecules. The factor FA in eq 7 accounts for the gauge independency of the dipole moment of the fragment and is defined as
(9) in the electronic reaction field tensor for fragments with a total nuclear charge ZAand with a negative net charge q A (FA= 1 in
the nuclear reaction field tensor in that case) or r)
in the nuclear reaction field tensor for fragments with a positive net charge q A (FA = 1 in the electronic reaction field tensory in that case). The quantum-mechanical extension of eq ( 6 )has the following form:
where H, is the solvent-unperturbed Hamiltonian of the solute molecule, p~ and p~represent the dipole moment operators of the respective molecular fragments, and PAand PBare the projection operators that yield the partial charge on fragments A and B. The PA and PB operators can be defined following one of the possible LCAO MO charge-partitioning schemes in the molecule (Mulliken, ZDO, etc.), and different schemes may be used in the framework of the present theory. For a totally consistent theory, the total charge of the fragments should equal the total charge of the molecule, and the total dipole should equal the (vector) sum of the fragment dipoles. The former is easy to ascertain, as all schemes for charge partitioning ascribe electrons to nuclei, and it is unlikely that dividing a molecule into fragments will include a single atomic nucleus to both fragments. The same is not true for the dipole or higher moments, which include (in addition to one-center terms) bond terms. In most ZDO schemes (as CNDO, INDO, AM1, etc.), these two-center terms are not considered. In ab initio methods, and in spectroscopic INDO,5*?59 these two-center terms are included. In such a case, the bond contribution that lies between fragments must be divided between the two fragments (or treated as a separate "bond fragment") in such a way that the sum of the fragment moments equals the total solute moment. By reordering the terms and defining (for brevity) new reaction field factors as r A ( l ) = rA/RAB, rA(') = rA/RAB*, = FA/ RAB',FA0(')= rAo/RAB, and r A d 2 ) = rAo/RAB2,One iS led to a more compact representation of eq 11:
Karelson et al.
11904 The Journal of Physical Chemistry, Vol. 97, No. 46, 1993
A#B PA(
= (6$1H,
($Pel$) - JW)+ C.C.
+ H , - wl$)
+c.c. = 0
(14)
where HS is the solvation Hamiltonian as defined in the last equation and c.c represents the Hermitian conjugate. This generates a many-electron SchrMinger equation
( H , + H , ) 9 = W$
+
6G = ( S $ l H l $ ) = (S$IH,, Hs - 1/2H, -El$) = 0 (22) and the one-electron Fock operator (cf. eq 17) is given as
(15)
with W,in this case (see below, method A), playing the role of the solute energy plus the solute-solvent interaction energy, not including the changes in solvent energy. By separating the nuclear and electronic variables, the following Hartree-Fock-type equations can be derived from the variation (14):
f& = ti$i
gives a correction for the electrostatic terms, doubly accounted for in model A. However, one can define a variational functional under the condition that the resulting total energy of the solute molecule would contain exactly the solvation energy (model B). In that case, the variational condition has the following form (cf. eq 13):
wheref, is again the Fock operator for the solvent-unperturbed solute molecule. The terms for the nuclear part of the fragment dipole moment
i = 1, n
(16) where 9, are the molecular orbitals in the n-electron system and
A
mEA
+ + (rA(Z)+ r;(2) + rAr;(2))($lpB~$)])p
($lpBI$)
(17)
In the SCRF procedure, ($Ip~l$)= FA and ($IP&) = q B are, respectively, the dipole moment and the partial charge on the corresponding fragment A or B, calculated from the charge distribution of the previous self-consistent-field cycte. Analogous terms for the nuclear part of the fragment dipole moment are accounted for classically as follows: A
+ r ~ r ~ "( J/IPB~$)])Z,R, ~)) (24)
and the terms corresponding to the interactions considering in eq 19 for model A
+
E,, = c 2 r o A q 2 A
mfA
where Z, denotes the nuclear (core) charge of the nucleus m with the radius vector R,. Also, the interaction of the net charge of each fragment with its own reaction field image charge, the interaction between the total dipole moment and net charge of each fragment, and the reaction field image dipoles and charges on other fragments, as well as the interaction between different reaction field image charges and dipoles, are treated according to classical electrostatics:
Therefore, the total energy E,,, of the solute-polarizable medium system is represented by the following formula:
(20) = E(qc) + ENUC + ECL+ Ecorrcction where E(qc) is the electronic energy of the molecule in the Etot
3. Results and Discussion
In order to test the model proposed in this work, we calculated the energies of transfer from the gas phase (g) to the aqueous solution (w) as the difference between the calculated heats of formation of compounds in these two media: AEaq= AHO,(w) - AH0Xg) (26) using AMIm and PM361 single-cavity (SCa) and multicavity (MCa) SCRF Hamiltonians for a set of small organic compounds.63 The respective experimental data, AG(g+w), were obtained from the compilation by Hine and Mookerjee.64 The list of compounds and their molecular weights, densities, and cavity radii in aqueous solution as estimated from the mass densities are given in Table I. The cavity radii for fragments were calculated starting from the total volume of the molecule and from the volume for the methyl group, which was obtained as one-half of the molecular volume of ethane. The fragment cavity radii for compounds are also presented in Table I. In Tables I1 and 111, the AM1 calculated heats of formation AHof in the gas phase and in the dielectric medium with e = 80 (aqueous solution), and the corresponding transfer energies AEaq,are given
Solvent Effect in Flexible Molecular Systems
The Journal of Physical Chemistry, Vol. 97, No. 46, 1993 11905
TABLE I: Calculated Total Cavity Radii and Fragment Cavity Radii for the Molecules of This Study CH3CH 30.07 0.572 2.752 2.184 34.03 0.8428 2.520 2.184 CH3F CH3C1 50.49 0.9159 2.796 2.184 94.94 1.6755 2.822 2.184 CH3Br CH,I 141.94 2.279 2.912 2.184 CH3OH 32.04 0.7914 2.523 2.184 CH3NH2 31.06 0.6628 2.649 2.184 CH3COOH 60.05 1.0492 2.831 2.184
(CH3) (CHI) (CH3) (CHI) (CH3) (CHI) (CH3) (CH3)
2.184 (CHI) 1.775 (F) 2.253 (Cl) 2.292 (Br) 2.426 (I) 1.779 (OH) 2.013 (NH2) 1.87 (CO)
TABLE III: AM1 Si le-Cavity and Multicavity SCRF (Model B) Calculated fiats of Formation AIPf (kcal/mol) in the Gas Phase and in the Aqueous Solution (e = 80), Corresponding Energies of Transfer A&, and Experimental Free Energies of Transfer of Compounds from the Gas Phase to the Aqueous Solution A6lg-w)"
-
-
LY0r(c=80)
1.779 (OH)
CH3CN 41.05 0.7857 2.779 2.184 (CHI) 2.226 (CN) CH3CHN 44.05 0.7834 2.815 2.184 (CH3) 2.281 (CHO) CH3CONH2 59.07 0.9986 2.862 2.184 (CH3) 1.879 (CO) 2.013 CH3COCH3 CH$H CH3SCH3
-
58.08 0.7899 3.078 2.184 (CH3) 2.026 (CO) 48.11 0.8665 2.803 2.184 (CH3) 2.264 (SH) 62.13 0.8483 3.074 2.184 (CH3) 2.018 (S)
C H I C H ~ O H 46.07 0.7893 2.850 2.184 ICH?) . -, 1.922 (CH,) .
-I
(CH~OH)I (CH2NH2h
62.07 1.1088 2.810 1.757 (OH) 1.783 (CH2) 60.11 0.8995 2.981 1.892 (NH2) 1.863 (CH2)
(NH2)
2.184 (CH3) 1.779 (OH)
-
Molecular weight. Density.73e Cavity radius for a single cavity. Cavity radius for a fragment (given in parentheses).
TABLE Ik AM1 Single-Cavity and Multicavity SCRF (Model A) Calculated Heats of Formation AHDr (kcal/mol) in the Gas Phase and in the Aqueous Solution (c = 80), Corresponding Energies of Transfer A&, and Expenmental Free Energies of Transfer of Compounds from the Gas Phase to the Aqueous Solution AC(g+w)" Wr(c-80)
compd
W d g )
SCa
MCb
-17.41 -17.41 -17.51 -66.85 -61.02 -65.22 -21.23 -18.95 -22.10 -6.19 -9.25 -7.43 5.67 3.26 3.42 -67.07 -57.03 -61.31 -6.11 -13.85 -10.55 -102.99 -125.52 -124.42 19.28 6.28 13.14 -50.17 -41.56 -51.84 -73.05 -3 1.76 -49.36 -57.47 -49.10 -57.60 -3.61 -8.15 -8.73 -9.34 -1 1.82 -13.13 -62.63 -77.58 -66.43 -108.82 -1 19.72 -134.56 -1 1.01 -16.73 -29.18
SCo
MCb
0.00 -4.20 -3.15 -3.06 -2.41 -4.27 -4.44 -22.53 -13.00 -10.28 -41.11 -8.37 -4.54 -2.48 -3.80 -10.90 -5.72
-0.10 -5.83 -2.21 -1.24 -2.25 -10.04 -7.74 -21.43 -6.14 -8.61 -17.39 -8.50 -5.12 -3.79 -14.95 -25.84 -12.45
compd
W d g )
CH3CH3 CH3F CH3CI CH3Br CH31 CH3OH CH3NH2 CHsCOOH CH3CN CH3CHO CH3CONH2 CH3COCHj CH3SH CHjSCHs CH3CH2OH (CHIOH)* (CH2NH2)2
-17.41 -61.02 -18.95 -6.19 5.67 -57.03 -6.11 -102.99 19.28 -41.56 -3 1.76 -49.10 -3.61 -9.34 -62.63 -108.82 -11.01
a
-0.80 -0.90 -5.07 -4.56 -6.70 -3.89 -3.55 -9.70 -3.81 -1.26 -1.56 -4.96 -7.75 -9.88
Single cavity. Multicavity. for models A and B, re~pectively.6~In Tables IV and V are presented the same data as obtained by using the PM3 method. First of all, we notice that model A in most cases gives largely overestimated AM1- and PM3-calculated transfer energies both for the single cavity and multicavity representations. In some cases (CH31, CH~CONHZ),PM3 SCRF does not converge at all. The use of model A has recently been criticized by Wang and Ford,lg who point out that the effective Hamiltonian, eq 15, does not commute with the Hamiltonian that yields the total energy. The relevance of this observation, however, is not clear to us. We point out rather that the best model will be that which is consistent with a systematic choice of cavity radii. For our choice of radii, we observe that model A can, indeed, lead to unrealistic results. Therefore, we limit ourselves in further discussion to model B. In this case (model B), there is a good linear relationship between the MCa SCRF calculated and experimental transfer energies using the AM1 method (cf. also Figure 1). For these studies, we have broken down the solute molecule into several rotationally flexible fragments with individual reaction fields. Thus, one group of the molecules studied (CH3C1, C H J B ~CH31, ,
MCb
-17.41 -17.51 -66.06 -62.67 -19.79 -2 1.02 -6.98 -6.89 5.05 4.43 -58.29 -62.93 -10.77 -8.31 -104.26 -108.80 15.96 13.98 -45.22 -44.24 -31.97 -42.56 -53.48 -51.27 -5.52 -5.48 -1 1.45 -10.02 -66.63 -63.47 -1 10.38 -118.25 -20.29 -11.99
SCa
MCb
0.00 -0.25 -0.84
-0.10 -5.04 -2.07 -0.70 -1.24 -5.90 -4.66 -5.79 -5.30 -3.66 -10.80 -4.38 -1.91 -2.11
-0.n -0.62 -1.26 -2.20 -1.27 -3.32 -2.68 -0.21 -2.17 -1.87 -0.58 -0.84 -1.56 -0.98
-4.00 -7.87 -9.28
AG(g+w) 1.81 -0.22 -0.56
-0.80 -0.90 -5.07 -4.56 -6.70 -3.89 -3.55 -9.70 -3.81 -1.26 -1.56 -4.96 -7.75 -9.88
Single cavity. Multicavity.
TABLE IV: PM3 Single-Cavity and Multicavity SCRF (Model A) Calculated Heats of Formation AIPf (kcal/mol) in the Gas Phase and in the Aqueous Solution (c = 80), Corresponding Energies of Transfer A&, and Expenmental Free Energies of Transfer of Compounds from the Gas Phase to the Aqueous Solution AC(g+w)"
AG(g+w) 1.81 -0.22 -0.56
SCa
compd
Wr(g)
-17.63 -53.80 -14.68 -1.98 9.45 -51.88 -4.00 -102.00 23.29 -44.20 -49.33 -53.31 -5.53 -10.96 -55.66 -95.14 -5.01
LY0(c=80) SCa MCb
hE,p SCa
MCb
-17.63 -18.29 0.00 -0.66 -56.98 -58.60 -3.18 -4.80 -17.43 -16.18 -2.75 -1.50 -6.58 -4.28 -4.60 -2.30 7.72 n.c.c n.c. -1.73 -55.42 -6 1.29 -3.54 -9.41 -9.74 -7.81 -3.81 -5.74 -125.41 -128.75 -23.41 -26.75 7.17 14.50 -16.12 -8.79 -53.10 -52.83 -8.90 -8.63
n.c.
n.c.
-60.56 -10.28 -14.58 -60.10 -104.88 -10.53
-63.17 -10.26 -16.07 -72.19 -111.63 -17.22
n.c.
n.c.
-7.25 -9.86 -4.75 -4.73 -3.62 -5.11 -4.44 -16.53 -9.26 -16.49 -5.52 -12.21
AG(g+w) 1.81 -0.22 -0.56
-0.80 -0.90 -5.07 -4.56 -6.70 -3.89 -3.55 -9.70 -3.81 -1.26 -1.56 -4.96 -7.75 -9.88
Single cavity. Multicavity SCRF not converging. CHsOH, CH3NH2, CH3CN, CH3CH0, CH3SH) is divided into two fragments, one of which is the methyl group and the other involves the substituent attached toit. Another groupof molecules (CH&OOH, CH3CONH2, CH3COCH3, CH3SCH3, CH3CH2OH) has three fragments, where two terminal substituents (-CH3, -OH,-NH2) areconnected with a bridgegroup ( > C = O , - C H r , 4%). Ethylene glycol and ethylenediamine are divided into four fragments (two methylene groups and two terminal substituents, either -OH or -NH2). The results of the calculation indicate that methyl fluoride is better represented by a single cavity, and we used the respective value in the correlation (see later). Also, for acetonitrile, which has a rigid CN group connected to the methyl group, the single cavity representation gives a value closer to the experimentally observed transfer energy. However, for the rest of the molecules studied one has obtained remarkably good agreement between experimental and AM1 MCa SCRF energies. The linear regression of the calculated energies of transfer against the experimental values AG(g+w) gives the following relationship:
Karelson et al.
11906 The Journal of Physical Chemistry, Vol. 97, No. 46, 1993
TABLE V: PM3 Single-Cavity and Multicavity SCRF (Model B) Calculated Heats of Formation A I P t (kcalhol) in the Gas Phase and in the Aqueous Solution (c = 80), Corresponding Energies of Transfer A&,, and Experimental Free Energies of Transfer of Compounds from the Gas Phase to the Aqueous Solution A C ( ~ - W ) ~
-8
-6
-17.63 -17.63 CHpCH3 -54.77 -53.80 CH3F -15.40 CH3Cl -14.68 -2.92 CHpBr -1.98 8.67 CH31 9.45 -52.53 CHjOH -51.88 -5.98 CH3NHz -4.00 CHsCOOH -102.00 -106.04 19.21 CH3CN 23.29 CHoCHO -44.20 -46.56 CH3CONH2 -49.33 -54.02 CHaCOCH3 -53.31 -55.03 -6.89 CHBH -5.53 -12.00 CHoSCH3 -10.96 CH3CHzOH -55.66 -51.53 (CHZOH)~ -95.14 -98.26 -12.30 (CH2NHz)z -5.01 a Single cavity. Multicavity
-18.19 -57.73 -15.71 -4.09 7.79 -57.41 -7.33 -1 11.52 15.95 -47.87 -64.79 -57.34 -7.04 -13.42 -64.96 -105.14 -10.80
0.00 -0.56 -0.97 -3.93 -0.72 -1.03 -0.94 -2.11 -0.78 -1.66 -1.05 -5.53 -1.98 -3.33 -4.04 -9.52 -4.08 -7.34 -2.36 -3.57 -4.71 -15.46 -1.72 -4.03 -1.04 -1.51 -3.87 -2.46 -3.12 -9.30 -3.12 -10.00 -7.29 -5.79
1.81 -0-22 -0.56 -0.80 -0.90 -5.07 -4.56 -6.70 -3.89 -3.55 -9.70 -3.81 -1.26 -1.56 -4.96 -7.75 -9.88
-4
-
-2
-
lBj!j?$8Ob
10
0 -
2
17
15
2
I
I
I
11
I
I
AEaq -1 0
-8
-6
-4
-2
0 2
0
-2
-4
-6
-8
-10
AG(g-4 Figure 1. AM1 MCa SCRF calculated energies of transfer AE., (kcal/ mol) vs the experimental free energies of transfer AG(g+w) (kcal/mol) from the gas phase to the aqueous solution for the molecules of Table I. Model B is assumed in the calculations. The numbering of points corresponds to Table I; the data are those of Table 111.
+
AEaq= (0.33 f 0.30) (0.953 f 0.058)AC(g-w)
(27)
correlation coefficient R = 0.975 standard deviation u = 0.71 (kcal/mol) The slope of this relationship is close to unity, and the intercept is a small positive number. One must certainly not forget the other terms required in the total solvation energy, namely, the cavity formation, dispersion, and solvent restructurization effects. An effective method for a semiempirical treatment of these effects in combination with the general Born model for electrostatic interactions has been proposed by Cramer and Truhlar, and the results obtained by them for a wide range of compounds are very promising.6668 Notably, they also used the AM1 Hamiltonian to calculate the quantum-chemical part of the energy. However, the estimate of the absolute values of the cavitation and dispersion interactions is only a few kcal/mol for the set of compounds studied in this work, and, in addition, these terms are of opposite
sign and tend to ~ance1.6’~~The solvent restructurization term, which is still smaller, is also in part accounted for in the MCa SCRF model and is expected to be reasonably constant for this series of molecules all dissolved in water. Therefore, the existence of the relationship shown in eq 27 is not unexpected if the solvation model is good. However, in the case of a single-cavity model, which uses only the total dipole moment of the molecule, the corresponding relationship is quite poor, as shown in Figure 2, with the exceptions of those cases in which a single cavity is appropriate (Le., CHSF, CH3C1, CH3Br, CH31, and even CH3CN). In these cases, there is no real internal rotation of one group with respect to another. Consequently, the correct local representation of the solute molecule charge distribution and the respective reaction field lead to a more adequate description of the solvation interactions, particularly in the case of molecules with rotational flexibility. The results of the calculations using the PM3 parametrization are of considerably lower quality in comparison with the AM1 results, as shown in Figure 3. This is obviously due to the less satisfactory form of the corresponding semiempirical wave function for the description of the local molecular properties and charge distribution. Hence, one must be very careful in choosing the quantum-chemical method used for the calculation of solvation energies as this, in general, requires an accurate local representation of the charge distribution, not only a global one. In a b initiocalculations, particular attention must be paid to the choice of the basis set and to the method of accounting for electron correlation effects. The results of the present semiempirical treatment demonstrate, however, the possible usefulness of the MCa SCRF approach for testing the quality of different basis sets and various post-Hartree-Fockapproaches for the description of bulkintermolecular interactions in condensed media (solutions). More importantly, the practical applicability of the method present in this work would cover the study of the reaction paths in solutions-a major topic of contemporary chemistry. This work divides a molecule into several spherical regions in which the charge and dipolar terms are developed to represent the molecular charge density. This can be clearly extended to involve higher moments in each region, giving a more accurate representation of the electrostatic potential the solvent realizes. However, by breaking the molecule into smaller regions we not only better represent the internal rotational flexibility of the molecule, but we also need fewer moments to describe the charge
Solvent Effect in Flexible Molecular Systems
The Journal of Physical Chemistry, Vol. 97, No. 46, 1993
*%q
n
-10
0 15
-a -6
-4
-2
0
2
0
-2
-4
-6
-8
-10
WS-W Figure 3. PM3 MCa S C R F calculated energies of transfer AEaq(kcal/ mol) vs the experimental free energies of transfer AG(g-w) (kcal/mol) from the gas phase to the aqueous solution for the molecules of Table I. Model B is assumed in the calculations. The numbering of points corresponds to Table I; the data are those of Table V.
density to a given accuracy than would be required by a singlecavity treatment. The essential idea behind the MCa SCRF model presented here is that rotationally free fragments of the molecule each create their own reaction field. Thus, a molecule like CH3-OH, with free rotation about the C - O bond, requires two fragments, CH30-CH3, three fragments, etc. Molecules such as CH3F or CH3C N require only one fragment. The evidence for this is rather remarkable from our results. This is different in principle from other models that divide a molecule into overlapping (atomic) spheres to better represent the shape and electrostatic potential of a given solute molecule. Although, as mentioned above, our multicavity model also better represents the shape and potential, in general, than a single cavity model, it accounts in an average way for some of the solute internal dynamics.
Acknowledgment. Authors from the University of Tartu (M.K. and T.T.) thank the Estonian National Foundation of Pure and Applied Chemistry for financial support of this work. This work has been supported in part through a grant from the Office of Naval Research. References and Notes (1) Onsager, L. J . Am. Chem. Soc. 1936,58,1486. (2) BBttcher, C. J. F.; Bordewijk, P. Theory of Electric Polarization, 2nd ed.; Elsevier: Amsterdam, 1978;Vol. 11. (3) Hylton, J.; Christoffersen, R. E.; Hall, G. G. Chem. Phys. Lett. 1974, 26, 501. (4) Hoijtink, G.J.; de Boer, E.; van Meji, P. H.; Weijland, W. P. Recl. Trav. Chim. Pays-Bas. 1956,75,487. (5) Chalvet, 0.;Jano, I. C. R . S6anc6s Acad. Sci. 1964,259, 1867. (6) Fischer-Hjalmars, I.; Hendriksson-Enflo, A,; Hermann, C. Chem. Phys. 1977,24, 167. (7) Constanciel, R.;Tapia, 0. Theor. Chim. Acta 1980, 55, 177. (8) Constanciel, R.; Contreras, R. Theor. Chim. Acta 1984,65, 1. (9) Tucker, S. C.; Truhlar, D. G. Chem. Phys. Lett. 1989,157, 164. (10) Angyin, J. G. J. Math. Chem. 1992, 10, 93. (11) Karelson, M. M.; Zerner, M. C. J. Phys. Chem. 1992,92,6949. (12) Tapia, 0.; Goscinski, 0. Mol. Phys. 1976, 29, 1653. (13) KarlstrBm, G. J. Phys. Chem. 1987,92, 1315. (14) Pappalardo, R. R. J. Chem. Res. 1989,290. (15) Wong, M. W.; Frisch, M. J.; Wiberg, K. B. J . Am. Chem. SOC.1990, 112,4776. (16) Hylton McCreery, J.; Christoffersen, R. E.; Hall, G. G. J . Am. Chem. Soc. 1976,98,7191. (17) Rivail, J. L.; Rinaldi, D. Chem. Phys. 1976,18, 233. (18) Karelson, M. M.; Katritzky, A. R.; Zerner, M. C. Znt. J . Quantum. Chem. 1986,S20, 521.
11907
(19) Wang, B.; Ford, G. P. J. Chem. Phys. 1992,97,4162. (20) The analysis of the differences in two models as compared to experimental data will be the subject of a separate communication (Tamm, T.; Karelson, M., in preparation). See also ref 11. (21) Tapia, 0.In Molecular Interactions; Ratajczak, H., Orville-lhomas, W. J., Eds.; J. Wiley & Sons: Chichester, 1983;Vol. 3, p 183. (22) Rivail, J.-L.; Terryn, B. J. Chim Phys. 1982, 79, 1. (23) Rivail, J.-L.; Terryn, B.; Rinaldi, D.; Ruiz-Lopez, M. F. J. Mol. Struct. (THEOCHEM) 1985,120,387. (24) Ford, G. P.; Wang, B. J. Comput. Chem. 1992,13,229. (25) Mietrus, S.;Scrocco, E.; Tomasi, J. Chem. Phys. 1981,55, 117. (26) Hoshi, H.; Sakurai, M.; Yoshio, Y.; Chujo, R. J. Chem. Phys. 1987, 87, 1107. (27) Drummond, M. L. J. Chem. Phys. 1987,88,5014. (28)Drummond, M. L. J. Chem. Phys. 1987.88, 5021. (29) Frecer,V.;Majekova,M.;Miertus,S.J. Mol.Struct. (THEOCHEM) 1989,183,403. (30) Fox, T.; RBtsch, N. J. Mol. Struct. (THEOCHEM) 1992,276,279. (31) Negre, M.; Oroszko, M.; Luque, F. J. Chem. Phys. Left. 1992,196, 27. (32) Olivares del Valle, F. J.; Bonaccorsi, R.; Cammi, R.; Tomasi, J. J. Mol. Struct. (THEOCHEM 1991. 76. 295. (33) Olivares del Valle, F. J.; Bonakrsi, R.; Cammi, R.; Tomasi, J. J. Mol. Struct. (THEOCHEM) 1991,78,349. (34) Gould, I. R.; Green, D. V. S.; Young, P.; Hillier, I. H. J . Org. Chem. 1992.57.4434. (35) Parchment, 0.G.; Hillier, I. H.; Green, D. V. S.;Burton, N. A.; Morley, J. 0.;Schaefer,H. F., I11 J . Chem.Soc., Perkin Trans. 2,1992,1682. (36) Woodcock, S.;Green, D. V. S.; Vincent, M. A.; Hillier, I. H.; Guest, M. F.; Sherwood, P. J. Chem. Soc., Perkin Trans. 2 1992,2151. (37) Tamm, T.; Karelson, M. Organic Reactivity (Tartu) 1989,26,211. (38) Karelson, M. M.; Katritzky, A. R.; Szafran, M.; Zerner, M. C. J . Org. Chem. 1989,61,6030. (39) Karelson, M. M.; Katritzky, A. R.; Szafran, M.; Zerner, M. C. J . Chem. SOC.,Perkin Trans. 2 1990,161. (40) Karelson, M. M.; Tamm, T.; Katritzky, A. R.; Szafran, M.; Zerner, M. C. Int. J. Quantum Chem. 1990,37, 1. (41) Karelson, M.; Zerner, M. C. J. Am. Chem. Soc. 1990,112,9405. (42) Katritzky, A. R.; Karelson, M. J. Am. Chem. SOC.1991,113,1561. (43) Tamm, T.; Karelson, M. Proc. Est. Acad. Sci. Chem. 1991,40,122. (44) Agren, H.; Medina-Llanos, C.; Mikkelsen, K. V. Chem. Phys. 1987, 115,43. (45) Mikkelsen,K. V.; Agren, H.; Jensen, H. J. Aa.; Helgaker, T. J . Chem. Phys. 1988,89,3086. (46) Agren, H.; Mikkelsen, K. V. J . Mol. Struct. 1991. 80. 425. (47) Agren, H.; Knuts, S.; Mikkelsen, K. V.; Jensen, H. J: Aa. Chem. Phys. 1992,159,211. (48) Wong, M. W.; Frisch, M. J.; Wiberg, K. B. J . Am. Chem. SOC.1991, 9 A716 -I -I -, . . . -.
(49) Wong, M. W.; Wiberg, K.B.; Frisch, M. J . Chem. Phys. 1991,95, 8991. (50) Noell, J. 0.; Morokuma, K. J . Phys. Chem. 1976,80,2675. (51) Constanciel, R.; Tapia, 0. Theor. Chim. Acta 1978,48, 75. (52) Su,S.-G.; Simon, J. D. J . Phys. Chem. 1987, 91,2693. (53) Castner, E. W., Jr.; Maroncelli, M.; Fleming, G. R. J . Chem. Phys. 1987,86,1090. (54) Castner, E. W.,Jr.; Fleming, G. R.; Bagchi, B. Chem. Phys. Lett. 1988,143,270. (55) Bagchi, B.; Chandra, A. J. Chem. Phys. 1990,93, 1956. (56) The extension of the theory for using higher terms in the expansion shown in eq 1 as well as the use of ellipsoidal cavities instead of spherical is straightforward and will omitted here. From the technical point of view, this may not be trivial and will be discussed in further communications. (57) BBttcher, C. J. F. Theory ofElectric Polarization, 2nd ed.; Elsevier: Amsterdam, 1973;Vol. 1, Chapter IV. (58) Ridley, J. E.; Zerner, M. C. Theor. Chim. Acta 1973,32, 111. (59) Bacon, A. D.; Zerner, M. C. Theor. Chim. Acta 1979,53,21. (60) Dewar, M. J. S.; Zoebisch, E. G.; Healy, E. F.; Stewart, J. J. P. 1985, 107,3902. (61) Stewart, J. J. P. J . Comput. Chem. 1989, 10, 209. (62) Stewart, J. J. P. MOPAC 6.0. QCPE 1989, No. 455. (63) MOPAC version 6.062was modified with the inclusion of the MCa SCRF model discussed in the present work. (64) Hine, J.; Mookerjee, P. K. J . Org. Chem. 1975,40, 292. (65) The data by Hine and Mokerj& correspond to room temperature (25 “C), where thedielectric constant of water is 78.3.” However, the results of the calculations with e = 80 will be the same with the precision given in Tables 11-V. (66) Cramer, C. J.; Truhlar, D. G. J . Am. Chem. Soc. 1991,113, 8305, 8552. (67) Cramer, C. J.; Truhlar, D. G. J . Comput. Chem. 1992,13, 1089. (68) Cramer, C. J.; Truhlar, D. G. Science 1992, 256,213. (69) Pierotti, R. A. Chem. Rev. 1976,76,717. (70) Roesch, N.; Zerner, M. C., in preparation. (71) Karelson, M. M. Theor. Eksp. Khim. 1979,15, 80. (72) Karelson, M. M. Organic Reactivity (Tartu) 1984,21, 160. (73) CRC Handbook of Chemistry and Physics, 57th ed.; CRC Press: Boca Raton, FL, 1976.
