Methane in aqueous solution at 300 K

Volume 82, number I

METHANE

IN AQUEOUS

15 August 1981

CHEMICAL PHYSICS LETTERS

SOLUTION

AT 300 K

G. BOLIS and E. CLEMENT1 Department B28, I&?&Corporation. Poughkeepsie. New York 12602, USA

Received 21 April 1981; in final form 19 May 1981

A Monte Carlo simulation for a methane molecule within a cluster of 202 water molecules is presented and compared with previous studies by Scheraga and Beveridge. Whereas the structural information obtained agrees with previous fmdings, the energy points to the need for a very large number of Monte Carlo steps and for a different algorithm to compute the partial molar internal energy.

1. Introduction Methane in water solution can be considered as the prototype of hydrophobic systems. In this work we present a simulation of one molecule of methane fured at the center of a sphere ftied with 202 water molecules constrained within a sphere of 113 1 _&radius. The water molecuIes as well as methane are considered as vibrationless; the statistical positions of water molecules relative to methane are analyzed with the Metropolis [l] variant of the Monte Carlo (MC) method. The simulated temperature is 300 K. The water-water interaction is represented by the two-body potential previously reported [2] and tested in its ability to reproduce accurately the experimentally determined X-ray and neutron scattering for liquid water at 300 K [3]. The methane-water interaction is represented by two different two-body atom-atom potentials, derived from a fitting of computed intermolecular interaction energies for CH4-HZ0 complexes (with H,O at many positions and orientations relative to CH,). These potentids have been derived [4] by fitting interactions energes between CH, and Hz0 obtained (a) from-ab initio SCF MO computations, (b) by correcting the basis set superposition error with the counterpoise method [S] , and (c) by including the dispersion energy correction [6]. Details of the potentials are given elsewhere [4] ; here we point out only that we have selected the following fitting form to repre-

sent the atom-atom

interactions:

B/r12(i, j) - A/r6(i, j) + D/r4(i, j) + Ccl(i) s(i)lrG, i)

(1)

and the often discussed fitting form [73 B’/r12(& j) - Al/&i,

j) + C’q(i)q( j)/r(i, j),

(2)

where A, B, C, D and A’, B’, C’ are fitting constants, r(i, j) is the internuclear distance between atom i on CH4 and atomj on H20, q(i) and s(j) are the net changes on atoms i andj, obtained from SCF LCGO computations for the separated molecules [8] _ We shall refer to these potentials as the four-term (4T) and the three-term (3T) potentials, respectively. The quantum-mechanical computations have been performed with the IBMOL program * and the Monte Carlo simulation with an adaptation of a program originally written by W. Watts ’ and later extensively modified by Roman0 and Clementi [lo] _ It is noted that the 3T and the 4T potentials yield very similar energy surfaces [4] _ By locating the oxygen atom of a water molecule on a sphericaI surface (of radius 4.5 A, and with the methane’s carbon atom at its center) and by minimiz* IBMOL version 7: The latest documentation is available from L. Gianolio, Istituto G. Donegani, S.p.A., 28100 Novara, Italy. * For an early documentation, see ref. [9]. 147

Volume 82, number 1

CHEMICAL

PHYSICS LETTERS

ing the energy for the orient&ion of the water’s hydrogen atoms, with the above potentials one obtains eight shallow minima_ At four minima (with an enerkcal/mol) a water molecule points Its %y of x-0.64 hydrogen atoms towards CH,, for the remaining four minima (with an energy of x-O.69 kcal/mol) the water molecule points its hydrogen atoms away from CH, . The MC sunulation was performed with a forced bias technique to accelerate the convergence of the solvent-solute interaction_ Despite this, the convergence of the small solute-solvent interaction was notably slow, since the solvent-solvent interaction is very dominant. In fig. 1, we report the computed watermethane mteraction energy (solid line) and its standard deviation (dashed line) as a function of the number of the random “moves”, that is, the number of times a water molecule has been randomly selected and displaced; we have not reported the energy for the initial conformations since it is statistically meaningless. (We shall not discuss the water-water interaction as a function of the “moves” number, since it converges much faster than the water-methane interaction.) One can observe that a very large number of “moves” are needed_ For the 4T potential we have used as the starting configuration, the one obtained from the 3T potential after 3.75 X 106 configurations_ In the following we shall denote the average watermethane interaction energy as U(S’). The partial molar internal energy U(S) is given by the sum of U(S’) and AU(W, W’), the latter term being the computed energy difference of the water-water interaction be-

Fig. 1. Water-methane interaction energy as function of the number of “moves”_ The dotted lines report the mean skndard deviation (see text).

148

15 August 1981

tween pure water and solvent water for a sample with the same number of water molecules.

2. Computational

results

Let us indicate the number of oxygen and hydrogen atoms enclosed in two spheres of radius r and r + dr by N(0) and N(H), respectively. in fig. 2 (top inserts) we report dN(O)/dr and dN(H)jdr obtained with the 4T and the 3T potentials. We observe a number of peaks; the first peak for dN(O)/dr starts at e3.20 a and ends at x5.20 A; and for dN(H)/dr the first peak starts at ~2.70 A and ends at x5.30 A. For these peaks the integrals of dN(0) and dN(H) between the above limits yield a value of N(0) and N(H) corresponding to 19.5 oxygen atoms and 38.0 hydrogen atoms. Thus we can state that the first solvation shell contains 19-20 water molecules. The second peak contains =70 water molecules (from ~5.20 to ~8.8 A); thereafter it is impossible to identify additional well-defined shells. (The abrupt ending of dN(O)/dr and dN(H)/dr around 12 II. is due to the spherical boundary conditions imposed on the cluster of 202 water molecules.) At the middle and bottom inserts of fig. 2, we have added the dN/dr values for the half of the spherical cluster of water molecules which contains the -CH, group and for the second half containing the C-H group (for the 3T and 4T potentials). The hemisphere containdng the methyl group is the one with negative values for the 2 axis; therefore the short notation ‘negative” and L‘positie” is used to identify quantities related to the two hemispheres. We postulate that in general, a methyl group has a distribution similar to the one reported in fig. 2 (middle and bottom). Clearly in a molecule CHs-R, the nature of R can be vastly different: for example, when R is a hydroxyl group, the dN(O)fdr for CH, and the dN(H)/dr are expected to be less similar to -CH3 in CH, than for the case when R is another methyl group. As the water-methane interaction is about ten times smaller than the water-water interaction, it seems reasonable to expect only one shell. On the other hand the clothrate hydrogen-bonded structure is ixpetted to be different from the bulk water hydrogenbonded structure. The sharp separation after the second shells indicates that stmcturaliy one can talk of

Volume 82, number 1

CHEMICAL

PHYSICS LETTERS

15 August 1981

9

6 5

5

i

3T-Total

Hydrogen -’ I -J --.,-’

r

c

3T Negative Axis

7

7

6

6

.x.

4T Negative Axis

8

8

5

I

I

3T Positive Axis Hydrogen

*

A

5

4T Positive Axis

4

Fig. 2. Top: dN/dr for the oxygen and hydrogen atoms for the 3T and the 4T potentials. Middle: same quantities considering the half clusters containing the CH3 group, Bottom: same as middle insert, but for the half cluster containing CH.

two soivation she&, even if energetically, we observed only one shell. This can be seen very readily by considering the water-methane statistical distribution of the energy. In fig. 3, we report the energy of those water molecules having the oxygen atom enclosed within two concentric spheres of radius r and r + dr. it is interesting to note that the minimum is very shallow and comparable to the one determined for the 3T
the hydrogen-bonding network of the clathrate; the CH,-H,O attraction being about one-tenth of the H20-H20 attraction, the building up process of the water-water hydrogen-bonded structure will be insensitive to the attraction of CH4 but very sensitive to its kard-core repzdshr Thus, the population of the water molecules in the frost shell is expected to be optbnaily oriented to m aximize the water-water attraction and to avoid the methrme’s hard core. At the middle and bottom inserts of fig. 3, we report the energy distribution for the hemisphere containing either the +H3 or the CH group: we propose the energy profde obtained for the -CH3 group, name149

Volume 82, number 1

CHEMICAL PHYSICS LETTERS

for the CH group (positive axis). In the case of the 4T potential the two curves are more differentiated. These differences are however small. By considering the computed interaction energy between Ar and Hz0 [ 1 l] , we predict that the same distribution and energy will be found by considering a solution of water with Ar as solute. In fig. 4, we present the radial distribution functions for the carbon atom (at the origin of the R axis) with the hydrogen and oxygen atoms of the water molecules_ in fig. 5, we present a not-uncommon water distribution containing 22 oxygen atoms in the clathrate. We find 7 pentagonal and 2 hexagonal structures; in addition we find 1 trimer and an irregular structure with 9 water molecules. The latter has 2 water molecules in common with the trimer; this total of 10 water molecules is interpreted as an inremediare structure in a dynamical process that could yield for example, two hexamers (upon rearrangement)_ Alter-

2' -

1.

Positive

3-r

(B 0

15 August 198 1

1

%__ L Negative

)

I

H 2

3

i

E

6

7

8

s

3-r

IO

, -Hydrogen

-.

-..,

2

I .

-\ ,... \

0

._

._

-

,' .---.._,-*

H

, 0

1

2

A 3

4

F

6

7

8

3

3 I

.

Fig. 3. Top: energyfor water moleculescontainedin the volume dermedby two concentricspheresof radiusr and r + dr. Center: energydistribution(3T potential) for the half-sphere containingeitherthe methyl or the CH group. Bottom: same as centerinsert,but from the 4T potenti&

ly a very shallow minimum with a depth equal very small percent of the water-water energy, acteristic of any hydrophobic interaction. One tice that the 3T potential yields essentially the energy curve for the CH, group (negative axis)

to a as charcan nosame and

Fig. 4. Radial distribution functions for the carbon atom (origin) and the hydrogen and oxygen atoms of the water molecules; top insert, 3T potential; bottom insert, 4T potential.

10

15 August 1981

CHEMICAL PHYSICS LETTERS

Volume 82, number 1

more common structure of the clatbrate; witb regard to fig_ 5 we wish to stress that the clatbrate is not a rigid structure.

3. Discussion Three Monte Carlo simulations have been reported for methane in water solution at room temperature;

the first one is by Darshevsky and Sarkisov 1121, the second by Ovicki and Scheraga 1131, the tbird by Swarninatban et al. [14] _ In table 1, we compare the last two simulations with those reported here; we have excluded the first simulation [12] , since it was obtained with a numerically unstable algorithm and thus the reported data are statisticshy not too reliable. The simulations in refs. 113,141 have been performed with different methane-water two-body potentials, but with the same water-water two-body potential [2,3] . In table 1 we report (in order from top to bottom) some specifications on the CH,+-H20 potentiaLs(identification, reference, lowest energy at the minimum), computational details on the Monte Carlo simulations (number of simulated water molecules, en-

Fig. 5. Representationof the oxygen atoms for a clathrate structure(see text). natively, the process could yield two pentamers by

lransfem?zg two water molecules to the second shell With this example we wish to emphasize the dynamical nature of the clathrate. Below we shah discuss a Table 1 Comparison between Monte Carlo simulations a)

CH4-HsO potentials OS [12] energy niium N water ensemble boundaries

-0.42

RDF (CO) inten. RDF (CHI

inten.

A U(S’) AU(W, w’) US) NW @dhrate)

-0.75

100

124

NPT

NVT

per. 1O-6 N eonfigurations 0.67 RDF (CO) mex 4.2 RDF (CH) max 3.9 RDF (CO) inten. 2.5 RDF (CH) inten. 2.0 RDF (CO) min. RDF (CH) min.

MP [13]

6.0 6.3 0.7

2.0 -0.87 f 0.08 -11 f 15 -11 f 15 23

MO 1131

HM [13]

3T h)

4T h)

-1.17 124 NVT

0.0 124 NVT

-0.70 202 NVT

-0.65 202 NVT

spher. 5.8 3.90

spher. 3.4 3.8

per. 0.75 3.6 c)

per. 0.57 3.5 c)

per. 0.90 3.6 c)

3.70

3.8

2.2 =)

2.8 c)

2.6 c)

2.14 3.45

1.97 2.49

5.6 C)

5.0 c)

0.6 C)

0.6 =)

0.3 =I

0.6 =)

5.40 5.30 0.51

5-6 5.8 0.34

l-71+ 0.29

1.73 + 0.29

5.19 + 0.27

-25.0 i 6.6 -23 56.6 19.35

-235 -21.8 20

* 6.6 zt 6.6

-20.9 -15.7 20

+ 6.6 + 6.6

1.22 -1.02 * 0.06 -10.82 -11.84 20

f 2.87 f 2.83

1.12 -0.03 f 0.05 -6.38 -6.41 20

i 2.87 f 2.99

a) All energies in kcal/mol; all distances in A. b) This work. Cl Radial distribution function reported for carbon and the water molecule center of mass.

151

Volume 82, number 1

CHEMICAL

semble characterization, type of boundaries, number of ‘cmoves”, position and intensity of the radial distribution functions, RDF values at the first maximum and at the first minimum), and the average computed values for U(S), U(S’) and AU(W, W’). One could conclude that even using different potentials and a different number of solvent molecules, substantial agreement in all the reported simulations is found for stmctural predictions, that is, around CH, in water solution and at room remperature, one predicts a clathrate structure with ~20 water molecules. For a nice and “average” representation of the clathrate, see ref. [ 141. However, the next point to note is that the determination of U(S’) was previously reported [ 13,141 as a small positive quantity; in our work, by carrying out the MC simulation more extensively, a small negative value is obtained. In computing U(S) we make use, as previously done [ 13,141, of a precariously small difference AU(W, W’) between two large numbers of nearly comparable value. In addition, U(S’) is dependent on the number of water molecules considered in the simulation. Our simulation has a standard deviation smaller than those previously reported [13,14]. The computed value U(S) = -6.41 -t 2.99 kcal/mol is not far from the experimental result, namely -2.6 kcal/mol [IS] ; we stress however, that by considering intermediate results in our very extended MC simulation, values of AU(W, W’) such that U(S) becomes near to -2.6 kcal/mol have been determined_ Therefore, the direct determination of U(S) as here and previously [13,14] obtained, is not too reliable because of the somewhat unstable algorithm [small difference between large numbers and very slow convergence for W’)] In addition, it is well known from theoretical analyses of many-body perturbations [ 16,171 that threebody correction packs the solvent molecules somewhat differently from what is obtained by considering only two-body interactions. (See ref. [14] for additional comments on the three-body correction.) It is equally known that, in dealing with water, the packing brought about by the three-body interaction tends to lower the energy of the system [18-201. VIith this in mind, we conclude that our results, corroborated by the previous Monte Carlo studies [13,14], are near to what one would obtain by using acact two-body potential; in addition we feel confident in stating that the stntcfuraZ results presented here will change only slightly by including many-body corrections_ Other 152

15 August 1981

PHYSICS LETTERS

properties, however will depend on many-body corrections either notably (for example, pressure) or moderately [for example U(,Sj] . For this reason in this and previous papers [7] we have deferred simulations of those quantities which are known to be strongIy dependent on three-body corrections 1211.

Acknowledgement We acknowledge the contribution of Dr. G. Corongiu with suggestions and discussions. One of us (GB) wishes to acknowledge the IBM Corporation for a post-doctoral fellowship. References [1] N. Metropolis, SW. Rosenbluth, M.N. Rosenbluth, A.H. Teller and E. Teller, J. Chem. Phys. 21 (1953) 1087. [2] 0. Matsuoka, E. Clementi and M. Yoshimine, J. Chem. Phys. 64 (1976) 1351. [3] G-C. Lie, E. Clementi and M. Yoshimine. _J. Chem. Phvs_

64 (1976) 2314. 141 G. Bobs, E. Clementi, H.A. Scheraga, H. Wertz and

--e

--

C. Tosi, J. Am. Chem. Sot., submitted for publication.

151 SF. Boys and F. Bernardi, Mol. Phys. 19 (1970) 553. 161 W. Kofos, Theoret. Chim. Acta 51 (1979) 219. 171 E. Clementi, Lecture notes in chemistry, Vol. 19. Com181 I91 1101

putational aspects for large chemical systems (Springer, Berlin, 1980). R.S. Mulliken, J. Chem. Phys. 23 (1955) 1833,1841. J. Fromm, IBM Technical Report SJR (1973). S. Roman0 and E. Clementi. Gazz. Chim. Ital. 108 (1978) 319.

1111 W_ Kofos, E_ Clementi and G. Corongiu, Intern. J. Quantum Chem. 17 (1980)

775.

WI V-G. Dashevsky and G.N. Sarkisov, Mol. Phys. 21(1974) 1271.

1131 J-C. Owicki and H.A. Scheraga, I. Am. Chem. Sot. 99 (1977)

7403.

1141 S. Swaminathan, S-W. Harrison and D-L. Beveridge, J. Am. Chem. Sot. 100 (1978)

5705.

WI M. Yaacobi and A. Ben-Naim, J. Phys. Chem. 78 (1974) 175.

WI AD. Buckingham, Advan. Chem. Phys. 12 (1967) 107. r171 P. Claverie, in: From diatomics to biopolymers, ed. B. PuBman (Wiley, New York, 1978).

1181 E. Clementi, H. Kktemnacher, W. Koios and S. Romano, T’heoret. Chirn. Acta 5.5 (1980)

237.

WI E. Clementi, W. Kilos, G.C. Lie and G. RanghJno, Intern. J. Quantum Chem. 17 (1980) 377. 1201 J.W. Kress, E. Clementi, J.J. Kozak and MB. Schwartz, J. Chem. Phys. 63 (1975)

3907.

1211 A. Dalgamo, Advan. Chem. Phys. 12 (1967) 143.