- Email: [email protected]
An approach to the solvation free energy in terms of the distribution functions of the solute–solvent interaction energy
Journal of Molecular Liquids 119 (2005) 23 – 29 www.elsevier.com/locate/molliq
An approach to the solvation free energy in terms of the distribution functions of the solute–solvent interaction energy Nobuyuki Matubayasi*, Masaru Nakahara Institute for Chemical Research, Kyoto University, Uji, Kyoto, 611-0011, Japan Available online 23 November 2004
Abstract The energy representation of the molecular configuration in a dilute solution is introduced to express the solvent distribution around the solute over a one-dimensional coordinate specifying the solute–solvent interaction energy. On the basis of the energy representation, an approximate functional for the solvation free energy of a solute in solution is constructed by adopting the Percus-Yevick-type approximation in the unfavorable region of the solute–solvent interaction and the hypernetted-chain-type approximation in the favorable region. The solvation free energy is then given exactly to second order with respect to the solvent density and to the solute–solvent interaction. It is demonstrated that the solvation free energies of nonpolar, polar, and ionic solutes in water are evaluated accurately and efficiently from the single functional over a wide range of thermodynamic conditions. The extension to a flexible solute molecule is straightforward. The applicability of the method is illustrated for solute molecules with a stretching or torsional degree of freedom. D 2004 Elsevier B.V. All rights reserved. Keywords: Solvation free energy; Energy representation; Distribution function; Computer simulation
1. Introduction The most fundamental quantity to describe a process in solution is the free energy change. Indeed, it governs the equilibrium and rate constants of the process. The free energy change corresponding to the insertion process of a solute in solution is the chemical potential or the solvation free energy. Once the chemical potentials are known for the species present in the initial and final states of a process of interest, the free energy change for the process can be readily evaluated. Therefore, it is of primary importance in statistical mechanics of solutions to establish a scheme to determine the solvation free energy of a solute in solution. A molecular description of the solvation free energy can be implemented by formulating a functional which expresses the solvation free energy in terms only of distribution functions in the solution and pure solvent systems of interest. An approximate functional need to be
* Corresponding author. Tel.: +81 774 38 3073; fax: +81 774 38 3076. E-mail address: [email protected] (N. Matubayasi). 0167-7322/$ - see front matter D 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.molliq.2004.10.005
constructed in practice, however, since the exact functional involves an infinite series of many-body distribution functions [1]. When the distribution functions constituting the approximate functional can be obtained by computer simulation, the solvation free energy is readily determined by utilizing the exact structure of the solution. Therefore, it is desirable to develop a simple and accurate functional for the solvation free energy with a well-defined approximation. When a set of potential functions is given for the solution system of interest, the bexactQ solvation free energy of the solute molecule can be calculated by the free energy perturbation and thermodynamic integration methods [2,3]. These methods are, however, expensive in computation and practical only for small molecules. For the purpose of analyzing the solvation free energy on the molecular level, furthermore, it is necessary to express the solvation free energy in terms of distribution functions of the systems of interest. Thus, an accurate functional for the solvation free energy not only provides an efficient route of computation, but also sets a basis for the molecular understanding. In the present paper, we introduce a functional for the solvation free energy of a solute in solution. We adopt the
24
N. Matubayasi, M. Nakahara / Journal of Molecular Liquids 119 (2005) 23–29
energy representation introduced in Refs. [4–6], and provide an approximate but accurate functional which is expressed in terms of energy distribution functions in the solution and pure solvent systems of interest. When the energy representation is adopted, the coordinate of a solvent molecule around the solute molecule is the solute–solvent interaction energy and the solvent distribution around the solute is expressed over the one-dimensional coordinate for any type of solute–solvent interaction potential. In the energy representation, a straightforward construction of a functional for the solvation free energy is possible which is exact to second order in the solvent density and in the solute–solvent interaction. The method is thus suitable to study supercritical fluid. It describes the low- to medium-density region accurately and is useful to explore a wide range of thermodynamic conditions. The method of energy representation is also pertinent to the treatment of flexible molecules. Indeed, the energy distribution functions constituting the functional for the solvation free energy can be defined without referring to whether the molecules are structurally rigid or flexible. The organization of the present paper is as follows. In Section 2, a set of distribution functions are introduced to formulate a reduced description of solvent configuration relative to the solute molecule, and an approximate functional is presented for the solvation free energy. In Section 3, the procedures of computation are described. In Section 4, the performance of the functional for the solvation free energy is assessed with solvent water over a wide range of thermodynamic conditions for nonpolar, polar, and ionic solutes with and without structural flexibility.
2. Theory The system of our interest is a dilute solution containing a single solute molecule. The intermolecular interaction is supposed to be pairwise additive. The notations and developments in this paper are then parallel to those adopted in Refs. [4–6]. The complete set of the position and orientation of a solvent molecule is called the full coordinate and is denoted collectively by x. When the solute molecule involves structural flexibility, its intramolecular degrees of freedom are collectively written as w. Of course, when the solute is rigid, the following developments are valid simply by omitting the dependence on w of any function. The solute–solvent interaction potential of interest is v and is fixed at the outset in our developments. It may be expressed as v f (w,x) in the full coordinate representation, where a superscript f is attached to emphasize that a function is represented over the full coordinate. The solvation free energy Dl is the free energy change corresponding to the gradual insertion process of the solute molecule. In Dl, only the contribution from the potential energy is involved and the ideal (kinetic) contribution is excluded. When the intramolecular energy of the solute is
C(w) and the total solvent–solvent interaction energy is U(X), Dl is expressed as expð bDlÞ Z
( d wdX exp b WðwÞþ Z
¼
X
)! f
v ðw; xi Þ þ U ðXÞ
i
d wdX expð bfWðwÞ þ U ðXÞgÞ ð1Þ
where X represents the solvent configuration collectively and xi is the full coordinate of the ith solvent molecule. b is the inverse of k BT, as in the usual notational convention, with the Boltzmann constant k B and the temperature T. A restriction of attention to a certain set of solute intramolecular state can be made simply by the corresponding alteration of the domain of integration over w. The energy representation is introduced by adopting the value of the solute–solvent interaction v of interest as the coordinate e for the distribution of the solvent molecule around the solute molecule. The instantaneous distribution qˆe is then defined as X d v f ðw; xi Þ e ð2Þ qˆ e ðeÞ ¼ i
where the sum is taken over the solvent molecules and a superscript e is attached to emphasize that a function is represented over the energy coordinate. Note that the specification of v is necessary in Eq. (2). In the energy representation, the solvation free energy Dl is approximately expressed in terms of distribution functions constructed from qˆ e in the solution and pure solvent systems. In our treatments, the solution system refers to the system in which the solute molecule interacts with the solvent under the solute–solvent interaction v of interest at full coupling. In the solution, the average distribution q e of the v value is relevant in the approximate construction of Dl and is given by qe ðeÞ ¼ hqˆ e ðeÞi
ð3Þ
where h...i represents the ensemble average in the solution system of interest. On the other hand, the pure solvent system denotes the system in which no interaction is physically present between the solute and solvent molecules. At an instantaneous configuration of the pure solvent system, qˆ e is constructed by placing the solute molecule in the neat solvent system as a test particle. The average ˆ distribution q 0e and the correlation matrix v 0e then appear in the approximate functional for Dl and are expressed, respectively, as qe0 ðeÞ ¼ hqˆ e ðeÞi0
ð4Þ
and ve0 ðe; gÞ ¼ hqˆ e ðeÞqˆ e ðgÞi0 hqˆ e ðeÞi0 hqˆ e ðgÞi0
ð5Þ
N. Matubayasi, M. Nakahara / Journal of Molecular Liquids 119 (2005) 23–29
where h...i0 represents the ensemble average in the pure solvent system. In this case, the solute molecule is placed as a test particle in the neat solvent system and the solute and solvent degrees of freedom are uncoupled from each other in the probability distribution. An approximate functional for Dl is derived in Ref. [5]. The functional is constructed by adopting the PercusYevick-type approximation in the unfavorable region of the solute–solvent interaction and the hypernetted-chaintype approximation in the favorable region. Dl is then given by a set of equations listed as
qe ðeÞ w ðeÞ ¼ kB T log e q0 ðeÞ e
e
ð6Þ
Z
dðe gÞ e 1 v0 ðe; gÞ qe0 ðeÞ
e e
q ðgÞ q0 ðgÞ
we0 ðeÞ ¼ kB T
Dl ¼ kB T
dg
Z de
e
q ðeÞ
qe0 ðeÞ
ð7Þ
þ bwe ðeÞqe ðeÞ
e
faðeÞF ðeÞ þ ð1 aðeÞÞF0 ðeÞg q ðeÞ
qe0 ðeÞ
ð8Þ
F ðeÞ ¼
8 > < bwe ðeÞ þ 1 þ
e
bw ðeÞ expð bwe ðeÞÞ 1
> : 1 bwe ðeÞ 2
ðwhen we ðeÞV0Þ ðwhen we ðeÞz0Þ
ð9Þ 8 log 1 bwe0 ðeÞ > > < log 1 bwe0 ðeÞ þ 1 þ e bw0 ðeÞ F0 ðeÞ ¼ > > 1 bwe ðeÞ : 2 0
when we0 ðeÞV0
when we0 ðeÞz0
ð10Þ
að e Þ ¼
8 <1 :1
qe ðeÞ qe0 ðeÞ qe ðeÞ þ qe0 ðeÞ
2
when qe ðeÞzqe0 ðeÞ when qe ðeÞVqe0 ðeÞ ð11Þ
3. Procedures In the present work, the solvent is water. The solution is dilute and contains a single solute molecule. The thermodynamic state of each system of interest is specified by the water density and temperature. Four thermodynamic states are examined. One is an ambient state of 1.0 g/cm3 and 25 8C and the others are supercritical states of 1.0, 0.6, and 0.2 g/cm3 and 400 8C. The water molecule is treated as rigid and nonpolarizable, and the SPC/E model is adopted as the intermolecular potential function between water molecules [7].
25
All the nonpolar, polar, and ionic species are examined for rigid solutes. The nonpolar solutes treated are methane and ethane. Water, methanol, and ethanol are employed as polar solutes, and Na+ and Cl as ionic solutes. The interaction potential and the scheme for calculating the solvation free energy are described in Ref. [5]. The flexible solute molecules employed in the present work are comprised of two or four interaction sites. The two- and four-site molecules involve stretching and torsional degree of freedom, respectively, for intramolecular motion. The potential function for the stretching motion of the two-site solute is harmonic and is given by K ðr re Þ2
ð12Þ
where r is the distance between the two sites and r e is taken to be 3.0 2 [6]. Jorgensen’s parameter set for nbutane is adopted to model the four-site solute [8]. The bond lengths and angles are fixed at the values listed in Table 1 of Ref. [8] and the torsional degree of freedom is described by the potential energy expressed as Eq. (2) and Table 2 of Ref. [8]. The intermolecular interaction between the flexible solute and water molecules consists of the Lennard-Jones and Coulombic terms, as usual, and is given by Eq. (23) of Ref. [5] with the truncation factor S(x) set to unity. The Lennard-Jones parameters for the two-site molecule are taken from Ref. [9]. The values for the sodium ion in Table 1 of Ref. [9] are assigned to one of the sites, and those for the chloride ion to the other. For the four-site molecule, Jorgensen’s set for n-butane in Table 3 of Ref. [8] is employed as the Lennard-Jones parameters. The LennardJones parts of the solute-water potential functions are then constructed by the standard Lorentz-Berthelot combining rule [2]. For the Coulombic interaction, two cases are examined. One is the nonpolar case, in which no charge is given to any of the sites. In the other case, the charge of +1 in the unit of elementary charge is carried by the Na+-like site of the two-site molecule and by one of the methyl-like sites of the four-site molecule. The Cl like site of the twosite molecule and the other methyl-like site of the four-site molecule involve the corresponding negative charges, so that the molecules are neutral in total. In the rest of the paper, the two-site molecule with K=2 kcal/mol/22 in Eq. (12) is called HM-2-0 and HM-2-1, respectively, when the charge on the Na+-like site is 0 and +1 in the unit of elementary charge. Similarly, the two-site molecule with K=5 kcal/mol/22 is denoted as HM-5-0 and HM-5-1, respectively, at the charge of 0 and +1. The four-site molecule is called Bu-0 and Bu-1, respectively, when one of the methyl-like sites has a charge of 0 and +1. In any case, no charge is placed on the nonterminal methylene-like sites of the four-site solutes. To obtain the energy distribution function q e, a Monte Carlo simulation was carried out for the solution system of interest. The distribution functions q 0e and v 0e were also
26
N. Matubayasi, M. Nakahara / Journal of Molecular Liquids 119 (2005) 23–29
calculated from a simulation of the pure solvent system. In this case, the solute molecule was inserted as a test particle at a random position with a random orientation. The intramolecular degree of freedom of the solute molecule was sampled simply through the intramolecular potential (or through the modified scheme given in Ref. [6] when the modified form of functional in Ref. [6] was employed for
the solvation free energy). It should be noted that when the solute molecule is inserted into the pure solvent system, it often overlaps with solvent molecules. The overlapping configurations contribute to q 0e and v 0e at large energy coordinate and account for the excluded volume effect in the solvation free energy. The system size and the simulation lengths correspond to those in Refs. [5,6].
Fig. 1. Graphical comparison of the solvation free energy Dl for the rigid solute molecules. The approximate values are provided by Eqs. (6)–(11) and the exact values are obtained from the free energy perturbation calculations in Ref. [5].
N. Matubayasi, M. Nakahara / Journal of Molecular Liquids 119 (2005) 23–29
27
4. Results and discussion In Fig. 1, the approximate values of the solvation free energy Dl for the rigid solute molecules are graphically compared to the corresponding exact values obtained from Table 1 Solvation free energies in the unit of kcal/mola Solute
Thermodynamic Exact state
HM-2-0 1.0 g/cm3 and 25 8C 1.0 g/cm3 and 400 8C 0.6 g/cm3 and 400 8C 0.2 g/cm3 and 400 8C HM-5-0 1.0 g/cm3 and 25 8C 1.0 g/cm3 and 400 8C 0.6 g/cm3 and 400 8C 0.2 g/cm3 and 400 8C HM-2-1 1.0 g/cm3 and 25 8C 1.0 g/cm3 and 400 8C 0.6 g/cm3 and 400 8C 0.2 g/cm3 and 400 8C HM-5-1 1.0 g/cm3 and 25 8C 1.0 g/cm3 and 400 8C 0.6 g/cm3 and 400 8C 0.2 g/cm3 and 400 8C Bu-0 1.0 g/cm3 and 25 8C 1.0 g/cm3 and 400 8C 0.6 g/cm3 and 400 8C 0.2 g/cm3 and 400 8C Bu-1 1.0 g/cm3 and 25 8C 1.0 g/cm3 and 400 8C 0.6 g/cm3 and 400 8C 0.2 g/cm3 and 400 8C
Approximate Eqs. (6)–(11) Modified formb
3.2F0.2
3.5F0.1
3.4F0.1
15.1F0.3
12.4F0.1
12.5
3.5F0.1
3.9F0.1
3.9
0.7F0.1
0.7
0.7
3.3F0.2
3.6F0.1
3.4F0.2
15.3F0.3
12.6F0.1
12.6F0.1
3.5F0.1
3.9
3.9
0.7F0.1
0.7
0.7
104.8F0.9
–c
108.5F1.8
86.2F0.5
–c
87.2F0.5
86.4F0.3
–c
82.9F0.4
78.6F1.3
–c
70.9F0.4
92.5F0.7
–c
97.5F1.7
74.1F0.7
–c
74.3F0.2
74.7F0.5
–c
70.5F0.3
67.4F1.0
–c
60.8F0.7
3.1F0.4
4.4F0.2
4.5F0.3
21.6F0.2
17.9F0.1
17.9F0.1
3.8F0.2
4.8F0.1
4.8F0.1
0.7F0.2
0.6
0.6
63.5F0.7 63.2F1.0
63.3F1.0
41.1F0.7 43.9F0.3
44.7F0.3
47.2F0.4 46.4F0.1
47.2F0.1
39.5F0.4 38.4F0.4
38.7F0.4
a Each value is rounded to a multiple of 0.1 kcal/mol. The error is expressed as 95% confidence level and is not shown when it is rounded to 0. b A modified form of approximation given by Eqs. (29) and (35)–(41) of Ref. [6]. c The value calculated is numerically unstable or involves a large discretization error typically of a few kcal/mol.
Fig. 2. The probability distribution functions p in solution and p 0 at isolation of the dihedral angle / of the Bu-0 solute at an ambient state of 1.0 g/cm3 and 25 8C and at supercritical states of 1.0, 0.6, and 0.2 g/cm3 and 400 8C. The solute is at the trans conformation when /=p, and it is at the gauche when /=p/3 or 5p/3. At convergence, p and p 0 are symmetric against the reflection with respect to /=p. The asymmetry is present due to the statistical error. It should be noted that p 0 is common to the three supercritical states.
the standard free energy perturbation method. The good agreement is observed between the approximate and exact values. The agreement is particularly notable at the low- and medium-density states of 0.2 and 0.6 g/cm3 and 400 8C. When the solute is ionic, the density at the state of 0.2 g/cm3 and 400 8C is not yet low in the sense, for example, that the hydration number at that state is comparable to the numbers at ambient states [10,11]. Even in this case, our approximate procedure is effective in determining the solvation free energy. The solvation free energies of water at 1.0 g/cm3 and 400 8C and of methanol and ethanol at 0.6 g/cm3 and 400 8C are rather small in magnitude. These behaviors are caused by the balance between the favorable and unfavorable contributions of the solute–solvent interactions, and are well reproduced by our approximate method. Therefore, the single functional given by Eqs. (6)–(11) provides an accurate and efficient route to the solvation free energy for various types of solutes over a wide range of thermodynamic conditions. The Dl values for the flexible solute molecules are presented in Table 1. As noted in Ref. [6], the approximate expression listed as Eqs. (6)–(11) is not successful for the HM-2-1 and HM-5-1 solutes. Their intramolecular structures are largely changed through the solute–solvent interaction, and the modified scheme of approximation described in Ref. [6] is needed for good performance. The approximate Dl calculated by Eqs. (6)–(11) are in good agreement with the exact for the HM-2-0, HM-5-0, Bu-0,
28
N. Matubayasi, M. Nakahara / Journal of Molecular Liquids 119 (2005) 23–29
Fig. 3. The probability distribution functions p in solution and p 0 at isolation of the dihedral angle / of the Bu-1 solute at an ambient state of 1.0 g/cm3 and 25 8C and at supercritical states of 1.0, 0.6, and 0.2 g/cm3 and 400 8C. The solute is at the trans conformation when /=p, and it is at the gauche when /=p/3 or 5p/3. At convergence, p and p 0 are symmetric against the reflection with respect to /=p. The asymmetry is present due to the statistical error. It should be noted that when the Ewald method is employed, p 0 at a fixed temperature depends slightly on the density due to the interaction of the solute with its own images.
and Bu-1 solutes. For the HM-2-0 and HM-5-0 solutes, it is shown that the solvent effect on the solute intramolecular structure is weak. In this case, the modified form of approximation provided in Ref. [6] does not lead to an improvement and the original form given by Eqs. (6)–(11) performs well. In Table 1, the Dl values calculated by the modified form of approximation in Ref. [6] are also listed for the Bu-0 and Bu-1 solutes. They are observed to be coincident with the corresponding values from Eqs. (6)–(11). To see the origin of this observation, we examine the probability distribution function of the torsional angle / of the four-site molecule. In the following, the probability distribution function of the intramolecular coordinate / is denoted by p in solution and by p 0 at isolation (absence of the solute– solvent interaction). In Fig. 2, we show p and p 0 for the Bu-0 solute. Actually, the Bu-0 solute models n-butane and is often employed as a prototypic solute to study the hydrophobic effect [12,13]. As expected from the common notion of hydrophobicity, it is seen at each thermodynamic state that the trans conformation becomes less populated through the solute–solvent interaction. In other words, a compact conformation is preferable in water over a wide range of thermodynamic conditions. When the temperature is elevated from 25 to 400 8C at a fixed density of 1.0 g/cm 3 , p/p 0 at the typical trans and gauche conformations is observed to vary within 20%. Actually,
this observation is consistent with a previous finding that the solvent effect on the association and dissociation of a pair of nonpolar solutes is weakly dependent on the temperature when the density of solvent water is fixed [14]. The density reduction at a fixed supercritical temperature of 400 8C weakens the effect of solvation. At a low density of 0.2 g/cm3, the solvent effect is barely appreciable. In Fig. 3, p and p 0 are shown for the polar four-site solute, Bu-1. The solute molecule is extended due to the solute–solvent interaction, and the trans conformation is exclusively occupied in the solution systems examined. This observation is parallel to that for HM-2-1 and HM-5-1 reported in Ref. [6]. While p is broader at the higher temperature, the solvent effect is still bstrongQ even at a lowdensity state of 0.2 g/cm3 and 400 8C in the sense that the isothermal density reduction at 400 8C from 1.0 to 0.2 g/cm3 causes a weaker effect on p than that from 0.2 g/cm3 to the zero density limit. It is seen in Figs. 2 and 3 that the main domain for p overlaps significantly with the peak domain for p 0. In other words, the typical conformations of the solute in solution are frequently realized for the solute at isolation. In this case, the modification provided in Ref. [6] does not lead to an improvement of the approximate scheme, in agreement with the arguments presented in Ref. [6].
Acknowledgments This work is supported by the Grant-in-Aid for Scientific Research (Nos. 11740322, 13640509, and 15205004) from the Japan Society for the Promotion of Science and by the Grant-in-Aid for Scientific Research on Priority Areas (No. 15076205) and the NAREGI (National Research Grid Initiative) Project from the Ministry of Education, Culture, Sports, Science, and Technology. N.M. is also grateful to the Supercomputer Laboratory of Institute for Chemical Research, Kyoto University, for generous allocation of computation time.
References [1] N. Matubayasi, L.H. Reed, R.M. Levy, J. Phys. Chem. 98 (1994) 10640. [2] J.P. Hansen, I.R. McDonald, Theory of Simple Liquids, Academic Press, London, 1986. [3] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Oxford University Press, Oxford, 1987. [4] N. Matubayasi, M. Nakahara, J. Chem. Phys. 113 (2000) 6070. [5] N. Matubayasi, M. Nakahara, J. Chem. Phys. 117 (2002) 3605, 118 (2003) 2446. [6] N. Matubayasi, M. Nakahara, J. Chem. Phys. 119 (2003) 9686. [7] H.J.C. Berendsen, J.R. Grigera, T.P. Straatsma, J. Phys. Chem. 91 (1987) 6269. [8] W.L. Jorgensen, J.D. Madura, C.J. Swenson, J. Am. Chem. Soc. 106 (1984) 6638.
N. Matubayasi, M. Nakahara / Journal of Molecular Liquids 119 (2005) 23–29 [9] J. Chandrasekhar, D.C. Spellmeyer, W.L. Jorgensen, J. Am. Chem. Soc. 106 (1984) 903. [10] P.B. Balbuena, K.P. Johnston, P.J. Rossky, J. Phys. Chem. 99 (1995) 1554. [11] P.B. Balbuena, K.P. Johnston, P.J. Rossky, J. Phys. Chem. 100 (1996) 2706.
29
[12] R.O. Rosenberg, R. Mikkilineni, B.J. Berne, J. Am. Chem. Soc. 104 (1982) 7647. [13] W.L. Jorgensen, J. Gao, C. Ravimohan, J. Phys. Chem. 89 (1985) 3470. [14] N. Matubayasi, M. Nakahara, J. Phys. Chem., B 104 (2000) 10352.
An approach to the solvation free energy in terms of the distribution functions of the solute–solvent interaction energy Nobuyuki Matubayasi*, Masaru Nakahara Institute for Chemical Research, Kyoto University, Uji, Kyoto, 611-0011, Japan Available online 23 November 2004
Abstract The energy representation of the molecular configuration in a dilute solution is introduced to express the solvent distribution around the solute over a one-dimensional coordinate specifying the solute–solvent interaction energy. On the basis of the energy representation, an approximate functional for the solvation free energy of a solute in solution is constructed by adopting the Percus-Yevick-type approximation in the unfavorable region of the solute–solvent interaction and the hypernetted-chain-type approximation in the favorable region. The solvation free energy is then given exactly to second order with respect to the solvent density and to the solute–solvent interaction. It is demonstrated that the solvation free energies of nonpolar, polar, and ionic solutes in water are evaluated accurately and efficiently from the single functional over a wide range of thermodynamic conditions. The extension to a flexible solute molecule is straightforward. The applicability of the method is illustrated for solute molecules with a stretching or torsional degree of freedom. D 2004 Elsevier B.V. All rights reserved. Keywords: Solvation free energy; Energy representation; Distribution function; Computer simulation
1. Introduction The most fundamental quantity to describe a process in solution is the free energy change. Indeed, it governs the equilibrium and rate constants of the process. The free energy change corresponding to the insertion process of a solute in solution is the chemical potential or the solvation free energy. Once the chemical potentials are known for the species present in the initial and final states of a process of interest, the free energy change for the process can be readily evaluated. Therefore, it is of primary importance in statistical mechanics of solutions to establish a scheme to determine the solvation free energy of a solute in solution. A molecular description of the solvation free energy can be implemented by formulating a functional which expresses the solvation free energy in terms only of distribution functions in the solution and pure solvent systems of interest. An approximate functional need to be
* Corresponding author. Tel.: +81 774 38 3073; fax: +81 774 38 3076. E-mail address: [email protected] (N. Matubayasi). 0167-7322/$ - see front matter D 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.molliq.2004.10.005
constructed in practice, however, since the exact functional involves an infinite series of many-body distribution functions [1]. When the distribution functions constituting the approximate functional can be obtained by computer simulation, the solvation free energy is readily determined by utilizing the exact structure of the solution. Therefore, it is desirable to develop a simple and accurate functional for the solvation free energy with a well-defined approximation. When a set of potential functions is given for the solution system of interest, the bexactQ solvation free energy of the solute molecule can be calculated by the free energy perturbation and thermodynamic integration methods [2,3]. These methods are, however, expensive in computation and practical only for small molecules. For the purpose of analyzing the solvation free energy on the molecular level, furthermore, it is necessary to express the solvation free energy in terms of distribution functions of the systems of interest. Thus, an accurate functional for the solvation free energy not only provides an efficient route of computation, but also sets a basis for the molecular understanding. In the present paper, we introduce a functional for the solvation free energy of a solute in solution. We adopt the
24
N. Matubayasi, M. Nakahara / Journal of Molecular Liquids 119 (2005) 23–29
energy representation introduced in Refs. [4–6], and provide an approximate but accurate functional which is expressed in terms of energy distribution functions in the solution and pure solvent systems of interest. When the energy representation is adopted, the coordinate of a solvent molecule around the solute molecule is the solute–solvent interaction energy and the solvent distribution around the solute is expressed over the one-dimensional coordinate for any type of solute–solvent interaction potential. In the energy representation, a straightforward construction of a functional for the solvation free energy is possible which is exact to second order in the solvent density and in the solute–solvent interaction. The method is thus suitable to study supercritical fluid. It describes the low- to medium-density region accurately and is useful to explore a wide range of thermodynamic conditions. The method of energy representation is also pertinent to the treatment of flexible molecules. Indeed, the energy distribution functions constituting the functional for the solvation free energy can be defined without referring to whether the molecules are structurally rigid or flexible. The organization of the present paper is as follows. In Section 2, a set of distribution functions are introduced to formulate a reduced description of solvent configuration relative to the solute molecule, and an approximate functional is presented for the solvation free energy. In Section 3, the procedures of computation are described. In Section 4, the performance of the functional for the solvation free energy is assessed with solvent water over a wide range of thermodynamic conditions for nonpolar, polar, and ionic solutes with and without structural flexibility.
2. Theory The system of our interest is a dilute solution containing a single solute molecule. The intermolecular interaction is supposed to be pairwise additive. The notations and developments in this paper are then parallel to those adopted in Refs. [4–6]. The complete set of the position and orientation of a solvent molecule is called the full coordinate and is denoted collectively by x. When the solute molecule involves structural flexibility, its intramolecular degrees of freedom are collectively written as w. Of course, when the solute is rigid, the following developments are valid simply by omitting the dependence on w of any function. The solute–solvent interaction potential of interest is v and is fixed at the outset in our developments. It may be expressed as v f (w,x) in the full coordinate representation, where a superscript f is attached to emphasize that a function is represented over the full coordinate. The solvation free energy Dl is the free energy change corresponding to the gradual insertion process of the solute molecule. In Dl, only the contribution from the potential energy is involved and the ideal (kinetic) contribution is excluded. When the intramolecular energy of the solute is
C(w) and the total solvent–solvent interaction energy is U(X), Dl is expressed as expð bDlÞ Z
( d wdX exp b WðwÞþ Z
¼
X
)! f
v ðw; xi Þ þ U ðXÞ
i
d wdX expð bfWðwÞ þ U ðXÞgÞ ð1Þ
where X represents the solvent configuration collectively and xi is the full coordinate of the ith solvent molecule. b is the inverse of k BT, as in the usual notational convention, with the Boltzmann constant k B and the temperature T. A restriction of attention to a certain set of solute intramolecular state can be made simply by the corresponding alteration of the domain of integration over w. The energy representation is introduced by adopting the value of the solute–solvent interaction v of interest as the coordinate e for the distribution of the solvent molecule around the solute molecule. The instantaneous distribution qˆe is then defined as X d v f ðw; xi Þ e ð2Þ qˆ e ðeÞ ¼ i
where the sum is taken over the solvent molecules and a superscript e is attached to emphasize that a function is represented over the energy coordinate. Note that the specification of v is necessary in Eq. (2). In the energy representation, the solvation free energy Dl is approximately expressed in terms of distribution functions constructed from qˆ e in the solution and pure solvent systems. In our treatments, the solution system refers to the system in which the solute molecule interacts with the solvent under the solute–solvent interaction v of interest at full coupling. In the solution, the average distribution q e of the v value is relevant in the approximate construction of Dl and is given by qe ðeÞ ¼ hqˆ e ðeÞi
ð3Þ
where h...i represents the ensemble average in the solution system of interest. On the other hand, the pure solvent system denotes the system in which no interaction is physically present between the solute and solvent molecules. At an instantaneous configuration of the pure solvent system, qˆ e is constructed by placing the solute molecule in the neat solvent system as a test particle. The average ˆ distribution q 0e and the correlation matrix v 0e then appear in the approximate functional for Dl and are expressed, respectively, as qe0 ðeÞ ¼ hqˆ e ðeÞi0
ð4Þ
and ve0 ðe; gÞ ¼ hqˆ e ðeÞqˆ e ðgÞi0 hqˆ e ðeÞi0 hqˆ e ðgÞi0
ð5Þ
N. Matubayasi, M. Nakahara / Journal of Molecular Liquids 119 (2005) 23–29
where h...i0 represents the ensemble average in the pure solvent system. In this case, the solute molecule is placed as a test particle in the neat solvent system and the solute and solvent degrees of freedom are uncoupled from each other in the probability distribution. An approximate functional for Dl is derived in Ref. [5]. The functional is constructed by adopting the PercusYevick-type approximation in the unfavorable region of the solute–solvent interaction and the hypernetted-chaintype approximation in the favorable region. Dl is then given by a set of equations listed as
qe ðeÞ w ðeÞ ¼ kB T log e q0 ðeÞ e
e
ð6Þ
Z
dðe gÞ e 1 v0 ðe; gÞ qe0 ðeÞ
e e
q ðgÞ q0 ðgÞ
we0 ðeÞ ¼ kB T
Dl ¼ kB T
dg
Z de
e
q ðeÞ
qe0 ðeÞ
ð7Þ
þ bwe ðeÞqe ðeÞ
e
faðeÞF ðeÞ þ ð1 aðeÞÞF0 ðeÞg q ðeÞ
qe0 ðeÞ
ð8Þ
F ðeÞ ¼
8 > < bwe ðeÞ þ 1 þ
e
bw ðeÞ expð bwe ðeÞÞ 1
> : 1 bwe ðeÞ 2
ðwhen we ðeÞV0Þ ðwhen we ðeÞz0Þ
ð9Þ 8 log 1 bwe0 ðeÞ > > < log 1 bwe0 ðeÞ þ 1 þ e bw0 ðeÞ F0 ðeÞ ¼ > > 1 bwe ðeÞ : 2 0
when we0 ðeÞV0
when we0 ðeÞz0
ð10Þ
að e Þ ¼
8 <1 :1
qe ðeÞ qe0 ðeÞ qe ðeÞ þ qe0 ðeÞ
2
when qe ðeÞzqe0 ðeÞ when qe ðeÞVqe0 ðeÞ ð11Þ
3. Procedures In the present work, the solvent is water. The solution is dilute and contains a single solute molecule. The thermodynamic state of each system of interest is specified by the water density and temperature. Four thermodynamic states are examined. One is an ambient state of 1.0 g/cm3 and 25 8C and the others are supercritical states of 1.0, 0.6, and 0.2 g/cm3 and 400 8C. The water molecule is treated as rigid and nonpolarizable, and the SPC/E model is adopted as the intermolecular potential function between water molecules [7].
25
All the nonpolar, polar, and ionic species are examined for rigid solutes. The nonpolar solutes treated are methane and ethane. Water, methanol, and ethanol are employed as polar solutes, and Na+ and Cl as ionic solutes. The interaction potential and the scheme for calculating the solvation free energy are described in Ref. [5]. The flexible solute molecules employed in the present work are comprised of two or four interaction sites. The two- and four-site molecules involve stretching and torsional degree of freedom, respectively, for intramolecular motion. The potential function for the stretching motion of the two-site solute is harmonic and is given by K ðr re Þ2
ð12Þ
where r is the distance between the two sites and r e is taken to be 3.0 2 [6]. Jorgensen’s parameter set for nbutane is adopted to model the four-site solute [8]. The bond lengths and angles are fixed at the values listed in Table 1 of Ref. [8] and the torsional degree of freedom is described by the potential energy expressed as Eq. (2) and Table 2 of Ref. [8]. The intermolecular interaction between the flexible solute and water molecules consists of the Lennard-Jones and Coulombic terms, as usual, and is given by Eq. (23) of Ref. [5] with the truncation factor S(x) set to unity. The Lennard-Jones parameters for the two-site molecule are taken from Ref. [9]. The values for the sodium ion in Table 1 of Ref. [9] are assigned to one of the sites, and those for the chloride ion to the other. For the four-site molecule, Jorgensen’s set for n-butane in Table 3 of Ref. [8] is employed as the Lennard-Jones parameters. The LennardJones parts of the solute-water potential functions are then constructed by the standard Lorentz-Berthelot combining rule [2]. For the Coulombic interaction, two cases are examined. One is the nonpolar case, in which no charge is given to any of the sites. In the other case, the charge of +1 in the unit of elementary charge is carried by the Na+-like site of the two-site molecule and by one of the methyl-like sites of the four-site molecule. The Cl like site of the twosite molecule and the other methyl-like site of the four-site molecule involve the corresponding negative charges, so that the molecules are neutral in total. In the rest of the paper, the two-site molecule with K=2 kcal/mol/22 in Eq. (12) is called HM-2-0 and HM-2-1, respectively, when the charge on the Na+-like site is 0 and +1 in the unit of elementary charge. Similarly, the two-site molecule with K=5 kcal/mol/22 is denoted as HM-5-0 and HM-5-1, respectively, at the charge of 0 and +1. The four-site molecule is called Bu-0 and Bu-1, respectively, when one of the methyl-like sites has a charge of 0 and +1. In any case, no charge is placed on the nonterminal methylene-like sites of the four-site solutes. To obtain the energy distribution function q e, a Monte Carlo simulation was carried out for the solution system of interest. The distribution functions q 0e and v 0e were also
26
N. Matubayasi, M. Nakahara / Journal of Molecular Liquids 119 (2005) 23–29
calculated from a simulation of the pure solvent system. In this case, the solute molecule was inserted as a test particle at a random position with a random orientation. The intramolecular degree of freedom of the solute molecule was sampled simply through the intramolecular potential (or through the modified scheme given in Ref. [6] when the modified form of functional in Ref. [6] was employed for
the solvation free energy). It should be noted that when the solute molecule is inserted into the pure solvent system, it often overlaps with solvent molecules. The overlapping configurations contribute to q 0e and v 0e at large energy coordinate and account for the excluded volume effect in the solvation free energy. The system size and the simulation lengths correspond to those in Refs. [5,6].
Fig. 1. Graphical comparison of the solvation free energy Dl for the rigid solute molecules. The approximate values are provided by Eqs. (6)–(11) and the exact values are obtained from the free energy perturbation calculations in Ref. [5].
N. Matubayasi, M. Nakahara / Journal of Molecular Liquids 119 (2005) 23–29
27
4. Results and discussion In Fig. 1, the approximate values of the solvation free energy Dl for the rigid solute molecules are graphically compared to the corresponding exact values obtained from Table 1 Solvation free energies in the unit of kcal/mola Solute
Thermodynamic Exact state
HM-2-0 1.0 g/cm3 and 25 8C 1.0 g/cm3 and 400 8C 0.6 g/cm3 and 400 8C 0.2 g/cm3 and 400 8C HM-5-0 1.0 g/cm3 and 25 8C 1.0 g/cm3 and 400 8C 0.6 g/cm3 and 400 8C 0.2 g/cm3 and 400 8C HM-2-1 1.0 g/cm3 and 25 8C 1.0 g/cm3 and 400 8C 0.6 g/cm3 and 400 8C 0.2 g/cm3 and 400 8C HM-5-1 1.0 g/cm3 and 25 8C 1.0 g/cm3 and 400 8C 0.6 g/cm3 and 400 8C 0.2 g/cm3 and 400 8C Bu-0 1.0 g/cm3 and 25 8C 1.0 g/cm3 and 400 8C 0.6 g/cm3 and 400 8C 0.2 g/cm3 and 400 8C Bu-1 1.0 g/cm3 and 25 8C 1.0 g/cm3 and 400 8C 0.6 g/cm3 and 400 8C 0.2 g/cm3 and 400 8C
Approximate Eqs. (6)–(11) Modified formb
3.2F0.2
3.5F0.1
3.4F0.1
15.1F0.3
12.4F0.1
12.5
3.5F0.1
3.9F0.1
3.9
0.7F0.1
0.7
0.7
3.3F0.2
3.6F0.1
3.4F0.2
15.3F0.3
12.6F0.1
12.6F0.1
3.5F0.1
3.9
3.9
0.7F0.1
0.7
0.7
104.8F0.9
–c
108.5F1.8
86.2F0.5
–c
87.2F0.5
86.4F0.3
–c
82.9F0.4
78.6F1.3
–c
70.9F0.4
92.5F0.7
–c
97.5F1.7
74.1F0.7
–c
74.3F0.2
74.7F0.5
–c
70.5F0.3
67.4F1.0
–c
60.8F0.7
3.1F0.4
4.4F0.2
4.5F0.3
21.6F0.2
17.9F0.1
17.9F0.1
3.8F0.2
4.8F0.1
4.8F0.1
0.7F0.2
0.6
0.6
63.5F0.7 63.2F1.0
63.3F1.0
41.1F0.7 43.9F0.3
44.7F0.3
47.2F0.4 46.4F0.1
47.2F0.1
39.5F0.4 38.4F0.4
38.7F0.4
a Each value is rounded to a multiple of 0.1 kcal/mol. The error is expressed as 95% confidence level and is not shown when it is rounded to 0. b A modified form of approximation given by Eqs. (29) and (35)–(41) of Ref. [6]. c The value calculated is numerically unstable or involves a large discretization error typically of a few kcal/mol.
Fig. 2. The probability distribution functions p in solution and p 0 at isolation of the dihedral angle / of the Bu-0 solute at an ambient state of 1.0 g/cm3 and 25 8C and at supercritical states of 1.0, 0.6, and 0.2 g/cm3 and 400 8C. The solute is at the trans conformation when /=p, and it is at the gauche when /=p/3 or 5p/3. At convergence, p and p 0 are symmetric against the reflection with respect to /=p. The asymmetry is present due to the statistical error. It should be noted that p 0 is common to the three supercritical states.
the standard free energy perturbation method. The good agreement is observed between the approximate and exact values. The agreement is particularly notable at the low- and medium-density states of 0.2 and 0.6 g/cm3 and 400 8C. When the solute is ionic, the density at the state of 0.2 g/cm3 and 400 8C is not yet low in the sense, for example, that the hydration number at that state is comparable to the numbers at ambient states [10,11]. Even in this case, our approximate procedure is effective in determining the solvation free energy. The solvation free energies of water at 1.0 g/cm3 and 400 8C and of methanol and ethanol at 0.6 g/cm3 and 400 8C are rather small in magnitude. These behaviors are caused by the balance between the favorable and unfavorable contributions of the solute–solvent interactions, and are well reproduced by our approximate method. Therefore, the single functional given by Eqs. (6)–(11) provides an accurate and efficient route to the solvation free energy for various types of solutes over a wide range of thermodynamic conditions. The Dl values for the flexible solute molecules are presented in Table 1. As noted in Ref. [6], the approximate expression listed as Eqs. (6)–(11) is not successful for the HM-2-1 and HM-5-1 solutes. Their intramolecular structures are largely changed through the solute–solvent interaction, and the modified scheme of approximation described in Ref. [6] is needed for good performance. The approximate Dl calculated by Eqs. (6)–(11) are in good agreement with the exact for the HM-2-0, HM-5-0, Bu-0,
28
N. Matubayasi, M. Nakahara / Journal of Molecular Liquids 119 (2005) 23–29
Fig. 3. The probability distribution functions p in solution and p 0 at isolation of the dihedral angle / of the Bu-1 solute at an ambient state of 1.0 g/cm3 and 25 8C and at supercritical states of 1.0, 0.6, and 0.2 g/cm3 and 400 8C. The solute is at the trans conformation when /=p, and it is at the gauche when /=p/3 or 5p/3. At convergence, p and p 0 are symmetric against the reflection with respect to /=p. The asymmetry is present due to the statistical error. It should be noted that when the Ewald method is employed, p 0 at a fixed temperature depends slightly on the density due to the interaction of the solute with its own images.
and Bu-1 solutes. For the HM-2-0 and HM-5-0 solutes, it is shown that the solvent effect on the solute intramolecular structure is weak. In this case, the modified form of approximation provided in Ref. [6] does not lead to an improvement and the original form given by Eqs. (6)–(11) performs well. In Table 1, the Dl values calculated by the modified form of approximation in Ref. [6] are also listed for the Bu-0 and Bu-1 solutes. They are observed to be coincident with the corresponding values from Eqs. (6)–(11). To see the origin of this observation, we examine the probability distribution function of the torsional angle / of the four-site molecule. In the following, the probability distribution function of the intramolecular coordinate / is denoted by p in solution and by p 0 at isolation (absence of the solute– solvent interaction). In Fig. 2, we show p and p 0 for the Bu-0 solute. Actually, the Bu-0 solute models n-butane and is often employed as a prototypic solute to study the hydrophobic effect [12,13]. As expected from the common notion of hydrophobicity, it is seen at each thermodynamic state that the trans conformation becomes less populated through the solute–solvent interaction. In other words, a compact conformation is preferable in water over a wide range of thermodynamic conditions. When the temperature is elevated from 25 to 400 8C at a fixed density of 1.0 g/cm 3 , p/p 0 at the typical trans and gauche conformations is observed to vary within 20%. Actually,
this observation is consistent with a previous finding that the solvent effect on the association and dissociation of a pair of nonpolar solutes is weakly dependent on the temperature when the density of solvent water is fixed [14]. The density reduction at a fixed supercritical temperature of 400 8C weakens the effect of solvation. At a low density of 0.2 g/cm3, the solvent effect is barely appreciable. In Fig. 3, p and p 0 are shown for the polar four-site solute, Bu-1. The solute molecule is extended due to the solute–solvent interaction, and the trans conformation is exclusively occupied in the solution systems examined. This observation is parallel to that for HM-2-1 and HM-5-1 reported in Ref. [6]. While p is broader at the higher temperature, the solvent effect is still bstrongQ even at a lowdensity state of 0.2 g/cm3 and 400 8C in the sense that the isothermal density reduction at 400 8C from 1.0 to 0.2 g/cm3 causes a weaker effect on p than that from 0.2 g/cm3 to the zero density limit. It is seen in Figs. 2 and 3 that the main domain for p overlaps significantly with the peak domain for p 0. In other words, the typical conformations of the solute in solution are frequently realized for the solute at isolation. In this case, the modification provided in Ref. [6] does not lead to an improvement of the approximate scheme, in agreement with the arguments presented in Ref. [6].
Acknowledgments This work is supported by the Grant-in-Aid for Scientific Research (Nos. 11740322, 13640509, and 15205004) from the Japan Society for the Promotion of Science and by the Grant-in-Aid for Scientific Research on Priority Areas (No. 15076205) and the NAREGI (National Research Grid Initiative) Project from the Ministry of Education, Culture, Sports, Science, and Technology. N.M. is also grateful to the Supercomputer Laboratory of Institute for Chemical Research, Kyoto University, for generous allocation of computation time.
References [1] N. Matubayasi, L.H. Reed, R.M. Levy, J. Phys. Chem. 98 (1994) 10640. [2] J.P. Hansen, I.R. McDonald, Theory of Simple Liquids, Academic Press, London, 1986. [3] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Oxford University Press, Oxford, 1987. [4] N. Matubayasi, M. Nakahara, J. Chem. Phys. 113 (2000) 6070. [5] N. Matubayasi, M. Nakahara, J. Chem. Phys. 117 (2002) 3605, 118 (2003) 2446. [6] N. Matubayasi, M. Nakahara, J. Chem. Phys. 119 (2003) 9686. [7] H.J.C. Berendsen, J.R. Grigera, T.P. Straatsma, J. Phys. Chem. 91 (1987) 6269. [8] W.L. Jorgensen, J.D. Madura, C.J. Swenson, J. Am. Chem. Soc. 106 (1984) 6638.
N. Matubayasi, M. Nakahara / Journal of Molecular Liquids 119 (2005) 23–29 [9] J. Chandrasekhar, D.C. Spellmeyer, W.L. Jorgensen, J. Am. Chem. Soc. 106 (1984) 903. [10] P.B. Balbuena, K.P. Johnston, P.J. Rossky, J. Phys. Chem. 99 (1995) 1554. [11] P.B. Balbuena, K.P. Johnston, P.J. Rossky, J. Phys. Chem. 100 (1996) 2706.
29
[12] R.O. Rosenberg, R. Mikkilineni, B.J. Berne, J. Am. Chem. Soc. 104 (1982) 7647. [13] W.L. Jorgensen, J. Gao, C. Ravimohan, J. Phys. Chem. 89 (1985) 3470. [14] N. Matubayasi, M. Nakahara, J. Phys. Chem., B 104 (2000) 10352.