Abstract

The pair correlation function g₂(r) provides a basic geometric descriptor for many-particle systems. It must obey two necessary conditions: (i) non-negativity for all distances r, and (ii) non-negativity of its associated structure factor S(k) for all k. Here we utilize an improved stochastic construction algorithm for particle configurations to establish conditions in which (i) and (ii) are also sufficient, i.e., g₂(r) is in fact realizable. Two types of target pair correlation functions have been investigated in one, two, and three dimensions for hard-core particles, specifically a unit step function, and a contact δ plus step pair correlation function. Results indicate that the former target function is realizable up to a terminal density set by necessary condition (ii), at which the particle core packing fraction equals 2^-d in d dimensions. Furthermore, results are consistent with the proposition that for d>1 the contact δ plus step function is realizable up to a terminal density due to condition (ii) at which the packing fraction of cores is (d+2)/2^d+1 [Torquato and Stillinger, J. Phys. Chem. B 106 (2000) 8354, 11406].

Keywords: Pair correlation function; Realizability; Hard spheres and disks

PACS: 05.20.-y; 05.10.-a; 02.70.Uu

Article Outline

1. Introduction

2. Numerical procedures

2.1. Construction algorithm

2.2. Bin width selection algorithm

2.3. Ensemble averaging of data

3. Results

3.1. Unit step function—three dimensions

3.2. Contact δ plus step function—one dimension

3.3. Contact δ plus step function—two dimensions

3.4. Contact δ plus step function—three dimensions

4. Discussion and conclusions

References

1. Introduction

The generation of particle configurations consistent with a limited amount of microstructural information (lower-order correlation functions) constitutes an example of an inverse problem. Computer simulations that accomplish this goal utilize reconstruction algorithms [1], [2], [3] and [4]. Such reconstructions can potentially identify the appropriate correlation functions that effectively characterize a class of random structures [4]. An intriguing extension of the reconstruction algorithm is the creation of particle configurations conforming to a model or hypothetical statistical correlation function [2], [3] and [5]. In this instance, the numerical protocol is termed a construction algorithm. This paper investigates the realizability of certain pair correlation functions by d-dimensional particle systems using the construction algorithm.

Applications of our methodology include the construction of three-dimensional structures using information obtained from two-dimensional images [2]. This ability has high value in a wide variety of fields, such as petroleum engineering, biology, and medicine, where often only two-dimensional images are available for analysis [4]. In addition, our techniques can be used for construction of spatial point patterns based on correlation functions employed by ecological and cosmological models. An attractive fundamental problem would entail probing the non-uniqueness of configurations specifically constructed for a determined pair correlation function [4].

In order for a given pair correlation function g₂(r) to be realizable by a statistically homogeneous many-particle system, it must obey certain conditions [6]. Firstly, as g₂(r) represents the probability of locating particles separated by the radial distance r, it cannot be negative

(1.1)

The structure factor S(k) is related to the Fourier transform of the total correlation function h(r) = g₂(r)-1 as follows [4]:

(1.2)

where k is the wave vector, k is its magnitude, and ρ is the number density of the particle system. Secondly, this structure factor is constrained to nonnegative values for all values of k

S(k)0

(1.3)

on account of its relation to mean-square fluctuations of collective density variables.

A third condition forces a lower bound on the variance associated with the number of particles within a d-dimensional spherical window of radius R [7]:

(1.4)

where θ is the fractional (non-integer) part of the average number of particle centers contained within the window N(R) = ρυ₁(R), υ₁(R) is the volume of the spherical window given by [4]

(1.5)

and Γ(x) is the gamma function. For a statistically homogeneous and isotropic system, the number variance σ²(R) is given by [8]

(1.6)

where α(r;R) is a scaled intersection volume of two identical d-dimensional spheres of radius R whose centers are separated by a distance r, and dr is a d-dimensional volume element. It should be noted that Eq. (1.4) is a condition that becomes progressively weaker for higher dimensions because the variance grows at least as rapidly as R^d-1 for sufficiently large R [8]. However, the three inequalities (1.1), (1.3), and (1.4) necessarily must be satisfied for realizability of any hypothetical pair correlation function [6]. Nonetheless, the above conditions constitute a small subset of the infinite number of conditions that may be necessary for correlation function realizability [9].

In this paper, we investigate the realizability of two distinct types of pair correlation functions, namely the unit step function, and the contact δ plus step function. Both investigated target functions are defined for statistically homogeneous and isotropic systems of identical d-dimensional hard particles of diameter D. The unit step function case [5] and [10] is defined as

(1.7)

where Θ(x) is the Heaviside step function,

(1.8)

The contact δ plus step function [6] is defined as the following:

(1.9)

Z is the mean number of contacts experienced by particles, δ(x) is the Dirac delta function, and s₁(R) is the surface area of a d-dimensional sphere

(1.10)

Fig. 1 shows graphical representations of the above target functions. The unit step pair correlation function was examined by Stillinger et al. [10] during their investigation of “iso-g₂” processes. One of our goals is to determine whether the unit step pair correlation function, Eq. (1.7), is realizable for hard-particle configurations with packing fraction φ in the range 0φφ_c, where the packing fraction is defined by [4]

(1.11)

Note that the provisional terminal density for the unit step function φ_c=2^-d [10]. The provisional terminal density has been defined as the highest possible packing fraction that can be achieved without violating the non-negativity condition on S(k). Condition (1.4) is satisfied for this choice of pair correlation function. The realizability of the unit step pair correlation function in one and two dimensions has been investigated by Crawford et al. [5]. In this paper, one of the problems that we explore is the realizability of the unit step function g₂ for three-dimensional hard-particle systems. Our approach involves making use of a construction algorithm to generate configurations in finite systems that provide the best estimate to the unit step pair correlation function.

(16K)

Fig. 1. Graphical representations of the target functions in Eqs. (1.7) and (1.9). Left panel: unit step pair correlation function. Right panel: contact δ plus step pair correlation function.

The realizability of the contact δ plus step pair correlation function was first studied by Torquato and Stillinger [6]. We attempt construction of d-dimensional hard-particle configurations with packing fractions up to and beyond a provisional terminal density φ_c. For Eq. (1.9), a provisional terminal density for dimensions d>1 is given by the formula [6]

(1.12)

For this case, we have identified the provisional terminal density as the highest possible density subject to optimization with respect to Z such that the non-negativity condition on S(k) is observed. It should be noted that Eq. (1.12) is obtained for d-dimensional systems with coordination number Z = d/2 [6]. Thus, the provisional terminal density φ_c=0.5 in two dimensions and φ_c=0.3125 in three dimensions. Application of Eq. (1.12) to one-dimensional systems formally yields a terminal density of 0.75. However, it should be noted that one-dimensional systems with packing fraction φ=0.75 violate inequality (1.4). In reality, for the choice of Z=1/2, the terminal density φ_c=0.7185 is the highest possible attained packing fraction for one-dimensional hard-particle systems. We stress that this reduction in terminal density below the provisional terminal density in Eq. (1.12) does not apply for d2.

The rest of this article is organized as follows. In Section 2, we present our construction algorithm for generation of particle configurations for the investigated target correlation functions. In addition, we introduce a bin width selection procedure for use in calculating the pair correlation function of finite-sized rigid particle systems. In Section 3, we discuss results in three dimensions for the unit step pair correlation function and results in one, two, and three dimensions for the contact δ plus step pair correlation function. Concluding remarks and discussion of future work are contained in Section 4.

2. Numerical procedures

In Section 2.1, we present a construction algorithm for the investigation of the realizability of a given target function. The construction algorithm will be slightly modified for the second investigated pair correlation function, Eq. (1.9). In Section 2.2, we discuss the procedure for calculation of the pair correlation function, including a bin width selection algorithm for hard-particle systems. Finally, in Section 2.3 we discuss our approach to data analysis.

2.1. Construction algorithm

In all of our simulations, all N system particles are contained within a d-dimensional cube of length L = 1 subject to periodic boundary conditions. We use the method of simulated annealing [11] to construct configurations of particles for a given target function. It should be noted that the pair correlation function g₂(r) is termed the radial distribution function (RDF) for the statistically homogeneous and isotropic systems that are explored in our study. We require that our results match the target RDF within the maximum radial distance selected in the binning procedure discussed below.

We begin our simulation by initializing a random configuration of particles at the packing fraction of interest as defined by Eq. (1.11). A fictitious energy is introduced as follows:

(2.1)

where is the target RDF, defined by Eq. (1.7) or (1.9), g₂(r_i) is the calculated RDF at any time step in the simulation, and the sum operates over all bin distances r_i up to a specified limit termed the sampling distance. The energy is calculated for the initial configuration of particles, then a new configuration is generated by adhering to the set of rules outlined below.

In previous work [5], the energy was computed after movement of any single particle. In contrast, we will employ a collective movement of particles to generate new configurations. In particular, sets of particles (previously defined fractions of the number of particles in the system [12]) are moved by displacement along each axis by amounts randomly and uniformly distributed in the interval [-τ,τ], where τ is the maximum step size. The energy is computed for each new configuration and the energy associated with the previous configuration is saved for comparison. It should be noted that the system RDF g₂(r_i) must only be partially recomputed after each collective movement. This time-saving measure is implemented as the system RDF is not significantly altered for each new configuration. In other words, as the binned number of particles varies for a selected number of bins after each collective movement, the binned particles and hence, the system RDF remains unchanged for the remaining bins. The bins affected by the collective movement involve those bins that have lost or gained particle pairs of separation distance that fall within the specified bin distance. The collective movement is accepted or rejected with probability P(ΔE) according to the Metropolis acceptance rule [4]

(2.2)

where ΔE is the energy difference between the present and most recent particle configuration and T is a fictitious temperature [4]. Modifying the simulated annealing algorithm to include the collective movement of particles and the partial recalculation of the system RDF reduces the running time by a factor of approximately 10.

At each temperature T, the system is thermalized until the particles have experienced multiple displacements from their previous positions. The acceptance of uphill moves with a Boltzmann factor probability enables the particle configuration to navigate the energy landscape without inevitably getting caught in local minima for which E>0. Repeated application of the collective movement of particles allows the configuration RDF the opportunity to converge to the target RDF. A cooling schedule, which governs the value and the rate of change of T, is chosen to allow the system to evolve to the desired state as quickly as possible. In our simulations, we adopt a cooling schedule in which the temperature is reduced by a factor of 0.9 for each temperature cycle. Eventually, we approach the global minimum of the energy E, Eq. (2.1). At any stage of the calculation, the configuration energy E can be viewed as a least-squares error.

The construction algorithm can be applied to an investigation of the realizability of the unit step RDF without any modifications. However, the algorithm must be adapted for construction of particle systems that match the contact δ plus step RDF. In particular, we begin with a random configuration of contact n-mers i.e., monomers, dimers, trimers, etc. at a fixed mean contact value Z. In order to preserve the contact between particles in an n-mer, we ensure that all particles in that n-mer simultaneously participate in the same collective movement and thus are moved by the same displacement and direction. In essence, the clusters move as an integrated unit during the simulation. This contrasts with an earlier study [1] in which contacting n-mers were not included in the initial conditions but formed as a result of the subsequent optimization procedure. Inclusion of n-mers at the outset of the algorithm makes the task of achieving the target pair correlation function appreciably quicker. The calculated system RDF will involve contributions between particles in different n-mers and contributions from particles that are within the same n-mer. It should be noted that particle pairs that are in direct contact within the same n-mer contribute to the Dirac delta function in Eq. (1.9), whereas all other particle pairs provide binned contributions to the Heaviside function term in the contact δ plus step pair correlation function.

At the onset of our simulations, we employ a variety of system sizes N, packing fractions φ, and particle coordination numbers Z to obtain desired d-dimensional configurations. In particular, our system sizes have varied from 100 to 500 particles while investigated packing fractions have ranged from 0.05 to densities above the calculated terminal density. Likewise, for simulations concerning the contact δ plus step pair correlation function, the selected particle contact values have fallen in the range 0.10Z0.5 in one dimension, 0.10Z1.0 in two dimensions, and 0.10Z1.5 in three dimensions. Note that there exists an infinite variety of cluster distributions consistent with a given Z value. This point will be discussed in greater detail in the upcoming sections.

2.2. Bin width selection algorithm

After each collective movement, the RDF is partially recalculated from a histogram of the average number of particle centers n(r) contained in a concentric shell of finite thickness Δr at radial distance r from an arbitrary reference particle center [4]. The radial distance r is defined as halfway between the inner radius (r-Δr/2) and the outer radius (r+Δr/2) of each shell. The shell thickness is the bin width. Let n_k(r) represent the accumulated count of particles placed in bin k for the entire system corresponding to a radial distance r. Consequently, all pairs are enumerated twice and by virtue of its definition, n_k(r) is an even integer. Then

(2.3)

and the RDF in one, two, and three dimensions is given by the following set of equations in which we have fixed L=1:

(2.4)

(2.5)

(2.6)

For both investigated pair correlation functions, we select r and Δr to ensure that the target number of pairs in each bin is an even integer such that g₂(r)=1 for r>D.

Considering all possible particle configurations consistent with Eq. (1.7) or (1.9), we sum over all particles to obtain the expected number of particle centers that lie in the annular region r-Δr/2r<r+Δr/2 around any chosen particle as follows:

(2.7)

where υ₁(x) is as defined in Eq. (1.5). Using Eq. (1.11), we rewrite Eq. (2.7) as the following expression:

(2.8)

Using an approximate bin width Δr = 0.1D, we calculate B(D,1.1D) from Eq. (2.8). In general, B(D,1.1D) is not an even integer, thus, the function ν_e[B(D,1.1D)] yields its closest even integer

(2.9)

Here, Q^* is the number of pairs required to fall in the bin of interest. We can obtain the outer radius by application of Eq. (2.8). For example, in two dimensions, the outer radius is given by

(2.10)

Then, the bin width Δr and radius r are easily obtained as

(2.11)

The above procedure is repeated for each successive bin until the radial distance is greater than the chosen sampling distance for the investigated particle configuration. Our bin width selection algorithm is consistent with g₂(r)=1 for r>D in the infinite system limit. It should be noted that the particle pairs binned by the selection algorithm excludes contact pairs in a given n-mer.

2.3. Ensemble averaging of data

In the analysis of the collected data, we depart from the approach used by Crawford et al. [5] which required that a single configuration match the target RDF within the sampling distance. Such a requirement is both unrealistic and overly stringent. As we intend to equal the target RDF in the thermodynamic limit, we employ the ergodic hypothesis i.e., the result of averaging over all realizations of an ensemble is equivalent to averaging over the volume for one realization in the infinite-volume limit [4]. This approach enables us to replace volume averaging in the limit that the volume tends to infinity with ensemble averaging. Thus, we average system RDFs of constructed configurations (of fixed φ) to obtain a match with the target RDF.

3. Results

Our results indicate that the unit step function is indeed realizable for three-dimensional hard-particle configurations below the terminal density φ_c=0.125. In addition, the realizability of the contact δ plus step RDF for hard-particle systems in one, two, and three dimensions below the provisional terminal density is supported by the results of our investigation. Provided that the system density is not too large, satisfying the non-negativity conditions on g₂(r) and S(k) suffices to ensure that both target RDFs are realizable by hard-particle configurations. However, in order for the contact δ plus step RDF to be realizable at the terminal density φ_c, the mean number of contacts Z must be equal to d/2 where d is the system dimension.

3.1. Unit step function—three dimensions

Simulations were run for three-dimensional particle systems of packing fractions up to and beyond the terminal density. In particular, simulations of packing fraction φ=0.05, 0.1, 0.125, 0.15, 0.2, and 0.3 and system size N=100 and 256 were executed by application of the algorithms discussed in Section 2. It should be noted that multiple cases were executed for each investigated packing fraction for systems of N=100.

For systems of φ0.15, we obtained a perfect match to the unit step RDF up to a sampling distance of 3 particle diameters. Fig. 2 and Table 1 display a constructed particle configuration at the terminal density and associated parameters for the simulation of the constructed configuration, respectively. As indicated in Table 1, the bin widths Δr used in all simulations of the unit step RDF were of order 0.1D.

(28K)

Fig. 2. A three-dimensional configuration of 256 particles that has been constructed for the unit step pair correlation function at the terminal density φ_c=0.125.

Table 1.

Inner radius (r-Δr/2)/D, outer radius (r+Δr/2)/D, and the number of particles n_k(r) in each bin k for r>D for construction of the unit step function for a three-dimensional hard-particle system displayed in Fig. 2

k	(r-Δr/2)/D	(r+Δr/2)/D	nk(r)	k	(r-Δr/2)/D	(r+Δr/2)/D	nk(r)
1	1	1.09952	84	11	2.00069	2.10072	322
2	1.09552	1.20024	102	12	2.10072	2.20069	354
3	1.20024	1.30048	120	13	2.20069	2.30075	388
4	1.30048	1.40071	140	14	2.30075	2.40097	424
5	1.40071	1.50121	162	15	2.40097	2.501	460
6	1.50121	1.60109	184	16	2.501	2.60098	498
7	1.60109	1.70077	208	17	2.60098	2.70099	538
8	1.70077	1.8005	234	18	2.70099	2.8011	580
9	1.8005	1.90044	262	19	2.8011	2.90107	622
10	1.90044	2.00069	292	20	2.90107	3.00098	666

Based on an earlier investigation on one- and two-dimensional systems [5], the realizability of the unit step function for φφ_c=0.125 was not unexpected. However, we were surprised to achieve the target RDF for a system of packing fraction φ=0.15. We hypothesize that the particles are able to conform to the unit step RDF for the investigated sampling distance but with some deviation beyond that distance. In fact, it should be noted that the RDF deviates uniformly and randomly from unity beyond the sampling distance for all systems that result in perfect matches to the unit step function. Upon averaging the RDFs for the constructed realizations, these deviations are greatly diminished beyond the sampling distance.

For packing fractions that are substantially above the terminal density, we are unable to construct perfect matches to the unit step RDF. At a density of φ=0.3, the constructed RDF has an initial peak at r=D that is followed by a trough (see Fig. 3). A comparison of the RDFs for constructed and equilibrium systems at φ=0.3 (above the terminal density) is displayed in Fig. 3. Note that the RDF value at contact g₂(D⁺)≈1.5 differs considerably from the contact value for the corresponding equilibrium system RDF, g₂(D⁺)≈2.3. In addition, for increasing r, the fluctuations in the constructed RDF appear damped relative to the fluctuations in the equilibrium RDF indicating the efforts of the construction algorithm to flatten the RDFs of generated systems.

(20K)

Fig. 3. Radial distribution functions for constructed particle configurations for the unit step function, and for the equilibrium configurations. All sampled configurations have packing fraction φ=0.3 in three dimensions. The constructed RDF is averaged over ten hard-particle configurations.

3.2. Contact δ plus step function—one dimension

For d=1, simulations were run for packing fraction φ=0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.71, and 0.75 at different average contact values Z=0.1, 0.34, and 0.5. In one dimension, the denominator of Eq. (2.4) is independent of r but dependent on the bin width Δr. As a result, we are able to bypass application of the bin width selection algorithm and simply select appropriate bin widths of order 0.1D. Our constructed hard-particle systems included configurations in which half the number of particles were present as dimers (Z=0.5) and others in which 34% of the particles occurred as dimers (Z=0.34).

For systems with φ0.7, we obtained particle configurations of varied mean contact numbers with RDFs that matched the contact δ plus step pair correlation function up to a sampling distance of 3.5 particle diameters. It should be noted that for φ<0.4, we derived perfect matches to the target function for an extended sampling distance of 5 particle diameters. Beyond the sampling distance, these constructed RDFs uniformly and randomly deviate from unity. An ensemble average of the constructed RDFs reveals that the deviations are greatly diminished relative to a single constructed RDF for distances greater than the sampling distance (see Fig. 4).

(25K)

Fig. 4. An illustration of the effect of ensemble averaging of constructed RDFs beyond the sampling distance. Left panel: A single constructed RDF. Right panel: An ensemble averaged RDF for 400 one-dimensional constructed configurations. All sampled configurations have been generated for the contact δ plus step pair correlation function with packing fraction φ=0.1 and mean contact value Z=0.1.

However, for φ=0.71, we were unable to achieve the target RDF for average contact values other than Z=0.5. Furthermore, particle configurations of packing fraction φ=0.75 were unrealizable for the contact δ plus step pair correlation function for any mean particle contact number. This is due to the fact that at φ=0.75 the one-dimensional system violates Eq. (1.4). Indeed, the lower bound on the number variance σ²(R) is obeyed over the range 0<φ0.7185 for Z=0.5. Our lack of success in realizing the investigated RDF for φ=0.7185 may be attributed to a violation of a yet undiscovered necessary condition. Fig. 5 compares the constructed RDFs for systems of packing fraction φ=0.71 and 0.75. We note that the attempted construction of the hard-particle system of φ=0.75 indicates an inherent robustness in our construction algorithm.

(19K)

Fig. 5. RDFs for one-dimensional configurations that have been constructed for the contact δ plus step pair correlation function. Left panel: Sampled configuration at packing fraction φ=0.71. It should be noted that we succeeded in obtaining a perfect match to the target RDF prior to ensemble-averaging our results. Right panel: Sampled configuration at packing fraction φ=0.75. Both sampled configurations have mean contact value Z=0.5.

3.3. Contact δ plus step function—two dimensions

Simulations were run for packing fraction φ=0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.4, and 0.5. For each investigated packing fraction, we varied our system size (N=100, 200, and 500) and the average contact value (Z=0.1, 0.25, 0.5, 0.6, and 1.0). During each simulation, the bin width selection algorithm was employed to select bin widths of order 0.1D. The constructed configurations consisted of a distribution of monomers, dimers, and trimers. Trimers occurred in both extended (linear) and contracted (triangular) forms. The distribution of monomers, dimers, and trimers is dependent on the desired mean contact number. For example, for N=200 and Z=1.0, two possible configurations are one in which all particles appear as dimers and another in which 7.5% of the particles occur as trimers, 90% occur as dimers, and the remaining particles are monomers. As indicated, a specific distribution of n-mers need not be used in simulations to obtain agreement with the investigated pair correlation function.

For systems of packing fractions below 0.2, we obtained a perfect match to the target RDF up to a sampling distance of 5 particle diameters. These perfect results were achieved for all tested values of Z. Similarly, we achieved perfect results for systems of φ=0.25, 0.3, and 0.4 up to a sampling distance of 3 diameters for varying mean contact values. At the terminal density φ_c=0.5, we achieve a perfect match to the contact δ plus step function up to a sampling distance of 3 diameters for Z=1.0. Fig. 6 displays a constructed particle configuration at the terminal density. It should be noted that the RDF uniformly and randomly deviates from unity beyond the sampling distance for all systems that result in perfect matches to the contact δ plus step function.

(16K)

Fig. 6. A two-dimensional configuration of 500 particles that has been constructed for the contact δ plus step pair correlation function. The configuration consists of only dimers at the terminal density φ_c=0.5 with an average contact value Z=1.0.

3.4. Contact δ plus step function—three dimensions

Simulations were run for packing fraction φ=0.05, 0.1, 0.125, 0.15, 0.2, and 0.3125. For each investigated packing fraction, we varied our system size (N=100 and 256) and the mean contact value (Z=0.1, 0.3, 0.6, 0.7, and 1.5). During each simulation, the bin width selection algorithm was employed to select bin widths of order 0.1D. Initial configurations for the construction algorithm consisted of a distribution of monomers, dimers, trimers, and tetramers. As before, the distribution of contact n-mers was dependent on the desired mean particle contact value. Both trimers and tetramers were initialized as combinations of extended (linear) and compact forms. For example, a configuration in which all particles appear in extended tetramers (two particles per tetramer have two contacts and the other two particles have a single contact) has a mean contact value of Z=1.5. The above configuration was used in our investigation of systems of Z=1.5 below and at the terminal density.

For systems of φ=0.05, 0.1, 0.125, 0.15, and 0.2, we obtained a perfect match to the target RDF up to a sampling distance of 3 particle diameters. These perfect results were achieved for all tested values of Z. At the terminal density φ_c=0.3125, we achieve a perfect match to the contact δ plus step function up to a sampling distance of 2.75 diameters for an average contact value Z=1.5. Fig. 7 displays a constructed configuration at the terminal density. Beyond the sampling distance, the constructed system RDF randomly deviates from unity in a uniform fashion for all systems that result in perfect matches to the contact δ plus step function.

(27K)

Fig. 7. A three-dimensional configuration of 256 particles that has been constructed for the contact δ plus step pair correlation function at the terminal density φ_c=0.3125 and mean contact value Z=1.5. All particles in the configuration appear as tetramers (yellow).

4. Discussion and conclusions

During the course of our investigation, we have modified the generalized simulated annealing algorithm [11] to include two algorithm-speeding measures: the collective movement of particles and the partial recalculation of the system RDF. Both time-saving implementations have added a substantial utility to our construction algorithm on account of increased likelihood of obtaining converged results within a reasonable time. This added utility is also evident in our investigation of the near approach of the pair correlation function to realizability for systems above the terminal density. In particular, the algorithm finds a close compromise when an exact solution is not realizable. The right panel of Fig. 5 displays the system RDF for a configuration above the terminal density. Clearly, the program attempts to match the contact δ plus step RDF as closely as possible. In addition, we note that the results obtained for all investigated systems are consistent with the necessary Yamada condition defined in Eq. (1.4) indicating the inherent robustness of our construction algorithm.

The unit step RDF was achieved for three-dimensional systems below or at the terminal density φ_c=0.125. This is consistent with previous results in one and two dimensions [5] for systems below the terminal density. In contrast, no simulation studies had previously been executed for the contact δ plus step pair correlation function in any dimension. We obtained perfect matches to the contact δ plus step RDF for one-dimensional systems with packing fractions φ<0.7185. However, it should be noted that at densities close to the terminal density, only systems with a mean particle contact value Z=0.5 yielded matches to the target RDF. Similarly, the contact δ plus step function was achieved for all investigated values of Z below the provisional terminal density defined in Eq. (1.12) for two- and three-dimensional hard-particle systems. At the terminal density, matches to the target RDF were only obtained for systems in two and three dimensions with Z=1.0 and 1.5, respectively.

It is useful to remark on the choice of the cluster distribution for realizability simulations of the contact δ plus step pair correlation function at the terminal density. Our results indicate that matches to the target RDF are only achieved for d-dimensional systems with mean contact value Z=d/2 at densities close to the terminal density. However, there are an infinite number of different cluster distributions that would be consistent with Z=d/2. It is remarkable that our selected cluster distributions result in perfect matches to the pair correlation function, and thus indicates an insensitivity or non-uniqueness to the cluster distribution.

We note that it has been shown that the two- and three-dimensional systems constructed for the contact δ plus step RDF as well as the three-dimensional systems generated by the unit step RDF are hyperuniform systems at the terminal density [8]. Hyperuniformity is concerned with behavior of density fluctuations at large length scales. For a very large class of particle systems, the number variance within a specified window (see Eq. (1.6)), asymptotically approaches proportionality to the window volume for large window size. In contrast, for hyperuniform systems the number variance grows asymptotically as the surface area of the d-dimensional window. Thus, infinite-wavelength density fluctuations vanish for hyperuniform systems.

More elaborate sets of target functions defined for amorphous sphere packings qualify for an analysis such as presented in this paper. Particularly, the deviations of the target functions from unity for distances beyond contact are evidently important in determining associated properties i.e., jamming density. Applications of the methods proposed in this paper cannot fail to produce a deeper understanding of the statistical geometry involved in jammed sphere packings. It should be noted that the hyperuniformity of such sphere packings is dependent on the construction protocol [6].

As a final subject, we note that pair correlation functions are of central importance in the thermodynamics of polymer melts, solutions, and glasses because they not only contain information on the global effects of intermolecular (polymer–polymer) interaction but also give detailed information on the consequences of the intermolecular potential of interaction. In addition, calculation of the intermolecular RDF for polymeric liquids assists in the analysis of polymer equilibrium structure and provides a direct route to the calculation of all thermodynamic properties, including an equation of state improving upon those of the solid-like lattice or cell models [13]. The reference interaction site model (RISM) pioneered by Chandler et al. [14] in their investigation of the local structure of molecular liquids has served as an underlying framework in investigations on polymer systems. In particular, variations on the RISM integral equation approach have been used in the analysis of equilibrium theory of polymer liquids in several studies of linear polymer chains [15] and [16].

A desirable expansion of our present analysis involves investigation of the constraints concerning pair correlation function realizability for polymer chains under very general circumstances. In the intended investigation, one would employ a hypothetical polymer model of chains of contacting monomers in either branched or unbranched forms. We stress that such an implementation is a model in the venerable tradition of simplifying details of chemical bonding structure. In closing, we note that an extension to the study of the realizability of the three-body distribution function g₃ presents an attractive opportunity for future work. For example, one could study the realizability of g₃'s that satisfy the Kirkwood superposition approximation [17].

References

[1] M.D. Rintoul and S. Torquato, J. Colloid Inter. Sci. 186 (1997), p. 467. Abstract | PDF (348 K)

[2] C.L.Y. Yeong and S. Torquato, Phys. Rev. E 57 (1998), p. 495. Abstract-INSPEC   | MathSciNet | Full Text via CrossRef | APS full text

[3] D. Cule and S. Torquato, J. Appl. Phys. 86 (1999), p. 3428. Abstract-INSPEC   | OJPS full text | Full Text via CrossRef

[4] S. Torquato, Random Heterogeneous Materials: Microstructure and Macroscopic Properties, Springer, New York (2002).

[5] J. Crawford, S. Torquato and F.H. Stillinger, J. Chem. Phys. 119 (2003), p. 7065. Abstract-INSPEC | Abstract-Compendex   | Full Text via CrossRef

[6] S. Torquato, F.H. Stillinger, J. Phys. Chem. B 106 (2002) 8354, 11406.

[7] M. Yamada, Prog. Theor. Phys. 25 (1961), p. 579 Note that the inequality (1.4) applies to windows of arbitrary shape. MathSciNet | Full Text via CrossRef

[8] S. Torquato and F.H. Stillinger, Phys. Rev. E 68 (2003), p. 041113. Full Text via CrossRef

[9] O. Costin and J.L. Lebowitz, J. Phys. Chem. B 108 (2004), p. 19614. Abstract-Compendex | Abstract-INSPEC   | Full Text via CrossRef

[10] F.H. Stillinger, S. Torquato, J.M. Eroles and T.M. Truskett, J. Phys. Chem. B 105 (2001), p. 6592. Abstract-Compendex | Abstract-INSPEC   | Full Text via CrossRef

[11] P.J.M. van Laarhoven and E.H.L. Aarts, Simulated Annealing: Theory and Applications, D. Reidel Publishing Co., Holland (1987).

[12] During each simulation, a fixed fraction of the number of particles in the system was displaced by each collective movement. However, this fraction was dependent on the system of interest. In some cases, 10 percent of the particles were displaced per collective movement; in other cases, twenty percent of the particles were displaced during each such cycle.

[13] K.S. Schweizer and J.G. Curro, Phys. Rev. Lett. 58 (1987), p. 246. Abstract-INSPEC   | Full Text via CrossRef

[14] D. Chandler In: E.W. Montroll and J.L. Lebowitz, Editors, Studies in Statistical Mechanics VIII, North-Holland, Amsterdam (1982), p. 275 and references cited therein.

[15] J.G. Curro and K.S. Schweizer, J. Chem. Phys. 87 (1987), p. 1842. Abstract-INSPEC   | Full Text via CrossRef

[16] J.P. Donley, J.J. Rajasekaran and A.J. Liu, J. Chem. Phys. 109 (1998), p. 10499. Abstract-INSPEC   | OJPS full text | Full Text via CrossRef

[17] J.G. Kirkwood, J. Chem. Phys. 3 (1935), p. 300. Full Text via CrossRef

Corresponding author. Princeton University, Department of Mechanical Engineering, Princeton, NJ 08544, USA.
¹ Supported by the U.S. Department of Energy Computational Science Graduate Fellowship.
² Supported by the Office of Basic Energy Sciences, DOE, under Grant no. DE-FG02-04ER46108.