Calculation of pair and triplet correlation functions for a Lennard - Jones fluid at and

On the basis of Monte Carlo simulations, configurations and pair and triplet correlation functions for a Lennard - Jones fluid are reported over a large range of densities and temperatures (0.002 ≤ ρ ∗ ≤ 1.41 and 0.45 ≤ T ∗ ≤ 25, dimensionless). In total, data for 27 615 ρ ∗ – T ∗ -state points including 750 configurations for each state point are used for the calculation of the pair and triplet correlation functions. For the pair configuration function, an approximation over the whole set of state points is provided, which reproduces the pair configuration with high accuracy (3 σ < 0.0075). The results for the triplet configuration functions are compared to the Kirkwood superposition approximation. With the exception of low ρ ∗ ( < 0.3), the application of the Kirkwood superposition approximation seems not to be a proper choice, showing errors > 20%. The configurations and pair and triplet correlation functions of all simulations are openly accessible and can be used as a reference for future theoretical developments of Lennard - Jones fluids and, especially, the liquid state.


I. INTRODUCTION
The Lennard-Jones potential 1 u LJ is one of the prototypes of an intermolecular interaction between particles, where r is the distance between particles.The potential u LJ is pair additive and it depends on two parameters ε and σ.The parameter σ is the solution for u LJ (r) = 0 and ε is the absolute depth of the trough of the u LJ -potential, as illustrated in Fig. 1.
Only two parameters are involved in Eq. ( 1) (compare with the Duhem theorem), so ε and σ can be used to make all variables and functions dimensionless, where T is the temperature, ρ is the number density, p is the pressure, N is the number of particles, kB is the Boltzmann constant, U is the internal energy, and F is the Helmholtz free energy.Dimensionless properties are designated by an asterisk * , and f * and u * are energies per particle.This feature of dimensionless parameters makes the Lennard-Jones potential especially attractive because simulations can be carried out with dimensionless variables ρ * and T * without knowing ε and σ in advance.Real ρ and T can be inserted after the simulations, using a proper choice of ε and σ, to make interpretations.Gottschalk 2 performed 27 579 Monte Carlo (MC) simulations using 1372 LJ-particles (N) in the stable and metastable single fluid phase field considering a large number of densities ρ * and temperatures T * (state points) providing an extensive set of f * -, u * -, and p * -values.These were used to derive an equation of state (EOS) for a Lennard-Jones fluid (LJ-EOS).Here, 36 additional MCsimulations are reported so the total amount increases to 27 615.In addition to f * , u * , and p * , these simulations did provide a large number of particle configurations, 750 per state point.We have derived pair and triplet correlation functions for all 27 615 simulations.The motivation for this simulation effort was to derive the ingredients for a general EOS, applicable for fluid mixtures at high pressures p and temperatures T in the geological and planetary settings.The general EOS will be based on perturbation theory 3-6 using a LJ-EOS 2 as a reference EOS.][9][10][11][12][13] In this article, the configurations and pair and triplet correlation functions are reported.The J-, K-, L-, M-, and N-integrals 11,14 required by perturbation theory will the subject of future communication.

ARTICLE
Particle configurations and pair correlation functions (i.e., radial distributions) have been the subject of numerous publications, starting from the theoretical work by Kirkwood 15 and the computational study by Rosenbluth and Rosenbluth. 16The present study provides an extensive and consistent database of calculated pair and triplet correlation functions over a large set of ρ * -T * state points.Such a database could be also of further use in the Kirkwood-Buff theory 17 to link the microscopic details to macroscopic properties.
In the range of 0 ≤ ρ * ≤ 1.41 and 0.45 ≤ T * ≤ 25, a continuous function for the pair correlation is presented, which is applicable to stable and metastable single fluid state points, and here, r * is not restricted to any upper value.
The direct calculation of the triplet correlation functions is scarce, [18][19][20][21][22][23] and the available range is limited to a few ρ * -T * state points.For applications of triplet correlation functions, therefore, the Kirkwood superposition approximation SA 15 is often used.The validity of the SA is discussed below in detail and is compared to the direct triplet calculations.
The here calculated configurations and pair and triplet correlation functions are also made available (see the Data Availability section).

II. MONTE CARLO SIMULATIONS
Gottschalk 2 performed 27 579 NVT MC-simulations in the stable and metastable single fluid phase field.This article adds 36 more simulations in this region for a total of 27 615.The simulations did provide u * -and p * -values.Values for f * were evaluated by integrating p * from simulations using the relationship However, integration is only possible when the integrand is continuously known starting from ρ * = 0.This is only the case for T * ≥ 1.3 or for vapor outside the two-phase field (ρ * < 0.3, T * < 1.3).For ρ * = 0, the Lennard-Jones fluid behaves as an ideal gas (p * = T * ρ * ) and the numerator is therefore 0 and so is the integrand.Simulation state points did range from 0.002 ≤ ρ * ≤ 1.41 and 0.45 ≤ T * ≤ 25 using a r * cutoff of 5.0 for ρ * ≤ 1.37 and 4.5 for ρ * > 1.37.At each state point, 1372 LJ particles (N = 1372) were equilibrated using 5 × 10 4 cycles and, subsequently, sampled for 0.75 × 10 6 cycles.One cycle did consist of 1372 trial moves, one for each particle.Every 1000 cycles, the thermodynamic properties u * and p * as well as the current configuration of particles were recorded.
Subsequently, the 750 stored configurations per state point were used to calculate the pair and triplet correlation functions using a bin width of 0.005 for pair and 0.04 in all directions for triplet correlations.For further details regarding correlation calculations, see below.
Stable and metastable fluid field boundaries, a table of state points, as well as additional details of the MC-simulations were also given in the work of Gottschalk. 2

III. THE PAIR CORRELATION FUNCTION A. Calculation of the pair correlation function
The pair correlation function g (2) (r * ) was calculated using the following recipe: 24,25 For the evaluation of J-integrals 11 required in perturbation theory, a continuous approximation of the pair correlation function would be convenient.Morsali et al. 26  According to de Boer and Michels, 27 the pair distribution function g (2) (r * 12 ) can be expressed in a virial-like expression (see also Kirkwood 15 ), In the limit of zero density, the pair correlation g ′ 0 function for a Lennard-Jones fluid is with For g ′ 1 and g ′ 2 , de Boer and Michels, 27 Rowlinson, 28 and Dyer et al. 29 provided the following expressions: where f ij is the Mayer function.In the special case of the Lennard-Jones potential, f ij is The variables ri are the vectors from molecule 1 to molecule i.Since for spherical molecules f ij is a function of r * ij only, it is convenient to use these scalar distances as the variables for integration. 28The Jacobian for the change from Cartesian to polar coordinates can be factorized into three Jacobians for molecules 2, 3, and 4 and can be taken independently.The details for this procedure can be found in the work of Rowlinson. 28he integrals in Eqs. ( 9) and ( 10) have been determined numerically in the range of 0.45 ≤ T * ≤ 25, resulting in 3 × 10 6 and 47 880 evaluations for g ′ 1 and g ′ 2 , respectively.The results have been fitted to the following equations: and using Mathematica. 30The respective coefficients for these equations are provided in the supplementary material.
In a brute-force attempt, an additional total empirical function g ′ 3 is fitted to the residuals left behind by Eq. ( 6) using g ′ 0 , g ′ 1 , and g ′ 2 .The fit involves 7.5 × 10 6 ρ * -T * -r data triples using the following function: The fitting was also performed with Mathematica.While over 4 × 10 7 triples are available from simulations, the memory requirement of Mathematica was extensive, so even with 128 GB RAM, only 7.5 × 10 6 triples could be handled.The reduced amount of data was selected randomly from the larger data pool.The fit procedure was repeated several times with different randomly chosen datasets, showing that the obtained solution was independent from the selection process.
The respective coefficients for Eq. ( 14) are provided in the supplementary material.

The complete (2)
12 approximation is then Figure 6 shows a histogram of the deviations of Eq. ( 15) for all simulation results.The deviations (3σ) are below 0.0075 and Eq. ( 15) is therefore valid for the whole range of 4 × 10 7 triples.

FIG. 6. Deviation of g
(2) 12 calculated by Eq. ( 15) to simulation results for all 40 502 166 simulation triples.Only a few triplet correlation functions for a few state points are published by Krumhansl and Wang, 18 Raveché et al., 19 Haymet et al., 20 McNeil et al., 21 Attard, 22 and Fushiki. 23Their limited results are consistent with the findings here.
One application for triplet correlation functions is the development and calculation of Kand L-integrals 11 required by the Pople 32 expansion of perturbation theory.

B. Calculation of the quadruplet correlation function
The calculation of the quadruplet correlation function g (4) 1234 needs sufficiently small values of Δr * , and as a consequence, a large amount of data storage is needed for its computation, and additionally, the computing time would be extensive.Therefore, the calculations were not performed here; however, they could be done using the available configurations.

FIG. 1 .
FIG. 1.The parameter σ is the solution for u LJ (r) = 0 and ε is the absolute depth of the trough of u LJ .