Topological generalization of the rigidnonrigid transition in softsphere and hardsphere fluids
Abstract
A fluid particle changes its dynamics from diffusive to oscillatory as the system density increases up to the melting density. Hence, the notion of the Frenkel line was introduced to demarcate the fluid region into rigid and nonrigid liquid subregions based on the collective particle dynamics. In this work, we apply a topological framework to locate the Frenkel lines of the softsphere and the hardsphere models relying on the system configurations. The topological characteristics of the ideal gas and the maximally random jammed state are first analyzed, then the classification scheme designed in our earlier work is applied. The classification result shows that the fraction of solidlike atoms increases from zero to one in the rigid liquid region. The dependence of the solidlike fraction on the bulk density is understood based on the theory of fluid polyamorphism. The percolation behavior of solidlike clusters is described based on the fraction of solidlike molecules in an integrated manner. The crossover densities are obtained by examining the percolation of solidlike clusters. The resultant crossover densities of softsphere fluids converge to that of hardsphere fluid. Hence, the topological method successfully highlights the generality of the Frenkel line.
I Introduction
Hardsphere and softsphere models have been widely used to simulate the dynamics and structure of the fluid phase. The hardsphere model is one of the most extensively studied model in statistical physics. The softsphere model, in which the penetrability of a sphere depends on the slope of the repulsive wall, has also been used as a simple model. In this model, the pair potential is given by
(1) 
where is the interatomic pair potential between two spheres of diameter separated by a distance , is the energy parameter, and is the repulsive exponent. In these repulsive systems, no firstorder gasliquid transition occurs because the attractive interaction is absent. Hence, it is believed that the hardsphere and the softsphere systems would follow the same dynamics scheme in the whole fluid region Bacher et al. (2014).
Brazhkin et al. questioned this continuous picture of dynamics Brazhkin et al. (2012). Based on the phonon theory Bolmatov et al. (2012, 2015a) and on experimental validations Bolmatov et al. (2015b), they proposed the notion of the Frenkel line. They argued that the particle dynamics of supercritical fluid, a state of matter beyond the gasliquid critical point, changes from diffusive (gaslike) to oscillatory (solidlike) across the Frenkel line. They proposed to use the heat capacity criterion and the velocity autocorrelation function to locate the Frenkel line of supercritical fluids Brazhkin et al. (2013). They further applied these thermodynamic and dynamic criteria to more general classes of fluid models such as the slightly softsphere model and the hardsphere model.
However, a series of recent works on the phonon theory revealed that these thermodynamic and dynamic criteria could not be directly extended to locate the Frenkel line of all fluid models. Bryk et al. noted that the Frenkel lines of the softsphere fluids located from these thermodynamic criteria did not converge to the density where the anomalous behavior of the hardsphere models was observed Bryk et al. (2017a). The nonconvergence of the Frenkel lines of the softsphere fluids to that of the hardsphere fluid was critical since the transport properties of the softsphere fluids including selfdiffusion, shear and bulk viscosity, and thermal conductivity coefficients converge to those of the hardsphere fluid Borgelt and Hoheisel (1989); Heyes and Brańka (2005).
The nonconvergence of the Frenkel lines mainly originates from the quasicrystalline approximation (QCA), which is the basis of the conventional thermodynamic and dynamic criteria. Khrapak et al. examined the validity of the QCA for the softsphere potentials and reported that the QCA fails when the repulsive exponent of the softsphere potential is higher than 20 Khrapak et al. (2017). Hence, Brazhkin et al. recently proposed to use the anomalous transport properties as criteria to locate the rigidnonrigid transition density of the hardsphere model Brazhkin et al. (2018). Yoon et al. proposed a dynamic criterion based on the notion of the solidicity from the twophase thermodynamic (2PT) model Yoon et al. (2018a). They observed that the solidicity, which is defined as the ratio of the diffusivity of a system to the diffusivity of the hardsphere system at the zero pressure limit, shows an inflection behavior near the dynamic crossover density. Since the solidicity does not depend on the QCA, they demonstrated that the Frenkel lines of the soft spheres located from the new criterion converge to that of the hardsphere model Brazhkin et al. (2018).
In contrast to the dynamicsbased approaches, there have been only a few works that attempted to locate the Frenkel line based not on the collective particle dynamics but on the system configuration. Bolmatov et al. proposed that the third maximum of the pair correlation function be related to the structural crossover across the Frenkel line Bolmatov et al. (2014). Ghosh et al. applied this observation to characterize the Frenkel line of the confined fluids Ghosh and Krishnamurthy (2018). However, Bryk et al. criticized the use of the pair correlation function as a geometrical evidence of the Frenkel line since the intensity of the third maximum was so small that it cannot be distinguished from either the thermal noise or numerical errors Bryk et al. (2017b). Fomin et al. discovered that the packing fraction of effective hard spheres reaches the percolation threshold near the Frenkel line Fomin et al. (2014). Ryltsev et al. proposed that the concentration of tetrahedral clusters is higher than the percolation threshold above the Frenkel line Ryltsev and Chtchelkatchev (2013). Unfortunately, they could explain how the solidlike structure evolves in the rigid liquid region only approximately.
In our last work, thus, we suggested a topological classification procedure to analyze the dynamic crossover of supercritical argon across the Frenkel line Yoon et al. (2018b). In this procedure, geometric details of the Voronoi cells are not used to classify a particle as either gaslike or solidlike. Instead, the topological information, or the local connectivity of an atom to its neighbors, was used to classify whether an atom resembles the ideal gas or the maximally random jammed (MRJ) state. Based on this procedure, we discovered that the fraction of the solidlike particles steeply increases over the rigid liquid region, which is enclosed by the Frenkel line and the freezing line. This result provided the physical meaning of the Frenkel line as a percolation transition line from the ideal gas to the MRJ state.
In this work, we apply the designed method to analyze the rigidnonrigid crossover of the softsphere and hardsphere fluids. We first examine the topological characteristics of the ideal gas and the MRJ state, respectively. The topological classification method is then applied to locate the Frenkel lines of softsphere and hardsphere models. The dependence of the solidlike fraction on the bulk density is explained from the scope of the fluid polyamorphism theory and the isomorph theory. We further represent the generality of the dynamic crossover in terms of the solidlike fraction. The rigidnonrigid crossover densities of softsphere fluids converge to that of the hardsphere model as the repulsive exponent increases. These results substantiate that the topological framework can be used to locate the Frenkel line of general types of potentials.
Ii Methods
ii.1 Molecular Dynamics (MD) simulations
We performed timedriven NVT simulations Plimpton (1995) of fluids modeled with the repulsive potentials [Eqn. (2)]
(2) 
where is given as:
(3) 
Each system contained 2,000 particles. The repulsive exponent was chosen from . The size parameter () was 3.405 Å, the same as that of argon modeled with the LennardJones potential. Energy parameters () were determined so that the coefficients become equal to where is equal to , the energy parameter of argon. The potentials were shifted and truncated at the cutoff radius where the potential has its minimum. The simulation temperatures were selected as . The timestep was varied from to depending on the repulsive exponent () and the simulation conditions. The systems were equilibrated for 100,000 steps. After the equilibration, they were run for additional 100,000 steps to obtain the configurations and the system pressures. We additionally performed MD simulations of 16,000 particles modeled with the repulsive (WeeksChandlerAndersen Weeks et al. (1971), WCA) potential at to examine the influence of the finite size effect.
For the hardsphere systems, the eventdriven Molecular Dynamics (EDMD) simulations Bannerman et al. (2011) were conducted. For each simulation condition, two hundred configurations of 2,048 hard spheres (, ) were collected with a fixed dimensionless time interval of 1.00. The packing fractions () were from 0.06 to 0.54 where is the number density of particles ().
ii.2 Topological framework for the local structure analysis
The dynamics of a particle in a dense system is dominated by the relative positions of its surrounding neighbors. Hence, we conjectured that the topological framework for local structure analysis proposed by Lazar et al. Lazar et al. (2015) would build a robust connection between the dynamics and the geometry. In this framework, the arrangement of neighbors surrounding a central particle is described based on the Voronoi tessellation, the partitioning of space into regions which consist of all points closer to a central particle than to any other. A Voronoi cell provides the connectivity information as the geometric feature of a face and its adjacent faces. For instance, the number of faces of a Voronoi cell is equal to the number of nearest neighbors. An sided polygonal face of a Voronoi cell implies that the two atoms which share that face have common neighbors. The connectivity information of faces wrapping around a Voronoi cell can be encoded as a series of the integers called the Weinberg vector Lazar et al. (2012). Hence, the Weinberg vector can be viewed as a name of the type of a Voronoi cell.
This topological description of the local configuration is in line with the isomorph theory. The isomorph theory Schrøder and Dyre (2014) states that the following simple relation is satisfied for any two configurations and if they are isomorphic to each other:
(4) 
Here, is the location vector of the particle , is the number of atoms in the system, and is the Boltzmann statistical weight of the given configuration. In other words, two configurations are regarded to be isomorphic if their configurations in the reduced unit are the same. Similarly, the topological description does not depend on the distances between the central particle and its neighbors. Two local configurations are considered the same if their Voronoi cells are topologically identical regardless of their volumes.
ii.3 Topological characterization of the dynamic limits
Since the Frenkel line was originally defined as the thermodynamic states where the particle dynamics change from diffusive to oscillatory, it was required to select the configurations that are opposite to each other from the viewpoint of the dynamics. A manybody system in which particles only translate without the interference of their neighbors would be the ideal gas. In contrast, systems in which most particles are randomly distributed and only vibrate in their place would be the MRJ state Klatt and Torquato (2014). Therefore, we constructed 50 samples of 500,000 particles each of both systems. The configurations of the ideal gas were generated by distributing points randomly, whereas those of the MRJ state were produced by compressing the hardsphere system using the LubachevskyStillinger algorithm Lubachevsky and Stillinger (1990). We then analyzed these samples using the opensource VoroTop software Lazar (2017), which computed the distribution of topological features of the model systems.
ii.4 Topological classification strategy
Fig. 1 schematically describes the topological classification strategy used in this work. The first step of the classification scheme is based on the list of the Weinberg vectors discovered in the dynamic limits. First, we compute the list of the Weinberg vectors of soft (hard) spheres obtained from the MD simulations using VoroTop. We then compare the likelihood of finding the Weinberg vector of a single atom in the list of the Weinberg vectors of the ideal gas () to that of the MRJ state (). If , the atom is classified as solidlike ( where is a state number). Otherwise, it is classified as gaslike (). After this initial classification, a weighted meanfield strategy Yoon et al. (2018c) is used to reclassify an atom based on its state number and those of its neighbors. The averaged state number is defined as
(5) 
where is the number of the nearest neighbor atoms of the particle. This procedure makes it possible to remove the influence of small local fluctuation on the classification result. Note that solidlike local structures might appear in ideal gas by pure chance, since the ideal gas essentially includes every possible configurations: if the topological framework alone is used, a nonzero fraction of particles in ideal gas are always classified as solidlike. By applying the weighted mean field procedure, the probability to find a solidlike molecule in the ideal gas converges to zero. As a last step, the averaged state number is rounded to decide the state of molecule .
ii.5 Percolation analysis
Percolation theory Stauffer and Aharony (2014) has frequently been used to analyze the structural characteristics of the connected molecules (clusters) discovered in the fluid phase. According to percolation theory, an infinite (spanning) cluster appears when the particle concentration is higher than a particular concentration of the particles called the percolation threshold (). Here, a cluster is defined as a connected assembly of atoms that are classified as the same state, either solidlike or gaslike. From the viewpoint of the Voronoi tessellation, two particles can be regarded to be connected if they share a face with each other. Hence, two atoms belong to the same cluster if they are connected through the atoms whose classification results are identical to theirs. A clustering algorithm proposed by Stoll Stoll (1998) is used to detect the cluster structure. In the first step of the algorithm, the list of clusters is obtained without consideration of the periodic boundary conditions. In the second step, the algorithm determines whether a cluster is infinite or not by examining the connectedness of two Voronoi cells which are located at opposite sides of the simulation box and assigned to the same cluster. If they are Voronoi neighbors across the boundary, the cluster is regarded to be a spanning (infinite) cluster. After this test, the independent clusters which are connected through the periodic boundaries are assigned to a single cluster if two atoms in these clusters are connected to each other.
Iii Results and Discussion
iii.1 Topological characterization of the ideal gas and the MRJ state
Ideal gas  










MRJ  








As the topology of each Voronoi cell encodes the manner in which neighbors of a particle are arranged Lazar et al. (2015), the distribution of Voronoi topologies provides a meaningful statistical description of local ordering in a system. In crystals, the set of possible Voronoi types is finite, restricted by the symmetry of the crystal and by the manner in which unstable vertices of an ideal structure can resolve Lazar et al. (2015). In contrast, the set of possible Voronoi types in an ideal gas is infinite, as almost all arrangements of particles are possible with some finite probability. Because the number of topological types is infinite, it is clear that not all arrangements are equally likely. Indeed, in prior work we documented the distribution of Voronoi topologies in the ideal gas and found that certain ones appear more frequently than others Lazar et al. (2013). In this section we report data from the topological analysis of the ideal gas and MRJ systems.
The average number of faces in the ideal gas and MRJ systems were 15.54 () and 14.28 (), respectively. The higher concentration of the distribution about the mean in the MRJ system suggests more order in that system than in the ideal gas; this similarly expresses itself in terms of the Voronoi indices and Voronoi topologies. Table 1 shows the 10 most common sets of Voronoi indices in the two system. Each set of indices, also called a vector Lazar et al. (2012), records the number of faces with a given number of edges, beginning with 3. The ten most frequent vectors in the MRJ system account for almost a third of all constituent atoms; in contrast, the ten most common vectors in the ideal gas account for less than three percent. The relative dispersion of vectors in the ideal gas case can be seen as reflecting the relative lack of order in that system, as compared with the MRJ one.
Fig. 2 illustrates the 8 most common Voronoi topologies in both the ideal gas and the MRJ systems, along with their estimated frequencies. The distribution of Voronoi topologies in the MRJ system is more concentrated among a smaller number of types, with each of the eight most frequent types accounting for at least 1.4% of all atoms. In contrast, in the ideal gas, no single Voronoi topology accounts for more than 0.27%. This can again be understood as reflecting the relative disorder of the ideal gas as compared with the MRJ system.
Of additional potential interest in studying the relationship between the local neighborhood of a particle and its dynamics is the symmetry group associated with its Voronoi cell Weinberg (1966). The average symmetry group order is 1.161 in the ideal gas, and 3.365 in the MRJ system. Table 2 lists the number of atoms in each data set with nontrivial symmetry groups. The MRJ system has many more Voronoi cells with highorder symmetries than does the ideal gas. In particular, the number of atoms with icosahedral symmetry () is three orders of magnitude larger in the MRJ than in the ideal gas.
The distribution of symmetry orders in the two systems might be correlated with their dynamics. In particular, the more symmetric a particle’s Voronoi cell, the more confined that particle will be due to the influence of its symmetricallyarranged neighbors. For reference, in bodycentered cubic crystals, the order of the symmetry group of each Voronoi cell is 48; other crystals have similarly highorder symmetries. Atoms with highlysymmetric arrangements of neighbors, and thus highlysymmetric Voronoi cells, might be thought of as trapped in a potential “cage”.
Topological features of the Voronoi cells in the ideal gas and MRJ systems thus provide a quantitative description of local structure, and provide reference states against which to compare more general fluids. In a future paper we will provide a statistical description of a more general class of fluids.
iii.2 Frenkel lines of softsphere and hardsphere fluids
Fig. 3 shows the fraction of solidlike molecules () of the softsphere fluids and the hardsphere fluid. For all repulsive exponents and simulation temperatures, ’s show a sigmoidal dependence on the density [Eqn. (6)].
(6) 
They start to steeply increase near the dynamic crossover densities obtained from the 2PT model and converge to unity near the freezing densities.
The sigmoidal dependence of on the density can be understood from the viewpoint of fluid polyamorphism Anisimov et al. (2018); Yoon et al. (2018c); Ha et al. (2018). According to the theory of fluid polyamorphism, the interconversion of gaslike and the solidlike states can be expressed as the chemical reaction (). The equilibrium constant of the interconversion is represented as
(7) 
where is the Gibbs energy difference between the two states. Assuming that the Gibbs free energy of the interconversion () is proportional to , Eqn. (7) can be transformed into Eqn. (6).
Figure 3 demonstrates that the structural evolution of softsphere fluids approaches that of the hardsphere fluid as the repulsive exponent is increased, where the hardsphere fluid is understood as the limit. When the repulsive exponent is small (Fig. 3a), curves largely depend on the temperatures. As the repulsive exponent increases (Fig. 3b and c), different isothermal curves of more closely resemble that of the hardsphere fluid.
Figure 3 also shows that the crossover density is smaller in systems with larger . As Brazhkin et al. stated Brazhkin et al. (2018), the rigidnonrigid crossover densities would be determined by the cage effect of surrounding molecules. The extent of the cage effect is determined by the competition between the relative kinetic energy of the central particle and the softness of the repulsive wall of its neighbors. When the repulsive exponent is high or the simulation temperature is low, particles are easily arrested by their nearest neighbors. Hence, the crossover density of the hardsphere fluid should be lower than any softsphere fluids. When is small and the local cage is softer, more neighbors are required to trap the central particle at the high temperature.
Interestingly, while the 2PT model and the theory of collective phonon dynamics based on the Frenkel frequency understand this cage effect based on the dynamics of a particle, the topological classification method quantifies it based on the fraction of solidlike molecules. This scenario implies that the attractive interaction would not play an important role in determining the location of the Frenkel line. Fig. 4 demonstrates this idea; the crossover densities of the LJ fluid from our earlier work, where the interatomic potential includes attraction, are almost consistent with those of the repulsive 126 fluid. When the system temperature is low, they slightly disagree with each other, but the discrepancy between these systems decreases as the temperature increases. Hence, the rigidnonrigid crossover can be well explained by the hardsphere paradigm, which states that the excluded volume effects dominate the liquid behavior Weeks et al. (1971); Dyre (2016).
Next, we analyze the solidlike structures of the softsphere and the hardsphere fluid systems from the viewpoint of the percolation theory. When the probability of finding an infinite cluster in a configuration () is expressed as a function of the , curves at different conditions collapse to a single line obtained from the LJ simulations Yoon et al. (2018b) regardless of the repulsive exponents and temperatures (Fig. 5a). Moreover, both curves of LJ and WCA fluid show the same dependence on the system size (see Fig. A1b in the Appendix).
This result demonstrates that the rigidnonrigid transition across the Frenkel line is quasiuniversal; it does not depend on the types of interactions. In line with the isomorph theory, it can be restated that two configurations with the same are isomorphic, showing the same percolation behavior. Provided that the percolation of solidlike clusters is related to the thermodynamic and transport properties of general fluid systems, the generality of the Frenkel line would have a deep relation with the excess entropy scaling Rosenfeld (1977); Jakse and Pasturel (2016), which will be dealt with in our future studies.
Since the percolation behaviors of the softsphere and hardsphere fluids are equal to that of the LJ fluid, the percolation threshold obtained from our previous work () can be used to locate the dynamic crossover densities. Fig. 5b shows the crossover densities of softsphere and hardsphere fluids. As shown in the 2PT method Yoon et al. (2018a), the dynamic crossover densities of softsphere fluids converge to that of hardsphere fluid, which corresponds to . This result again demonstrates the advantage of the topological classification method. The thermodynamic and dynamic criteria by Brazhkin et al. are based on QCA, and the crossover densities from these criteria do not converge to that of the hardsphere systems where QCA breaks down Bryk et al. (2017a). The solidicity criterion from the 2PT model is free of QCA, yet it relies on CarnahanStarling equation of state Carnahan and Starling (1969), which is an approximate model. On the other hand, the topological criterion from this work does not rely on hypotheses that are constrained to the repulsive exponents. Moreover, it does not require a vast amount of the postprocessing procedure and data to obtain the thermodynamic properties of a system; it only requires the location of the particles and the simulation box length. Overall, these results show that the topological framework successfully generalizes the notion of the rigidnonrigid transition in the fluid models.
Iv Conclusions
The topological framework sheds light on the generalization of the notion of the rigidnonrigid crossover. It not only offers physical insight into the relationship between the dynamics and the geometry of particles but also overcomes a limit of the conventional methods based on the thermodynamics and dynamics of the system. The dynamic limits of fluid particle systems can be clearly characterized based on their topological characteristics. The topological framework deduced from these dynamic limits provides a classification scheme that can locate the rigidnonrigid transition of softsphere and hardsphere systems. The fraction of solidlike molecules () from the classification method can be used as an order parameter to describe the rigidnonrigid transition in an integrated manner. The dynamic crossover densities of the softsphere particles converge to that of the hardsphere particles, which was also observed in the 2PT model. Hence, it would be advantageous to expand our understanding of the fluid physics as well as to calculate the thermodynamic properties of the fluid systems.
Acknowledgements
E. A. L. gratefully acknowledges the generous support of the US National Science Foundation through Award DMR1507013. M.Y.H. acknowledges the support of the National Research Foundation of Korea Grants funded by Korean Government (NRF2017H1A2A1044355).
Appendix: Finitesize effect
Fig. A1a compares the fraction of solidlike molecules modeled with the repulsive potential at . The fraction of solidlike molecules does not change significantly when the number of molecules in the system increases. Hence, the topological classification results do not depend on the size of the system. In contrast, the percolation behavior of the system depends on the system size as shown in our earlier work Yoon et al. (2018b). Compared to the system with 2,000 molecules, of the system shows more abrupt increase when the system density increases. As shown in Fig. A1b, the dependence of on of the repulsive 126 (WCA) fluid is the same as that of the LennardJones fluid. Again, this result substantiates that the percolation of solidlike structures in the rigid liquid region is universal.
References
 Bacher et al. (2014) A. K. Bacher, T. B. Schrøder, and J. C. Dyre, Nat. Commun. 5, 5424 (2014).
 Brazhkin et al. (2012) V. Brazhkin, Y. D. Fomin, A. Lyapin, V. Ryzhov, and K. Trachenko, Phys. Rev. E 85, 031203 (2012).
 Bolmatov et al. (2012) D. Bolmatov, V. Brazhkin, and K. Trachenko, Sci. Rep. 2, 421 (2012).
 Bolmatov et al. (2015a) D. Bolmatov, D. Zavâyalov, M. Zhernenkov, E. T. Musaev, and Y. Q. Cai, Ann. Phys. 363, 221 (2015a).
 Bolmatov et al. (2015b) D. Bolmatov, M. Zhernenkov, D. Zavâyalov, S. N. Tkachev, A. Cunsolo, and Y. Q. Cai, Sci. Rep. 5, 15850 (2015b).
 Brazhkin et al. (2013) V. Brazhkin, Y. D. Fomin, A. Lyapin, V. Ryzhov, E. Tsiok, and K. Trachenko, Phys. Rev. Lett. 111, 145901 (2013).
 Bryk et al. (2017a) T. Bryk, A. Huerta, V. Hordiichuk, and A. Trokhymchuk, J. Chem. Phys. 147, 064509 (2017a).
 Borgelt and Hoheisel (1989) P. Borgelt and C. Hoheisel, J. Chem. Phys. 91, 7872 (1989).
 Heyes and Brańka (2005) D. Heyes and A. Brańka, Phys. Chem. Chem. Phys. 7, 1220 (2005).
 Khrapak et al. (2017) S. Khrapak, B. Klumov, and L. Couëdel, Sci. Rep. 7, 7985 (2017).
 Brazhkin et al. (2018) V. Brazhkin, Y. D. Fomin, V. Ryzhov, E. Tsiok, and K. Trachenko, Physica A (2018).
 Yoon et al. (2018a) T. J. Yoon, M. Y. Ha, W. B. Lee, and Y.W. Lee, arXiv preprint arXiv:1806.07608 (2018a).
 Bolmatov et al. (2014) D. Bolmatov, D. Zavâyalov, M. Gao, and M. Zhernenkov, J. Phys. Chem. Lett. 5, 2785 (2014).
 Ghosh and Krishnamurthy (2018) K. Ghosh and C. Krishnamurthy, Phys. Rev. E 97, 012131 (2018).
 Bryk et al. (2017b) T. Bryk, F. A. Gorelli, I. Mryglod, G. Ruocco, M. Santoro, and T. Scopigno, J. Phys. Chem. Lett. 8, 4995 (2017b).
 Fomin et al. (2014) Y. D. Fomin, V. Ryzhov, E. Tsiok, V. Brazhkin, and K. Trachenko, Sci. Rep. 4, 7194 (2014).
 Ryltsev and Chtchelkatchev (2013) R. Ryltsev and N. Chtchelkatchev, Phys. Rev. E 88, 052101 (2013).
 Yoon et al. (2018b) T. J. Yoon, M. Y. Ha, E. A. Lazar, W. B. Lee, and Y.W. Lee, arXiv preprint arXiv:1807.02761 (2018b).
 Plimpton (1995) S. Plimpton, J. Comput. Phys. 117, 1 (1995).
 Weeks et al. (1971) J. D. Weeks, D. Chandler, and H. C. Andersen, J. Chem. Phys. 54, 5237 (1971).
 Bannerman et al. (2011) M. N. Bannerman, R. Sargant, and L. Lue, J. Comput. Chem. 32, 3329 (2011).
 Lazar et al. (2015) E. A. Lazar, J. Han, and D. J. Srolovitz, Proc. Natl. Acad. Sci. U.S.A. 112, E5769 (2015).
 Lazar et al. (2012) E. A. Lazar, J. K. Mason, R. D. MacPherson, and D. J. Srolovitz, Phys. Rev. Lett. 109, 095505 (2012).
 Schrøder and Dyre (2014) T. B. Schrøder and J. C. Dyre, J. Chem. Phys. 141, 204502 (2014).
 Klatt and Torquato (2014) M. A. Klatt and S. Torquato, Phys. Rev. E 90, 052120 (2014).
 Lubachevsky and Stillinger (1990) B. D. Lubachevsky and F. H. Stillinger, J. Stat. Phys. 60, 561 (1990).
 Lazar (2017) E. A. Lazar, Modelling Simul. Mater. Sci. Eng. 26, 015011 (2017).
 Yoon et al. (2018c) T. J. Yoon, M. Y. Ha, W. B. Lee, and Y.W. Lee, J. Chem. Phys. 149, 014502 (2018c).
 Stauffer and Aharony (2014) D. Stauffer and A. Aharony, Introduction to Percolation Theory: Revised Second Edition (CRC press, 2014).
 Stoll (1998) E. Stoll, Comput. Phys. Commun. 109, 1 (1998).
 Lazar et al. (2013) E. A. Lazar, J. K. Mason, R. D. MacPherson, and D. J. Srolovitz, Phys. Rev. E 88, 063309 (2013).
 Weinberg (1966) L. Weinberg, SIAM J. Appl. Math. 14, 729 (1966).
 Anisimov et al. (2018) M. A. Anisimov, M. Duška, F. Caupin, L. E. Amrhein, A. Rosenbaum, and R. J. Sadus, Phys. Rev. X 8, 011004 (2018).
 Ha et al. (2018) M. Y. Ha, T. J. Yoon, T. Tlusty, Y. Jho, and W. B. Lee, J. Phys. Chem. Lett. 9, 1734 (2018).
 Dyre (2016) J. C. Dyre, J. Phys. Condens. Matter 28, 323001 (2016).
 Rosenfeld (1977) Y. Rosenfeld, Phys. Rev. A 15, 2545 (1977).
 Jakse and Pasturel (2016) N. Jakse and A. Pasturel, Sci. Rep. 6, 20689 (2016).
 Carnahan and Starling (1969) N. F. Carnahan and K. E. Starling, J. Chem. Phys. 51, 635 (1969).