J . Phys. Chem. 1988, 92, 3280-3284
3280
arise naturally in a variety of physical contexts, including the numerical path integral study of quantum dynamics. The development of the present approach and the related method, the coordinate rotation t e c h n i q ~ e , suggests ~.~ that progress is being made toward the development of a general Monte Carlo theory of quantum dynamics. It is likely that the methods described here will find use in other problem areas where there are difficulties associated with phase oscillations. It is important to recognize that there are a number of potential pitfalls associated with applications of any of the approximate forms of the damping function. Perhaps the most serious is the fact that the results of integrations using an approximate damping function are not exact. There is an additional difficulty in applications of any of the second-order damping functions [eq 13 or 151. For multidimensional integrations the second derivatives required in the evaluation of the second-order damping functions take the form of determinants. The generation of such deter-
minants for many-dimensional systems may prove to be difficult. A more useful procedure may be one requiring only first derivatives which can be corrected to the exact result. Such a procedure is given in ref 12. Acknowledgment. We thank J. 0.Hirschfelder, R. Wyatt, and R. D. Coalson for organizing the Lasers, Molecules and Methods Conference last summer at Los Alamos. That conference served as an initial sounding board for many of the ideas presented here. The timely quantum dynamics meetings organized this spring at Courant by M. Kalos and J. Lebowtiz and this summer at CECAM by M. Gillan and P. Madden were likewise of substantial assistance. Work at the University of Rhode Island was supported in part by grants from Research Corporation, the donors of the Petroleum Research Fund, administered by the American Chemical Society, and the University of Rhode Island Academic Computer Center.
Structure of Hard-core Models for Liquid Crystals Daan Frenkel FOM Institute for Atomic and Molecular Physics, P.O. Box 41883, 1009 DB Amsterdam, The Netherlands (Received: July 16, 1987; In Final Form: October 28, 1987)
The results of recent computer simulations on fluids of nonspherical hard-core particles are discussed. New data are presented on the structure and dynamics of a system of hard spherocylinders with length-to-width ratio LID = 5 . These data show that such spherocylinderscan occur in at least four stable phases, viz., isotropic fluid, nematic liquid crystal, crystalline solid, and, surprisingly, a smectic A phase.
Introduction Computer simulations of classical many-body systems can be used to gain insight in the microscopic behavior of real liquids and solids. Two distinct, and often complementary, approaches may be distinguished. On the one hand, one may carry out simulations on a realistic models in order to assist the interpretation of real experiments. On the other hand, computer simulations on idealized models of dense phases may be used to test theoretical concepts. Occasionally, computer simulations on simple model systems have yielded results that were qualitatively different from what was expected on the basis of the theories current at the time. A prime example of such a "computer discovery" is the observation by Alder and Wainwright of the long-time tails in the velocityautocorrelation function of hard spheres. Our current theoretical understanding of simple liquids is, to a large extent, based on the results of such simulations. The situation is different for more complex fluids, such as liquid crystals. For the latter class of materials, direct comparison between experiment and simulation is difficult, because such simulations are very time-consuming (although not impossible; see ref 2). Simulation of idealized models for liquid crystals is also less than straightforward because there is no consensus as to what constitutes an "ideal" liquid crystal. For atomic liquids it is well-known that the structure of the fluid is almost completely determined by the short-range repulsive forces acting between the atoms. In fact, the success of the hard-sphere fluid as a reference system in thermodynamic perturbation theories for simple liquids3v4is largely a consequence of this fact. In contrast, it is at present not known whether the structure of more complex liquids, such as liquid crystals, is also (1) Alder, B. J.; Wainwright, T. E. Phys. Rev. A 1970, 1, 18-21. (2) Picken, S.J. Internal Report; University of Groningen, 1984. (3) Barker, J. A.; Henderson, D. J . Chem. Phys. 1967, 48, 4714. (4) Weeks, J . D.; Chandler, D.; Andersen, H. C. J . Chem. Phys. 1971, 54, 5237.
0022-3654/88/2092-3280$01.50/0
primarily determined by excluded volume effects. From the theoretical work of Onsager5 we know that a fluid of (infinitely) thin spherocylinders with length L and diameter D must undergo a transition from the isotropic to the nematic phase at a number density of order l/(L*D). At this density the fraction of the volume occupied by the spherocylinders is still vanishingly small (of order D/L). Recent computer simulations on hard ellipsoids-of-revolution with more realistic shapes6!' have shown that a stable nematic phase is possible for this class of hard-core molecules if the length-to-breadth ratio is either larger than 2.5 or less than 0.4.8 These results do not yet imply that nonspherical hard-core interactions are the cause of orientational order in real nematic liquid crystals. In fact, two additional factors are often invoked to explain the stability of nematic liquid crystals, namely, (1) long-range anisotropic forcesg which tend to induce orientational order and (2) the presence of flexible tails attached to the rigid molecular core.1° The effect of the flexible tails is to stabilize the liquid phase with respect to the crystalline solid. However, now that we know that hard spheroids, that have neither longrange interactions nor flexible tails, form a nematic phase over a rather wide range of length-to-breadth ratios, we can begin to test thermodynamic perturbation theories. Such tests should enable us to ascertain whether nonspherical hard-core fluids can serve as "reference" systems for nematic liquid crystals in the same way that the hard-sphere fluid is a reference system for, say, liquid argon. Onsager, L. Ann. N . Y. Acad. Sci. 1949, 51, 627. Eppenga, R; Frenkel, D. Mol. Phys. 1984, 52, 1303. Frenkel, D; Mulder, B. M. Mol. Phys. 1985, 55, 117 1. Pioneering simulations of hard spherocylinders with LID = 2 were carried out by Vieillard-Baron Mol. Phys. 1974, 28, 809. However, in that system Vieillard-Baron did not observe nematic ordering. (9) Maier, W.; Saupe, A. Z. Naturforsch., A: Astrophys., Phys., Phys. Chem. 1958, 13A, 564. (10) Dowell, F.; Martire, D. E. J. Chem. Phys. 1978,68, 1088. Ibid. 1978, 68, 1094. Ibid. 1978, 68, 2322. Dowell, F. Phys. Rev. A 1983, 28, 3526. (5) (6) (7) (8)
0 1988 American Chemical Society
The Journal of Physical Chemistry, Vol. 92, No. 11, 1988 3281
Structure of Hard-core Models for Liquid Crystals For all liquid-crystalline phases other than the nematic, the role of a hard-core reference fluid is even more uncertain for the very simple reason that is is far from obvious that such phases can be reproduced at all by a purely repulsive hard-core model. In fact, the possibility is not even discussed in the existing textbooks on liquid crystals (see, e.g., ref 11). It was therefore a great surprise (see, however, ref 12) when Stroobants et al.I3 found evidence for stable smectic phases in computer simulations on systems of parallel spherocylinders with L (length) over D (diameter) ratios between 0.5 and 5. This surprise was heightened when subsequent simulations14 revealed the existence of stable columnar phases in fluids of parallel spherocylinders with LID > 4. However, the fact that the molecular orientation of parallel spherocylinders is fixed makes this model less suitable as reference system for real smectics. It is clearly important to establish whether or not the smectic order is destroyed if the molecular orientations are left free. This paper describes a computer study of hard spherocylinders with a length-to-width ratio LID = 5 and full translational and orientational freedom. I shall present new results on the structure and phase transitions of this model system. These results clearly demonstrate that such hard spherocylinders can undergo spontaneous ordering as the density is increased, first to a nematic phase and next to a smectic liquid crystal. This smectic phase is stable in a density range where the crystalline solid melts spontaneously. The remainder of this paper is organized as follows. First I briefly describe the techniques used to carry out the simulations. For more details the reader is referred to ref 15. Next I will present data on the equation of state and the first five virial coefficients for spherocylinders. Evidence for the isotropic-nematic transition is derived from the behavior of the orientation-dependent radial distribution function. In addition I show evidence for pronounced pretransitional fluctuations in the vicinity of this transition. The nematic-smectic transition is studied by using similar tools. In addition, snapshots of the molecular arrangements are reproduced. These pictures illustrate the nature of the molecular ordering more directly than any correlation function. Computer Simulations
The starting configuration for the computer simulations was a crystalline lattice of parallel spherocylinders in a almost cubic box. The initial configuration was generated starting from an fcc crystal of hard spheres. The [ 11 11 direction was chosen parallel to the z axis. After expansion of the distance between the close-packed [ 1111 planes by a factor 1 + L / D , the spheres are replaced by spherocylinders with length-to-width ratio LID. It is clear that this procedure will, in general, not yield a cubic box. However, by a judicious choice of the number of crystal unit cells in the x, y , and z directions, the box can be made very nearly cubic. The number of spherocylinders in a periodic box was chosen to be 576. This number of particles was needed to guarantee that a given particle would only interact with the nearest periodic image of any other particle. This requirement implies that the linear dimensions of the periodic box in any direction are at least 2(L D ) . The initial configuration was instantaneously expanded to a 20% of the density of regular close packing. Subsequently, a constant-volume Monte Carlo simulation was carried out during which the original lattice rapidly melted while the distribution of molecular orientations became isotropic. After that, the system was compressed stepby-step to higher densities, with density steps of 5% of the density at regular close packing. At each density the system was equilibrated again. At low densities (35% of close
+
~~
~~
~
~~
~
~~
(1 1) de Gennes, P. G. Physics of Liquid Crystals;Oxford University Press:
New York, 1974. (12) Hosino, M.;Nakano, H.; Kimura, H. J . Phys. SOC.Jpn. 1979,46, 1709. (13) Stroobants, A.; Lekkerkerker, H. N. W.; Frenkel, D. Phys. Rev.Lett i9a6,57,1452. (14)Stroobants, A,; Lekkerkerker, H. N. W.; Frenkel, D. Phys. Rev.A 19a7.36.2929. (15) Allen, M. P.; Frenkel, D.; Talbot, J., manuscript in preparation.
1 L
45
L/D
=
5.00
P Figure 1. Equation of state of hard spherocylinders with length-to-width ratio L I D = 5. The density is expressed in units of the density of regular close packing. Open circles indicate the MD results for the fluid branch and open triangles for the solid branch. The statistical error in the data points is typically a few tenths of a percent and is therefore less than the size of the symbols. The filled circles indicate the prediction of the five-term virial series computed with the data in Table I. The dashed lines are fits to the MD data points. They have been drawn as a guide to the eye.
packing, and less), equilibration was always achieved in 20 000 trial moves per particle. For the higher densities, where spontaneous ordering occurred, equilibration was achieved by alternating molecular dynamics simulations and (constant stress) M C simulations. After equilibration, M C and M D production runs were carried out. At a density above 45% of regular close packing the system transformed into an orientationally ordered fluid. Once the stability of this phase had been established, further compression was carried out on an orientationally ordered fluid that had previously been aligned parallel to one of the edges of the periodic box. The duration of the M C runs was typically 20000 trial moves per particle. The length of the MD run varied from 2 X lo5 collisions at the lowest density, to almost 2 X lo6 collisions at higher densities. In order to speed up the simulations, the M C and MD programs were vectorized to run on a CYBER 205 supercomputer. The vectorized M D program ran at a speed of some 4 X lo5 collisions/h of CPU time. At densities beyond 60% of close packing, the compression procedure for the 576 particle system became too time-consuming and was abandoned. Instead, a solid of 144 spherocylinders was slowly expanded from 90% of close packing. It melted irreversibly, and with appreciable hysteresis, at 70% of close packing. The reason why a system of 144 rather than 576 particles could be used at high densities is that, after equilibration, the system remains strongly aligned parallel to the longest edge of the periodic box. As a consequence, the probability of an overlap with a next nearest periodic image becomes negligible. Results
The equation of state obtained from the M D simulations is shown in Figure 1 . The pressure of P* is expressed in units kT/vo, where uo is the molecular volume of a spherocylinder: u,, = (7rD3/6)(1 + ( 3 / 2 ) L / D ) . The volume per spherocylinder at regular close packing is uq = (03/2'/*)(1 + ( 3 / 2 ) 1 / 2 L / D )Hence, . for spherocylinders with LID = 5, the packing fraction at close packing equals 88.35%, which should be compared with the corresponding number for spheres (74.05%) on the one hand and for cylinders (90.69%) (=2-D hard disks) on the other. Figure 1 clearly shows the density gap between the fluid branch and the solid branch. In the same figure the virial equation of state is also indicated, using the first five virial coefficients for this system. The virial coefficients were computed by using the method of Ree
3282 The Journal of Physical Chemistry, Vol. 92, No. 11, 1988
Frenkel
TABLE I: Vinal Coefficients for Spherocytinders with L I D = 5‘ B,lB,2 0.4346 + 0.0005 B,IB,~ 0.0669 + 0.0005 0.01 15
+ 0.0008
“The higher order virial coefficients have been expressed in units of B2. The statistical error was estimated from the variance of block averages in the Monte Carlo sampling method of ref 15. The number of independent trial configurations was 5 X lo5. TABLE II: Equation of State of Hard Spherocylinders with L I D = 5’ (PVlnkT) - 1 P* = PIP,, fluid solid N 2.67 0.200 576 5.33 576 0.300 576 7.18 0.350 9.26 0.400 576 1 1.29 576 0.450 12.45 576 0.500 576 13.76 0.569 14.09 0.604 576 576 15.89 0.650 18.11 144 0.700 144 18.49 0.720 144 20.18 0.740 22.62 144 0.759 144 24.63 0.770 144 26.94 0.778 144 33.64 0.800 144 0.820 41.15 15.40 144 0.700 20.00 144 0.760 21.11 144 0.780 23.20 144 0.800 26.08 144 0.820 33.47 144 0.860 47.49 144 0.900
1 -
0.75 -
0.50
-
0.25
o
2
4
6
a
io
Figure 2. Radial distribution function g ( r ) (r = distance between the centers of mass) for hard spherocylinders with L I D = 5. Drawn line: 20% of close packing. Dashed line: 30%of cp. Dotted line: 35% of cp. Long dash-short dash: 40% of cp. Long dash-dot: 45% of cp.
‘The pressure was obtained during constant-volume molecular dynamics runs. The statistical error in the data points is estimated to be less than 1%. and HooverI6 and have been collected in Table I. It should be noted that the virial coefficients for the same model system have been computed earlier by Monson and Rigby.” However, the present values for B3 through BS are slightly, but systematically, larger than the data of ref 17. This same trend is observed when we try to reproduce the virial coefficients of spherocylinders with lower LID ratios ( L I D = 1 and 2) reported in ref 17. However, our data for L / D = 1.5 are in perfect agreement with the corresponding results of Nezbeda (see ref 18). As can be seen in Figure 1 , the five-term virial equation of state accurately reproduces the molecular dynamics data up to a reduced density of 0.4. From that moment on the M D pressure drops below the value computed with the virial series. It is interesting to look at the changes in the local structure of the spherccylinder fluid as the density is increased. One obvious measure of the degree of local order is the radial distribution function g(r). This quantity is shown in Figure 2. At low densities (20% of close packing), g(r) is almost featureless: it increases smoothly from 0 at r = D to 1 for r > 4 0 . As density is increased to 45% of close packing, a peak develops at short distance ( r = 1.2 D ) . But otherwise the correlation function remains featureless. More interesting is the behavior of the distribution function that measures the decay of angular correlation: g2(r) (P2(cose ( r ) ) ) , where e ( r ) is the angle between the axes of two molecules separated by a distance r. The density dependence of g2(r)is shown in Figure 3. Clearly, the local orientational order increases rapidly with increasing density. More importantly, the range of angular (16) Ree, F. H.; Hoover, W. G.J . Chem. Phys. 1964, 40, 939. (17) Monson, P. A.; Rigby, M . Mol Phys. 1978, 35, 1337. (18) Boublik, T.;Nezbeda, I. Collecr. Czech. Chem. Commun. 1986, 51, 2301.
Figure 3. Density dependence of the orientational correlation function g2(r)= ( P 2 (cos e ( r ) ) )(see text), for hard spherocylinders with L I D = 5 . Meaning of symbols as in Figure 2.
correlations grows. In fact, at 45% of close packing, gz(r) does not decay to zero within half the diameter of the periodic box. This suggests that the system at this density is just at the transition point to a nematic phase. Further evidence for the occurrence of spontaneous orientational ordering is provided by the behavior of the nematic order parameter S, defined as -2 times the middle eigenvalue of the tensor Q:6
where u,1 is the arth component of the unit vector along the symmetry axis of molecule i. At lower densities S is always, within the statistical accuracy of the simulations, equal to zero. At 45% of close packing S fluctuates strongly (and slowly; see ref 19): its average value is equal to 0.3. This is a typical value for a nematic order parameter at the transition to the isotropic phase. As the density is increased even more (Figure 4), g2(r)becomes truly long ranged. In the nematic phase, g2(r) S2, as r m. The increase in nematic order is reflected in the behavior of the nematic order parameter S, which grows from 0.3 to more than 0.9, as the density is increased to 60.4% of close packing. Actually, the point at 45% of close packing is possibly still just isotropic: this can be seen from the behavior of the correlation function g4(r) (P4(cos 9 ( r ) ) ) . Figure 5 shows that at 45% of close packing, unlike g2(r),g4(r) still decays to zero within half a box length. However, at higher densities g4(r)shows the same long-ranged
-
(19) Frenkel, D. J . Phys. Chem. 1987, 91, 4912
-
The Journal of Physical Chemistry, Vol. 92, No. 11, 1988 3283
Structure of Hard-core Models for Liquid Crystals 1
0.75 30
n
k
W N
0.50
M
20 0.25
u
0
2
4
6
8
1
0.00
0
Figure 4. Density dependence of the orientational correlation function g2(r)2 ( P 2 (cose(r)))at densities above the isotropic-nematic transition. Drawn line: 50% of cp. Dashed line: 56.9% of cp. Dotted line: 60.4% of cp.
l
i
-
v
0.25
l
0.50
0.75
Figure 6. Density dependence of the static orientational correlation function 1 + gze C(P2(cos e,)) (closed triangles) and of the amplitude of the first peak in the structure factor, S(k,), which is a measure for the amplitude of smectic precursor fluctuations (open triangles). For both 1 g2 and S(kmax),there are higher density points which are off-scale. The drawn lines indicate the steepness of the divergence.
+
2 0.75 n
L
Y
0.50
M
n
p:
, 1
M
0.25
Figure 5. Density dependence of the orientational correlation function g4(r) (P4(cos e@))) (see text), for hard spherocylinders with L I D = 5. Meaning of symbols as in Figure 2.
behavior as g2(r). Incidentally, this is the first direct evidence from computer simulation for nematic ordering in a fluid of hard spherocylinders. The fact that spherocylinders form a nematic phase is gratifying, but it is not the central point of this work. The question we set out to answer is whether the smectic ordering observed in parallel spherocylinders can still occur in a spherocylinder fluid with full orientational freedom. In a previous publication18 I have already indicated that there appears to be a critical increase in the amplitude and decay time of smectic precursor fluctuations as the density of the spherocylinder fluid is increased to 60% of close packing. More direct evidence for the growth of smectic ordering is presented in Figure 6 . This figure shows the growth in the first peak of the structure factor for k values along the direction of the nematic director. The value of the structure factor at a wave vector k measures the mean-squared amplitude of density fluctuations with a wavelength 27rlk. Hence the onset of smectic order should cause a peak in the structure factor at a value of k = (2?r/d)n,where d is the smectic-layer spacing and n the nematic director. The peak of the structure factor clearly diverges a t a density around 60% of close packing. Also indicated in the same figure is the static orientational correlation factor 1 g2 (see, e.g., ref 20). g2 is the "quadrupolar" equivalent of the well-known Kirkwood g factor. g2 is defined through the relation
+
(20) For a discussion, see Kivelson, D.; Madden, P.A. Annu. Reu. Phys. Chem. 1980, 31, 523.
0
4
8
12
16
20
r/D Figure 7. Behavior of the longitudinal density correlation function gll(r), below (drawn line, 50% of close packing) and above (dashed line, 60.4% of close packing) the nematic-smectic transition.
1.25 1 n
k
3 M
0.75
0.25 0.50
I
0
2
4
6
8
1
0
r/D Figure 8. Density dependence of the transverse density correlation function g J r ) (see text), for hard spherocylinders with L I D = 5 , at densities above the isotropic-nematic transition. Drawn line: 50% of cp. Dashed line: 56.9% of cp. Dotted line: 60.4% of cp.
and it measures the degree of orientational correlation in a system. It is of order 1 in an isotropic fluid and of order N (number of particles) in a nematic liquid crystal. Figure 6 shows that g2 diverges at a density well below the point where the peak i n t h e structure factor diverges. The density at which g2 appears to
3284
Frenkel
The Journal of Physical Chemistry, Vol. 92, No. 11, 1988 0.10
0.08
I1
i
0 06 k %Y
tlD
004
1
~
(I II
002
I'
Figure 9. Density dependence of the bond-orientational correlation function g&) (see text), for hard spherocylinders with LID = 5, at densities above the isotropic-nematic transition. Drawn line: 50% of cp. Dashed line: 56.9% of cp. Dotted line: 60.4% of cp.
Figure 11. Snapshot of a configuration of 576 hard spherocylinders at 56.9% of close packing. This density is inside the nematic range
However, smectic precursor fluctuations are already clearly visible Viewpoint as in Figure 10
I
Figure 10. Snapshot of a configuration of 576 hard spherocylinders at 40% of close packing. For the sake of cIarity the spherocylindersare
indicated by a line of length L along the molecular axis. The view shows a projection of all molecules on the x-z plane of the periodic box. At this density the system is in the isotropic phase. diverge should be very close to the isotropic-nematic transition. A direct measure for the degree of smectic order is given by the longitudinal density correlation function g l l ( r ) . This function measures static density correlations along the nematic director. Figure 7 shows that this function is almost structureless at 50% of close packing, while at 60% it is strongly modulated. Evidence that the observed phase is not a solid comes from the transverse correlation function gxy(r). As can be seen in Figure 8, this function develops some short-range structure with increasing density, but solid order is not observed. This is shown more convincingly by the bond-order correlation function g6(r),which measures transverse bond-order correlations. The local in-plane bond order around a given molecule q6(ri)is defined by q6(ri) C w ( r i j )exp(i68,) (3) J
where eijis the angle between the xy projection of rij and the x axis. The weight function w(rij)has been chosen such that effectively only nearest neighbors contribute to the sum in eq 3. In the present case w(ru) was chosen to be a Gaussian, but it should be stressed that, as long as w(rij)is short-ranged, its precise form is unimportant. g6(r) is defined as g6(r) E ($6(0)*$6(r)) / d r ) (4) where the asterisk denotes complex conjugate and g(r) is the radial distribution function. In a solid (with triangular close-packed
Figure 12. Snapshot of a configuration of 576 hard spherocylinders at 65% of close packing. At this density the system is clearly smectic
Viewpoint as in Figure 10 planes) g6(r)should be long-ranged. In fact, Figure 9 shows that this function is short-ranged for all densities shown. For the denser smectic phase (65% of close packing) there are indications that g6(r) becomes long-ranged. This would suggest that at higher densities there may be a transition from a smectic A phase to a hexatic B phase (note that the transition to the solid cannot take place below the limit of mechanical stability of the solid, i e , at 70% of close packing). A picture is worth a thousand words. Figures 10-12 show the change in molecular ordering that takes place as the system is compressed from isotropic (40% of close packing; Figure 10) to nematic (56.9% of close packing; Figure 11) to smectic (65% of close packing; Figure 12) Acknowledgment. It is a pleasure to acknowledge the collaboration with Mike Allen and Alain Stroobants. This work was performed as part of the research program of the Stichting voor Fundamental Onderzoek der Materie (FOM) with financial support from the Nederlandse Organisatie voor Zuiver Wetenschappelijk Onderzoek (ZWO).
