Mapping diffusivity of narrow channels into one-dimension
MMapping diffusivity of narrow channels into one-dimension
Mahdi Zarif and Richard K. Bowles Department of Physical and Computational Chemistry,Shahid Beheshti University, Tehran 19839-9411, Iran ∗ Department of Chemistry, University of Saskatchewan, Saskatoon, S7N 5C9, Canada † (Dated: revised version – January 14, 2020)The diffusion of particles trapped in long narrow channels occurs predominantly in one dimension.Here, molecular dynamics simulation is used to study the inertial dynamics of two-dimensional harddisks, confined to long, narrow, structureless channels with hard walls in the no-passing regime. Weshow that the diffusion coefficient obtained from the mean squared displacement can be mappedonto the exact results for the diffusion of the strictly one-dimensional hard rod system throughan effective occupied volume fraction obtained from either the equation of state or a geometricprojection of the particle interaction diameters. I. INTRODUCTION
With recent advances in nanotechnology and our abil-ity to produce nano-structured materials, there is nowconsiderable interest in understanding the transport offluids inside them [1–3]. Fluids in highly confined sys-tems such as carbon nanotubes [4], zeolites [5–7], metal–organic frameworks [8] and ion channels in biologicalmembranes [9, 10] exhibit single-file diffusion (SFD) inwhich particles cannot pass each other and the diffusionis constrained in one dimension. This geometric restric-tion, combined with the nature of the fundamental par-ticles dynamics, can have a significant effect on the waythe particles move. For example, the transport coeffi-cient of a fluid along an arbitrary axis, D x , can generallybe calculated from the long time behavior of the meansquared displacement (MSD), which is given by the Ein-stein relation: (cid:68) ( x ( t ) − x (0)) (cid:69) = (cid:10) ∆ x ( t ) (cid:11) ∝ αl ( D x t ) γ , (1)where α depends upon the distribution of jumps in thebasic motion, l is the “free volume” along the x axisper particle and γ is a variable [11]. In a fluid perform-ing deterministic dynamics, the particles exhibit ballisticmotion at short times, where γ = 2, before caging andcage hopping leads to a crossover to diffusive motion atlong times where γ = 1, and Eq. 1 becomes, (cid:10) ∆ x ( t ) (cid:11) = 2 D x t . (2)Equation 2 is valid in the bulk at long times, independentof the nature of the particle dynamics, and remains truein single file systems when the dynamics is determinis-tic [12]. However, if the single file particles are subject toa Brownian background [11], or independent stochasticforces [13], then the MSD of a tagged particle is describedby an Einstein-like relation, (cid:10) ∆ x ( t ) (cid:11) = 2 F x t / , (3) ∗ m [email protected] † [email protected] where γ = 1 /
2, indicating the system exhibits anoma-lous sub-diffusion, and the diffusion coefficient has beenreplaced by the mobility factor F x [14, 15].One of the most challenging tasks in the field of con-densed matter is to find a relationship between the trans-port coefficients of a fluid and its equilibrium thermody-namics properties. Rosenfeld [16] originally proposed ascaling relation, D R ≈ .
58 exp( As ex ) between a reducedtransport coefficient, D R = DT − / ρ / ( ρ is numberdensity), and the excess entropy per particle, s ex , where A is a material dependent constant. This was later ex-panded to consider dilute fluids [17]. Dzugutov [18] de-veloped a similar scaling relation, D D ≈ .
078 exp( s ex ),linking the excess entropy to a diffusion coefficient, D D = Dρ / Γ − , that is reduced in terms of an effective Enskoginter-particle collision frequency, Γ E , that includes con-tributions obtained from microscopic structural quanti-ties such as the radial distribution function. These stud-ies have subsequently formed the basis for scaling typeapproaches focused on connecting the excess entropy todynamics in a variety of systems [19–32]. Furthermore,it is important to note that entropic scaling laws, such asthe Adam-Gibbs relation [33], have also successfully de-scribed diffusion and structural relaxation in terms of theconfigurational entropy in supercooled liquids [34, 35], atleast at temperatures above the observed break down inthe Stokes-Einstein relation [36].In the context of confined fluids, Mittal et. al. [37]showed that the scaling relationships between the self-diffusivity, excess entropy and density remain valid forboth quasi-two dimensional (particles confined betweenplates) and quasi-one-dimensional fluids. The same au-thors also found that self diffusivity of the confined sys-tems could essentially be mapped onto the behaviour ofbulk fluids through either an effective packing fraction orthe excess entropy [38]. Similarly, it has been shown thata slightly modified version of the Dzugutov scaling couldaccount for the diffusivity in water and water-methanemixtures both in the bulk and when confined to carbonnanotubes of different diameters [39]. This suggests thatthe properties of bulk materials could be used to predictthe dynamics of highly confined fluids. a r X i v : . [ c ond - m a t . s o f t ] J a n However, rather than use the bulk as a reference,Percus and co-workers [40–42] have explored the possi-bility of mapping the dynamic properties of quasi-one-dimensional fluids to those of a truly one-dimensional(1 d ) system where it is often possible to obtain exactresults, or is relatively easy to perform highly accuratenumerical simulations. In particular, Percus [41] devel-oped a heuristic approach that suggests, (cid:10) ∆ x ( t ) (cid:11) = K (cid:104)| ∆ x ( t ) |(cid:105) , (4)where (cid:104)| ∆ x ( t ) |(cid:105) contains the details of the particle mo-tion in the strictly one-dimensional system, so is propor-tional to t for inertial dynamics and proportional to t / in the case of anomalous sub-diffusion. For a 1 d systemof N hard rods with diameter, σ , on a line of length L , K = ( L/N − σ ), but more generally K could be approx-imated by kT /P , or ( − kT ∂L/∂N P ) / , where the P isthe pressure in the quasi-1 d system, T is the temperatureand k is the Boltzmann constant. The proposed strategyis then to identify characteristics of the distribution ofnext nearest neighbour separations, to which P is con-nected, that can be used to map the dynamics betweenthe system of interest and the 1 d reference system.The goal of the current paper is to simply demonstratethat an effective density and effective particle diameter,obtained either from the equation of state or a direct ge-ometric measure, is able to provide such a map betweenthe long time dynamics of a strictly 1 d hard rod systemand the dynamics of a quasi-1 d system. Here, we willfocus on the study of two dimensional hard discs con-fined to a structureless channel exhibiting deterministicdynamics as an example. The remainder of the paper isorganized as follows: Section II describes the details ofthe model and the simulation method used to examinethe dynamics. Section III outlines the calculation of theeffective density and particle size, and the mapping of thedynamics of the 1 d hard rod system. We also show thatthe same scaling appears through the use of the excessentropy. Our conclusions are summarized in Section IV. II. MODEL AND METHODS
We consider a system of two-dimensional hard disks ofdiameter, σ , confined between parallel hard lines sepa-rated by a distance H d . The channel length in the lon-gitudinal direction is L and the two ends obey periodicboundary conditions. The occupied volume fraction forthe system is then given by, φ = N πσ H d L , (5)where N is the number of particles. The particle–particleand particle–wall interaction potentials are given by, V ( r ij ) = (cid:26) r ij (cid:62) σ ∞ r ij < σ , (6) t < ∆ x > Slope ~ 1Slope ~ 2 φ =0.1 φ = 0.5 FIG. 1. The MSD along the direction of the pore axis, asa function of t for H d = 1 . φ . The short linesegments have slopes as indicated and are included as a guideto the eye. V w ( r i ) = (cid:26) | r y | (cid:54) h / ∞ otherwise , (7)respectively, where r ij = | r j − r i | is the distance betweenparticles, r y is the component of the position vector fora particle perpendicular to the wall and h = ( H d − σ ),is the reduced channel diameter.In the current work, we restrict our analysis to sys-tems performing inertial dynamics, so in the long timelimit the diffusion coefficient is given by Eq. 2. TheMSD is obtained from event-driven molecular dynamics(EDMD) simulations [43] in the canonical ( N, V, T ) en-semble where N = 30000 particles, V = h L is the totalavailable volume of system accessible to the center of theparticles and T is the temperature. The units of time inthe simulation are σ (cid:112) m/kT , where m is the mass of theparticles, which is set equal to unity, and all lengths arein units of σ . We study channels with widths H d = 1 . (cid:112) /
4) which ensures that only nearest neigh-bors can interact and prevents the particles from passingeach other.At the start of each run, the particles are placedon a lattice at packing fraction φ = 0 . kT = 1. At each φ studied, 200 N collisions are usedto reach equilibrium and the MSD in the longitudinaldirection of the channel ( x -axis) was measured over thenext 400 N collisions in which particle coordinates weresaved 80 times, separated by 5 N collisions. After col-lecting data at a given φ the system is compressed to ahigher occupied volume fraction using a modified versionof the Lubachevsky and Stillinger [44] (LS) algorithmthat ensures H d /σ remains constant as the diameter ofthe discs is changed ( L fixed), with a compression rate ds = dσ/dt = 0 . ln( t ) γ FIG. 2. Exponent, γ , as a function of ln t for H d = 1 . φ . We monitor the self-diffusivity of the fluid by fittingthe long-time behavior of the MSD in the longitudinaldirection of the channel ( x -axis) for the particles to theEinstein equation (Eq. 1). Figure 1 shows the MSD asa function of time for the case with H d = 1 . φ = 0 . φ = 0 . γ by taking the numericalderivative of the data presented in Fig. 1 to determinewhere we should evaluate the slope to obtain the diffu-sion coefficient. Figure 2 shows γ evolves over time. Atlow densities, γ decreases slowly from its ballistic valuebecause the particles are well separated. At higher densi-ties, the expected diffusive behaviour is obtained quickly,but at longer times, correlations build up as a result offinite size effects [12] and a non-monotonous behavior isobserved. Figure 3 shows the value of D x , obtained bymeasuring the slope of the MSD over the time periodwhere γ ≈
1, as a function of occupied volume fractionfor different channel diameters.The equation of state for this quasi-1 d system is ob-tained following the transfer matrix analysis developedby Kofke and Post [45] and is shown in Figure 4. III. SCALING DYNAMICS WITH AONE-DIMENSIONAL EFFECTIVE DENSITYAND PARTICLE SIZE
The exact equation for diffusion of hard rods ina strictly one-dimensional (1 d ) system was solved byJepsen [46] and is given by: D x σ = (1 − φ ) φ (2 πβm ) / , (8)where β = ( kT ) − , m is the mass of a particle, and theoccupied volume fraction becomes φ = N σ/L in 1 d . Fig-ure 3 compares D x for the quasi-1 d systems with the φ -2 -1 D x FIG. 3. Diffusivity, D x as a function of occupied volumefraction, φ , for various channel widths. f P V / N k T f q* f P V / N k T FIG. 4.
P V /NkT vs. φ for different channel H d , obtainedusing the transfer matrix method described in Ref [45]. Inset:The equilibrium fraction of defects, θ ∗ , in the jammed statesassociated with the inherent structure basins sampled by theequilibrium fluid as a function φ for channel widths H d = 1 . H d = 1 .
86 (dashed line). strictly 1 d system as a function of occupied volume frac-tion for different channel diameters. At low densities, D x for the quasi-one-dimensional channels generally fol-lows the trend of the strictly 1 d system, but decreaseswith increasing H d at fixed φ . At higher densities, thediffusion coefficients begin to exhibit behaviour that re-flects the differences between the two. For the narrowerquasi-1 d channels, D x decreases more rapidly than thestrictly 1 d system because D x necessary goes to zero as φ approaches φ J max <
1, the most dense packing for thesystem. For the wider channels, D x exhibits a plateauat intermediate densities that coincides with the plateauobserved in the equation of state (see Fig 4). Whilethe short ranged hard disc interaction and quasi-one-dimensional nature of the system rule out the possibilitythat the system has a phase transition, the fluid under-goes a significant, but continuous, structural rearrange-ment at intermediate densities [45, 47]. The low densityfluid is characterized by linear arrangements of the discswhere the collisions between discs occur predominantlyin the longitudinal direction, and the high density fluid isdominated by zig-zag arrangements of the discs, with discinteractions occurring across the channel. In the transi-tion region, the effects of increasing φ are somewhat offset by these structural rearrangements so the amount ofspace the discs have to move in the longitudinal directionchanges slowly, giving rise to the plateaus in the diffusioncoefficient and the longitudinal pressure. The diffusioncoefficient would be expected to decrease again at higherdensities where the pressure begins to increase again, butour simulations becomes trapped in a glassy state whereit is not possible to measure the equilibrium propertiesof the fluid.The goal of this work is to identify a scaling thataccounts for these differences and collapses D x for thequasi-1 d systems onto Eq. 8. As the diameter of channelchanges, the kinetic term in Eq. 8, (2 πβm ) − / remainsunchanged, so the focus of our task is to provide a mapbetween the quasi-1 d and strictly-1 d systems for the den-sity dependent term, (1 − φ ) /φ , which can be thoughtof as the mean clearance between neighbouring parti-cles [12]. One possibility is to explore this term in thecontext of the inherent structure paradigm [48, 49], whereconfigurations that map to the same mechanically stablejammed structure, or inherent structure, are grouped to-gether into local basins of attraction so that the thermo-dynamics and dynamics of the fluid can be described interms of the basins sampled at equilibrium. The amountof vibrational space the particles have when samplinga given basin then goes as ( φ J − φ ). The strictly-1 d hard rod system has an extremely simple inherent struc-ture landscape consisting of a single basin that jams in astate with all the rods in contact end-to-end and φ J = 1.The inherent structure landscape is more complicated inhigher dimensions. In two dimensions, a disc is locally lo-cally jammed if it has at least three rigid contacts, not allin the same hemisphere. However, a configuration madeup of locally jammed discs is not necessarily a mechan-ically stable inherent structure, or collectively jammedstate, because a collective motion of the particles canlead to unjamming [50], highlighting the global nature ofjamming phenomena.Confining the particles within narrow quasi-1 d chan-nels reduces the number of possible disc-disc contactsand prevents collective rearrangements of the particles.As a result, the complete inherent structure landscapeof the current model, and how it is sampled as functionof φ , can be constructed by considering local packing ar-rangements [51–53]. The most dense structure consistsof a zig-zag packing of discs contacting their two neigh- H d L h FIG. 5. A configuration of hard discs in a jammed packingcontaining discs in the most dense (red) and defect (green)local structures. bours across the channel as well as the wall and the lowerdensity jammed states are formed by introducing defectswhere discs form a contact with another disc along thechannel (see Fig. 5). The jamming density of a packingis then given by [47], φ J = π H d (cid:16) θ + (1 − θ ) (cid:112) H d (2 σ − H d ) (cid:17) , (9)where θ is the fraction of defects so that φ J → φ J max as θ →
0. Noting that two defects cannot appear next toeach other in a jammed configuration because the parti-cles would not satisfy the local jamming condition, it hasalso been shown that the lowest density jammed struc-ture consists of alternating dense and defect states, with θ = 0 .
5. The number of inherent structure basins is alsoa function of θ , with the maximum number of jammedstructures occurring when θ = 1 / − √ /
10 and decreas-ing to unity as θ approaches its upper and lower limits.At equilibrium, the fluid samples the set of basins onthe inherent structure landscape that maximizes the totalentropy so that at a given φ the system samples basinswith the equilibrium fraction of defects θ = θ ∗ . Theinsert to Figure 4 shows θ ∗ as function of φ , obtained bymaximizing the total entropy of the fluid [52, 53]. At low φ , the system samples the region of the inherent structurelandscape where there is a large number of basins with ahigh defect fraction and low φ J . With increasing φ , thesystem samples a smaller number of basins, and smallerdefect fraction, that have a greater vibrational entropydue to their larger φ J , effectively trading a reduction ofconfigurational entropy for increased vibrational entropy.In the context of the current problem, this means that φ J is a function of density, which in turn influences theamount of space the particles have move. Equation. 8can then be rewritten to account for density dependentjamming density as, D x σ = ( φ J − φ ) φ (2 πβm ) / , (10)where ( φ J − φ ) /φ is the effective mean clearance betweenadjacent particles, and φ J at each density is given byEq. 9 using θ ∗ . Figure 6 shows the predictions of theEq. 10 for different channel widths, where we have usedthe equilibrium value of θ to obtain φ J . While the data issomewhat linearized the slopes of the lines decrease with ( f J - f )/ f ( D x / s ) ( pb m ) / FIG. 6. Test of Eq. 10 plotting ( D x /σ N )(2 πβm ) / as a func-tion of ( φ J − φ ) /φ . increasing H d , which may be a result of the increasedmotion of the particles perpendicular to the longitudi-nal axis of the channel that is not accounted for in ourmeasurements of the MSD.The exact EOS for quasi-1 d systems was solved byKofke and Post [45] using a transition matrix approach.Their approach uses the idea that when the y positions ofthe particles are fixed the system can be characterized asa 1 d mixture of hard rods on a line with different contactlengths. This allow the integration in x to be performedindependently of the y integration when solving the par-tition function. The results of the integrals can be solvedas an eigenvalues problem where the largest eigenvalueis used, which essentially identifies an average effectivediameter of the projected hard rods.Motivated by this, we define an effective one-dimensional interaction diameter, σ i , between two neigh-bouring particles that simply depends on their relativeseparation in the y -direction, ∆ y = | y i +1 − y i | , σ i = (cid:112) σ − ∆ y , (11)and reduces to σ when ∆ y = | y i +1 − y i | = 0 (see Fig. 7).The average effective 1 d interaction diameter is thengiven by, σ N = 1 N (cid:88) i σ i , (12)and we can define an effective occupied volume fraction, φ eff = N σ N L , (13)that can be extracted from the equation of state of thequasi-one dimensional system and used to obtain the dif-fusion coefficient from Eq 8. To achieve this, we assumethe exact 1 d equation of state holds for φ eff : Z = 11 − φ eff , (14) σσ i (a)(b)(c) σ i FIG. 7. Effective diameter, a) purely 1 d system σ = σ N , b) H d = 1 . H d = 1 . where Z = P V /N kT is the compressibility factor for thequasi-one-dimensional system. Rearranging, we obtain φ eff as, φ eff = Z − Z . (15)Sububstituting eq. 13 into 8, gives the diffusion coeffi-cient in terms of the effective volume fraction, D x σ N = (1 − φ eff ) φ eff (2 πβm ) / . (16)Figure 8 shows Eq. 16 using the simulated data forthe diffusion coefficient at different channel widths com-pared to exact analytical results for 1 d system. The plotshows, the collapse of the data is good, which suggests φ eff provides a good thermodynamic connection throughthe equation of state to 1 d diffusion, as suggested by Per-cus [41]. The deviation at low densities may be relatedto a broad distribution of effective diameters at low den-sities, where the particles move freely in the transversedirection and can adopt many different effective diame-ters. At high densities, the distribution of effective diam-eters is restricted more closely to the average, except inthe cases where the particles are appear in defect statesand σ i ≈ σ . In addition, the measurements of D x atthe very lowest densities can be difficult because of thelow frequency of particle-particle collisions and the slowconvergence of γ → φ eff is obtained from the equation of state, it canbe measured geometrically. Figure 9 compares the valuesof σ N obtained using eq. 14 with the average values of theeffective diameter measured during our MD simulationsat different channel widths. This shows the usefulness ofthis quantity especially, in the cases that the exact EOSis not available.Finally, we simply demonstrate that φ eff also providesan appropriate link between the excess entropy and dif-fusion in these confined fluid systems. Here, the excess (1- f eff )/ f eff ( D x / s N )( pb m ) / FIG. 8. Scaled self-diffusion coefficient data represented byusing Eq. 16 by using effective volume fraction for variouschannel width. s N s N ( M D ) FIG. 9. Effective diameter calculated from Eqs. 13 and 14vs. those calculated from MD simulations using Eq. 11 fordifferent H d . entropy of a strictly one-dimensional hard rod system,relative to ideal gas, can be obtained, s ex /N k = − (cid:90) φ Z − φ (cid:48) d φ (cid:48) , (17)which yields, s ex /N k = ln(1 − φ ). (18)Again, replacing φ with φ eff in Eq. 18, and using Eq. 16 yields, D x = ( σ N /φ eff ) (2 πβm ) − / exp ( s ex ) , (19)which connects the excess entropy to a reduced diffusioncoefficient, D x φ eff /σ N , in a similar form to that obtained -S ex -1 ( D x / s ) f ( pb m ) / FIG. 10. Scaled self-diffusion coefficient data versus negativeexcess entropy calculated using Eq. 17 for different H d . by Rosenfeld [16] and Dzugutov [18]. Figure 10 showsthe relationship between the diffusion coefficient and ex-cess entropy in the quasi-one-dimensional system beforerescaling. Then use of φ eff necessarily reproduces thesame collapse observed in Fig. 8. IV. CONCLUSION
Simulation studies of both bulk and confined fluidshave shown the effectiveness of scaling relationships thatmap the transport properties of the system to a refer-ence state through a thermodynamic property such asthe excess entropy. In the case of highly confined, quasi-one-dimensional fluids, the strictly one-dimensional fluidis an ideal reference state because results are often knownexactly. This work shows that the proposed scaling canbe achieved for a system of quasi-one-dimensional harddiscs through an effective occupied volume fraction thatis obtained directly from the equation of state or canbe measured geometrically. However, the channel widthsconsidered here are still narrow and it is likely the scal-ing will begin to breakdown for wider channel. It alsoremains to be seen how well the scaling works for longerrange interactions.
ACKNOWLEDGMENTS
We thank NSERC for financial support. We also thankWestGrid and Compute Canada for providing computa-tional resources. [1] R. M. Robertson, S. Laib, and D. E. Smith, Proc. Natl.Acad. Sci. U.S.A. , 7310 (2006).[2] K. E. Gubbins, Y.-C. Liu, J. D. Moore, and J. C. Palmer,Phys. Chem. Chem. Phys. , 58 (2011).[3] R. Rodr´ıguez Schmidt, J. G. Hern´andez Cifre, andJ. Garc´ıa de la Torre, Euro. Phys. J. E, , 9806 (2012).[4] B. Mukherjee, P. K. Maiti, C. Dasgupta, and A. K. Sood,ACS nano , 985 (2010).[5] J. K¨arger, M. Petzold, H. Pfeifer, S. Ernst, andJ. Weitkamp, J. Catal. , 283 (1992).[6] A. V. A. Kumar, Mol. Phys. , 1306 (2014).[7] J. K¨arger, ChemPhysChem , 24 (2015).[8] H. Jobic, Phys. Chem. Chem. Phys. , 17190 (2016).[9] M. Ø. Jensen, S. Park, E. Tajkhorshid, and K. Schulten,Proc. Natl. Acad. Sci. U.S.A. , 6731 (2002).[10] J. C. Rasaiah, S. Garde, and G. Hummer, Annu. Rev.Phys. Chem. , 713 (2008).[11] D. G. Levitt, Phys. Rev. A , 3050 (1973).[12] K. Hahn and J. K¨arger, J. Phys. Chem. , 316 (1996).[13] J. K. Percus, Phys. Rev. A , 557 (1974).[14] K. Hahn and J. K¨arger, J. Phys. Chem. B (1998),10.1021/jp981039h.[15] B. Lin, M. Meron, B. Cui, S. A. Rice, and H. Diamant,Phys. Rev. Lett. , 216001 (2005).[16] Y. Rosenfeld, Phys. Rev. A , 2545 (1977).[17] Y. Rosenfeld, J. Phys.: Condens. Matter , 5415 (1999).[18] M. Dzugutov, Nature , 137 (1996).[19] J.-L. Bretonnet, J. Chem. Phys. , 9370 (2002).[20] I. Yokoyama, High Temp. Mat. Pr. , 341 (2001).[21] Q.-L. Cao, X.-S. Kong, Y. Li, X. Wu, and C. Liu, PhysicaB , 3114 (2011).[22] Q.-L. Cao, J.-X. Shao, P.-P. Wang, and F.-H. Wang, J.Appl. Phys. , 135903 (2015).[23] Q.-L. Cao, P.-P. Wang, J.-X. Shao, and F.-H. Wang,RSC Adv. , 84420 (2016).[24] N. Jakse and A. Pasturel, Scientific Reports , 20689(2016).[25] A. Samanta, S. M. Ali, and S. K. Ghosh, Phys. Rev.Lett. , 145901 (2004).[26] M. Wautelet, Euro. J. Phys. , 601 (2001).[27] M. J. Colbrook, X. Ma, P. F. Hopkins, and J. Squire,Monthly Notices of the Royal Astronomical Society ,2421 (2017).[28] G. Gioia, G. Lacorata, E. P. M. Filho, A. Mazzino, andU. Rizza, Boundary-Layer Meteorology , 187 (2004). [29] J.-M. Bomont and J.-L. Bretonnet, Physica A , 23(2015).[30] L. Guang-Xu, L. Chang-Song, and Z. Zhen-Gang, Chin.Phys. Lett. , 2489 (2004).[31] S. Ramachandran, S. Komura, K. Seki, and G. Gompper,Eur. Phys. J. E , 46 (2011).[32] A. Bhattacharya, W. H. Morrison, K. Luo, T. Ala-Nissila, S. C. Ying, A. Milchev, and K. Binder, Eur.Phys. J. E , 423 (2009).[33] G. Adam and J. H. Gibbs, J. Chem. Phys. , 139 (1965).[34] S. Sastry, Nature , 164 (2001).[35] R. J. Speedy, J. Chem. Phys. , 9069 (2001).[36] A. D. S. Parmar, S. Sengupta, and S. Sastry, Eur. Phys.J. E , 90 (2018).[37] J. Mittal, J. R. Errington, and T. M. Truskett, J. Phys.Chem. B , 18147 (2006).[38] J. Mittal, J. R. Errington, and T. M. Truskett, J. Chem.Phys. , 4708 (2007).[39] P. Sahu, S. M. Ali, and K. T. Shenoy, J. Phys. Chem. C , 11968 (2017).[40] P. Kalinay and J. K. Percus, J. Chem. Phys. , 204701(2005).[41] J. K. Percus, J. Stat. Phys. , 40 (2010).[42] P. Kalinay and J. K. Percus, Phys. Rev. E , 031109(2011).[43] B. J. Alder and T. E. Wainwright, J. Chem. Phys. ,1439 (1960).[44] B. Lubachevsky and F. H. S. Jr, J. Stat. Phys. , 561(1990).[45] D. A. Kofke and A. J. Post, J. Chem. Phys. , 4853(1993).[46] D. W. Jepsen, Journal of Mathematical Physics , 405(1965).[47] R. K. Bowles and I. Saika-Voivod, Phys. Rev. E ,011503 (2006).[48] F. H. Stillinger, E. A. DiMarzio, and R. L. Kornegay, J.Chem. Phys. , 1564 (1964).[49] F. H. Stillinger and T. A. Weber, Phys. Rev. A ,978989 (1982).[50] S. Torquato and F. H. Stillinger, J. Phys. Chem. B ,1184911853 (2001), doi: 10.1021/jp011960q.[51] M. Zaeifi Yamchi, S. S. Ashwin, and R. K. Bowles, Phys.Rev. Lett. , 225701 (2012).[52] S. S. Ashwin, M. Zaeifi Yamchi, and R. K. Bowles, Phys.Rev. Lett. , 145701 (2013).[53] M. Zaeifi Yamchi, S. S. Ashwin, and R. K. Bowles, Phys.Rev. E91