7378
J . Phys. Chem. 1990, 94, 7378-7385
Molecular Dynamics Simulation of Liquid Mixtures of AcetonitrHe and Chloroform H. Kovacs, J. Kowalewski, and A. Laaksonen* Division of Physical Chemistry, Arrhenius Laboratory, University of Stockholm, S - 106 91 Stockholm, Sweden (Received: December 19, 1989: In Final Form: March 30, 1990)
The molecular dynamics simulation technique is employed to simulate liquid chloroform and liquid mixtures of chloroform and acetonitrile at four molar fractions. A comparison with an earlier NMR investigation at this laboratory is carried out in order to verify the validity of the simulation prerequisites, whereafter changes in the structural and dynamic properties are examined as a function of mixture composition. An analysis of the atomatom radial distribution functions and interaction energies provides evidence of increased structural order and enhanced interactions between the unlike components upon dilution of acetonitrile. The acetonitrile molecules do not show any clustering tendency, although their long-ranging electrostatic interactions are strong. Chloroform seems to function more like an inert solvent. The translational and reorientational motions of acetonitrile are highly anisotropic, whereas the chloroform motions are fairly isotropic and can be described reasonably well by the small-step diffusion model. The dynamic properties of the two mixture components exhibit only minor variation with concentration.
Introduction TABLE I: Simulation Parameters mol%CH3CN 0 10 22 50 78 Only a few molecular dynamics (MD) simulations on mixtures time step, fs 1 1 1 1 1 of simple organic compoundsI4 have been reported to date, in cutoff radius, A 16.55 16.14 15.98 15.26 14.56 contrast to aqueous mixtures which have attracted considerable density, g/cm3 1.40 1.40 1.33 1.20 1.00 interest. However, in a system of polar organic substances one temp, K 288 f 15 302 f 24 299 & 10 293 f 8 301 f 6 might expect to trace weak or intermediate intermolecular interactions. Moreover, the detailed dynamic information obtained TABLE 11: Potential Parameters Used for Chloroform in the during the simulation can be used to select among different Simulation motional models for simple liquids. site Cii(kB/K) UiiI2, A 4ikl The acetonitrile-chloroform system has been subject to several C 51.0 1.6 0.32 thermodynamic s t u d i e P and quantum mechanical calc~lations.~~l~ CI 175.0 1.75 -0.14 Binary solutions of acetonitrile and chloroform have also been H 13.4 1.375 0.10 investigated by NMR spectroscopy"-16and Raman spectroscopy1' the present work we extract new information in the form of radial in order to elucidate the association phenomena and dynamic distribution functions and time correlation functions. We have properties of the mixture components. We chose the present liquid used the mixture concentration as a variable in order to reveal mixture system in order to complement our previous NMR work.I8 trends and subtle changes. Conclusions presented in that article were based on measurements Both pure acetonitrile and pure chloroform properties were of the translational diffusion coefficients and rotational correlation analyzed by Evans19*20 in 1982 by computer simulations using times both around and perpendicular to the symmetry axis in these effective pair potentials. Another simulation on chloroform with two symmetric top molecules. Viscosities were also determined. an effective pair potential was performed a couple of years later If the measured dynamic properties are reproduced reasonably by Dietz and Heinzinger?I while Bohm and Ahlrichs22published well through simulation, there is good reason to believe that the simulation results for chloroform acquired through a potential model liquids are comparable to the real mixtures to the extent derived in part from ab initio data. that, based on the MD data, conclusions may be drawn, for At about the same time Bohm, McDonald, and Maddenz3 instance, about the liquid structure. Such detailed molecular published an article in which they presented an effective pair information cannot be inferred from the NMR experiments. In potential for acetonitrile. That potential has been shown to yield structural and dynamic properties in good agreement with experimental data. Lynden-Bell and co-workers have used it later ( I ) Schoen, M.; Hoheisel, C.; Beyer, 0. Mol. Phys. 1986, 58, 699. in a study of molecular motion24 and several other investiga(2) Fincham, D.; Quirke, N.; Tildesley, D. J. J . Chem. Phys. 1986, 84, t i o n ~ * of ~ -liquid ~ ~ acetonitrile. We have used the same potential 4535; Fincham, D.; Quirke, N.; Tildesky, D. J. J. Chem. Phys. 1987,87,6117. for acetonitrile, and thus considered it superfluous to repeat the (3) Maliniak, A.; Laaksonen, A. Mol. Phys. 1987,62,489. Maliniak, A.; pure acetonitrile simulation despite the somewhat lower temLaaksonen, A.; Kowalewski, J.; Stilbs, P. J. Chem. Phys. 1988, 89, 6434. (4) Mittag, U.; Samios, J.; Dorfmiiller, Th. Mol. Phys. 1989, 66, 51. perature of the Lynden-Bell studies (291 K versus 298 K). Re( 5 ) Kreglewski, A. Bull. Acad. Pol. Sci. 1965, 13, 723. cently Jorgensen and BriggsB devised a simpler three-site potential
Handa, Y. P.; Jones, D. E. Can. J. Chem. 1977, 55, 2977. Handa, Y. P. J. Chem. Thermodyn. 1977, 9, 117. Murray, F. E.; Schneider, W. G.Can. J . Chem. 1955, 33, 797. Remko, M.; Polcin, J. Z . Phys. Chem. (Munich) 1976, 102, 161. (IO) Figeys, H. P.; Geerlings, P.; Berckmans, D.; Van Alsenoy, C. J . Chem. Soc.. Faraday Trans. 2 1981, 77, 721. ( I I ) Berkeley, P. J.: Hanna, M. W. J . Phys. Chem. 1963, 67, 846. (12) Howard, B. B.; Jumper, C. F.; Emerson, M. T. J . Mol. Spectrosc. (6) (7) (8) (9)
1963, 10, 117. (13) Lin, W.-c.; Tsay, S.-j. J. Phys. Chem. 1970, 74, 1037. (14) Lichter, R. L.; Roberts, J. D. J . Phys. Chem. 1970, 74, 912. (15) Woo. K. W.Ph.D. Thesis, University of Detroit, 1976; Diss. Abstr. Int. B. 1977, 37, 3986. (16) Tokuhiro, T. J . Chem. Sot., Faraday Trans. 2 1988, 84, 1793. ( I 7) Moradi-Araghi, A.; Schwartz, M. J . Chem. Phys. 1978, 68, 5548. (18) Kovacs, H.; Kowalewski, J.; Maliniak, A,; Stilbs, P. J . Phys. Chem. 1989, 93, 962.
0022-3654/90/2094-7378$02.50/0
(19) Evans, M. J . Mol. Liq. 1983, 25, 149. (20) Evans, M. W. Adv. Mol. Relax. Interaction processes 1982, 24, 123. Evans, M.W. J . Mol. Liu. 1983.25.21 1. Janik, B.; Sciesibski, J.; Evans, M. W.; Kluk, E.; Grochulsk', T. Chem. Phys. 1982, 70, 77. (21) Dietz, W.; Heinzinger, K. Ber. Bunsen-Ges. Phys. Chem. 1985.89, 968. (22) Bohm, H. J.; Ahlrichs, R. Mol. Phys. 1985, 54, 1261. (23) Bohm, H. J.; McDonald, I. R.; Madden, P. A. Mol. Phys. 1983,49, 347. (24) Bohm, H. J.; Lynden-Bell, R. M.; Madden, P. A,; McDonald, I. R. Mol. Phys. 1984, 51, 761. (25) Lynden-Bell, R. M.; Madden, P. A,; Stott, D. T.; Tough, R. Mol. Phys. 1986, 58, 193. (26) Westlund, P.-0.; Lynden-Bell, R. M. Mol. Phys. 1987, 60, 1189. (27) Westlund, P.-0.; Lynden-Bell, R. M. J . M a p . Reson. 1987, 72, 522. (28) Lynden-Bell, R. M.; Westlund, P.-0. Mol. Phys. 1987, 61, 1541. Westlund, P.-0.; Lynden-Bell, R. M. Chem. Phys. Lett. 1989, 154, 67.
0 1990 American Chemical Society
Liquid Mixtures of Acetonitrile and Chloroform
The Journal of Physical Chemistry, Vol. 94, No. 19, 1990 7379
function for acetonitrile, enabling a 4 times faster simulation performance while essentially reproducing the equilibrium results of BBhm et al.
TABLE It(: Self-Diffusion Coefficients, D mol %
Physical Models and Computational Methods Constant-volume molecular dynamics simulations have been carried out for pure chloroform and for four different chloroform/acetonitrile mixtures, consisting of IO, 22, 50, and 78 mol % of acetonitrile. All the simulated systems contained totally 256 molecules in a cubic box. Temperatures during the simulations were close to room temperature and experimental densities were used for the different binary mixtures and pure chloroform. A summary of the physical states and characteristic simulation parameters is given in Table I. Both molecular species, Le., chloroform and acetonitrile, were kept rigid throughout the simulation and the equations of motions were solved by using the fifth-order Nordsieck-Gear predictorcorrector integration s ~ h e m e . ~ ~The , ~ ' rotational motion was treated by using the second-order quaternion method of Sonnen~chein.~~ Pair potentials used are so-called site-site type and given in standard Lennard-Jones and electrostatics 12-6-1 form. The acetonitrile-acetonitrile potentials, taken from the work of Bohm et are based on atomic interaction parameters taken from studies of other liquids and partial charges derived by a b initio calculations. Further details on the acetonitrile potentials are found in the above-mentioned papers. A chloroform-chloroform potential set was constructed by using standard Lennard-Jones parameters for short-range interactions and placing fractional charges on the atoms. Potential parameters for chloroform used in the simulation are gathered in Table 11. We preferred to use the same type of potential (12-6-1) for both acetonitrile and chloroform in order to facilitate the construction of cross potentials between the mixture components. The potential of Bohm and AhlrichsU is of another type (exp6-1). Our potential is the same as that of Evans20 or Dietz and Heinzinger,21except the charges. The charges used in the present simulation give the chloroform dipole moment a value about 35% higher than the gas-phase value of 1.04 D.33 The reason for the enhanced dipole moment is to include, in a rough way, some of the many-body character present in the condensed state due to interactions between molecules at close distance. Strong inductive effects are to be expected in the present solutions owing to the considerable polarizability of chloroform and highly dipolar character of acetonitrile. Figeys et a1.I0 have reported a dipole moment of 1.61 D for chloroform from a quantum mechanical calculation on a molecule pair in a linear arrangement with the chloroform-hydrogen directed toward the acetonitrile-nitrogen. Also, the potential model of Bohm and Ahlrichs22 for simulating chloroform in chloroform yielded a dipole moment of p = 1.23 D. They determined the charge distribution by ab initio calculations. Presently, the same set of fractional charges was used in all of the simulations, including pure chloroform. The high dipole moment leads to strong electrostatic interactions, which in turn can bring about the strong fluctuations in temperature. The chloroform-acetonitrile potential was constructed by employing Berthelot-Lorentz combination rules. The long-range electrostatic interactions were treated by using Ewald summation techniques. Also, minimum image convention was used together with periodic boundary conditions. The number of steps in all of the simulations was 50 000, while the production was accomplished during the last 30000 steps (Le., 30 ps). All the simulations were performed using a modified version of McMOLDYN package34 on a CONVEX C220 supercomputer.
MD
(29) Jorgensen, W. L.; Briggs, J. M. Mol. Phys. 1988, 63, 547. (30) Nordsieck, A. Math. Comput. 1962, 16, 22. (31) Gear, C. W. In Numerical Initial Value Problems in Ordinary Dif ferential Equations; Prentice-Hall: Princeton, NJ, 1971. (32) Sonnenschein, R. Comput. Phys. 1985, 59, 347. (33) Gray, C. G.; Gubbins, K. E. Theory of Molecular Fluids, Volume I : Fundamentals; Clarendon Press: Oxford, U.K.,1984. (34) Laaksonen, A. Comput. Phys. Commun. 1986, 42, 271.
X
lo9 (m3 s-')
acetonitrile
acetonitrile 0 10 22 50 78 100
chloroform
NMR
3.81 3.85 3.88 4.00 4.3"
2.66 2.60 (15%) 2.88 3.28 (70%) 4.04
MD
NMR
3.38 3.38 3.43 3.37 3.63
2.45 2.40 2.33 (15%) 2.42 2.68 (70%)
"291 K (ref 24).
TABLE I V Rotational Correlation Times,
T
~
(p) , ~
acetonitrile
mol 5% acetonitrile 0
MD
NMR
IO
1.25 1.25 1.20 1.oo 0.95'
1.87 1.83 1.52 1.21 1.02
22 50 78 100
chloroform
MD
NMR
0.96 0.94 0.95 1.02 0.96
1.92 2.09 2.23 2.38 2.09
'291 K (ref 24)
All the time correlation functions were calculated upto 0.50 ps. This was long enough for the translational and rotational velocity correlation functions to decay to zero. The procedure to obtain the reorientational correlation times is discussed in the following section.
Comparison with Experiment To begin with, a comparison of the self-diffusion coefficients and correlation times for the top axis rotation with the experimental data is presented in order to illustrate the extent to which the model liquids are compatible with real solutions. The reasons for certain discrepancies are also discussed. The self-diffusion coefficient is a sensitive indicator of the accuracy of a molecular potential. There are several ways to obtain the diffusion coefficient from simulations; it can be determined, for instance, from the mean-squared displacement of the molecule, or from the integral of the linear velocity autocorrelation function. The simulated self-diffusion coefficients in Table I11 are calculated from the integral of the linear velocity correlation function through the relation D =
1 J m ( a ( t ) O(0)) dt = kLl T 7, 3
0
m
For the sake of comparison, the experimental self-diffusion coefficients recorded by the Fourier transform pulse gradient spin-echo technique in NMR35are also given in Table 111. It can be seen that the experimental concentration dependence is properly reproduced in the computer liquid. As for the absolute values, the simulated diffusion coefficients are consistently 30% larger than the experimental ones. Thus, the agreement is reasonably good, particularly if it is borne in mind that, although the experiments were carried out with deuteriated substances, no correction for the viscosity has been performed. It seems that the diffusion coefficients, as well as other results, are reproduced more accurately for acetonitrile than for chloroform. The explanation for this difference must lie in the simulation model. The correlation times for the reorientation of the symmetry axis, T ~ , were ~ , extracted from the corresponding time correlation functions for second-rank spherical harmonics. This was done in two steps: first by integrating the simulated part of the correlation function up to 0.5 ps, and second by fitting the remainder of the correlation function to an exponential expression of type A exp(B/t) and integrating. The sum of these two integrals is the simulated correlation time, the values of which are compared in Table IV with the experimental values determined by N M R (35) Stilbs, P. Prog. N M R Spectrosc. 1987, 19, I .
7380 The Journal of Physical Chemistry, Vol. 94, No. 19, 1990
Kovacs et al.
TABLE V: Comparison of Simulated Dynamic Dot. of Chloroform dipole temp, D X IO9,' moment, D K m2 s-I 71) PS present work 1.40 288 3.38 I,0.96 I,1.3 Evans (ref 20) 1.04 293 4.00
11, 1.5 Bohm and Ahlrichs (ref 22) Dietz and Heinzinger (ref 21) a
1.23
293
2.0
I , 1.65
II, 1.3 1.1
295
2.6
I ,1.2 11, 1.4
Diffusion coefficient. bRotational correlation time.
TABLE Vi: Interaction Energies (kJ mol-') CHCI3CHCI3mol % acetonitrile CHC13 CH3CN 0 (Coulomb) -455 -71 I7 0 (Lennard-Jones) I O (Coulomb) -402 -316 I O (Lennard-Jones) -764 -5997 22 (Coulomb) -214 -623 22 (Lennard-Jones) -4715 -1328 -787 50 (Coulomb) -6 8 50 (Lennard-Jones) -2276 -2203 -45 1 -14 78 (Coulomb) -530 -1692 78 (Lennard-Jones)
CH3CNCH3CN
-92 -24 -312 -104 -1271 -596 -2583 -1600
relaxation measurements. Here again, the concentration dependence is reproduced accurately, especially for acetonitrile, and the agreement of the numerical values is reasonable. We should recall that the experimental chloroform data do not follow the viscosity, as was stated in the previous N M R study.'* The simulated chloroform T~~ actually varies very little with concentration. A comparison of some dynamic results from the various available simulations of chloroform are presented in Table V. The rotational correlation times appear to be quite insensitive to the different potential parameters, all of the values being around unity. Besides, the chloroform reorientation is, to a very good approximation, isotropic according to all the various models. The diffusion coefficient of Bohm and AhlrichsZ2and Dietz and Heinzinger2' are almost equal to the observed values. Dietz and Heinzinger reported also that the total removal of the partial charges did not affect the translational mobility of chloroform in a detectable way. The present potential model, being constructed to compensate for the polarization of chloroform in acetonitrile, gives somewhat too high molecular mobility in pure chloroform. In conclusion, we find that the agreement between the MD simulations and the NMR experiments is good enough to motivate further analysis of the simulated trajectories. Static Properties We begin the discussion on the time-averaged properties by investigating the interaction energies obtained during the simulation. The Lennard-Jones and Coulomb contributions to the pairwise interaction energies between either like or unlike components in the different mixtures are reported in Table VI. It can be observed that, for the chloroform-chloroform interactions in the simulated liquids, the Coulombic contributions are of minor size, not more than 6% of the Lennard-Jones term. The chloroform potential is thus-also thanks to the fairly globular shape of the molecule-quite isotropic, contrary to the acetonitrile potential, which is strongly dominated by the long-ranging electrostatic contribution. The motional anisotropy of acetonitrile depends partly on the high ratio between its components of the = 9.8, but it is the electrostatics moment of inertia tensor; 11/111 that makes a major difference in the motional behavior of these two molecules. Concerning the acetonitrileacetonitrile interaction energies reported in the third column of Table VI, their Coulombic part is much larger than the Lennard-Jones part. It can be observed that not only do the long-range electrostatic interactions persist in the diluted solution but also they have more than doubled in proportion to the Lennard-Jones part. Thus, there is a qual-
CIIJCN 10%
2'5
1
2.0
-
1.5
-
1.0
-
.-''-. CHsCN
JN-CZ
J
22% 50% 18%
1
4
0.0
1, O
I
,
,
,
I,
2
, ,
--4
6
a
IO
(.i) Figure 1. Some of the a t o m a t o m radial distribution functions between two acetonitrile molecules: (top) g c z a ( R ) .(middle) gNK2(R),(bottom) gN-H(R).
itative correlation with the experimental dielectric constant, which diminishes in this system from 37 of pure acetonitrile to 4 of pure chloroform. Next we shall analyze the radial distribution functions (rdf) and coordination numbers. Throughout the following analysis, the differences and variations in molecular size and shape of the two mixture components have not been corrected for. This is because we regard the two types of molecules as sufficiently similar to neglect such effects; for instance, the molar volume of acetonitrile is 68% of the molar volume of chloroform. Acetonitrile-Acetonitrile. In Figure 1 we present 3 of the 10 different atom-atom pair distribution functions for acetonitrile. (Concerning notations, C1 signifies the methyl carbon and C2
.
2.0 1
'
1
The Journal of Physical Chemistry, Vol. 94, No. 19, 1990 7381
Liquid Mixtures of Acetonitrile and Chloroform
. CEsCN
OX 10%
9c-c
22%
'*' 1.6
1.5
50%
IC-8
7ax 1.0
I
0.5
0.5
0.0 I
0.0
l " " l " " l " " l " " 1
O
2
8
6
4
10
0
2
4
6
8
10
(1) OX
C E s C N OX
10%
CEfCN 1 0 %
'Ci-E
22%
CEfCN 2 2 %
f.ol
50%
I
78%
0.5
O
2
4
6
E
10
1(I)
O
2
4
6
a
IO
(1)
Figure 2. Some of the atom-atom radial distribution functions bctween two chloroform molecules, from top to bottom: (top left) gCx(R), (bottom left) gCI-H(R), (top right) gC-H(R)9 (bottom right) gcl-cdR). on NMR measurements. Previously, we suggested that the the nitrile carbon in acetonitrile; C denotes carbon in chloroform.) structuring of acetonitrile grows upon dilution until 20 mol %, First, it should be noted that the shapes of these functions closely resemble those of Bohm et al.;23in particular, the positions of the whereafter further dilution leads to a breakdown of the ordered arrangement. maxima and minima are identical with the results reported by Bohm. This implies that the acetonitrile orientations in the The integral of the radial distribution function between atoms i and j in different molecules, multiplied by the number density, mixtures are similar to the molecular arrangement in pure acetonitrile. Bohm et al.23suggested that there are about 10 acegives the coordination number, zij, at a certain distance. In Table tonitrile molecules around the central molecule in the first coVIIa are given the coordination numbers as a function of mixture ordination shell. These neighboring molecules are of two kinds, composition derived from the radial distribution functions for the three closest ones being oriented predominantly in a antiparallel chloroform-carbon and nitrile-carbon. These two atoms were chosen because they nearly coincide with the centers of masses. fashion while the orientational correlations for the rest of the The coordination numbers were obtained by integrating to 7.02 molecules are significantly weaker. A, a compromise for the first minima in the three rdfs of current In all of the acetonitrile rdfs, the amplitude increases upon interest, gCZC2(R), dilution, and the increase is particularly sharp in the range of gc2c(R),and gcC(R). In the third column of Table VIIa we have calculated the local molecular fraction of strong orientational correlation close to the central atom. For acetonitrile, zczcz/(zc2c2 + zCz& when the central atom is the instance, the examination of the gcZx2(R)function, which corresponds approximately to the center of mass radial distribution nitrile carbon. The local molecular fraction is actually always function, shows that the shoulder a t about 3.8 A which was ina little higher than the mole fraction, but the difference is not large enough to indicate clustering in the model liquid. terpreted by Bohm as originating from the three closest, antiparallel neighbors grows into a peak upon dilution. Thus the Chloroform-Chloroform. In Figure 2 are plotted 4 of the 6 nearest-neighbor orientational correlation becomes much stronger atom-atom pair correlation functions for chloroform. A quick as the acetonitrile concentration decreases. The same thing is comparison with the rdfs published previously for pure chloroapparent in the gNC2(R)function, and especially in the gN+,(R) forms22 confirms that the contours of the curves are reproduced. Particularly, the simulations by Evans20 and Dietz and Heinfunction, where the first peak has gained very much in amplitude. zingerz1 performed with potentials similar to ours seem to yield That peak corresponds to two parallel neighbors with their nia structure identical with the present result. The simulation trogens directed toward the methyl protons of the central molecule. reported by Bohm and Ahlrichs,22 performed with an a b initio, The orientational correlation seems to increase even for the repure pair potential, leads to somewhat closer distances of the first maining molecules in the first coordination sphere. maxima in the radial distribution functions involving hydrogen. It may be pointed out that none of the peaks has shifted sigTheir distribution of charges differs from the other potential nificantly, which shows that the mutual positions or orientations of the molecules are not altered. The increasing rdf amplitude models by having a large negative partial change at the carbon is a sign of enhanced acetonitrile interactions in diluted solutions, site. A controversial point has been how extensive the structure which agrees with the observations we have made earlier based might be in this fairly nonpolar liquid.
7382 The Journal of Physical Chemistry, Vol. 94, No. 19, I990 TABLE VII: Analysis of Coordination Numbers at 7.02 A for Acetonitrile and Chloroform
Kovacs et al. C B j C N 10% 2.5
,-'...C B j C N
,,
O1-N
a. Acetonitrile , e ' . .
mol % acetonitrile
gC*
