Cluster pair correlation function of simple fluids: energetic connectivity criteria

We consider the clustering of Lennard-Jones particles by using an energetic connectivity criterion proposed long ago by T.L. Hill [J. Chem. Phys. 32, 617 (1955)] for the bond between pairs of particles. The criterion establishes that two particles are bonded (directly connected) if their relative kinetic energy is less than minus their relative potential energy. Thus, in general, it depends on the direction as well as on the magnitude of the velocities and positions of the particles. An integral equation for the pair connectedness function, proposed by two of the authors [Phys Rev. E 61, R6067 (2000)], is solved for this criterion and the results are compared with those obtained from molecular dynamics simulations and from a connectedness Percus-Yevick like integral equation for a velocity-averaged version of Hill's energetic criterion.


I. INTRODUCTION
The concepts of clustering and percolation have been widely used in order to explain several phenomena in very diverse areas including Physics, Chemistry, Biology, Geology, Sociology and Economics. In particular, with reference to chemical-physics, phenomena such as nucleation, 1 hydrogen bonding, 2 insulator-conductor, sol-gel and glass transitions 3,4,5,6,7,8,9 as well as bridging in granular materials 10 are currently studied from this point of view. In all these cases, the system under study can be thought of as a collection of individuals (atoms, molecules, grains, etc.) that, with generality, we call particles. Most of the efforts have been based on lattice representations of the systems. The relative simplicity of lattice models allows for a wide variety of treatments, which extend from almost heuristic 11 to quite rigorous. 12 Whatever the treatment is, the concept of connectivity between the particles plays an important role.
Sometimes, however, a continuous description-where particles can occupy any point in a continuum phase space-is needed to reach a more realistic picture of the phenomena under consideration. For this context, the concept of connectivity has been generalized and adapted to describe clustering and percolation in continuum systems. The main ideas have been established in the pioneering works of Hill 13 and Coniglio et al. 14 Hill considers a partition of the whole system into subsystems of particles (the clusters) that satisfy some linking properties. The concept of cluster is thus directly related to the idea of bonded pairs. A bonded pair is a set of two particles that are linked by some direct mechanism. A cluster is then defined as a set of particles such that any pair of particles in the set is connected through a path of bonded pairs. We call these clusters "chemical clusters" to distinguish them from the non-pair-bonded clusters we have introduced in a previous work 15 -note, however, that this does not mean that clusters are necessarily formed through a true chemical bonding. A system is said to be in a percolated configuration if it contains a cluster that spans the system volume.
From Hill's theory, we see that a connectivity criterion is needed in order to decide whether two particles are bonded or not. This connectivity criterion has to be defined in accordance with the phenomenon under study. 4,16,17 In the search for stable atomic clusters, which mark the onset of a phase transition in a monatomic gas, Hill proposed a simple energetic criterion (HE): two particles are bonded if their relative kinetic energy is less than the negative of their relative potential energy. 13 In principle, this criterion takes into account the relative positions and velocities of the relevant pair of particles. For molecular fluids, instead of just their distance, the potential energy could in general depends on the direction and magnitude of the vector position of each of the two involved particles and their relative orientations.
From a theoretical point of view, a criterion that involves the velocity of the particles prevents the straightforward integration of the momenta in the partition function, which is the great advantage of classical statistical mechanics. To avoid this obstacle, Hill 13 himself has proposed a velocity-averaged (VA) version of his criterion giving effective potentials between bound and unbound particles.
The VA and the complete HE criteria have been used in computer simulations as well as in integral equations studies. For the VA criterion only the particle positions come into account, so it is suitable to both, Monte Carlo (MC) and molecular dynamics (MD) calculations. With respect to the integral equations approach, Coniglio et al. 14 have obtained a connectedness Ornstein-Zernike (OZ) relationship for the pair connectedness function g † (r 1 , r 2 ) (see also a review from Stell). 18 This function is proportional to the joint probability density of finding two particles belonging to the same cluster and at positions r 1 and r 2 , respectively. Therefore, by integrating g † (r 1 , r 2 ), the mean cluster size S and the percolation density ρ p -i.e. the value of ρ for which S(ρ) diverges-can be obtained. Since Coniglio's theory deals only with the positions of the particles, the HE criterion cannot be implemented. Instead, the VA criterion was used by Coniglio et al. 14 to analytically calculate the percolation loci, for a potential made up of a hard core plus an attractive interaction, in a crude mean-field approximation.
It is worth mentioning that most of the theoretical studies 19,20,21,22 on connectivity and percolation in continuum systems based on Coniglio's type equations were focused in the rather simple Stillinger's connectivity criterion. 23 This criterion states that two particles are bonded if they are separated by a distance shorter than a given connectivity distance d. In this case, d is an ad hoc parameter, which must be chosen on physical grounds. Although this criterion might be sensible in the study of certain insulator-conductor transitions, it is unrealistic regarding clustering in saturated vapors.
A general theory which is appropriate for bonding criteria involving both, the momenta and positions of a pair of particles, has been developed by two of us. 24 The main object in our theory is the pair connectedness function g † (r 1 , r 2 , p 1 , p 2 ) which is proportional to the joint probability density of finding two particles at positions r 1 and r 2 with momenta p 1 and p 2 , respectively, and belonging to the same cluster. This function verifies also an OZ like relationship. In a previous paper 25 (thereafter denoted as I) we applied our general theory to study the complete HE criterion for the same model fluid considered by Coniglio et al. 14 under the VA criterion. We also used the same simple closure relation proposed by Coniglio et al. More recently, we have reported 26 the solution of our generalized connectedness OZ type relation for a Lennard-Jones fluid closed with a connectedness Percus-Yevick (PY) condition. We implemented a connectivity criterion which generalizes Stillinger's criterion in that a life time τ for the bonds is required. 15 In Ref. I we have also performed MD simulations of the Lennard-Jones fluid and have used both criteria (HE and VA) to define the clusters. We concluded that the VA criterion strongly overestimates percolation densities. We will partially revise these results here and will discuss some subtleties related to the identification of percolating clusters. Notice that MD simulations are convenient when the HE criterion is used to identify clusters since MC algorithms do not provide per se the velocities of the particles. 27 It should be mentioned that the HE criterion has been considered in MD studies of small clusters and the critical percolation behavior of Lennard-Jones fluids by several authors. 28,29,30 It has been suggested that the percolation line-the line that separates the temperature-density phase diagram into percolated and non-percolated states-might be experimentally observable. 30,31 Moreover, cluster analysis based on this criterion seems to be useful in locating the gas-liquid coexistence curve. 30 The main purpose of this work is to apply the generalized connectedness OZ type relationship closed with a connectedness PY condition for the Lennard-Jones fluid in order to obtain the pair connectedness function g † HE (r 1 , r 2 , p 1 , p 2 ) for HE clusters and thus the related cluster pair correlation function: g † HE (r 1 , r 2 ) = ρ(r 1 , p 1 )ρ(r 2 , p 2 )g † HE (r 1 , r 2 , p 1 , p 2 )dp 1 dp 2 , where ρ(r 1 , p 1 ) is N times the probability density of finding a particle at the phase space configuration (r 1 , p 1 ). We compare g † HE (r 1 , r 2 ) so obtained with the function g † VA (r 1 , r 2 ) for the VA criterion calculated using the integral equation that results when the OZ type relationship of Coniglio et al. 14 is closed with a PY like conditions. Both functions-g † HE (r 1 , r 2 ) and g † VA (r 1 , r 2 )-are compared with the corresponding curves given by MD simulations.
The paper is organized as follow. In Sec. II we present the model system and the two connectivity criteria, i.e. HE and VA, we will work on. Also we use this section to discuss some aspects about the MD simulations. The continuum clustering theories suitable for each criteria will be sketched in Sec. III. There, we briefly describe the integral equation for g † HE (r 1 , r 2 , p 1 , p 2 ) and its solution following Lado's orthogonal polynomials method. 26,32 Finally in Sec. IV we compare the theoretical results with those obtained from simulations. We then summarize and give the conclusions.

II. MODEL SYSTEM AND ENERGETIC CRITERIA
We consider a system of N particles whose configurations are given by their positions and momenta (r i , p i ) (i = 1, ..., N). The canonical (NV T ) ensemble will be used throughout. We assume that particles interact via the Lennard-Jones pair potential where The clustering criteria are expressed in terms of the bond conditional probability density P (r i,j , p i,j ), say the probability density that two particles i and j are bonded under the condition that their positions and momenta are (r i , p i ) and (r j , p j ), respectively.

A. HE criterion
The original Hill's criterion (HE) identifies clusters by defining: 13 with p i,j the relative momentum: p i,j = p j − p i . A maximum connectivity distance d has been added to the criterion in order to avoid unrealistic bonding.

B. VA criterion
By integrating the relative momenta weighted by the Maxwell distribution in the region where the relative kinetic energy is lesser than minus the pair potential, the momenta are eliminated and the VA bond conditional probability density is obtained: 13,14 where Γ[a] is the gamma function and γ[a, x] is the incomplete gamma function.

C. Molecular dynamics calculations
We consider a system of Lennard-Jones particles in a cubic box with periodic boundary conditions in the NV T ensemble and use a leap-frog algorithm with velocity correction. 33 The time step is chosen as ∆t * = ∆tσ −1 k B T /(εm) = 0.01. Quantities are averaged over 10 3 configurations chosen every 100 ∆t after equilibration. A cut off distance equal to 2.5σ was used in the pair potential.
In the VA case, for each pair of particles that satisfies v(r i,j ) ≤ 0 and r i,j ≤ d, we generate a random number z, between 0 and 1. If [3/2] we consider that the particles form a bonded-pair, otherwise we do not. Note that this criterion can also be used in MC simulations because it does not require information on the particle velocities. To identify the clusters from the list of bonded-pairs we use the Stoddard's algorithm. 33,34 A system is said to be in a percolated state if a cluster that spans the replicas is present 50 percent of the time. 35 It is known that this criterion yields results that are marginally affected by finite size effects. 36 Then, a percolation transition curve, which separates the percolated from the non-percolated states of the system, can be drawn above the coexistence curve in the T − ρ phase diagram.
In Fig. 1, the percolation loci for HE and VA connectivity criteria are presented for d = 3σ. These simulations were performed with N = 1372 particles. The gas-liquid coexistence curve obtained by Panagiotopoulos 37 and the MC liquid-solid coexistence curve of Hansen and Verlet 38 are also shown. The MD results for the HE criterion are similar to those obtained by Campi et al. 30 These authors consider that the system is on the percolation line if the second moment of the cluster size distribution that excludes the largest cluster n ′ (s) reaches its maximum. In Fig. 2 we show the dependence of the calculated percolation density with system size and connectivity distance d. Extrapolation to infinite systems can be obtained by fitting a straight-line to a plot of ρ p versus L −1/ν . 35 We have used the universal value of ν = 0.88 ± 0.02 reported by Gaunt and Sykes 39 for three-dimensional systems. The error due to finite size effects in the calculated value of ρ p for the 1372-particle system is of 1.0%. This correction is smaller than the size of our symbols in Fig. 1. Also from Fig. 2, we see that the effect of the connectivity distance d on ρ p is negligible for d > 2.5.
As we can see from Fig. 1, the VA criterion is a very good approximation to the full HE criterion as far as the percolation loci is concerned. In Ref. I, we reported a VA percolation line that was located at much larger densities, and concluded that the approximation was rather poor. The revised results reported here show that this is not the case. The source of the error in Ref. I comes from the way percolating clusters are detected according to the Seaton-Glandt prescription. 35 All clusters in a given configuration are first identified by the Stoddard's algorithm, then each cluster is analyzed separately to determine if its replicas are connected with one another. Since the VA criterion implies the use of random numbers to decide whether two particles are connected, the second step where each separated cluster is analyzed for percolation needs to reuse the same random numbers generated when it was first identified. This subtlety was overseen in Ref. I, which led to the incorrect identification of actual percolating clusters as disconnected replicas.

III. CLUSTER PAIR CORRELATION FUNCTIONS
In the remainder of the paper we restrict our attention to the cluster pair correlations for the HE and VA energetic criteria. We calculate them by using the above mentioned integral equations and MD simulations. Thus, this section will be devoted to pose the integral equations for the cluster correlation functions g † HE (r 1 , r 2 ) and g † VA (r 1 , r 2 ) and briefly discuss their solutions.

A. VA criterion
In order to study clustering in a system composed of N classical particles interacting via a pair potential v(r 1 , r 2 ), Hill separated the Boltzmann factor e(r 1 , r 2 ) = exp[−βv(r 1 , r 2 )], into bonded ( †) and unbounded ( * ) terms: 13 e(r 1 , r 2 ) = e † (r 1 , r 2 ) + e * (r 1 , r 2 ). As usual β = 1/k B T , being k B the Boltzmann constant. Since e † (r 1 , r 2 ) represents the basic probability density of finding two particles bonded and at positions r 1 and r 2 , this separation yields a diagrammatic expansion for the partition function in terms of "chemical" clusters. We express Hill's separation as follows where P (r 1 , r 2 ) = P VA (r 1,2 ) is given by Eq. (4) in the case of the VA energetic criterion. Fugacity and density expansions have been found, within Hill's formalism, by Coniglio and co-workers 5 for the pair connectedness function g † (r 1 , r 2 ) ≡ g † VA (r 1 , r 2 ). As it was already mentioned, this function is proportional to the joint probability density of finding two particles belonging to the same cluster and at positions r 1 and r 2 , respectively. Moreover, by collecting nodal and non-nodal diagrams in these expansions an OZ-type relationship is obtained The function c † (r 1 , r 2 ) is the direct pair connectedness function. By posing a closure relation, an integral equation for g † (r 1 , r 2 ) is obtained. Here, we use the more reliable closure available, i.e. the PY-like relation 5 In Eq. (8), f * (r 1,2 ) = e * (r 1,2 ) − 1 = exp[−βv(r 1,2 )][1 − P VA (r 1,2 )] − 1 is the unbound Mayer function and g(r 1,2 ) is the thermal pair distribution function (PDF). In order to solve the integral equation given by Eqs. (7) and (8), for the Lennard-Jones potential we have implemented Labik's numerical algorithm. 40

The integral equation
We summarize here the basic theory that we have developed 24 to describe the clustering and percolation for clusters whose bond definition depends on the positions and momenta of the two particles under consideration.
For a system of N classical particles that interact through a pair potential v(r i , r j ), we define a density correlation function ρ(r 1 , r 2 , p 1 , p 2 ) that is N(N − 1) times the probability density of finding two particles at the phase space configurations (r 1 , p 1 ) and (r 2 , p 2 ) respectively: Here h is Planck's constant and Q(N, V, T ) the canonical partition function of the system. Then, in the same spirit of Hill and Coniglio et al., 5,13 we separate exp[−βv(r i , r j )] into connecting and blocking parts, Here f † (r i , r j , p i , p j ) represents the basic probability density that two particles in configuration (r i , r j , p i , p j ) are bonded. We will sometimes use the shorthand notation where the sum is carried out over all possible arrangements of products of functions f † i,j and f * k,l . We note that the functions f † i,j and f * i,j can depend on the momenta as well as on the positions of the two particles, but the sum of f † i,j and f * i,j must be momentum independent in order to conform to Eq. (10). Except for this last condition, the functions f † i,j and f * i,j are otherwise arbitrary for thermodynamic purposes. Of course, we choose them in such a way that the desired definition of bonded particles for HE clusters is achieved, i.e., where P (r i,j , p i,j ) = P HE (r i,j , p i,j ) is given in Eq. (3). Each term in the integrand of Eq. (11) can be represented as a diagram consisting of two white e 1 and e 2 points, N − 2 black e i points and some f † i,j and f * i,j connections except between the white points. Here we take e i ≡ exp[−β p 2 i 2m ]. White points are not integrated over whereas black points are integrated over both their positions and momenta. All the machinery normally used to handle standard diagrams in classical liquid theory 41 can now be extended to treat these new type of diagrams. By following Coniglio's recipe to separate connecting and blocking parts in the PDF, g(r 1 , r 2 ) = g † (r 1 , r 2 , p 1 , p 2 ) + g * (r 1 , r 2 , p 1 , p 2 ), we obtain an OZ-like integral equation for g † (r 1 , r 2 , p 1 , p 2 ), Here ρ(r 1 , p 1 )ρ(r 2 , p 2 )g † (r 1 , r 2 , p 1 , p 2 ) is N(N − 1) times the joint probability density of finding two particles at positions r 1 and r 2 with momenta p 1 and p 2 , respectively, and belonging to the same cluster, where the bonding criterion is given by Eqs. (12), (13) and (3), while The function c † (r 1 , r 2 , p 1 , p 2 ) denotes the sum of all the non-nodal diagrams in the diagrammatic expansion of g † (r 1 , r 2 , p 1 , p 2 ). We recall here that a nodal diagram contains at least one black point through which all paths between the two white points pass. For a homogeneous system, we have To obtain a closed integral equation with Eq. (14) or Eq. (16), we need a closure relation between g † (r 1 , r 2 , p 1 , p 2 ) and c † (r 1 , r 2 , p 1 , p 2 ). Here we will use the PY approximation g(r 1 , r 2 ) exp[βv(r 1 , r 2 )] = 1 + N(r 1 , r 2 ), where the function N(r 1 , r 2 ) is the sum of the nodal diagrams in the expansion of g(r 1 , r 2 ). Separation into connecting and blocking parts, g(r 1 , r 2 ) = g † (r 1 , r 2 , p 1 , p 2 ) + g * (r 1 , r 2 , p 1 , p 2 ) and N(r 1 , r 2 ) = N † (r 1 , r 2 , p 1 , p 2 ) + N * (r 1 , r 2 , p 1 , p 2 ), yields or, for a homogeneous system, Equation (14) joined with Eq. (17) or Eq. (16) joined with Eq. (18) give a closed set of equations for g † (r 1 , r 2 , p 1 , p 2 ). From the function g † HE (r 1 , r 2 , p 1 , p 2 ) ≡ g † (r 1 , r 2 , p 1 , p 2 ) we define the pair correlation function for energetic clusters g † HE (r 1 , r 2 ) according to Eq. (1).

Solution of the integral equation
Equivalence with an integral equation for polarizable fluids Our problem consists in solving Eq. (16) for g † (r 12 , p 1 , p 2 ) closed by the connectedness PY relation (18) with f † (r i , r j , p i , p j ) and f * (r i , r j , p i , p j ) given by Eqs. (12) and (13). In the closure relation (18), g(r 12 ) is the thermal PDF of the system. In this work we take g(r 12 ) from the solution of the thermal OZ equation in the PY approximation. 41 An equation mathematically equivalent to Eq. (16) has been previously solved by Lado 32 in the study of nonpolar polarizable molecules. Explicitly, the equation considered there, which is a generalized OZ equation, relates the fluid total correlation function (TCF) h(r 12 , p 1 , p 2 ) = g(r 12 , p 1 , p 2 ) − 1 (with g(r 12 , p 1 , p 2 ) the PDF) and the direct correlation function (DCF) c(r 12 , p 1 , p 2 ), where p i denotes the instantaneous dipolar moment induced on molecule i by the remaining molecules of the system. The function f (p) gives the instantaneous dipolar moment thermal distribution which is assumed to have a Gaussian form where α is the effective polarizability of the molecules. We observe that Eqs. (16) and (19) are the same if we identify h with g † , c with c † , the induced dipolar moment p i with the kinetic momentum p i and the polarizability α with the particle mass m. There are, however, some differences between the connectivity problem and the polarizable-molecule problem. The form of f (p) does not need to be Gaussian for polarizable molecules; moreover, f (p) is coupled to the TCF. Therefore, the value of the effective polarizability α depends on the density and temperature of the system. In the connectivity problem, however, the equivalent of f (p), ρ(r, p)/ρ, is intrinsically Gaussian and independent of the thermodynamic macrostate of the system.
Another difference between the connectivity problem here and the problem described by Lado is that our closure relation must be complemented with the condition given by Eqs. (12) and (13). In addition, the closures are different. Here we consider the connectedness version of PY whereas an almost exact relation between DCF and TCF (van Leeuwen-Groeneveld-De Boer 42 exact relation with approximate bridge function) is used by Lado. 32 Nevertheless, these differences do not affect the general method of solution developed by Lado and we can apply the same principle of expansions in orthogonal functions.
Thus, following Lado, 26,32 we start by reassigning the unknown function to be the indirect correlation function γ † (r 12 , p 1 , p 2 ) = g † (r 12 , p 1 , p 2 ) − c † (r 12 , p 1 , rather than g † (r 12 , p 1 , p 2 ), and rewriting Eq. (16) in Fourier representation, The closure given by the PY relation [Eq. (18)] together with the conditions (12), (13) and (3) yield The connectivity part of the PDF is then computed from γ † as g † (r 12 , p 1 , p 2 ) = g(r 12 ) The Fourier transform in Eq. (21) and its inverse are defined as The standard method for solving Eqs. (21) and (22) is to explicitly break out the angular dependence of all functions in the form of expansions in spherical harmonics. 43 Expansion of the pair functions in orthogonal polynomials The essential point in the integral equation solution method 32 is the expansion of all the pair functions, like γ † (r 12 , p 1 , p 2 ), in terms of orthogonal polynomials. First we expand where ω 1 and ω 2 are the directions of the momenta p 1 and p 2 , m = −m, and m = −l, −l + 1, ..., l. In this and similar expressions, the vector r 12 has been implicitly chosen as the z direction in the specification of the Euler angles ω = (θ, φ). The spherical harmonics satisfy the orthogonality condition so that the coefficients of the expansion (26) are immediately obtainable as Similarly, we can break out the kinetic momentum dependence in the form of expansions in polynomials of p, γ † l 1 l 2 m (r, p 1 , p 2 ) = n 1 ,n 2 γ † n 1 n 2 l 1 l 2 m (r) Q n 1 l 1 (p 1 ) Q n 2 l 2 (p 2 ) , which are constructed to be orthogonal with Gaussian weight function namely, The coefficients of the expansion are then again obtainable by quadratures, Given the Gaussian form of the weight function f (p), the associated polynomials are 44 where L b n (t) are the associated Laguerre polynomials 45 and Γ (z) is the gamma function. Accordingly, all the functions in r-space are expanded in the form γ † (r, p 1 , p 2 ) = 4π where the z axis is along r and the summation indices satisfy the constraints n = 0, 1, 2, ..., l = n, n − 2, n − 4, ..., 1 or 0, (35) m = 0, ±1, ±2, ..., ±l.

1/2
Np (t), k for x k = cos θ k , the kth root of P Np (x), and j for y j = cos φ j , the jth root of T Np (y), where L 1/2 Np (t), P Np (x) , and T Np (y) are the associated Laguerre, Legendre, and Chebyshev polynomials, respectively, all of order N p ; here the associated Legendre functions P lm (x) are normalized to 2. The w are the corresponding Gaussian weights,

IV. RESULTS
Firstly, as a complement to Fig. 1, we show in Fig. 3 the cluster pair correlation functions g † HE (r 12 ) and g † VA (r 12 ) obtained from MD using N = 4000 particles for (ρ * ,T * ) = (0.24,1.4), (ρ * ,T * ) = (0.42,1.4) and (ρ * ,T * ) = (0.429,1.4) where ρ * = ρσ 3 and T * = k B T /ε. The two last points correspond, respectively, to the VA and HE percolation loci for T * = 1.4. The main peak in the cluster correlation functions is higher for the VA criterion than for the HE criterion. This implies that there is a larger tendency to consider as directly connected two neighbor particles by this criterion. However, as it can be appreciated in Fig. 3, g † VA (r 12 ) falls faster than g † HE (r 12 ) for larger r. This contrasting behavior can be understood by analyzing the cluster size distribution function n(s), which gives the number of clusters in the system consisting of s particles. Figure 4 shows n(s) for the same three state points as Fig. 3. We can see here that the VA criterion identifies a larger amount of clusters than the HE criterion up to a certain size-which depends on the density. However the HE criterion always identifies some clusters larger than the largest clusters identified by the VA criterion. This means that the cluster correlation function is more long ranged for the HE clusters.
Figs. 5 and 6 show the theoretical cluster correlation functions g † HE (r 12 ) and g † VA (r 12 ), calculated by solving the corresponding integral equations [Eqs. (16), (18) and (7), (8), respectively] following the methods indicated in the previous section, together with the corresponding MD simulation results. We show the curves at temperatures T * = 1.4 and 3.0 and densities ρ * = 0.24 and 0.55, respectively. This values are rather far from the percolation loci since we have been unable so far to obtain convergence of the numerical algorithms at higher densities when the HE criterion is considered.
With generality, Figs. 5 and 6 say that the theoretical results follow quite well the trends of the corresponding simulations for each criterion.

V. CONCLUSIONS
We have shown that, contrary to our previous results (Ref. I), the VA energetic criterion is, in general, a good approximation to the full HE criterion for estimating the percolation loci in a Lennard-Jones fluid. However, the cluster correlation functions are somewhat different in the VA case. We have obtained the cluster pair correlation functions for both energetic criteria through the numerical integration of connectedness OZ integral equations. In particular, we have used a generalization of the integral equations that allows the implementation of the HE criterion. The theoretical results agree rather well with the simulations.  Open diamonds correspond to the gas-liquid coexistence. 37 Open squares correspond to the fluid-solid coexistence. 38 The percolation loci for the HE criterion (open circles) and for the VA criterion (open triangles) are compared. Solid circles correspond to the percolation loci for the HE criterion from Campi et al. 30 Lines are only to guide the eye. Figure 2. Percolation density for the HE criterion as a function of L −1/ν . L is the simulation box length and ν = 0.88. 39 The largest system corresponds to 4000 particles. The system size used for the calculation of the percolation curves in Fig. 1 is indicated by an arrow. The inset shows the percolation threshold as a function of d for the 1372-particle system. The arrow shows the value of d used in the rest of the paper.