J. Phys. Chem. 1982,86,838-847
838
body effects has been fortuitously compensated by errors in the quantum-mechanical energies. There are some limitations in fitting to experimental data. It is necessary that the experimental data adequately represent all important regions of configuration space. This was the difficulty that was found in constructing an EPEN potential for water. The original EPEN/ 1function was fitted to data for only ice and steam, and hence was inadequate for reproducing the important repulsive regions of the water-dimer-pair energy surface. At this point, the concept of fitting to experimental data was abandonedbut only for water.’ Although high-quality, quantummechanical energies were used in the parameterization of EPEN/2,’ the resulting potential did not provide a good simulation of liquid water. Recently, Marchese et al.42have constructed a successful EPEN/ 2 water potential using nearly 300 dimer energies derived from quantum-mechanical calculations and an empirical fit thereof. Besides the previously discussed problems in fitting to quantummechanical energies, it would be excessively time con-
suming to compute such a large number of quantum-mechanical energies for every type of molecule under consideration. The answer to this dilemma may lie in the present simulation. The parameters for saturated amines were obtained by fitting to data which included rotational barriers. Such data would reflect to some extent the type of repulsive interactions missing in the data used for the EPEN/1 water potential. Perhaps, if the known rotational barriers in hydrogen peroxide& were included, it might be possible to obtain a good EPENI2 water potential from only experimental data. It is encouraging that results for systems other than water support the validity of the original EPEN procedure in which parameterization is carried out with readily available experimental data.
(42) F. T. Marchese, P. K. Mehrotra, and D. L. Beveridge, J. Phys.
(43) R. H. Hunt, R. A. Leacock, C. W. Peters, and K. T. Hecht, J . Chem. Phys., 42, 1931 (1965).
Chem., 85, 1 (1981).
Acknowledgment. We are indebted to Drs. W. J. Peer and D. C. Rapaport for many helpful discussions and Dr. D. L. Beveridge for details of the new EPEN/2 water potential prior to publication.
Revised Empirical Potential for Conformational, Intermolecular, and Solvation Studies. 7. Testing of Parameters by Application to Liquid Methane‘ Robert H. Klncaidl and Harold A. Scheraga’ &kef Laboratory of Chemistry, Cornell University, Ithaca, New York 14853 (Received: July 14, 1981; In Final Form: September 7, 1981)
The EPEN/2 potential for methane, parameterized previously on experimental data, is used here in a Monte Carlo simulation of liquid methane in the NVT ensemble at 130 K at a number density of 0.015 A-3 (corresponding to an experimental pressure of 43.4 atm). Simple cubic boundary conditions are used with 125 methane molecules. Considering that the parameters for the methane potential were derived entirely from systems other than methane, the calculated thermodynamic and structural properties of liquid methane are in reasonable agreement with the results of experiment and previous simulations.
Introduction In the accompanying paper,3 we discussed the background for the development of the EPEN/2 potential and applied it to a simulation of liquid ammonia. In the present paper, we carry out a similar simulation for liquid methane. Unlike ammonia and water (which are representative of the parameters for the NH2 and OH groups, respectively, in the EPEN formalism), methane (which is representative of the parameters for the CH2 group) is nonpolar and does not participate in hydrogen bonding. In addition, the parameters for CHI can be used to describe HzO. .CHI intermolecular interactions which (together with the CH4-.-CH4pairwise interaction) are important for studying the nature of aqueous solutions of methane and hydrophobic interactions. The simulation of liquid
-
(1) This work was supported by research grants from the National Science Foundation (PCM79-20279 and PCM79-18336) and from the National Institute of General Medical Sciences of the National Institutes of Health (GM-14312 and GM-25138). (2) NIH Trainee 1978-1981. (3) R. H. Kincaid and H. A. Scheraga,J. Phys. Chem., preceding paper in this issue.
0022-3654/82/2086-0838$0 1.25/0
TABLE I: EPENI2 Parameters for Hydrocarbonsa parameter
value
distance of C-H bonding electron pair from C nucleus
0.8035~
A b
218.44858 (kcali m o l ) I ’2 2.11985 A I 10.550446 (kcal A6/m01)”2
B b C b
Adapted from ref 4 . tron site.
For a saturated carbon elec-
methane provides a test of the parameters for the methane potential. This is especially important since the EPEN/2 parameters for saturated hydrocarbons were obtained4 entirely from systems other than methane. Hence, this simulation provides a test not only of the parameters for methane but also of the extent to which the parameters are transferable from one molecule to another. (4) R. A. Nemenoff, J. Snir, and H. A. Scheraga, J.Phys. Chem., 82, 2504 (1978).
0 1982 American Chemical Society
The Journal of Physical Chemistry, Vol. 86, No. 5, 1982 839
EPENIS Study of Methane
H
\
TABLE 11: Thermodynamic Properties of Liquid Methane a t a Number Density of 0.015 K 3 mean potential energy, kcalimol
\H
H
Flgure 1. Structure of minimum-energy methane dimer.
Potential Function The EPEN/2 potential is described in the accompanying paper,3 and the parameters used for methane are given in Table 1. They pertain to eq 5 of paper 6.3 Values of 1.094 A (ref 5) and 109.47’ were used for the C-H bond length and the H-C-H bond angle, respectively. There are no “lone-pair” sites in the methane potential and, as is customary with ESPNI2, ”bonding-electron” sites are located along the C-H bond. Snir et a1.6 discovered regions on the EPEN/1 waterdimer energy surface, where the pair energy went anomalously to negative infinity. In order to determine whether there are any such singularities on the methane-dimer energy surface, we minimized the energies of 200 randomly selected dimer configurations by using the algorithm of Powell.’ Many different zero-gradient points were found. We have not ascertained whether each of these is a true minimum or saddle point. In Figure 1,however, we report the structure of the global minimum-energy dimer of the methane potential surface. This dimer has a C...C distance of 3.70 A and an intermolecular potential energy of -0.48 kcal/mol. The important result of these minimizations is that no singularities were found in regions of the potential energy surface that are retained in the Monte Carlo selection procedure. On the EPEN/1 water potential surface, the singularities that were found rendered the potential useless for computer simulations of liquid water.4 In contrast, the EPEN/2 potentials for methane, ammonia,3 and water4 appear to be free of such singularities. Monte Carlo Procedure The Monte Carlo calculations were carried out on a Floating Point Systems 120B array processor utilizing a Prime 350 minicomputer as a host machine.* A brief discussion of this system is given in the accompanying pape? and elsewhere? The NVT ensemble was used with 125 methane molecules at 130 K and a number density of 0.015 A-3 (the associated experimental pressure being 43.4 atm).1° These conditions were chosen for comparison with a previous simulation by Murad et a1.,l1 which included calculated radial distribution functions. At the time that our calculations were carried out, no experimental radial distribution functions were available, and the paper of Habenschuss et a1.12 had not yet appeared. The potential was truncated at a distance of 10.1A from a given carbon to any other carbon. The sampling pro(5)G. Herzberg, “Electronic Spectra of Polyatomic Molecules”,Van Nostrand-Reinhold, Princeton, NJ, 1966,p 619. (6)J. Snir, R. A. Nemenoff, and H. A. Scheraga,J.Phys. Chem., 82, 2497 ~.(1978). ~~- --,(7)M. J. D. Powell, Comput. J.,7, 155 (1964). (8)C. Pottle, M. S. Pottle, R. W. Tuttle, R. J. Kinch, and H. A. Scheraga, J . Comput. Chem., 1, 46 (1980). (9)D. C. Rapaport and H. A. Scheraga, Chem. Phys. Lett., 78, 491 (1981). (10)R. D. Goodwin, NBS Tech. Note (US.), No. 653 (1974). (11)S. Murad, D. J. Evans, K. E. Gubbins, W. B. Streett, and D. J. Tildesley, Mol. Phys., 37, 725 (1979). (12)A. Habenschuss,E. Johnson, and A. H. Narten, J. Chem. Phys., 74, 5234 (1981). ~
EPEN/2 ( T = 1 3 0 K ) exptlC ( T = 1 3 0 K ) Murad e t aLe ( T = 127 K)
--1.56‘ - 1.64d -1.66
internal C,,, cal/ energy, (mol kcal/mol deg) .- 0.78’ 8.4’ - 0.87 7.7 -0.85 1.6 I
A correction of -0.1 kcal/mol was added t o the calculated mean potential energy to correct for the neglect of interactions beyond the 10.1-A cutoff. The calculated internal energy was obtained by adding 3 R T t o the mean potential energy, t o include an estimate of the kinetic energy. The calculated heat capacity includes a similar 3R term. Reference 10. The experimental mean potential energy was obtained by subtracting 3 R T from the experimental internal energy for a comparison with the quantity calculated directly by the Monte Carlo procedure. e Reference 1 1 .
’
cedure of Metropolis et al.13 with Euler angles 8, 4, and $ (described by Marion14)was employed. As in the previous paper,36 was sampled as cos 6’ to avoid the necessity of calculating the Jacobian of the configurational integral. The Euler angles 4 and $ were sampled uniformly. Maximum step sizes of 0.25 A, 0.35 rad, and 0.25 were used to alter the Cartesian coordinates, 4 and $, and cos 8, respectively. These step sizes yielded an average acceptance ratio of -50%. Equilibration was carried out over 500 000 attempted moves, and averages were taken over the next lo6 moves. Because of limitations in program memory of the array processor, averages were computed in the same manner as in the accompanying simulation of liquid a m m ~ n i a . ~ Energy averages were taken from each configuration. The carbon-carbon radial distribution function was averaged from every 500th configuration. The C...H and H...H radial distribution functions were obtained from configurations that were saved on magnetic disk or tape after every 10000 attempted moves.
Results A summary of calculated and experimentallo values of the thermodynamic quantities is presented in Table 11. The calculated internal energy was obtained by correcting the mean potential energy and the heat capacity by 3RT and 3R, respectively. These corrections correspond to the classical kinetic energy contributions due to rotation and translation. We do not include quantum corrections to the vibrational components of the energy and the heat capacity, since these contributions were found to be relatively small for both liquid water and liquid ammonia when compared to the typical accuracy of the simulation results. A correction of -0.1 kcal/mol was added to the internal energy to compensate for the neglect of pair energies beyond the 10.1-A cutoff; eq 6 of the previous paper3 was used for this purpose, with a value of 3444 kcal A6/mol for the coefficient C ((obtained empirically for methane15). Since methane is nonpolar, no consideration of neglected dipolar interactions need be made. Corrections to higher multipole interactions are neglected. Using control functions described by Wood,16we estimate the accuracy of the ~~~
~
~
_
_
_
(13)N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, J. Chem. Phys., 21, 1087 (1953). (14)J. B. Marion, “Classical Dynamics”,Academic Press, New York, 1970,p 384. (15)J. 0. Hirschfelder, C. F. Curtiss, and R. B. Bird, “Molecular Theory of Gases and Liquids”,Wiley, New York, 1954,p 1111.
_
840
Kincaid and Scheraga
The Journal of Physical Chemistty, Vol. 86, No. 5, 7982 I .6
3.2 EPEN/2
2.4
M u r a d et.
1
01.
------
~
€PEN12
-
M u r a d et. 0 1 .
------
A
1.6
t I
t I/i
i ,I
O ' 2 1 0 " ' 4.0
O 4 I
" 6.0
8 1. 0
" IO0
4.0
R(i)
Flgure 2. Comparison of the calculated carbon-carbon radial distribution functions from the EPENIP and moleculardynamics" simulations.
60
8.0
10.0
R(i)
Flgure 3. Comparison of the calculated carbon-hydrogen radial distribution functions from the EPEN/2 and moleculardynamics'' simulations. I .6
calculated energies as f3%. It is difficult to make a good EPEN/2 estimate of the accuracy of the calculated heat capacity M u r a d et. 0 1 . - - - - - since it is known that there are as yet convergence problems in the computation of this q ~ a n t i t y . ' ~ J ~ Considering the data of Table 11, and the fact that the potential function was constructed entirely of parameters derived from systems other than methane, we feel that the agreement between the EPEN/2 and experimental energies and heat capacities is reasonable. It is important to note that, while the Williams potentiallg used by Murad 04 et al." was intended to be a transferable one, the functional form chosen by Williams was intended to apply only to hydrocarbon molecules (e.g., no point charges were included). In contrast, the EPENI2 function is applicable 40 60 8.0 I .o to all types of molecules, both polar and n ~ n p o l a r . ~ ~ ~ R(i1 Consequently, point charges are included to provide a consistent formalism; this introduces Coulombic interacFlgure 4. Comparison of the calculated hydrogen-hydrogen radial tions, even in a neutral, nonpolar molecule such as methdistribution functions from the EPEN/P and moleculardynamics" simane. In addition, the EPEN/2 formalism is applicable to ulations. intramolecular interactions such as those that determine rotational barriers, as well as to intermolecular interact i o n ~ .Great ~ ~ ~versatility is demanded of the EPENIP function. It is for this reason that we feel that some small error in quantities such as the energy and the heat capacity can be tolerated at this stage of development of transferable potentials. We chose to compare the calculated radial distribution functions from our simulation (130 K) with those of Murad et al.ll (127 K). Recently, an experimental carbon-carbon P radial distribution function has been obtained by Haben0 shuss et a1.12from X-ray diffraction data. Since the exl6r periment was carried out at 97 K, only qualitative comparisons with our results can be made. Habenshuss et al.,12 however, compared their results with correlation functions calculated from the reference-interaction-site model (RISM),20and good agreement was found. They also compared the RISM results with the molecular-dynamics Flgure 5. Comparison of the calculated coordination numbers from results of Murad et al." Again, good agreement was found. the EPEN/P and moleculardynamics" simulations. This implies, indirectly, that one would expect the molecular-dynamics results of Murad et a1.l1 to be a good representation of liquid methane at 127 K. Further, the molecular-dynamics results yielded the carbon-hydrogen (16)W. W. Wood in 'Physics of Simple Liquids", H. N. V. Temperley, and hydrogen-hydrogen radial distribution functions, J. S. Rowlinson, and G. S. Rushbrooke, Eds., North-Holland Publishing which have not yet been determined by experiment. Co., Amsterdam, 1968,Chapter 5. Figure 2 compares our carbon-carbon radial distribution (17)C. Pangali, M. Rao, and B. J. Berne, Chem. Phys. Lett., 55,413 (1978). ~ - -,. . function with that from the molecular-dynamics calcula(18)M. Mezei, S.Swaminathan, and D. L. Beveridge, J. Chem. Phys., tion. Good overall qualitative agreement is found. The 71,3366 (1979). major discrepancy is that the EPENIB function is some(19)D.E.Williams, J. Chem. Phys., 47,4680 (1967). what too structured, since the first maximum appears to (20) D.Chandler, Mol. Phys., 31,1213 (1976).
:
EPEN/P Study of Methane
be somewhat exaggerated in comparison with the molecular-dynamics results; the same type of discrepancy was observed for the EPENI2 ammonia ~imulation.~ Figures 3 and 4 compare the EPENI2 and moleculardynamics results for the carbon-hydrogen and hydrogenhydrogen radial distribution functions, respectively. The agreement among these functions is quite good. We also computed the coordination number, n(r), by eq 7 of the previous paper: for both our simulation and that of Murad et al.ll Figure 5 depicts n(r) for the EPEN/2 and molecular-dynamicssimulations. Again, we find good agreement, especially near the first minimum of the carbon-carbon radial distribution function. This value is typically used as the average number of neighbors in the first solvation shell. Both our simulation and that of Murad et al." lead to a position of the first minimum in the carbon-carbon radial distribution function near 5.9 A. At this distance, we find n(r) = 12.5 and 12.3 for the EPEN/2 and molecular-dynamics simulations, respectively.
Conclusions EPENI2 has adequately reproduced the internal energy, heat capacity, and radial distribution functions of liquid methane. Considering the early stage of development of the construction of versatile transferable potentials, we feel that the overall agreement between calculated and experimental properties is satisfactory. The only intermolecular properties used in determining the parameters of the EPEN/2 methane potential were the lattice energy and unit cell dimensions of solid n~ e n t a n e .The ~ only other properties used to obtain the parameters for saturated hydrocarbons were intramolecular rotational energies. In the work of Nemenoff et a1.4 the properties of n-pentane were fitted very well, but some internal rotational energies (viz., in ethane and n-butane) were not fitted perfectly. Thus, the disagreement in the calculated internal energy of liquid methane may be a manifestation of the imperfection of the fitting of the internal rotational energies of small hydrocarbon molecules. In contrast, the results of Nemenoff et al.4seem to indicate that these parameters should give better agreement for CH2 fragments in relatively large molecules. It is the latter that would be most important in constructing potential functions for polypeptides and proteins. If experimental intermolecular properties of smaller hydrocarbons were included in the data used to determine the parameters for saturated hydrocarbons, then the resulting potentials would probably be more representative of these smaller hydrocarbons.
The Journal of Physical Chemistry, Vol. 86,
No. 5, 1982 041
From the agreement between the calculated radial distribution functions and those derived from experiment and from previous simulations, it appears that the structural properties of liquid methane are well represented by the EPEN/2 parameters. Such structural information is of importance in the study of aqueous solutions. I t was in the computation of the structure of liquid water that the most severe problems were encountered in the construction of EPEN potentials. Fortunately, methane seems to be free of similar structural inadequacies. The methane-pair energy surface is free from disastrous singularities such as those found on the EPEN/l water-dimer energy surface: The results of this simulation along with the results of the accompanying ammonia simulation are encouraging. The ammonia simulation3shows that the saturated amine parameters of EPEN/2, which were obtained entirely from experimental data, yield good thermodynamic and structural results for the liquid state. The present methane simulation shows that the EPEN/2 function can adequately represent nonpolar liquids as well. Even though EPEN/2 could not simulate liquid water,6 using either experimental data or a limited set of quantum-mechanical data, an EPEN/ 2-type water-water-pair potential (based on a more extensive set of quantum-mechanically derived pair energies) was found to yield good thermodynamic and structural results for liquid water.21 Thus, the three major molecular fragments, NH2, CH2, and OH, all seem well described by EPENI2 functions. Our results for liquid methane and liquid ammonia support the validity of the original EPEN formalism in which experimental rather than quantum-mechanical data are utilized to obtain parameters for potential functions. Fitting to experimental data has several advantages that are discussed in the accompanying paper describing our ammonia ~imulation.~ It is noteworthy that the same functional form can describe such vastly different systems as liquid water, ammonia, and methane. While we have focused attention only on a few systems thus far, EPEN/2 appears to be a useful and versatile formalism capable of describing a variety of molecular systems in the liquid state. Acknowledgment. We are indebted to Drs. W. J. Peer and D. C. Rapaport for many helpful discussions, Dr. D. L. Beveridge for details of the new EPENI2 water potential prior to publication, and Dr. K. E. Gubbins for providing the data from the molecular-dynamics calculations on liquid methane. (21) F. T. Marcheae, P. K. Mehrotra, and D. L. Beveridge, J.Phys. Chem., 86, 1 (1981).
