- Email: [email protected]
Benzene solubility in water: A reassessment
Chemical Physics Letters 429 (2006) 114–118 www.elsevier.com/locate/cplett
Benzene solubility in water: A reassessment Giuseppe Graziano
*
Dipartimento di Scienze Biologiche ed Ambientali, Universita` del Sannio, Via Port’Arsa 11 – 82100 Benevento, Italy Received 12 July 2006; in final form 24 July 2006 Available online 5 August 2006
Abstract It is shown that the results of molecular dynamics simulations on the hydration thermodynamics of benzene at room temperature [Schravendijk and van der Vegt, J. Chem. Theory Comput. 1 (2005) 643] are in line with a former theoretical analysis [Graziano and Lee, J. Phys. Chem. B 105 (2001) 10367]. In fact: (a) the benzene–water van der Waals interaction energy proves to be larger in magnitude than the work of cavity creation and is able to account for the experimental finding that the hydration of benzene is a spontaneous process under the Ben-Naim standard conditions around room temperature; (b) the weak benzene–water H-bonds do not provide a significant contribution to benzene solubility in water because the favorable enthalpic component is almost entirely compensated for by an unfavorable entropic component. This enthalpy–entropy compensation occurs because the H-bonding potential of benzene is not strong. 2006 Elsevier B.V. All rights reserved.
1. Introduction It is well established that the solubility of benzene in water is markedly larger than that of an alkane of similar size [1]. At room temperature, the Ben-Naim [2] standard Gibbs energy change DG associated with the hydration (i.e., gas-to-water transfer) of benzene is a negative quantity, whereas it is a large positive quantity for n-alkanes [3]. Clearly, it is important to provide an explanation at a molecular level of the relatively large solubility of benzene in water. Makhatadze and Privalov [3] suggested that the weak H-bonds that benzene forms with water provide a favorable energetic contribution that determines the solubility enhancement [4,5]. Graziano and Lee [6], G&L, reached a different conclusion on the basis of a precise division of the hydration thermodynamic quantities grounded on statistical mechanics. They concluded that: (a) the relatively large solubility of benzene in water is caused by stronger solute–solvent van der Waals attractive interactions with respect to those of alkanes; (b) the weak benzene–water H-bonds do not afford a significant contribution to DG because the favorable enthalpic term is to a *
Fax: +39 0824 23013. E-mail address: [email protected].
0009-2614/$ - see front matter 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2006.08.006
large extent compensated for by an unfavorable entropic term. Graziano [7] strengthened this conclusion by pointing out that the polarizability of benzene is markedly larger than that of an alkane of similar size. In a recent article Schravendijk and van der Vegt [8], S&vdV, performed molecular dynamics, MD, simulations of two benzene models in simple point charge, SPC, water [9] to try to settle down the matter. In order to address the role played by the weak benzene–water H-bonds, S&vdV studied a real benzene model that favors the formation of such weak H-bonds as well as a van der Waals, vdW, benzene model obtained by removing all partial charges of the first model and keeping the dispersion interactions unaffected. The main point of the S&vdV work was the calculation of DG by means of a thermodynamic integration procedure for the two models of benzene, using a constant pressure–temperature ensemble within the GROMACS simulation package [10]. They obtained, at 302 K, DG = 4.8 kJ mol1 for the real benzene model and +1.0 kJ mol1 for the vdW benzene model. On the basis of these two qualitatively different DG values, S&vdV [8] concluded that the partial charges on the benzene ring leading to the existence of the weak H-bonds with two water molecules located over the two faces of the ring are necessary to account for the relatively large solubility of benzene in
G. Graziano / Chemical Physics Letters 429 (2006) 114–118
water. Note that the experimental value at 25 C is DG = 3.6 kJ mol1 for benzene [3]. In the present Letter I would like to show that the numbers coming from a careful analysis of the S&vdV work are not in contrast with the conclusion reached by G&L [6]. Specifically, the MD simulations of S&vdV confirm that: (a) the benzene–water van der Waals interaction energy overwhelms in magnitude the work of cavity creation; (b) the enthalpic contribution of the weak benzene–water Hbonds is to a large extent compensated for by a corresponding entropic contribution. Since some confusion seems to be caused by an incomplete understanding of the approach devised by G&L, I start by spelling out in detail the fundamental points of their statistical mechanical theory. 2. Theory of solvation Solvation corresponds to the transfer of a solute molecule from a fixed position in the ideal gas phase to a fixed position in a solvent at constant temperature and pressure [2]. The process can be treated as the insertion of an external perturbing potential W(X), where X is a vector representing a single configuration of the N molecules of the solvent [11,12]. By using the Widom’s potential distribution theorem [13], DG is given by: DG ¼ RT lnhexp½WðXÞ=RT ip
ð1Þ
where the subscript p means that the ensemble average is performed over the pure solvent configurations, assuming an NPT ensemble (i.e., keeping fixed the number of molecules, the pressure and the temperature), whose probability density function is: Z qp ðXÞ ¼ exp½H ðXÞ=RT = exp½H ðXÞ=RT dX ð2Þ where H(X) = U(X) + P Æ V(X) is the enthalpy function of a configuration, U(X) and V(X) are the corresponding internal energy and volume, and the denominator is the isobaric partition function of the system. A liquid is a condensed state of the matter and the insertion of a solute molecule requires the exclusion of the solvent molecules from a region of space that is suitable to host the solute. Therefore, the perturbing potential can be expressed as [11]: exp½WðXÞ=RT ¼ fðXÞ exp½wðXÞ=RT
ð3Þ
where f(X) can assume only the values of one or zero depending on whether or not there is a suitable cavity to host the solute molecule in the given solvent configuration; and w(X) represents the attractive potential between the solute molecule and the surrounding solvent molecules. Insertion of Eq. (3) into Eq. (1) leads to: DG ¼ RT lnhfðXÞip RT lnhexp½wðXÞ=RT ic
ð4Þ
where the subscript c means that the ensemble average is performed over the solvent configurations possessing a
115
suitable cavity to host the solute whose probability density function is: Z qc ðXÞ ¼ fðXÞ exp½H ðXÞ=RT = fðXÞ exp½H ðXÞ=RT dX
ð5Þ
Such configurations are only a small fraction of the total solvent configurations (i.e., those for which the function f(X) is equal to one). Eq. (4) indicates that DG is the sum of two terms: the work to create the cavity and the work to turn on the attractive solute–solvent potential. However, this does not mean additivity of contributions: the attractive solute–solvent potential is switched on given that the cavity has already been created, and the average is calculated by using the conditional probability density function qc(X). The work to create the cavity is: DGc ¼ RT lnhfðXÞip RT ln pinsertion
ð6Þ
It is the reversible work to single out the configurations containing the cavity from the ensemble of pure solvent configurations [14,15] and, according to a theorem of statistical mechanics [16], it is exactly related to the probability that the desired region of space in the liquid is devoid of liquid molecules (i.e., the so-called insertion probability). The work to turn on the attractive solute–solvent potential is: DGa ¼ RT lnhexp½wðXÞ=RT ic
ð7Þ
By setting j = w(X) Æw(X)æc, expanding in power series the exponential function, and keeping in mind that Æjæc ” 0, one obtains [6,7]: DGa ffi hwðXÞic ½hj2 ic =2RT
ð8Þ
if higher order terms are ignored. It is simple to show that the fluctuation term Æj2æc/2RT is of entropic origin [11,12]. When the attractive potential w(X) is weak, the fluctuations in the value of Æw(X)æc are small, and the second term on the right-hand side of Eq. (8) can be neglected. Under this condition, the solvation Gibbs energy change is given by: DG ¼ DGc þ hwðXÞic
ð9Þ
For the hydration of alkanes and noble gases, which are unable to form H-bonds with water, the fluctuation term Æj2æc/2RT is small because w(X) represents the solute– water van der Waals interactions which are weak compared to the water–water H-bonds [17]. In such cases it should be safe to assume that Æw(X)æc is equal to the direct solute– water interaction energy DUsw (i.e., the subscript sw stands for solute–water interaction). When the solute molecule can form H-bonds with water, the condition of small fluctuations for the solute–water interaction energy could not hold [17]. However, in the ensemble of the water configurations possessing a suitable cavity to host the solute, the quantity Æw(X)æc should not account for the solute–water H-bond energy. A negligible
116
G. Graziano / Chemical Physics Letters 429 (2006) 114–118
fraction of water molecules in the qc(X) distribution would have the correct orientation to form a H-bond with the solute molecule inserted in the cavity because water molecules have not yet reorganized in response to the attractive potential of the solute, since the latter acts as a ghost [6,7]. Thus Æw(X)æc should account only for the van der Waals component of the solute–water interaction energy, Æw(X)æc = DUsw(vdW), the fluctuation term is expected to be small, and the contribution of solute–water H-bonds should be included in the solvent reorganization term. By adopting this line, G&L found that [6,7]: (a) the benzene–water van der Waals interaction energy is larger in magnitude than the work of cavity creation and is able to account for the experimental negative DG value; (b) the fluctuation term proves to be small for the hydration of benzene in view of the weakness and negligible anisotropy of the van der Waals component of the benzene–water potential; (c) the reorientation of water molecules to form H-bonds with the aromatic ring is characterized by enthalpy–entropy compensation because the H-bonding potential of benzene is not strong. Further light can be shed on the matter by using an inverse relationship proved by Widom [13], so that the work to turn on the solute–solvent attractive potential is: DGa ¼ RT lnhexp½wðXÞ=RT ic ¼ RT lnhexp½wðXÞ=RT ia
ð10Þ
where the subscript a means that the ensemble average is performed over the solution configurations possessing the cavity occupied by the solute molecule that interacts with the surrounding solvent molecules. The probability density function of such configurations is: Z qa ðXÞ ¼ fðXÞ expf½wðXÞ þ H ðXÞ=RT g= fðXÞ expf½wðXÞ þ H ðXÞ=RT gdX
ð11Þ
In other words, in such ensemble the solute molecule does not act as a ghost but fully interacts with and influences the surrounding solvent molecules. In general, it should be noted that the Widom’s inverse relationship is valid only when the external potential W(X) is not infinite, and that W(X) to insert a cavity at a fixed position is infinite for many solvent configurations if the position is occupied by solvent molecules. This is why cavity creation is a special process and is separated out first, before applying the Widom’s inverse relationship. By setting c = w(X)Æw(X)æa, expanding in power series the exponential function in the last expression of Eq. (10), and keeping in mind that Æcæa” 0, one obtains: 2
DGa ffi hwðXÞia þ ½hc ia =2RT
ð12Þ
if higher order terms are ignored. Even though Eq. (12) looks like Eq. (8), its physical meaning is entirely different: (a) Æw(X)æa accounts for the full solute–water interaction energy, including possible H-bonds, Æw(X)æa = DUsw(full); (b) the fluctuation term Æc2æa/2RT, always of entropic ori-
gin, should be a large quantity for a solute able to form H-bonds with water since the H-bonding potential is strong and is characterized by a marked anisotropy due to the geometric requirements for the occurrence of H-bonds. Actually, for strongly H-bonding solutes such as n-alcohols, the truncation of the power series expansion cannot be a correct procedure [18]. If the fluctuation term Æc2æa/2RT is not a negligible quantity, the solvation Gibbs energy change is given by: DG ¼ DGc þ hwðXÞia þ ½hc2 ia =2RT
ð13Þ
Furthermore, since both Eqs. (8) and (12) are correct from the statistical mechanical point of view, by equating the two expressions, one obtains [6]: hj2 ic =2RT ¼ ½hc2 ia =2RT þ ðhwðXÞia hwðXÞic Þ
ð14Þ
that should allow a reliable evaluation of the fluctuation term Æj2æc/2RT. Now, the numerical values calculated by S&vdV can be inserted in the above equations in order to test the validity of the G&L analysis of benzene solubility in water. 3. Analysis of the S&vdV work Beyond to calculate DG, S&vdV [8] computed directly from the MD simulation runs the benzene–water interaction energy DUsw. It resulted that, at 302 K, DUsw = 59.4 kJ mol1 for the real benzene model and 45.6 kJ mol1 for the vdW benzene model, respectively. S&vdV assumed also that the water reorganization process is characterized by a perfect enthalpy–entropy compensation [19] (i.e., the subscript ww stands for water–water interaction): DH ww ¼ T DS ww
ð15Þ
Therefore DG was solely given by the enthalpic and entropic contributions due to solute–water interactions (note, in general, that in a condensed phase of the matter, such as a liquid, the difference between energy and enthalpy is very small): DG ¼ DU sw T DS sw
ð16Þ
S&vdV, by using Eq. (16), obtained the solute–water entropy contribution at 302 K: T Æ DSsw = DUsw DG = 54.6 kJ mol1 for the real benzene model and 46.6 kJ mol1 for the vdW benzene model, respectively. In order to provide an interpretation for these large and negative values of the benzene–water entropy change, S&vdV used a theoretical relationship derived by Sanchez and colleagues [20]: T DS sw ¼ RT ln pinsertion ½hc2 ia =2RT
ð17Þ
where the first term on the right-hand-side is minus the work of cavity creation, DGc, and the second term on the right-hand-side is a fluctuation term identical to that of Eq. (12) – it accounts for the fluctuation of the attractive solute–solvent potential energy w averaged over the ensemble of solution configurations in which there are no solute–
G. Graziano / Chemical Physics Letters 429 (2006) 114–118
solvent repulsions and the attractive solute–solvent potential is turned on. S&vdV performed MD trajectories 90 ns long for the two models of benzene, and found that the fluctuation term amounts to 12.8 kJ mol1 for the real benzene model and 3.9 kJ mol1 for the vdW benzene model [8]. This means that the work to create in SPC water at 302 K a cavity suitable to host benzene is: (a) DGc = RT Æ ln pinsertion = 54.6 12.8 = 41.8 kJ mol1 for the real benzene model; (b) DGc = RT Æ ln pinsertion = 46.6 3.9 = 42.7 kJ mol1 for the vdW benzene model. It is important to note that these two values should be identical because DGc is a property of the pure liquid and cannot depend on the potential energy of the solute molecule [14,15]; in the following the corresponding mean value, DGc = 42.3 kJ mol1 will be used. The two numbers derived from the simulations of S&vdV are close to each other and are reliable. Indeed (a) from the DGc values calculated by van Gunsteren and colleagues [21] in SPC water at room temperature, DGc 43 kJ mol1 for a cavity suitable to host a benzene molecule; (b) according to scaled particle theory [22], DGc = 41.3 kJ mol1, using the experimental density of water at 25 C, and the hard sphere diameters r(ben˚ and r(water) = 2.80 A ˚ [6,23]. zene) = 5.26 A Clearly, it is now possible to apply Eqs. (9) and (13) to the two benzene models of S&vdV. By using Eq. (9) for the vdW benzene model and considering that Æw(X)æc = DUsw(vdW), it results: DG ffi DGc þ DU sw ðvdWÞ ¼ 42:3 45:6 ¼ 3:3 kJ mol1 ð18Þ By using Eq. (13) for the real benzene model and considering that Æw(X)æa = DUsw(full), it results: DG ¼ DGc þ DU sw ðfullÞ þ ½hc2 ia =2RT ¼ 42:3 59:4 þ 12:8 ¼ 4:3 kJ mol1
ð19Þ
These calculations indicate unequivocally that both Eqs. (9) and (13) are right and provide a similar numerical result and the same physical explanation for the relatively large solubility of benzene in water. Eq. (18) indicates that the experimental negative DG value of benzene in water is accounted for by simply considering the balance between the work of cavity creation and the contribution of van der Waals interactions. But this is the same explanation emerging from Eq. (19), because the favorable enthalpic contribution of the weak benzene–water H-bonds is compensated for by a corresponding unfavorable entropic contribution: DU sw ðfullÞ þ ½hc2 ia =2RT ¼ 59:4 þ 12:8 ¼ 46:6 kJ mol1 DU sw ðvdWÞ ð20Þ The entropic fluctuation term practically cancels the difference in the benzene–water interaction energy calculated for the two different benzene models by S&vdV. This is exactly the manifestation of the enthalpy–entropy compensation characterizing the weak benzene–water H-bonds,
117
as indicated by G&L [6]. The magnitude of the fluctuation term Æc2æa/2RT expresses the extent to which the attractive solute–solvent interactions bias positions and orientations of the vicinal solvent molecules, causing a reduction of the accessible configuration space that, in turn, leads to an entropy decrease. It is also possible to obtain an estimate of the fluctuation term Æj2æc/2RT in Eq. (8) using Eq. (14) and the values calculated by S&vdV: hj2 ic =2RT ffi ½hc2 ia =2RT þ DU sw ðfullÞ DU sw ðvdWÞ ¼ 12:8 59:4 þ 45:6 ¼ 1:0 kJ mol1
ð21Þ
This estimate indicates that such fluctuation quantity should be small [6]. It is also important to note that Æj2æc/2RT is a negative quantity implying that the reorganization of water molecules to form the weak H-bonds with the benzene ring is a spontaneous process, even though characterized by an almost complete enthalpy–entropy compensation. S&vdV determined the fluctuation term also for the vdW benzene model and the value was 3.9 kJ mol1, that cannot be considered negligible. However, according to the statistical-mechanical approach of G&L, the quantity Æj2æc/2RT has to be calculated over the ensemble of solvent configurations possessing the suitable cavity, which is different from the ensemble of solution configurations used by S&vdV to calculate the fluctuation term for the vdW benzene model. 4. Discussion The performed analysis indicates that the values obtained from the MD simulations of S&vdV are not in contrast with the G&L approach: the benzene–water van der Waals interaction energy is larger in magnitude than the work of cavity creation and is able to account for the experimental negative DG value. It is worth noting that the van der Waals benzene–water interactions are stronger than those of an aliphatic hydrocarbon of similar molecular size with water, due to the larger polarizability of the aromatic ring [7]. Consider, for instance, that, even though ˚ 3, is similar to that of the polarizability of c-hexane, 10.78 A 3 ˚ benzene, 10.32 A [22], the c-hexane molecule is significantly larger than that of benzene. The effective hard ˚ for c-hexane and 5.26 A ˚ for sphere diameters are 5.63 A benzene [12], so that the molecular volume of the latter is 82% of that of the former. Such a value should correspond also to the polarizability ratio because this physical quantity is proportional to molecular volume. However, experimental data indicate that the polarizability of benzene is 96% of that of c-hexane [22]. In addition, van Gunsteren and colleagues [24] obtained, at room temperature, DUsw = 35.0 kJ mol1 for i-butane and 36.4 kJ mol1 for n-butane; these values are about 10 kJ mol1 smaller in magnitude than the benzene–water van der Waals interaction energy, DUsw(vdW) = 45.6 kJ mol1. This finding
118
G. Graziano / Chemical Physics Letters 429 (2006) 114–118
is in line with the analysis by G&L [6,7], but was not underscored by S&vdV. A further point should be considered: if the contribution of the weak benzene–water H-bonds is small due to enthalpy–entropy compensation, the thermodynamic integration procedure should produce the same DG estimate for the two models of benzene devised by S&vdV. Indeed, thermodynamic integration takes into account solely the effective contributions to the DG quantity. The difference found by S&vdV, DG = 4.8 kJ mol1 for the real benzene model and+1.0 kJ mol1 for the vdW benzene model, is not small, but cannot be considered to be large. In order to gain perspective, one has to remind that, for an alkane having a size comparable to that of benzene (i.e., i-butane and nbutane), DG amounts to 9–10 kJ mol1 at 25 C [2,3]. The latter numbers indicate that the vdW benzene model of S&vdV is markedly less hydrophobic than an alkane of similar size. This result was neglected by S&vdV, even though it emerged directly from their MD simulations. In addition, S&vdV showed that a change of GROMOS parameters produces a marked effect on the DG value calculated for the real benzene model: DG = 4.8 kJ mol1 for the first-set, and 6.7 kJ mol1 for the second-set [8]. S&vdV did not show the DG value obtained for the vdW benzene model using the second-set of GROMOS parameters. One could speculate that, by using the second-set of GROMOS parameters, DG would be a negative quantity also for the vdW benzene model. In conclusion, the performed analysis points out that the numerical values obtained from MD simulations by S&vdV indicate that: (a) the benzene–water van der Waals interaction energy overwhelms the work of cavity creation, accounting for the experimental negative DG value; (b) the weak benzene–water H-bonds do not really influence the benzene solubility in water because the favorable enthalpic component is almost entirely compensated for by an unfavorable entropic component.
Acknowledgement I thank Dr. B. Lee (Center for Cancer Research, NCI, NIH, Bethesda, MD) for reading earlier drafts of the manuscript and providing useful suggestions. References [1] C. McAuliffe, J. Phys. Chem. 70 (1966) 1267. [2] A. Ben-Naim, Solvation Thermodynamics, Plenum, New York, 1987. [3] G.I. Makhatadze, P.L. Privalov, Biophys. Chem. 50 (1994) 285. [4] S. Suzuki, P.G. Green, R.E. Bumgarner, S. Dasgupta, W.A. Goddard, G.A. Blake, Science 257 (1992) 942. [5] P. Linse, J. Am. Chem. Soc. 112 (1990) 1744. [6] G. Graziano, B. Lee, J. Phys. Chem. B 105 (2001) 10367. [7] G. Graziano, Biophys. Chem. 110 (2004) 249. [8] P. Schravendijk, N.F.A. van der Vegt, J. Chem. Theory Comput. 1 (2005) 643. [9] H.J.C. Berendsen, W.F. van Gunsteren, J.P.M. Postma, J. Hermans, in: B. Pullman (Ed.), Intermolecular Forces, Reidel, Dordrecht, 1981, p. 331. [10] H.J.C. Berendsen, D. van der Spoel, R. van Drunen, Comput. Phys. Commun. 91 (1995) 43. [11] B. Lee, Biopolymers 31 (1991) 993. [12] G. Graziano, Can. J. Chem. 80 (2002) 401. [13] B. Widom, J. Phys. Chem. 86 (1982) 869. [14] B. Lee, J. Chem. Phys. 83 (1985) 2421. [15] G. Graziano, J. Phys. Soc. Jpn. 69 (2000) 1566. [16] R.C. Tolman, The Principles of Statistical Mechanics, Oxford University Press, London, 1938. [17] B. Lee, Biophys. Chem. 51 (1994) 271. [18] G. Graziano, J. Phys. Chem. B 109 (2005) 12160. [19] H.A. Yu, M. Karplus, J. Chem. Phys. 89 (1988) 2366. [20] I.C. Sanchez, T.M. Truskett, P.J. In’t Veld, J. Phys. Chem. B 103 (1999) 5106. [21] T.C. Beutler, D.R. Beguelin, W.F. van Gunsteren, J. Chem. Phys. 102 (1995) 3787. [22] R.A. Pierotti, Chem. Rev. 76 (1976) 717. [23] G. Graziano, Chem. Phys. Lett. 396 (2004) 226. [24] D. Trzesniak, N.F.A. van der Vegt, W.F. van Gunsteren, Phys. Chem. Chem. Phys. 6 (2004) 697.
Benzene solubility in water: A reassessment Giuseppe Graziano
*
Dipartimento di Scienze Biologiche ed Ambientali, Universita` del Sannio, Via Port’Arsa 11 – 82100 Benevento, Italy Received 12 July 2006; in final form 24 July 2006 Available online 5 August 2006
Abstract It is shown that the results of molecular dynamics simulations on the hydration thermodynamics of benzene at room temperature [Schravendijk and van der Vegt, J. Chem. Theory Comput. 1 (2005) 643] are in line with a former theoretical analysis [Graziano and Lee, J. Phys. Chem. B 105 (2001) 10367]. In fact: (a) the benzene–water van der Waals interaction energy proves to be larger in magnitude than the work of cavity creation and is able to account for the experimental finding that the hydration of benzene is a spontaneous process under the Ben-Naim standard conditions around room temperature; (b) the weak benzene–water H-bonds do not provide a significant contribution to benzene solubility in water because the favorable enthalpic component is almost entirely compensated for by an unfavorable entropic component. This enthalpy–entropy compensation occurs because the H-bonding potential of benzene is not strong. 2006 Elsevier B.V. All rights reserved.
1. Introduction It is well established that the solubility of benzene in water is markedly larger than that of an alkane of similar size [1]. At room temperature, the Ben-Naim [2] standard Gibbs energy change DG associated with the hydration (i.e., gas-to-water transfer) of benzene is a negative quantity, whereas it is a large positive quantity for n-alkanes [3]. Clearly, it is important to provide an explanation at a molecular level of the relatively large solubility of benzene in water. Makhatadze and Privalov [3] suggested that the weak H-bonds that benzene forms with water provide a favorable energetic contribution that determines the solubility enhancement [4,5]. Graziano and Lee [6], G&L, reached a different conclusion on the basis of a precise division of the hydration thermodynamic quantities grounded on statistical mechanics. They concluded that: (a) the relatively large solubility of benzene in water is caused by stronger solute–solvent van der Waals attractive interactions with respect to those of alkanes; (b) the weak benzene–water H-bonds do not afford a significant contribution to DG because the favorable enthalpic term is to a *
Fax: +39 0824 23013. E-mail address: [email protected].
0009-2614/$ - see front matter 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2006.08.006
large extent compensated for by an unfavorable entropic term. Graziano [7] strengthened this conclusion by pointing out that the polarizability of benzene is markedly larger than that of an alkane of similar size. In a recent article Schravendijk and van der Vegt [8], S&vdV, performed molecular dynamics, MD, simulations of two benzene models in simple point charge, SPC, water [9] to try to settle down the matter. In order to address the role played by the weak benzene–water H-bonds, S&vdV studied a real benzene model that favors the formation of such weak H-bonds as well as a van der Waals, vdW, benzene model obtained by removing all partial charges of the first model and keeping the dispersion interactions unaffected. The main point of the S&vdV work was the calculation of DG by means of a thermodynamic integration procedure for the two models of benzene, using a constant pressure–temperature ensemble within the GROMACS simulation package [10]. They obtained, at 302 K, DG = 4.8 kJ mol1 for the real benzene model and +1.0 kJ mol1 for the vdW benzene model. On the basis of these two qualitatively different DG values, S&vdV [8] concluded that the partial charges on the benzene ring leading to the existence of the weak H-bonds with two water molecules located over the two faces of the ring are necessary to account for the relatively large solubility of benzene in
G. Graziano / Chemical Physics Letters 429 (2006) 114–118
water. Note that the experimental value at 25 C is DG = 3.6 kJ mol1 for benzene [3]. In the present Letter I would like to show that the numbers coming from a careful analysis of the S&vdV work are not in contrast with the conclusion reached by G&L [6]. Specifically, the MD simulations of S&vdV confirm that: (a) the benzene–water van der Waals interaction energy overwhelms in magnitude the work of cavity creation; (b) the enthalpic contribution of the weak benzene–water Hbonds is to a large extent compensated for by a corresponding entropic contribution. Since some confusion seems to be caused by an incomplete understanding of the approach devised by G&L, I start by spelling out in detail the fundamental points of their statistical mechanical theory. 2. Theory of solvation Solvation corresponds to the transfer of a solute molecule from a fixed position in the ideal gas phase to a fixed position in a solvent at constant temperature and pressure [2]. The process can be treated as the insertion of an external perturbing potential W(X), where X is a vector representing a single configuration of the N molecules of the solvent [11,12]. By using the Widom’s potential distribution theorem [13], DG is given by: DG ¼ RT lnhexp½WðXÞ=RT ip
ð1Þ
where the subscript p means that the ensemble average is performed over the pure solvent configurations, assuming an NPT ensemble (i.e., keeping fixed the number of molecules, the pressure and the temperature), whose probability density function is: Z qp ðXÞ ¼ exp½H ðXÞ=RT = exp½H ðXÞ=RT dX ð2Þ where H(X) = U(X) + P Æ V(X) is the enthalpy function of a configuration, U(X) and V(X) are the corresponding internal energy and volume, and the denominator is the isobaric partition function of the system. A liquid is a condensed state of the matter and the insertion of a solute molecule requires the exclusion of the solvent molecules from a region of space that is suitable to host the solute. Therefore, the perturbing potential can be expressed as [11]: exp½WðXÞ=RT ¼ fðXÞ exp½wðXÞ=RT
ð3Þ
where f(X) can assume only the values of one or zero depending on whether or not there is a suitable cavity to host the solute molecule in the given solvent configuration; and w(X) represents the attractive potential between the solute molecule and the surrounding solvent molecules. Insertion of Eq. (3) into Eq. (1) leads to: DG ¼ RT lnhfðXÞip RT lnhexp½wðXÞ=RT ic
ð4Þ
where the subscript c means that the ensemble average is performed over the solvent configurations possessing a
115
suitable cavity to host the solute whose probability density function is: Z qc ðXÞ ¼ fðXÞ exp½H ðXÞ=RT = fðXÞ exp½H ðXÞ=RT dX
ð5Þ
Such configurations are only a small fraction of the total solvent configurations (i.e., those for which the function f(X) is equal to one). Eq. (4) indicates that DG is the sum of two terms: the work to create the cavity and the work to turn on the attractive solute–solvent potential. However, this does not mean additivity of contributions: the attractive solute–solvent potential is switched on given that the cavity has already been created, and the average is calculated by using the conditional probability density function qc(X). The work to create the cavity is: DGc ¼ RT lnhfðXÞip RT ln pinsertion
ð6Þ
It is the reversible work to single out the configurations containing the cavity from the ensemble of pure solvent configurations [14,15] and, according to a theorem of statistical mechanics [16], it is exactly related to the probability that the desired region of space in the liquid is devoid of liquid molecules (i.e., the so-called insertion probability). The work to turn on the attractive solute–solvent potential is: DGa ¼ RT lnhexp½wðXÞ=RT ic
ð7Þ
By setting j = w(X) Æw(X)æc, expanding in power series the exponential function, and keeping in mind that Æjæc ” 0, one obtains [6,7]: DGa ffi hwðXÞic ½hj2 ic =2RT
ð8Þ
if higher order terms are ignored. It is simple to show that the fluctuation term Æj2æc/2RT is of entropic origin [11,12]. When the attractive potential w(X) is weak, the fluctuations in the value of Æw(X)æc are small, and the second term on the right-hand side of Eq. (8) can be neglected. Under this condition, the solvation Gibbs energy change is given by: DG ¼ DGc þ hwðXÞic
ð9Þ
For the hydration of alkanes and noble gases, which are unable to form H-bonds with water, the fluctuation term Æj2æc/2RT is small because w(X) represents the solute– water van der Waals interactions which are weak compared to the water–water H-bonds [17]. In such cases it should be safe to assume that Æw(X)æc is equal to the direct solute– water interaction energy DUsw (i.e., the subscript sw stands for solute–water interaction). When the solute molecule can form H-bonds with water, the condition of small fluctuations for the solute–water interaction energy could not hold [17]. However, in the ensemble of the water configurations possessing a suitable cavity to host the solute, the quantity Æw(X)æc should not account for the solute–water H-bond energy. A negligible
116
G. Graziano / Chemical Physics Letters 429 (2006) 114–118
fraction of water molecules in the qc(X) distribution would have the correct orientation to form a H-bond with the solute molecule inserted in the cavity because water molecules have not yet reorganized in response to the attractive potential of the solute, since the latter acts as a ghost [6,7]. Thus Æw(X)æc should account only for the van der Waals component of the solute–water interaction energy, Æw(X)æc = DUsw(vdW), the fluctuation term is expected to be small, and the contribution of solute–water H-bonds should be included in the solvent reorganization term. By adopting this line, G&L found that [6,7]: (a) the benzene–water van der Waals interaction energy is larger in magnitude than the work of cavity creation and is able to account for the experimental negative DG value; (b) the fluctuation term proves to be small for the hydration of benzene in view of the weakness and negligible anisotropy of the van der Waals component of the benzene–water potential; (c) the reorientation of water molecules to form H-bonds with the aromatic ring is characterized by enthalpy–entropy compensation because the H-bonding potential of benzene is not strong. Further light can be shed on the matter by using an inverse relationship proved by Widom [13], so that the work to turn on the solute–solvent attractive potential is: DGa ¼ RT lnhexp½wðXÞ=RT ic ¼ RT lnhexp½wðXÞ=RT ia
ð10Þ
where the subscript a means that the ensemble average is performed over the solution configurations possessing the cavity occupied by the solute molecule that interacts with the surrounding solvent molecules. The probability density function of such configurations is: Z qa ðXÞ ¼ fðXÞ expf½wðXÞ þ H ðXÞ=RT g= fðXÞ expf½wðXÞ þ H ðXÞ=RT gdX
ð11Þ
In other words, in such ensemble the solute molecule does not act as a ghost but fully interacts with and influences the surrounding solvent molecules. In general, it should be noted that the Widom’s inverse relationship is valid only when the external potential W(X) is not infinite, and that W(X) to insert a cavity at a fixed position is infinite for many solvent configurations if the position is occupied by solvent molecules. This is why cavity creation is a special process and is separated out first, before applying the Widom’s inverse relationship. By setting c = w(X)Æw(X)æa, expanding in power series the exponential function in the last expression of Eq. (10), and keeping in mind that Æcæa” 0, one obtains: 2
DGa ffi hwðXÞia þ ½hc ia =2RT
ð12Þ
if higher order terms are ignored. Even though Eq. (12) looks like Eq. (8), its physical meaning is entirely different: (a) Æw(X)æa accounts for the full solute–water interaction energy, including possible H-bonds, Æw(X)æa = DUsw(full); (b) the fluctuation term Æc2æa/2RT, always of entropic ori-
gin, should be a large quantity for a solute able to form H-bonds with water since the H-bonding potential is strong and is characterized by a marked anisotropy due to the geometric requirements for the occurrence of H-bonds. Actually, for strongly H-bonding solutes such as n-alcohols, the truncation of the power series expansion cannot be a correct procedure [18]. If the fluctuation term Æc2æa/2RT is not a negligible quantity, the solvation Gibbs energy change is given by: DG ¼ DGc þ hwðXÞia þ ½hc2 ia =2RT
ð13Þ
Furthermore, since both Eqs. (8) and (12) are correct from the statistical mechanical point of view, by equating the two expressions, one obtains [6]: hj2 ic =2RT ¼ ½hc2 ia =2RT þ ðhwðXÞia hwðXÞic Þ
ð14Þ
that should allow a reliable evaluation of the fluctuation term Æj2æc/2RT. Now, the numerical values calculated by S&vdV can be inserted in the above equations in order to test the validity of the G&L analysis of benzene solubility in water. 3. Analysis of the S&vdV work Beyond to calculate DG, S&vdV [8] computed directly from the MD simulation runs the benzene–water interaction energy DUsw. It resulted that, at 302 K, DUsw = 59.4 kJ mol1 for the real benzene model and 45.6 kJ mol1 for the vdW benzene model, respectively. S&vdV assumed also that the water reorganization process is characterized by a perfect enthalpy–entropy compensation [19] (i.e., the subscript ww stands for water–water interaction): DH ww ¼ T DS ww
ð15Þ
Therefore DG was solely given by the enthalpic and entropic contributions due to solute–water interactions (note, in general, that in a condensed phase of the matter, such as a liquid, the difference between energy and enthalpy is very small): DG ¼ DU sw T DS sw
ð16Þ
S&vdV, by using Eq. (16), obtained the solute–water entropy contribution at 302 K: T Æ DSsw = DUsw DG = 54.6 kJ mol1 for the real benzene model and 46.6 kJ mol1 for the vdW benzene model, respectively. In order to provide an interpretation for these large and negative values of the benzene–water entropy change, S&vdV used a theoretical relationship derived by Sanchez and colleagues [20]: T DS sw ¼ RT ln pinsertion ½hc2 ia =2RT
ð17Þ
where the first term on the right-hand-side is minus the work of cavity creation, DGc, and the second term on the right-hand-side is a fluctuation term identical to that of Eq. (12) – it accounts for the fluctuation of the attractive solute–solvent potential energy w averaged over the ensemble of solution configurations in which there are no solute–
G. Graziano / Chemical Physics Letters 429 (2006) 114–118
solvent repulsions and the attractive solute–solvent potential is turned on. S&vdV performed MD trajectories 90 ns long for the two models of benzene, and found that the fluctuation term amounts to 12.8 kJ mol1 for the real benzene model and 3.9 kJ mol1 for the vdW benzene model [8]. This means that the work to create in SPC water at 302 K a cavity suitable to host benzene is: (a) DGc = RT Æ ln pinsertion = 54.6 12.8 = 41.8 kJ mol1 for the real benzene model; (b) DGc = RT Æ ln pinsertion = 46.6 3.9 = 42.7 kJ mol1 for the vdW benzene model. It is important to note that these two values should be identical because DGc is a property of the pure liquid and cannot depend on the potential energy of the solute molecule [14,15]; in the following the corresponding mean value, DGc = 42.3 kJ mol1 will be used. The two numbers derived from the simulations of S&vdV are close to each other and are reliable. Indeed (a) from the DGc values calculated by van Gunsteren and colleagues [21] in SPC water at room temperature, DGc 43 kJ mol1 for a cavity suitable to host a benzene molecule; (b) according to scaled particle theory [22], DGc = 41.3 kJ mol1, using the experimental density of water at 25 C, and the hard sphere diameters r(ben˚ and r(water) = 2.80 A ˚ [6,23]. zene) = 5.26 A Clearly, it is now possible to apply Eqs. (9) and (13) to the two benzene models of S&vdV. By using Eq. (9) for the vdW benzene model and considering that Æw(X)æc = DUsw(vdW), it results: DG ffi DGc þ DU sw ðvdWÞ ¼ 42:3 45:6 ¼ 3:3 kJ mol1 ð18Þ By using Eq. (13) for the real benzene model and considering that Æw(X)æa = DUsw(full), it results: DG ¼ DGc þ DU sw ðfullÞ þ ½hc2 ia =2RT ¼ 42:3 59:4 þ 12:8 ¼ 4:3 kJ mol1
ð19Þ
These calculations indicate unequivocally that both Eqs. (9) and (13) are right and provide a similar numerical result and the same physical explanation for the relatively large solubility of benzene in water. Eq. (18) indicates that the experimental negative DG value of benzene in water is accounted for by simply considering the balance between the work of cavity creation and the contribution of van der Waals interactions. But this is the same explanation emerging from Eq. (19), because the favorable enthalpic contribution of the weak benzene–water H-bonds is compensated for by a corresponding unfavorable entropic contribution: DU sw ðfullÞ þ ½hc2 ia =2RT ¼ 59:4 þ 12:8 ¼ 46:6 kJ mol1 DU sw ðvdWÞ ð20Þ The entropic fluctuation term practically cancels the difference in the benzene–water interaction energy calculated for the two different benzene models by S&vdV. This is exactly the manifestation of the enthalpy–entropy compensation characterizing the weak benzene–water H-bonds,
117
as indicated by G&L [6]. The magnitude of the fluctuation term Æc2æa/2RT expresses the extent to which the attractive solute–solvent interactions bias positions and orientations of the vicinal solvent molecules, causing a reduction of the accessible configuration space that, in turn, leads to an entropy decrease. It is also possible to obtain an estimate of the fluctuation term Æj2æc/2RT in Eq. (8) using Eq. (14) and the values calculated by S&vdV: hj2 ic =2RT ffi ½hc2 ia =2RT þ DU sw ðfullÞ DU sw ðvdWÞ ¼ 12:8 59:4 þ 45:6 ¼ 1:0 kJ mol1
ð21Þ
This estimate indicates that such fluctuation quantity should be small [6]. It is also important to note that Æj2æc/2RT is a negative quantity implying that the reorganization of water molecules to form the weak H-bonds with the benzene ring is a spontaneous process, even though characterized by an almost complete enthalpy–entropy compensation. S&vdV determined the fluctuation term also for the vdW benzene model and the value was 3.9 kJ mol1, that cannot be considered negligible. However, according to the statistical-mechanical approach of G&L, the quantity Æj2æc/2RT has to be calculated over the ensemble of solvent configurations possessing the suitable cavity, which is different from the ensemble of solution configurations used by S&vdV to calculate the fluctuation term for the vdW benzene model. 4. Discussion The performed analysis indicates that the values obtained from the MD simulations of S&vdV are not in contrast with the G&L approach: the benzene–water van der Waals interaction energy is larger in magnitude than the work of cavity creation and is able to account for the experimental negative DG value. It is worth noting that the van der Waals benzene–water interactions are stronger than those of an aliphatic hydrocarbon of similar molecular size with water, due to the larger polarizability of the aromatic ring [7]. Consider, for instance, that, even though ˚ 3, is similar to that of the polarizability of c-hexane, 10.78 A 3 ˚ benzene, 10.32 A [22], the c-hexane molecule is significantly larger than that of benzene. The effective hard ˚ for c-hexane and 5.26 A ˚ for sphere diameters are 5.63 A benzene [12], so that the molecular volume of the latter is 82% of that of the former. Such a value should correspond also to the polarizability ratio because this physical quantity is proportional to molecular volume. However, experimental data indicate that the polarizability of benzene is 96% of that of c-hexane [22]. In addition, van Gunsteren and colleagues [24] obtained, at room temperature, DUsw = 35.0 kJ mol1 for i-butane and 36.4 kJ mol1 for n-butane; these values are about 10 kJ mol1 smaller in magnitude than the benzene–water van der Waals interaction energy, DUsw(vdW) = 45.6 kJ mol1. This finding
118
G. Graziano / Chemical Physics Letters 429 (2006) 114–118
is in line with the analysis by G&L [6,7], but was not underscored by S&vdV. A further point should be considered: if the contribution of the weak benzene–water H-bonds is small due to enthalpy–entropy compensation, the thermodynamic integration procedure should produce the same DG estimate for the two models of benzene devised by S&vdV. Indeed, thermodynamic integration takes into account solely the effective contributions to the DG quantity. The difference found by S&vdV, DG = 4.8 kJ mol1 for the real benzene model and+1.0 kJ mol1 for the vdW benzene model, is not small, but cannot be considered to be large. In order to gain perspective, one has to remind that, for an alkane having a size comparable to that of benzene (i.e., i-butane and nbutane), DG amounts to 9–10 kJ mol1 at 25 C [2,3]. The latter numbers indicate that the vdW benzene model of S&vdV is markedly less hydrophobic than an alkane of similar size. This result was neglected by S&vdV, even though it emerged directly from their MD simulations. In addition, S&vdV showed that a change of GROMOS parameters produces a marked effect on the DG value calculated for the real benzene model: DG = 4.8 kJ mol1 for the first-set, and 6.7 kJ mol1 for the second-set [8]. S&vdV did not show the DG value obtained for the vdW benzene model using the second-set of GROMOS parameters. One could speculate that, by using the second-set of GROMOS parameters, DG would be a negative quantity also for the vdW benzene model. In conclusion, the performed analysis points out that the numerical values obtained from MD simulations by S&vdV indicate that: (a) the benzene–water van der Waals interaction energy overwhelms the work of cavity creation, accounting for the experimental negative DG value; (b) the weak benzene–water H-bonds do not really influence the benzene solubility in water because the favorable enthalpic component is almost entirely compensated for by an unfavorable entropic component.
Acknowledgement I thank Dr. B. Lee (Center for Cancer Research, NCI, NIH, Bethesda, MD) for reading earlier drafts of the manuscript and providing useful suggestions. References [1] C. McAuliffe, J. Phys. Chem. 70 (1966) 1267. [2] A. Ben-Naim, Solvation Thermodynamics, Plenum, New York, 1987. [3] G.I. Makhatadze, P.L. Privalov, Biophys. Chem. 50 (1994) 285. [4] S. Suzuki, P.G. Green, R.E. Bumgarner, S. Dasgupta, W.A. Goddard, G.A. Blake, Science 257 (1992) 942. [5] P. Linse, J. Am. Chem. Soc. 112 (1990) 1744. [6] G. Graziano, B. Lee, J. Phys. Chem. B 105 (2001) 10367. [7] G. Graziano, Biophys. Chem. 110 (2004) 249. [8] P. Schravendijk, N.F.A. van der Vegt, J. Chem. Theory Comput. 1 (2005) 643. [9] H.J.C. Berendsen, W.F. van Gunsteren, J.P.M. Postma, J. Hermans, in: B. Pullman (Ed.), Intermolecular Forces, Reidel, Dordrecht, 1981, p. 331. [10] H.J.C. Berendsen, D. van der Spoel, R. van Drunen, Comput. Phys. Commun. 91 (1995) 43. [11] B. Lee, Biopolymers 31 (1991) 993. [12] G. Graziano, Can. J. Chem. 80 (2002) 401. [13] B. Widom, J. Phys. Chem. 86 (1982) 869. [14] B. Lee, J. Chem. Phys. 83 (1985) 2421. [15] G. Graziano, J. Phys. Soc. Jpn. 69 (2000) 1566. [16] R.C. Tolman, The Principles of Statistical Mechanics, Oxford University Press, London, 1938. [17] B. Lee, Biophys. Chem. 51 (1994) 271. [18] G. Graziano, J. Phys. Chem. B 109 (2005) 12160. [19] H.A. Yu, M. Karplus, J. Chem. Phys. 89 (1988) 2366. [20] I.C. Sanchez, T.M. Truskett, P.J. In’t Veld, J. Phys. Chem. B 103 (1999) 5106. [21] T.C. Beutler, D.R. Beguelin, W.F. van Gunsteren, J. Chem. Phys. 102 (1995) 3787. [22] R.A. Pierotti, Chem. Rev. 76 (1976) 717. [23] G. Graziano, Chem. Phys. Lett. 396 (2004) 226. [24] D. Trzesniak, N.F.A. van der Vegt, W.F. van Gunsteren, Phys. Chem. Chem. Phys. 6 (2004) 697.