Hydrophobic Solvation: Aqueous Methane Solutions

May 5, 2007 - Dartmouth College. Hanover ... (2, 3) published recently in this Journal are more closely re- ... second provides a brief review of Henr...
1 downloads 0 Views 185KB Size
Research: Science and Education edited by

Advanced Chemistry Classroom and Laboratory

Joseph J. BelBruno Dartmouth College Hanover, NH 03755

Hydrophobic Solvation: Aqueous Methane Solutions Oliver Konrad Universität Hamburg, Institut für Physikalische Chemie, Martin-Luther-King Platz 6, D-20146 Hamburg, Germany Timm Lankau* Department of Chemistry, Room 303, National Tsing Hua University, 101 KuangFu Road Sec. 2, Hsinchu 300, Taiwan; *[email protected]

An interesting application of Henry’s law to noisy knuckles was provided by Kimborough (1) in 1999. Two articles (2, 3) published recently in this Journal are more closely related to the ideas presented here. The first describes experiments on the hydrophobic solvation of toluene (2) and the second provides a brief review of Henry’s law (3). This article focuses on the solvation of hydrophobic gases and on the application of Henry’s law. Most undergraduate students are familiar with the concept of solvation shells around ionic solutes and with basic electrochemical methods to detect and analyze them. The idea of a solvation shell around an apolar solute as well as its detection seems to baffle them, as they expect electrostatic interactions among solute and solvent to be the source of such a structural effect. This article is based on a short lecture, that gives a basic introduction to this subject without being too specific about the applied methods. The first section will be used to develop an expression for Henry’s law (eq 14) as it is frequently used today (4). Since this explanation is based solely on phenomenological thermodynamics, it is possible to discuss the thermodynamics of hydrophobic solvation in detail even in undergraduate courses. An application of Henry’s law to methane solvation is provided in the second section. Results from computational simulations will be used to discuss structural aspects of hydrophobic solvation and how these structural effects relate to the thermodynamic properties of hydrophobic solvation. Henry’s Law The dissolution process of methane in water can be represented as a chemical reaction,

CH4(gas)

CH4(sol)

(1)

where CH4(gas) describes “methane in the gas phase” as the educt and CH4(sol) describes “methane in aqueous solution” as the product of the reaction. Owing to the low concentration of methane in the solution, we can assume that the vast majority of water molecules are not affected by the methane and that it is possible to neglect them in the following discussion. The chemical drive, ∆µ, of this reaction can therefore be computed by the difference of the chemical potentials of methane in solution and in the gas phase: ∆µ = µ CH 4 (sol) − µ CH4 (gas) 864

Journal of Chemical Education



(2)

To simplify the argument, we assume further an ideal behavior of the solution and the gas, n ∆µ = µ° CH4 (sol) + RT ln

p c − µ° CH4 (gas) − RT ln c ° p ° (3)

where µC H4(sol) and µC H4(gas) are the chemical potentials of methane in the solvent and in the gas phase under standard conditions, c is the methane concentration in solution, p is the methane pressure, R is the universal gas constant, T is the absolute temperature, c represents the standard concentration of 1 mol兾L, and p the standard pressure of 1 bar. Chemical equilibrium is established when the drive of the reaction vanishes, ∆µ = 0. Hence, it is possible to rewrite eq 3 as p p ° ∆ µ ° c = exp = KH (4) RT c c ° where ∆µ = µC H4(sol) − µC H4(gas) is the difference in the standard potentials and KHc is the Henry constant based on the solute’s concentration. Equation 4 suggests that the value of KHc and therefore the solubility of methane in water depends solely on the difference of the standard chemical potentials ∆µ and the absolute temperature T. Hydration effects are difficult to observe directly with eq 4, as the number of molecules in solution and in the gas phase differ significantly. One liter of a standard methane solution contains 1 mol of methane and 1 L of methane gas under standard conditions contains only (1兾22.4) mol = 0.04 mol. It is therefore helpful to remove differences in the chemical potential arising from standard definitions to observe directly the consequences of methane solvation. This can be done by defining a new reference state. To distinguish between the new reference and the commonly used standard state, the new reference state is labeled with an open circle  while the commonly used standard definition is marked by  . The definition of the new reference state for the calculation of the Henry constant KH starts from a sample of volume V  taken from the saturated methane solution. This sample defines the standard for the chemical potential of methane in aqueous solution µCH4(sol), which is fixed by the temperature of the sample T , its volume V , and the number n2 of methane molecules in the sample. As we will see later, the inclusion of the temperature T into the new standard (T ) is not necessary for the argument. Most students would define a standard by temperature, pressure, and concentration. To stay close to the student’s perception

Vol. 84 No. 5 May 2007



www.JCE.DivCHED.org

Research: Science and Education

of a standard and to maintain the idea to define a standard by the experimental environment, T was included into the discussion. Finally, the new reference concentration can be computed from the properties of the sample,

n2° (5) V ° The second reference state (methane gas, µCH4(gas) ) is created from the first by removing all water molecules n1 from the sample while keeping T  and V  constant. If we assume further ideal behavior for the dilute methane gas, the new reference pressure p can be calculated from the ideal gas law,

Although chemists generally express concentrations by the number of particles per unit volume, it is easier to use mole fractions such that X2 is the mole fraction of methane,

p° V ° = n2° RT °

(6a)

and therefore n° RT ° p ° = 2 V °

(6b)

Repeating the initial calculations with the new reference states replaces eq 4 by a new set of reference variables: p p° ∆µ° = exp c c ° R T °

(7)

Additional physical insight can be drawn from the difference of the new reference potentials, ∆µ = µCH4(sol) − µCH4(gas). As the reference states differ only in the number of water molecules n1 it is possible to identify ∆µ as the excess chemical potential of hydration µex = µex(T  ), which summarizes the influence of the methane–water interactions on the chemical potential of methane. The second part of the excess free enthalpy of solvation, the excess chemical potential of water, is not covered by eq 1 and can be neglected in the discussion owing to the low methane concentration. (The ratio between water molecules with and without a direct contact to a methane molecule is approximately 1:1111, if we assume a methane concentration of 2.5 mmol兾L and that each methane molecule is surrounded by 20 water molecules.)

p p° µ° ex = exp p c c° RT °

(8)

Since the volume of the sample V  is the same in both reference states (gas phase and aqueous solution), it is possible to use eqs 5 and 6b to simplify the pre-exponential factor p° n ° R T ° = 2 = R T ° c ° V ° c°

(9a)

such that p µ° ex = RT ° exp c RT °

(9b)

This calculation is valid for any given temperature T  and therefore it is possible to use the symbol T instead of T  in the following, where µex becomes a function of the temperature µex = µex(T ): p µ = RT exp ex c RT www.JCE.DivCHED.org

(10)



n2 n ≈ 2 n1 + n2 n1

X2 =

c ° =

(11a)

since (11b)

n1 >> n2

where n1 is the number of water molecules in moles in solution and n2 is the same for methane. Since the solubility of methane in water is low, X2 is small. Hence, the density of the solution ρsol may be replaced by the density of the solvent ρH2O, Vsol =

mH2 O + m CH 4 n1 MH2 O m sol = ≈ ρsol ρsol ρH2 O

(12)

where mH2O is the mass of water in the sample, mCH4 is the mass of methane in the sample, and MH2O is the molar mass of water. Since mH2O is much larger than that of methane mCH4, the latter term can be neglected. Also note that mH2O = n1MH2O. Applying eqs 11a and 12 it is possible to write p = c

p Vsol n2

=

p n1 M H2 O n2 ρH2 O

=

p MH2 O X2 ρH2 O

(13)

and consequently eq 10 becomes R T ρH2 O p µ = exp ex RT X2 M H2 O

X

= KH

(14)

where KHX is the Henry constant basing on the solute’s mol fraction. The ratio ρH2O兾MH2O is equal to the number of moles solvent per volume unit, thus ρ H 2O 兾M H 2O = mH2O兾(VH2OMH2O) = (nH2OMH2O)(VH2OMH2O) = nH2O—兾VH2O, which is equal to the inverse of the molar volume VH2O of the solvent and can be used to simplify eq 14 further: X KH =

p RT µ = exp ex X2 VH2O RT

(15)

This expression of KHX has the advantage that it shows clearly that the unit of KHX is equal to that of a pressure on both sides of eq 14. — Since most students are more familiar with ρH2O than with VH2O and tables of solvent densities as a function of temperature are more easily found than those of molar volumes, we will retain eq 14 in the following discussion. It is tempting to argue that eq 14 is more difficult to use for teaching purposes than eq 4, because eq 14 is more complicated and it takes much longer to derive eq 14 from the chemical equilibrium conditions than eq 4. On the other hand, two important issues can be raised in favor of eq 14. 1. µex describes purely solvation effects, whereas the discussion of ∆µ is complicated by effects originating from the definition of the thermodynamic standard states. 2. Equation 14 is a function of the density of the solvent ρH2O. This can be used to move the focus from thermodynamical to structural aspects of hydrophobic solvation.

Vol. 84 No. 5 May 2007



Journal of Chemical Education

865

Research: Science and Education

Figure 1. Simulated (6) and experimental values (7–9) for KHX as a function of temperature T. Experimental values (crosses) are connected by a spline interpolation and values from the simulation (diamonds) by a polynomial fit. A better agreement between both curves can be obtained, if selected experimental measurements were excluded from the smoothing procedure. The plot shows the fit if all points are used.

Figure 2. The excess chemical potential of solvation from simulation (solid line and diamonds) and experiment (cross).

Methane Solvation Methane is a standard example for the hydrophobic solvation of small molecules. A comparison of our molecular dynamics (MD) (5) simulation results with experimental data is shown in Figure 1. A brief summary of the computational setup is provided in the Supplemental MaterialW while details of the simulations are published in ref 6. Experimental data for KHX have been taken from the literature (7–9). A typical feature of hydrophobic gas solvation is the maximum in KHX. The analysis of the simulations in ref 6 showed that the value of KHX at the maximum is controlled by the methane– water interaction potential while its position on the temperature scale is dominated by the water–water interaction potential used for the simulation. These observations suggest that the amount of methane in the solution is controlled by direct solute–solvent interactions, while the temperature of minimum methane solubility is controlled by the inner structure of the solvent (6, 10). The shape of both curves in Figure 1 are in good agreement and differ only in the position of the maximum. Hence, either curve can be used for the discussion of hydrophobic solvation. We will use the results from the computer simulations, since the simulations give rise to all thermodynamic data simultaneously plus additional information on the microscopic structure of the methane solution. Figure 2 shows the value of µex obtained from computer simulations (diamonds) and experiments in the temperature range of liquid water (cross). The solid line is a polynomial fitted to results of the simulations,

µ ex = a + bT + cT

2

+ dT

3

+ eT

4

(16)

which will be used later for the calculation of derivatives.1 The agreement between simulated and experimental values for µex is much better than that observed for KHX. The misalignment between theory and experiment in Figure 1 is caused by the poor agreement between theory and experiment for the density of the solvent ρH2O (Figure 3). Such a dependence of KHX on the solvent’s density is another clue that the inner structure of the solvent is an important factor in the discussion of hydrophobic solvation. 866

Journal of Chemical Education



Figure 3. Density of water obtained from computer simulations.

The maximum in µex suggests that the usual linear approach for the temperature dependence of µex T

µ ex (T ) = µ° ex −

∫ Sex dT

° (T − T ° ) ≈ µ° ex − S ex (17)

T °

(Sex is the excess entropy of solvation) cannot be used to describe hydrophobic solvation accurately and that the temperature dependency of Sex has to be taken into account. Figure 4 depicts Sex as a function of temperature. Since all simulations were done at the same pressure (6), we were able to use the expression Sex = −

∂ µex ∂ T

p

(

)

= − b + 2 cT + 3 d T 2 + 4 eT 3 (18)

to calculate Sex directly from µex via eq 16. In the temperature range of liquid water, Sex is negative. The solvation of a hydrophobic gas reduces the entropy of the system. Frank and Evans proposed the “iceberg model” (11) to explain this reduction. Using the interpretation of entropy as a measure for order and disorder, Frank and Evans suggested that the hydrophobic solute stiffens the solvent in its direct neighborhood, which increases the crystallinity of the solvent and therefore its order. They did not propose a structure of the “iceberg” and they stressed that the term iceberg was chosen for convenience and was not intended to suggest a structure similar to that of ice.

Vol. 84 No. 5 May 2007



www.JCE.DivCHED.org

Research: Science and Education

Figure 4. Excess entropy of solvation Sex as a function of temperature T calculated from the simulation.

Figure 5. Excess enthalpy of solvation Hex as a function of temperature T calculated from the simulation.

With increasing values of T the excess entropy Sex becomes more positive as shown in Figure 4. This increase can be rationalized by the simultaneous increase of molecular motions within the iceberg, which correlates with a lower degree of order and therefore with a smaller entropy decrease. The chemical work associated with the formation of the iceberg can be computed from µex by subtracting the entropy term TSex from µex. The chosen simulation conditions allow us to identify this chemical work as the excess enthalpy Hex of the process, H ex = µex + T Sex = a − cT 2 − 2 dT 3 − 3 eT 4 (19) Hex varies between 20 kJ兾mol at low temperatures and +20 kJ兾mol at the other end of the temperature scale with an average value close to zero. This is much smaller than expected for a large rearrangement of the solvent structure. The chosen model for the water–water interactions predicts the strength of the hydrogen bond in the water dimer (H2O)2 to be close to 7.17 kcal兾mol ≈ 30 kJ兾mol with a bond length of 2.75 Å, whereas experiments suggest a value of approximately 23 kJ兾mol for the energy and a bond length of 3 Å (12–14). It is therefore possible to say that Hex represents 0 ± 1 hydrogen bond! Today, the iceberg effect is commonly explained by the formation of an inverse solvation shell around the hydrophobic solute. A detailed analysis of the simulation at 300 K shows that 20 water molecules surround the methane molecule at an average carbon–oxygen distance of 3.81 Å. This www.JCE.DivCHED.org



Figure 6. Idealized model of the (H2O)20 cage with a methane molecule inside.

information can be used to compute the methane–water interaction energy with the given methane–water interaction potential as 18 kJ兾mol, which is again in the energetic range of an isolated hydrogen bond. Figure 5 shows that Hex (300 K) is approximately 6 kJ兾mol, which suggests that the energy to form a cage in water for the methane molecule is about +12 kJ兾mol—approximately half the binding energy of an isolated hydrogen bond. The small quantity of energy necessary for the formation of the cavity can be rationalized through its structure. Figure 6 depicts an idealized image of a methane molecule and its first solvation shell. MD simulations suggest that about 20 water molecules form the solvation shell around the methane molecule (6, 15). A first approximation to its structure is the water dodecahedron observed in gas hydrates (16, 17). Experimental evidence for and against this simple model has been published previously. Bowron et al. (18) reported EXAFS (extended X-ray absorption fine structure) experiments on the formation of krypton hydrates that suggest that the structure of the first solvation shell and that of the cage are very similar, whereas Koh et al. (19) suggest that the first solvation shell around a methane molecule is considerably smaller than the water dodecahedron in the gas hydrate. Snapshots from MD simulations show that the solvation cage around the methane is less symmetrical than the dodecahedron in Figure 6 and its shape is closer to that of a potato tuber or a battered golf ball than to that of a sphere. Further, the number of molecules in the first solvation shell can change by one or two water molecules as time passes. Despite these problems, the dodecahedron can be used to highlight a principal feature of hydrophobic solvation. First, we need an estimate of the number of hydrogen bonds per water molecule in liquid water. A water molecule is tetrahedrally coordinated in regular ice (ice Ih), which floats on liquid water. The water molecules in the liquid must therefore be packed more densely than those in ice, which again suggests that the number of hydrogen bonds per water molecule is slightly larger than four. Such a small increase in the coordination number can be rationalized by thermal fluctuations that can move water molecules closer to each other than in ice and thus enables penta-coordinated water molecules to exist for a very short time. The water molecules in a hydrophobic solvation shell do not build direct bonds with the solute but form a hydrogen bonded network. The molecules are oriented in such a

Vol. 84 No. 5 May 2007



Journal of Chemical Education

867

Research: Science and Education

way that their molecular plane points to the solute and their “sticky ends” are either tangentially oriented or point away from the methane molecule. The number of hydrogen bonds per water molecule within the network is three, as Figure 6 suggests. Small deviations from this number are possible owing to thermal fluctuations. Each water molecule can form another hydrogen bond with water molecules surrounding the shell either as hydrogen donor (10 oxygen–hydrogen bonds point away from the methane molecule) or as hydrogen acceptor (10 lone electron pairs from the water molecules are oriented similarly than the dangling hydrogen atoms). The number of hydrogen bonds per water molecule in the first solvation shell is therefore close to four, but likely to be slightly smaller. It is a unique feature of water as solvent that the number of hydrogen bonds per water molecule hardly changes during the formation of the cavity for the hydrophobic solute. It has been suggested that the hydrogen bonds that stabilize the cavity are stronger than those in bulk water (20). Such an increase in hydrogen bond stability can be rationalized by the decrease of translational and rotational mobility of the cage building water molecules. As the temperature of the system rises, such a stabilization of the cage becomes less likely and the formation of the cavity is energetically more expensive, which can be observed in the increasing value of Hex. The fragility of this argument can be seen in Figure 3. With increasing values of temperature T the density of the solvent ρsol decreases. The formation of a cavity can be regarded as an artificial increase in the distance among a subset of water molecules. A decrease in the density of the solvent should therefore result in a reduction of the energetic costs of cavity formation and Hex should become more negative, which is not seen in Figure 5. Figure 7 can be used to analyze these effects in greater detail. For this analysis eq 14 for KHX was split into two parts: the pre-exponential factor F(T ) and the exponential function E(T ): K HX

= F (T ) E (T )

F (T ) =

R T ρH2 O

(20)

M H2 O

E (T ) = exp

µ ex RT

Journal of Chemical Education

If we assume the excess chemical potential µex to be independent of the temperature, the value of E(T ) should decrease monotonically with increasing T. Figure 7 shows an unexpected maximum in E(T ) at approximately 325 K. For T ≤ 325 K the increase in µex (Sex < 0, Figures 2 and 4) compensates for the decrease caused by the increase of T in the denominator. After passing through a plateau between 325 and 350 K, E(T ) decreases for T > 350 K. The value of µex does not change significantly (Sex ≈ 0) in this region until it finally decreases (Sex > 0, Figures 2 and 4) and the value of E (T ) is therefore dominated by the T 1 term in the argument of the exponential function. The shape of E(T ) is therefore the result of the nonlinearity of µex (The changes in Sex are too large for Sex to be assumed constant.) and a simple approach such as eq 17 cannot be used satisfactorily for the analysis of the hydrophobic solvation of small molecules. The shape of the function E(T ) is similar to that of KHX(T ) as shown in Figure 1. Multiplication by F(T ) simply increases the value of KHX and moves the maximum to larger values of T, which suggests that the changes in the solvent’s density ρH2O have a stronger influence on the work necessary to form a cavity in the solvent than on volume effects in the chosen measurement2 for the gas concentration in the solvent. Conclusions

If we assume ρH2O to be independent of temperature, the preexponential factor F should be a linear function of temperature T (indicated by a dotted line in Figure 7). Figure 3 displays the decrease of ρH2O with the temperature. The slope of F(T ) should therefore decrease as the temperature increases until it finally becomes negative (Figure 7). But, F(T ) increases with temperature in the chosen temperature range, because the decrease in ρH2O is smaller than the increase in T. The comparison between the solid and the dotted line for F(T ) shows that the decrease of ρH2O causes a decrease in KHX, which is equivalent to an increase in the methane mole fraction in the solution. Low density solvents should therefore dissolve methane more easily than high density ones.

868

Figure 7. Analysis of eq 14 as shown in eq 20. F(T) is the preexponential factor and E(T) the exponential function. The dotted line is a linear extrapolation of F(T) with a constant value for ρH2O constructed as a tangent to F(T) for small values of T.



The standard form of Henry’s law (eq 4) has been modified (eq 14) to show the connections among the Henry constant KHX, the density of the solvent ρH2O, and the excess chemical potential of solvation µex. The methane solubility in water as a function of temperature shows a typical minimum at approximately 350 K. The corresponding maximum in KHX(T ) correlates with a maximum of µex, which again is the result of the change of sign in Sex. These results, as well as the zero value of Hex, can be rationalized by the internal structure of water, which is governed by hydrogen bonding. The hydrophobic solvation process is therefore a good teaching example, which connects macroscopic, phenomenological thermodynamic results with an atomistic point of view.

Vol. 84 No. 5 May 2007



www.JCE.DivCHED.org

Research: Science and Education

This connection becomes even more important, as the standard linear approach (eq 17) can be used only with difficulty to explain the maximum in KHX. Further, the discussion can be easily extended towards computational simulation techniques such as molecular dynamics or Monte Carlo methods. Acknowledgments The authors would like to thank K. Nagorny (general support of this work), I. L. Cooper (discussion of eq 14), and the Job-Foundation (financial support). Notes 1. The values of the constants in the polynomial for the description of µex(T ) between 240 and 460 K are as follows: a = 117.173, b = 1.7742, c = 0.00421301, d = 7.01704 × 106, and e = -4.5751 × 109. 2. One should keep in mind that ρsol in eq 14 originates from the volume Vsol of the sample taken from the methane solution.

Literature Cited 1. Kimborough, Doris R. J. Chem. Educ. 1999, 76, 1509–1510. 2. Serafin, J. M. J. Chem. Educ. 2003, 80, 1194–1196. 3. Rosenberg, R. M.; Peticolas, W. L. J. Chem. Educ. 2004, 81, 1647–1652. 4. Boulougouris, G. C.; Voutsas, E. C.; Economou, I. G.; Theodorou, D. N.; Tassios, D. P. J. Phys. Chem. B 2001, 105, 7792–7798. 5. Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids, 1st ed.; Oxford University Press: New York, 1987. 6. Konrad, O.; Lankau, T. J. Phys. Chem. B. 2005, 109, 23596– 23604. 7. Rettich, T. R.; Handa, Y. P.; Battino, R.; Wilhelm, E. J. Phys. Chem. 1981, 85, 3230–3237. 8. Cramer, S. D. Ind. Eng. Chem. Process Des. Dev. 1984, 23, 533–538. 9. Crovetto, R.; Fernández-Prini, R.; Japas, M. L. J. Chem. Phys. 1982, 76, 1077–1086. 10. Paschek, D. J. Chem. Phys. 2004, 120, 6674–6690. 11. Frank, H. S.; Evans, M. W. J. Chem. Phys. 1945, 13, 507– 532. 12. Kim, K. S.; Mhin, B. J.; Choi, U.-S.; Lee, K. J. Chem. Phys. 1992, 97, 6649–6662. 13. Curtiss, L. A.; Frurip, D. J.; Blander, M. J. Chem. Phys. 1979, 71, 2703–2711. 14. Dyke, T. R.; Muenter, J. S. J. Chem. Phys. 1974, 60, 2929– 2930. 15. Konrad, O. Molekulardynamische Simulationen zur Solvatation von Methan in Wasser. Ph.D. Thesis, University of Hamburg, Hamburg, Germany, 2006. 16. Sloan, E. D., Jr. Clathrate Hydrates of Natural Gase, 2nd ed.; Marcel Dekker Inc.: New York, 1998.

www.JCE.DivCHED.org



17. Gas Hydrate Research at Geomar. http://www.gashydrate.de/ (accessed Feb 2007) 18. Bowron, D. T.; Filipponi, A.; Roberts, M. A.; Finney, J. L. Phys. Rev. Lett. 1998, 81, 4164–4167. 19. Koh, C. A.; Wisbey, R. P.; Wu, X.; Westacott, R. E.; Soper, A. K. J. Chem. Phys. 2000, 113, 6390–6397. 20. Silverstein, K. A. T.; Haymet, A. D. J.; Dill, K. A. J. Am. Chem. Soc. 2000, 122, 8037–8041. 21. Lindahl, E.; Hess, B.; van der Spoel, D. J. Mol. Model 2001, 7, 306–317. 22. Nosè, S.; Klein, M. L. Mol. Phys. 1983, 50, 1055–1076. 23. Parrinello, M.; Rahman, A. J. Appl. Phys. 1981, 52, 7182– 7190. 24. Nosè, S. Mol. Phys. 1984, 52, 255–268. 25. Hoover, W. G. Phys. Rev. A 1985, 31, 1695–1697. 26. Miyamoto, S.; Kollman, P. A. J. Comput. Chem. 1992, 13, 952–962.

Appendix All molecular dynamics calculations presented in this communication were done using the GROMACS software package (21). Molecular dynamics simulations for the calculation of the Henry constant KH started with a box containing 216 water and one methane molecule at a pressure of 1 bar at a fixed temperature. Thermodynamic equilibrium was assumed to be reached if no further drift was observed for 400 ps in the mean value of the thermodynamic variables density, pressure, temperature, and potential energy. Each of these equilibrated systems was subsequently used to run 21 different simulations for the thermodynamic integration. Soft core potentials for Lennard-Jones and electrostatic interactions have been used (α = 1.51 kJ兾mol, σ = 0.3 nm) within the thermodynamic simulation. The results from these 21 simulations were numerically integrated with the trapezoid rule to obtain finally the value of the excess chemical potential µex. All simulations were 500 ps long. The leapfrog algorithm (5) (time step 1 fs) in a cubic box with periodic boundary conditions combined with the minimum image convention was used. The cutoff radius was set to 0.9 nm for the LennardJones and electrostatic forces. The long range electrostatic forces were calculated with the reaction field (RF) (5) (εr = 78.5) approach. Additionally the dispersion part of the Lennard-Jones interaction beyond the cutoff was corrected for pressure and energy. The neighbor list with a radius of 0.9 nm was updated every 10 time steps. The NpT-ensembles (1 bar) have been obtained with Parrinello–Rahman pressure (22, 23) (1 ps) and Nosé–Hoover temperature (24, 25) (0.1 ps) coupling. The internal degrees of freedom of the water molecules were constrained with the SETTLE algorithm (26). The densities of pure liquid water (Figure 3) at various temperatures were obtained from the last 300 ps of simulations containing 216 water molecules.

Vol. 84 No. 5 May 2007



Journal of Chemical Education

869