J. Phys. Chem. 1993, 97, 9470-9477
9470
Computer Simulation of Liquid Methyl Chloride F. F. Martins Freitas,+ B. J. C . Cabral,* and F. M. S. Silva Fernandes'9t Departamento de Quimica, Faculdade de Cibncias da Universidade de Lisboa and CECUL (INIC), Rua Ernest0 Vasconcelos, Bloc0 C l , Campo Grande 1700 Lisboa, Portugal Departamento de Quimica, Faculdade de Cibncias da Uniuersidade de Lisboa and CFMC (INIC) 2, Av. Prof. Gama Pinto 1699 Lisboa, Portugal Received: March 1 . 1993
Molecular dynamics and Monte Carlo methods were applied, in the isobaric-isothermal NPT ensemble to study several thermodynamic states of CH3Cl in liquid phase. An intermolecular site-site (CH3 and C1) model was used. The thermodynamical properties were studied in the 175 K (mp)-250 K (bp) range. The calculated densities and vaporization enthalpies for the liquid are in good agreement with experimental results. The analysis of the site-site distribution functions indicates the presence in the liquid of some energetically favored relative configurations as antiparallel and linear dimers. The study of the dynamical properties was carried out from the calculation of several time-dependent functions as center-of-mass velocities CV(t)= ( v ( t ) v ( O ) ) , angular Cw(t) = ( J ( t ) . J ( O ) ) , and reorientational c ~ ( = t )(PL(u,(t).uj(t))>( L = 1, 2) correlation functions. These functions were used for the calculations of the respective correlation times. The contributions of the dipolar interactions to the dynamics of liquid CH$l were analyzed, and it was found that they are modest for the translational but are relevant for the molecular angular velocity correlations. The behavior of rWfrom the mp to the bp shows a small decrease with the temperature. We have observed a fair agreement between simulational and experimental results for 71, the correlation time related to C l ( t ) . We also provide a comparison of the computer simulation results with experimental data for the dielectric relaxation of liquid methyl chloride.
1. Introduction Computer simulationsare extremely important to study manybody systems in condensed phases and are essential to discuss properties of liquids.' Molecular dynamics (MD) and Monte Carlo (MC) computer simulations can provide not only relevant information directly comparable with experimentsbut also original information not accessible to conventional experimental methods. In addition, computer simulations are very convenient to assess the adequacy of theoretical approaches and interaction models. In this sense, the validity of a particular intermolecular interaction model should be, in principle, related to the determination of thermal as well as dynamical properties over a set of different thermodynamic states. The discussion on the adequacy of intermolecular potentials is of special importance in the case of molecular liquids. In fact, although some well-known effective intermolecular potentials for liquid simulations have been proposed,Z systematic analysis on their general validity over a wide range of thermodynamic states and its adequacy to correctly reproducedifferentproperties as thermal, dielectric, and dynamics, as well as interfacial properties, are relatively scarce. Another aspect which is crucial for computer simulations concerns the dependence of the results on the sum conventions and on the number of molecules in the system. This point is very important for systems with long-ranged interactions as is usually the case of polar liquids. Thus, in this case, the simulation of systems with different sum conventions and different number of particles is very important in order to discuss finite size effects and to allow any realistic comparison with experimental data. We are reporting in the present study computer simulations of a site-site (CH3 and Cl) intermolecular model for liquid methyl chloride (CHsCI). This system has been the subject of several computer simulation studies. Bbhm et al.3 carried out MD simulations using an ab initio potential. Bigot et al.4 performed
* To whom correspondence should be addressed. t Faculdade de Ciancias da Universidade de Lisboa and CECUL (INIC). t Faculdade de Cikncias da Universidade de Lisboa and CFMC (INIC) 2.
MC simulationsfor CH3Cl using a site-site (TIPS format) model. The simulation of a polarizable model of CH3Cl has been also reported.5 More recently, MD studies of the static dielectric constant6 and of correlation functions' have been carried out. Very recently, the liquid-vapor interface of CH3Cl has been investigated by Monte Carlo simulations.* The resultsofthecomputer simulationsfor CH3C1, in particular those based on the site-site model,4were in good agreement with experimental data9 for some thermal properties, as the density and the vaporization energy. However, except for the study of Millot et a1.,6 all the simulations were limited to one thermodynamic state, which is not necessarily conclusive on the general validity of the intermolecular model. In addition, we are not aware of any study of dynamical properties of CH3Cl based on the TIPS model.4 The present study intends to provide information on the thermal and also on time-dependent properties of CH3Cl over a wide set of thermodynamical conditions and to allow a more detailed discussion on the adequacy of the intermolecular model. The present results are based on Monte Carlo and molecular dynamics computer simulations in the isothermal-isobaric ensemble (MC-NPT and MD-NPT). The interest in using both methods, apart from obvious verification of the agreement between them, may be relevant for this study, since the implementation and practice of MD-NPT for molecular liquids are relatively recent.lO We have selected four thermodynamic states in the liquid including the experimental melting (mp = 175.4K) and boiling points (bp = 249.4 K). We have compared our results with the available experimental information on this system in the liquid phase. The experimental structure for the liquid has been reported." Thermal properties for the liquid are known.9 Absorption-dispersion spectroscopyinfrared (FIR) measurements have been carried out by Gerschel et a1.12 Dielectric relaxation for this system has been discussed by Gerschel et aI.13 2. Computational Details
2.1. Intermolecular Potential. The intermolecular site-site (TIPS format) model has been described in detail elsewhere.*
QO22-3654/93/2Q97-947Q.SQ4.QQ/Q0 1993 American Chemical Society
The Journal of Physical Chemistry, Vol. 97, No. 37, I993 9471
Computer Simulation of Liquid Methyl Chloride
TABLE I: Vaporization Enthalpies for Liquid CH3CI AHH,f kJ mol-' MD T/K 175.4 200.0 220.0
-p
......... .__: _ U _
-1.0
-8.0,
-1.5 0.0
2.0
4.0
6.0
10.0
12.0
249.39
14.0
a
MC
N=108
Ne256
N=500
N=128
exptb
26.0k0.3 26.2 f 0.3O 25.3 A 0.4 25.5 k 0.3" 24.7k0.4 24.8 k 0.4' 23.8 k 0.5 24.0 k 0.4"
26.4f0.2
26.5 f O . 1
25.9 k0.1"
24.3
25.6 k 0.2
25.7 k 0 . 2
25.0k 0.1"
23.5
25.0k0.2
24.8f0.3
24.1fO.la
22.6
24.2f 0.2
24.3 f 0.4
23.3 k0.2'
21.5
Calculated values without Ewald summation. * From ref 9.
r / a 1.5
I
I
I
I
CHJ -CH,-Ci
CI
......... .__u _ t./_ , ,, , , , , , , , , , , , , ,, , , ,, , , , , , , , , , ,
-1.5 -r , , , , , , , , , , , , , ~ , , , , , , , , , , , , , ,, , ,-,
I
--
CH,-CCI-
+
CH,-Cl
1.0
'
0.5
0
-1.5
............................................................
2.0
4.0
6.0
8.0
10.0
12.0
configurationsafter the equilibration phase. Stochastic parameters have been chosen in order to give acceptance rates of about 30-50% for molecular trial moves (translational and rotational) and for volume trial moves at each of the 2500 configurations (a configurationcorresponding to the perturbation of one randomly selected molecule). The cooling process and equilibration/ averaging procedures have been repeated for the set of temperatures reported in Table I. 2.3. Molecular Dynamics Simulations. The MD-NPT calculations were performed with 108, 256, and 500 molecules in a cube with periodic boundary conditions. It is known that the slow decay of coulombic interactions with the intermolecular distance imposes the use of correction procedures in computer simulationsof condensed phases. These procedures are based on lattice summation, as the Ewald sum,1616 or reaction field methods." These corrections seem more important for purely dipolar liquids, for which the orientational order is due essentially to dipolar correlations, as, for example, dipolar hard spheres or Stockmayer fluids. For systems with dipolar quadrupolar interactions, it was arguedl8 that the quadrupolar interactions tend to destroy the structure due to dipole forces and long-range corrections were assumed to play a minor role. For anisotropic site-site models, as the present one for CHsCl, it seems reasonable to assume that the orientational order is determined by a competition between anisotropic site-site and coulombic intera c t i o n ~ .However, ~~ as it will be verified, long-range coulombic correctionsare not very important for thermal properties but are essential in the study of thedynamics. In the present calculations thecoulombicinteractions were treated by a practicalandaccurate procedure, proposed by Adams and Dubey,zothe modified Ewald sum. The method of Andersenzl for constraining the average pressure coupled with Hoover's method22 for constraining the instantaneous temperature was used. Our program has been constructed by adapting a NVE program by Finchamz3which uses the SHAKE algorithm" for constrainingthe intersite distances. The effective "mass of the piston" was chosen as 2 X l e 3 (in reduced units), and the times step was 4 X s. The initial configuration for each run was always a faced-centered-cubic lattice. For each number of particles and for each temperature, two sets of calculations were carried out: one with the minimum image convention and the other with the modified Ewald sum.20 Equilibration over 50 000 time steps was followed by averaging over 40 000 time steps. The calculations were performed on CONVEX machines (C220 and C3410) where a run of 10 000 steps takes 18 h for a system with 500 molecules (C220).
1
14.0
Figure 1. CHzCl TIPS potentials for a parallel (top), perpendicular (center), and linear (bottom) dimer. However, in order to illustrate the Occurrence of some typical relative configurationsof CHpCl dimers in the liquid, we present in Figure 1 the intermolecular potentials for the antiparallel, perpendicular, and linear dimers. The TIPS potential is an effective intermolecularpotential with two interaction sites (CH3 and C1) at a distance of 1.781 A, which is the experimental bond length. The charges on each site (+0.25 and 4 . 2 5 au, respectively) give a dipole moment of 2.13 D, which is higher then the observed gas-phase dipole moment for CH3Cl (1.9 D), taking into account, in an average way, the polarization of the molecule in condensed phases.s.6 The MC and MD procedures are described in detail in the following sections. 2.2. Monte Carlo Simulrtio~. MC-NPT calculations have been carried out with 128 molecules in a cube with periodic boundary conditions. The sum convention for the set of interactions of one central molecule with the others was the minimum image convention. The initial configuration has been taken from a previous study at the boiling point (bp = 249.39 K).s This system has been cooled to 220 K and equilibrated over 5 X 10s configurations. Reported averages are over 106-2 X lo6
3. Structure
The site-site radial distribution functions (RDF's) for CH$l are presented in Figures 2 4 . They reflect the presence of some preferential relative configurationsfor the molecules in condensed phases. Although experimental results on the structure of liquid CHaCl are reported," we are not aware of data on the partial structures. We will discuss in a sequence the RDF's for CH3Cl.
Martins Freitas et al.
9472 The Journal of Physical Chemistry, Vol. 97, No. 37, 1993 TABLE Ik Densities for Liquid C H 3 C l
3.0
dig cm-3
T/K 175.4
N = 108 1.06 f 0.01 1.08f 0.01' 200.0 1.04 f 0.01 1.05f 0.01" 220.0 1.01 f 0.01 1.02 f 0.01' 249.39 0.98f 0.02 0.99f 0.02'
2.0
4
Q
w
1.0
0.0
0.0
2.0
4.0
6.0
r
/Y
10.0
12.0
14.0
Figure 2. CH3-CIradial distribution functions.
,
2.0
I
-
T / K
w.wa 175.4
200 220 249.39
1
l l . o i
0.0 0.0
2.0
4.0
6.0
r
/*f
10.0
12.0
14.0
Figure 3. C H d H 3 radial distribution functions. 2.0
-
T / K
1.5
1
Cl-Cl
249.39
-220
0.5
0.0
4 0.0
I
I 2.0
4.0
6.0
r
/Y
10.0
12.0
*
1.04 0.01 1.04i 0.01 1.05 f 0.03" 1.075 1.02f 0.01 1.02 i 0.01 1.03 f 0.03' 1.04 0.99f 0.01 0.99 f 0.01 0.98f 0.04' 0.985
Calculated values without Ewald summation. From ref 9. orientational freedom between these configurations. It is reduced at higher densities, illustrating again that the orientational motion is more hindered. A clear reinforcement of the role played by antiparallel configurations at higher densities is reflected in the relative height of the two maxima. Our CH3-CH3 site-site RDF is very similar to that one reported by BBhm et ala3 However, for the Cl-Cl RDF we have observed some small differences in the height of the second maxima which is higher than the first one in the study of BBhm et al.3 This result indicates that the presence of linear dimers is enhanced in their simulation. The analysis of the site-site correlation functions suggests that the structure of CH3Cl indicates the presenceof some energetically favorable relative configurationsof two molecules, Le., antiparallel, linear, and perpendicular dimers. The importance of the role played by dipolar interactions at higher densities is evident by a reinforcementof the first peak of the CH3-CI RDF's. In addition, we have verified that the structure is not sensitive to the number of molecules and to the sum convention.
4. Thermai Properties
175.4 200
e.-+a
MD MC N=256 N = 500 N = 128 exptb 1.07f 0.01 1.07 f 0.01 1.09 t 0.04" 1.12
14.0
Figure 4. CI-Cl radial distribution functions.
3.1. C H d . CH3-Cl RDF's (Figure 2) are characterized by a first maximum at 3 . 8 4 . 0 A. This maximum is related with the Occurrence of antiparallel and linear configurations, but the pronounced sharpness of this first maximum at higher densities comes essentially from the contribution of antiparallel configurations. The broad shoulder at 4.8-5.4 A reflects some orientational freedom between antiparallel and linear configurations, and it is related to the Occurrence of perpendicular configurations. This shoulder is reduced at lower temperatures, indicating orientational hindrance. Our CH3-Cl site-site RDF in the 2 W 2 4 9 K range is very similar to the atom-atom C-Cl RDF reported by BBhm et ala3at 273 K. The height and positions of the first (2.0 at 4.0 A) and second (1.2 at 8.0 A) maxima as well as of the shoulder (0.9 at 4.7 A) and first minimum (0.8 at 6.1 A) are in very good agreement. 3.2. CH3-CH3 and CI-CI. The CH3-CH3 (Figure 3) and C1C1 (Figure 4) RDF's are both similar and bimodal, i.e., present two maxima. These maxima are related respectively with the Occurrence of antiparallel and linear configurations. A plateau between the two maxima reflects (as referred above) some
Our results for thermal properties and the correspondingRMS deviations are presented in Tables I and 11, where they are compared to the experimental data? The agreement between NPT-MC and NPT-MD results supports the validity of our code for NPT-MD simulations of molecular systems. The enthalpies of vaporization (Table I) are in good agreement with experiment and reproduce the smallvariation with temperature experimentally observed for CH&1.9 The computer simulation results for the densities (Table 11) are also in very good agreement with the experimental data for the liquid phase. The deviations from experimental values are, in the worst case (1 75.4 K), less than 5%. It can be observed that thermodynamic properties for this system are, within error bars, not sensitive to the the sum convention and to the number of molecules. This conclusion is based on the comparison between results with and without Ewald sum and also on the comparison between results with different numbers of molecules (Tables I and 11).
5. Dynamical Properties 5.1. Translational Mu8ion. The center-of-mass velocity correlation functions c,.(t)are shown in Figure 5 . We can observe that for higher densitiesthese functions present a negativeportion reflecting the tendency to have a component of velocity in the reversal direction after some time (0.4 ps). This reversal is usually interpreted as a rebound from a cage of neighbors.2526 Table I11 shows the translational diffusion coefficientscalculated from the integration of the translational center-of-mass velocity correlation functionsand from the mean-squaredisplacementsvia the Einstein relation. Our results show that the translational diffusion is not significantly sensitive to the number of molecules, but it is dependent on the sum convention. Thus, the results from the simulations with the Ewald sum are systematically higher than the values with the minimum image convention. In order to discuss
The Journal of Physical Chemistry, Vol. 97, No. 37, 1993 9413
Computer Simulation of Liquid Methyl Chloride 1.0
1.0
-
T / K
0.8
*mt
0.6
175.4 200 220 249.39
-
175.4 200 220
0.6
0.5 0.4 v u
t' 0.2
0.0
0.5
1.0
1.5
2.0
t /P. Figure 5. Center-of-mass translational correlation functions for liquid CH3Cl. CV(t) (v(r).v(O)/v(O).v(O)).
TABLE Ilk Translational Diffusion Coefficients D/1 W9mz s-l
N = 108 T/K 175.4
200.0 220.0 249.39
a
1.3 0.8' 1.7 1.ZC 2.3 1.7c 3.3 2.8'
b 1.5 0.8' 1.9 1.3c 2.5 1.8' 3.5 2.9c
N = 256 a
b
1.2
1.4
N = 500 a 1.2
b
ref 3
N = 256
1.3
1.7
1.9
1.9
2.0
2.4
2.6
2.5
2.6
3.4
3.6
3.6
3.7
2.8 (237 K)
From the translational correlation functions. From the mean-square displacement. Calculated values without Ewald summation.
the importance of the coulombicinteractions on the translational diffusion, we have also performed NVT simulations for a system with the charges equal to zero which will be denoted by LJ system. In these simulations we have used the density from the NPT simulation of CH3C1. The diffusion coefficient for a LJ system with 108 molecules at 175.4K is 1.4X 1V m2s-I, a value which is close to that one for CH3Cl with Ewald sum. This result suggests that although the translational dynamics is dependent on the sum convention for the coulombic interactions, for the present intermolecular model, the differences between the translational dynamics of the LJ and dipolar systems are not significant. Unfortunately, we are not aware of experimental values for the diffusion coefficients of CH3Cl. However, our results at 220 K, which are in the 2.3 X 10-9-2.6X lO-9 m2 s-l range, are in good agreement with the value at 237 K (2.8X m2 s-l) predicted by BBhm et al.3 5.2. Angular Correlation Functions. The angular velocity correlation functions C,(t) are presented in Figure 6. These functions are very similar for the set of states in the liquid phase. All of them show a Gaussian decay to zero (0.154.18ps) and present a negative anticorrelation part, which indicates the occurrence of partial reversal of the angular momentum (cage effect) and a longtime negative tail with small oscillations.Some theoretical models have been proposed to describe both orientational and reorientational time correlations in liquids.2'-3O More recently, a Gaussian cage model" has been proposed to describe angular velocity and reorientational correlation functionsof dense liquids. The computer simulation data for CH3Cl at 175.4 K (mp) and 249.4 K (bp) have been fitted to the Gaussian cage model (Figure 6). The model describes correctly the short-time behavior of the angular velocity correlation functions but fails to fit the computer simulation data from 0.3 ps to longer times. Our results for the rotational correlation times are reported in Table IV. We can observe that these times are not sensitive to number of molecules in the system but are dependent on the extension of the correlations. Thus, if we evaluate the correlation
N = 500
T/K
N = 108
a
b
a
175.4
0.0137 0.01 04c 0.0140 0.01 24' 0.0139 0.0136' 0.0119 0.0128'
0.0139
0.00856
0.0129
200.0 220.0 249.39
0.0141
b 0.008 55
0.0143
0.0137
0.00555
0.0136
0.006 35
0.0124
0.00424
0.0125
0.003 75
Correlation extension 2 ps. b Correlation extension 4ps. Calculated values without Ewald summation.
times from the integration of C,(t) up to 2 ps, very similar results are obtained for the liquid phase (0.013ps). Integration up to 4 ps illustrates that orientational correlations are not dead at 2 ps but that small negative correlations persist out to longer times. Our results show a small reduction of 7, with the increase of the temperature. This seems to be related to the oscillatory behavior of the negative long-time tail in Cw(r). As it is illustrated in Figure 6 (inset), comparison between CJt) at 175.4and 249.4 K shows that the first one has a deeper anticorrelation well but decays faster to zero. 7, can be compared with the spin-rotation correlation time T~~ determined by NMR. Unfortunately, for CH3C1, NMR measurements have been carried out only at the solid pha~e,3~-33 and we are not aware of experimental measurements for the liquid. However, NMR measurements of T,, for a polar liquid as HC134showed a similar behavior of rSras function of the temperature. This effect for HC1 has been associated with hydrogen bonding in the liquid. Deviations from the Arrhenius behavior for 7,have been attributed tothecomplexity of molecular orientational motion in liq~ids.3~ Comparison between the rotational correlation times for the LJ system and CHsCl with 108 molecules at 175.4K shows that coulombic interactions contribute to anticorrelate orientational motions. Thus, 7, (from the integration of the rotational correlation function up to 2 ps) for the LJ system is 0.017 ps while it is 0.013 ps for CHsC1. This result seems to reflect the importance of dipolar coulombic contributions to enhance the cage effect in dense liquids, and it is consistent with other results31~35on the role played by quadrupolar contributions to this effect. 5.3. Reorientational Correlation Functions. The total reorientational correlation functions ctt(t)(L= 1,2)are of particular importance for the discussion of the dynamics in dense phases. They are related to absorption and light scattering in condensed
9474
The Journal of Physical Chemistry, Vol. 97, No. 37, 1993
TABLE V
Correlation Times 71' and T + (ps) N = 108 without Ewald sum with Ewald sum
T/K
71' a
71ta 71'
175.4 200.0 220.0 249.39
3.7 2.5 1.9 1.4
5.0 2.3 2.2 1.9
71t
3.4 4.7 2.5 2.3 1.9 2.2 1.5 1.9
TI"
71ta
71'
71t
5.2 3.5 2.3 1.7
5.9 2.6 1.7 1.2
4.5 3.1 2.3 1.7
4.6 2.4 1.8 1.2
7mptc
5.5 (177 K) 4.0(199K) 2.95 (223 K) 2.0(247 K)
with Ewald sum
N = 500
N = 256 T/K 175.4 200.0 220.0 249.39
~
1
3.8 2.6 2.0 1.5
' 71ta ~
4.6 2.5 2.7 1.5
7lSb
71tb
71"
71t"
~ l ' b
~ l t b
3.5 2.6 2.0 1.5
4.7 2.7 2.7 1.5
4.0 2.7 2.1 1.5
5.1 2.3 2.3 1.8
3.6 2.6 2.1 1.5
4.4 2.4 2.3 1.8
From the integration of Cli(t) and Clt(r). From d In Cli(t)/dt and d In Clt(t)/dt. From ref 12. pha~es.3~38CLt(t) is defined as c;(t) = V & ( ~ h j ( O ) ) ) (1) where u is a unit vector along the molecular axis. Clt(t) and CZt(r)can be written as
Martins Freitas et al. 8. From the results with Ewald sum we can observe that 7 1 s < ~ 1 t f othe r thermodynamic states in the liquid phase but also that thevalues at 249.4 K (bp) arevery close. Thedifferences between T ~ Sand ~~t reflect the importance of collective effects.42943 Our results for T I * in the 220-249 K range indicate a good agreement with the value predicted by B6hm et al.3at 237 K (1.9 ps). Absorption-dispersion spectroscopyin the far-infrared (FIR) measurements has been reported by Gerschel et a1.I2over a wide range of thermodynamic states including the solid at 129 K. The experimental correlation times are reported in Table V. We observe a reasonable agreement between our values for and experiment.12 Although we can interpret this agreement as an indication on the adequacy of the TIPS model to describe reorientational correlations in liquid CH3C1, the differences at 200 and 220 K suggests that this model could be improved. In addition, as is illustrated in Figure 8, the statistical errors in Clt(t)arecertainlyreflected on the estimationof the Tltcorrelation times. 5.3.2. Dielectric Relaxation. Experimental data for the dielectric constant of methyl chloride have been reported by Gerschel et al.13 We report in Table VI the Kirkwood factor and the dielectric constant €0 that has been calculated through the relation43 (6)- 1)/3 'YgK
(6)
where
and N
N
From these functions the corresponding correlation times, given by (4) or d In CL(t) dt can be calculated. Our results for Crs(t)and CLt(t) are shown in Figures 7,8, 11, and 12. The respective correlation times are reported in Tables V and VI. To calculate T L from ~ expression 4, we have numerically integrated the bare data up to 2 ps and fitted the long-time behavior of C&) (from 2 ps to infinity) to an exponential function. We shall discuss in the sequence the behavior of these function for the TIPS model of CH3Cl. 5.3.1. Cl(t). Cls(t) and Clt(r) are presented in Figures 7 and 8. The corresponding correlation times are reported in Table V. After a common Gaussian short-time behavior, the reorientational correlation functions decay slowly and persist out for long times. Several points are worth remarking: (1) These functions are not significantlydependent on the number of particles, and in general a good agreement between the correlation times evaluated using expressions 4 and 5 is observed. (2) The reorientationalcorrelation functions are extremely sensitive to the sum convention, namely, minimum image or Ewald summation. This result is related to the pronounced dependence of the total dipole moment correlation function of the system on the boundary conditions, and it has been extensively discussed in the (3) Although the averages have been carried out along 160 ps, C,t(r) has been calculated with significant error bars, as is indicated in Figure 7;
=
y = 47rp~~/9k,T (7 ) p is the density and p is the dipole moment. The Kirkwood factor has been estimated from the fluctuations of the mean-square dipole moment of the sample. Although the agreement with the experimental data is reasonable, we observe that our theoretical dielectric constants are in general underestimated when compared with the experimental ones. However, it is known that many-body induction contributions can play a significant role in the evaluation of this property.44 The real ( ~ ' ( u ) )and imaginary ( ~ " ( v ) parts ) of the dielectric constant as a function of u are compared with the experimental data in Figures 9 and 10, respectively. A fairly good agreement with experiment is observed, and the small differences suggest that the intermolecular model for methyl chloride could be improved. This conclusion is supported by the well-known importance of many-body induction interactions to describe the dielectric relaxation of molecular liquids. 5.3.3. C2(t). CzS(t)and C+(t) are presented in Figures 11 and 12,respectively. The correlation times related to C+(t)and Czs(t) are presented in Table VII. The behavior of these functions indicates that they decay in a different way which is reflected in the differences between the correlation times ~2~ and 72'. These functions are not significantly dependent on the number of molecules but are very sensitive to the sum convention. We can also observe that the agreement between the correlation times calculated from expressions 4 and 5 is not very good. For the set of states considered T+ > 72s. but the differences between them are very small at higher temperatures. Once again, our results for ~ 2 in s the 220-249 K range suggest a good agreement with ~ 2 ~237 a t K (0.65-0.8 ps), predicted by B6hmet ala3Weobserve that, for Cz(t), antiparallel and linear relative configurations between the dimers play the same role. However, perpendicular configurations contribute with negative values and thus reduce the correlation times. This could explain why for this model of CH3C1, at higher temperatures, when the number of relative configurations between the antiparallel and linear ones increases, indicating some reorientational freedom (see the discussion on the RDF's), the correlation times 72t decay rapidly to zero. The reorientational correlation functions CzS(t)may, in principle, be compared with the corresponding correlation function obtained
Computer Simulation of Liquid Methyl Chloride
The Journal of Physical Chemistry, Vol. 97, No. 37, 1993 9415
TABLE VI: Kirkwood Factors and Dielectric Constants ~~
~
aa T/K 175.4 200.0 220.0 249.39 a
~
~
~~
N = 108
N = 256
N=500
N = 108
N = 256
N = 500
€Ob
2.2 1.5 1.4 1.8
2.1 1.4 1.8 1.3
3.3 1.4 1.6 1.3
22.6 13.4 11.8 12.6
21.8 12.7 15.0 9.7
16.8 13.3 13.1 9.5
20.2 (176 K) 18.0 (193 K) 15.8 (213 K) 12.5 (253 K)
Calculated values from MD with Ewald sum. Experimental values from ref 13. 17.5
1.0
0.a
4
I
*
T / K
-2 0 0
*****
193
;;
47 **L 1 1 I 52. .O5
h
~
Q'
-
0.6
-
v .c\
-
10.0
7.5
.1 4
0.4
0.2
1
-
175.4
e+e-a+
2.5
200
220 249.39
5'0
'7
0.5
0.0
\
1.5
2.0
t ps Figure 7 . Cll(t) reorientational correlation functions for liquid CHsCl. 17.5 1
I
1.0
T / K
15.0 -
0.a
12.5 y
0.6
-w 10.0 -
-
7.5 y 0.4
5.0
0.2 -
1
.
0.0
-
-
175.4 200 220 249.39
1
~
1
/
1
1
1
,
/
~
1
,
,
1
1
,
1
,
,
,
,
,
1
1
1
,
,
,
,
,
,
,
,
,
,
,
,
1
,
,
17.5
by a Fourier transform of the vibrational Raman spectrum. We do not know of any Raman experiments for liquid CH3Cl.
1
12.5 I5.O
- 249.34 (4
* * * * * 253 *
i
6. Conclusions Computer simulationsfor a TIPS model of CH3C1have allowed a discussion of relevant aspects for this system in the liquid phase, namely, its structural thermal, and dynamical properties. The behavior of the RDF's has been discussed, and the relationship between the structure of the liquid and the presence of some preferential configurations for the dimers has been clarified. It was also verified that some structural features were in agreement with MD results of B6hm et al.3 Comparison between theoretical results and experimental data shows that this model for CH3Cl provides a correct description of the thermal properties, and in particular, the densities are in very good agreement with experiment over a wide set of thermodynamic states in the liquid phase. Thedynamics of CH3Cl has been extensively discussed. Several time correlation functions and their corresponding correlation times have been determined. We have also reported a detailed discussion about the dependence of the results on the number of particles and on the sum convention. It has been verified that structural and thermal properties are not sigficantly sensitive to the size of the system and to the coulombiclong-range corrections. On the other hand, dynamic properties are very dependent on the
1
T / K
7.5
1
*.\
2.5 1
0.0
I I t 8 1 1 1 1 ,
10'
io'
I81111111
io'
b
II81111,
111118181
iofa io!l ' u / s
I l # \ > # 3 # ,
io"
8
8 ~ 8 8 ~ 1 1 , I
io"
I I # I # , /
io"
Figure 9. Real part of the dielectric constant ~ ' ( u ) (a) : calculated values (N= 500); (b) experimental data." long-range corrections. The simulations of a system with zero charges (LJ system) provided us with some elements to discuss the role played by dipolar interactions in the dynamics of CH3Cl. We have concluded that the translational correlation functions are only modified in a modest way by dipolar coulombic contributions while orientational correlations are more dependent on these interactions. In addition, T~ is very similar in the liquid phase and shows a small decrease when we move from 175.4 K (mp) to 249.4 K (bp). Our results for some dynamical properties such as the translational diffusioncoefficientand ~ 1 ~ a ~z~correlation nd times
9476 The Journal of Physical Chemistry, Vol. 97, No. 37, 1993
1
7.0 J 6.0
Martins Freitas et al.
1
0.8
-
0.6
-
'
v u "
-
0.4
0.2
-
DP880
175.4 200 220 249.39
0.5
0.0
7.0
e.0
1.5
t
ps
Figure 12. CZt(t)reorientational correlation functions for liquid CH3Cl.
,
1
TABLE VII: Correlation Times r i and
(ps)
N = 108 with Ewald sum
TfK 175.4 200.0 220.0 249.39
1
2.0 3.0
72"
72'
7zt
72"
72' a
72'b
1.2 0.9 0.7 0.5
1.7 1.2 1.0 0.7
1.4 1.0 0.8 0.6
1.9 1.4 1.1 0.8
2.1 1.1 0.8 0.6
3.7 2.0 1.0 0.8
1.9 1.3 1.0 0.8
10'
1 1 1 1 1 1 1 /
10'
,
& ,,,,, , io' 10'0 ,
,
4
, , , , ,,
,
N = 256
, , , , , , , ,,
lotl
io"
10"
- 249.94
* * * * * 253 T / K
J
4.0
-
3.0
1
2.0
-
T/K 175.4 200.0 220.0 249.39
1
1
w
10"
s
7.0 7
5.0 6.0
72t
3.1 2.3 1.3 1.0
with Ewald sum ,
0.0
:
72"
:
1.0
without Ewald sum
I
(4
N = 500
72"
7zta
72#
7ztb
~
1.3 0.9 0.7 0.5
1.8 2.1 0.9 0.7
1.4 1.1 0.9 0.7
1.9 1.6 1.0 0.8
1.4 0.9 0.7 0.5
2
' 7zt ~
2.2 1.5 1.2 0.7
7zSb
7zt
1.5 1.1 0.9 0.7
2.1 1.5 1.4 0.9
From the integration of Czl(t) and C2t(t). From d In C2*(t)/dtand d In Czt(r)/dt.
agreement supports the adequacy of the TIPS model to describe reorientational correlations in liquid CHSCl. The differences between theoretical and experimental results for the enthalpies of vaporization and correlation times related to clt(t) suggest that a refinement of this effective model or the use of a nonpairwise additive intermolecular potential for CHsC15 could be of interest.
1.0 1
O.O 10 ~o
8.
v
/
s-l
Figure 10. Imaginary part of thedielectricconstant c"(v): (a) calculated values (N = 500); (b) experimental data." 1.0
II
Acknowledgment. The authors thank Prof. William Steele for interesting and helpful discussions. The Foundation for National Scientific Computing (FCCN, Convex C220) and the Faculty of Sciences Computer Center at the University of Lisboa (FCUL, Convex C3410) are gratefully acknowledged for providing us with computational resources.
0.8 -
0.6
-
0.4
-
References and Notes (1) Allen, M. P.;Tildesley, D. J. In Computer Simulation of Liquids;
-
Clarendon Press: Oxford, 1987. (2) Jorgensen, W. L. J. Phys. Chem. 1986, 90,1276. (3) BBhm, H. J.; Meissner, C.; Aldrichs, R. Mol. Phys. 1984, 53, 651. (4) Bigot, B.; Cabral, B. J. C.; Rivail, J. L. J. Chem. Phys. 1985, 83,
3083. ( 5 ) Cabral, B. J. C.; Rivail, J. L.; Bigot, B. J. Chem. Phys. 1987, 86, 1467. 200 220 (6) Millot, C.; Rivail, J. L.;Diguet, R. Chem. Phys. Lett. 1989,160,228. 249.39 44, 27. (7) Evans, M. W.; Heyes, D. M. J. Mol. Liq. 1989, 0.0 I , (8) Amaral, L. A. N.; Cabral, B. J. C. J. Phys: Condens. Matter 1993, 0.0 0.5 1.5 2.0 5, 1919. t ly p s (9) Gallant, R.In Physical Properties ofHydroearbons; Gulf: Houston, Figure 11. C2*(t)reorientational correlation functions for liquid CHaCl. 1968. (10) Nose, S.;Klein, M . L. Mol. Phys. 1983,50, 1055. (11) Zeidler, M. D. Z . Phys. Chem. (Munich) 1982, 133, 1. are consistent with data from MD simulations carried out using (12) Gerschel, A.; Dimicoli, I.; Jaffre, J.; Riou, A. Mol. Phys. 1976, 32, an ab initio intermolecular potential.' We have observed a fair 619. agreement between theoretical and experimentalcorrelation times (1 3) Gerschel, A.; Grochulski, T.; Kisiel, Z.; Pszczolkowski, L.; Leibler, related to Cl(t) which have been measured by Gerschel.12 This K.Mol. Phys. 1985, 54, 91. 0.2 -
175.4
Computer Simulation of Liquid Methyl Chloride (14) De Leeuw, S.W.; Perram, J. W.; Smith, E. R. Proc. R. Soc. London 1980, A373, 27. (15) De Leeuw. S. W.: Perram. J. W.: Smith. E. R. Proc. R. Soc. London 1980, A373, 57. ' (16) Morrow, T. J.; Smith, E. R. J. Srur. Phys. 1990, 61, 187. (17) Hertzner, A. W.; Schoen, M.; Morgner, H. Mol. Phys. 1991, 73, 1011. (18) Murad, S.;Gubbins, K. E.; Powles, J. G. Mol. Phys. 19%0,40,253. (19) Lomba, E.; Femandez, A.; Martin, C.; Lombardero, M. Phys. Scr. 1991. T38. 79. ~ 0 Adams, ) D.; Dubey, G. J. Compur. Phys. 1987, 72, 39. (21) Andersen, H. C. J. Chem. Phys. 1980, 72, 2384. (22) Hoover, W. G.; Ladd, A. J. C.; Moran, B. Phys. Reu. Lett. 1982,48, 1818. (23) Fincham, D. CCP5 Program Library. (24) Rychert, J. P.; Cicotti, G.; Berendsen, H. J. C. J. Compur. Phys. 1977, 23, 327. (25) Lynden-Bell, R. M.; Hutchinson, D. J. C.; Doyle, M. J. Mol. Phys. 1986,58, 307. (26) Madden, P. Molecular Dynamics Simulationof Statistical Mechanics System. In Proceedingsofthe InrernutionalSchoolofPhysics'Enrico Fermi"; Cicotti. Hoover. Eds.; North-Holland Amsterdam, 1986. (27) Gordon, R. G. J. Chem. Phys. 1966,44, 1830. ~
The Journal of Physical Chemistry, Vol. 97, No. 37, 1993 9477 (28) Fixman, M.; Rider, K. J. Chem. Phys. 1%9,51, 2425. (29) Detyna, E.; Singer, K.; Singer, J. V. L.; Taylor, A. J. Mol. Phys. 1980.41. 31. LGden-Bell, R. M.; McDonald, I. R. Mol. Phys. 1981, 43, 1429. Lynden-Bell, R. M.; Stale, W. A. J. Phys. Chem. 1984,88,6514. Watton, A.; Pratt, J. C. J. Mugn.Reson. 1976, 64, 296. Eguchi, T.; Chihara, H. J. Mu@. Reson. 1988, 76, 143. Krynicki, K.; Powles, J. G. J. Mugn. Reson. 1972, 6, 539. Steele, W. A.; Street, W. B. Mol. Phys. 1980, 39, 299. van Konynenburg, P.; Steele, W. A. J. Chem. Phys. 1972,56,4776. Brot, C . In Molecular Motions in Liquids; Lascombe, J., Ed.; D. Boston, 1974. Rosenthal, L.C.; Strauss, H. L. J. Chem. Phys. 1975,64,282. Neumann, M.; Steinhauser, 0.;Pawley, G. S . Mol. Phys. 1984,52, Kusalik, P. G. Mol. Phys. 1991, 73, 1349. Belhadj, M.; Alper, H. E.; Levy,R. M.Chem. Phys. Leu. 1991,179, Pollock, E. L.; Alder, B. J. Phys. Reu. Lerr. 1981, 46, 950. Leveaque, D.; Weir, J.-J. J. Chem. Phys. 1983, 79, 917. Caillot, J. M.; Levespue, D.; Weis. J.; Kusalik, P. G.; Patev, G. N. Mol. Phys. 1985, 55, 65.
