Model for Aqueous Solvation Based on Class IV Atomic Charges and

J. Phys. Chem. 1996, 100, 16385-16398

16385

Model for Aqueous Solvation Based on Class IV Atomic Charges and First Solvation Shell Effects Candee C. Chambers,† Gregory D. Hawkins, Christopher J. Cramer,* and Donald G. Truhlar* Department of Chemistry, Supercomputer Institute, and Army High Performance Computing Research Center, 207 Pleasant Street SE, UniVersity of Minnesota, Minneapolis, Minnesota 55455-0431 ReceiVed: April 10, 1996; In Final Form: May 31, 1996X

We present a new set of geometry-based functional forms for parametrizing effective Coulomb radii and atomic surface tensions of organic solutes in water. In particular, the radii and surface tensions depend in some cases on distances to nearby atoms. Combining the surface tensions with electrostatic effects included in a Fock operator by the generalized Born model enables one to calculate free energies of solvation, and experimental free energies of solvation are used to parametrize the theory for water. Atomic charges are obtained by both the AM1-CM1A and PM3-CM1P class IV charge models, which yield similar results, and hence the same radii and surface tensions are used with both charge models. We considered 215 neutral solutes containing H, C, N, O, F, S, Cl, Br, and I and encompassing a very wide variety of organic functional groups, and we obtained a mean unsigned error in the free energy of hydration of 0.50 kcal/mol using CM1A charges and 0.44 kcal/mol using CM1P charges. The predicted solvation energies for 12 cationic and 22 anionic solutes have mean unsigned deviations from experiment of 4.4 and 4.3 kcal/mol for models based on AM1 and PM3, respectively.

1. Introduction The standard-state free energy of solvation ∆GS° is a fundamental quantity characterizing the interaction of a solute molecule with a given solvent. In addition to its fundamental interest, ∆GS° may be combined with other thermodynamic data to predict a variety of equilibrium constants, the most important of which are solubility1 and the partitioning of a solute between immiscible phases.2 In previous work,3-11 an approach to calculating ∆GS° has been presented based on a two-step calculation: (1) an estimate of the electrostatic contribution based on treating the solvent as a dielectric continuum; (2) a semiempirical estimate of the additional contribution to the free energy from the first solvation shell. We believe that several aspects of such an approach are essential for a realistic treatment of the physical situation.12,13 (a) Both steps should treat the shape of the solute realistically (as opposed, for example, to assuming a spherical or ellipsoidal cavity). (b) In step 1, the description of the solute charge distribution should allow, at a minimum, for a partial charge at each atomic center (as opposed, for example, to including only the overall dipole moment or some other limited single-center multipole description). (c) Also in step 1, the solute should be polarizable.12-14 (d) Because the concepts of solute radius and of a boundary between the solute and the dielectric medium are intrinsically ill-defined,15 the parameters in the second step must be consistent with the treatment of the dielectric boundary assumed in the first step (for example, if the solute radius is made larger, the electrostatic contribution from the first solvation shell included in step 1 will be smaller, so the first solvation shell contributions that are included in step 2 must account for the difference). Step 2 must account not only for the deviation of the electrostatics in the first solvation shell from the continuum estimate of step 1 but also for other specific effects of which † Address as of August 1, 1996: Departments of Physics and Chemistry, Mercyhurst College, 501 East 38th St., Erie, PA 16546. X Abstract published in AdVance ACS Abstracts, September 15, 1996.

S0022-3654(96)01077-5 CCC: $12.00

the most prominent are solute-solvent dispersion interactions, solvent cavitation, nonelectrostatic components of hydrogenbonding, solute-induced changes in solvent structure and entropy, and so forth. It is reasonable to assume that such effects are proportional to the average number of solvent molecules in the first solvation shell, and the essence of our approach is to assume that this average number is proportional to the solventaccessible surface area.16,17 The proportionality constant (which is a microscopic surface tension18-22) depends on the nature of a particular part of the solute with which the solvation shell is in contact. For example, the surface tension of a fluorine atom at the surface of the solute is different from the surface tension of a chlorine atom. By estimating the surface tensions semiempirically, we compensate for the theoretical ambiguity of what values to use for the radii of the solute atoms in accounting for electrostatic effects. To the extent that our assumption about the nonelectrostatic contributions being proportional to solventaccessible surface area is correct, our parameters will be physical if the functional form of the surface tensions captures the chemical variation of first solvation shell interactions. In light of the above discussion, we emphasize that the precise division of physical effects into “electrostatics” and “first solvation shell effects” is ambiguous for at least three reasons. First, free energies of systems with interacting parts are never, in either theory or experiment, decomposable into contributions from individual parts; due to the nature of entropy itself, they are thermodynamic properties of the whole system.23 Second, the electrostatic contribution of the first solvation shell is partly contained in step 1 and partly contained in step 2, and there is no unique way to get around this because a discrete boundary with a discontinuity in the dielectric constant between the solute and the first solvation shell is an ill-defined concept. Third, due to the semiempirical nature of step 2, any systematic errors in step 1 (for example, due to the partial charge model of the solute or the algorithmic solution of the electrostatic problem) will be compensated to some extent by the parametrization of step 2, again spreading electrostatic effects over both steps. Having recognized these ambiguities, we stress two important © 1996 American Chemical Society

16386 J. Phys. Chem., Vol. 100, No. 40, 1996 points. First, the overall free energy of solvation need not be inaccurate due to the difficulty of separating effects of individual solvation shells, provided that the electrostatics and first solvation shell effects are treated consistently. Second, when interpreted with appropriate attention to the dependence on solute atomic radii, the separate evaluation of electrostatic and first solvation shell effects can be a source of considerable physical insight. Nevertheless, one should make every effort to make the radii as physical as possible because the radii modulate the influences of the solvent not only on the solute energetics but also on the solute wave function and hence on the change in molecular properties induced by the solvent. For example, chemical reactivity, molecular recognition, and biological activity are known to correlate with electrostatic potentials, and it may be critical to achieve a physical allocation of solvation effects between the electrostatic and surface tension terms in the effective Hamiltonian if one is to calculate such molecular properties as well as solvation energies. Both physical radii and realistic partial charges can contribute to achieving such a useful allocation. The reader is referred to more elaborate treatments for further attempts to identify individual contributions to solvent effects.24 Although it need not be an intrinsic choice for the approach described, we have chosen to solve the electrostatic problem by the generalized Born approximation25-27 (GBA) with solute descreening effects included by integrating the dielectric energy density over the region exterior to the solute.22,28-30 One obvious alternative would be a numerical solution of the Poisson equation, but we believe that the extra effort involved in such an approach is not warranted because of the intrinsic uncertainty in the location of the dielectric boundary. That is, since the electrostatic formulation of the physical situation is intrinsically approximate, the seeming rigor of a full numerical solution in step 1 need not improve the oVerall model. Having opted for the GBA, we gain an advantage in interpretation since the GBA permits a decomposition of electrostatic effects into contributions from individual atomic sites and pairs thereof. The purpose of the present paper is to give a new parametrization for aqueous solvation along the above lines. The new parametrization is called Solvation Model 5.4 or SM5.4, and it allows the treatment of a wide variety of organic (and some inorganic) solutes in water. Furthermore, it is designed to be extendable to other solvents (e.g., alcohols, ethers, hydrocarbons, and halocarbons) in later work. We note that the first part of the name, SM5, denotes what will eventually become a whole suite of models based on the functional form presented here for the surface tensions, and the second part of the name, the 4 after the point, denotes the use of class IV charges in the present paper. At this point it is useful to distinguish the SM5 approach from the earlier models, called SM1 through SM4. SM1,3 SM1a,3 SM2,4,6 SM2.1,7 SM2.2,31 SM3,5,6 and SM3.17 were parametrized only for water, and the partial charges for the electrostatics were based on Mulliken analysis32 of semiempirical molecular orbital wave functions, in particular of wave functions obtained by the popular Austin Model 133-35 (AM1) and Parameterized Model 336 (PM3). These charges, like those calculated from ab initio minimum-basis-set Hartree-Fock calculations, tend to underestimate charge separations and polarities. Since the errors are systematic, the surface tensions do a good job of compensating when the methods are applied to “typical” solutes. However, one’s confidence in predictions for transition states10 and other difficult cases would be justifiably higher with more accurate partial charges. Models SM1 through SM3.1 differ from one another in terms of the functional forms used for surface tensions and the

Chambers et al. numerical treatment of dielectric descreening; however, all of these models treat dielectric descreening by a volume integral over the free energy density of electric polarization in the space occupied by solvent: this approach3,22,28 is called density descreening. In the density descreening algorithms for offdiagonal Coulomb integrals, a volume integral is used to estimate descreening of one part of the solute by the presence of another part of the solute. Our current suite of descreening algorithms consist of density decreening by converged trapezoidal rule quadrature7 (DD, highest level), rectangular rule approximation to DD3,6,22 (RD, middle level), and pairwise descreening31 (PD, lowest level). In the present paper, we will use only the DD algorithm. The SM4 model8-11 was based on more accurate atomic charges, in particular on semiempirical class IV charges,37 which, for the purposes of calculating dipole moments or electrostatic potentials, appear from validation tests to be comparable to or better than the best partial charges available from any method that is currently practical for large molecules. The SM4 model has a new form for bond-order-dependent surface tensions and is again based on density descreening. It was parametrized for all alkane solvents for a wide variety of solutes containing H, C, N, O, F, S, Cl, Br, and I.8,9 For water it was parametrized solely on the basis of a limited range of solutes (hydrocarbons, ethers, aldehydes, alcohols, and diols) containing H, C, and O.10,11 The present paper reports a more general parametrization based on class IV charges suitable for a broad range of solutes containing H, C, N, O, F, S, Cl, Br, and I in water. The present data set provides a very demanding test of theory, since aqueous electrostatic solvation effects are very large. We note that not only are class IV charges more accurate than class III charges38-40 but when based on semiempirical molecular orbital theory as used here, they are very economical. An important feature of SM2 through SM4 is the use of calculated bond orders41 to identify the “nature” of the part of the solute in contact with the first solvation shell and hence to determine the local surface tension (this avoids the need for atomic “types”, as widely employed in molecular mechanics). In the course of developing the new SM5 general parametrization, we changed the functional forms of the surface tensions from bond-order-based to geometry-based because this yields a more stable parametrized model. This shift in functionalization is deemed significant enough to warrant designation as a distinct SMx model, i.e., SM5 instead of, for example, SM4C.10 In particular, we label the models presented in this paper SM5.4, where the 4 after the point denotes class IV charges. In future work, the present model will be seen as the first in a suite of SM5.x models combining various levels of atomic charges and varying degrees of “rigor” in the dielectric descreening algorithm. In the work presented here we use converged density descreening, and since this is the default choice for SM5-type models, we will not indicate this explicitly in the names of the new models presented here. In addition to combining class IV charges with better stability and a broader range of solutes, the SM5.4 model has three other advantages over earlier models: (i) the set of compounds used for parametrization is larger and more diverse, (ii) the performance of the method for various functional groups was analyzed in greater detail to try to identify the most physical functional form for the parametrization, and (iii) we present a single set of parameters that can be used with both AM1 and PM3 wave functions, rather than having separate parameters for each. Features i and ii will be true for the entire suite of SM5 models that will follow this one. Feature iii is possible because class

Model for Aqueous Solvation IV charges derived from AM1 and PM3 wave functions are reasonably similar. Thus it is tempting to postulate that the unspecific parameters (see below) of the present solvation model could also be used with ab initio calculations if they are of a high enough level (e.g., MP2/6-31G* class III charges) to yield comparably accurate charges. We have established the following notation for using the three sets of surface tension coefficients that will be presented in this article: SM5.4/U (with U for unspecific) will refer to the first set of coefficients developed, applicable to both AM1 and PM3. If we wish to denote that these are being used in a given calculation or set of calculations with the AM1 solute Hamiltonian and CM1A charge model, we will denote this usage as AM1-SM5.4/U. Similarly, PM3-SM5.4/U will denote usage of SM5.4/U with the PM3 solute Hamiltonian and the CM1P charge model. The surface tension coefficients derived specifically for use with the AM1 and CM1A models will be denoted SM5.4/A, and those derived specifically for use with the PM3 and CM1P models will be called SM5.4/P. We believe that the SM5 functional forms for the surface tensions account for the most significant geometry dependence of the first solvation shell effects that can be represented by simple functions of bond distances. Thus, further studies in our group will be the development of a range of levels of the treatment of solvation based on using these functional forms and the generalized Born approximation but with simpler prescriptions for the charges or descreening algorithm. One could also envision using these functional forms with other treatments of the electrostatics. We refer to the final models of this paper as SM5.4/A and SM5.4/P, and we refer to unspecific models obtained at an intermediate stage of the parametrization as AM1-SM5.4/U and PM3-SM5.4/U. Additionally, we will refer to the unspecific models simply as SM5.4/U when we are referring to both AM1SM5.4/U and PM3-SM5.4/U. Finally, we say SM5.4 to refer to all of these at once, and we say SM5 to discuss the functional form of the surface tensions or the whole suite of models employing them. The SM5.4 models presented here were developed using a locally modified version of version 5.0 of the AMSOL program,42 and they will be included in version 6.0 of AMSOL. Section 2 presents the list of solutes used for the parametrization of SM5 and the sources of the experimental data. Section 3 describes the parametrization of the SM5.4 aqueous models. Section 4 contains discussion, and Section 5 offers concluding remarks. 2. Experimental Data Neutrals. The set of neutral molecules that we used for the parametrization of the present aqueous solvation model is called the SM5.4 aqueous models neutral training set. The origin of this training set is most clear if we first define a larger set of molecules, called the meta set, that we defined for potential use to parametrize both aqueous and nonaqueous models. The meta set comes from several sources. First we included all neutral solutes that were used in creating any of the prior solvation models4-7 for water or the SM4 solvation model for n-hexadecane.8 The merger of these sets provides a broad coverage of commonly encountered functional groups. Since the meta set is designed to be used to parametrize not only aqueous and n-hexadecane solvation models but also solvation models for a wide variety of other solvents, it also includes some additional solutes that were chosen on the basis of the large amounts of data available in other solvents,43 in particular, 1-hexanol, o-cresol, m-cresol, p-cresol, 1-heptanol, 3,3-dimeth-

J. Phys. Chem., Vol. 100, No. 40, 1996 16387 ylbutanone, 2-heptanone, pentanoic acid, hexanoic acid, propyl ethanoate, methyl pentanoate, and methyl hexanoate. Finally, as we developed the SM5.4 aqueous model, we found it desirable to add more molecules to the meta set for various reasons. For instance, diethyl disulfide was added because the union of the previously used sets contained only one disulfide: dimethyl disulfide. Methylamine, 2-ethylpyrazine, and Nmethylmorpholine were added to aid in understanding trends in primary, secondary, and tertiary amines. Hydrazine was added to balance the otherwise strong effect that ammonia had on the optimization. Ammonia and hydrazine are the only compounds in the set that include nitrogen but not carbon. The original set contained morpholine and substituted piperazines; piperidine was added for balance. Our original set also did not include any amides; therefore we added acetamide, (E)-Nmethylacetamide, (Z)-N-methylacetamide, and N,N-dimethylacetamide. Our original set did not contain any fivemembered-ring aromatic amines; thus we added imidazole to the meta set. In order to obtain a greater understanding and also to predict experimental results of the effect of adding -CH2- groups to molecules containing a variety of functional groups, we added 1-octanol, octanal, methyl octanoate, and octanone. Finally, we added H2 to the meta set to improve the parametrization for H2 and hence for H; the latter may be important in kinetics. We used several sources for the experimental free energies of aqueous solvation. We took data from these references in the following order of preference: Suleiman and Eckert,44 Cabani et al.,45 Abraham and Whiting,46 Hine and Mookerjee,47 Wagman,48 Wolfenden,49 Wauchope and Haque,50 and Han and Bartles.51 All experimental data were converted to a consistent standard state, namely 1 mole per liter ideal gas for the vapor phase and 1 mole per liter ideal solution for the solution state. Now we can define the SM5.4 aqueous models neutral training set; this set includes all molecules in the meta set that satisfy all three of the following criteria: (i) an experimental free energy of aqueous solvation is available in the references listed in the previous paragraph, (ii) the molecule contains no elements other than H, C, N, O, F, S, Cl, Br, and/or I, and (iii) the molecule appears to be modelable in both the gas phase and solution in terms of a single conformation that is correctly predicted by AM1. The third (conformational) criterion is discussed further in section 3.1. The number of molecules in the SM5.4 aqueous models neutral training set is 215. Ions. All experimental data for ions were taken from Pearson.52 Data for several ions are “quite uncertain” according to Pearson;52 thus, these ions were not taken into consideration when creating the SM5.4 aqueous models ionic training set. The resulting number of molecules in the SM5.4 aqueous models ionic training set is 34. The most critical issue in assessing the reliability of Pearson’s results concerns the free energy of hydration assumed for the proton. Pearson used a value of -261.4 kcal for this quantity (this value, like all other values in the present paper, is converted to a standard state of 1 M in the gas phase as well as in solution). Ford and Wang,53 however, recommended the more recent value54 of -260.9 kcal for this quantity. Switching to this value changes Pearson’s free energies of solvation by +0.5 kcal for cations and -0.5 kcal for anions, which is within our rounding error for ions (since we round experimental solvation energies for ions to the nearest kilocalorie). Of course, there are additional sources of uncertainty as well, and some workers53,55 have made other estimates of solvation energies which may be compared to Pearson’s. Such comparisons are provided in Table 1. Comparing the Ford-Wang values and the corrected Klots

16388 J. Phys. Chem., Vol. 100, No. 40, 1996

Chambers et al.

TABLE 1: Experimental Free Energies of Solvation (kcal) for Some Ions correcteda

published ion

Pearsonb

Klotsc

Ford-Wangd

Pearson

Klots

NH4+ CH3NH3+ (CH3)2NH2+ (CH3)3NH+ (CH3)2OH+ CH3C(OH)CH3+ C5H5NH+ OHCH3COCH2CH3OC6H5OCH3CO2FClBrI-

-79 -70 -63 -59 -70 -64 -59 -106 -81 -95 -72 -77 -107 -77 -72 -63

-84

-79.5 -72 -64.5 -57 -73 -70 -57 -107.5 -76 -95 -76 -82

-78 -69 -62 -58 -69 -63 -58 -106 -81 -95 -72 -77 -107 -77 -72 -63

-80.5

-104

-104 -76 -72.5 -61

-107

-107 -79 -76 -65

a These values are corrected to a free energy of hydration of -260.9 kcal for H+, as recommended by Trasatti (ref 54) and Ford and Wang (ref 53). b Reference 52. c Reference 55. d Reference 53.

values to the uncorrected values of Pearson, which are the values we used, shows that in 17 of the 18 cases the deviation is less than or equal to the 5 kcal uncertainty that we estimated previously6 for ions. In fact, in 13 of the cases the deviations are less than or equal to 3 kcal, which is smaller than the mean accuracy of our method for ions. Thus, we believe that it is adequate to use the values of Pearson for fitting and comparison, and we shall do so in the rest of the paper. 3. Parametrization This section presents some of the issues arising in the interpretation of the experimental data, such as conformational issues for individual solutes, and it describes the functional forms used in the parametrizations. 3.1. Conformational Averaging. Many of the compounds in the meta set have multiple conformational minima on their respective potential energy surfaces. In principle, accurate prediction of the absolute free energy of solvation requires statistically averaging (in a Boltzmann sense) over all minima both in the gas phase and in solution,56 i.e.,

exp[-∆GS°/RT] ) ∑PC exp[-∆GS°(C)/RT]

(1)

C

where PC is the equilibrium mole fraction of conformation C in the gas phase. Evaluation of eq 1 can be tedious, especially for solutes characterized by a very large number of conformational minima. Furthermore, it is inconvenient to use during parametrization because it makes the solvation energies nonlinear functions of the surface tensions. Therefore, we avoided using eq 1 during parametrization. This is possible because the situation is simplified when a single conformer A dominates the equilibrium population in both the gas phase and solution or when all appreciably populated conformers have essentially the same solvation free energy. In the former case, it is justifiable to calculate the solvation free energy using only the dominant conformer. In the latter case, one may use any appreciably populated conformer. We took advantage of the first consideration to simplify the calculations in many instances by applying standard precepts of conformational analysis57 to eliminate consideration of highenergy conformers in the parametrization. For example, chair conformations were used for saturated six-membered rings

without consideration of twist-boats, equatorially substituted conformations of monosubstituted saturated six-membered rings were used instead of axially substituted ones, gauche interactions were minimized in branched hydrocarbon conformations, etc. We took advantage of the second consideration as well. For example, we assumed that relative conformations of substituents were unimportant in disubstituted aromatics where the substituents are in a 1,3 or 1,4 relationship, e.g., m-hydroxybenzaldehyde. We have also assumed this to be the case for conformers differing only by methyl group rotations. For such solutes, we simply used the lowest energy conformer at the AM1 and PM3 levels. In the same spirit, we have assumed all alkyl chains four carbon atoms or longer to be in a fully extended conformation. More rigorous studies employing full statistical averaging on n-alkanes in solution have indicated that “balling up” of the alkyl chain reduces the total solvent-accessible surface area by only a small amount.17,58 Although a nearby functional group may perturb the average by comparison to the alkane case, we retain this approximation for aliphatic substituents because of its simplicity. As another general issue, it is important to note that semiempirical methods are not always successful in predicting gas-phase conformational equilibria, so our choice of the dominant configuration is not always the AM1 or PM3 predicted global minimum. In prior studies, we have identified serious deficiencies for either or both AM1 and PM3 in predicting aspects of the conformational equilibria for cyclohexylammonium,59 glucose,11,60a and 1,2-ethanediol.60c Sometimes, semiempirical minima do not have overall geometries that are even comparable to more accurate ones; for example, it is well established that AM1 and PM3 predict saturated five-membered rings to be much closer to planar than they actually are.61-63 For these latter cases, we have made the assumption that the free energy of solvation will be only minimally affected by the flattening of the ring. Moreover, this flattening prevents us from addressing issues like envelope vs twist conformations or pseudoequatorial vs pseudoaxial substitution at the semiempirical level; we simply assume that their effects in a more rigorous conformational averaging would be small. For cyclopentanol, AM1 and PM3 predict that the endo orientation of the hydroxyl group is preferred, but we used an exo orientation that is presumably more similar to the actual equatorial orientation in the more puckered ring. All that being said, several solutes remain for which conformational issues were addressed as special cases or at least explicitly checked for potential problems. Some examples are discussed in the rest of this section. For 1,2-dimethoxyethane and 2-methoxyethanol, we explicitly checked that the trans orientations of the oxygens are lower than the gauche in both gas phase and solution for both AM1 and PM3. We found this is indeed the case, by 0.6-2.8 kcal. We checked that in propylamine, the C-2 carbon is trans to the nitrogen lone pair in both the gas phase and solution; AM1 and PM3 each predict it to be lower in energy in both phases by 2.9-3.9 kcal. In 2,2,2-trifluoroethyl vinyl ether, AM1 and PM3 predict that the s-cis conformation is more favorable than the s-trans by 1.5-2.1 kcal both in the gas phase and in solution (conformation about the H2C-O bond is trans as expected). In diethylamine, we based our analysis on the gas-phase global minimum found using both AM1 and PM3; the minimum obtained by both of these methods places both methyl groups anti to the nitrogen lone pair. Although other low-energy structures were found, the solvation energies of all of these structures are reasonably similar so conformational averaging

Model for Aqueous Solvation was omitted. The situation for dipropylamine was also checked thoroughly and found to be similar to that for diethylamine. For 2-ethylpyrazine, the structures with CH3-CH2-C-N dihedral angles of 60° and 120° are very similar in energy in both the gas phase and solution. We used the 60° structure. For hydrazine we used a structure with gauche lone pairs, which in AM1 is lower than trans in both the gas phase and solution and which also agrees with the experimental64 gasphase conformation, which has an H-N-N-H dihedral angle of 91°. For PM3, the gauche conformation does not exist in the gas phase or in solution. Thus, hydrazine was included in the parametrization only for the AM1 solute description, whereas every other molecule in the training set was included for both AM1 and PM3. Finally, some molecules in the meta set were not included in the SM5.4 aqueous models neutral training set because checks like this indicated a more complicated situation. Specifically, 1,2-dichloroethane, 1,2-dibromoethane, 2,2-dichloro-1,1-difluoroethyl methyl ether, 3-ethyl-2-methoxypyrazine, trimethyl orthoacetate, trimethyl orthobenzoate, and triethylamine were not used in the parametrization of the surface tensions. 3.2. Atomic Radii. The importance of the choice of effective atomic radii for the calculation of solvation free energies is controversial. In the present work, considerable effort was expended not only in searching for appropriate electrostatic radii but also in examining the sensitivity to our choices. Although the electrostatic free energies for individual ions and molecules show some sensitivity to the intrinsic Coulomb radii for all atom types, if one uses a large parametrization set and finds the surface tension coefficients by linear regression for each set of intrinsic Coulomb radii, most of the sensitivity vanishes in the final root-mean-square (rms) errors (the major exceptions being the radii for N and S). Nevertheless, as mentioned in the Introduction, one should try to keep the radii as realistic as possible in order to get physical wave functions. Thus, except for N and S, we took the intrinsic Coulomb radii for nonhalogen atom types to have a systematic physical relationship to the actual sizes of the atoms; for the halogen radii, the values of the radii were selected solely on the basis of the predicted free energies of solvation for the monatomic ions. For C and O we set the radii 0.08 Å larger than Bondi’s values65 for the van der Waals radii in an attempt to make the partition between electrostatic and nonelectrostatic contributions as meaningful as possible. We started with the same 0.08 Å criterion for H, but as will be shown in greater detail below, we found it necessary to “switch” to a smaller radius for H when it is bonded to O or N. Since the effects of the electrostatics are much greater in the ions, they were carefully studied in conjunction with the neutrals to aid in our selection of atomic radii. Both S and N radii needed to be more than 0.08 Å larger than Bondi’s van der Waals radii. When changing the electrostatic radius for S from 1.78 to 1.92 to 2.1 Å, the mean signed error, mean unsigned error, and rms error over the subset of the SM5.4 aqueous models neutral training set that contains sulfur but no halogens all changed by 0.01 kcal or less. Similarly, when changing the electrostatic radius for N from 1.78 to 1.85 to 1.92 Å (we did not want the electrostatic radius for N to be larger than S, thus we did not use radii larger than 1.92 Å), the mean signed error, mean unsigned error, and rms error over the nitrogen-containing subset of the SM5.4 aqueous models neutral training set all changed by 0.01 kcal or less. However, the ions showed clear preferences for particular choices of the radii. Mean unsigned and rms errors for the above listed S radii over all sulfur containing ions varied by 1.4 and 1.9 kcal, respectively. Mean

J. Phys. Chem., Vol. 100, No. 40, 1996 16389 TABLE 2: Elctrostatic Radii (Å), van der Waals Radii (Å), and Microscopic Surface Tension Coefficients (cal mol-1 Å-2) for SM5.4 by Atom Type and Atom Pair SM5.4/A k

σk(0)

H C O N F S Cl Br I

20.00 63.70 -80.35 3.74 18.89 -43.28 -2.29 -6.84 -7.96

SM5.4/P σk(0)

σHk′ -21.21 38.66 -59.05 50.22

σHk′

20.00 59.37 -97.46 18.35 18.35 -37.08 -1.91 -6.82 -9.39

k

k′

(1) σkk′

O O O C S N

O N C C S C

38.19 187.05 112.90 -36.19 35.87 -33.15

-20.72 61.54 -50.76 28.15

(2) σkk′

-23.21

SM5.4/A and SM5.4/P Rk

Fk(0)

FHk′

1.2 1.7 1.52 1.55 1.47 1.8 1.75 1.85 1.98

1.28 1.78 1.6 1.92 1.5 1.92 2.13 2.31 2.66

-1.28 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0

(1) σkk′

25.16 161.73 146.31 -44.56 38.09 -42.03

(2) σkk′

-13.88

unsigned and rms errors for the above listed N radii over all nitrogen-containing ions varied by 0.7 and 1.0 kcal, respectively. The final values of the electrostatic radii chosen on the basis of these considerations are listed in Table 2. A comment is also in order on the choice of 1.7 Å as the effective solvent radius. This number may be “justified” in various ways. First of all, it is an average of the values 1.4 and 2.0 Å that we previously used for water and alkanes. Second, it may be calculated from the density of water by assuming that water consists of close packed spheres (which occupy 68% of the volume). Third, we parametrized a large set of compounds in the SM5.4 aqueous models training set containing H, C, N, and O using solvent radii of 1.4, 1.65, 1.7, and 1.9 Å and, as long as the surface tension coefficients are found by linear regression using areas calculated with a given solvent radius, the change in results is negligible. In particular the change in rms error is on the order of 0.02 kcal, and the solvation free energy of 99% of the molecules changes by less than 0.15 kcal. Thus 1.7 Å is the center of a stable range for water, and we saw no reason to retain the traditional value of 1.4 Å,16 which is based on the O-O pair correlation function in water. The change from the traditional value helps to emphasize that in a continuum model we make no assumptions about the orientation of the O atoms in the first-shell water molecules. Our interpretation of the present model is that the effective solvent radius is closely related both to the range of dispersion interactions and to the size of a solvent molecule. In this regard it is important to remember that free energies of solvation represent averages over a myriad of solvent configurations, and the temptation to associate statistically averaged quantities with easily visualized pictures of individual solvent molecules with specific sizes, shapes, and orientations should be strongly resisted. Even such appealingly named features as the “solvent-separated” dimer or “solvent-separated” pair result from a complex mix of configurations, as shown for example by Jorgensen et al.66 We have concluded, therefore, that within reasonable bounds and with due attention to keeping the wave functions physical consistency is more important than a specific molecular model in choosing the solvent radii; the same considerations apply to the solute van der Waals radii, for which we used the standard values of Bondi,65 and to some extent the intrinsic Coulomb radii, which have already been discussed in greater detail. It is worth repeating here a point made in the Introduction: the

16390 J. Phys. Chem., Vol. 100, No. 40, 1996

Chambers et al.

boundary between the solute and the dielectric medium is an intrinsically artificial concept. The surface tensions will give a different contribution from the first solvation shell depending on how much of that shell’s electrostatics are included in the electrostatic term, which depends on the intrinsic Coulomb radii chosen. Any attempt to interpret the microscopic surface tensions that does not take this into account misses the point, except possibly for alkane solutes, where electrostatic effects are negligible. 3.3. Functional Forms. We arrived at the new functional forms by an iterative process of trial and error and by comparing trends in predicted free energies of solvation to experiment for various functional groups and substitution patterns. We learned in this way which parametric dependencies need to be explicitly included in the functional forms and which dependencies will result directly from the form of the electrostatics and solventaccessible surface areas alone. We also learned which explicit dependencies were unstable and hence to be avoided or were highly correlated with others and hence unnecessary if the others were included. In our opinion this exhaustive process (as opposed to some ideal artificial intelligence method such as a neural network or a genetic algorithm) is still the best way to develop forms for parametrizing semiempirical theories when the domain of intended applicability is broad, such as the present work, where we include a very wide variety of organic functionalities. (In fact we did try to use genetic algorithms67 for several aspects of the problem, and in later work we found them to be quite useful. We will not discuss the wide variety of parametrization functionalities that did not survive the winnowing process or attempt to rationalize the optimality of the final choices, but will simply present the final functional forms and the results we obtained with them as their justification. The rest of this section presents these forms. In all SMx models, the standard-state free energy of solvation is written as3-11,68

∆GS° ) ∆GENP + GCDS°

(2)

where ∆GENP includes the change in the electronic and nuclear internal energy of the solute and the electric polarization free energy of the solute-solvent system upon insertion of the solute in the solvent and GCDS° is the contribution of local (e.g., first solvation shell) effects to the standard-state free energy of transfer. Note that S in ∆GS° stands for solvation, ENP stands for electronic-nuclear polarization, and CDS stands for cavitation-dispersion-solvent structure, as explained elsewhere.4-9 The ∆GENP term can be written as

∆GENP ) ∆EEN + GP

(3)

where ∆EEN is the change in the electronic and nuclear energy of the solute in going from the gas phase to solution and GP is the polarization free energy. The standard-state free energy in solution is approximated as

GS° ) E(g) + ∆GS°

(4)

where E(g) is the Born-Oppenheimer solute energy in the gas phase. We note that in principle GS° should include contributions from zero-point and thermal vibrational motions and from thermal population of electronically excited states and gas-phase rotationally excited states, but in practice these are included only in applications where very high accuracy is attempted,11,60 and they will not be discussed further in the present paper.

A key element of our parametrization of ∆GENP and GCDS° is the use of a new form of switching function. The function, denoted a “cutoff tanh function” (or COT) because of its general shape, has the form

{

T(R|R h ,∆R))

[ (∆R -∆RR + Rh )]

exp -

ReR h + ∆R

0

otherwise

(5)

where R is a distance, and R h and ∆R are parameters. In all cases R will be equal to the distance Rkk′ between two atoms k and k′. In the SM5.4 models, GP has the form22

GP ) -

( )∑

1 2

1-

1 

qkqk′(rkk′2 +

k,k′

RkRk′ exp[-rkk′2/dkk′RkRk′])-1/2 (6) where rkk′ is the interatomic distance between atoms k and k′, Rk is an effective atomic radius for atom k in the electrostatic part of the calculation and will be called a Coulomb radius (more precisely, this is the effective Coulomb radius, and it depends on the intrinsic Coulomb radius, which is parametrized), and dkk′ is an empirically optimized constant equal to 4.0 for all k and k′ except when one of k and k′ is H and the other is C, in which case the value for dkk′ is 4.2. Two major methodological changes were made in the SM5 formalism as compared to the formulation used in models SM1,3 SM2,4,6 SM2.1,7 SM3,5,6 SM3.1,7 and SM4.8,9 These changes predominately affect only molecules containing C, H, N, and O. The first change is the removal of the two cutoff Gaussians3,6 in the electrostatics. These cutoff Gaussians primarily affected systems with two geminal oxygens or with a vicinal N and H pair. The second change in the model involves C, H pairs. Previous models produced positive ∆GENP values, which are unphysical, for longer unbranched alkanes and some branched and cycloalkanes. Although these values were small, they were clearly wrong and so we have removed this problem. We previously had set dkk′ (eq 6) equal to 4 in all cases. In the SM5 approach we retain this value except when k and k′ correspond to a C, H pair. This eliminated all positive ∆GENP values. The values of the effective atomic radii are calculated using the intrinsic Coulomb radii that are called F.3 For C, N, O, F, S, Cl, Br, and I, Fk is a constant, called FC, FN, etc. For H,

Fk|k)H ) FH(0) + FHk′

∑

T[-T(Rkk′|-0.4,0.4)|R h Hk′, ∆RHk′]

k′)O,N

(7)

This functional form ensures a somewhat smaller radius for hydrogen atoms bonded to nitrogen or oxygen than for other hydrogen atoms. Similar ideas, namely, that hydrogens bonded to heteroatoms have smaller radii than hydrogens bonded to carbon, have also been widely used in molecular mechanics, and this idea has been systematically validated in a solvation context by Luque, Negre, Orozco, and Bachs.69 In our own earlier models this same effect was accomplished more indirectly by making the effective radius an inverse function of atomic charge.3,6 The first solvation shell term has the form

GCDS° ) ∑σkAk(RSCDS)

(8)

k

where k denotes an atom, Ak(RSCDS) is the solvent-accessible surface area of atom k, RSCDS is the effective radius used for a

Model for Aqueous Solvation

J. Phys. Chem., Vol. 100, No. 40, 1996 16391

solvent molecule (the radius used for water is 1.7 Å, as discussed in section 3.2), and σk is the microscopic surface tension of the atom. The electrostatic cutoff Gaussians for O-O interactions which were present in the previous models were inserted because otherwise agreement with experiment was poor. Similarly poor results for compounds with geminal O-O interactions were found in early stages of the development of the SM5 model (if the cutoff Gaussians were not included). To obtain better agreement with experiment without using cutoff Gaussians in the electrostatic term, a microscopic surface tension, σOO, was added, which comes into effect for systems with geminal O-O interactions, such as carboxylic acids, esters, and orthoesters. Additionally, microscopic surface tension terms have been included for N-C, O-N, and S-S interactions. We also added a microscopic surface tension for H-bondedto-O which affects the alcohols and carboxylic acids, and we added an O-C microscopic surface tension that predominantly affects the aldehydes, ketones, carboxylic acids, and esters. The O-C microscopic surface tension also makes a small contribution to GCDS° for a few other molecules in the training set, notably 1,4-dioxane and some alcohols. We added two C-C microscopic surface tension terms; the first affects all molecules with C-C bonds. The second C-C surface tension term affects mainly carbon-carbon triple bonds, but some carbon-carbon double bonds have contributions that are approximately 1/40th of the amount for the triple bonds. Cavitation-dispersion-structure terms were also computed on the basis of oxygen-bonded-to-nitrogen, hydrogen-bondedto-nitrogen, and nitrogen-bonded-to-carbon distances. The oxygen-bonded-to-nitrogen term only affects the nitrohydrocarbons. We also found it necessary to include a surface tension term of a more complicated form for nitrogen bonded to carbon. We found it necessary to include a microscopic surface tension for sulfur bonded to sulfur. There are only two disulfides in the test set; however, without the sulfur-bondedto-sulfur surface tension term, the SM5.4 models predict free energies of solvation that are too negative for the disulfides by 2.0-2.7 kcal. As a consequence of the considerations just enumerated, the forms used σk for in eq 8 are

σk|k)H ) σH(0) +

∑

σHk′T(RHk′|R h Hk′, ∆RHk′)

(9)

k′)C,O,N,S 2

(i) (i) (i) , ∆RCC ) σk|k)C ) σC(0) + ∑σCC ∑ T(Rkk′|Rh CC i)1

σk|k)O ) σO(0) + σOC

(10)

k′)C k′*k

∑ T(Rkk′|Rh OC, ∆ROC) +

k′)C

σOO ∑ T[-T(Rkk′|-0.4, 0.4)|R h OO, ∆ROO] + k′)O k′*k

∑ T(Rkk′|Rh ON, ∆RON)

σON

(11)

k′)N

h SS, ∆RSS) σk|k)S ) σS(0) + σSS ∑ T(Rkk′|R

(12)

k′)S k′*k

σk|k)N ) σN(0) +

h NC, ∆RNC)[ ∑ T(Rk′k′′|R h NC, ∆RNC)]2}1.3 σNC { ∑ T(Rkk′|R k′)C

k′′*k k′′*k′

(13)

TABLE 3: COT Parameters (Å) k

k′

R h kk′

∆Rkk′

H H H C

O or N C S C

O O O N S

C O N C S

1.7 1.85 2.14 1.84a 1.27b 1.33 2.75 1.5 1.84c 2.75

0.3 0.3 0.3 0.3a 0.07b 0.1 0.3 0.3 0.3 0.3

(1) (1) b (2) (2) R h CC′ , ∆RCC′ . R h CC′, ∆RCC′ . 1.55 Å. a

c

For k′′ ) H in eq 13 this value is

The value for σH(0) was adjusted so that SM5.4 correctly predicts the free energy of solvation for H2.51 The centers, R h kk′, of the COTs were chosen on the basis of bond-order functions and distances between the atoms, whereas the widths, ∆Rkk′, of the COTs were determined by systematic searches. For each trial choice of these parameters we fit the microscopic surface tension coefficients by linear regression to minimize the rms errors for the neutrals. We developed the model using a single set of parameters for both AM1 and PM3. We shall refer to this model as either AM1-SM5.4/U or PM3-SM5.4/U, depending on which Hamiltonian is used in the gas phase. The “U” refers to the unspecific nature of the surface tension coefficients parametrized using both Hamiltonians. Table 3 contains the COT parameters R h kk′ (1) and ∆Rkk′. The SM5.4/U parameters for σ(0) k , σHk′, σkk′ , and are given in the supporting information. We then refit the microscopic surface tension coefficients for AM1 and PM3 individually, resulting in two more refined models that we shall refer to as SM5.4/A and SM5.4/P. The values for F(0) k , FHk′, (1) (2) σ(0) , and σ for the SM5.4/A and SM5.4/P models , σ , σ Hk′ k kk′ kk′ are in Table 2. 4. Discussion 4.1. General Remarks. We begin the discussion by making a few remarks on the surface tension terms. First we note that the microscopic surface tensions σγ are named on the basis of their units, and in no way should this be regarded as indicating that they are the same as macroscopic surface tensions. In particular, the structural reordering of solvent near a plane interface with either a dilute gas or an immiscible fluid would be expected to be different from the structural reordering near a solute, and this has a very significant effect on the local entropic contributions. In fact, one would expect that the free energy per unit of area would depend significantly on the nature of the exposed surface. The essence of semiempirical modeling is to reduce this to a manageable functionality. In the SM2SM4 series, we modeled this dependence using atomic numbers and ground-state bond orders, the latter being functions of molecular geometry. In the SM5 approach we use atomic numbers and molecular geometry directly. We have examined a large number of possible bond-order and geometric dependencies, and we settled on the functional forms used here as the simplest and most stable that capture the systematic deviations of the experimental free energy of solvation from our electrostatic estimates. The concept of an environment-dependent surface tension is intuitively obvious, but it can also be backed up by simulations. For example, Wallqvist and Berne70 have studied the free energy of a spherical probe that interacts with water by an r-12 repulsive potential, and they studied probes with various surface curva-

16392 J. Phys. Chem., Vol. 100, No. 40, 1996

Chambers et al.

TABLE 4: Mean Unsigned Errors (kcal/mol) in Aqueous Free Energies of Solvation by Functional Group for Selected SMx Models functional group unbranched alkanes branched alkanes cycloalkanes alkenes alkynes arenes alcohols ethers aldehydes ketones carboxylic acids esters bifunctional water, H2 subtotal

number of molecules 8 5 5 9 5 8 16 9 6 12 5 12 5 2 107

AM1-SM2.1

PM3-SM3.1

SM5.4/A

SM5.4/P

experimental dispersion

Compounds Containing at Most C, H, and/or O 0.55 0.52 0.63 0.61 0.32 0.37 0.22 0.19 0.52 0.67 0.69 0.90 0.56 0.89 1.46 1.44 0.37 1.22 0.80 0.33 0.44 0.61 0.54 1.16 0.49 0.66 1.24 1.23 0.62 0.77

0.60 0.65 0.21 0.45 0.25 0.14 0.51 0.93 0.29 0.37 0.77 0.44 0.32 1.59 0.49

0.56 0.69 0.13 0.33 0.21 0.18 0.40 1.03 0.36 0.37 0.75 0.41 0.34 1.17 0.46

3.06

0.73 0.39 0.49 0.49 1.84 0.50 2.20 0.72

0.74 0.56 0.48 0.14 0.50 0.68 0.75 0.57

1.85

aliphatic amines aromatic amines nitriles nitrohydrocarbons amides bifunctional ammonia, hydrazine subtotal

15 10 4 6 3 3 2a 43a

Compounds Containing N 1.08 2.21 1.00 1.35 0.27 1.85 0.50 1.24 0.57 0.53 1.28 2.38 3.15 0.18 0.98 1.67

thiols organic sulfides, H2S organic disulfides subtotal

4 5 2 11

Compounds Containing S, H, and/or C 0.68 0.37 0.78 0.48 2.20 4.37 1.00 1.15

0.30 0.40 0.20 0.33

0.24 0.35 0.05 0.26

0.58

3 8 5 3 10 8 17 54

Compounds Containing Halogens 0.62 0.73 0.23 0.50 0.21 0.61 0.30 0.90 0.44 0.44 0.38 0.58 0.60 0.80 0.43 0.64

0.54 0.27 0.65 0.13 0.25 0.23 0.48 0.37

0.47 0.23 0.55 0.19 0.24 0.15 0.41 0.32

1.56

0.50

0.44

2.86

fluorinated hydrocarbons chloroalkanes chloroalkenes chloroarenes brominated hydrocarbons iodinated hydrocarbons other halo molecules subtotal

Entire Set total

215a

0.67

0.93

a

Hydrazine was not calculated using either the PM3-SM3.1 or the SM5.4/P method (see text). Therefore, the number of elements for PM3 methods is one less than the value shown in the number of molecules column.

tures and shapes. They found that the free energy change per first-shell solvent molecule depends (although weakly) on the characteristics of the solute surface. This dependence is consistent with our use of environment-dependent surface tensions. Their model did not include attractive dispersion forces or solute internal energy, which further complicate real systems. One of the most severe physical simplifications of the present models, and hence the entire suite of SM5 models, is that the surface tensions depend on local charge density only through the atomic number and selected geometrical variables and notsfor examplesexplicitly on partial charges.71 Aside from these dependencies, this corresponds to postulating that first-solvation shell effects of varying partial charges are all entirely accounted for by ∆GENP. The SM5 functional forms are strictly functions of geometry, requiring no assignment of atom types, hybridization states, or bond orders. As a consequence, one never encounters a structure that suffers from the “missing parameter” problem of molecular mechanics. An eloquent statement of the disadvantages of atom typing and the assignment of chemical functionality, especially in the context of the richness of structures encountered in the pharmaceutical research environment, has

been provided by Gerber and Mu¨ller.72 An even more important motivation in our own case was to develop a parametrization that remains valid as the system proceeds along a reaction path. 4.2. Specific Comments on Parametrization. Tables 4 and 5 contain the mean unsigned and mean signed errors, respectively, over all functional groups included in the SM5.4 aqueous models neutral training set for four different solvation models: SM2.1, SM3.1, SM5.4/A, and SM5.4/P. To make the errors more meaningful, we have also included the number of molecules in the SM5.4 aqueous models neutral training set that fall into each group. Finally, in each of these two tables, we have also listed the dispersion of the experimental data for the subgroups of molecules containing various heteroatoms and for the entire set; for this purpose the dispersion of a set is defined as the rms deviation of the members of that set from their mean value. Table 6 gives SM5.4/A and SM5.4/P values for the ENP and CDS portions of the free energy of solvation for several individual molecules from each subgroup. Table 7 contains the ions used in the parametrization. Table 8 contains neutrals and ions that were not used in the parametrization of the SM5.4 aqueous models for various

Model for Aqueous Solvation

J. Phys. Chem., Vol. 100, No. 40, 1996 16393

TABLE 5: Mean Signed Errors (kcal/mol) in Aqueous Free Energies of Solvation by Functional Group for Selected SMx Models functional group unbranched alkanes branched alkanes cycloalkanes alkenes alkynes arenes alcohols ethers aldehydes ketones carboxylic acids esters bifunctional water, H2 subtotal aliphatic amines aromatic amines nitriles nitrohydro carbons amides bifunctional ammonia, hydrazine subtotal thiols organic sulfides, H2S organic disulfides subtotal fluorinated hydrocarbons chloroalkanes chloroalkenes chloroarenes brominated hydrocarbons iodinated hydrocarbons other halo molecules subtotal

number of molecules 8 5 5 9 5 8 16 9 6 12 5 12 5 2 107 15 10 4 6 3 3 2a 43a

AM1-SM2.1

PM3-SM3.1

Compounds Containing at Most C, H, and/or O -0.55 -0.52 -0.63 -0.61 0.32 0.37 -0.18 -0.04 -0.52 -0.67 0.69 0.90 0.38 0.86 1.15 1.19 -0.36 -1.22 0.75 0.06 0.01 -0.61 -0.23 -1.16 -0.01 -0.21 -1.24 -1.23 0.13 -0.04 Compounds Containing N 0.43 1.33 0.61 0.81 -0.09 -1.85 -0.50 1.24 0.45 -0.24 0.89 1.60 -3.15 -0.18 0.16 0.76

4 5 2 11

Compounds Containing S, H, and/or C 0.31 0.23 0.09 -0.21 -2.20 -4.37 -0.25 -0.81

3 8 5 3 10 8 17 54

Compounds Containing Halogens 0.62 0.73 -0.07 0.44 0.00 -0.01 -0.30 -0.90 0.11 -0.27 -0.32 0.37 0.32 -0.16 0.08 0.01

SM5.4/P

experimental dispersion

-0.60 -0.65 -0.15 0.33 -0.04 0.10 0.29 0.10 -0.29 0.29 0.77 -0.44 0.01 -1.56 -0.02

-0.53 -0.69 -0.06 0.25 -0.01 0.12 0.16 0.30 -0.36 0.30 0.75 -0.41 0.01 -1.13 -0.01

3.06

0.18 0.05 0.29 0.21 1.44 -0.10 0.17 0.23

0.24 0.03 0.13 0.10 0.42 -0.68 0.75 0.12

1.85

0.03 0.09 0.03 0.06

0.05 0.04 -0.01 0.03

0.58

0.54 -0.24 0.63 -0.10 -0.14 -0.07 -0.04 0.00

0.47 -0.20 0.52 -0.19 -0.13 -0.04 0.00 0.00

1.56

0.04

0.02

2.86

SM5.4/A

Entire Set total

215a

0.10

0.09

a

Hydrazine was not calculated using either the PM3-SM3.1 or the SM5.4/P method (see text). Therefore, the number of elements for PM3 methods is one less than the value shown in the number-of-molecules column.

reasons. We present these molecules for discussion purposes since they each involve some aspect of interest. 4.3. Accuracy Obtained for Neutrals. The SM5.4 models predict the free energy of methane more accurately than do previous models; see Table 4. In addition, there is a systematic improvement for all molecules containing methyl groups (such as methanol and methyl methanoate). This gain comes without sacrificing the accuracy of the model for the hydrocarbons; see Tables 4 and 5. The SM5.4 models yield lower mean unsigned, mean signed, and rms errors than the SM2.1 or SM3.1 models for the CHO test set. An additional experimentally observed trend that we wanted the models to reproduce is the change in free energy of solvation due to the insertion of a -CH2- group. To study this trend, we selected the following pairs of molecules (experimental differences of the butyl minus the octyl are given parenthetically in kilocalories after each set): n-butane and n-octane (-0.8), 1-butanol and 1-octanol (-0.6), butanal and octanal (-0.9), butanone and octanone (-0.7), and methyl butanoate and methyl octanoate (-0.8). SM5.4/A predicts these differences to be -0.7, -0.7, -0.7, -0.9, and -0.7 kcal, and SM5.4/P predicts them to be -0.7, -0.7, -0.7, -0.8, and -0.7 kcal.

The SM5.4 models incorporate systematic improvements over previous models for several functional groups, especially alkynes, ethers, and aliphatic amines. In addition, the treatment of molecules containing S, Br, and I is significantly improved. The rms errors over the entire SM5.4/A and SM5.4/P predictions for the neutral training sets are 0.72 and 0.62 kcal/ mol, respectively. The rms errors for subsets containing individual functional groups are given in the supporting information. 4.4. Accuracy Obtained for Ions. Table 7 gives the calculated values for ∆GENP, GCDS°, and ∆GS° for all 34 ions used in the parametrization of both SM5.4/A and SM5.4/P. The dispersion of the experimental data over the entire set is 13.5 kcal, whereas the mean unsigned errors for SM5.4/A and SM5.4/P predictions are 4.4 and 4.3 kcal, respectively. The rms errors over the SM5.4/A and SM5.4/P predicted sets are 5.7 and 5.4 kcal, respectively. 4.5. Molecules Not used in the SM5.4 Aqueous Models Parametrization Set. For discussion purposes, we present in Table 8 results predicted by SM5.4 models for some molecules that were not included in the training set.

16394 J. Phys. Chem., Vol. 100, No. 40, 1996

Chambers et al.

TABLE 6: Calculated and Experimentala Free Energies of Solvation (kcal/mol) for Selected Neutral Solutes Used in the Parametrization of SM5.4/A and SM5.4/P functional groups and molecules

∆GENP

SM5.4/A GCDS°

∆GS°

∆GENP

SM5.4/P GCDS°

∆GS°

experiment ∆GS°

methane n-butane n-octane

-0.1 -0.2 -0.1

2.1 1.7 2.2

Unbranched Alkanes 2.0 -0.0 1.4 -0.1 2.1 -0.0

2.1 1.6 2.2

2.1 1.5 2.2

2.0 2.1 2.9

2,2-dimethyl propane 2,2,4-trimethyl pentane

-0.2

2.0

Branched Alkanes 1.8 -0.1

1.8

1.8

2.5

-0.1

2.5

2.4

-0.0

2.3

2.3

2.9

cyclopentane cis-1,2-dimethylcyclohexane

-0.6 -0.2

1.4 1.9

Cycloalkanes 0.8 1.8

-0.2 -0.1

1.3 1.8

1.1 1.8

1.2 1.6

1-butene cyclopentene (E)-2-pentene

-0.8 -1.4 -0.9

2.4 1.6 2.0

Alkenes 1.6 0.2 1.2

-0.4 -0.8 -0.5

2.0 1.1 1.7

1.5 0.3 1.2

1.4 0.6 1.3

ethyne 1-hexyne

-2.3 -2.4

2.8 2.6

Alkynes 0.5 0.2

-1.9 -2.0

2.3 2.2

0.5 0.2

0.0 0.3

benzene m-xylene naphthalene

-3.1 -3.2 -4.9

2.1 2.4 2.6

Arenes -1.0 -0.8 -2.3

-1.9 -2.0 -3.1

0.8 1.4 0.7

-1.1 -0.7 -2.4

-0.9 -0.8 -2.4

1,2-ethanediol 1-butanol cyclopentanol p-cresol 1-octanol

-6.4 -3.7 -3.1 -6.7 -3.5

-2.3 -0.7 -0.8 0.2 -0.2

Alcohols -8.8 -4.4 -3.9 -6.5 -3.7

-6.3 -3.6 -3.0 -5.1 -3.5

-2.9 -1.0 -1.2 -1.1 -0.5

-9.2 -4.6 -4.2 -6.2 -4.0

-9.3 -4.7 -5.5 -6.1 -4.1

tetrahydrofuran 1,4-dioxane methyl propyl ether 1,2-dimethoxyethane anisole

-2.9 -3.9 -2.1 -3.3 -4.4

-0.6 -2.3 0.7 -0.3 1.4

Ethers -3.6 -6.2 -1.4 -3.6 -3.0

-2.2 -2.9 -1.6 -2.4 -2.8

-1.2 -3.4 0.4 -0.8 0.0

-3.4 -6.3 -1.2 -3.2 -2.8

-3.5 -5.1 -1.7 -4.8 -1.0

butanal benzaldehyde octanal

-4.3 -6.1 -4.2

1.0 1.7 1.6

Aldehydes -3.3 -4.4 -2.6

-4.5 -5.3 -4.5

1.2 0.7 1.8

-3.3 -4.6 -2.6

-3.2 -4.0 -2.3

butanone cyclopentanone 3-pentanone methyl phenyl ketone 2-octanone

-4.6 -4.3 -4.3 -6.3 -4.3

1.0 0.4 1.3 1.6 1.6

Ketones -3.6 -3.9 -3.0 -4.7 -2.7

-4.6 -4.5 -4.3 -5.4 -4.4

1.1 0.6 1.3 0.7 1.7

-3.6 -3.9 -3.0 -4.7 -2.7

-3.6 -4.7 -3.4 -4.6 -2.9

ethanoic acid pentanoic acid

-8.0 -7.0

1.5 1.8

Carboxylic acids -6.5 -7.6 -5.2 -6.6

1.1 1.4

-6.6 -5.3

-6.7 -6.2

methyl methanoate methyl butanoate butyl ethanoate methyl octanoate

-6.0 -4.8 -5.3 -4.7

2.1 1.9 2.2 2.5

2-propen-1-ol 2-methoxyethanol butenyne m-hydroxybenzaldehyde

-4.3 -5.4 -3.2 -9.5

methylamine dimethylamine propylamine trimethylamine pyrrolidine diethylamine N-methylpiperazine N,N′-dimethylpiperazine dipropylamine

-2.5 -2.0 -1.8 -1.8 -1.9 -1.1 -3.0 -2.6 -1.1

Esters -3.8 -2.9 -3.1 -2.2

-5.5 -4.4 -4.9 -4.3

1.6 1.5 1.9 2.1

-3.9 -2.9 -3.0 -2.3

-2.8 -2.8 -2.6 -2.0

-0.1 -1.8 3.1 -0.2

Bifunctional CHO -4.4 -7.2 -0.1 -9.7

-3.7 -4.9 -2.3 -8.2

-0.8 -2.4 2.2 -1.4

-4.5 -7.3 -0.2 -9.6

-5.1 -6.8 0.0 -9.5

-3.2 -2.7 -2.6 -1.3 -3.1 -1.6 -4.1 -2.5 -1.8

Aliphatic Amines -5.8 -4.7 -4.4 -3.2 -5.1 -2.8 -7.1 -5.1 -2.9

-1.7 -1.0 -1.4 -0.3 -1.1 -0.8 -1.4 -0.5 -0.7

-3.4 -3.9 -2.7 -3.0 -4.2 -2.4 -5.5 -4.1 -2.5

-5.1 -4.9 -4.1 -3.3 -5.3 -3.2 -6.9 -4.6 -3.2

-4.6 -4.3 -4.4 -3.2 -5.5 -4.1 -7.8 -7.6 -3.7

Model for Aqueous Solvation

J. Phys. Chem., Vol. 100, No. 40, 1996 16395

TABLE 6 (Continued) functional groups and molecules

∆GENP

SM5.4/A GCDS°

∆GS°

Aromatic Amines -0.1 -5.0 -1.0 -5.7 0.7 -4.4

pyridine 2-ethylpyrazine 2,5-dimethylpyridine

-4.8 -4.7 -5.1

ethanonitrile benzonitrile

-6.1 -5.1

1.8 2.2

nitroethane 2-nitropropane 2-methyl-1-nitrobenzene

-9.2 -8.6 -8.1

Nitrohydrocarbons 5.0 -4.1 5.2 -3.4 5.2 -2.8

(E)-N-methylacetamide

-6.9

-0.8

Nitriles -4.3 -3.0

Amides -7.8

∆GENP

SM5.4/P GCDS°

∆GS°

experiment ∆GS°

-3.4 -3.7 -3.5

-1.5 -2.7 -0.4

-4.9 -6.3 -3.9

-4.7 -5.5 -4.7

-6.9 -4.9

2.3 1.7

-4.5 -3.2

-3.9 -4.1

-7.4 -6.7 -7.0

3.8 4.0 3.3

-3.7 -2.7 -3.7

-3.7 -3.1 -3.6

-7.2

-1.5

-8.7

-10.0

Bifunctional N and Hydrazine -3.5 -6.9 -3.0 -4.2 -7.8 -3.5 -2.5 -5.7 -2.6 -3.2 -6.9 b

-3.9 -5.2 -3.8 b

-6.9 -8.7 -6.4 b

-6.5 -7.2 -6.3 -9.3

Thiols -1.0 -3.1

-1.0 -2.2

0.0 -0.8

-1.0 -2.9

-1.3 -2.6

2-methoxyethanamine morpholine N-methylmorpholine hydrazine

-3.3 -3.6 -3.3 -3.7

ethanethiol thiophenol

-1.3 -3.5

hydrogen sulfide dimethyl sulfide thioanisole

-1.1 -1.3 -3.9

Organic Sulfides, H2S 0.3 -0.8 -0.4 -1.7 0.7 -3.3

-0.3 -1.4 -3.0

-0.5 -0.3 -0.3

-0.8 -1.6 -3.3

-0.7 -1.5 -2.7

dimethyl disulfide diethyl disulfide

-2.4 -2.3

Organic Disulfides 0.4 -2.0 0.9 -1.4

-2.8 -2.9

1.0 1.2

-1.8 -1.7

-1.8 -1.6

fluoromethane 1,1-difluoroethane fluorobenzene

-1.8 -2.8 -3.1

Fluorinated Hydrocarbons 2.5 0.7 -1.6 2.8 0.0 -2.6 2.9 -0.2 -2.3

2.5 2.7 1.6

0.9 0.0 -0.6

-0.2 -0.1 -0.8

dichloromethane 1,1,2-trichloroethane 2-chloropropane

-1.8 -1.8 -1.5

Chloroalkanes 0.1 -1.6 -0.1 -1.9 1.1 -0.4

-1.8 -1.9 -1.5

0.2 0.0 1.1

-1.6 -1.9 -0.4

-1.4 -2.0 -0.3

3-chloropropene trichloroethene

-1.5 -0.7

Chloroalkenes 1.7 0.2 0.3 -0.4

-1.3 -0.7

1.3 0.2

-0.1 -0.5

-0.6 -0.4

chlorobenzene p-dichlorobenzene

-2.8 -2.0

Chloroarenes 1.6 -1.2 1.1 -1.0

-1.7 -1.3

0.5 0.2

-1.3 -1.2

-1.1 -1.0

bromomethane 2-bromopropane bromobenzene

-1.5 -1.6 -2.8

Brominated Hydrocarbons 0.5 -1.0 -1.3 0.7 -0.9 -1.6 1.1 -1.7 -1.8

0.5 0.7 0.0

-0.8 -1.0 -1.8

-0.8 -0.5 -1.5

diiodomethane 1-iodobutane iodobenzene

-0.7 -0.9 -2.6

Iodinated Hydrocarbons -1.3 -1.9 0.4 -0.5 0.9 -1.6

-0.9 -0.7 -1.6

-1.3 0.3 -0.3

-2.2 -0.4 -1.9

-2.5 -0.3 -1.7

1-bromo-1-chloro-2,2,2-trifluoroethane 1-bromo-1,2,2,2-tetrafluoroethane 1-chloro-2,2,2-trifluoroethane 1,1,1-trifluoropropan-2-ol bis(2-chloroethyl)sulfide 2,2,2-trifluorethyl vinyl ether p-bromophenol

-1.4 -1.8 -2.8 -6.0 -2.0 -3.7 -6.5

Other Halo Molecules 1.8 0.4 2.8 1.0 2.7 -0.1 1.4 -4.6 -0.8 -2.9 3.8 0.0 -0.7 -7.2

-1.7 -1.9 -2.8 -5.6 -2.3 -3.2 -5.0

1.7 2.8 2.6 1.1 -0.7 3.2 -2.0

0.1 0.8 -0.2 -4.5 -3.0 0.1 -7.0

-0.1 0.5 0.1 -4.2 -3.9 -0.1 -7.1

0.3 0.4

a Sources of experimental data and the order of preference used when more than one source contained a free energy of solvation for a single solute are discussed in section 2. b Hydrazine was not calculated using SM5.4/P (see text).

SM5.4/A predicts free energies of solvation for the ions CO32- and HCO2- to be -273.8 and -76.9 kcal, respectively, whereas SM5.4/P predicts -275.0 and -79.2 kcal, respectively. For comparison, we note that the values presented by Marcus73 are -314 and -94 kcal, respectively.

Table 8 contains SM5.4/A and SM5.4/P free energies of solvation for imidazole and N,N-dimethylacetamide. These two molecules are in the meta set, but not in the SM5.4 aqueous models neutral training set. Although no experimental values are available for these molecules, Pearson52 gives an “estimated”

16396 J. Phys. Chem., Vol. 100, No. 40, 1996

Chambers et al.

TABLE 7: Calculated and Experimentala Free Energies of Solvation (kcal/mol) for Ionic Solutes Used in the Parametrization of SM5.4 SM5.4/A HC2CH3OH2+ (CH3)2OH+ CH3C(OH)CH3+ H3O+ CH3OCH3CO2CH3COCH2C6H5OC6H5CH2OHHO2O2CH3NH3+ CH3C(OH)NH2+ (CH3)2NH2+ (CH3)3NH+ C5H5NH+ NH4+ CNCH2CNNH2NO2NO3N3CH3SH2+ (CH3)2SH+ HSn-C3H7SC6H5SFClBrIa

SM5.4/P

∆GENP

GCDS°

∆GS°

∆GENP

GCDS°

∆GS°

experiment ∆GS°

-80.8 -93.4 -72.7 -67.4 -117.6 -82.0 -77.5 -72.1 -66.8 -59.8 -104.3 -94.5 -88.0 -74.4 -68.2 -65.9 -58.8 -56.4 -84.3 -78.2 -67.4 -93.0 -75.3 -68.9 -70.0 -67.7 -60.2 -84.4 -75.5 -67.5 -109.3 -76.9 -70.9 -61.6

3.6 2.4 2.4 1.7 2.7 -2.6 0.9 0.4 0.1 2.9 -5.6 -4.9 -8.7 -1.6 -0.7 0.3 1.7 1.3 -3.7 6.3 2.5 -1.7 5.5 8.6 0.7 2.5 1.7 -3.1 -2.5 -1.8 2.4 -0.3 -1.1 -1.4

-77.2 -91.0 -70.3 -65.7 -114.9 -84.5 -76.6 -71.7 -66.7 -56.8 -109.9 -99.5 -96.7 -76.0 -68.9 -65.6 -57.2 -55.1 -88.0 -71.9 -64.9 -94.7 -69.8 -60.3 -69.4 -65.2 -58.5 -87.5 -78.1 -69.3 -106.9 -77.3 -72.0 -63.0

-82.0 -92.5 -73.3 -67.7 -113.6 -83.6 -78.9 -72.9 -67.1 -59.1 -105.3 -92.9 -87.2 -77.2 -67.9 -69.1 -61.2 -57.4 -84.5 -78.4 -69.0 -87.9 -76.5 -69.4 -66.1 -67.9 -60.3 -83.7 -76.7 -69.6 -109.3 -76.9 -70.9 -61.6

2.8 2.8 2.2 2.1 4.2 -4.6 0.2 -0.1 -1.2 1.5 -6.5 -8.0 -12.7 -1.2 0.1 0.5 1.6 0.4 -2.7 6.9 2.7 -0.2 4.5 6.1 3.3 1.8 1.5 -3.0 -2.1 -2.5 2.3 -0.3 -1.1 -1.6

-79.2 -89.7 -71.1 -65.7 -109.4 -88.2 -78.6 -73.0 -68.3 -57.6 -111.8 -101.0 -99.9 -78.4 -67.8 -68.7 -59.5 -57.0 -87.1 -71.6 -66.3 -88.1 -72.0 -63.2 -62.7 -66.2 -58.8 -86.8 -78.8 -72.1 -106.9 -77.2 -72.0 -63.2

-73 -85 -70 -64 -104 -95 -77 -81 -72 -59 -106 -101 -87 -70 -66 -63 -59 -59 -79 -77 -75 -95 -72 -65 -74 -74 -61 -76 -76 -67 -107 -77 -72 -63

Reference 48.

TABLE 8: Calculated and Experimental Free Energies of Solvation, ∆GS°, (kcal/mol), for Selected Solutes Not Used in the Parametrizations of SM5.4 3-ethyl-2-methoxypyrazine 1,2-dichloroethane 1,2-dibromoethane 2,2-dichloro-1,1-difluoroethyl methyl ether triethylamine trimethyl orthoacetate trimethyl orthobenzoate imidazole N,N-dimethylacetamide HCO2CO32H3S+

number of conformations

SM5.4/A

SM5.4/P

3 2 2 5 8 4 2 0 0 0 0 0

-7.0 -1.8 -2.7 -1.2 -1.1 -2.9 -2.4 -11.4 -5.8 -76.9 -273.8 -74.8

-7.7 -1.9 -2.6 -2.2 -1.9 -5.2 -2.9 -13.5 -8.2 -79.2 -275.0 -77.2

experimenta -4.4 -1.7 -2.1 -1.1 -3.0 -4.4 -4.0 -94b -314b -87c

a Sources of experimental data (for the neutral solutes) and the order of preference used when more than one source contained a free energy of solvation for a single solute are discussed in Section 2. b Reference 73. c Reference 52. Uncertainty in the pKa (upon which the experimental ∆GS° is based) led us to not include H3S+ in the SM5.4 aqueous models ionic training set.

free energy for imidazole of -6 kcal/mol (as usual, values are converted to our standard state); he references Hine and Mookerjee,47 but he does not say which of their several schemes he used. Actually, ∆GS° for imidazole can be estimated in more than one way using the theoretical schemes presented by Hine and Mookerjee.47 Using their bond order scheme and assuming that the imidazole bond contributions are 3 [Car-H], [CardCar], 2 [CardNar], 2 [C-N], and [N-H], the free energy of aqueous solvation for imidazole is -9.7 kcal/mol. Yet, the carbonnitrogen bonds in imidazole are not of purely single- or doublebond character, and they are all members of an aromatic ring. Thus imidazole could also be calculated using the following scheme: 3 [Car-H], [CardCar], 4 [CardNar], and [N-H]. This

predicts a free energy of aqueous solvation for imidazole of -10.5 kcal/mol. Of these values, we consider that -10.5 kcal may be the most reasonable, or perhaps the average of the SM5.4 calculations should be the preferred estimate. The remaining neutral results in Table 8 are for several molecules in the meta set for which we have experimental data, but which were not included in the SM5.4 aqueous models neutral training set because they would require conformational averaging, as discussed earlier. 4.6. AM1-SM5.4/U and PM3-SM5.4/U Aqueous Models. Although in general, we recommend using SM5.4/A and SM5.4/P over AM1-SM5.4/U and PM3-SM5.4/U, we believe that there are occasions when one might wish to use AM1-

Model for Aqueous Solvation SM5.4/U or PM3-SM5.4/U. Therefore, we present parameters and results for these model that use a single set of microscopic surface tensions that are equally applicable to either Hamiltonian, SM5.4/U. The microscopic surface tensions for SM5.4/U and the free energies of aqueous solvation for all molecules in the SM5.4 aqueous models neutral training set for AM1SM5.4/U and PM3-SM5.4/U are given in the supporting information. The only differences between AM1-SM5.4/U and SM5.4/A or between PM3-SM5.4/U and SM5.4/P are changes in the microscopic surface tension coefficients, many of which differences are quite small. Only six of the microscopic surface tension coefficients change by more than 20%, and only one of them changes by more than 35%. The mean unsigned errors for AM1-SM5.4/U and PM3SM5.4/U predictions are 0.56 and 0.52 kcal/mol, respectively, the mean signed errors are -0.16 and 0.23 kcal/mol, respectively, and the rms errors are 0.77 and 0.71 kcal/mol, respectively. The improvements in the SM5.4/A and SM5.4/P models over SM5.4/U are mainly in the alkynes, arenes, aromatic amines, nitrohydrocarbons, and polyfunctional aromatics. 5. Concluding Remarks The SMx solvation models have proven to be very robust and useful methods for calculating free energies of solvation in aqueous media. These models are based on estimating electrostatic effects by the generalized Born model8 with molecularshape-dependent dielectric screening and descreening effects and estimating first solvation shell effects from neighborhooddependent atomic surface tensions combined with solventaccessible surface areas for every atom in the solute. Parameters are obtained by empirical fitting to experimental data. A critical element in the SMx models is the incorporation of the electrostatic effects into a Fock operator through the use of a self-consistent reaction field. This allows the calculation of electronic structures optimized in the presence of solvents and the calculation of molecular geometries including the effects of solvent geometry-dependent electronic reorganization. Although the SMx models have been quite successful, we recognized several areas for improvement in previous models, and in this paper we presented new models, called SM5.4, that improve the following aspects of previous models. (1) The previous aqueous solvation models were based on semiempirical partial charges that typically underestimated bond polarity. Although the resulting effect on predicted free energies of solvation is apparently minimized by the semiempirical nature of the parametrization, the cancellation of systematic errors in the electrostatics by systematic compensation in surface tension terms makes the SM1-SM3.1 models less physical than SM4 and SM5.4 and makes interpretation of factors contributing to solvation energies much less reliable than the solvation energies themselves. (2) In previous parametrizations through SM4, we employed Coulomb radii and surface tensions that depend on partial charges and bond orders of the atoms involved. While this provides a very physical way to take into account the neighborhoods and functionalities of the atoms, it also has at least two disadvantages: (i) it has the potential to misclassify the functionality; that is, the results of unpublished attempts to extend the theory occasionally exhibited dependencies on partial charges or bond orders that do not appear to be present in the experimental data; and (ii) such dependencies introduce computational complications into the Fock operator and into algorithms for optimizing geometries. In the present article we introduced a new solvation model that removes these and other disadvantages. In particular,

J. Phys. Chem., Vol. 100, No. 40, 1996 16397 we use the CM1A and CM1P class IV charge models, which are characterized by a remarkable combination of low computational cost and high accuracy. Furthermore all charge and bond-order dependencies of the parameters are replaced by purely geometrical dependencies, which are both more stable and more computationally convenient. In addition, the parameters are chosen such that the calculated separation of first solvation shell effects from bulk electrostatic effects is not counterintuitive. The new models that we present, SM5.4/A, SM5.4/P, AM1SM5.4/U, and PM3-SM5.4/U, have electrostatic and nonelectrostatic components that are more physical than our previous models, can easily be extended to other solvents, have been parametrized on the basis of a large (215) set of molecules, and have rms errors lower than our previous models. Acknowledgment. The authors are grateful to David Giesen, Xavier Assfeld, and Michael Gu for many helpful discussions and for other parametrization efforts in our group that contributed in parallel to the success of the effort presented here. We also acknowledge Joey W. Storer for work on an abandoned SM4 water model and Ivan Rossi for work on an abandoned genetic algorithm parametrization strategy; both of these projects contributed to our appreciation of the important factors in the parametrization and hence indirectly to the final SM5.4 model. The authors are grateful to Thomas Halgren for an informative discussion of the experimental data for ions, which led to Table 1. This work was supported in part by the National Science Foundation under Grant No. CHE94-23927 and by the Army High Performance Computing Research Center under the auspices of the Department of the Army, Army Research Laboratory. The content does not necessarily reflect the position or the policy of the government and no official endorsement should be inferred. Supporting Information Available: Tables of rms errors over functional groups for SM5.4/A and SM5.4/P and microscopic surface tension coefficients for the SM5.4/U model as well as full tables of SM5.4/A and SM5.4/P results (∆GENP, GCDS°, and ∆GS°), AM1-SM5.4/U and PM3-SM5.4/U results (∆GS°), and experimental results (∆GS°) for all solutes in the SM5.4 aqueous models neutral training set (20 pages). Ordering information is given on any current masthead page. References and Notes (1) Grant, D. J. W.; Higuchi, T. Solubility BehaVior of Organic Compounds; John Wiley & Sons: New York, 1990. (2) (a) Leo, A.; Hansch, C.; Elkins, D. Chem. ReV. 1971, 71, 525. (b) Leo, A. J. Chem. ReV. 1993, 93, 1281. (3) Cramer, C. J.; Truhlar, D. G. J. Am. Chem. Soc. 1991, 113, 8305. Erratum: 1991, 113, 9901. (4) Cramer, C. J.; Truhlar, D. G. Science 1992, 256, 213. (5) Cramer, C. J.; Truhlar, D. G. J. Comput. Chem. 1992, 13, 1089. (6) Cramer, C. J.; Truhlar, D. G. J. Comput.-Aided Mol. Design 1992, 6, 629. (7) Liotard, D. A.; Hawkins, G. D.; Lynch, G. C.; Cramer, C. J.; Truhlar, D. G. J. Comput. Chem. 1995, 16, 442. (8) Giesen, D. J.; Storer, J. W.; Cramer, C. J.; Truhlar, D. G. J. Am. Chem. Soc. 1995, 117, 1057. (9) Giesen, D. J.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem. 1995, 99, 7137. (10) Storer, J. W.; Giesen, D. J.; Hawkins, G. D.; Lynch, G. C.; Cramer, C. J.; Truhlar, D. G.; Liotard, D. A. In Structure, Energetics, and ReactiVity in Aqueous Solution: Characterization of Chemical and Biological Systems; Cramer, C. J., Truhlar, D. G., Eds.; American Chemical Society: Washington, DC, 1994; p 24. (11) Barrows, S. E.; Dulles, F. J.; Cramer, C. J.; Truhlar, D. G.; French, A. D. Carbohydr. Res. 1995, 256, 219. (12) Cramer, C. J.; Truhlar, D. G. In ReViews in Computational Chemistry; Lipkowitz, K. B., Boyd, D. B., Eds.; VCH Publishers: New York, 1995; Vol. 6, p 1.

16398 J. Phys. Chem., Vol. 100, No. 40, 1996 (13) Cramer, C. J.; Truhlar, D. G. In SolVent Effects and Chemical ReactiVity; Tapia, O., Bertra´n, J., Eds.; Kluwer: Dordrecht, in press. (14) Tapia, O. In Quantum Theory of Chemical Reactions; Daudel, R., Pullman, A., Salem, L., Viellard, A., Eds.; Reidel: Dordrecht, 1980; Vol. II, p 25. (15) van Duijnen, P. Th. J. Chem. Soc., Faraday Trans. 1994, 90, 1611. (16) Lee, B.; Richards, F. M. J. Mol. Biol. 1971, 55, 379. (17) Hermann, R. B. J. Phys. Chem. 1972, 76, 2754. (18) Amidon, G. C.; Yalkowsky, S. H.; Anik, S. T.; Valvani, S. C. J. Phys. Chem. 1975, 879, 2239. (19) Rose, G. D.; Geselowitz, A. R.; Lesser, G. J.; Lee, R. H.; Zehfus, M. H. Science 1985, 229, 834. (20) Eisenberg, D.; McLachlan, A. D. Nature 1986, 319, 199. (21) Ooi, T.; Oobatake, M.; Nemethy, G.; Scheraga, H. A. Proc. Natl. Acad. Sci. U.S.A. 1987, 84, 3086. (22) Still, W. C.; Tempczyk, A.; Hawley, R. C.; Hendrickson, T. J. Am. Chem. Soc. 1990, 112, 6127. (23) van Gunsteren, W. F.; Luque, F. J.; Timms, D.; Torda, A. E. Annu. ReV. Biophys. Biomol. Struct. 1994, 23, 847. (24) See, for example, De Vries, A. H.; van Duijuen, P. T. Int. J. Quantum Chem. 1996, 57, 1067. (25) Hoijtink, G. J.; de Boer, E.; van der Meij, P. H.; Weijland, W. P. Recl. TraV. Chim. Pays-Bas 1956, 75, 487. (26) Peradejordi, F. Cah. Phys. 1963, 17, 393. (27) Jano, I. C. R. Acad. Sci. Paris 1965, 261, 103. (28) Bucher, M.; Porter, T. L. J. Phys. Chem. 1986, 90, 3406. (29) Ehrenson, S. J. Comput. Chem. 1989, 10, 77. (30) Schaefer, M.; Froemmel, C. J. Mol. Biol. 1990, 216, 1045. (31) Hawkins, G. D.; Cramer, C. J.; Truhlar, D. G. Chem. Phys. Lett. 1995, 246, 122. (32) Mulliken, R. S. J. Chem. Phys. 1935, 3, 564. (33) Dewar, M. J. S.; Zoebisch, E. G.; Healy, E. E.; Stewart, J. J. P. J. Am. Chem. Soc. 1985, 107, 3902. (34) Dewar, M. J. S.; Zoebisch, E. G. J. Mol. Struct. (THEOCHEM) 1988, 180, 1. (35) Dewar, M. J. S.; Yate-Ching, Y. Inorg. Chem. 1990, 29, 3881. (36) Stewart, J. J. P. J. Comput. Chem. 1989, 10, 221. (37) Storer, J. W.; Giesen, D. J.; Cramer, C. J.; Truhlar, D. G. J. Comput.Aided Mol. Des. 1995, 9, 87. (38) (a) Chirlian, L. E.; Francl, M. M. J. Comput. Chem. 1987, 8, 894. (b) Breneman, C. M.; Wiberg, K. B. J. Comput. Chem. 1990, 11, 361. (39) (a) Singh, U. C.; Kollman, P. A. J. Comput. Chem. 1984, 5, 129. (b) Besler, B. H.; Merz, K. M., Jr.; Kollman, P. A. J. Comput. Chem. 1990, 11, 431. (40) Thole, B. T.; van Duijuen, P. T. Theor. Chim. Acta 1983, 63, 209. (41) Armstrong, D. R.; Perkins, P. G.; Stewart, J. J. P. J. Chem. Soc., Dalton Trans. 1973, 838. (42) Cramer, C. J.; Hawkins, G. D.; Lynch, G. C., Giesen, D. J.; Rossi, I.; Storer, J. W.; Truhlar, D. G.; Liotard, D. A. QCPE Bull. 1995, 15, 41. (43) Masterfile (1994) from MedChem Software, BioByte Corp., P.O. 517, Claremont, CA 91711-0157. (44) Suleiman, D.; Eckert, C. A. J. Chem. Eng. Data 1994, 39, 692.

Chambers et al. (45) Cabani, S.; Paolo, G.; Mollica, V.; Lepori, L. J. Solution Chem. 1981, 10, 563. (46) Abraham, M. H.; Whiting, G. S. J. Chem. Soc., Perkin Trans. 1990, 2, 291. (47) Hine J.; Mookerjee, P. K. J. Org. Chem. 1975, 40, 292. (48) Wagman, D. D. J. Phys. Chem. Ref. Data 1982, 11, Suppl. 2. (49) Wolfenden, R. Biochemistry 1978, 17, 201. (50) Wauchope, R. D.; Haque, R. Can. J. Chem. 1972, 50, 133. (51) Han, P.; Bartels, D. M. J. Phys. Chem. 1990, 94, 7294. (52) Pearson, R. G. J. Am. Chem. Soc. 1986, 108, 6109. (53) Ford, G. P.; Wang, B. J. Am. Chem. Soc. 1992, 114, 10563. (54) Trasatti, S. Pure Appl. Chem. 1986, 58, 955. (55) Klots, C. E. J. Phys. Chem. 1981, 85, 3585. (56) Ben-Naim, A. Statistical Thermodynamics for Chemists and Biochemists; Plenum: New York, 1992; p 421. (57) Eliel, E. L.; Wilen, S. H. Stereochemistry of Organic Compounds; John Wiley and Sons: New York, 1994. (58) Tun˜o´n, I.; Silla, E.; Pascual-Ahuir, J.-L. Chem. Phys. Lett. 1993, 203, 289. (59) Cramer, C. J. J. Org. Chem. 1992, 57, 7034. (60) (a) Cramer, C. J.; Truhlar, D. G. J. Am. Chem. Soc. 1993, 115, 5745. (b) Cramer, C. J.; Truhlar, D. G. J. Am. Chem. Soc. 1993, 115, 8810. (c) Cramer, C. J.; Truhlar, D. G. J. Am. Chem. Soc. 1994, 116, 3892. (61) Ferguson, D. M.; Gould, I. R.; Glauser, W. A.; Schroeder, S.; Kollman, P. A. J. Comput. Chem. 1992, 13, 525. (62) Csonka, G. I. J. Comput. Chem. 1993, 14, 895. (63) Dobado, J. A.; Molina, J. M.; Rodriguez Espinosa, M. J. Mol. Struct. (THEOCHEM) 1994, 303, 205. (64) Kohata, K.; Fukuyama, T.; Kuchitsu, K. J. Phys. Chem. 1982, 86, 602. (65) Bondi, A. J. Phys. Chem. 1964, 68, 441. (66) (a) Jorgensen, W. L.; Buckner, J. K.; Huston, S. E.; Rossky, P. J. J. Am. Chem. Soc. 1987, 109, 1981. (b) Jorgensen, W. L. Computer Simulation of Biomolecular Systems; van Gunsteren, W. F., Weiner, P. K., Eds.; ESCOM Science Publishers: Leiden, 1989. (c) Jorgensen, W. L.; Severance, D. L. J. Am. Chem. Soc. 1990, 112, 4768. (67) Schraudolph, N. N.; Grefenstette, J. J. GAucsd 1.4, 1992. (68) Cramer, C. J.; Truhlar, D. G. In QuantitatiVe Treatments of Solute/ SolVent Interactions (Theoretical and Computational Chemistry, Vol. 2); Politzer, P., Murray, J. S., Eds.; Elsevier: Amsterdam, 1994; p 9. (69) (a) Luque, F. J.; Negre, M. J.; Orozco, M. J. Phys. Chem. 1993, 97, 4386. (b) Bachs, M.; Luque, F. J.; Orozco, M. J. Comput. Chem. 1994, 15, 446. (c) Orozco, M.; Bachs, M.; Luque, F. J. J. Comput. Chem. 1995, 16, 563. (70) Wallqvist, A.; Berne, B J. J. Phys. Chem. 1995, 99, 2885. (71) Rick, S. W.; Berne, B J. J. Am. Chem. Soc. 1994, 116, 3949. (72) Gerber, P. R.; Mu¨ller, K. J. J. Comput.-Aided Mol. Des. 1995, 9, 251. (73) Marcus, Y. Biophys. Chem. 1994, 51, 111.

JP9610776