J . Phys. Chem. 1987, 91, 4334-4341
4334
A Molecular Dynamics Study of Liquid Methanol with a Flexible Three-Site Model G . Pilinkis,+ E. Hawlicka,: and K. Heinzinger* Max-Planck-Institut f u r Chemie (Otto-Hahn-Institut),Mainz, Federal Republic of Germany (Received: February 3, 1987)
A new potential is presented which describes the methanol-methanol interactions on the basis of a flexible three-site model. The intramolecular part of the potential has been derived from spectroscopic data. A molecular dynamics study has been performed with this potential at 286 K. The structural properties of liquid methanol calculated from the simulations are in good agreement with X-ray measurements. The average geometrical arrangement of nearest neighbors and their hydrogen bonding are discussed. The potential describes correctly the gas-liquid frequency shifts of the intramolecular vibrations. Several thermodynamic properties calculated from the simulation compare favorably with experimental results.
I. Introduction In recent years several computer simulations of liquid methanol have been reported.'-5 In all cases the pair potentials employed in the simulations were based on a rigid methanol molecule, which excluded the investigation of changes in molecular geometry and the internal vibrational motions caused by intermolecular interactions. Molecular dynamics (MD) simulations with flexible water models proved already to be capable of reproducing the gas-liquid frequency shift of the intramolecular vibrations and allowed to study the effect of hydrogen bonding on the water molecule ge~ m e t r y .Also, ~ ~ the effects of single ions on these intramolecular properties of water have been calculated from MD simulations.10*" In this paper results of an M D simulation of liquid methanol with a flexible three-site model are reported for the first time. The total potential is separated into an intra- and an intermolecular contribution. The structure of liquid methanol is discussed on the basis of radial distribution functions, the orientation of the molecules, and the geometrical arrangement of nearest neighbors. The gas-liquid frequency shift of the intramolecular vibrations and the self-diffusioncoefficient have been calculated from velocity autocorrelation functions. The differences in the geometry of the methanol molecule between gas phase and liquid are reported. The interaction energies and the heat capacity have been calculated and quantum corrections applied. The results of the simulation are compared with experimental data. Simulations of methanol-water mixtures with the potential presented here have been performed, and the results will be reported in a forthcoming paper.
11. Potential Model The potential describing the interaction between two methanol molecules is based on a flexible three-site model, consisting of an oxygen atom, a hydrogen atom, and the methyl group as a whole. The total potential consists of an intra- and an intermolecular part J'(re,p,pi) = vintra(pi) + V,nter(ru,p) where the p i are the internal coordinates and the ru,Bthe distances between sites in different molecules. The intramolecular potential has been chosen in the form proposed by Carney, Curtiss, and Langhoff (CCL) for water'* vintra(p,)
= CLijpipj
+ CLij,pipjp, + ZLijk/pip,pkp/
with p1 = OH - r & d / r O H I ~ 2 =, (rco - r$o)Ir$o9and ~3 = ( ~ C O H - f f $ o H ) / y t O y . The quantities r&,, re and a;OHare the gas-
2,
phase equilibrium values being 0.9451 1.425 A, and 108.53', respectively. The potential constants used in the simulation are given in Table I. The strong reduction of the number of parameters is based on the approximation that the C-0 stretching and the bending vibrations are harmonic, the interactions between
'Permanent address: Central Research Institute for Chemistry, Hungarian Academy of Sciences, Budapest, Hungary. *Permanent address: Institute of Applied Radiation Chemistry, Technical University (Politechnika). Lodz, Poland.
0022-3654/87/2091-4334$01.50/0
TABLE I: Potential Constants for the Intramolecular Potential Energy of the Methanol Molecule (in kJ/mol) P:
2138.0
p5 PI
2241.7 452.6 -45 22.5
p3
PIp3
450.3 5383.7 661.3
the 0-H and the C-0 stretching vibrations can be neglected, and only interactions between the 0 - H stretching and the C-0-H bending vibrations which are linear in p3 need to be taken into account. The nonquadratic potential constants are taken from the BJH potential for water.6 In order to determine the potential constants for the quadratic terms, several simulations with 32 methanol molecules without intermolecular interactions were performed. In the first simulation the constants given by Timidei and Zerbi13 were used. They were modified in the following simulations until the experimentally known gas-phase frequencies for methanol were reproduced. In addition, a further constraint was imposed in such a way that at an 0-H distance of about 2.2 AI4 the dissociation energy of 457 kJ/molI5 is reproduced. The resulting intramolecular potential energy surface is shown in Figure 1 in the form of a contour plot in rOHand rco for a C-0-H angle of 108.3' (left side). The right half of the figure results if only harmonic contributions are taken into account. The potentials describing the intermolecular interactions can be separated in a Coulombic and non-Coulombic term:
The charges assigned to the three sites are qo = -0.61e1, qH = 0.351e1, and qCH3 = 0.251el. The charges on 0 and H are very similar to the ones in the C F model for water.I6 They yield a
Jorgensen, W. L. J . Am. Chem. SOC.1980, 102, 543. Jorgensen, W. L. J . Am. Chem. SOC.1981, 103, 341. Jorgensen, W. L. J . Phys. Chem. 1986, 90, 1276. Fujimori, H.; Matsumoto, T.; Katayama, M . Mem. Fac. Eng., Hokkaido Univ. 1981, 15, 395. ( 5 ) Haughney, M.; Ferrario, M.; McDonald, I. R. Mol. Phys. 1986, 58, (1) (2) (3) (4)
849. (6) Bopp, P.; JancsB, G.; Heinzinger, K. Chem. Phys. Lett. 1983, 98, 129. (7) Jancsb, G.; Bopp, P.; Heinzinger, K.Chem. Phys. 1984, 85, 377. (8) Toukan, K.; Rahman, A. Phys. Rec. B: Condens. Matter 1985, 31. 2643. (9) Lie, G. C.; Clementi, E. Phys. Rev. A 1986, 33, 2679. (10) Probst, M. M.; Bopp, P.; Heinzinger, K.; Rode, B. M. Chem. Phys. Lett. 1984, 106, 317. (11) Bopp, P. Chem. Phys. 1986, 106, 205. (12) Carney, G. D.; Curtiss, L. A,; Langhoff, S. R. J . Mol. Spectrosc. 1976, 61, 371. (13) Timidei, A,; Zerbi, G. 2.Naturforsch. A , 1970, 25, 1729. (14) Sorbie, K. S.; Murell, J. N. Mol. Phys. 1975, 29, 1387. (1 5) Atkins, P. W. Physical Chemistry; University Press: Oxford, 1982. (16) Stillinger, F. H.; Rahman, A. J . Chem. Phys. 1978, 68, 666.
0 1987 American Chemical Society
The Journal of Physical Chemistry, Vol. 91, No. 16, 1987 4335
Molecular Dynamics Study of Liquid Methanol I
rOH/a-' I .20
1.10
I .oo
'
I
'
'
'
I
'
I
I
'
'
I
' -
I
'OO/A
- ,v
Vh 2.95
i
2.90
i
0.95
t-
0.90
-
0.85
-
0.80
-
2.85
2.83
2.75
L
0.75
-
'
'
I
'
I
'
L
I
'
I
'
'
'
'
I
2.70 9.4
9.6
9.5
9.7
'OH
9.9.10-'
9.8
/A
Figure 3. Contour plot of the potential energy surface of the methanol dimer for 6' = 30'. The energy difference between two adjacent contour lines is 0.25 kJ/mol.
Figure 2. Definition of the angle 6'. TABLE II: Optimized Geometry, Energy, Dipole Moment, and Association Enthalpy for the Methanol Dimer from Various Potential Models and ab Initio Calculations 6:' -AH? ~~
model this work TIP2 (ref 2) OPLS (ref 3) STO-3G(ref 18) 6-31G* (ref 18)
ki:fh
24.69 23.71 28.43 23.30 23.68
"r'
deg
DM, D
kJ/mol
2.85 2.77 2.73 2.74 2.95
30 27 22 48 48
2.78 (1.93)c 3.18 (2.20) 3.60 (2.22) 2.90 (1.51) 3.02 (1.94)
16.67 15.69 20.41 15.28 15.66
-301
2
3
L
5
6
8r/A
7
1
Figure 4. Potential energy for the linear methanol dimer as a function of 0-0 distance with a fixed 0-H distance of 0.954 8, (gas-phase geometry) and 6' = 30'. The dots result from the TIP2 potential.*
"6' is defined in Figure 2. bExperimentalvalues for -AH range from 12.11 to 18.83 kJ/mol." CDipolemoment of the monomer is given in
parentheses. dipole moment of 1.93 D for the gas-phase geometry and lead to approximate agreement with the experimental dimerization energy. The non-Coulombic parts of the 0-0, 0-H, and the H-H potentials are taken from the modified CF model for water' and given by V" = 11 1889/raa6 - 1.045(exp[-4(r - 3.4)*] exp[-1.5(r - 4.5)2]}
+
V'OH
= 26.07/rg2 - 41.79/(1 V'HH = 418.33/(1
+ exp[40(r - l.05)]] 16.74/{1 + exp[5.439(r - 2.2)])
+ exp[29.9(r - 1.968)])
While the non-Coulombic term for the CH3-H interaction is assumed to be zero, the ones for CH3-O and CH3-CH, are described by a (12;6) Lennard-Jones potential with parameters taken from Jorgensen2 (CH3-CH3:a = 3.86 A, c = 0.1883 kJ/mol; CH3-O:a = 3.47 A, c = 0.3713 kJ/mol). This set of potentials is very convenient for the simulations of methanol-water mixtures with the CF model for water. Such investigations are in progress. The potentials described above lead to a hydrogen-bonded linear dimer as the configuration with the lowest energy minimum (Figure 2). The characteristic data for this configuration are given in Table I1 and compared with other methanol potentials as well as with data from ab initio calculations. They agree
-90
-60
-30
0
30
60
90
Figure 5. Potential energy for the linear methanol dimer as a function of 6' with a fixed 0-0 distance of 2.85 A. The dots result from the TIP2
potential.* reasonably well with the TIP2 potential* and the ab initio calculations while they differ significantly from the OPLS p ~ t e n t i a l . ~ The intramolecular 0-H distance is 0.01 1 8,larger than for the monomer. The potential surface for the dimer is shown in Figure 3 in the form of a contour plot as a function of the 0-0 distance and the 0-H bond length for 8 = 30'. In Figure 4 the potential energy for the linear dimer with an 0-H bond length of 0.945 8, (gas-phase geometry) and 8 = 30' is shown as a function of 0-0 distance and compared with the TIP2 potential. The most significant difference exists for the position of the minimum which is found a t a 0.08 8, larger distance compared with the TIP2 (Table 11). Figure 5 shows the dependence of the dimerization energy on 8 for the optimal dimer distance. The curve is very flat around 8 = 30°, which is the most favorable dimer angle and is very similar to the one for TIP2 in this range. The experimental
4336 The Journal of Physical Chemistry, Vol. 91, No. 16, 1987
PZlinkSls et al.
TABLE III: Characteristic Values of the Intermolecular Parts of the Site-Site Radial Distribution Functions g&) XY
R,
rM I
drMI)
00
2.67 1.72 2.16 3.3 1
2.85 1.95 2.45 3.62 2.87 4.15
3.25 3.36 2.98 1.94 0.93 2.06
3H HH
co CH cc
3.63
R2 3.12 2.20 2.88 5.06
rml
3.41 2.65 3.31 4.50 3.27 6.00
g(rml) 0.20 0.07 0.30 1.02 0.81 0.68
for Liquid Methanol" rM2
n(rnl1)
4.9 3.4 5.0 4.9 4.2 7.8
2.0 2.0 2.3 6.1 2.1 12.8
" R i ,rMi, and rmiare the distances (in A) where for the ith time g J r ) is unity, has a maximum, and has a minimum, respectively. The uncertainty is at least =kO.O2A.
Figure 7. Separate contributions of the first (full), second (dashed), and third (dotted) nearest neighbors to the oxygen-oxygen radial distribution function and running integration number.
r/A
Figure 6. Site-site radial distribution functions and running integration numbers for liquid methanol from a simulation at 286 K.
values for the dimerization enthalpy in the gas phase estimated from thermodynamic and spectroscopic data range from -1 2 to -19 kJ/mol." All potentials except the OPLS one lead to values in this range (Table 11). 111. MD Simulation
The system studied consisted of 200 methanol molecules. With an experimental density of 0.7866 g/cm3 a side length of the basic cube of 23.82 8, results. The program which includes the intramolecular three-body interactions is a modified version of the one employed for the simulation of water with a central force modeL7 The Ewald summation was used for all Coulombic interactions and the shifted force potential method19 for the non-Coulombic parts of the potential. The simulation extended over 14000 time steps equivalent to a total ellapsed time of 3.5 ps. The average temperature was 286 K without rescaling. The stability of the total energy of the system was better than 0.1%.
IV. Results and Discussions ( A ) Radial Distribution Functions (RDF).The six site-site RDFs are denoted by gxyand are shown in Figure 6 together with the running integration numbers as defined by
where pm is the number density of molecules. Some characteristic values of the RDFs are given in Table 111. The positions of the intermolecular peaks and the running integration numbers of all six RDFs are consistent with a pref(17) Weltner, W.; Pitzer, K. S. J. Am. Chem. Soc. 1951, 73, 2606. Miller, G. A. J . Chem. Eng. Data 1964, 9, 418. Schuster, Y.A.; Kozlova, V. A.; Grauzhan, V. A . Zh. Fir. Khim. 1978, 52, 799. (18) Tse, Y . C.; Newton, M. D.; Allen, L. C. Chem. Phys. Lett. 1980, 75, 350. (19) Streett, W. B.; Tildesley, D. J.; Saville, G. ACS Symp. Ser. 1980, No. 102, 543.
-0.21
I
Figure 8. Average value of the cos odd as a function of 0-0 distance, where odd is the angle between the dipole moments of two methanol molecules.
erence for a hydrogen-bonded chainlike structure in liquid methanol as deduced from various experimental investigations and computer simulations. There are small differences between the RDFs presented in Figure 6 and the ones reported by Jorgensen simulated with different methanol The positions of the first peaks in goo(r), goH(r), and gHH(r)generally appear at smaller distances. The differences range up to 0.1 8,. This is not unexpected on the basis of the 0-0 distances for the dimers as given in Table 11. Differences in the peak heights may have to be attributed at least partly to the temperature difference of about loo. An average coordination number of about two can be read from noo(r) if integrated up to the first minimum of goo(r) at 3.4 8, (Figure 6 and Table 11). In order to get more detailed information on the nearest-neighbor coordination, the 0-0 RDFs have been calculated separately for the first, second, and third nearestneighbor oxygen atoms. The result is shown in Figure 7 together with the corresponding running integration numbers. For the first, second, and third neighbor the running integration numbers taken again at 3.4 8, are 1.0, 0.8, and 0.2, respectively. This means the existence of branching points, as on the average every fifth oxygen atom has three nearest neighbors. This analysis shows that the simple formula which relates the coordination number noo(rml) = 2(n - l ) / n with the chain length and which has been employed in analysis of X-ray dataz0is incorrect as it neglects the existence (20) Magini, M.; Paschina, G.; Piccaluga, G. .I. Chem. Phys. 1982, 77, 205 1.
Molecular Dynamics Study of Liquid Methanol
The Journal of Physical Chemistry, Vol. 91, No. 16, 1987 4331 1
1
P ( COS e,, 1
kHlki
1.0
05 cos
I
e,,
- 0.5 0 0.5 1 Figure 9. Distribution of cos Odd for nearest-neighbor methanol molecules up to 0-0distances of 3.4 A. -1
of branching points. Figure 7 also shows that the third nearest-neighbor R D F is responsible for the sharp increase of goo(r) around 4 A (Figure 6 ) . ( B ) Orientation of the Methanol Molecules. The average value of cos o d d , where odd is the angle between the dipole moment directions of two methanol molecules, is shown in Figure 8 as a function of the 0-0distance. For two methanol molecules close to each other is about 0.35, in accordance with the formation of a linear hydrogen bond (see Figure 2 ) . Over the range of the first peak in goo(r) the preferential orientation of the dipole moment decreases rapidly with increasing distance. Unlike the water case (see, e.g., Figure 6 in ref 21), there also exists a preferential orientation between 4 and 6 A-the range of the second peak in goo(r). This long-range correlation of the dipole moments is a consequence of the chain formation in liquid methanol. The distribution of cos o d d for nearest-neighbor methanol molecules up to 0-0distances of 3.4 A-the position of the first minimum in gw(r)-is shown in Figure 9. From this distribution an average value for odd of 72" results. The most probable angle is about 4O0, but as the distribution is rather flat around cos odd = 0.5, this angle is not really significant. Beyond 90" the distribution strongly decreases with increasing angle. (C)X-ray Structure Function. Neutron and X-ray diffraction studies have been performed on liquid methanol at room temperature20$22-24 which can be used to check the reliability of the potentials employed in the simulation as far as the structure of the liquid is concerned. The hydrogens in the methyl group give a significant contribution to the neutron scattering pattern. As these hydrogens are not explicitly included in the three-site model, the neutron diffraction results will not be discussed here, but in a forthcoming paper where the results for a six-site model of methanol will be presented. For the comparison with the simulation results the most recent X-ray measurements of Narten and HabenschussZ3will be employed. As the X-ray measurements do not permit the separate calculation of the site-site RDFs, the comparison is restricted to the total structure function and the total RDF. The total X-ray structure function is the weighted sum of the partial site-site structure functions
H ( k ) = E , d k ) h,pW with k = 41r sin O/X where X is the wavelength of the incident radiation and 20 the scattering angle. The weighting functions
c,,(k) = ( 2 - S,,)f,(k) fB(k)/[Ufu(k)I2 depend only on the scattering amplitudes fa(k) of the constituent sites, which vary slightly with k . For methanol the major contributions to the X-ray structure function result from the 0-0, (2-0, and C-C interactions with the approximate average values for c,@(k)of 0.3,0.45, and 0.2, respectively. The partial structure (21) Szlsz, Gy. I.; Heinzinger, K.; Ride, W. 0. Z . Naturforsch. A . 1981, 36, 1067. ( 2 2 ) Montague, D. G.; Dore, J. C.; Cummings, S. J . Mol. Phys. 1984, 53,
1049. (23) Narten, A. H.; Habenschuss, A. J . Chem. Phys. 1984, 80, 3387. (24) Tanaka, Y . ;Olitomo, N.; Arakawa, K. Bull. Chem. SOC.Jpn. 1985,
58, 270.
0
'
2
L
6
8
10
12
k/&-'
Figure 10. Total (below) and distinct (above) structure function for liquid methanol from the simulation (full line) and from an X-ray
measurement (open circles).23
1
,
,
3
,
5
7
,
riA
9
Figure 11. Total intermolecular radial distribution function for liquid methanol from the simulation (full line) and an X-ray measurement
(open circles).23 functions can be calculated from the simulation by Fourier transformations of the site-site RDFs by
kh,,o(k) = 4np,Jmr(g,,,(r) - 1) sin kr d r The experimental structure function published by Narten and H a b e n ~ c h u s swas ~ ~ derived from the measured scattering cross section on the basis of a two-site model (OH and CH,) and normalized to the scattering from a single average site in a methanol molecule. This representation of the structure function depends, however, on the number of sites used to approximately describe the molecule. In order to get a three-site model representation, their measured cross section Z(k) had to be recalculated and renormalized to the scattering of a single methanol molecule using the equation 3
H(k) = [cl(k)-
Cf,2(k)l/tUa(k)l2 Ol=1 a
The total structure functions from the X-ray experiment and the simulation are compared in the lower part of Figure 10. The comparison shows that the flexible three-site model described above leads to a basically correct description of the structure of liquid methanol including intramolecuar geometry and vibrations. A more sensitive check of the intermolecular structure is possible by the comparison of the distinct structure functions where the intramolecular contributions are not taken into account. It can be seen from the upper part of Figure 10 that for the distinct structure function good agreement also exists between simulation and experiment except for a small phase shift at larger k values. In order to better understand the reason for this phase shift, the distinct RDF has been calculated by Fourier transformation from the distinct structure function
Gd(r) = 1 + (1/2a2p,r)
r
k H ( k ) sin kr dk
The Gd(r) from experiment and simulation are compared in Figure 11. It can be seen from the comparison that the phase shift has
4338 The Journal of Physical Chemistry, Vol. 91, No. 16, 1987 TABLE IV: Average Values of the Bond Lengths, Bond Angles, and Dipole Moments of the Methanol Molecules in the Liquid and in the Gas Phase' OlCOHt
liquid gas phase I?
ROH,8, 0.964 (0.028) 0.945
RCO,8, 1.426 (0.021) 1.425
RCH.8, 1.945 (0.050) 1.944
deg 107.4 (3.8) 108.5
DM,D 1.97 (0.1) 1.93
Root-mean-square deviations are given in parentheses.
to be attributed to a shift of the first peak which mainly reflects the 0-0 distance spectrum. Thus, the experiment leads to a most probable 0-0distance of 2.81 A, which is smaller by 0.04 A than the value found from the simulation (Figure 6 and Table 111). ( D ) Geometrical Arrangement of Nearest Neighbors. In order to calculate the particle and time-averaged geometrical arrangement of the nearest-neighbor methanol molecules around a central one, a molecule fixed coordinate system has been introduced where the oxygen atom of the methanol molecule defines the origin, the molecular plane coincides with the yz plane, and the dipole moment points in the negative z direction. The positions of the oxygen atoms of the three nearest-neighbor methanol molecules in this coordinate system have been collected at several hundred time steps spread over the whole simulation time. The resulting projection densities are shown in the left column of Figure 12. In the right column the result of a similar investigation of water is shown for comparison. In the water case the four nearest neighbors are included and the coordinate system is also defined in the insertion. The result for water shows a strong preference for the occupation of tetrahedral sites with a narrower distribution in the hydrogen atom directions than in the lone-pair direction^.^^ In the case of methanol a strong preference exists only for the site opposite of the hydrogen atom. The narrow distribution around the hydrogen atom direction indicates strong hydrogen bonding. In a small but significant number of cases oxygen atoms can be found in the dipole moment direction as can be seen from the small peak in the projection onto the yz plane. Different from water the lone-pair directions in methanol do not lead to sites with preferential occupation. The y z projection shows that the neighboring oxygen atoms prefer the direction between the lone-pair orbitals. They are almost uniformly distribued as far as their x and y coordinates are concerned. ( E ) Intramolecular Geometry and Vibrations. The intermolecular interactions in the liquid will change the gas-phase geometry of the flexible methanol molecule. The average bond lengths, bond angle, and dipole moment together with their ' ~ been root-mean-square deviations I, = ( ( x 2 ) - ( x ) ~ ) ] have calculated from the simulation and are compared in Table IV with the equilibrium gas-phase values. Substantial differences exist for the 0-H bond lengths and the angle which in turn lead to a change in the dipole moment. The distribution of the intramolecular 0-H distances can be seen from the first peak in goH(r) in Figure 6 and the one for the dipole moments is shown in Figure 13. The total spectral density of the intramolecular vibrations in wavenumbers is given by 3 SCH,OH(t)
=
CS"(i4
a=1
with the Fourier transformation
S*('v)= ( 2 m , c / k n J m c u ( t ) cos 2rc5t dt
PBlinkQset al. TABLE V Intramolecular Vibrational Frequencies of Methanol in Wavenumbers" MD expt vib mode gas liquid gas' liquidc CO stretch 1030 1055 1033 1029 COH bend 1331 1407 1345 1420 OH stretch 3693 3342 3679 3337 "The statistical uncertainty of the MD values is at least f 1 0 cm-'. 'Reference 13. Reference 26.
TABLE VI: Average Inter- and Intramolecular Potential Energies, Total Energy, and Constant-Volume Heat Capacity at 286 K As Calculated from the Simulations of Gas-Phase and Liquid Methanol with the Three-Site Model" liquid gas
Ui"W -31.17
K,,. 3.95 3.57
E,,, -16.53 14.27
C,
69 f 10 50 k 2
"The energies are given in kJ/mol and the heat capacity in J / ( K
mol).
number of time averages, and P j ( t ) the velocity of site j of kind a at the time t . m , is the mass of the site and c the velocity of light. The spectra are normalized in such a way that
J P S a ( t )dij = 3 the number of degrees of freedom for each site. The total spectral density for liquid methanol as calculated from the simulation is shown in Figure 14. It is compared with the one for gas-phase methanol (dashed line) which has simply been obtained from a simulation where all intermolecular interactions had been switched off. The positions of the band maxima for both spectra are compared in Table V with spectroscopic data.13-26The band around 500 cm-' which appears only in the spectral density of liquid methanol has to be attributed to librational motions. The band at low frequencies corresponds to hindered translations and is not resolved because of the limited correlation length of 0.25 ps used for cv(t). It can be seen from Table V that the frequencies as well as the gas-liquid frequency shifts calculated from the simulations are in good agreement with the experimental data if it is taken into account that the statistical uncertainty of the MD values is estimated to be at least f 1 0 cm-l. These results show that the simplifications introduced in the intramolecular part of the potential are justified. A more detailed analysis of the intra- and intermolecular motions in liquid methanol will be reported in a subsequent paper. (F)Self-Diffusion Coefficient. The normalized center-of-mass velocity autocorrelation function together with the running selfdiffusion coefficient calculated according to the Green-Kubo relation is shown for liquid methanol in Figure 15. The resulting self-diffusion coefficient of 2.62 X cm2/s seems to have reached its limiting value. It is about 30% higher than the experimental value of 1.97 X 10-j cmz/s at 286 K.27 The discrepancy might be due to the approximation introduced by using a three-site model, as preliminary results of a simulation with a six-site model indicate much better agreement with the experimental value. ( G ) Interaction Energies and Heat Capacity. The intra- and intermolecular potential energies, the total energy, and the constant-volume heat capacity, which has been calculated from the temperature fluctuationsz8by
and the velocity autocorrelation function are reported in Table VI for gas-phase and liquid methanol. The where N , denotes the number of sites ( a = 0, H, CH,), NT the ( 2 5 ) Palinkis, G.; Bopp, P.; Jancsb, G.; Heinzinger, K. Z.Naturforsch. A 1984, 39. 79.
~~
~
( 2 6 ) Falk, M.; Walley, E. J . Chem. Phys. 1961, 34, 1555. (27) Pratt, K. C.; Wakeham, W.A. J . Chem. Soc., Faraday Trans. 2 1977, 73, 997. (28) Lebowitz, J. L.; Percus, J. K.; Verlet, L. Phys. Rec. 1967, (53.250.
Molecular Dynamics Study of Liquid Methanol
The Journal of Physical Chemistry, Vol. 91, No. 16, 1987 4339
2
Figure 12. Densities of the projections of the oxygen atom positions of the three nearest-neighbor methanol molecules (left) and of the four nearest-neighbor water molecules (right) around central ones onto the three planes of coordinate systems as defined in the insertions.
intramolecular potential energy and the heat capacity calculated from the gas-phase simulation give values of 3 / 2 R Tand 6R, respectively. This shows that the anharmonicities in the intramolecular potential energy do not contribute to these properties in the gas phase. The result is different for liquid methanol. The intramolecular energy is 10% higher than the harmonic value, while for the heat capacity the statistical uncertainty of about 15% resulting from the temperature fluctuations does not permit a definite statement. In order to compare the data from the simulation with the measured heat capacity and the heat of vaporization, quantum corrections have to be applied: E:"' = E,
+ A@ + EQ(CH3)
A@ corrects for the intra- and intermolecular vibrations of the, three-site model. It is approximately calculated by treating the
frequencies in the spectral densities ScH,oH(S)(Figure 14) as quantum harmonica1 oscillators (see, e.g., Berens et al.29): AE? = R T
s-
Sc,,,H(>)(Q(S)
- 1) dS
where ~ ( 5 = ) u/2 + u/(eu 1 ) and = h;,+/kT. E Q ( c H ~ ) corrects for the missing nine intramolecular modes neglected by using a three-site model. They are calculated in the Same proximation from experimental gas-phase frequencies13 by 9
EQ(CH3) = R T X Q i ( > ) i=l
(29) Berens, P. H.; Mackay, H. J.; White, G. M.; Wilson, K. R. J . Chem. Phys. 1983, 79, 2375. (30) Ivash, E. V.; Li, J. M.; Pitzer, K. S. J . Chem. Phys. 1955, 23, 1814. (31) Wilhoit, R. C.; Zwolinski, B. J. J . Phys. Chem. Ref.Data, Suppi. 1973, 2.
4340 The Journal of Physical Chemistry, Vol. 91, No. 16, 1987
Pglinkgs et al.
TABLE VII: Comparison of the Corrected Values of the Total Energy and the Heat Capacity from the Simulations of Gas-Phase and Liquid Methanol at 286 K with Experimental Data' AE? EQ(CH3) P A@ G(CH3) cy p t l A C AZp" liquid 29.3 95.2 108.0 -26.1 10.0 5 2 f 10 66.5b 34.5 31.4d gas 30.7 95.2 140.1 -23.1 10.0 31 f 2 35.Ic
-
uFor notations see text. The energies are -given in kJ/mol and the heat capacity in J / ( K mol). * A t 298 K (ref 8). C A t 280 K (ref 30). dAt 298
K (ref 31).
30
' 0
I 00 16
1
18
20
22
21,
26
Figure 13. Distribution of the molecular dipole moments in liquid methanol.
Figure 14. Total spectral densities for liquid (full) and gas-phase (dashed) methanol as calculated from M D simulations.
I
I
2
,
6
4
8 r,
d/A
Figure 16. Average potential energy of two methanol molecules as a function of the 0-0 distances and running integrated interaction energies.
V / k J mol-'
Figure 17. Pair interaction energy distribution for all methanol molecules ~ ~ (dashed) and separately for pairs with hydrogen bond angles I J (for definition see Figure 18) smaller than 5O, lo", 20°, 30°, and 60'.
for the liquid than for the gas phase. The average potential energy of two methanol molecules as a function of 0-0distance is shown in Figure 16 together with the integrated interaction energies as defined by
Figure 15. Normalized center-of-mass velocity autocorrelation function and running self-diffusion coefficient of liquid methanol as calculated from an M D simulation at 286 K.
The correction terms for energy and heat capacity are given in Table VII, and the resulting values are compared with experimental data, where the heat of vaporization has been calculated from the simulations by AH? = EY(g) - Ey(1)
+ RT
The quantum corrections to the total energy are very large, but their contribution to the heat of vaporization is very small as they are almost the same for the gas phase and liquid. The resulting difference between simulation and experiment is less than 10%. The heat capacity for the gas phase calculated from the simulation agrees in the limits of uncertainty with the experimental value while there is a significant difference for liquid methanol. It is difficult to discuss this discrepancy because of the large uncertainty caused by the slow convergence of the temperature fluctuations. Also, the quantum harmonic oscillator approximation employed for the calculation of the correction terms might be less reliable
The energy minimum at 2.85 A coincides with the maximum of goo(r) as can be seen from Figure 6. The two nearest-neighbor methanol molecules contribute 37 kJ/mol or about 2/3 to the total interaction energy. About 80% of these 37 kJ/mol result from Coulomb interactions, while from the total interaction energy of 6 2 kJ/mol only about 60% are Coulombic. Almost all of the remaining 1/3 of the total interaction energy comes from the 10 second nearest neighbors (Figure 6) as indicated by the shallow second minimum in (V(r)) around 5 A. It results from the second peak in gm(r) together with the energetically favorable preferential orientation of the dipole moments in this distance range (Figure 8). The pair interaction energy distribution a s defined by .
P(v) = L N
r
k
c ( 6 ( V - vi,),
i
