13544
J. Phys. Chem. 1993,97, 13544-13549
Analysis of Orientational Order in Molecular Clusters. A Molecular Dynamics Study Shimin Xu and Lawrence S. Bartell' Department of Chemistry, University of Michigan, Ann Arbor, Michigan 48109 Received: July 28, 1993; In Final Form: October 20, I 9 9 P
In Monte Carlo and molecular dynamics investigations of the structures of, and transformations between, phases of clusters of plyatomic molecules, the diagnostic functions conventionally applied to atomic clusters have proven to be helpful. These functions include the pair-correlation function g(r), angular distribution function P(B), and Lindemann index The systems of hexafluorides we have been examining exhibit a much richer diversity of phases than do atomic clusters, however, and molecular orientations play crucial roles. It has been fruitful to introduce orientational variants of the above translational functions, which we designate g&), Q(B),and so( 7'). It has been found that, for rapid, subjective assessments of the orientational characteristics of molecules in the various phases, the conventional, two-dimensional Pawley-Fuchs projection is particularly convenient. In analyses of the degree of completion of a phase change in a given cluster, however, several of the other diagnostic functions have proven to be more amenable to quantification. Illustrations of the various proposed functions are presented as applied to seven crystalline phases of small clusters of TeF6.
S(n.
Introduction
Cluster P
The study of atomic and molecular clusters as intermediates between free monomers and the bulk has become increasingly popular.'-3 Molecular dynamics (MD) simulations have been illuminating in the explorationof how properties of systemsevolve as their size increases. Molecular clusters have been particularly attractive systems to investigate in research on structures and transformations of condensed matter because their behavior closely parallels that of bulk matter, even when clusters are very small, and the behavior can be intriguingly complex.67 Techniques for characterizing the various types of structures seen which, for simplicity of expression, we shall refer to as phases, are only now being developed. A number of the diagnostic procedures designed for atomic clusters are proving useful for molecular clusters. These include the pair-correlation function g(r),8 the angular distribution function P(8),8-10and the Lindemann index 6( 2")8J1-13which correlates intermolecular amplitudes of vibration with the onset of melting. Molecular clusters, however, exhibit a much wider variety of phases than atomic clusters, and the orientations of molecules must be addressed, as well. For example, in the hexafluoride systems we have been examining, we have already encountered well over a dozen different solid-state packings,I4 and rotational melting can occur before translational melting. One powerful technique for portraying molecular orientations is that used to good advantage by Pawley and co-workers,5 which we shall call the Pawley-Fuchs projection. While it gives valuable and quickly interpretable graphic information, other representations reducibleto onedimensional numerical functionsare worth considering. The aim of the present paper is partly to introduce some of the new phases found for the hexafluorides, as typified by TeF6, that will be treated in forthcoming MD studies of the phase behavior of quasispherical molecules. Such studies have already been of substantial assistance in the interpretation of experimental electron diffraction experiments on molecular cluster^.'^ They have also provided a fuller understanding of how molecules organize into crystalline aggregates. The primary aim of this presentation, however, is to examine the performance of several diagnostic functions alternative to the Pawley-Fuchs projection, each in some manner a rotational analog of one of the functions g(r), P(8), or 6( T ) mentioned above.
All of the clusters of AFs molecules examined experimentally to date (A = S,Se, Te, Mo, and W) freeze into a body-centered cubic phase (bcc) below their melting points, as do the bulk systems.16 All also transform into a monoclinic phase, space group C2/m,Z = 3 (monol) when cooledcomparativelyrapidly.17 The stable low-temperature phase of TeFs in large clusters and in the bulk is orthorhombic, space group Pnma, Z = 4 (orth).l* This phase is readily prepared in simulations for small clusters, of course, but does not spontaneously form from warmer phases on the time scale of simulations. Other phases of TeF6 that can exist at least metastably in small clusters of comparable stability include trigonal, space group P j m l , Z = 3 (trigl) known experimentally for SFa,I9 WCL2O and UC16,21rhombohedral, space group Rg, Z = 3 (rhmb) known experimentally for WCL22 another trigonal form, space group Pglc, Z = 2 (trig2), and another monoclinic form, space groupP21lc. Z = 2 (mono2). We are unaware of experimental observations of the latter two structures among the hexahalides although trig2 has been postulated elsewhere.23 Because phase mono2 has not been described in the literature to our knowledge, we list the atomic coordinates as calculated by a crystal-packing analysis in Table 1. This phase is of interest in view of its particulary dense packing. Several other crystal organizations to be characterized elsewhere have also been found, but the above seven phases will suffice for the present purposes.
e
Abstract published in Aduonce ACS Abstrocts. December 1, 1993.
0022-3654/93/2091-13544304.00/0
h Of TeF6
Computational Detnils Molecular dynamics runs were based on a modified version of the program MDMPOL.24 Simulations were performed on approximately spherical clusters of tellurium hexafluoride containing 100, 128, 150, 250, and 350 molecules. The molecules were taken to be rigid octahedra with effective bond lengths of 1.8 15 A.25 Pairwise-additive atom-atom intermolecular interactionsbased on Lennard-Jones potential functionswere employed for the clu~ters.~ In runs to be presented elsewhere,16 partial charges on atoms were included in the model potential function in an attempt to make it more realistic. Results for structures and dynamic properties were not changed enough to alter conclusions presented here. Therefore,all illustrativecalculations shown will be based on the model of neutral atoms. As borne out el~ewhere,~ this ~ . ~model ~ gives a good account both of the structures and the transformations of condensed TeF6. Q 1993 American Chemical Society
The Journal of Physical Chemistry, Vol. 97, No. 51. 1993 13545
Orientational Order in Molecular Clusters
TABLE 1: Fnctional Coordinates of Atoms in Monoclinic Cell,' Pzdc, from Packinp Analysis atom u U W 0 0 Te 0 -0.1104 FI 0.2971 -0.0092 Fz 0.2017 -0.3986 0.0371 F3 0.1353 0.7377 0.1866 "The remaining coordinates are generated by the Wyckoff e site symmetry, choice 2. Lattice parameters are u = 5.361 A, b = 5.075 A, c = 8.256 A, 7 = 61.95O.
0.035 0.03
0.025
s
0.02
PI 0.015 0.01 0.005 0 0
30
80
90
e
120
150
180
(0)
a
o.025
7
b
a
0
30
Bo
90
e
120
160
180
(0)
b Figure 3. (a) Angular distribution functions P(0) and (b) orientational angular distribution functions Q(0) of Te atoms in a 150-moleculecluster
ofTeF6at lSOK(solidcurve,bcc)andat50K(dashedcurve,monoclinic). The cutoff distance for nearest neighbors is taken to be 5.25 A at both temperatures for P(0). d
C
Figure 1. Pawley-Fuchs projections of a 350-molecule cluster of TeF6 at various stages of cooling: (a) 140 K,bcc, (b) 100 K,bcc, (c) 90 K, monoclinic; (d) 50 K,monoclinic.
where N is the number of molecules in the cluster, Nb is the number of bonds in the molecule, and is the angle between bond a in molecule i and bond @ in molecule j , or
a
b
C
Figure 2. Pawley-Fuchs projections of various phases of a 128-molecule cluster of TeF6 at 50 K: (a) rhombohedral; (b) orthorhombic; (c) a second monoclinic (mono2, see text).
Diagnostic Fuactiom Each of the different phases of TeFs imposes its own unique signature upon the conventional diagnostic functions g(r) and P(9). Because molecular rotations play a more important role than do molecular translations in some of the phase transitions, changes in g(r) and P(9) can be more subtle than is desirable for the identification of phases. Pawley-Fuchs two-dimensional projections5 display molecular orientations in a particularly accessible way for visual inspection as shown in Figures 1 and 2. Several alternative one-dimensional functions that explicitlytake molecular orientation into account are described below and applied to results of MD simulations of TeF6 clusters. The examples given are illustrative, not exhaustive. One natural choice is to examine the orientational angular distribution function which we designate as Q(0). For a molecule as simple as TeF6 in which all the bondsare equivalent,the simplest expression of such a function is a normalized histogram of all of the angles between directions of bonds in different molecules, or
in which (n& and (np), are the unit vectors along the bonds involved. For example, the function Q(0) for a perfectly ordered bcc crystal of TeF6 would have all covalent bonds parallel to the crystallographic axes and Q(0) would display peaks only at Oo, 90°, and 1 8 0 O . Phases with several distinct orientations would yield more peaks. In realistic simulations, librational (angular) amplitudes of motion characteristically broaden the peaks and broaden them very considerably in the plastically crystalline bcc phase as shown in Figure 3, which compares Q(0) with P(9) for the same phase. Orientational angular distributions for other phases are shown in Figures 4-6. It has proven to be helpful to definea function Si,(r,n) to provide an index to describe how closely matched are the orientations of any two molecules i and j whose centers of mass are separated by distance r. For hexafluorides such a function can be taken to be
where the unit vectors are as defined above and parameter n controls the latitude tolerated in the matching. If parameter n is greater than 6, say, the contribution Sij(a,B) from any angle B a , ~is essentially a Gaussian function of Ba,B with ug2 = 2/n. In our illustrative calculations, n was taken to be 265, corresponding to ug = So. A less severe requirement for alignment might be more appropriate for general application. The function Sij will
Xu and Bartell
13546 The Journal of Physical Chemistry, Vol. 97, No. 51, 1993 0.1
a
0.08
-
0.06
-
W
b
-
0.04 0.02
I
h
. .p,,;; .* :.,... , .i; ',i.. ..\ r.,:; .
:
0 0
30
60
90 8 (")
150
120
10
12
14
16
12
14
16
(A) a
0.04 I
CD
8
r
Figure 4. Angular distribution functions of Te atoms in 128-molecule clusters of TeF6 in the rhombohedral (solid curve) and second monoclinic (mono2, dashed curve) phases at 50 K. Cutoff distance = 5.25 A.
n
6
4
180
2 1.61
j/
I
0.03
3 8
0.02
n
CD
k
4
0.01
60
30
90
e
l20
160
18
correlation functions g&) of Te atoms in a 15o-mo~ccu~e cluster of TeF6 at 150 K (solid curve, bcc) and 50 K (dashed curve, monoclinic).
(0)
Figure 5. Angular distribution functions of Te atoms P(0) (solid curve) and orientational angular distribution functions Q(0) (dashed curve) of an orthorhombic, 128-molecule cluster of TeF6 at 50 K. Cutoff distance = 5.25 A.
6r
1
5t1 i 4
0.045
:3 h
bo
0.035
A
2
n
k
10
8
r (A) b Figure 7. Pair-correlation functions g(r) and (b) orientational pair-
0 0
6
0.026
0 4
0.016
6
8
10
12
14
16
10 12 r (A) b
14
16
r 0.005
I
0
a
I 55
110
165
(A)
22 2.5
T (K) Figure 6. Peak heights of Q(0) = 90' for clusters of TeF6 at various temperatures as they are being cooled: circles, 100-molecule; triangles, 150-molecules; squares, 250-molecule; diamonds, 350-molecule clusters.
be shown elsewhere to be useful in recognizing critical nuclei in phase changes. For the present paper we incorporate it into an orientational pair-correlation function ge(r), which is a variant of the conventional pair-correlation function g(r),modulated by suppressing certain peaks. The orientational pair-correlation functionge(r),then, is defined as (4)
where u is the volume per molecule of the cluster (taken to be the same as that of the bulk in our calculations). If all molecular orientatiOns in the cluster were identical, Sij(r,n) would be equal to unity and g&) would reduce to theconventional pair-correlation function g(r). If several well-ordered orientations exist, g&) disappears for those pairs which are orientated differently. Illustrative examples of ge(r) are compared with g(r) in Figures 7-10.
I
I
t jlEi 4
6
8
figure 8. (a) Pair-correlation functions g(r) and (b) orientational paircorrelation functions g&) of Te atoms in 128-moleculeclusters of TeF6 in the rhombohedral phase (solid curve) and a second monoclinic phase (mono2, dashed curve) at 50 K.
For molecular systems, a more complex function called the angular pair-correlation function had already been intr0duced,2~*~8 for a purpose somewhat analogous to that of ge(r). It is based on threevariables,namely,thedistancebetweenapairofmolecules and two orientational parameters. Although it includes information about both translation and orientation, the result is not so simply plotted because of its multidimensionality. Although
Orientational Order in Molecular Clusters 8
The Journal of Physical Chemistry, Vol. 97, No. 51, 1993 13547
6
0.11
1
I I ‘
1
I
0.00
A
5
bo”
4
0.07 b o
b 3
0.06
h
5
M
2
0.03
1
0.01
0
4
6
8
10
12
14
16
r (A) Figure 9. Pair-correlation functions g(r) (solid curve) and orientational pair-correlation functions g&) (dashed curve) of Te atoms in an orthorhombic, 128-moleculecluster of TeF6 at 50 K.
0
8
10
r
12
14
166
22
Figure 11. Lindemann indices 6 ( T ) for adjacent molecules, averaged over all molecules in clusters of ?‘e&as they are being heat& circlcs,
100-molecule;triangles, 150-molecule;squares 250-molecule;diamonds, 350-molecule clusters.
0‘43 0.33
6
110
T (K)
1
4
66
4
16
(A>
Figure 10. Pair-correlation functions g(r) of Te atoms in trigonal, 128molecule clusters of TeF6 in the first (trigl, solid curve) and second (trig2, dashed curve) trigonal variants at 50 K.
the orientational pair-correlation function we proposed is onedimensional, it contains information about the ordering of both translations and orientations. In addition, an index is introduced to assess the magnitude of the amplitudes of orientational displacementsof pairs of molecules relative to each other. Although its analogy to the traditional Lindemann index’’ 6 of translational vibrations is imperfect, we refer to it as the “orientational Lindemann index” 68. In its construction the variance in the nearest-neighbor distances that are incorporated into 6 is replaced by the variance in the angles between (n& and (n&, and the average linear separation of molecules in 6 is replaced by an angular separation taken to be the angle between two indistinguishable molecular orientations. Assigning a value of 90’ to the angular separation in A&, we have
6, =
T (K) Figun 12. Orientational Lindemann indices de( 7‘)averaged over all the molecules in clusters of TeF6 as they are being heated: circles, 100-
molecule; triangles, 150-molecule;squares,250-molecule;diamonds, 350molecule clusters. commonly adopted variant for small clusters is the definition
in which rijis the separation between any two atoms (or molecules in molecular clusters). Obviously analogous to eq 5 , such a definition includes all internuclear distances in a cluster instead of just those distances correspondinp to atomic contacts. Only the latter were intended by Lindemann to be considered. Therefore, in calculations leading to Figure 11 we included in 6 only those rij corresponding to molecular contacts (that is, for TeF6, only those within a cutoff distance of 7.5 A). On the other hand, since the variance in 8 4 does not appear to change markedly with the distance between i and j, we retain eq 5 for 68, summing over all pairs in a cluster to reduce noise.
Discussion
When amplitudes of librational motions of molecules are small and molecules all have the same mean orientation, 68 is small. As will be discussed elsewhere, 68 provides a satisfactory criterion for rotational melting. The functions 6 and 68 are compared in Figures 11 and 12 for clusters of tellurium hexafluoride that are undergoing phase changes on heating. Because several different definitions of the translational Lindemann index 6 have been introduced in the literature, we mention the most common. Some authors use the ratio of the rms amplitude of vibration of atoms about their equilibrium positions to the mean distance separating the atoms,29a definition closely following Lindemann’s original usage. * A perhaps more
It is not the purpose of the present illustrative examples to characterize the various structures of TeF6 clusters which heve been encountered, for these are already known from a detailed analysisof thecoordinatesof theconstituent atoms in theclusters. Instead, the aim is to find which, if any, of the proposed functions provide readily calculated indices permitting a rapid, routine recognition of transitions from one phase to another and, less importantly, a diagnosis of which phases are present. As confirmed by Figure 1, Pawley-Fuchs projections offer a very convenient and intuitively appealing way to monitor the progress of a phase change. The view down one of the crystallographic axes of the stable high-temperature bcc phase shows that all bonds tend to be aligned along the crystallographic axes but with considerable disorder in orientation. At about 90 K the cluster transforms to monoclinic (monol) by reorienting one-third of its molecules, as can be seen in the figure. In Figure 2 are shown
13548
Xu and Bartell
The Journal of Physical Chemistry, Vol. 97, No. 51. 1993
projections of the rhombohedral, orthorhombic, and second monoclinic (mono2) phases. In the rhombohedral structure, all molecules have the same orientation so that, if the projection were viewed down one of the bond directions instead of approximatelyalong a molecular 3-fold axis, the projection would be indistinguishablefrom that of the bcc phase shown in Figure la,b except for the much smaller thermal disorder. In Figure 2b,c it is evident that sites corresponding to molecules with two distinct orientations exist in both the orthorhombic and the monoclinic (mono2) phases. Because the molecules fit together in such different ways in these two phases, there is no difficulty in distinguishingbetween the phases by means of the projections. As illustrated in Figures 3-6, the angular distribution function P(0) and orientationalangular distribution function Q(0) provide additional information. Figures 3a and 7a, considered together, indicate how closely the phase monol resembles bcc in spatial arrangements of molecular centers of mass. The broad peaks of bcc in P(0) split as the lowering of symmetry characterizing the monoclinic phase develops, but the overall envelope changes only modestly. It is, of course, the similarity of organization of the two phases that accounts for the transition being so facile that it can occur on the time scale of MD runs. The orientational angular distribution functionsQ(0) of the two phases are strikingly different, however, because the transition, as mentioned above, involves the reorientation of one-third of the molecules. The function Q(0) also demonstrates how great is the orientational disorder in the bcc phase, which phase would display peaks only at Oo, 90°, and 180°, if it were perfectly ordered. If molecular orientations were completelyrandom,thedistribution Q(0) would become proportional to sin 6, a function whose shape the orientational angular distribution of the bcc phase has begun to approach. The more conventional function P(0) distinguishes reasonably well between all of the phases, although monol and mono2 are somewhat similar. Function Q(0), while clearly discriminatingbetween bcc and monoclinic, is less successful with some of the others. For example, the curves for monol (Figure 3b) and orthorhombic (Figure 5 ) are nearly the same. As suggested in Figure 6,however, the temperature dependence of Q(0) at I9 = 90° can provide a convenient monitoring of the transition from bcc to monoclinic. Selecting a threshold value of the index will be shown elsewhere to provide a useful measure of the characteristic temperature of the transition.jO Figures 7-10 show that each of the seven phases can be distinguished by a careful examination of the g(r) function. Such examinations are not particularly convenient for the purpose of identifying phases, however. Introducing the modulation by Sij(r,n) adds additional discrimination such that the attenuation of certain peaks by disorder or by reorientation affords simple, alternative indices of transformation. For example, the distribution of intermolecular distances in the two monoclinic phases monol and mono2 is similar according to the g(r) functions in Figures 7a and 8a. Nevertheless, the two phases differ considerably in their orientational pair correlation functions (Figures 7b and Sb). For some phases, e.g., rhombohedral among the present examples, g(r) and g&) are the same except for the modest motional attentuation in the latter function, because all of the molecules have identical orientations. All of the phases described in this paper except bcc are metastable at low temperatures, and all ultimately transform to bcc on warming. Even though orientations of all molecules in bcc are identical over long time averages, the motional attenuation (Figures 7a and 8a) is so great that the transition from any of the other phases to bcc can be readily followed by monitoring the decay of a strong feature in g&). In addition to providing a way to probe the course of a transition, the g&) function reveals at a glance which of the intermolecular distances in g(r) for a given phase correspond to molecules with common orientations. An illustrative example is shown in Figure 9.
From the foregoing it can be seen that the monitoring of each of the proposed diagnostic functions at a suitable value of the argument can yield a working criterion, not only for melting but also for other phase changes encountered. Perhaps the most popularly applied criterion for the phase change of melting in atomic clusters has been the Lindemann index 6. It is instructive to examine the behavior of this index in transitions of molecular clusters, both according to the conventional variant, 6, and according to the angular analog, 68. This is done in Figures 11 and 12. When small clusters are heated, their melting is diffuse and extended over a substantial range of temperatures. No sharp melting point exists. Cluters of small, rigid polyatomic molecules appear to become thoroughly molten when 6 approaches 0.09 as deduced from the coefficients of diffusion of their constituent molecules.25 Figure 11 indicates the melting and shows that smaller clusters melt (as expected) at lower temperatures than do larger clusters. An additional inflection at a lower value of 6 also occurs which has no counterpart in the published studies of atomicclusters. It correspondstothetransitionfrommonoclinic to bcc. The behavior of the index 60 is markedly different from that of 6, as is illustrated in Figure 12 where 68 reveals only a single transition for each cluster. This transition corresponds to rotational melting. When molecular orientations become completely random over the time steps monitored, the index 60 approaches its limiting value of (1 or 0.4352. A somewhat lower value of about 0.33 has been found to be a more practical threshold to be associated with the onset of rotational melting over intervals of 4000 time steps.30 On this basis it has been inferred that the rotational melting of TeF6 clusters is indistinguishable from the monoclinic to bcc solid-state transition from an ordered crystalline array to a plastically crystalline, orientationally disordered phase. Although each of the above diagnostic functions can play a role in helping to identify crystalline phases, their most useful application is probably in the monitoring of the course of a phase transition. For a quick, easily interpreted assessment of the state of a system, the Pawley-Fuchs projection is indispensable. For deriving estimates of the fraction of a cluster to have undergone a transformation, however, it has been found that several of the other functions are more amenable to quantification. Further details of their performance will be presented in forthcoming papers in this series on MD simulations of the structures and transformations of small molecular clusters.
Acknowledgment. This research was supported by a grant from the National Science Foundation. We thank Dr.F. Dulles and Mr. K. Kinney for their valuable assistance with computations. References and Notes (1) Jena, P., Rao, B. K., Khanna. S.N. Eds. Physics and Chemistry of Small Clusters; Plenum: New York, 1987. (2) Echt, O., Recknagel, E., Eds. Proceedings of the Fifrh International Symposium on Small Particles and Inorganic Clusters. In 2.Phys. D 1991, 19dt20.
(3) Jena, P., Rao, B. K., Khanna, S.N., Eds. Physics and Chemistry of Finite Systems: From Clusters to Crystals; Kluwer Academic: Dordrecht, 1992. (4) Boyer, L. L.; Pawley, G. S.J . Comput. Phys. 1988, 78, 405. ( 5 ) Fuchs, A. H.; Pawley, G . S. J . Phys. (Paris) 1988,19, 41. (6) Bartell. L. S.; Dibble, T. S.;Hovick, J. W.; Xu, S. In Physics and Chemistry of Finite Systems: from Clusters to Crystals; Jena, P., Rao. B. K., Khanna, S. N., Eds.; Kluwer Academic: Dordrecht, 1992; Vol. I, p 71. (7) Bartell, L. S.; Xu, S.J. Phys. Chem. 1991, 95, 8939. (8) Berry, R. S.;Beck, T. L.; Davis, H. L.; Jellinek, J. In Euolution of Size Effects in Chemical Dynamics, Part 2; Prigogine, I., Rice, S.A., Eds. (Adu. Chem. Phys. 1988, 70, (2). 75). (9) Haymet. A. D. J. Chem. Phys. Irtt. 1984, 107, 7 7 . (IO) Quirke, N.;Sheng, P. Chem. Phys. Lett. 1984.110, 63. ( 1 1) Lindemann, F. A. Phys. 2.1910, I I , 609; Engineering 1912,94515. (12) Kaclberer, J. B.; Etters, R. D. J . Chcm. Phys. 1977.66, 3233. (13) Etters, R. D.; Kaelberer, J. B. Phys. Reo. A 1975, 11, 1068. (14) Xu, S.;Kinney, K. E.; Hovick, J. W.; Bartell, L. S. Unpublished research.
Orientational Order in Molecular Clusters (15) Dibble, T. S.; Bartell, L. S . J. Phys. Chem. 1992, 96, 8603. (16) Siegel, S.; Northrup, D. A. Inorg. Chem. 1966,5, 2187. (17) Bartell, L. S.; Hovick, J. W.; Dibble, T. S.; Lennon, P. J. J . Phys. Chem. 1993, 97, 230. (18) Bartell, L. S.;Powell, B. M.Mol. Phys. 1992, 75, 689. (19) Raynerd, G.; Tatlock, G. J.; Venables, J. A. Acto Crystullogr. 1982, 838, 1896. (20) Taylor, J. C.; Wilson, P. W. Acto Crystallogr. 1974, 830, 1216. (21) Taylor, J. C.; Wilson, P. W. Acta Crystullogr. 1974, 830, 1481. (22) Ketelaar, J. A. A.; van Oosterhout, G. W. Red. Trau. Chim. PuysBus 1943.62, 197.
The Journal of Physical Chemistry, Vol. 97, No. 51, 1993 13549 (23) Willing, W.; MOller, U. Acta Crystallogr. 1988, ,344, 1. (24) Smith, W.; Fincham, D. CCPS Program Documentation: The Program MDMPOL; SERC Daresbury Laboratory: Dareebury, UK, 1982. (25) Gundersen, G.; Hedberg, K.J . Chem. Phys. 1978, 68, 3548. (26) Xu,S.; Bartell, L. S . J . Phys. Chem., following paper in this issue. (27) Haile, J. M.;Gray, C. G. Chem. Phys. Lcrr. 1980, 76, 583. (28) Gray, C. G.;Gubbins, K.E. Theory of Moleculur Fluids,Clarendon: Oxford, 1984; Vol. 1. (29) Stillinger, F. H.; Stillinger, D. K.J . Chem. Phys. 1990, 93, 6013. (30) Xu,S.; Bartell, L. S. 2.Phys. D, submitted.
