Solubility of Methane in Water: The Significance of the Methane

Solubility of Methane in Water: The Significance of the Methane−Water ... Department of Chemistry, National Tsing Hua University, 101 KuangFu Road S...
0 downloads 0 Views 135KB Size
Download PDF
23596

J. Phys. Chem. B 2005, 109, 23596-23604

Solubility of Methane in Water: The Significance of the Methane-Water Interaction Potential Oliver Konrad UniVersita¨t Hamburg, Institut fu¨r Physikalische Chemie, Martin-Luther-King Platz 6, D-20146, Germany

Timm Lankau* Department of Chemistry, National Tsing Hua UniVersity, 101 KuangFu Road Sec. 2, HsinChu 30013, Taiwan ReceiVed: August 4, 2004; In Final Form: September 26, 2005

The influence of the methane-water interaction potential on the value of the Henry constant obtained from molecular dynamics simulations was investigated. The SPC, SPC/E, MSPC/E, and TIP3P potentials were used to describe water and the OPLS-UA and TraPPE potentials for methane. Nonbonding interactions between unlike atoms were calculated both with one of four mixing rules and with our new methane-water interaction potential. The Henry constants obtained from simulations using any of the mixing rules differed significantly from the experimental ones. Good agreement between simulation and experiment was achieved with the new potential over the whole temperature range.

1. Introduction The calculation of chemical potentials is an active research topic in computational chemistry, as various properties of the simulated system can be derived from the chemical potential. Molecular dynamics (MD) simulations provide an elegant way to calculate both the chemical potential and the corresponding properties.1,2 MD simulations solve Newton’s equations of motion to analyze the development in time of the system. The necessary information about the system’s potential energy is provided through a force field. The choice of the force field is therefore essential for the reliability of the results. A force field is described by a set of potential energy functions and their associated parameters. In the majority of cases, the interactions can be subdivided into four groups (electrostatic, bonded, nonbonded, special) with different potential energy functions. Notable problems arise from the parametrization of the nonbonding interactions between chemically different species, as a large number of parameters is required to address each possible nonbonding interaction. Therefore, mathematical formulas, termed combination or mixing rules, are used to derive these parameters from the force fields of pure substances. Usually, the Lorentz Berthelot rules3 are used for this task, but often, they seem to be an inappropriate choice.3-9 Consequently, a large selection of different mixing rules has appeared in the literature,3-7,10-17 and it seems therefore necessary to verify the applicability of these rules. The solubility of methane in water was chosen because of its relative simplicity. Although numerous MD simulations on methane solvation have been carried out,2,18-28,52,56,57,65 only few focus on solvation energies.2,22-28,56,57,65 Two of these publications28,65 address the calculation of the Henry constant directly from MD simulations as we do here. These two publications differ from our work by the chosen algorithm for the calculation * T.L.: email [email protected] (preferred to phone and fax), phone (+)886-3-5715131-3414, fax (+)886-3-5711082. O.K.: email [email protected].

of the excess chemical potential. While Paschek28 and Guillot and Guissani65 rely on Widom's particle insertion algorithm,68 the thermodynamic integration29 algorithm will be used in this work. In the present paper, the influence of four different mixing rules on the excess chemical potential of methane solvation and the value of the Henry constant was investigated over the temperature range from 250 to 450 K. The value of the Henry constant has been calculated using the thermodynamic integration (TI)29 method with soft core potentials and is compared with results from Monte Carlo (MC) calculations,30,66 MD calculations,28,65 and experimental data.32-34 2. Computational Setup Following a short introduction to the Henry constant and its computation, we will focus briefly on computational details and end this section with a discussion of the various methane-water interaction potentials used in this work. 2.1. Henry Constant. The solubility of a solute (methane) in an ideal dilute solution (water) can be calculated from Henry’s law

KH ) lim

xsolutef0

( ) fsolute xsolute

(1)

where KH is the Henry constant, xsolute is the mole fraction of solute, and fsolute is the fugacity of the solute. It is possible to calculate the absolute value of KH, provided that the excess chemical potential at infinite dilution (µ) and the molar density of the liquid are known.30,31

KH ) lim

xsolutef0

exp (RTF M )

(µ/RT)

(2)

F is the density of the pure solvent, and M is its molar mass. The excess chemical potential µ of a species in a mixture is defined as the difference between the chemical potential of the solute in the solvent and that of the solute as an ideal gas under

10.1021/jp0464977 CCC: $30.25 © 2005 American Chemical Society Published on Web 11/17/2005

Solubility of Methane in Water

J. Phys. Chem. B, Vol. 109, No. 49, 2005 23597

TABLE 1: Water Potentialsa SPC43

SPC/E43

TABLE 2: Methane-Water Interaction Potentialsa MSPC/E44

TIP3P45

r(OH) (nm) 0.1 0.1 0.09839 0.09572 ∠(HOH) (deg) 109.47 109.47 109.47 104.52 σ (nm) 0.316557 0.316557 0.3116 0.315061  (kJ/mol) 0.650629 0.650629 0.61942 0.636812 qO (e) -0.82 -0.8476 -0.8216 -0.834 qH (e) 0.41 0.4238 0.4108 0.417 rdf (nm) 0.277 0.275 0.271 0.278 rcav (nm) 0.363 0.360 0.355 0.364

eq 1 σ  rmin σ  rmin

eq 2

identical conditions. For brevity, µ will be described as the excess chemical potential in the following, always implying infinite dilution. µ equals the difference in the Gibbs free energy between a system containing the solute in an arbitrary amount of solvent (state A) and solute and solvent completely separated (state B).29 The Gibbs free energy difference between the states A and B was calculated using the thermodynamic integration (TI) algorithm1,2,29 as implemented in GROMACS version 3.1.4.35 It was shown in ref 36 that the TI algorithm can be used successfully for the analysis of hydrophobic solvation, because the reaction pathway λ shows no energetic fine structure, which may obscure the result of the thermodynamic integration. 2.2. Computational Details. The calculation of KH started with the MD simulation of a system containing 216 (820) water and 1 methane molecule until thermodynamic equilibrium was reached at a pressure of 1 bar at the fixed temperature. (Note: The parameters used in the second set of simulations are given brackets.) Thermodynamic equilibrium was assumed to be achieved if no further drift was observed for 400 ps in the mean values of the thermodynamic variables density, pressure, temperature, and potential energy. Each of these equilibrated systems was used to run 21 simulations for the thermodynamic integration. Soft core potentials for Lennard-Jones (LJ) and electrostatic interactions have been used (soft core parameter R ) 1.51 kJ/mol; interaction radius σ ) 0.3 nm) for the TI. The results from these 21 simulations were numerically integrated with the trapezoid rule to obtain finally the value of the excess chemical potential µ. Arguments for28 and against36 the particle insertion algorithm within solvation simulations have been published. The data listed in Tables 1-3 can be used to justify the selection of the TI algorithm instead of Widom's particle insertion method. Table 3 shows that the first solvation shell of the methane molecule contains about 20 water molecules. A simple geometrical model of such a sphere would be a dodecahedron. Table 1 lists the positions of the first maximum of the oxygen-oxygen radial distribution function obtained from MD runs for pure water (section 3.1). If we now assume that the oxygen-oxygen distance in the solvation shell is about the same as the oxygenoxygen distance in liquid water, we can estimate the size of the cavity from the midradius of the dodecahedron.69 Table 1 lists the results from these calculations. The optimal cavity radius can be estimated from the minimum of the LJ function, which describes the methane-water interaction energy. Table 2 lists these minima for the various combination rules. The comparison of the radii listed in Tables 1 and 2 as well as the fact that the water molecules within the virtual water dodecahedron have to be replaced by the methane molecule show that substantial changes in the water structure are to be expected during the

σ  rmin σ  rmin σ  rmin

eq 4

0.344031 0.885361 0.386

0.349915 0.782811 0.393

0.344031 0.885228 0.386

TIP3P TraPPE 0.342809 0.348875 0.885228 0.861647 0.385 0.392

0.349915 0.782693 0.393

0.344778 0.894932 0.387

SPC OPLS-UA 0.343621 0.349368 0.894932 0.872577 0.386 0.392

0.350361 0.796477 0.393

0.344778 0.894932 0.387

SPC/E OPLS-UA 0.343621 0.349368 0.894932 0.872577 0.386 0.392

0.350361 0.796477 0.393

0.342300 0.873347 0.384

MSPC/E OPLS-UA 0.340921 0.347763 0.873347 0.848129 0.383 0.390

0.348911 0.760003 0.392

a

σ, : Lennard-Jones parameters. qO: charge on oxygen. qH: charge on hydrogen. rdf: first maximum in the oxygen-oxygen radial distribution in pure water. rcav: midsphere radius [1/4(3 + x5)‚a ≈ 1.3‚a] of a dodecahedron69 with an edge length a equal to the rdf value of water.

eq 3

TIP3P OPLS-UA 0.342809 0.348875 0.885361 0.861756 0.385 0.392

optimized potential σ 

0.356000 1.013125 0.400

rmin

6 σ in nm;  in kJ/mol; rmin ) x2‚σ in nm (minimum of the resulting Lennard-Jones function and therefore the radius of the optimal methane cavity). a

TABLE 3: Excess Chemical Potential at 300 K for Various Water Potentialsa pot.

rule

SPC SPC SPC SPC SPC SPC/E SPC/E SPC/E SPC/E SPC/E MSPC/E MSPC/E MSPC/E MSPC/E MSPC/E TIP3P TIP3P TIP3P TIP3P TIP3P exp.50 exp.2

a b c d opt. a b c d opt. a b c d opt. a b c d opt.

µ [kJ mol-1] 9.29 9.03 9.55 11.48 8.84 10.01 9.81 10.72 11.74 9.39 9.69 9.73 11.39 12.49 8.95 9.31 9.28 9.84 10.99 8.19 8.09 8.38

N 19.4 ( 0.1 18.3 ( 0.2 20.1 ( 0.2 20.3 ( 0.2 21.0 ( 0.1 19.9 ( 0.1 19.4 ( 0.3 20.0 ( 0.2 20.2 ( 0.1 21.5 ( 0.3 19.9 ( 0.1 19.5 ( 0.1 20.9 ( 0.3 20.7 ( 0.2 21.5 ( 0.2 18.0 ( 0.2 18.4 ( 0.1 19.9 ( 0.2 19.7 ( 0.2 19.7 ( 0.2

a Pot: water-water interaction potential. Rule: combination rule according to eq 4. N: number of water molecules in the first solvation shell. Opt.: optimized methane-water interaction potential.

simulation. The TI algorithm, which allows rearrangements of the solvent structure as the solute molecule blends into the solvent, was therefore preferred by us. All simulations were 500 ps long. The leapfrog algorithm37 (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 (1.2 nm) for the Lennard-Jones and electrostatic forces. The long-range electrostatic forces were calculated with the reaction field (RF)37 (r ) 78.5) approach. Additionally, the dispersion part of the Lennard-Jones interaction

23598 J. Phys. Chem. B, Vol. 109, No. 49, 2005

Konrad and Lankau

beyond the cutoff was corrected for pressure and energy. The neighbor list with a radius of 0.9 nm (1.2 nm) was updated every 10 time steps. The NpT ensembles (1 bar) have been obtained with Parrinello-Rahman pressure38,39 (1 ps) and Nose´-Hoover temperature40,41 (0.1 ps) coupling. The internal degrees of freedom of the water molecules were constrained with the SETTLE algorithm.42 The densities of pure liquid water (eq 2) at various temperatures were obtained from the last 300 ps of simulations containing 216 (820) water molecules. The statistical error within the molecular simulations was computed by repeating the 300 K simulation 17 times and calculating the excess chemical potential µ each time. The mean value of µ was taken as the result of the simulation, while the standard deviation of the mean was used to draw the error bar. The corresponding error bars of the remaining points of the graph were taken from the 300 K simulations and scaled according to the size of µ. Further, to remove this statistical scattering from the plots to enhance their readability, the results from the molecular dynamics simulations were not connected by straight lines, but a smoothing spline interpolation was used to represent the data. The original data will be supplied by the authors upon request. 2.3. Mixing Rules. The difficult part in the computation of the excess chemical potential µ with MD simulations is the selection of a suitable force field. The water molecules were described with four different models using rigid monomer geometries with three interaction sites (Table 1): SPC,43 SPC/ E,43 MSPC/E,44 and TIP3P.45 Two Lennard-Jones type potential functions are commonly used within MD simulations containing methane. The TraPPE (σ ) 0.373 nm,  ) 1.2306 kJ/mol)46 and the OPLS-UA (σ ) 0.373 nm,  ) 1.2309 kJ/mol)47 potentials have nearly identical parameters. An analysis of the MD simulations using either of these potentials showed that both simulations produced the same thermodynamical results within the statistical error margin (subsection 3.4). Therefore, results obtained with the OPLSUA potential are mainly reported here. The nonbonding interaction energies are described in the force fields for methane and water by a Lennard-Jones (6,12) function.

[( ) ( ) ]

u(rij) ) 4ij

σij rij

12

-

σij rij

6

(3)

It is therefore possible to construct a methane-water interaction potential using the mixing rules and the Lennard-Jones parameters of the force fields for the pure substances. The intermolecular Lennard-Jones parameters were calculated either with one of the published mixing rules3,5,11

σii + σjj 2

ij ) xii jj

(ref 37)

(b) σij ) xσii σjj

ij ) xii jj

(ref 10)

(a) σij )

(c) σij )

(d) σij )

σ3ii + σ3jj σ2ii

+

ij )

σ2jj

( ) σ6ii + σ6jj 2

1/6

4ii jj (xii + xjj)2

( )

ij ) 2xii jj

σ3iiσ3jj

σ6ii + σ6jj

(ref 11)

(ref 5)

or our new methane-water interaction potential (subsection 3.5).

Figure 1. Density of SPC/E water as a function of box size. The error bars indicate the rms values of the calculated densities (216 water molecules) taken from the GROMACS output.

The Lennard-Jones parameters of the methane-water interaction potentials are summarized in Table 2. 3. Results and Discussion The first two subsections contain tests of the chosen computational method. This method was applied to the solvation problem in the following two subsections. The development of the optimized methane-water interaction potential was included in this section, because its development was the direct consequence of the results in subsection 3.4. Another subsection is devoted to the results obtained with the optimized interaction potential, and this section closes with the discussion of our results. 3.1. Pure Water. Figure 1 shows the density F of SPC/E water as a function of box size (216 and 820 water molecules, 1 bar ambient pressure) and temperature. The error bars shown are the rms values taken from the GROMACS output from the simulations containing 216 water molecules. The corresponding error bars from simulations containing 820 molecules are approximately half as large as those obtained from the simulations with 216 molecules but still cover both curves. The plots shown in Figure 1 show that a larger computational setup for the simulation does not change the results of the simulations significantly. For comparison, reference values for the density of liquid water (1 bar ambient pressure) are given in Figure 1. The reference values for the density are reproduced by the simulation at approximately room temperature, while the slope of the reference F(T) curve is not reproduced by the experiment. Further, no temperature of maximum density (TMD) can be observed in the calculated F(T) curve. This lack of a TMD has to be partially attributed to the chosen computational method. The same calculation using the TIP5P interaction potential48 (data not shown here) instead of the SPC/E one did not show a TMD in the expected temperature range (∼277 K) but at a much lower temperature (∼250 K), although the TIP5P potential has been parametrized to reproduce the TMD. These problems seem to be caused by the original parametrization of the TIP5P potential as suggested by Lı´sal et al.67 A successful reproduction of the TMD using the TIP5P potential and an MD simulation package has been published by Paschek,28 although the small number of simulations (at 275, 300, 325, 350, and 375 K) makes it difficult to decide whether the observed TMD is the result of the MD simulation or that of the chosen fit algorithm for his data. The F(T) curve in Figure 1 extends far beyond the boiling point of water (373 K) at an ambient pressure of 1 bar. MD simulations are known to describe phase transitions inaccurately37 if periodic boundary conditions are used within the simulation. Consequently, the final density of the system is

Solubility of Methane in Water

J. Phys. Chem. B, Vol. 109, No. 49, 2005 23599

Figure 2. The influence of the simulation pressure on the value of the excess chemical potential µ. The continuous lines represent the results from calculations with a simulation pressure of 20 bar, and the points with error bars represent those obtained for 1 bar.

preset by the density of the initial setup, and it is therefore possible to analyze liquid water at temperatures beyond the boiling point without increasing the pressure. The boiling point of water at an ambient pressure of 10 bar is approximately 453 K.49 The density of water under these conditions is about 887 kg/m3, which has to be compared with 848 kg/m3 (216 water molecules) and 858 kg/m3 (820 molecules). 3.2. Influence of Simulation Pressure on µ. Since phase transitions from liquid to vapor are usually not observed within MD simulations as discussed in the previous section, it might be possible to calculate the Henry constant over the complete temperature range also at an ambient pressure of 1 bar. This assumption is tested in Figure 2. The solid line shows the excess chemical potential calculated at an ambient pressure of 20 bar, while the dots with error bars represent the results from the corresponding calculations at a pressure of 1 bar. The data in Figure 2 demonstrate very well that the value of the excess chemical potential µ is nearly independent of the chosen simulation pressure. Both simulation runs produced the same values for µ within the error margin of the chosen computational method. This result and those shown in Figure 1 demonstrate that it is possible to analyze the complete range of experimental values of KH with a single computational setup. 3.3. Excess Chemical Potential. The majority of simulations at 300 K (Table 3) yielded excess chemical potentials for methane solvation in the range of 9.0 to 9.8 kJ/mol, which is in reasonable agreement with values from experiments ranging from 7.31 to 8.37 kJ/mol2,50,51 and agree well with other simulations (8.04-12.14 kJ/mol).2,22,23,25-28,30,31,52-54,57,62-64 This range of possible results can be reduced further (8.4212.14 kJ/mol) if the selection of simulations concentrates on those with a pressure of approximately 1 bar or a comparable density has been predefined.2,22,23,25-28,31,52-54,57,62-64 The combination of the SPC and OPLS-UA potential with the geometric mixing rules (rule b) produced the best result (9.03 kJ/mol) at a temperature of 300 K compared with the other simulations. Consequently, the SPC and the OPLS-UA potentials have been used to test the effect of the mixing rules on the Henry constant over a larger temperature range. The data in Table 3 show that the value of µ varies only little with the first three combination rules (a, b, c), while results for µ obtained with Halgren's rule (rule d) are always larger than those computed with the first three mixing rules. This effect is also observed on the change of the methane force field from OPLS-UA to TraPPE, as the difference between these two OPLS-UA calculations is usually small (|∆TraPPE µ| < 0.05 kJ/mol). The data for µ in Table 3 show that the geometric mixing rules perform slightly better than the popular Lorentz-Berthelot rules, although the results obtained with both rules are usually very similar.

Figure 3. Influence of the chosen force field on the Henry constant KH as a function of temperature. (a) Variation of the water potential. (b) Variation of the methane and water potentials. (c) Variation of the methane-water potential.

The value of the Henry constant KH depends on the excess chemical potential µ and on the density of the solvent F as shown in eq 2. As F is also obtained from MD simulations, the value of KH is much more sensitive to the simulation’s quality in terms of the water structure than that of µ. 3.4. Henry Constant. Figure 3 summarizes the influence of the chosen force field on the calculated values for the Henry constant KH. Each force field can be subdivided into three groups of interaction potentials (H2O-H2O, CH4-CH4, and CH4H2O), which will be addressed separately in the following discussion. The majority of molecules within the simulation are provided by the solvent. Figure 3a illustrates the influence of the chosen water-water interaction potential of the solvent on the value of KH. Four different water-water interaction potentials (MSPC/ E, SPC/E, SPC, TIP3P) have been used in a series of simulations using the same methane-methane interaction potential (OPLSUA) and the Lorentz-Berthelot mixing rules. The curvatures of all calculated KH(T) curves are nearly the same for all chosen water-water interaction models. Also, the slopes of the curves in the ascending and descending branches of the curves are nearly the same. MD simulations seem to capture the trend in KH(T) regardless of the chosen model, but fail to reproduce the maximum of the curve or the precise value of KH. The maxima of the four calculated curves differ not only in the position of Tmax but also in their values. The primary goal of a further refinement of the simulations has therefore to be the reproduction of this maximum.

23600 J. Phys. Chem. B, Vol. 109, No. 49, 2005 Figure 3b shows the effect of an exchange of the methanemethane interaction (TraPPE and OPLS-UA) on the result of the simulation (TIP3P and SPC/E potentials, Lorentz-Berthelot mixing rules). Both curves obtained with the TIP3P potential and either of the methane-methane interaction potentials are nearly indistinguishable, and the difference between both SPC/E curves is very small. If the statistical noise, which has been removed from the plot by the spline interpolation, in the value of KH is added to Figure 3b, both curves for the SPC/E water potential are also indistinguishable. As the simulation contains only one methane molecule, direct methane-methane interactions between different methane molecules can be ruled out as the source of the small difference between the TraPPE and the OPLS-UA curve. However, an exchange of the methane-methane interaction potential will affect the methane-water interaction potential via the mixing rules. This effect is likely to be small, as the similarity of the TraPPE and the OPLS-UA curves suggest. Figure 3c displays the influence of the mixing rules on the value of the Henry constant obtained with the SPC and OPLSUA potentials. The SPC potential has been chosen for this set of computations, since Berendsen et al.58 showed in contrast to Figure 3a that the SPC potential produces better results for solvation effects than the SPC/E one. It was not possible to improve the quality of the MD runs for KH with SPC water significantly. As observed previously in Figure 3c, the MD simulations reproduce the curvature of the KH(T) curve obtained from experiments but fail to reproduce the maximum. A comparison of Figure 3b and 3c shows that a change of the mixing rules has a much stronger influence on the simulation’s quality than a change in the methane potential. SPC water performs very well in MD simulations of interfaces, which contain water,70,71 whereas SPC/E water yields better results for pure liquid water.43 The chosen computational setup contains far more water than methane molecules, and the properties of the bulk water seem to be more important for the quality of the simulations than the properties of the methanewater interface. Hence, the changes in KH on the variation of the mixing rules are smaller than the changes associated with the change of the water potential from SPC/E to SPC. A comparison of all graphs in Figure 3 shows that the position of the maximum in KH(T) seems to be dominated by the chosen water potential, and it might be therefore possible to create a water potential, which allows a perfect description of the methane solvation with a given methane potential. On the other hand, well-established water-water interaction potentials are provided with many simulation packages, and the creation of yet another water interaction potential seems to be of little use. Another possibility to improve the simulations is given by the mixing rules. The curves shown in Figure 3c indicate that it is possible to move the maximum of KH(T) to a small extent. However, none of the KH(T) curves obtained with the commonly used mixing rules was close to that obtained from experiments. Mixing rules are used to generate interaction potentials between different chemical species from the properties of the pure substances. This approach is likely to work for chemically similar substances but might be prone to error for chemically different species. The Lennard-Jones potential used for the description of the methane-methane interactions has to reproduce long-range attractive as well as short-range repulsive interactions within pure methane. The Lennard-Jones interaction potential in liquid water serves effectively as a short-range repulsive potential keeping the molecules apart to avoid unphysically strong

Konrad and Lankau Coulomb interactions. The TIP3P potential will serve as an example to this interpretation: The oxygen-oxygen distance in the TIP3P water dimer is about 0.274 nm,45,59 which is significantly shorter than the Lennard-Jones potential suggests (σ ) 0.315 nm, Table 1). The Lennard-Jones potential is therefore of a repulsive nature (u(2.74 Å) ) 7.9 kJ/mol) in the (H2O)2 equilibrium geometry. The bonding energy of the TIP3P water dimer is equal to -27.2 kJ/mol, which suggests that attractive interactions are predominately described by the Coulomb interactions and not by the Lennard-Jones potential. Since the repulsive part of the Lennard-Jones potential varies with r-12, short-range water-water interactions will be dominated by the repulsive part while the attractive r-6 part is overpowered by the Coulomb term (r-1) at large oxygenoxygen distances. Calculations with small water clusters showed that the attractive part of the Lennard-Jones potential can be neglected completely in energy calculations as long as the repulsive part of the interaction potential is parametrized appropriately.59 It might be therefore possible that the mixing rules are likely to fail to produce a suitable methane-water interaction potential, as potentials of the pure substances differ significantly in their physical origin. This assumption is tested in the next subsection. 3.5. An Optimized Methane-Water Interaction Potential. It was shown in the previous section that the properties of the KH(T) curve are dominated by the properties of the chosen water-water interaction potential. Further, the analysis of Figure 3b showed that it is still possible to influence the KH(T) by the selection of the mixing rules. This raises two questions: Is it possible to reproduce the experimental KH(T) curve in MD simulations with a given water potential? And if so, how will such an optimized potential change the structure of the bulk water? The new methane-water interaction potential was optimized to reproduce the experimentally observed values of the Henry constant KH with MD simulations. For convenience, the program was written to reproduce the experimentally observed value of the excess chemical potential (µ ) 6.91 kJ/mol) at 277 K. The SPC/E potential was chosen for the description of the water-water interactions, because the SPC/E potential produced the most reliable description of the density of pure water (section 3.1) and the maximum of the KH(T) curve was closest to that observed in experiments (Figure 3a). The initial pair of LennardJones parameters for the optimization were calculated by applying the Lorentz-Berthelot rules to the SPC/E and OPLSUA parameters. Then, the initial parameter set was systematically optimized using a standard simplex algorithm60 with more than 200 sequential MD runs (125 water molecules, T ) 277 K, 1 bar, 50 ps). Since the trajectories calculated during the simplex optimization were too short for statistically reliable predictions of KH, the simplex algorithm did not converge well. This problem is illustrated in Figure 4a, which shows the value of the scoring function F(x) ) x(6.9 kJ/mol-x)2 (simplified to |6.9 kJ/mol - x| in the PERL script for the optimization) as a function of the optimization progress measured in computational steps. As the optimization proceeds, the scoring function F(x) does not become equal to zero, because the statistical error in the MD simulations prohibits the contraction of the simplex. The simplex fluctuates around the minimum and keeps its area constant. The three most promising sets of Lennard-Jones parameters from the simplex optimization were chosen for more precise MD simulations (216 water molecules, 500 ps) from 270 to 450 K. The set of parameters, which reproduced the experi-

Solubility of Methane in Water

Figure 4. Development of an optimized methane-water interaction potential. (a) Convergency of the simplex optimization. (b) Improvement of the KH(T) curve.

mentally observed values of KH best over the complete temperature range, was chosen for a final manual refinement. The value of σ was increased in small steps until further improvement of the simulated KH(T) curve was impossible (σopt ) 0.356 nm). It was not possible to improve the value of  further with this method (opt ) 1.013 kJ/mol). The development of the new methane-water interaction potential is summarized in Figure 4b. The dots show the value of the Henry constant KH at different temperatures obtained from MD simulations using the final methane-water interaction potential. 3.6. Calculations with the Optimized Potential. Figure 5a shows our results for KH(T) computed with the SPC/E potential for water and the OPLS-UA one for methane. The first curve shown was obtained for an ensemble with 216 water molecules and the second with 820 molecules. The methane-water interaction energy has been calculated using our new methanewater interaction potential instead of the mixing rules. Results published previously by other groups, as well as the values of KH(T) obtained from different experiments, are added for comparison. Our results with 216 water molecules reproduce the experimental data reasonably well, and those with 820 water molecules are a nearly perfect reproduction of the experimental curve. The values for KH(T) published by Paschek28 tend to be too large, and those published by Guillot and Guissani65 tend to be too small. Paschek used for his MD simulations the SPC/E waterwater interaction potential, as we did. The temperature Tmax at the maximum of KH(T) is the same in all his simulations with SPC/E water and agrees with our value. This observation is in agreement with our analysis in section 3.4, which suggested that the value of Tmax is dominated by the properties of the chosen water potential. Changing the system size from 216 to 820 water molecules also changed the position of the maximum in our KH(T) curve. This shift in the maximum cannot be observed in the results published by Paschek. The value of Tmax read from the data published by Guillot and Guissani reproduces very well the value of Tmax obtained from experiments. However, the number of published points for KH(T) is far too small to decide whether the observed

J. Phys. Chem. B, Vol. 109, No. 49, 2005 23601 position of the maximum is simply caused by the selection of simulation temperatures or if it is the true maximum. Despite the high quality in Tmax, the value of KH(Tmax) is too small. It should be mentioned also that Paschek28 failed in his attempts to reproduce the results of Guillot and Guissani. The second half of Figure 5a contains data obtained from MC simulations published by Bonifa´cio et al.30 and Errington et al.66 The value of Tmax published by Bonifa´cio et al. agrees with our result, but the values of KH(Tmax) are much smaller than the experimental one. The values for KH(T) published by Errington et al. agree with experiment very well for T > 375 K but fail to reproduce the experimental KH(T) curve for low temperatures. Finally, the comparison of the left and right sides of Figure 5a suggests that MC simulations reproduce the shape of the experimental KH(T) curve slightly less well than MD simulations. The calculations published by Paschek28 suggest that the quality of MD simulations on the methane solvation depends critically on the chosen water potential, which agrees with our results in Figure 3a. On the other hand, the calculations with our new interaction potential (Figure 5a) showed that it is possible to reproduce the experimentally observed values of KH with a less suitable water potential (SPC/E), if a suitable potential was chosen for the methane-water interaction. Paschek suggests in his paper that the structural properties of pure liquid water simulated with the chosen water potentials determine the final result of the simulation. Figure 5b shows the radial carbon-oxygen distribution functions g(rCO) obtained from simulations with the SPC/E and OPLS-UA potentials using the Lorentz-Berthelot rules and our optimized methane-water interaction potential. The most important change is the size of the water cavity. Calculations using the optimized methanewater interaction potentials predict a larger cavity for the methane molecule than the calculations using the LorentzBerthelot mixing rules. This structural change is restricted to the first solvation shell as shown in Figure 5b. The difference between both calculations vanishes for large values of the carbon-oxygen separation (r g 0.45 nm), if the statistical noise of the function is taken into account. The second peak in g(rCO) reflects the second solvation shell. Its structure is caused by interaction between the first solvation shell and the surrounding water molecules. Therefore, the position and shape of this peak are controlled by water-water interactions and consequently by structural effects of the chosen water potential. The structure of the liquid water itself seems to be of lesser importance than the structure of the cavity, which is generated by the watermethane interaction potential. Finally, Figure 5c shows the results for KH from simulations using the OPLS-UA potential for methane and our new methane-water interaction potential with various water-water interaction potentials. The wide variety of results suggests that the optimized methane-water potential cannot be used universally. An optimized potential, like that suggested here, can be used to create a cavity in liquid water simulated with a single force field, but its application to a wider range of problems is problematic. Since this reasoning applies also to mixing rules, any suggestion of a universal set of mixing rules obtained from a similar optimization process for this or other problems seems to be questionable. 3.7. Discussion. Figure 5a shows that it is possible to create a nearly perfect KH(T) with a less than optimal water-water interaction potential just by varying the methane-water interaction potential. To analyze this result in greater detail, a comparison of the optimized Lennard-Jones parameters with

23602 J. Phys. Chem. B, Vol. 109, No. 49, 2005

Konrad and Lankau

Figure 5. Calculations with the optimized interaction potential. (a) Comparison with experiment and published calculations. Left side: MC calculations. right side: MD calculations. P256, P500, P864: Paschek, 256, 500, and 864 H2O molecules.28 G256: Guillot 256 H2O molecules. W216: this work, 216 H2O molecules. W820: this work, 820 H2O molecules. Exp.: experiment.32-34 E250: Errington 250 H2O molecules.66 B256: Bonifa´cio 256 H2O molecules.30 (b) Structural effects. (c) Variation of the H2O potential.

those obtained from the mixing rules is shown in Table 2. The 6 position rmin ) x2‚σ ) 0.4 nm (optimized interaction potential) of the minimum of the Lennard-Jones potential energy curve moves further out during the optimization (σopt ) 0.356 nm, 0.344 nm e σmix e 0.350 nm, SPC/E water), whereas the well (opt) becomes significantly deeper compared with that calculated with the mixing rules (opt ) 1.013 kJ/mol, 0.796 kJ/mol e mix e 0.895 kJ/mol). Trout et al.61 published basis set superposition error (BSSE)-corrected MP2 calculations with large basis sets on the CH4-H2O cluster indicating much larger interaction energies between -2.4 and -4.7 kJ/mol at an oxygen-carbon distance of 0.35 nm. Basic BSSE corrections of the cluster geometry done at our lab suggest carbon-oxygen distances between 0.38 and 0.39 nm. A comparison of the quantum chemical results for an isolated CH4-H2O cluster with the optimized Lennard-Jones potential suggests that the effective methane-water interaction energies in the solvent are weaker than those predicted by the analysis of the isolated system, which might be interpreted as an uncooperative effect. All calculations for KH(T) produce a maximum for KH regardless of the chosen force field. This maximum can be observed in experiments, which suggests that it is not an artifact of the calculation. The value of KH depends according to eq 2 on the density of the pure solvent F and on the excess chemical potential µ. The density of the solvent decreases as the temperature increases (Figure 1), while the excess chemical potential increases (Figure 2). The increase of KH for small values of T can be understood by the fact that small changes in µ will increase the value of KH, because µ is the argument of the exponential function. This increase in KH correlates well with the general chemical experience that the solubility of any gas is reduced by heating the solvent. Figure 1 shows that a temperature increase beyond the solvent boiling point does not change the slope of the density decrease significantly. The maximum in KH(T) has therefore to be attributed to the curvature of µ(T) (Figure 2). An analysis of eq 2 showed that a plot of the exp(µ/RT) term shows the correct

curvature with a maximum. The position of the maximum is simply shifted by the preexponential factor. Table 3 displays the number of water molecules N in the first solvation shell of the methane molecule. The value of N was calculated by integrating the carbon-oxygen distribution functions g(rCO) from its origin to the first minimum. The number of water molecules in this shell is close to 20, which suggests the use of a water dodecahedron as a very first model of the solvation shell. The average carbon-oxygen distance in such a cluster is about 1.3 times larger69 than the oxygenoxygen distance. Table 1 lists the estimated cavity sizes in bulk water on the assumption that the length of the edge of the dodecahedron is equal to the average oxygen-oxygen distance in the liquid. An optimized carbon-oxygen distance of 0.4 nm as required by the optimized methane-water interaction potential requires therefore an oxygen-oxygen distance of 0.31 nm. This distance is significantly larger than the optimized oxygen-oxygen distance (0.27 nm) in the SPC/E water dimer or than the average oxygen-oxygen distance in bulk SPC/E water. An enlargement of the cavity in the solvent therefore increases the attractive methane-water interactions. This increase is only possible on the expense of the water-water interaction energy, and in order to compensate for this energy loss, the methane-water interaction energy has to become more attractive, which can be observed in the larger values for  in the optimized potential. Preliminary 6-31G*/HF calculations on the optimized CH4@(H2O)20 cluster (water dodecahedron with a methane molecule at its center) showed that the water-water bond is about 50 times stronger than the methane-water bond. In a dodecahedral geometry of the first solvation shell, 30 waterwater bonds compete with 20 methane-water bonds, which yields a geometry factor of 1.5. It is therefore possible to assume that the water-water interaction energy is about 75 times larger than the methane-water energy and that the methane-water interaction energy has only a small impact on the geometry of the CH4@(H2O)20 cluster. The small influence of the methanewater energy on the cluster size can be observed in the radius

Solubility of Methane in Water of the cage. The radius increased from 4.02 to 4.03 Å on the insertion of the methane molecule, but the strongest methanewater interaction energy has been observed in an artificially bloated CH4@(H2O)20 cluster at a radius of 4.11 Å. Thermal effects were not part of the quantum chemical simulation, and the comparison of Figure 5b with the 6-31G*/HF results suggests that the thermal motion of the water molecules softens the hydrogen bond network of the water cage and therefore facilitates its inflation. The decrease in the density F of the pure solvent reflects an increase in the average water-water distance as the temperature increases. This increase in the oxygen-oxygen distance increases the likelihood of large water cavities for the methane molecule and consequently enhances the methane-water bonding energy. This attractive interaction energy reduces the value of µ, as the energy to form a suitable cavity is effectively reduced. The decrease in density decreases therefore the value of µ, and KH decreases for large temperatures as shown in Figure 5a. Figure 5b shows the changes in the water structure induced by the new potential. An effective methane-water interaction potential can therefore be used to partially compensate for poor water-water interaction potentials. Since the main effect of the optimized water-methane interaction seems to be a correction of the cavity size computed with a given water potential, the possibilities to use such an optimized methane-water potential with other water-water potentials are limited as shown in Figure 5c. Conversely, the development of an optimized methanewater interaction potential with an ideal water-water potential offers the possibility to check the mixing rules. 4. Conclusions The systematic combination of various water-water interaction potentials with the OPLS-UA methane-methane potential via selected mixing rules showed that the resulting force field failed to reproduce in MD simulations the Henry constant observed in experiments. Starting from the SPC/E water-water and the OPLS-UA methane-methane interaction potentials, it was possible to create an optimized methane-water interaction potential which was successfully used in MD simulations to compute the Henry constant in a temperature range from 250 to 450 K. The analysis of the optimized potential showed that the main effect of the potential optimization is an enlargement of the cavity for the methane molecule and an enforcement of the methane-water interaction energy compared with the potentials obtained from the mixing rules. The influence of the optimized potential on the structure of the solvent is small compared to that on the cavity. The optimization process can therefore be used to compensate minor flaws in the structure of the solution caused by the chosen water potential. As the result of the optimization process, the resulting methane-water interaction potential works satisfactorily only with the originally chosen water-water interaction potential and cannot be used with other water potentials. The same argument has to be true for optimized mixing rules. The application of mixing rules to problems other than the original one seems to be troublesome. Acknowledgment. The authors would like to thank K. Nagorny in Hamburg (D) and I. L. Cooper in Newcastle upon Tyne (GB) for their help at all stages of this work and the JobFoundation for financial support.

J. Phys. Chem. B, Vol. 109, No. 49, 2005 23603 References and Notes (1) Shirts, M. R.; Pitera, J. W.; Swope, W. C.; Pande, V. S. J. Chem. Phys. 2003, 119, 5740-5761. (2) Wescott, J. T.; Fisher, L. R.; Hanna, S. J. Chem. Phys. 2002, 116, 2361-2369. (3) Good, R. J.; Hope, C. J. J. Chem. Phys. 1971, 55, 111-116. (4) DiPippo, R.; Kestin, J. J. Chem. Phys. 1968, 49, 2192-2198. (5) Ewig, C. S.; Thacher, T. S.; Hagler, A. T. J. Phys. Chem. B 1999, 103, 6998-7014. (6) Lin, H.-M.; Robinson, R. L., Jr. J. Chem. Phys. 1971, 54, 52-58. (7) Pen˜a, M. D.; Pando, C.; Renuncio, J. A. R. J. Chem. Phys. 1982, 76, 325-332. (8) White, A.; Tech. Rep. DSTO-TN-0302, Aeoronautical and Maritime Research Laboratory, Weapons System Division, DSTO Aeronautical and Maritime Research Laboratory, P. O. Box 4331, Melbourne Victoria 3001 Australia 2000. (9) Zhang, J.; Kwok, D. Y. Langmuir 2003, 19, 4666-4672. (10) Good, R. J.; Hope, C. J. J. Chem. Phys. 1970, 53, 540-543. (11) Halgren, T. A. J. Am. Chem. Soc. 1992, 114, 7827-7843. (12) Kong, C. L. J. Chem. Phys. 1973, 59, 2464-2467. (13) Kong, C. L. J. Chem. Phys. 1973, 59, 968-969. (14) Kramer, H. L.; Herschbach, D. R. J. Chem. Phys. 1970, 53, 27922800. (15) Kwok, D. Y.; Neumann, A. W. J. Phys. Chem. B 2000, 104, 741746. (16) Pen˜a, M. D.; Pando, C.; Renuncio, J. A. R. J. Chem. Phys. 1982, 76, 333-339. (17) Pen˜a, M. D.; Pando, C.; Renuncio, J. A. R. J. Chem. Phys. 1980, 72, 5269-5275. (18) Wallqvist, A. J. Phys. Chem. 1991, 95, 8921-8927. (19) Smith, D. E.; Zhang, L.; Haymet, A. D. J. J. Am. Chem. Soc. 1992, 114, 5875-5876. (20) Skipper, N. T. Chem. Phys. Lett. 1993, 207, 424-429. (21) Laaksonen, A.; Stilbs, P. Mol. Phys. 1991, 74, 747-764. (22) Slusher, J. T. J. Phys. Chem. B 1999, 103, 6075-6079. (23) Lyubartsev, A. P.; Førrisdahl, O. K.; Laaksonen, A. J. Chem Phys. 1998, 108, 227-233. (24) Lu¨demann, S.; Schreiber, H.; Abseher, R.; Steinhauser, O. J. Chem. Phys. 1996, 104, 286-295. (25) Lin, C.-L.; Wood, R. H. J. Phys. Chem. 1996, 100, 16399-16409. (26) Guillot, B.; Guissani, Y.; Bratos, S. J. Chem. Phys. 1991, 95, 3643-3648. (27) Chipot, C.; Millot, C.; Maigret, B.; Kollman, P. A. J. Phys. Chem. 1994, 98, 11362-11372. (28) Paschek, D. J. Chem. Phys. 2004, 120, 6674-6690. (29) Kirkwood, J. G. J. Chem. Phys. 1935, 3, 300-313. (30) Bonifa´cio, R. P.; Pa´dura, A. A. H.; Gomes, M. F. C. J. Phys. Chem. B 2001, 105, 8403-8409. (31) Boulougouris, G. C.; Voutsas, E. C.; Economou, I. G.; Theodorou, D. N.; Tassios, D. P. J. Phys. Chem. B 2001, 105, 7792-7798. (32) Cramer, S. D. Ind. Eng. Chem. Process Des. DeV. 1984, 23, 533545. (33) Crovetto, R.; Ferna´ndez-Prini, R.; Japas, M. L. J. Chem. Phys. 1982, 76, 1077-1086. (34) Rettich, T. R.; Handa, Y. P.; Battino, R.; Wilhelm, E. J. Phys. Chem. 1981, 85, 3230-3237. (35) Lindahl, E.; Hess, B.; van der Spoel, D. J. Mol. Model. 2001, 7, 306-317. (36) Beutler, T. C.; Be´guelin, D. R.; van Gunsteren, W. F. J. Chem. Phys. 1995, 102, 3787-3793. (37) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids, 1st ed.; Oxford University Press: New York, 1990. (38) Nose´, S.; Klein, M. L. Mol. Phys. 1983, 50, 1055-1076. (39) Parrinello, M.; Rahman, A. J. Appl. Phys. 1981, 52, 7182-7190. (40) Nose´, S. Mol. Phys. 1984, 52, 255-268. (41) Hoover, W. G. Phys. ReV. A 1985, 31, 1695-1697. (42) Miyamoto, S.; Kollman, P. A. J. Comput. Chem. 1992, 13, 952962. (43) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269-6271. (44) Boulougouris, G. C.; Economou, I. G.; Theodorou, D. N. J. Phys. Chem. B 1998, 102, 1029-1035. (45) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926-935. (46) Martin, M. G.; Siepmann, J. I. J. Phys. Chem. B 1998, 102, 25692577. (47) Jorgensen, W. L.; Madura, J. D.; Swenson, C. J. J. Am. Chem. Soc. 1984, 106, 6638-6646. (48) Mahoney, M. W.; Jorgensen, W. L. J. Chem. Phys. 2000, 112, 8910-8922. (49) Handbook of Chemistry and Physics, 72nd ed.; CRC Press: Boca Raton, 1991.

23604 J. Phys. Chem. B, Vol. 109, No. 49, 2005 (50) Ben-Naim, A.; Marcus, Y. J. Chem. Phys. 1984, 81, 2016-2027. (51) Yaacobi, M.; Ben-Naim, A. J. Phys. Chem. 1974, 78, 175-178. (52) Bash, P. A.; Singh, U. C.; Langridge, R.; Kollman, P. A. Science 1987, 236, 564-568. (53) Gallicchio, E.; Kubo, M. M.; Levy, R. M. J. Phys. Chem. B 2000, 104, 6271-6285. (54) Jorgensen, W. L.; Buckner, J. K.; Boudon, S.; Tirado-Rives, J. J. Chem. Phys. 1988, 89, 3742-3746. (55) Lazaridis, T. Acc. Chem. Res. 2001, 34, 931-937. (56) Rick, S. W.; Berne, B. J. J. Phys. Chem. B 1997, 101, 1048810493. (57) Smith, P. E. J. Phys. Chem. B 1999, 103, 525-534. (58) van der Spoel, D.; van Maaren, P. J.; Berendsen, H. J. C. J. Chem. Phys. 1998, 108, 10220-10230. (59) Lankau, T. A Computational Analysis of the Platinum-WaterVacuum Interface, Ph.D. Thesis, University of Hamburg, 2000; http:// www.sub.uni-hamburg.de/disse/398/Diss.pdf. (60) Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; Vetterling, W. T. Numerical Recipes in Pascal, 1st ed.; Cambridge University Press: Cambridge, 1989.

Konrad and Lankau (61) Cao, Z.; Tester, J. W.; Trout, B. L. J. Chem. Phys. 2001, 115, 2550-2559. (62) Shimizu, S.; Chan, H. S. J. Chem. Phys. 2000, 113, 4683-4700. (63) van der Vegt, N. F. A.; Trzesniak, D.; Kasumaj, B.; van Gunsteren, W. F. ChemPhysChem. 2004, 5, 144-147. (64) Villa, A.; Mark, A. E. J. Comput. Chem. 2002, 23, 548-553. (65) Guillot, B.; Guissani, Y. J. Chem. Phys. 1993, 8075-8094. (66) Errington, J. R.; Boulougouris, G. C.; Economouu, I. G.; Panagiotopoulos, A. Z.; Theodorou, D. N. J. Phys. Chem. B 1998, 88658873. (67) Lı´sal, M.; Kolafa, J.; Nezbeda, I. J. Chem. Phys. 2002, 117, 88928897. (68) Widom, B. J. Chem. Phys. 1963, 39, 2808-2812. (69) Weisstein, E. W. Dodecahedron. MathWorld-A Wolfram Web Resource, http://mathworld.wolfram.com/Dodecahedron.html. (70) van Buuren, A. R.; Marrink, S.-J.; Berendsen, H. J. C. J. Phys. Chem. 1993, 97, 9206-9212. (71) Tieleman, D. P.; Berendsen, H. J. C. J. Chem. Phys. 1996, 105, 4871-4880.