J . Phys. Chem. 1987, 91, 4907-4912 bonded pairs to link up in formation of chains and eventually rings as in the fully coupled model. It will be important to see whether manipulation of X in series of molecular dynamics simulations at low temperature can hasten the otherwise very slow attainment of chemical equilibrium.
VII. Summary and Conclusions Our first attempt1' at modeling liquid sulfur demonstrated that at least the gross features of the structural chemistry of that element could be reproduced with a combination of two-atom and three-atom potentials. The present study reports a refinement in choice of those potentials which improves the fidelity of description and in particular gives cyclic s8 mokcules special stability. No doubt further refinements are still possible, but probably any serious effort of this sort should await appearance of a reasonably complete set of ab initio quantum mechanical calculations for small sulfur clusters.18 In the long run it is desirable to have the model potential reproduce the known crystal structures for sulfur? though the van der Waals forces involved in these structures are quite weak and may thus be largely irrelevant for the chemical reactivity studied in this paper. The available experimental measurements of liquid sulfur pair correlation function^^^*^^ are few in number and too imprecise to
4907
allow instructive comparisons with molecular dynamics simulation results. We hope the present paper and its predecessor" will stimulate more experimental activity that takes advantage of recent advances in instrumentation. While the perturbation proposal outlined in section VI looks promising, it needs to be supported by (a) more simulation studies of g(r,X) at X values between 0 and 1, (b) development of analytical theory for g(r,O) and for fo(T,p),and (c) examination of the way that decreasing X below 1 affects the polymerization anomaly in the liquid, observed experimentally at 159 O C . I 0 Finally, we mention that the present class of interaction models could be used to investigate the possibility of divalent, but topologically distinctive, polymeric species. These include interlinked large rings (sulfur analogues of the catenanes) and knotted rings. Adroit selection of initial atom positions, followed by steepest descent motions on the potential energy hypersurface, should help to estimate the lower size limits for sulfur clusters which could exhibit such unusual structures without undue steric hindrance. Registry No. S , 7704-34-9; Sa, 10544-50-0. (24) TomDon. C. W.: Ginerich. N. S. J . Chem. Phvs. 1959. 31. 1598. (25j Poltivtskv, Yu. G.; Tytenko, Yu. V. Russ. J . j h p . Chkm.'(Engl. Transl.) 1975, 49, 178.
Computer Simulation of Mixtures of Hard Spheres G. Jackson,+ J. S. Rowlinson,* and F. van Swolt Physical Chemistry Laboratory, South Parks Road, Oxford, OX1 3QZ, U.K. (Received: December 3, 1986)
Computer simulationsof binary mixtures of hard spheres have been used to generate their equation of state, diffusion coefficients, and radial distribution functions. The size ratios used were 11/10, 5/3, 3, 5 , and 20. The results confirm the accuracy, on the fluid branch of the isotherm, of the equation of state of Mansoori, Carnahan, Starling, and Leland. Significant deviations were found only at the highest densities and/or for large size ratios. The pressure, diffusion coefficients, and radial distribution functions all tell us something of the conditions under which one or both of the components can solidify.
Introduction The properties of assemblies of hard spheres and mixtures of hard spheres are of interest in their own right and are important ingredients in perturbation theories of mixtures. Computer simulation is difficult at high densities for mixtures of spheres whose ratios of diameters differ greatly from unity. Here we obtain the compression factors, Z = p V / N k T , and the diffusion coefficients of binary mixtures with a ratio of diameters between 5/3 and 20, and at reduced densities E of up to 0.6, where
E
= (.rr/6V)(N1d13+ N2dZ3) (1) for N 1 spheres of diameter d l and N2 of diameter d2. The cross diameter is ( d , + d 2 ) / 2 . The maximum value of 5 for a system of spheres of uniform size occurs in crystalline close packing, for which E = ( ~ / 6 ) ( 2 ) = ' /0.7405. ~ The maximum values in mixtures are not known but are probably a little higher. For a few of the systems we have also obtained the pair distribution functions, gll(r), g12(r), and g z m . We have used the results, first, to test thoroughly the equation of state proposed by Mansoori, Carnahan, Starling, and Leland.' This is a generalization to mixtures of their earlier and well-tested proposal2 that Z ( t ) can be approximated by 9
1
Present address: School of Chemical Engineering, Olin Hall, Cornel1 University, Ithaca, NY 14853-5201.
0022-3654/87/2091-4907$01 S O / O
where Z, and Z, are calculated from the Percus-Yevick equations for g(r) via the compressibility and virial routes, respectively. The PY equations for mixtures of hard spheres were solved by Le-
bow it^.^ Our second purpose is to study the conditions under which these systems solidify. This is a difficult problem to which our results are a partial solution. Earlier simulations have used the Monte Carlo method in the canonical ensemble, which we abbreviate MC(NVT), or the molecular dynamic, MD, method. In the first method we have the results of Smith and Lea4 for (mainly) equimolar mixtures with a size ratio d 2 / d l = 5/3, the results of Rotenbergs for a ratio of 11/ 10, and the recent results of Fries and Hansen6 for ratios of 2 and 3 for mixtures dilute in the larger component. The molecular dynamic results of Alder' are for an equimolar mixture of size ratio 3. It is difficult to generate physically acceptable configurations for the MC(NVT) method if the density is high and the size ratio large, so we have used the isobaric ensemble, MC(NpT), in which (1) Mansoori, G. A,; Carnahan, N. F.; Starling, K. E.; Leland, T. W. J . Chem. Phys. 1971, 54, 1523. (2) Carnahan, N. F.; Starling, K. E. J . Chem. Phys. 1969, 51, 635. (3) Lebowitz, J. L. Phys. Rev. A 1964, 133, 895. (4) Smith, E. B.; Lea, K. R. Trans. Faraday SOC.1963, 59, 1535. (5) Rotenberg, A. J . Chem. Phys. 1965, 43, 4377. (6) Fries, P. H.; Hansen, J. P. Mol. Phys. 1983, 48, 891. (7) Alder, B. J. J . Chem. Phys. 1964, 40, 2724.
0 1987 American Chemical Society
4908
Jackson et ai.
The Journal of Physical Chemistry, Vol. 91, No. 19, 1987
25 '
20-
0 5000
~2
z
0
MCiNpT) MCiNVTi
0
+I
/
d2 i d , = 513
A MD
d21d,=513
2
15-
Z
/.'
10-
5-
/
f
/Q
AZ
0
0
Figure 1. The compression factor Z for a size ratio of 5 / 3 . The MC(NVT) results are those of Smith and Lea. The curve is the equation of state of Mansoori et al.
-1'
I
I
0.5
I
I
0.6
5
Figure 3. The difference, AZ, between the computer value of Z and that calculated from the equation of Mansoori et al. for a size ratio of 5 / 3 . Other conventions as in Figure 1.
Figure 2. The compression factor Z for a size ratio of 5 / 3 . The curve is the equation of state of Mansoori et al.
i OO
the difficulty is less, although there is still a danger of not sampling correctly all regions of configuration space. We therefore check the accuracy of these results by using also M D runs that start from one of the acceptable M C configurations. The choice of a common configuration does not entirely eliminate the sampling problem, but the two kinds of simulations evolve quite differently, and we hope that the good agreement we obtain between the two methods shows that both are accurate. We have also made two MC(NpT) runs with cells that are rectangular prisms of variable shape to see if this extra freedom allows a better sampling or eases the transition to a solid phase but found that there were only negligible changes of a few parts per thousand in the ratios of the lengths of the sides. Results The system comprised 108 molecules in a cubic box with the , usual repeating boundary conditions. The size ratios, d 2 / d l were 5/3, 3, 5, and 20. The number of molecules of species 2, the larger sphere, was 11, 54, or 97; that is, x2 = 0.1019, 0.5, or 0.8981. A few runs were made with seven molecules (x2 = 0.0648) for comparison with the results of Fries and Hansen. In the MC(NpT) runs the system was equilibrated at each pressure for (0.5-1.5) X lo6 configurations, depending on the density, and averages taken over runs of (1 5 2 . 0 ) X lo6 further configurations.
4
t
8
12
Figure 4. The mean-square displacements of the spheres in the system of size ratio 5 / 3 , x2 = 0.5000, [ = 0.5400. The time is in units of d,(m/kT)'f2.
The volume was changed randomly every 100 configurations. The MD runs started from a M C configuration, were equilibrated for 5 X lo4 collisions, and averages taken over runs of (1.5-3.0) X lo5 further collisions. A few longer runs were made to get good values of g(r). Size Ratio o f 5 / 3 . Figure 1 shows the values of 2 in the equimolar mixture from the MC(NVT) results of Smith and Lea at low densities, and the new values obtained by MC(NpT) and by MD at high densities. All are in good agreement with each other and with the equation of state of Mansoori et al. There is no sign of a transition to a solid. The agreement is equally good at x2 = 0.1019 and on the fluid branch at x2 = 0.8981 (Figure 2). This last system does, however, form a solid phase at densities above 0.52, as is shown by the sudden fall in the pressure. Figure 3 shows in more detail the departures from the equation of Mansoori et al. These are small but positive at x2 = 0.1019 and 0.5000, except at the highest density where the deviations increase sharply. The deviations at x2 = 0.8981 seem to be small but negative on the fluid branch and large and negative on the solid. The precision of the MD results is a little higher than those of
The Journal of Physical Chemistry, Vol. 91, No. 19, 1987 4909
Mixtures of Hard Spheres
1.0
d,ld,=513
p0'04 0.0'
I
I
I
I
- 0.02 2
;'30:
0,ZO-
0.00,
4
2
6
8
t
0.10
Figure 5. The mean-square displacements of the spheres in the system of size ratio 5/3, x2 = 0.8981, 6 = 0.6000. 0
O.OO
d, /d,=5/3
0.30
0.00 O ..
.
l
0
o=ho-&-
F
1
\o \
0.00'0145
\O
0150
J
0,60
5
Figure 6. The self-diffusioncoefficients,D,and D,,in units of d , ( k T / m)'/*in the system of size ratio 5/3, x2 = 0.1019.
I
I -
00
I
L
O
An.
04
r:
0.55
0.60
SI
Figure 7. The self-diffusion coefficients in the system of size ratio 5/3, x 2 = 0.5000. Conventions as in Figure 6 .
the MC, and both are poor for the solid. Figure 4 shows the mean square displacements for spheres of equal mass in the M D run at f = 0.52, x2 = 0.5; those for x2 = 0.1019 are similar. At x2 = 0.8981 the behavior is different and at densities above 0.52 there is little or no diffusion of either component, as is apparent, for example, in Figure 5 . The derived diffusion coefficients are shown in Figures 6-8. Size Ratio of 3. The results are similar to those for the size ratio of 5/3. In the equimolar mixture (Figure 9) the MC(NpT)
Figure 9. The compression factor for a system with a size ratio of 3, x 2 = 0.5000. Conventions as in Figure 1, but the MD results are those of
Alder. results agree well with Alder's MD results, and both agree with the equation of Mansoori et al. At x2 = 0.0648 there is similar agreement with the MC(NVT) results of Fries and Hansen, and with the equation. At x2 = 0.8981 there is again a phase transition at a density of about 0.53, but now the diffusion coefficients (Figure 10) show that it is only the larger component that has crystallized; the smaller continues to move through the lattice. Size Ratio of 5. The deviations from the equation of Mansoori et al. are again small in the fluid state (Figure 1l ) , but there is now apparently a solid branch also in the equimolar mixture, although the evidence from the diffusion coefficients (Figure 12) is less clear-cut. It seems (Figure 11) as if the solidification is a continuous change, that is, it may not be a first-order transition, in this simulation. Size Ratio of20. Figure 13 shows that the positive deviations from the equation of Mansoori et al. are now strong at x2 = 0.1019; they are clearly visible also for smaller size ratios in Figure 3 and Figure 11. Again there is a sharp transition to the solid at x2 = 0.8981, and a smoother transition at x2 = 0.5. However, the diffusion coefficients show that there is virtually no motion of the larger component at either mole fraction beyond a density of 0.52 or 0.54 (Figures 14 and 15). Radial distribution functions were calculated for x2 = 0.5, 5 = 0.45,0.5681 and x2 = 0.8981, E = 0.45,0.5651. The two lower
4910 The Journal of Physical Chemistry, Vol. 91, No. 19, 1987
Jackson et ai.
-
0'0 ' 0
2-
-
-
0 '
I
O\O' - ' .
100.5 -
\ O
b
Yo :
0
Con-
Figure 10.
ventions as in Figure 6 .
-1
Figure 12. The self-diffusion coefficients in the system of size ratio 5 , x2 = 0.8981. Conventions as in Figure 6 .
-I t
-11
x,=O 0981
-2
-t+f
A
x2
= 0'8981
I
-3
-4
+
AZ
-4)
-11
I
I
I
0.5
I
5 Figure 11. The difference, A Z , between the computer values of Z and that of the equation of Mansoori et al. for a size ratio of 5 . Conventions as in Figure 1.
densities show typical fluid curves of which Figure 16 is an ex) ample, while the two higher densities have values of g Z 2 ( r that are small compared with unity at certain separations (Figures 17 and 18) which is what might be expected in a solid. Again the evidence for solidification is not so strong at x2 = 0.5 as for x 2 = 0.8981 where g22(r)shows a particular strong peak at 1.85d2. The next nearest neighbor in a close-packed crystal lattice of spheres of uniform size is at a reduced separation of 4(2)'j2/3 = 1.886. Discussion The primary purpose of this paper was to obtain good values of the pressure and density on the fluid branches of the isotherms, and to use them to test the equation of state of Mansoori et al. The results show that the equation behaves well except at the
Figure 13. The difference, A Z , between the computer values of Z and that of the equation of Mansoori et al. for a size ratio of 20.
highest densities ( E > O S ) , and that the deviations that occur are worst for the largest size ratio of 20. The accuracy of our results appears, from the agreement of the MD and M C runs, to be about il% on the fluid branch, and the agreement with the equation of state is similar. The failure of the equation at very large size ratios was foreshadowed by the warning of Waisman, Henderson, and Lebowitz* from their studies of a system of infinite size ratio, that is, for hard spheres against a hard planar wall. Our secondary purpose was to study the solidification of these systems. A complete study of this problem is difficult since we do not know what crystal lattices can be formed by binary mixtures of arbitrary size ratio and cannot therefore choose in advance a cell that can become a multiple of the unit cell of the lattice. (8) 1373.
Waisman, E.; Henderson, D.; Lebowitz, J. L. Mol. Phys.
1976, 32.
The Journal of Physical Chemistry, Vol. 91, No. 19, 1987 491 1
Mixtures of Hard Spheres
I
I
I
I
I
I
8L
922(r)
4 5
Figure 14. The self-diffusion coefficients in the system of size ratio 20, x2 = 0.5000. Conventions as in Figure 6.
10
OO
rid,
30
20
Figure 17. The three radial distribution functions in the system of size ratio 20, x2 = 0.5000, f = 0.5681. This is probably a solid state. I
\O:::
20
-
0 I
0
-
Figure 15. The self-diffusioncoefficients in the system of size ratio 20,
x2 = 0.8981. Conventions as in Figure 6.
10
0'
20
30
40
rldl d21d,=20
S = O 4500
Figure 18. The three radial distribution functions in the system of size ratio 20, x2 = 0.8981, f = 0.5651. This is a solid state. as proposed by Parrinello and R a h ~ n a n . The ~ density-functional theory of Barrat, Baus, and Hansen'O suggests that even in the limited range of size ratios of between 1 and 1.18 there are complicated phase diagrams to be found. Our results suggest crystallization occurs, if at all, at a reduced density of about 0.50-0.53, which is a little higher than the value of 0.49 found for a one-component system. Systems with a small mole fraction of large spheres, x2 = 0.0648 and 0.1019, did not solidify, those with x2 = 0.5 solidified only if the size ratio was close to unity (11/10) or large ( 5 and 20). The intermediate range, 5/3 and 3, remained fluid. The systems with few small spheres, x 2 = 0.8981, solidified readily, but for a size ratio of 3 and larger the smaller component moved through the fixed array of larger spheres. If the larger spheres form a close-packed fcc structure then there are octahedral holes that can accommodate spheres smaller by a factor of 1 2lI2 = 2.41, and tetrahedral holes that can only accommodate spheres smaller by a factor of 2 6lI2 = 4.45. Our structures were not close packed and it is therefore reasonable to find that spheres smaller by a factor of 3 or more can diffuse through what was certainly a very imperfect crystal. Ermak, Alder, and Prattl' found that a size ratio no smaller than
4L
922 (r)
2
'0
10
30
20
40
rldl Figure 16. The three radial distribution functions in the system of size ratio 20, x2 = 0.8981, f = 0.4500. This is a fluid state. In this and the next two figures the lines are the results of the simulations, and are not related to the calculations of Mansoori et al. Ideally we should allow for random changes of cell lengths and of the angles between the faces, while keeping the stress constant,
+
+
(9) Parrinello, M.; Rahman, A. Phys. Reu. Lett. 1980, 45, 1196. (10) Barrat, J. L.; Baus, M.; Hansen, J. P. Phys. Reu. Lett. 1986,56, 1063.
J. Phys. Chem. 1987, 91, 4912-4916
4912
2.5 can be tolerated in a fcc lattice without either both components becoming localized or both fluid. We do not, however, know anything of the crystal structures formed in these runs, nor how they packed into the cells available to them. A full study of the (11) Ermak, D. L.; Alder, B. J.; Pratt, L. R. J . Phys. Chem. 1981, 85, 3221.
problem of the crystallization of mixtures of hard spheres has yet to be made.
Acknowledgment. We thank Professor B. J. Alder for useful discussions. G.J. thanks the SERC for the award of a CASE studentship, and F.V.S. thanks the Ramsay Trustees for a Fellowship.
Onsager’s Spherocylinders Revisited Daan Frenkel FOM Institute for Atomic and Molecular Physics, 1009 DB Amsterdam, The Netherlands (Received: December 4, 1986; In Final Form: February 9, 1987)
Since the work of Onsager, systems of hard spherocylinders have played a special role as the theoretician’s “ideal” nematic liquid crystals. It is, however, well-known that the Onsager theory is only a good description for extremely nonspherical particles. A quantitative measure for the range of validity of this theory is obtained from direct numerical calculation of the third through fifth virial coefficients of hard spherocylinders with length-to-width ratios, LID, between 1 and lo6. The phase diagram of spherocylinderswith “realistic” LID ratios can only be obtained by computer simulation. We report molecular dynamics and Monte Carlo simulations for spherocylinders with LID = 5 . These calculations suggest that the oldest model for a hard-core nematic may, in fact, have a smectic liquid-crygalline phase as well.
Introduction Liquid crystals are probably good candidates for the title ”nonsimple liquids”. They are complex both in the microscopic nature of their building blocks and in the diversity of partially ordered phases that may be found in nature. The implicit hope of theoreticians who propose models for liquid-crystalline mesophases is that only one or two aspects of the microscopic structure of mesogens (liquid-crystal-forming molecules) are essential for the understanding of the liquid-crystalline phases. Unfortunately, it appears that no single structural property can account for the wide diversity of known mesophases. What is worse, there is not even consensus about the choice of this “essential feature” for any given phase. Onsager was the first to propose that orientational ordering in a molecular fluid may be explained as an excluded-volume effect.’ The model considered by Onsager was an assembly of thin spherocylinders with length L and diameter D. Onsager showed that, in the limit L I D m, this model system undergoes a transition from the isotropic fluid phase to an orientationally ordered (nematic) liquid-crystalline phase. The physical reason for this phase transition is that, by forming a nematic, the system can gain translational entropy at the expense of some orientational entropy. The transition takes place at a density p* of order 1/(L2/D) (p*B2 = 3.29,2 where the second virial coefficient B, = rL2D/4). At the transition, the volume fraction occupied by the spherocylinders is still vanishingly small (of order D / t ) . Subsequently, alternative physical mechanisms have been invoked to describe the isotropic-nematic transition. For instance, the mean-field theory of Maier and Saupe3 attributes the orientational ordering to anisotropic intermolecular dispersion forces. Among experimentalists, the Maier-Saupe theory enjoys much greater popularity than the Onsager description. The reason is simple: the Onsager theory yields rather poor predictions for the properties of real thermotropic nematics (such as the value of the nematic order parameter S and the density change Ap at the
-
(1) Onsager, L. Proc. N . Y . Acad. Sci. 1949, 51, 627. (2) Lasher, G. J. Chem. Phys. 1970, 53, 4141. Kayser, R. F.; Ravechi, H.J. Phys. Rev. A 1978, 17, 2067. (3) Maier, W.; Saupe, A. Z . Naturforsch., A 1958, 13, 564.
0022-3654/87/2091-4912SOl S O / O
isotropic-nematic t r a n ~ i t i o n . ~Still, Onsager’s model plays a unique role in the theory of liquid crystals, because it is the only exactly solvable model with full translational and orientational degrees of freedom which exhibits a transition to the nematic phase. Virial Coefficients Unfortunately, the Onsager theory is only valid in the limit LID m, while most thermotropic liquid crystal have effective L I D ratios of 3-5.5 The reason why the Onsager theory cannot be used to describe molecular systems with such “small” LID values is the following: an essential assumption in its derivation is that, in the expansion of the free energy in powers of the density, all virial coefficients B, with n > 2 may be neglected. This condition is satisfied if the reduced nth virial coefficient, defined as
-
B,(reduced)
Bn/B2(”-
(1)
I)
satisfies the condition B,(reduced) > 1.
0 1987 American Chemical Society
