{'fraction_non_alphanumeric': 0.04209847771111221, 'fraction_numerical': 0.011900629741657162, 'mean_word_length': 5.022096148103912, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Predicting the structure of fluids with piecewise constant interactions: Comparing the accuracy of five efficient integral equation theories

We use molecular dynamics simulations to test integral equation theory predictions for the structure of fluids of spherical particles with eight different piecewise-constant pair interaction forms comprising a hard core and a combination of two shoulders and/or wells. Since model pair potentials like these are of interest for discretized or coarse-grained representations of effective interactions in complex fluids (e.g., for computationally intensive inverse optimization problems), we focus here on assessing how accurately their properties can be predicted by analytical or simple numerical closures including Percus-Yevick, hypernetted chain, reference hypernetted chain, first-order mean spherical approximation, and a modified first-order mean spherical approximation. To make quantitative comparisons between the predicted and simulated radial distribution functions, we introduce a cumulative structural error metric. For equilibrium fluid state points of these models, we find that the reference hypernetted chain closure is the most accurate of the tested approximations as characterized by this metric or related thermodynamic quantities.


I. INTRODUCTION
A common challenge in materials science is the "inverse design problem" [1,2], wherein one seeks to use theoretical models to discover the microscopic characteristics (e.g., the effective pair interactions) of a new system which, if fabricated or synthesized, would yield a targeted material property. Recent applications include designing materials that self-assemble into specific crystalline lattices [3][4][5][6], fluids that display optimized structural correlations and related transport properties [7][8][9], or solids that exhibit specific optical characteristics [10]. Inverse design problems are commonly addressed by stochastic optimization strategies like simulated annealing. Such approaches have the advantage of being general and easy to apply, and they can also be effective as long as material properties required for evaluating the objective function can be accurately and efficiently computed for large numbers of trial interactions during the optimization. This requirement typically means that "exact" yet computationally intensive methods for property determination (e.g., molecular simulations) are impractical for use within such calculations. Approximate theories with analytical or simple numerical solutions are attractive alternatives to molecular simulation in these contexts, provided that they can make sufficiently accurate predictions for a wide range of microscopic interaction types.
For bulk fluids, a key aim for property prediction is to discover the one-to-one link [11] between g(r), the radial distribution function (RDF) of a system at a given set of conditions, and ϕ(r), the interparticle pair potential. Knowledge of these functions of interparticle separation r allows for the direct calculation of the static structure factor, the energy, the pressure, and the isothermal compressibility [12]. Estimations of other properties can be * truskett@che.utexas.edu directly obtained from knowledge of the RDF as well. One example is the two-body excess entropy, which is often a good approximation of the total excess entropy [13] for simple liquids. Another is the information-theoretic estimate for the probability p n (Ω) of observing n particle centers in a molecular-scale subvolume Ω, a quantity which characterizes the fluid's density fluctuations [14]. Excess entropy, its two-body approximation, and p 0 have been shown to correlate with various dynamic properties of equilibrium fluids, e.g. diffusivity or viscosity [8,[15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31]. Mode-coupling theory also predicts that dynamic phenomena can be directly estimated from knowledge of the static structure factor [32].
With these considerations in mind, herein we use molecular simulations to test the accuracy of RDF predictions for five approximate integral-equation theory closures: Percus-Yevick, hypernetted chain and reference hypernetted chain [12], first-order mean spherical approximation (FMSA) [33], and a modified exponential version of FMSA [34]. Other more resource-intensive theories, like the Rogers-Young and hybrid mean-spherical approximations [35,36], self-consistent Ornstein-Zernike approaches [12], and thermodynamic perturbation theories [37][38][39] are not considered here. We apply the simpler five theories listed above to a diverse suite of eight pair potentials previously introduced by Santos et al [40], each composed of a hard core at r = σ plus two piecewise constant sections at larger r (i.e. wells or shoulders), that qualitatively mimic some of the features observed in the effective interactions of complex fluid systems. For each interaction, we investigate four thermodynamic state points with various combinations of low and high density and low and high temperature, and we compare the theoretical predictions for the RDF, the energy, and the two-body excess entropy to data from event-driven molecular dynamics simulations. To facilitate the RDF comparisons we introduce a "cumulative squared error" metric, which provides a quantitative characterization of the overall quality of each theoretical prediction. We also assess the accuracy of predictions for the potential energy and the two-body excess entropy.

A. Integral Equation Theory
Integral equation theories for uniform, isotropic fluids typically involve solving a system of two equations: the Ornstein-Zernike relation, which defines the direct correlation function c(r) in terms of the number density ρ and the total correlation function h(r) = g(r) − 1, and a closure, e.g., which introduces the link to the pair potential ϕ(r), where β = (k B T ) −1 , T is temperature, k B is Boltzmann's constant, and B(r) is the so-called bridge function. Two common approximations for B(r) are the Percus-Yevick (PY) closure, and the hypernetted chain (HNC) closure, Another is the so-called reference hypernetted chain approximation (RHNC), which assumes that the bridge function can be accurately approximated by that of a reference fluid, typically one of hard spheres at the same density: The hard-sphere fluid's bridge function B HS (r) has been calculated through careful molecular simulations, and multiple parameterizations for its density dependence exist [41][42][43]. For this work, we employ the analytical parameterization proposed by Malijevský and Labík [41] for the RHNC closure. With B(r) specified by these closures, we solve the coupled equations (1) and (2) using a rapidly-converging combination of Newton-Raphson and Picard root-finding methods developed by Labík et al. [44].
An alternative strategy is to replace the closure of Eq. 2 with separate expressions. For example, the mean spherical approximation (MSA) assumes the following relations hold, By further assuming first-order expansions in the characteristic dimensionless energy of the potential βε for both g(r) and c(r)-e.g., g FMSA (r) = g HS (r) + βεg 1 (r), where g HS (r) is the pair correlation function for a hard sphere system at the same density-Tang and Lu closed the equations analytically for several common pair interactions, including square wells [33]. We refer to this solution as the first-order mean spherical approximation (FMSA). In principle, FMSA can be applied to potentials with square shoulders as well. But for strong interactions, FMSA is known to incorrectly predict RDFs with negative values for some interparticle separations [34].
To resolve this, Hlushak et al. modified the FMSA to make it equally applicable to wells and shoulders by rearranging the terms in the series expansion, so that g EFMSA (r) = g HS (r) exp[−βεg 1 (r)] [34]. In this work, we refer to this analytical solution as the exponential firstorder mean spherical approximation (EFMSA).

B. Suite of Two-Step Potentials
Motivated by Santos et al. [40], we examine predictions for fluids from a set of pair interactions comprising a hard core and two piecewise-constant steps, where ε 1 and ε 2 are the energies of the first and second steps, respectively, and λ 1 and λ 2 are the outer edges of the first and second steps, respectively. Furthermore, as in Santos et al., we restrict the values of ε i to the set {−ε, −ε/2, 0, ε/2, ε}, where ε is a characteristic energy scale. Cases where ε 1 = ε 2 or ε 2 = 0 reduce to either single square wells or shoulders, or hard spheres, which have all been studied extensively elsewhere (see, e.g., refs. 1-41 in [45]) and are not considered here. We also exclude cases where max{|ε 1 |, |ε 2 |} = ε/2. Of the cases where ε 1 and ε 2 have opposite sign, we consider only combinations where ε 2 = −ε 1 = ±ε. We choose λ 1 = 1.5σ and λ 2 = 2σ in order to provide challenging perturbations to the bare hard sphere system that are still amenable to molecular simulation and theoretical treatment. After imposing these restrictions, the remaining eight pairwise interactions shown in Fig. 1, which we refer to as "Type A-H," form our test suite.
To explore how the accuracy of the various theories varies with density and temperature, we investigate each interaction at the four state points comprising combinations of packing fraction η = ρπσ 3 /6 = 0.15 or 0.45 and dimensionless temperature T * = k B T /ε = 0.67 or 2.0.

C. Molecular Simulations
We compare the theoretical predictions for the RDF, the energy, and the two-body excess entropy to the results of event-driven molecular dynamics simulations per-   [40], is topologically exhaustive (e.g., there are no other qualitative arrangements of two constant pairwise pieces that are not more appropriately labeled single wells or shoulders).
formed with the DynamO simulation engine [46]. Periodic boundary conditions were used, and the simulated systems were sized such that adequate RDF statistics could be collected for separations up to at least r = 10σ. In practice, this required N = 4000 particles when η = 0.15, and N = 8788 particles when η = 0.45. The "bins" for particle counts were 0.005σ wide. Temperatures were set and maintained using an Andersen thermostat [47].
Each simulation was initialized as an FCC lattice of the desired density at a high temperature, with randomly assigned particle velocities. After equilibrating for ten million events, the simulations were cooled to the desired temperature and re-equilibrated for a further ten million events. Then, the thermostat was removed, and the RDF was measured over the final five million events.

D. Quantifying Error in Predictions
To compare the various RDF theoretical predictions to simulations at a given state point, we define a metric we call the cumulative squared error, CSE(r): The integrand in the numerator characterize the squared deviation in the total correlation function between the prediction of a given theory h thy (r) and the result of the 'exact' simulation h sim (r); the power of two eliminates any possible cancellation of error, e.g. for cases where a theory both underpredicts and overpredicts the value of h(r) at different values of r. The denominator accumulates the total squared correlations in the simulated system, and thus normalizes the overall function to facilitate comparison between systems with different degrees of correlation (e.g., between low-density and high-density systems).
As r approaches infinity, all h(r) curves converge to zero and the CSE converges to a finite value, CSE ∞ : which is a measure of the summed squared correlations as a fraction of the total squared correlations in the system; thus, a larger value of CSE ∞ indicates that a theoretical prediction deviates more significantly from the "exact" simulation results. By construction, CSE ∞ has a defined minimum of 0 and, while it does not have a rigorous maximum, its value is typically less than 1 except in cases where the theoretical predictions are qualitatively very poor.
We also calculate the potential energy per particle U/ε, and the two-body contribution to excess entropy s (2) /k B , from simulations and theoretical predictions. Both quantities can also be directly computed from g(r) and thus, the normalized absolute deviation of the predicted versus simulated values can be used as an indication of the success of theoretical predictions. However, note that different RDFs can, in principle, give rise to the same value of U/ε or s (2) /k B . Moreover, U/ε only depends on correlations within the range of the pair interaction. As a result, we argue here that since the RDF is weighted differently for each thermodynamic quantity, the CSE metric we introduce-which tests the overall similarity between predicted and simulated RDFs-represents a more sensitive measure for the overall predictive quality of particular theory.  [12,41]; the first-order mean spherical approximation solution (FMSA) [33]; and the simple exponential first-order mean spherical approximation (EFMSA) [34], for the "type A" pair interaction. Shaded regions adjacent to each g(r) indicate the difference between the theory and simulation results.

III. RESULTS AND DISCUSSION
Structural predictions for the Type A pair interaction are compared to simulation results in Fig. 2, along with the corresponding cumulative squared errors as calculated via Eq. (8). For this interaction, the analytic solutions (FMSA and EFMSA) perform better at higher rather than at lower equilibrium fluid densities. As density increases, the effect of the excluded volume captured by the well-modeled hard-sphere RDF, g HS (r), overwhelm the energetic perturbations from the repulsive steps and dominate the resulting structure. Of the tested integral-equation theories with simple numerical closures, the PY closure tends to perform least well near contact, and for interaction Type A, the RHNC offers the best predictions at all four state points investigated. Analogous figures for each of the other interactions are presented for the interested reader in Appendix A.
It is tempting to conclude from a visual comparison of theoretical and simulated radial distribution functions that all of the theories perform similarly well, especially at the higher temperature (Figs. 2d and 2h). However, the resulting CSEs differ by nearly two orders of magnitude from most to least accurate (Figs. 2c and 2g), which underscores the utility and sensitivity of the CSE metric. As discussed below, these differences in the CSE become important when computing other quantities that depend on the RDF, especially when one considers that each thermodynamic quantity weights the RDF in a different way.
The total cumulative squared errors CSE ∞ for all interactions, state points, and theories are listed in Table  I. Six of the total thirty-two combinations of interaction type and state point considered did not produce singlephase, uniform fluids when simulated. Of the remaining twenty-six systems, the RHNC offered the most accurate structural predictions for all but four; however, at three of these four points, the CSE ∞ of the RHNC is still within ca. 65% of the most accurate theory (HNC). All four points are at low temperature (T * = 0.67) and high packing fraction (η = 0.45), and each of the pair interactions include attractions (types D, F, G, and H).
We also compare CSE ∞ against the absolute normalized errors for predictions of two example thermodynamic quantities, two-body excess entropy s (2) /k B and potential energy U/ε, in Fig. 3. Fig. 3a shows that CSE ∞ is generally a good predictor of s (2) /k B accuracy, although there are a handful of instances where the fractional error in the excess entropy is low while CSE ∞ is higher.The correlation between CSE ∞ and the potential energy is weaker, but still present; this is likely due to opportunities for fortuitous cancellation of error when pair interactions contain both positive and negative contributions (e.g., types D and H), when portions of the interactions are zero (types C and G), or when significant contribu-TABLE I. Total cumulative squared errors (CSE∞) for all theoretical approaches, thermodynamic state points, and interactions considered. "R," "H," and "P" are the RHNC [12,41], HNC [12], and PY [12] closures to the Ornstein-Zernike relation, respectively. "F" is the FMSA [33], and "E" is the EFMSA [34]. Italics indicate the lowest value of CSE∞ (and hence the theory with the most accurate structural prediction) at each combination of state point and interaction type.  tions to CSE ∞ occur beyond the range of the pair interaction. Overall, however, it is clear that the accuracies of both example thermodynamic quantity predictions correlate well with the cumulative squared error. For the interested reader, the values of s thy /s (2) sim − 1 and |(U thy /U sim ) − 1| are tabulated in Appendix B. If other thermodynamic quantities that depend on the RDF in a different way (e.g., the pressure or the isothermal compressibility) are also of interest, then the necessity to have an independent structural metric like CSE ∞ to assess the quality of the structural predictions is even more critical.

IV. CONCLUSION
In order to quantify the overall accuracy of theoretical predictions for fluid structure, we have introduced the total cumulative squared error (CSE ∞ ) metric, which accumulates squared discrepancies between a theoretical prediction and a reference "exact" result at all separation distances along the total correlation function and avoids any possible cancellation of error. We find that this CSE ∞ metric is very sensitive and tends to forecast the overall accuracy of structure-dependent thermodynamic calculations. As a result, it is an excellent tool for comparing accuracy between multiple theories, particularly when differences are difficult to discern by visual inspection.
We have used this metric to test the performance of five integral equation theory-based approaches for predicting equilibrium fluid structure in systems with pair interactions comprising a hard core plus two piecewise constant interactions, and we find that the reference hypernetted chain (RHNC) integral equation closure offers accurate and efficient predictions across a broad range of interactions and thermodynamic state points. This kind of analysis, i.e., considering the accuracy of various efficient theoretical methods for predicting the structure consistent with a broad range of possible interactions, will be particularly important for inverse design problems where the goal is to rather accurately predict which interaction is consistent with a targeted structure (or structurallyrelated property).

ACKNOWLEDGMENTS
The authors thank Anatol Malijevský for sharing his rapidly-converging integral equation theory code. T.M.T. acknowledges support of the Welch Foundation (F-1696) and the National Science Foundation (CBET-1403768). We also acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing HPC resources that have contributed to the research findings reported within this paper.

Appendix A: Extended Type B-H Structure Plots
The predicted radial distribution functions g(r) compared against simulation results, and the resulting cumulative squared errors CSE(r), are shown for interaction types B through H in Figs. 4-10, respectively.  [12,41]; the first-order mean spherical approximation solution (FMSA) [33]; and the simple exponential first-order mean spherical approximation (EFMSA) [34], for the "type A" pair interaction. Shaded regions adjacent to each g(r) indicate the difference between the theory and simulation results.      Radial distribution functions and cumulative squared errors for the "type H" interaction. Series and labeling are as in Fig. 4.
thy /s (2) sim ) − 1|, for all approaches, state points, and interactions considered. Labels are as in Table I. Italics indicate the value closest to zero (e.g., a perfect prediction) at each combination of state point and interaction type.