Free Volume in the Hard­Sphere Liquid 
Srikanth Sastry (1) , Thomas M. Truskett (1) , Pablo G. Debenedetti (1)\Lambda , 
Salvatore Torquato (2;3) , and Frank H. Stillinger (4;2) 
(1) Department of Chemical Engineering 
Princeton University, Princeton, NJ 08544 
(2) Princeton Materials Institute 
Princeton, University, Princeton, NJ 08544 
(3) Department of Civil Engineering and Operations Research, 
Princeton, University, Princeton, NJ 08544 
(4) Bell Laboratories, Lucent Technologies 
Murray Hill, NJ 07974 
\Lambda Electronic address: pdebene@pucc.princeton.edu 

Abstract 
We present a method for the efficient calculation of free volumes and cor­ 
responding surface areas in the hard­sphere system by extending a previous 
method for calculating, exactly, cavity volumes in sphere packings. Using 
this method, we evaluate, for the first time, the free­volume distribution of 
the hard­sphere liquid over a range of densities near the freezing transition. 
From the distribution of free volumes, the equation of state can be obtained 
from a purely geometric analysis, which permits the calculation of pressure in 
Monte Carlo simulations where the dynamic definition cannot be employed. 
Furthermore, we obtain the cavity­volume distributions indirectly from the 
free­volume distributions in a density range where direct measurement is in­ 
adequate. Direct measurement of the first moment of the cavity­volume dis­ 
tribution enables us to calculate the chemical potential in the vicinity of the 
freezing transition. 

I. INTRODUCTION 
It is well established that the structure of most dense liquids is dominated by repulsive 
interactions. The simplest model liquid that embodies this feature is the hard­sphere fluid, 
in which impenetrable particles interact solely via hard­core repulsions. The hard­sphere 
system has played a major role in liquid­state theory ever since seminal investigations by 
computer simulation [1--3] strongly suggested that it exhibits a first­order freezing transition. 
In the late 1950s, comparable strides occurred on the theoretical front with the introduction 
of the scaled­particle theory of fluids [4], which offered simple and accurate equations of state 
for systems comprising hard­core molecules. These contributions allowed Longuet­Higgins 
and Widom [5], and later Guggenheim [6], to extend the van der Waals theory by refining the 
repulsive contribution to the equation of state. In turn, the theories yielded quantitatively 
accurate predictions for the melting properties of argon. 
Following the high­temperature expansions of Zwanzig [7], Barker and Henderson [8,9] 
introduced their landmark perturbation theory as applied to a simple fluid in 1967. In their 
approach, the unperturbed reference system consisted of the positive part of the Lennard­ 
Jones potential, which in turn was related to an essentially equivalent hard­sphere system. 
The remarkable success of the theory of Barker and Henderson demonstrated quantitatively 
that simple liquids near their triple point are removed by a minor perturbation from a purely 
repulsive fluid. This result was reaffirmed by the first­order perturbation expansion of Weeks, 
Chandler and Andersen [10,11]. These methods inspired researchers [12,13] to explore in 
detail the structure of the equilibrated hard­sphere fluid. In particular, the perturbation 
techniques motivated Verlet and Weis [14] to develop a semi­empirical parameterization of 
the radial distribution function for the hard­sphere liquid, which has become the standard 
for numerical calculations. 
Just as the hard­sphere fluid provides a reference for understanding liquid structure, it 
also represents the simplest system which exhibits a fluid­solid transition and, possibly, a 
glass transition [15,16]. The fact that the properties of the hard­sphere liquid arise from 
1 

strictly entropic contributions, that is to say from purely geometric considerations, underlies 
the continued interest in this system, with recent emphasis directed towards understanding 
the statistical geometry of dense sphere packings [16--31]. In particular, quantities that 
describe the void space (volume available for insertion of an additional hard sphere), the 
free volume (volume within which a given hard­sphere center can move without requiring 
alteration of the other sphere positions), and the corresponding surface areas are directly 
related to thermodynamic quantities. 
A computational method has recently been presented by us [32] which permits exact 
determination of cavity volumes and surface areas in three­dimensional mono­ and polydis­ 
perse spherical packings. We utilize this method here to evaluate the chemical potential of 
the hard­sphere liquid over a range of densities near the freezing transition. Furthermore, 
we extend the above method for the calculation of free volumes and use this information to 
predict the equation of state of the fluid near the freezing density. We also determine, for 
the first time exactly, the free­volume distribution in the hard­sphere system over a range 
of densities. The sparcity of void space at high densities renders the cavity­volume distribu­ 
tion statistically inaccessible by the direct approach. However, it will be demonstrated that 
the cavity­volume distribution can be deduced indirectly from the free­volume distribution, 
which can itself be determined with a high degree of precision. 
In Sec. II we present definitions and background information. In Sec. III we describe 
the method for calculating free volumes in sphere packings. The results of our calculations 
of chemical potentials, pressures, and free­ and cavity­ volume distributions in the dense 
hard­sphere liquid are presented in Sec. IV. Sec. V contains a summary and concluding 
remarks. 
II. BACKGROUND 
In a system containing N hard spheres, a geometrical free volume v f can be defined [33] 
as the volume over which the center of a given sphere can translate, given that the other 
2 

N cavity volume 
v, which is the volume of a connected region of space available for the addition of another 
sphere. By definition, a point is inside of a cavity if it lies outside of the exclusion spheres 
surrounding each particle center, i.e., if it is separated from each particle center by at 
least one hard­core diameter oe. At low densities, the void space present in the system is 
connected, and hence the free volume approaches the void volume. As the density of the 
system is increased, the void space becomes disconnected, corresponding to the percolation 
of the exclusion spheres. This change in topography of the void space occurs when the 
exclusion spheres occupy approximately 30% of the space [34], which translates to a reduced 
density of ae \Lambda 
c ß 0:076 [35], where ae \Lambda = N oe 3 =V and V is the system volume. Here, we focus 
on the statistical geometry of the high­density liquid, i.e., systems well above the percolation 
threshold. 
Speedy and Reiss [36] have demonstrated that the free­volume and cavity­volume distri­ 
bution functions are related. For completeness, their arguments are reproduced below. In 
what follows p(v)dv is the probability that a cavity has a volume between v and v+dv. while 
f(v f )dv f is the probability that the free volume of a sphere lies between v f and v f + dv f . 
Analogous probability densities can be defined for the cavity surface p s (s) and the free 
surface f s (s). 
For a given configuration of spheres, the union volume of the cavities represents the 
available space V 0 . The available surface area S 0 comprises the surface areas of the individual 
cavitites. The average cavity volume and surface area are given by 
hvi = hV 0 i 
N c 
= 
Z 1 
0 
xp(x)dx hsi = hS 0 i 
N c 
= 
Z 1 
0 
yp s (y)dy (1) 
where N c represents the number of cavities in the system, which is averaged over all real­ 
izations of the particles. Speedy [37] has shown that the equation of state of an equilibrium 
hard­sphere fluid can be expressed in terms of the statistical geometry of the cavities 
fiP 
ae 
= 1 + oe 
2D 
hsi 
hvi 
(2) 
3 

where P is the pressure, ae is the number density, fi is (kT ) 
Boltzmann [39] may have 
been the first to derive this result. When an additional sphere is added to the system, it 
enters a cavity of size x which becomes its free volume. Since the sphere additions sample 
the available space uniformly, a cavity is visited with a frequency that is proportional to its 
size, i.e., 
f(x)dx = xp(x)dx 
R 1 
0 xp(x)dx 
= xp(x)dx 
hvi : (3) 
Strictly speaking, this equation relates the free­volume distribution of a system with N + 1 
spheres to the cavity­volume distribution of the N­sphere system. These two systems become 
equivalent in the thermodynamic limit. 
For any quantity g(v f ) that depends on the free volume, the following relationship holds 
hg(v f )i = 
R 1 
0 g(x)xp(x)dx 
hvi = 
hvg(v)i 
hvi : (4) 
In particular, if we choose g(v f ) = v n 
f , then (4) yields an equation relating the moments of 
the free­ and cavity­volume distributions 
hv n 
f i = hv n+1 i 
hvi : (5) 
From (5), it follows that 
hv 
free surface area s f . Choosing g(v f ) to be 
the ``free surface''­to­``free volume'' ratio reveals the following valuable relationship 
h s f 
v f 
i = hsi 
hvi 
= hS 0 i 
hV 0 i 
(7) 
which was originally proved by Speedy [40] using a slightly different argument. From (2) 
and (7), it is apparent that the pressure can be deduced from free­volume information alone 
4 

fiP 
ae 
= 1 + oe 
2D h s f 
v f 
i: (8) 
This equation was suggested 25 years ago by Hoover et al. [33] when they considered the 
dynamics of a light particle in a classical system. 
The chemical potential of the hard­sphere system is also directly related to its statistical 
geometry 
¯ = kT ln 
/ 
– D N 
hV 0 i 
! 
= kT ln 
/ 
– D Nh1=v f i 
N c 
! 
(9) 
where N is the number of particles in the system and – is the familiar thermal wavelength. 
The second equality follows from (1) and (5) and establishes the connection between the 
chemical potential and the free­volume distribution. Notice that the number of cavities 
N c appears in the second relationship, indicating that the chemical potential cannot be 
determined from free­volume information alone. It is conventional to separate the chemical 
potential ¯ into an ideal and an excess contribution 
¯ = ¯ id + ¯ ex = kT ln 
/ 
– D N 
V 
! 
+ kT ln 
/ 
V 
hV 0 i 
! 
: (10) 
The former term represents the chemical potential for an ideal gas, while the latter embodies 
the reversible work required to form a cavity of radius oe. 
Both analytical and numerical methods have been employed [40--43] to study the cavity­ 
volume and free­volume distributions in two dimensions (hard disks). The present work will 
focus on the exact determination of these quantities for the three­dimensional system. 
III. METHODOLOGY 
The algorithm described below is an extension of the method proposed by Sastry et 
al. [32] to calculate cavity volumes and surface areas in particle packings. For brevity, only 
a summary of the method is given; the interested reader can find more details in their 
original paper. 
5 

Given a configuration of hard spheres, the first step in the algorithm is the generation of 
Voronoi and Delaunay tessellations. Both of these constructions divide space into distinct, 
non­overlapping regions. The Voronoi tessellation divides the system into convex polyhedra 
which surround each atom. Specifically, a Voronoi polyhedron V i consists of points closer 
to atom i than any other atom. There will be several polyhedra V k which share a face 
with V i . The atoms corresponding to the neighboring polyhedra V k are termed ``geometric 
neighbors'' of atom i. The Delaunay tessellation is obtained by connecting all geometric 
neighbors, forming a ``primitive graph'' of Delaunay simplices (triangles in 2D, tetrahedra 
in 3D). A schematic of the dual construction is given in Figure 2. A cavity corresponds 
to a percolation cluster of Voronoi edges that lie entirely within the available space (i.e., 
outside of the exclusion spheres) [44,45]. Sastry et al. [32] have demonstrated that a cavity 
is enclosed entirely by the Delaunay simplices which are dual to the Voronoi vertices within 
the cavity. 
To calculate the volume and surface area of each cavity, the corresponding Delaunay 
simplices are subdivided into a set of (generally overlapping) subsimplices. The advantage 
of this division, which is described in detail elsewhere [32], is that for each subsimplex, only 
one exclusion sphere must be considered for determining its contribution to the surface area 
and volume. With an appropriate sign designation, the subsimplex contributions can be 
added directly to yield the total surface area and volume for each cavity (by design this 
method avoids the cumbersome calculation of multiple sphere overlaps). 
To determine the free volume and the free surface area for a given sphere, the sphere 
is simply removed from the system. The cavity that contains the center of the sphere so 
removed is then analyzed using the aforementioned algorithm. Of course when a sphere is 
removed from the configuration, a local region surrounding the resulting cavity must be re­ 
tessellated. The efficiency of this routine can be maintained if the region to be re­tessellated 
is minimal. Fortunately, some properties of the dual tessellation allow for the determination 
of such a minimal region: 
Theorem 1: If an atom is removed from the system, only the Voronoi polyhedra of its 
6 

geometric neighbors must be re­tessellated. 
To see this, consider the removal of atom i from the configuration. By definition, the 
only Voronoi polyhedra that are affected are those whose atoms are closer to a point in V i 
than any other atom. To prove the above theorem, we must demonstrate that at least one 
geometric neighbor k is closer than any non­neighboring atom j to an arbitrary point p in 
V i . 
Given a point p in V i , consider the nearest non­neighboring atom j (see Figure 3). If a 
vector, r jp , is drawn to connect the point of interest to atom j, then it will intersect V i at 
some point p 0 . The point p 0 will be on the face shared between V i and V k , and the definition 
of V k requires that jr jp 0 j ? jr kp 0 j. Furthermore jr jp j = jr jp 0 j + jr pp 0 j and jr kp j Ÿ jr kp 0 j + jr pp 0 j 
(by the triangle inequality). It then follows that jr jp j ? jr kp j, and the theorem is proved. 
Theorem 2: Pairs of geometric neighbors of atom i that share a Voronoi face continue 
to do so after atom i is removed. 
Consider atoms k and k 0 which are geometric neighbors of atom i and share a common 
face. Any point on the common face is closer to atoms k and k 0 than any other atom in the 
system. Clearly the removal of any other atoms (including atom i) will not change this fact, 
verifying the theorem. 
Theorem 1 and Theorem 2 indicate that the minimal region in which the tessellation 
must be reconstructed after the removal of atom i is the superpolyhedron S i composed of 
all Delaunay simplices that share atom i as a vertex. This region is illustrated in Figure 4. 
When atom i is removed, the tessellation is reconstructed inside of S i , and the surface area 
and volume of the cavity that contained the center of atom i can be calculated directly. 
The above procedure is applied to each atom in the system, for each configuration con­ 
sidered. The procedure described above is very efficient, consuming less than two minutes of 
CPU time per configuration (for 500 hard spheres on an HP 715/100 workstation). In Sec­ 
tion IV, results are presented for the equilibrated hard­sphere liquid and a modest extension 
into the metastable region. 
7 

IV. RESULTS FOR THE HARD­SPHERE LIQUID 
Systems of N = 500 spheres of diameter oe were simulated in a cubic cell of volume V using 
a standard molecular dynamics (MD) algorithm [31,46]. The initial configuration was chosen 
to be a face­centered cubic lattice at a reduced density ae \Lambda = 0:80, where ae \Lambda = N oe 3 =V . In each 
case, the lattice was melted by simulating for 5000N collisions to obtain the equilibrated 
fluid. Higher densities were achieved by allowing the diameter of the spheres to increase 
linearly in time via the prescription of Lubachevsky and Stillinger [46,47]. This compression 
protocol will, in general, create a non­equilibrium state at the density of interest. The 
properties of this state will be an extremely complex function of the system's history. To 
remove these effects, compressed packings were allowed to relax over a period of 3500N 
sphere collisions, which was sufficient to guarantee reproducible thermodynamic properties. 
In order to explore the statistical geometry of the dense liquid, several packing fractions 
were investigated in the vicinity of the freezing transition (ae \Lambda 
f ß 0:943). In particular, runs 
were performed at reduced densities ae \Lambda = 0.80, 0.85, 0.90, 0.91, 0.93, 0.943, 0.95, and 0.96. 
Strictly speaking, any amorphous packing with a density ae \Lambda ? ae \Lambda 
f exists in a state that 
is metastable with respect to formation of the crystalline phase. However, the entropic 
barriers to crystallization are large for modest extensions along the metastable branch. 
Indeed, Speedy [31] has found that metastable hard­sphere systems with ae \Lambda ! 1:03 will not 
crystallize even after 10 5 collisions per particle. 
As can be seen from (8), the equation of state for the hard­sphere fluid can be determined 
from free­volume considerations alone. The exact algorithm presented in Sec. III allows for 
the first direct test of (8) by computer simulation. For the calculation of the pressure, 500 
configurations (separated by 10 4 collisions each) were stored at each state point. Results 
for the 500­sphere system are shown in Figure 5. Given the relatively small number of 
configurations considered, the agreement between the pressure calculated from (8) and from 
the virial (collision rate) is remarkable. As a check, the pressure was also determined from 
the free­volume distributions of configurations generated by a series of Monte Carlo (MC) 
8 

simulations (see the discussion of the chemical potential results for details on the Monte 
Carlo runs). The equation of state produced by MD and MC routes were statistically 
indistinguishable. In principle, equation (2) provides yet another geometric route to the 
equation of state. However, the available space vanishes for many dense configurations of 
10 2 
:03 ! ae \Lambda ! 1:11) [31,48]. It has been recognized that the statistical­mechanical for­ 
malism used to desribe such states must be modified through the introduction of specific 
constraints [49--52]. Furthermore, Corti et al. [53] have demonstrated that geometrical con­ 
straints can be applied to simulate superheated liquids using a Monte Carlo algorithm. In 
their investigation, the liquid was prevented from boiling by constraining the size of the 
largest void in the system. For the metastable hard­sphere liquid, the corresponding con­ 
straint should prevent the formation of crystallites. Rintoul and Torquato [48] have shown 
that a bond­orientational order parameter [54] can be invoked to filter out configurations 
that contain significant crystallization. Due to the need for repeated enforcement of con­ 
straints, MC simulations have an obvious advantage over their deterministic counterpart 
(MD) for simulating metastable phases. However, efficient and accurate methods for cal­ 
culating the hard­sphere equation of state have been lacking in Monte Carlo simulations, 
where the dynamic defintion cannot be employed. The free­volume algorithm presented here 
provides one such efficient route to the pressure, requiring only static information. 
The distribution of free volumes f(v f ) is shown in Figure 6 for densities in the vicinity 
of the freezing transition. Notice that there seems to be a smooth change in the behavior of 
the free­volume distribution as the fluid enters the metastable region. It is not known if the 
free­volume distribution continues to vary in a regular way as the fluid is compressed along 
the metastable extension of the fluid branch. If a thermodynamic glass transition occurs 
in the metastable region, it will be a result of structural arrest. Such a profound signature 
of attenuated particle mobility should be evidenced by a change in form of the free­volume 
9 

distribution. Studies are underway to probe the statistical geometry of the dense, metastable 
fluid. 
It should be noted that in one dimension the free­volume distribution is known ex­ 
actly [55] and is given by 
f(v f ) = v f 
hv f i 2 
exp 
/ 
et al. [41] proposed 
the following form for the free­volume distribution for hard disks (above the percolation 
threshold) 
f(v f ) / v ff 
f exp( 
ff = 0:28 
volume diverges, 
i.e., p(0) = 1. We note that there is no fundamental problem with a probability density 
diverging, as long as the integrated probability is suitably defined for a given interval. More 
precisely, the probability density must be non­negative and normalize to unity. Hence, the 
form proposed by Hoover et al. [41] is both plausible and accurate for the three­dimensional 
hard­sphere system near its freezing point. 
Although the free volume is positive for every particle that is not rigidly jammed, the 
available space is identically zero for many configurations in the dense liquid. Therefore, 
direct measurement of the entire cavity­volume distribution is futile for densities near the 
freezing transition. Fortunately, equation (3) provides an indirect method for determining 
the cavity­volume distribution from the free­volume distribution. Results are presented 
in Figures 7 and 8 for the liquid at reduced densities of ae \Lambda = 0.8 and 0.943 respectively. 
For the system at ae \Lambda = 0:8, the free­volume data provides information about the tail of 
10 

the distribution where direct sampling fails. As can be seen in Figure 8, even less can be 
determined by direct measurement at the freezing transition. In simulations of reasonable 
size, the fluctuations which give rise to the tail of the cavity­volume distribution are so rare 
as to be virtually nonexistent. 
Although it is difficult to obtain the entire distribution of cavity volumes at high densities, 
reasonable statistics can be obtained for the average cavity size hvi and, through (1) and (10), 
the excess chemical potential. This is true, in part, because the large cavities that contribute 
to the tail of the distribution are indeed rare occurrences. To measure the excess chemical 
potential directly, 5000 configurations were generated by a standard NVT Monte Carlo 
algorithm for the densities ae \Lambda = 0.8, 0.85, 0.9, 0.91, 0.93, and 0.943. At each density, a face­ 
centered cubic lattice was melted for 10 5 Monte Carlo cycles (attempted moves per particle). 
The configurations were saved in intervals of 200 cycles during a 10 6 cycle production run. 
The available volume of each configuration was measured via the method of Sastry 
et al. [32], and the excess chemical potential ¯ ex was determined using equation (10). 
The results are presented in Figure 9 along with the excess chemical potentials consis­ 
tent with the accurate hard­sphere equations of state developed by Carnahan­Starling [56] 
and Sanchez [57]. For comparison, the precise calculations of Attard [58] and Labik and 
Smith [59] are also shown. The standard deviation was estimated by blocking the configura­ 
tions into 20 sub­sets. As can be seen, good agreement is obtained for all densities including 
the freezing transition. 
V. CONCLUSIONS 
A methodology for determining exactly the free volume and free surface area of a given 
particle in a configuration of hard spheres is presented. Using previously derived identi­ 
ties [36,39,40] that relate the statistics of the free­volume distribution to the hard­sphere 
equation of state, the pressure was determined in the vicinity of the freezing transition. The 
efficiency of this algorithm allows for the pressure to be determined precisely in a Monte 
11 

Carlo simulation, where the collision rate is inaccessible. 
Free­volume distributions for the dense, hard­sphere fluid are characterized for the first 
time. The distributions provide an indirect route to information about the statistics of cavity 
volumes. This is significant because the infrequent appearance of void space at high densities 
prevents the direct measurement of such quantites. Characterization of both cavity­ and 
free­volume distributions should prove to be interesting for metastable sphere systems, where 
the statistical geometry is poorly understood. 
It is shown that the first moment of the cavity­volume distributions, i.e., the average 
cavity size, can be obtained by direct measurement. This quantity was calculated for the 
hard­sphere liquid in the vicinity of the freezing transition. From this information the excess 
chemical potential was determined and was found to be in good agreement with previously 
tabulated results. 
ACKNOWLEDGEMENTS 
We thank Robin Speedy and David Corti for helpful discussions. PGD gratefully ac­ 
knowledges support of the U.S. Department of Energy, Division of Chemical Sciences, Office 
of Basic Energy Sciences (Grant DE­FG02­87ER13714), and of the donors of the Petroleum 
Research Fund, administered by The American Chemical Society. TMT gratefully acknowl­ 
edges The National Science Foundation for a graduate research fellowship. 
12 

REFERENCES 
[1] Alder, B. J., and Wainwright, T. E., 1957, J. chem. Phys., 27, 1208. 
[2] Alder, B. J., and Wainwright, T. E., 1962, Phys. Rev., 127, 359. 
[3] Hoover, W. G., and Ree, F. H., 1968, J. chem. Phys., 49, 3609. 
[4] Reiss, H., Frisch, H. L., and Lebowitz, J. L., 1959, J. chem. Phys., 31, 369. 
[5] Longuet­Higgins, H. C., and Widom, B., 1964, Molec. Phys., 8, 549. 
[6] Guggenheim, E. A., 1965, Molec. Phys., 9, 43. 
[7] Zwanzig, R. W., 1954, J. chem. Phys., 22, 1420. 
[8] Barker, J. A., and Henderson, D., 1967, J. chem. Phys., 47, 4714. 
[9] Barker, J. A., and Henderson, D., 1976, Rev. Mod. Phys., 48, 587. 
[10] Chandler, D., Weeks, J. D., and Andersen, H. C., 1983, Science, 220, 787. 
[11] Chandler, D., Weeks, J. D., and Andersen, H. C., 1976, Adv. chem. Phys., 34, 
105. 
[12] Barker, J. A., and Henderson, D., 1971, Molec. Phys., 21, 187. 
[13] Stillinger, F. H., 1971, J. comput. Phys., 7, 367. 
[14] Verlet, L., and Weis, J.­J., 1972, Phys. Rev. A, 5, 939. 
[15] Woodcock, L., 1981, Ann. N. Y. Acad. Sci., 371, 274. 
[16] Speedy, R. J., 1993, Molec. Phys., 80, 1105. 
[17] Hiwatari, Y., Saito, T., and Ueda, A., 1984, J. chem. Phys, 81, 6044. 
[18] Reiss, H. R., and Hammerich, A. D., 1986, J. Phys. Chem., 90, 6252. 
[19] Tobochnik, J., and Chapin, P., 1988, J. chem. Phys., 88, 5824. 
13 

[20] Zhou, Y., and Stell, G., 1988, J. statist. Phys., 52, 1389. 
[21] Torquato, S., Lu, B., and Rubinstein, J., 1990, Phys. Rev. A, 41, 2059. 
[22] Gonzalez, D. J., and Gonzalez, L. E., 1991, Molec. Phys., 74, 613. 
[23] Nezbeda, I., and Smith, W. R., 1992, Molec. Phys., 75, 789. 
[24] Bowles, R. K., and Speedy, R. J., 1994, Molec. Phys., 83, 113. 
[25] Yeo, J., 1995, Phys. Rev. E, 52, 853. 
[26] Torquato, S., 1995, Phys. Rev. E, 51, 3170. 
[27] Rintoul, M. D., and Torquato, S., 1996, Phys. Rev. Lett., 77, 4198. 
[28] Yuste, S. B., de Haro, M. L., and Santos, A., 1996, Phys. Rev. E., 53, 4820. 
[29] Dasgupta, C., and Valls, O. T., 1996, Phys. Rev. E, 53, 2603. 
[30] Reiss, H., Ellerby, H. M., and Manzanares, J. A., 1996, J. Phys. Chem., 100, 
5970. 
[31] Speedy, R. J., 1997, J. Phys.: condens. Matter, 9, 8591. 
[32] Sastry, S., Corti, D. S., Debenedetti, P. G., and Stillinger, F. H., 1997, 
Phys. Rev. E, 56, 5524. 
[33] Hoover, W. G., Ashurst, W. T., and Grover, R., 1972, J. chem. Phys., 57, 
1259. 
[34] Sevick, E. M., Monson, P. A., and Ottino, J. M., 1988, J. chem. Phys, 88, 
1198. 
[35] Lee, S. B., and Torquato, S., 1988, J. chem. Phys, 89, 3258. 
[36] Speedy, R. J., and Reiss, H., 1991, Molec. Phys., 72, 1015. 
[37] Speedy, R. J., 1980, J. chem. Soc. Faraday, Trans. II, 76, 693. 
14 

[38] Stell, G., 1966, Lecture Notes, Polytechnic Institute of Brooklyn, Brooklyn, New 
York. 
[39] Boltzmann, L., 1964, Lectures on Gas Theory (University of California Press), trans­ 
lated by S. Brush. 
[40] Speedy, R. J., 1981, J. chem. Soc. Faraday, Trans. II, 77, 329. 
[41] Hoover, W. G., Hoover, N. E., and Hanson, K., 1979, J. chem. Phys., 70, 1837. 
[42] Sturgeon, K. S., and Stillinger, F. H., 1992, J. chem. Phys., 96, 4651. 
[43] Rintoul, M. D., and Torquato, S., 1995, Phys. Rev. E, 52, 2635. 
[44] Kerstein, A. R., 1983, J. Phys. A.: Math. Gen., 16, 3071. 
[45] Medvedev, N. N., 1994, Doklady Akademii Nauk, 337, 767, (English translation: 
1994, Doklady Physical Chemistry, 337, 157). 
[46] Lubachevsky, B. D., and Stillinger, F. H., 1990, J. statist. Phys., 60, 561. 
[47] Lubachevsky, B. D., and Stillinger, F. H., 1991, J. statist. Phys., 64, 501. 
[48] Rintoul, M. D., and Torquato, S., 1996, J. chem. Phys., 105, 9258. 
[49] Schaaf, P., and Reiss, H., 1990, J. chem. Phys., 92, 1258. 
[50] Stillinger, F. H., 1995, Comput. Mater. Sci., 4, 383. 
[51] Stillinger, F. H., 1995, Science, 267, 1935. 
[52] Debenedetti, P. G., 1996, Metastable Liquids (Princeton University Press). 
[53] Corti, D. S., Debenedetti, P. G., Sastry, S., and Stillinger, F. H., 1997, 
Phys. Rev. E, 55, 5522. 
[54] Steinhardt, P. J., Nelson, D. R., and Ronchetti, M., 1983, Phys. Rev. B, 28, 
784. 
15 

[55] Speedy, R. J., and Reiss, H., 1991, Molec. Phys., 72, 999. 
[56] Carnahan, N. F., and Starling, K. E., 1969, J. chem. Phys, 51, 635. 
[57] Sanchez, I. C., 1994, J. chem. Phys., 101, 7003. 
[58] Attard, P., 1993, J. chem. Phys, 98, 2225. 
[59] Labik, S., and Smith, W. R., 1994, Molec. Sim., 12, 23. 
16