J. Phys. Chem. 1993,97, 10175-10185
10175
Solvation Thermodynamics of Polar Molecules in Aqueous Solution by the XRISM Method Pi1 H. Lee’ and Gerald M. Maggiora Computational Chemistry, Upjohn Laboratories, 301 Henrietta Street, Kalamazoo, Michigan 49007-4940 Received: December 17, 1992; In Final Form: April 12, 1993’
The aqueous solvation free energies and the average solute-solvent interaction energies of a diverse set of organic molecules were calculated by the XRISM (extended reference interaction-site model) method for two related potential-energy functions based upon the optimized potentials for liquid simulation (OPLS)parameter set. The results are compared with available data obtained from experimental and other theoretical methods. The XRISM method with the given potential-energy function parameter sets produces reasonable free-energy values for a number of molecules but not for all of the molecules studied. For some molecules, the results were also quite sensitive to the potential-energy function parameters used.
1. Introduction Most biomolecules contain polar and charged groups that can interact strongly with polar solvents such as water. Since many biological processes take place in an aqueous environment, it is important to be able to assess the role that solvent effects play on the thermodynamic properties and structure of biologically interesting molecules involved in these processes.’ Numerous theoretical studies have been carried out by a variety of methods to address this issue. For example, free-energy perturbation calculationsbased upon either Monte Carlo (MC)2,3or molecular dynamics (MD)67 simulations have been used extensively in studying the solvation properties of biological as well as organic molecules. These simulation methods, which treat water molecules explicitly, provide detailed information on the interactions among the molecules in system. The calculations, however, are very time consuming, and thus the size of systems that can be practically studied at this time is limited.8 In addition, the effect of truncating long-range Coulombic interactions on the thermodynamic properties of a system are uncertain.gJ0 In addition to simulation-basedmethods, analytical techniques can provide a computationallyless demandingalternative. Recent applications of dielectric continuum theory based on the finitedifference Poisson-Boltzmann (FDPB) method have shown that reasonable estimates of the electrostatic component of the free energy of solvation can be obtained for charged and polar molecules in water.” Two hybrid methods that employ a continuum model of the solvent within classical12and quantummechanical” frameworks show considerable promise. In the former case, Still and co-workers12developed a solvation model in which the solvent-solvent cavity and solute-solvent van der Waals terms are represented by the solvent-accessible surface area and the solute-solvent electrostatic term is represented by the generalized Born model. In the latter case, Cramer and Truhlarl3 developed a semiempirical quantum-mechanical procedure that extendsthe classical approach by replacing molecularmechanics-basedcharges with charges obtained self-consistently from a suitably modified AM1 Hamiltonian. Over the years, integral-equation theories have been used successfully in the study of pure liquids and solutions. In particular, the referenceinteraction-site model (RISM) developed by Chandler and Andersen14 has proved to be an appropriate means for treating the solution properties of a variety of nonpolar molecular systems. The XRISM method, an “extended”version of the RISM method, was developed by Rossky and co-workersll-17 for treating highly polar and charged molecular systems. The Author to whom all correspondence should be addressed. *Abstract published in Aduance ACS Abstracts. September 1, 1993.
0022-365419312097-10175$04.00/0
XRISM method was applied to polar solvents,16ionicsolutions,18J9 polar and nonpolar solutes,18and chemical reactions in aqueous solution,21-22with reasonable success. Recently, the XRISM method was applied by Yu, Pettitt, and Karplus23 in a study of the aqueous solvation of N-methylacetamide. As was shown by these authors, the XRISM method within the Gaussian fructuation (GF) approximation produces values for solvation thermodynamic properties in favorable agreement with Jorgensen’s MC3 and their own MD simulations, as well as with experimentaldata.24.2sThey also showed that the functional used to calculate the solvation thermodynamics without further approximation (“HNC approximation” in their paper) from the HNC closure provides a reasonable estimate of the relative but not the absolutesolvationthermodynamicsof N-methylacetamide. In essentially all of the studies cited above that deal with polar or ionicsolutesin water, thenatureof the potential-energy function used was critical, especially the electrostatic component. In particular, van Gunsteren51 noted that in molecular simulations the free energy of a system is rather sensitive to the interaction potential used. For example, in the relatively simple case of the free energy of solvation of methanol in water, application of parameters from established potential-energy functions yields a spread in free energies of approximately 2 kT about the experimental value. Also, the inability of currently available potential-energy functions to account for changes in charge distribution brought about by conformational changes is well illustrated by numerous studies3.5.23 on the free-energydifference between the cis and trans isomers of N-methylacetamide, where different charge distributionsare needed for each isomer to obtain agreement with experiment. Calculation of the free energy of solvation of a diverse, but related, set of molecules provides an opportunity for testing the performance of the potential-energy function in comparison with experimental data and for obtaining indications of possible improvements to the potential-energy functions. In the present work, we calculate the average solute-solvent interaction energy and the free energy of solvation of a diverse set of organic solute molecules, all of which were treated with the same potential-energyfunction by the XRISM method. We also compare the performances of the functionals to calculate the free energy of solvation derived from both the G F and HNC approximations on an extensive set of molecules. As the results of preliminary calculations were quite sensitive to the nature of the electrostatic parameters employed, two different but related potential-energy function parameter sets, OPLS and OPLSA, were investigated. The OPLS parameter set was developed by Jorgensen and Tirado-Rives,16 from simulation results of organic liquids, and the OPLSA26 set is a Q 1993 American Chemical Society
Lee and Maggiora
10176 The Journal of Physical Chemistry, Vol. 97, No. 39, 1993 combination of parameters obtained from the OPLS and AMBER2' potential-energy functions. In section 2, a brief discussion of the theoretical method for calculating solvation free energy from the radial-distribution function, g(r), is presented. This is followed by a discussion of the general features of the XRISM method and how it can be used to determine g(r). The nature of the potential-energy function parameters employed is also discussed. In section 3, the relevant radial-distribution functions are given along with the values of the solvation free energies obtained from both potentialenergy functions. A comparativediscussion of the results obtained for each potential-energy functions is given, along with a comparison of these results with data obtained by experiment and other theoretical methods. Section 4 summarizes the paper with an overview of the workand a discussionof several important, unresolved methodological issues.
2. Methodology 2.1. Free-Energy Calculations. In a canonical ensemble, the Helmholtz free energy of the system,A, is related to the partition function, Z, by the following equation -PA = In Z (1) where8 = l/kT, kistheBoltzmannconstant,and Tistheabsolute temperature. The partition function is given by
2=
s ...
sexp[-pU(r)] dr
2.2. XRISMTbeary. TheXRISM theo~y,l~~~asnotedearlier, is an extension of the RISM integral equation theory of Chandler and AndersenI4allowing it to be applied to polar and chargedmolecules. This is accomplished through a renormalization of the long-range Coulombicpotentials and appropriately modified closure relations. The theory consists of two equations in two unknowns, the total and direct correlation functions. The first equation isa site-siteomstein-Zernike (0Z)-likeequation which can be expressed in a matrix form as
h = W*C*W
(6)
where h = g - 1 represents the matrix of the total correlation function, c the matrix of the direct correlation function, w the matrix of the intramolecular correlation function, and p the diagonal number density matrix. The asterisk denotes a matrix convolution product. Each term in eq 6 is an N X N matrix where N = nu &, the sum of the numbers of interaction sites on the solute and solvent molecules, respectively. The total intermolecular potential energy is represented by a sum of spherically-symmetric,site-site pair potentials, U4(r), which can be divided into two terms,
+
-8UayW = - 8 V a y ( r ) + cb&
(7) where ub,, is a short-range potential and 44 is a long-range Coulombic interaction potential. Using this decomposition, the short-range total, hl, and direct, e,correlation functions can be expressed as
(2)
where the total potential energy of the system, U(r), is a function of the set of molecular coordinatesdenotedby r. In an interactionsite model, molecules are representedby sets of discrete interaction sites that are commonly located at the positions of the atomic nuclei. The total potential energy of a mixture of "interactionsite molecules" is then obtained as the sum over all pairwise, spherically symmetric interaction-site potentials. Thus, in a solute-solvent mixture at infinite dilution, the total potential energy of the system, assuming pairwise additivity of the site-site intermolecular interaction potential, is given by
(3) where the coupling parameter,28 A, is defined over the interval [0,1] such that X = 0 for pure solvent and X = 1 for the solutesolvent mixture. The index u refers to sites on the solute molecule and the indices u and u' to sites on the solvent molecules, respectively, and the sums are over all distinct solute-solvent sitesite pairs. The excess free energy of the system relative to pure solvent, AA, can then be expressed as an integral over the coupling parameter,
cs=c-4
Q = Ch [ ~ I P w I (10) The second equation in XRISM theory is a closure relation which represents an approximaterelation between hand c. Rossky and co-workers adopted the site-site hypernetted chain (HNC)like closure,1s which can be expressed in component form as hSay= exP(-BU",y
+ hs,y - p a y + Q q ) - Q a y - 1
(1 1)
Equations 6 and 11 form a set of coupled equations that can be solved numerically. The matrices in eq 6 can be expressed as submatrices that represent the solutesolute (uu),solute-solvent (uu), andsolventsolvent (UU) terms explicitly. Expansion of eq 6,assuming an ideal solution limit (Le., p, 0), yields a set of three matrix equations
-
h,, = w,*c,,*w,
eq 4b
(8)
h'=h-Q (9) where the renormalized Coulomb interaction, Q, is defined by the chain
Using eqs 1-3, eq 4a can be expressed as
From the definition of the radial distribution can then be expressed as
+ w*c*ph
+ w,*cw*p4,
(1 2 4
where the subscript u has been dropped from p for convenience. Q is the result that would be obtained forb if 4 alone were inserted for c in eq 12,17thus,
Q,, = wD*4,*wD + W , * ~ ~ * P Q ,
(134
UD
where p is the solvent density and gw(r,X) is the solute-solvent radial-distribution function for a system with interaction potential U(r,X). Radial-distribution functions can be obtained experimentally, theoretically by simulation, or by integral equation methods.
Q,, = W , * ~ ~ , * W , + W,*~,*PQ,,,
(13c)
Equation 12 can now be rewritten in terms of the renormalized functions as
The Journal of Physical Chemistry, Vol. 97, No. 39, 1993 10177
Solvation Thermodynamics by the XRISM Method
TABLE I: Potential Energy Parameters for the Water Molecule!
A
atom H
0.40
0.046
49 e 0.4 17
0
3.15
0.152
4834
0,
e,
kcal/mol
a ROH 0.9572 A, LHOH = 104.52O.
following equation: The functionsQ are obtained from the intramolecular correlation functions and long-range Coulombic interactions using eq 13 and then used to solve eq 14. It is w e l l - k n ~ w nthat ~ ~ all the RISM theories yield an underestimated solvent dielectric constant due to the implicit constraint in the closure that the direct correlation function cay(r) is asymptotically equal to I$&). There are two ways to correct for the underestimated solvent dielectric constant in the HNC-RISM formulation. One way19 is to remove ~ R I S Mfrom the potential-of-mean-force(PMF) between the solute sites and replace it with the experimental value, eo. The other way22,42is to modify the HNC-like closure so that the calculated solvent correlationfunctionsare consistent with the macroscopicdielectric constant of the solvent, CO. The second option was used here, and the modified closure is now given by
+ D9,, + hay - cay) - 1
ha, = exp(-BU,,
(1 5)
where D is a site-independent constant given by
D=
1
+ e0(3y-
The free energy of solvation can be obtained from two functionals. The first functionalis derived from the HNC closure, and according to Singer and Chandler,30the coupling parameter integration in eq 5 in section 2.1 can be performed in closed form within the XRISM-HNC formalism. Their result, for infinitely dilute solute, is given by ApHNC=
P -7 7 K4.2
dr [ 1/2ha:(r)
- cay(r)-
Pa=l y=l
'/2hay(r) cay(r)I (19) where A p = limp,+ (AAIN,) and Nu is the number of solute molecules. The second functional was derived from the assump tion of Gaussian fluctuations (GF) for the solvent by Chandler, Singh, and Richardson,20.M
1)
3y(cO - 1) withy = 4*/3p(d,2)/9, where ( 4 2 )is the average squared dipole moment of the solvent molecule and p is the solvent density. It should be noted that the modified HNC closure was used only for the calculations of solvent-solvent correlation functions and if either CY or y is a solute site in eq 15, D = 1, which makes it the actual HNC closure. Equations 14 and 15 can be solved by standard iterative methods; an initial guess for c(r) is used in the OZ equation to obtain h(r) which is used in the closure to obtain a new estimate for c(r). This process is repeated until c(r) has changed less than a given convergence criteria (6 = 0,001). To avoid divergence, the function, &(r) is mixed with &-l(r) before the next iteration according to cimix(r)= CYci-l(r)
+ (1 - CY)&)
(17) where i indicates the ith iteration and CY is the mixing parameter. The Fourier transformswere performed with Talman's algorithm52 as implemented for spherically symmetric transforms by Rossky and Friedman" on a logarithmically spaced grid. p = In r,
pn = pm
+ nAp,
r = ep
Ap = 0.02 and pm = -5.12
where n = 0, ..., N - 1. N = 512 grid points were used in all calculations in this work. An initial guess can be obtained from functions with similar potential-energy parameters. To obtain solutions to the XRISM equation in cases where no reasonable initial guess exists, the interactions of the solute with the solvent are slowly increased from small values up to the full potential through a series of computations. For a solute-solvent mixture, the thermodynamic quantities can be expressed in terms of solute-solvent correlation functions. The average solutesolvent interaction energy is an important term in solvation thermodynamics and can be obtained from
This functional appears to give better values of the solvation thermodynamic properties for N-methylacetamidez3than that based on the HNC closure (eq 19). The numerical integration in eqs 18-20 was obtained by the Gill and Miller algorithm which is implemented in the FORTRAN program, CUBINT in ref 54. The cut-off distance in the numerical integrations in eqs 18-20 was 21.3 A. Results obtained from both the HNC and G F functionalswill be compared for the diverse set of solute molecules studied here. An advantage of the interaction-site model is that it provides a means for calculating the contribution to the total energy from each interaction site on the solute molecule. The contribution from each site to the average solute-solvent interaction energy and the free energy of solvation can be calculated from eq 18 and eqs 19 and 20, respectively. Chemical group contributions then can also be obtained by combining appropriate sites into a chemical group. This point will be discussed further in section 3. 2.3. Potential-Energy Function Parameters and Thermodynamic Conditions. A three-site model similar to the TIP3P model developed by Jorgensen31 was used for water in all calculations reported here. The geometry and potential-energy function parameters for water are given in Table I. The structures of the solute molecules studied here are depicted in Figure 1. The intermolecular interaction energy between molecules [I and b is given by the sum of interactions between the na sites on molecule a and the nb sites on molecule b, respectively,
where CY and y are sites on molecules a and b, respectively. The intersite interactions were calculated using the standard mixing rules:
Lee and Maggiora
10178 The Journal of Physical Chemistry, Vol. 97, No. 39, 1993 %lute Molecules OPLSA
OPLS
OPLS
OPLSA 12, 13, 14. Alanine Dipeptide
1. Methanol
2. Ethanol
15. Benzene ?7
Me 1 4 ° 3 - H 4
3. 2-Propanol
16. Toluene
4. Acetone 17. Phenol
5. Methyl acetate Hll
18. Ammonium
OPLS
OPLSA 6. Acetic Acid
i:
H4-i1 -H2
H5
OPLS
H5
OPLSA 19. Acetate Ion
7. N,N-Dimethyl Acetamide
20. Ammonium Amtale 8. Acetamide
9, 10. N-Methyl Acetamide
11. cis N-Methyl Acetamide Dimer
p7
H23
Figure 1. Solute molecules. The potential energy parameters for the numbered sites are given in Table 11. The structures depicted are purely for the purpose of identifying individual atomic sites and do not reflect the actual conformationsused in the calculations.
Solvation Thermodynamics by the XRISM Method
The Journal of Physical Chemistry, Vol. 97, No. 39, 1993 10179
obtained for the set of solute molecules depicted in Figure 1 are compared in Table 111 with the relative percent error due to differences in the potential-energy functions. As is seen in the table, both parameter sets yield very similar values for (U,) in a majority of the moleculesstudiedhere. However, for 2-propano1, aceticacid, NJv-dimethylacetamide, benzene, toluene,and phenol, The OPLS and OPLSA parameter sets are listed in Table 11.The the differences between the results obtained from the OPLS and main differenca between the two potentials are in the CH,group. OPLSA parameter sets are rather large (relative % error 1 10%). In OPLS, these groups are treated as united atoms centered on For benzene, toluene, and phenol, the OPLS and OPLSA the carbon atom, while in OPLSA, each atom in the group is potential-energy parameters for the aromatic ring are quite treated as a separate interaction site. The parameters for OPLSA different, especially the partial charges on H and C (see Table are taken from MacroModeP developed in Clark Still’s lab at 11). Thus, the observed differences in (V,) are not surprising. Columbia University. It should be noted that the parameters of To understand the difference in (U,) for the other molecules, all hydrogen atoms bonded to nitrogen or oxygen atoms were appropriately modified to avoid the Coulombic ~ i n g u l a r i t y . ~ ~ ~due ~ ~ to the two parameter sets, it is useful to compare the contributionsto ( U,) from individualsites of the solute molecules. The geometries used in all calculations were obtained by In Table IV, the contributions from each site of several solute optimizing the structures in the gas phase. The conformation of molecules calculated with the two parameter sets are compared. the alanine dipeptide was determined by two dihedral angles: It is seen that both parameter sets produce contributions from 4(Cl-N&&) and $(NfC&-N9) (see Figure 1). The C7, polar sites that are generally large compared to the contributions C5,andaconformationshave (4,$) valuesof (82’,-67’), (177’,from nonpolar sites. In addition, contributions to (U,) from 178’), and (67’,33O), respectively. positively charged sites are positive and those from negatively All calculations were carried out at a temperature at 298.15 charged sites are negative. Significant variations also exist K and a pressure of 1 atm. The density of water under these between values obtained by the OPLSA and OPLS parameter conditions is p = 0.997 g/cm3, yielding 0.03334 molecules perA3 sets. For 2-propanol, it is obvious that the difference between the anda dielectricconstantof g = 78.4. Toobtain solvent correlation results from the two parameter sets is due to the contribution functions that are consistent with the solvent dielectric constant, from the oxygen atom. For acetic acid and N,N-dimethylaceta value of 0.959 is needed for D (see eq 16 in section 2.2). amide, the difference in negative charge on the oxygen and nitrogen atoms, respectively, makes about 5 kcal/mol difference 3. Results and Mscwion in (V,). The OPLS and OPLSA parameter sets do not make any difference in ( V,) for the Cs conformation, which is linear, 3.1. SolventSolvent Radial-Distribution Functions. The but do make a difference for the C7 and a conformationsof about solventdolvent radial-distribution functionsfor liquid water have 5 kcal/mol. For acetone and methyl acetate, the differences in been calculated earlier by the XRISM method using potentialthe contributions from each site obtained with the two parameter energy function parameters similar16 and identical32to those used sets cancel out. in the present study. The results obtained here are given in Figure 2. When compared with experimental data,33-35the locations of TheresultsinTables I11 and IValsoindicate that thedifferences the first and second neighbor peaks of &H, 2.32 and 3.9 A, are in ( U,) between cis- and trans-N-methylacetamide are almost in excellent agreement with the experimental values, 2.29 and the same regardless of the parameter set used, while in the case 3.9 A, respectively. For goo, the location of the first peak, 2.9 of the alanine dipeptide, the relative (U,) among different A, is in excellent agreement with the experimental value of 2.9 conformations depends upon the parameter set. A, but the location of the second peak, 4.22 A, is 0.28 A closer 3.3. Free Energy of Solvatim. Solvationfree energies at infinite to the origin, compared to the experimental value. The first and dilution calculated within the G F and HNC approximationsfrom second peaks of gOH are located at 1.49 and 3.53 A, respectively, the OPLS and OPLSA parameter sets are compared with and are 0.37 and 0.33 A closer to the origin than the experimental experimentaldata3639inTable V. As has been observed by other values. investigators, the G F approximation is superior to the HNC The number of atomic sites surrounding the oxygen atom (i.e., approximationfor predicting the solvationfree energy of nonpolar its coordination number) provides a useful means for evaluating solutes,Nmonovalent cations,41and N-methylacetamide,” where the accuracy of the solvent-solvent radial distribution functions both nonpolar and polar interactions are present, but not for and can be obtained from the following equation anions.23 In the present work, the G F approximation (eq 20) is seen to N&) = ~ ~ J , ‘ d r ’ 4 a r ’ ~ g & 3 provide considerably better estimates of the solvation free energy than the HNC approximation (eq 19) for every polar solute where y is either H or 0 and r is the distance to any minimum molecule studied here except methanol and acetate ion. Within of g*(r). The coordination number of H about 0 within the theGF approximation the resultsobtained vary with the potentialfirst coordinationshell, NOH(1.97 A) = 0.47, and within the second energy function used, although in most cases the values are quite shell, No~(4.57A) = 12.79, and the coordination number of 0 similar except for 2-propanol, cis-N-methylacetamide dimer, aboutOwithinthefirstshell,N~(3.82A)=7.41,wereobtained. benzene, toluene, and the C5 conformer of the alanine dipeptide, Even though the radial-distribution functions calculated here do where the differences are quite large ( 1 4 kcal/mol). In some not quantitativelyreproduce experimentaP3-3sresults, they exhibit cases, one parameter set produced values clearly in better qualitative features of the known water structure. agreement with experiment. For example, OPLS-based results 3.2. Average Soluted3olventInteractionEnergy. The average for 2-propanol and acetamide, are in better agreement with solute-solvent interaction energy, (U,),plays an important role experiment than OPLSA-based results, while the situation is in solvationthermodynamics and can be obtained from the solutereversed for acetic acid and N,N-dimethylacetamide. solvent interaction energy and the solute-solvent radial-distriAs discussed earlier for ( U,), the differences in calculated bution functions by eq 18 in section 2.2. Both the HNC and the solvation free energies can be attributed to, among other things, G F approximations give the same values of the averageinteraction differences in the OPLS and OPLSA potential-energy functions. energy. Yu and co-workers23 showed that XRISM calculations As seen in Table 11, the main difference between the two lies in of ( V,) for trans- and cis-N-methylacetamide determined with the CH. groups. In OPLS, these groups are treated as united the OPLS parameter set were in very good agreement with their atoms centered on the carbon atom, while in OPLSA, each atom own molecular dynamics results. In the present work, the results
Lee and Maggiora
10180 The Journal of Physical Chemistry, Vol. 97, No. 39, 1993
TABLE II: Potential-Energy Parameters for Solute Molecules' atom no. 1 2 3
atom type CH3 0 H
CH3 CH2
0 H
OPLS 4 0.265 -0.700 0.435
OPLSA c
Q
3.775 3.07 0.40
atom no.
atom type
4
Q
e
H C H H H
0.038 0.150 0.038 -0.700 0.038 0.435
2.50 3.80 2.50 3.07 2.50 0.40
0.05 0.05 0.05 0.17 0.05 0.046
1 2 3-6 7 8 9
C C H 0 H H
0.189 -0.115 0.038 -0.700 0.038 0.435
3.80 3.80 2.50 3.07 2.50 0.40
0.05 0.05 0.05 0.17 0.05 0.046
1 2 3 4 5 6-12
C
-0.115 0.227 -0.115 -0.700 0.435 0.038
3.80 3.80 3.80 3.07 0.40 2.50
0.05 0.05 0.05 0.17 0.046 0.05
1 2 3 4 5-10
C C 0
3.8 3.75 2.96 3.8 2.5
0.05 0.105 0.21 0.05 0.05
1. Methanol 0.207 1 0.17 2 0.046 3 4 5 6
0.000 0.265 -0,700 0.435
3.905 3.905 3.07 0.40
2. Ethanol 0.175 0.118 0.17 0.046
0.000 0.265 0.000 -0.700 0.435
3.91 3.85 3.91 3.07 0.40
0.16 0.08 0.16 0.17 0.046
0.062 0.300 -0,424 0.062
3.91 3.75 2.96 3.91
4. Acetone 0.16 0.105 0.210 0.160
0
3.2-Propanol C C
0 H H
H
-0.115 0.430 -0,430 -0.115 0.038
C C 0 0 C H
-0.115 0.600 -0.450 -0,400 0.135 0.038
3.8 3.75 2.96 3.00 3.80 2.5
0.05 0.105 0.21 0.17 0.05 0.05
C
0.05 0.55 -0.45 -0.40 0.25
3.91 3.75 2.96 3.00 3.80
5. Methyl Acetate 0.16 1 0.105 2 0.210 3 0.170 4 0.170 5 6-1 1
0.08 0.55 -0.50 -0.580 0.45
3.91 3.75 2.96 3.00 0.40
6. Acetic Acid 0.16 1 0.105 2 0.210 3 0.170 4 0.046 5 6-8
C C 0 0 H H
-0.115 0.600 -0.450 -0.585 0.435 0.038
3.8 3.75 2.96 3.00 0.40 2.5
0.05 0.105 0.21 0.17 0.046 0.05
1 2 3 4 5-6
-0.57 0.50 -0.50 0.00 0.285
3.25 3.75 2.96 3.91 3.80
7. N,N-Dimethylacetamide 0.17 1 0.105 2 0.210 3 0.160 4 0.170 5-7 8-9 10-15
N C 0 C H C H
-0.400 0.500 -0.500 -0,115 0.038 0.085 0.038
3.25 3.75 2.96 3.80 2.50 3.80 2.50
0.17 0.105 0.21 0.05 0.05 0.05 0.05
1 2 3 4 5-6
-0.85 0.50 0.00 -0.50 0.425
3.25 3.75 3.91 2.96 0.40
N C 0 C H H
-0.850 0.500 -0.500 -0.115 0.038 0.425
3.25 3.75 2.96 3.80 2.50 0.4
0.17 0.105 0.21 0.05 0.05 0.046
-0.57 0.20 0.50 0.00 -0.50 0.37
3.25 3.80 3.75 3.91 2.96 0.40
N C
0 C H H C H
-0.57 0.500 -0.500 -0.115 0.038 0.370 0.085 0.038
3.25 3.75 2.96 3.80 2.50 0.4 3.8 2.5
0.17 0.105 0.21 0.05 0.050 0.046 0.05 0.05
0.50 -0.57 0.00 -0.50 0.20 0.37
3.75 3.25 3.91 2.96 3.80 0.40
C N C 0 C H
0.500 -0.570 -0.115 -0.500 0.085 0.370
3.75 3.25 3.80 2.96 3.80 0.4
0,105 0.170 0.05 0.21 0.050 0.046
1 2 3 4 5 6
8. Acetamide 0.17 1 0.105 2 0.160 3 0.210 4 0.046 5-7 8-9 9, 10. N-Methylacetamide 0.170 1 0.170 2 0.105 3 0.160 4 0.210 5-7 0.046 8 9 10-12
11. cis-N-Methylacetamide Dimer 0.105 1 0.170 2 0.160 3 0.210 4 0.170 5 0.046 6
The Journal of Physical Chemistry, Vol. 97, No. 39, 1993 10181
Solvation Thermodynamics by the XRISM Method
TABLE II: (Continued) OPLSA
OPLS atom no. 7 8 9 10 11 12
1 2 3 4 5 6 7 8 9 10 11 12
atom type N
CH3 C H CHI, 0
C N 0
CH3 H C C CH3
N 0
H CH3
4 -0.57 0.2 0.5 0.37 0.0 -0.5
0.50 -0.57 -0.50 0.00 0.37 0.20 0.50 0.0 -0.57 -0.50 0.37 0.2
U
c
atom no.
11. cis-N-MethvlacetamideDimer 3.25 0.17 7 3.80 0.17 8 3.75 0.105 9 0.4 0.046 10 3.91 0.16 11 2.96 0.21 12 13-24 3.75 3.25 2.96 3.91 0.40 3.80 3.75 3.91 3.25 2.96 0.40 3.80
12-14. Alanine Dipeptide 0.105 1 0.170 2 0.210 3 0.160 4 0.046 5 0.080 6 0.105 7 0.16 8 0.170 9 0.210 10 0.046 11 0.17 12 13-22 15. Benzene 0.07 1-6 0.03 7-12
atom type
4
U
c
N
-0.57 0.085 0.5 0.37 -0.115 -0.5 0.038
3.25 3.8 3.75 0.4 3.8 2.96 2.5
0.17 0.05 0.105 0.046 0.05 0.21 0.05
H
0.500 -0.570 -0.500 -0.115 0.370 0.162 0.50 -0.115 -0.57 -0.5 0.37 0.085 0.038
3.75 3.25 2.96 3.80 0.40 3.8 3.75 3.80 3.25 2.96 0.40 3.8 2.5
0.105 0.170 0.21 0.05 0.046 0.05 0.105 0.05 0.17 0.21 0.046 0.05 0.05
C H
-0.15 0.15
3.75 2.5
0.1 1 0.05
C
-0.15 0.00 -0.15 0.15 -0.115 0.150 0.038
3.75 3.75 3.75 2.50 3.80 2.50 2.50
0.1 1 0.11 0.11 0.05 0.05 0.05 0.05
-0.70 0.265 -0.150 0.435 0.150
3.07 3.75 3.75 0.40 2.50
0.17 0.1 1 0.1 1 0.046 0.05
-0.32 0.33
3.25 0.40
0.17 0.046
-0.214 0.70 -0.80 0.038
3.80 3.75 2.96 2.50
0.05 0.105 0.21 0.05
-0.214 0.70 -0.80 0.33 -0.32 0.33 0.038
3.80 3.75 2.96 0.40 3.25 0.40 2.50
0.05 0.105 0.21 0.046 0.17 0.046 0.05
C C
H C 0
H C N 0
C H C C C N
0
H C
C H
-0.115 0.115
3.55 2.42
C C C H CH3
-0.115 O.OO0 -0.115 0.115 O.Oo0
3.55 3.55 3.55 2.42 3.91
1 2 3-7 8 9-13
0
-0.700 0.265 -0.115 0.435 0.1 15
3.07 3.55 3.55 0.40 2.42
17. Phenol 0.17 0.07 0.07 0.046 0.03
1 2-5
N
-0.400 0.350
3.25 0.40
18. Ammonium Ion 0.17 1 0.046 2-5
N
1 2 3-4
CHI C
-0.100 0.700 -0.800
3.91 3.75 2.96
19. Acetate Ion 0.16 1 0.105 2 0.210 3 4 5-7
C C 0
3.91 3.75 2.96 0.40 3.25 0.40
20. Ammonium Acetate 1 0.16 2 0.105 0.210 3 4 0.046 5 0.17 6 7-9 0.046 10-12
1-6 7-12 1 2 3-6 7-1 1 12
1 2 3-4 5 6 1-9
C C H H H
0
CH3
C 0
H N
H
-0.100 0.700 -0.800 0.350 -0.400 0.350
16. Toluene 0.07 1 4 0.07 5 0.07 6 0.03 7-10 0.16 11 12 13-15 1 2 3-7 8 9-1 3
C C
H C H H 0 C C
H H H
H C C 0
H N
H H
I, See Figure 1 for atom numbering. is treated as a separate interaction site. Although in OPLSA the sum of the partial atomic charges on the CH, group is the same as its total site charge in OPLS, the contribution of this group to the total solvation free energy is apt to differ for the two parameter sets-since each hydrogen atom interaction site contains a partial charge in OPLSA, a more polar interaction between the CH, group and the solvent is expected than is likely to be the case in OPLS, which has only a single interaction-site. Table VI illustrates, for four of the solutes studied here, the role that specific chemical functional groups play in determining the overall solvation free energy. Not unexpectedly, the results clearly show that it is more favorable to solvate polar groups (e.g., OH or C 4 ) in water than it is to solvate nonpolar ones
(e.g., CH2 or CH3). Moreover, in contrast to the situation for ( U w ) ,contributionsfrom CH3 groups are quite large. Generally, use of the OPLSA parameter set leads to more favorablesolvation of CH, groups than does use of the OPLS parameter set, accentuating the negative solvation free energy for molecules containing one or more of these groups. The solvation free energy cannot be rigorously separated into its electrostatic and hydrophobic contributions without making additional approximations as these two terms are correlated; solute-solvent radial distribution functions change in response to changes in molecular chargedistribution. Thus, the electrostatic contribution to the solvation free energy cannot be obtained from the free energy difference between a fully charged and a
10182 The Journal of Physical Chemistry, Vol. 97, No. 39, 1993 I
3.0-
I
I
4.0
6.0
Bw : :
2.0
0.0
8.0
r(A)
F i p e 2. Radial-distribution functions for pure water.
TABLE III: Average SoluteSolvent Interaction Energy, ( U,)’ solute
OPLS
1. methanol 2. ethanol 3. 2-propanol 4. acetone 5. methyl acetate 6. acetic acid 7. N,N-dimethylacetamide 8. acetamide 9. tram-N-methylacetamide 10. cis-N-methylacetamide 11. cis-N-methylacetamidedimer 12. ala-ala C7 13. ala-ala CS 14. ala-ala a 15. benzene 16. toluene 17. phenol 18. ammonium ion 19. acetate ion
-28.64 -30.40 -25.66 -19.04 -19.60 -29.30 -33.07 -36.72 -35.05 -29.98 -3 1.68 -58.55 -58.67 -71.83 -14.35 -13.43 -32.01 -143.2 -235.0 -109.1
20. ammonium acetate
OPLSA re1 % errorb -27.37 -29.61 -30.77 -17.30 -18.73 -25.09 -28.29 -36.42 -34.93 -29.79 -32.50 -52.81 -58.30 -67.02 -21.71 -20.89 -38.33 -143.50 -238.60 -112.10
4.4 2.6 -19.9 9.1 4.4 14.4 14.5 0.82 0.34 0.63 -2.6 9.8 0.63 6.7 -51.3 -55.5 -19.7 4.21 -1.5 -2.7
a In kcal/mol. Re1% error = [ (( UW)0prs- ( UW)0PrsA)/ ( ~ x 100.
W ) O P ~ I
hypothetically uncharged molecule, although such an approach does in many cases provide a reasonable estimate of this quantity. However, the hydrophobiccontribution to the solvation free energy can be obtained from simulations or XRISM calculations on hypothetically uncharged molecules. For example, in the case of benzene, the Lennard-Jones parameters and partial charges of the carbon and hydrogen atoms are larger in the OPLSA parameter set than in the OPLS set. If all charges are set to zero in both parameter sets, values for the “hydrophobic” free energy of 5.01 and 2.1 1kcal/mol are obtained for the OPLS and OPLSA parameter sets, respectively. This shows that the hydrophobic contribution to the free energy of solvation is also sensitive to the potential-energy parameter for this molecule. 3.4. SoluteSolvent Radial-Distribution Function. Although solvation thermodynamic propertiesdepend, among other things, upon the solute-solvent radial-distribution functions, gw(r),it is not possible in the present work to examine these functions for all of the molecules studied here. Nevertheless, it is instructive to consider, if only briefly, the effect that the environment surrounding a particular interaction site in a given molecule has on that site’s solute-solvent radial-distribution functions. Stated simply, the nature of the molecular framework surrounding the interaction site influences, to varying degrees, the ways in which solvent molecules can be distributed about that site. And the degree to which this environmental influence is felt relates strongly to the “transferability” of g,, values over a set of molecules containing a particular type of interaction site. For example, the oxygen atom in methanol, ethanol, and
Lee and Maggiora 2-propanol is located in slightly different environments for each molecule, although the potential-energy parameters describing the site are identical. Figures 3 and 4 depict the radial-distribution functions, goo, and goH, for the water 0 and H atoms about the 0 in the three alcohols, calculated with the OPLS and the OPLSA parameter sets. The OPLS-based results for goH show a strong peak at 1.46 A for all three alcohols, indicative of strong hydrogen bonding. However, the height of this peak differs in each case, being 3.25,3.08, and 2.59 for methanol, ethanol, and 2-propanol, respectively. The second maximum, due to the other H atom of the water molecule, is found a t 3.67, 3.67, and 3.6 A, with respective peak heights of 1.27, 1.15, and 1.04 for the three alcohols. For goo, the first maximum occurs at 2.95 A, with peak heights of 2.06, 1.94, and 1.86 for methanol, ethanol, and 2-propanol, respectively, while the location of the second peak varies from 4.57 A in methanol, to 4.76 A in ethanol, to 5.59 A in 2-propanol with corresponding peak heights of 1.13,1.05, and 1.13. Similar results are obtained with the OPLSA parameter set. However, the gOH for 2-propanol obtained from the OPLSA parameter set shows a stronger hydrogen bond between solute 0 and solvent H than is obtained for the OPLS parameter set, and this difference is reflected in the contribution from this site to (UlIu). Taken together, the solute-solvent distribution function results obtained with both parameter sets clearly show that the “molecular environment” surrounding the hydroxyhxygen interaction site influences the manner in which solvent (Le., water) molecules can distribute themselves about the interaction site. In addition to the effects produced by the “steric factors” considered here, electrostatic factors can also play a significant role: a detailed discussion of these issues will be presented in a subsequent publication. 3.5. Comparison with Results from Other Studies. A number of different techniques have been employed in the calculation of aqueous solvation free energies. Of the techniques in use today, MC- and MD-based simulation methods are the most popular. Comparison of the results produced by these methods with the results obtained here is useful but, in some cases, is difficult due to several factors. For example, in MC- or MD-based simulations, the solvation free-energy difference between two molecular species, and not the solvation free energy of the individual species, is generally calculated. Hence, direct comparison with the results of the XRISM calculations given here, which yield solvation free energies of individual molecular species, is not possible. Nevertheless, although it is generally less accurate, comparisons can be made by taking the difference between the calculated XRISM values, and this is the approach followed here. In addition, as also seen here, the nature of the potential-energy function employed in the calculations can have a dramatic effect on the results obtainedsmall differences in the potential-energy function can sometimeslead to rather large differences in the calculated values. And finally, the thermodynamic conditions under which the calculations were carried out can, in certain cases, also lead to significant differences. With these limitations in mind, the following brief comparisons will show that the XRISM method produces reasonable values of solvation free energy relative to a number of other methods. The alanine dipeptide has been studied extensively over the years by various m e t h o d ~ . l ” l ~ The - ~ ~ most ~ recent MD simulation studies by Tobias and Brooks7 of the solvation freeenergy differences between the cy and C7 conformations yielded a value of 8.7 kcal/mol (we compare their result for aLwith our result for cy conformation). The XRISM results obtained here for the OPLS and OPLSA parameter sets are, respectively, 6.88 and 7.86 kcal/mol, both of which are in reasonable agreement with the simulation results considering that different potentialenergy functions were used in two studies. They also obtained 19.2 kcal/mol for the relative (U,,”) difference for a C7
-
Solvation Thermodynamics by the XRISM Method
The Journal of Physical Chemistry, Vol. 97, No. 39, 1993 10183
TABLE I V Contributions to ( V,) from Individual Sites of Solute MolecuIes*b
1. CH3 2. CH 3. CH3 4.0 5. H total
OPLS 2-Propanol -1.79 2.81 -1.53 -27.24 2.10 -25.66
OPLSA
diffc
-2.21 4.12 -1.91 -32.52 1.75 -30.77
0.42 -1.31 0.38 5.28 0.35 5.1 1
1. CH3 2. c 3.0 4. CH3 total
-0.97 2.90 -20.01 -0.97 -19.04
-1.68 5.88 -19.81 -1.68 -17.30
0.71 -2.98
1. CH3 2. c 3.0 4.0 5. CH3 total
Methyl Acetate -0.72 -1.74 9.41 11.83 -18.81 -18.88 -10.98 -1 1.70 1S O 1.75 -19.60 -18.74
1-6. c 7-12. H total
-3.98 1.59 -14.35
-5.77 2.15 -21.71
1.9 -0.54 7.36
1. CH3 2. c 3.0 4.0 5. H
Acetic Acid -0.95 7.85 -20.67 -8.56 -6.97
-1.97 8.37 -15.56 -9.97 -5.95
1.02 -0.52 -5.1 1 1.41 -1.02
-29.30 -25.09 N,N-Dimethylacetamide -16.08 -10.98 14.73 14.12 -31.87 -31.21 -1.42 -1.87 2.17 1.91 -0.59 -0.26 -33.07 -28.29
-4.21
OPLS
Acetone
total 1. N 2. c 3.0 4. CH3 5. CH3 6. CH3 total
-0.2 0.71 -1.75 1.02 -2.42 0.07 0.72 -0.25 -0.86
Benzene
-5.1 0.61 -0.66 0.45 0.26 -0.33 -4.78
OPLSA
diffe
1. c 2. N 3.0 4. CH3 5. H 6. C 7. c 8 CH3 9. N 10.0 11. H 12. CH3 total
Alanine DiwDtide., CT , 6.36 7.10 -4.98 -5.77 -23.49 -21.31 -0.93 0.20 -7.51 -6.86 1.43 1.76 9.46 14.06 -1.17 -1.10 -7.64 -16.08 -27.38 -30.57 -3.45 4.91 0.75 0.85 -58.55 -52.81
474 0.79 -2.18 -1.13 -0.65 -0.33 -4.6 -0.07 8.44 3.19 -8.36 -0.1 -5.74
1. c 2. N 3.0 4. CH3 5. H 6. C 7. c 8. CH3 9. N 10.0 11. H 12. CH3 total
Alanine Dipeptide, C5 13.53 15.39 -13.33 -15.57 -3 1.72 -32.94 -0.97 -2.45 2.22 3.57 1.43 2.33 5.21 7.38 -1.22 -1.55 -1 -34 4.26 -20.77 -22.27 -10.91 -8.93 -0.80 1.oo -58.67 -58.30
-1.86 2.24 1.22 1.48 -1.35 -0.9 -2.17 0.33 2.92 1.5 -1.98 -1.8 -0.37
1. CH3 2. N 3.0 4. CH3 5. H 6. C 7. c 8. CH3 9. N 10. 0 11. H 12. CH3 total
Alanine Dipeptide, a 7.61 8.37 -1.99 -2.90 -29.68 -29.97 -0.92 -1.51 -12.39 -10.28 1.45 2.52 11.50 15.64 -1.33 -0.87 -5.43 -12.03 -30.76 -32.85 -11.12 -4.77 0.76 2.08 -71.83 -67.02
-0.76 0.91 0.29 0.59 -2.11 -1.07 -4.14 0.46 6.6 2.09 -6.35 -1.32 -4.81
..
In kcal/mol. See Figure 1 for site-numbering scheme. Diff = ( Uw)0pm - ( U w ) 0 p ~ ~ .
TABLE V
Free Enernv of Solvation (in kcal/mol) OPLS solute
1. methanol 2. ethanol 3.2-propanol 4. acctone 5. methyl acetate 6. acetic acid 7. N,N-dimethylacetamide 8. acetamide 9. tram-N-methylacetamide 10. cis-N-methylacetamide 11. cis-N-methylacetamide dimer 12. ala-ala C7 13. ala-ala C5 14. ala-ala a 15. benzene 16. toluene 17. phenol 18. ammonium ion 19. acetate ion 20. ammonium acetate
&OF
-10.55 -9.27 -4.97 -0.32 -0.72 -8.94 -6.27 -12.90 -9.96 -7.33 -3.20 -20.43 -19.43 -27.31 1.38 2.96 -8.32 -72.15 -117.30 -51.35
OPLSA APmc
A p
-2.72 1.36 8.72 11.91 14.01 3.82 12.90 1.22 7.06 8.93 32.93 16.94 16.38 10.91 25.30 30.59 18.95 -62.23 -99.89 -28.45
-9.80 -9.66 -10.54 -2.32 -2.49 -7.82 -8.23 -13.91 -1 1.84 -9.23 -8.92 -22.10 -24.13 -29.96 -2.76 -2.72 -12.02 -72.22 -120.10 -53.88
AP-c 1.90 8.10 15.14 19.61 21.50 9.24 26.50 4.48 4.47 16.62 47.76 32.23 29.49 25.82 23.82 3 1.99 17.71 -62.36 -98.58 -26.94
expt -5.1014b -5.00 -4.8' -3.90 -3.30 -6.706 -8.54' -9.7' -10.07'
41.9,~0.8br -0.90 -6.6' -796 -800
a Cabani, S.;Gianni, P.; Mollica, V.;Lepori, L. J. Solution Chem. 1981,10,563. b Pearson, R. G. J. Am. Chem. Soc. 1986,108,6109. Ben-Naim, A.; Marcus, Y . J. Chem. Phys. 1984,81, 2016.
transition. We obtained 13.28 and 14.21 kcal/mol from OPLS and OPLSA,respectively.
Based on MC simulations, Jorgensen e t ale4' obtained values of -4.0 f 0.2,3.9f 0.2, and 5.0 f 0.2 kcal/mol for the respective
Lee and Maggiora
10184 The Journal of Physical Chemistry, Vol. 9 7 , No. 39, 1993
TABLE M: Chemical Group Contributions of Solute Molecules to the Free Energy of Solvation (in kcal/mol)* OPLS OPLSA 2-Propanol 1. CH3 2.CH 3. CH3 4.0-H
4.59 -1.88 4.82 -12.51 -4.97
total
2.17 -0.57 2.46 -14.60 -10.54
Acetone 6.35 -13.03 6.35 -0.32
1. CH3 2. c=o 3. CH3
total
1.60 -5.52 1.60 -2.32
0.04 0.0
c 2.0
Methyl Acetate 6.64 -14.32 6.96 -0.72
1. CH3 2. -0c-o 3. CH3
total
6.0
8.0
6.0
I
r(A)4'0
2.26 -8.46 3.72 -2.48
Benzene 1. CH
0.23
-0.46 -2.76
total 1.38 See Figure 1 for relevant structures.
I
3.04
2.0
I I
0.0 0.0
J 2.0 r(A)4'0
Figure 4. Solute-solventradial distribution functions for methanol (solid line), ethanol (dotted line), and 2-propanol (dashed line) obtained from the OPLSA potential energy function: (a, top) 0-H, and (b, bottom) 0-0,.
TABLE VII: Comparison of Free Energies of Solvation (in kcal/mol) XRISM"
0.0
0.0
2.0
6.0
8.0 L
4.0-
molecule methanol ethanol 2-propanol acetone acetamide acetic acid methylacetate benzene toluene phenol
0.0 0.0
2.0
4.0
6.0
c
Figure 3. Solute-solventradial distributionfunctions for methanol (solid
line), ethanol (dotted line), and 2-propanol (dashed line) obtained from the OPLS potential energy function: (a, top) 0-H, and (b, bottom) 0-0,.
solvation free-energy difference of acetic acid with acetamide, acetone, and methyl acetate. These values are in reasonable agreement with the correspondingexperimental values -3.0.2.8, and 3.4 kcal/mol. The current XRISM calculations yield values of -3.96, 8.62, and 7.62 kcal/mol for the OPLS parameter set and -6.09,5.5, and 5.33 kcal/mol for the OPLSA parameter set. Essex et al.6 calculated the solvation free-energy differences for methanol ethanol and ethanol propanol using the freeenergy perturbation method within an MD framework for various potential-energy function parameter sets. For the methanol
-
OPLSA
-10.55 -9.27 -4.97 -0.32 -12.90 -8.94 -0.72 1.38 2.96 -8.32 -117.30 -72.15
-9.80 -9.66 -10.54 -2.32 -13.91 -7.82 -2.49 -2.76 -2.72 -12.02 -120.1 -72.22
QMb -2.5
-6.1 -4.6 -0.7 -0.7
GB/SAC exptd -6.2 -5.2 -4.3 -3.2 -10.6 -6.5 -2.2 -1.0 -0.1 -6.3 -82.9 -90.8
-5.1 -5.0 -4.8 -3.9 -9.7 -6.7 -3.3 -0.9 -0.9 -6.6 -77.0 -79.0
acetate ion -79.5 ammonium ion -78.2 Values obtained from the GF approximation.b See ref 50. Generalized Born equation/solventaccessible surfacearea. See ref 49. d See Table V for references.
8.0
r(4
-
OPLS
-
ethanol "transformation", these workers obtained values of -0.52 and 0.52 kcal/mol with the OPLS parameter set. In the latter case, the 6-1 2 van der Waals terms on the hydrogen atoms of the TIP3P water model used to represent the solvent were removed. The XRISM calculations yield 1.28 and 0.14 kcal/mol for the OPLS and OPLSA parameter sets, respectively. Experimentally, the difference is about 0.1 kcal/mol. In addition to the simulation methods, continuum-based classical12and quantum-mechanical13methods show considerable promise. In contrast to all of the simulation-based studies, the performance of both hybrid methods was evaluated over a considerable number and variety of different molecular species. Table VI1 provides a comparison of the XRISM results obtained here with those obtained by the two hybrid procedures and with experiment. From the table it is clear that the classical hybrid
Solvation Thermodynamics by the XRISM Method method of Still and co-workers12is superior to both the XRISM or quantum-mechanical methods.13 Interestingly, the quantummechanical method yields the best agreement with experiment for the acetate and ammonium ions. This is not totally surprising considering that the partial charges on the solute atoms are obtained in a self-consistent fashion in the quantum-mechanical method while they are fixed in both the XRISM and classicalhybrid methods. However, it should be noted that these hybrid methods involve the empirical parametrization of solventaccessible surface tensions and a solute-solvent electrostatic polarization term, unlike the XRISM method where no such parametrization was used.
4. Summary The aqueous solvation free energies and the average solutesolvent interaction energies of a series of organic molecules were calculated by the XRISM method within the HNC and GF approximations for both the OPLS and OPLSA parameter sets. The results show that solvation free energies obtained from the GF approximation are generally in reasonable agreement with experiment and in better agreement than those obtained from the HNC approximation, the only exceptions being methanol and acetate ion. No obvious trend was observed regarding the ability of either potential-energy function to provide reliable estimates of experimental solvation free energies. This lack of any trend was also observed in other commonly used potential-energyfunctionssuch as AMBER25 and CHARMM,4* although the data are less extensive49 than those for OPLS and OPLSA. The significant role played by electrostatic factors was clearly evident in the results reported here and points to the need for better descriptions of electrostatics in existing potential-energy functions. This point is evidenced by the fact that improved results were obtained for molecules using point charges fit to quantum-mechanicallycalculatedelectrostatic potential^;^ a fuller account of this work will be given in a subsequent publication. Moreover,the need to further investigate the underlying basis by which different approximations (e.g., the GF approximation) influence the results of XRISM calculations is clearly necessary if the XRISM method is to attain its full potential as a practical procedure for calculating the thermodynamic properties of solutions. Acknowledgment. We thank Professor B. M. Pettitt for many helpful discussions and sending us the program for logarithmic fast Fourier transformations. We also thank Dr. H. A. Yu for many interesting discussions and Professor B. Honig for sending coordinates of the solute molecules and also for many stimulating discussions. References and Notes (1) Brooks,C. L., 111;Karplus, M.; Pettitt, B. M. Proteins: A Theoretical Perspective of Dynamics, Structure, and Thermodynamics; John Wiley & Sons: New York, 1988. (2) Jorgensen, W. L.; Ravimohan, C. J . Chem. Phys. 1985, 83, 3050. (3) Jorgensen, W. L.; Gao, J. J. Am. Chem. SOC.1988, 110,4212. (4) Bash, P. A.; Singh, U. C.; Langridge, R.; Kollman, P. A. Science 1987, 236, 564.
The Journal of Physical Chemistry, Vol. 97, No. 39, 1993
10185
( 5 ) Cieplak, P.; Kollman, P. J. Comput. Chem. 1991, 12, 1232. (6) Essex, J. W.; Reynolds, C. A,; Richards, W. G.J . Am. Chem. Soc. 1992, 114, 3634. (7) Tobias, D. J.; Brooks, C. L.,111. J. Phys. Chem. 1992, 96, 3864. (8) Beveridge, D. L., Jorgensen, W. L., Eds. Computer Simulation of Chemical and Biomolecular Systems; Ann. N.Y. Acad. Sci. 1986. (9) Brooks, C. L., 111; Pettitt, B. M.; Karplus, M. J. Chem. Phys. 1985, 83, 5897. (10) Brooks, C. L., 111. J . Chem. Phys. 1987,86, 5156. (1 1) Jean-Charles, A.; Nicholls, A.; Sharp, K.; Honig, B.; Tempczyk, A.; Hendrickson, T. F.; Still, W. C. J . Am. Chem. Soc. 1991, 113, 1454. (12) Still, W. C.; Tempczyk, A.; Hawley, R. C.; Hendrickson, T. J. Am. Chem. SOC.1990, 112, 6127. (13) Cramer, C . J.; Truhlar, D. C. J . Am. Chem. SOC.1991, 113,8305. (14) Chandler, D.; Andersen, H. C. J . Chem. Phys. 1972.57, 1930. (15) Hirata, F.; Pettitt, B. M.; Rossky, P. J. J. Chem. Phys. 1982,77,509. (16) Pettitt, B. M.; Rossky, P.J. J. Chem. Phys. 1982, 77, 1451. (17) Pettitt, B. M.; Rossky, P. J. Israel J. Chem. 1986, 27, 156. (18) Pettitt, B. M.; Karplus, M. Chem. Phys. Lett. 1985, 121, 194. (19) Pettitt, B. M.; Karplus, M.; Rossky, P. J. J . Phys. Chem. 1986,90, 6335. (20) Chandler, D.; Singh, Y.; Richardson, D. J . Chem. Phys. 1984, 81, 1975. (21) Chiles, R. A.; Rossky, P. J. J. Am. Chem. Soc. 1984, 106, 6867. (22) Yu, H.; Karplus, M. J. Am. Chem. SOC.1990, 112, 5706. (23) Yu, H.; Pettitt, B. M.; Karplus, M. J . Am. Chem. SOC.1991, 113, 2425. (24) Wolfenden, R. Biochemistry 1978, 17, 201. (25) Spencer, J. N.;Berger,S. C.;Powell,C. R.; Henning, B. D.; Furman, G.S.; Loffredo,W. M.; Rydberg,E. M.;Neubert, R. A.;Shoop,C. E.; Blauch, D. N. J . Phys. Chem. 1981,85, 1236. (26) Jorgensen, W. J.; Tirado-Rives, J. J . Am. Chem. Soc. 1988, 110, 1657. (27) Weiner, S.J.; Kollman, P. A.; Case, D. A.; Singh, U. C.; Ghio, C.; Alagona, G.;Profeta, S.;Weiner, P. J. Am. Chem. Soc. 1984, 106, 765. (28) McQuarrie, D. A. Statistical Mechanics; Harper and Row: New York, 1973. (29) Rossky, P. J.; Dale, W. D. T. J . Chem. Phys. 1990. 73, 2457. (30) Singer, S.J.; Chandler, D. Mol. Phys. 1985, 55, 621. (31) Jorgensen, W. L. J . Am. Chem. SOC.1981, 103, 335. (32) Yu,H. A,; Karplus, M. J . Chem. Phys. 1988,89, 2366. (33) Narten, A. H. J. Chem. Phys. 1972, 56, 5681. (34) Egelstaff, A. Adv. Chem. Phys. 1983, 58, 1. (35) Narten, A. H.; Thiessen, W. E.; Blum, L. Science 1982, 217, 1033. (36) Ben-Naim, A.; Marcus, Y. J. Chem. Phys. 1984, 81, 2016. (37) Wolfenden, R.; Andersson, L.; Cullis, P. M.; Southgate, C. C. B. Biochemistry 1981, 20, 849. (38) Radzicka, A,; Wolfenden, R. Biochemistry 1988, 27, 1664. (39) Cabani, S.;Gianni, P.; Mollica, V.; Lepori, L. J. Solution Chem. 1981, 10, 563. (40) Ichiye, T.; Chandler, D. J . Phys. Chem. 1988, 92, 5257. (41) Yu, H. A.; Roux, B.; Karplus, M. J . Chem. Phys. 1990, 92, 5020. (42) Rossky, P. J.; Pettitt, B. M.; Stell, G.Mol. Phys. 1983, 50, 1263. (43) Rossky, P. J.; Karplus, M. J . Am. Chem. SOC.1979, 101, 1913. (44) Mezei, M.; Mehrotra, P. K.; Beveridge, D. L. J . Am. Chem. Soc. 1985, 107, 2239. (45) Brady, J.; Karplus, M. J . Am. Chem. SOC.1985, 107, 6103. (46) Anderson, A,; Hermans, J. Proteins 1988, 3, 262. (47) Jorgensen, W. L.; Briggs, J. M.; Contreras, M. L. J . Phys. Chem. 1990, 94, 1683. (48) Brooks, B. R.;Plafson, B. D.; States, D. J.; Swaminathan,S.; Karplus, M. J. Comput. Chem. 1983,4, 187. (49) Lee,P. H. Unpublished results. (50) Mohamadi, F.; Richards,N. G.J.;Guida, W. C.; Liskamp, R.; Lipton, M.; Caufield, C.; Chang, G.;Hendrickson, T.; Still, W. C. J . Comput. Chem. 1990, 11, 440. (5 1) van Gunstren, W. F. Computer Simulation of Biomolecular Systems, van Gunstren, W. F., Weiner, P. K., Eds.; ESCOM: Leiden, 1989; Vol. 27. (52) Talman, J. D. J . Comput. Phys. 1978, 29, 35. (53) Rossky, P. J.; Friedman, H. L. J. Chem. Phys. 1980. 72, 5694. (54) Davis, P. J., Rabinowitz, P., Eds. Methods ofNumericallntegration; Academic Press, 1984. ( 5 5 ) Lim, C.; Bashford, D.; Karplus, M. J . Phys. Chem. 1991,95,5610.
