The role of local bond-order at crystallization in a simple supercooled liquid
TThe role of local bond-order at crystallization in a simplesupercooled liquid
Søren Toxvaerd
DNRF centre “Glass and Time,” IMFUFA,Department of Sciences, Roskilde University,Postbox 260, DK-4000 Roskilde, Denmark (Dated: February 9, 2021)
Abstract
Large scale Molecular Dynamics simulations of sixty-five systems with N =80000 Lennard-Jonesparticles at two different supercooled liquid state points reveal, that the supercooled states containspatially heterogeneous distributed subdomains of particles with significant higher bond-order thanthe mean bond-order in the supercooled liquid. The onset of the crystallization starts in such anarea with relative high six-fold bond-order for a supercooled state, but low bond-order for a fcccrystal state, and the crystallization is initiated by a nucleus where all particles in the criticalnucleus have a significant lower bond-order than in a fcc crystal. The critical nucleus of N ≈ a r X i v : . [ c ond - m a t . s o f t ] F e b IG. 1. The liquid-solid phase diagram for a LJ system. The liquid states, which coexist with gasare with light green. The liquid states, which coexist with solid fcc are red, and the coexisting fcccrystal states are blue. The MD systems are quenched from the liquids (red points) down to thesupercooled liquids (blue points).
I. INTRODUCTION
Computer simulations have allowed us to determine the dynamics of crystallization. Alderand Wainwright’s [1] pioneering simulations of systems of hard spheres revealed, that a sys-tem of hard spheres at high packing fraction crystallizes into an ordered state with fcclattice structure. Later, computer simulations of a system of particles with the more real-istic Lennard-Jones (LJ) potential [2] not only confirmed the hard sphere result, but theywere also in agreement with crystallizations of noble gases at low temperatures. Many latersimulations of simple systems verified the result, that the minimum free energy of a con-densed system of simple spherical symmetrical particles is for the fcc lattice arrangement[3–6].The theory of crystallization is usually described by the classical nucleation theory (CNT)and its extensions [7–9]. But an exact analyse [10] and many computer simulations ofnucleation of liquids as well as crystals in supercooled gases exhibit a much more complexdynamics with polymorphism [11], than given by CNT. The critical nucleus in a gas phaseis not a compact object of the new phase, and it is not only initiated by density fluctuations,but also by temperature fluctuations [12, 13]. A recent review of theory and simulations ofcrystallizations is given in [14].The supercooled liquid exhibits bond orientational order [15], and here we show that the2
IG. 2. Energy per particle, u(t), as a function of the time after cooling to the supercooled state(
T, ρ )=(1.25,1.10) (left blue point in Figure 1). The figure shows six representative examples of u ( t ). The inset are the energies for all twenty five experiments in the time interval ∆ t ∈ [4000 , bond-order, given by Q [16], of the particles in a supercooled LJ liquid is heterogeneousdistributed. The supercooled liquid contains big areas with significant higher bond-orderthan the mean bond-order in the liquid, and the crystallization is initiated from such anarea. This result is consistent with the well known dynamic heterogeneities in supercooledliquids [17–19].The crystallization is very fast and accompanied by a growth of the supercooled areas withrelative high liquid bond-order. The growing crystal nuclei have fcc bond-order, but in manycases the nuclei also exhibit five fold symmetry with polycrystalline domains. The crystallinedomains are mainly with fcc crystal arrangements, but some are also with hcp structure. Hisbehavior confirms previous results for homogeneous crystallizations in hard sphere systems[20, 21] and in a LJ system [11, 22–24]. Many of the polycrystalline nucleations ended inlong time metastable polycrystalline states. II. HOMOGENEOUS CRYSTALLIZATION
Molecular Dynamics (MD) systems of 80000 LJ particles (see Appendix) at different liq-uid state points were cooled down to the supercooled states (Figure 1). We have performedtwenty five independent
N V E simulations of supercooling and crystallizations at the state3
IG. 3. Probability distribution of the bond-order parameter Q . The red curve is for the liquidbefore the cooling (left red point in Figure 1) and the blue curve is for the supercooled liquid (leftblue point in Figure 1). The green curve is for a perfect fcc crystal at the supercooled state point.The inset shows the log -distributions. point ( T, ρ )=(1.25,1.10), and twenty
N V E and twenty
N V T simulations at the state point(
T, ρ )=(2.80,1.30). The degree of supercooling given by the ratio between the supercooledtemperature T and the freezing point temperature of the liquid T f is 0.54 and 0.60, respec-tively and the systems crystallized at the time ∆ t ≈ −
150 after the quenches. (Units aregiven in LJ units, see Appendix 1, and the time unit ∆ t = 1 corresponds to 1000 MD times.jpg.)The crystallizations result in a decrease in pressure and energy and were essentiallycompleted within a crystallization time of ∆ t ≈
50. But for some of the simulations thecrystallizations were, however, first completed after longer times. The energies per particlefor the
N V E crystallizations at (
T, ρ )=(1.25,1.10) are shown in Figure 2, and the energiesfor the twenty
N V T crystallizations are shown in Figure 9.The description of the dynamics of the crystallization is divided into three subsections:The description of the supercooled liquid, the onset of crystallizations, and finally the crys-tallization and the description of the crystal states after the crystallizations were completed.
A. The supercooled liquid
The bond-order function Q (defined in Appendix 2) in the supercooled liquid is usedto detect the dynamics of crystallization. The liquid state is characterized by having a low4 IG. 4. Side view of the box with the N particles. value of bond-order compared with the bond-order in the crystal state. In [25] the authorsused bond-order in the Gaussian core model, which is a prototype for soft spheres, to analysethe onset of crystallization. They found, that the crystallization occurs in precursor regionsof high bond orientation order, and that the crystal which first nucleates is the one which hasthe closest symmetry to the ordered regions in the supercooled state. A later investigationof the bond-order in a compressed hard sphere fluid found, however, that the hexagonalordering appeared simultaneously with the density fluctuation at the onset of crystallization[26].The present investigation shows, that the supercooled liquid contains large subregions ofparticles with relative high Q bond-order, and the nucleus which initiates the crystalliza-tion has only a bond-order which is somewhat higher than the order in the heterogeneousdistributed bond-order domains, but on the other hand have a lower bond-order than in thecrystal.The distributions of the order parameter Q ( i ) for particles i in the liquid state ( T, ρ ) =(4 . , .
10) (left red point in Figure 1) and in the supercooled state (
T, ρ )=(1.25, 1.10)(blue point) are shown in Figure 3 together with the distribution in a fcc crystal state at(
T, ρ )=(1.25, 1.10). The distributions for the supercooled state (blue curve) and for the fcccrystal state (green curve) are separated, but the log distribution of Q ( i ) in the inset ofthe figure reveals, that there is an exponential decreasing probabilities for low order in thecrystal state and a high order in the supercooled state, and that there is an overlap of the5 IG. 5. Mean numbers of clusters n N in the liquid states for clusters with N particles withbond-order Q ( i ) > .
25 for all particles i in the clusters. The figure shows the log ( n N ( log (N))distributions and the inset is the log ( n N )(N) distributions. Red curves are for the liquid (left redpoint in Figure 1) and blue curves are for the supercooled liquid. The straight black line in thefigure is an algebraic fit, a × n b N , to the distribution in the supercooled liquid for clusters in theinterval n N ∈ [20 , t ), in one of the twenty five simulations, of particles i in the biggest clusterof crystal-ordered particles with 0 . < Q ( i ). The onset of crystallization at t cr =56 is marked byan arrow. From there the crystal nucleus grows quite monotonically. The green line is an estimateof the critical crystal cluster size of N cr. ≈
70. The inset shows the x-, y- and z-positions of thebiggest cluster. two distributions in the region 0 . < Q < .
45. Investigation of the crystal ordering atthe creation of the critical nucleus (next subsection) shows, that the successful nucleus havea mean order Q > .
35 in accordance with the distributions in Figure 3. The distributionsof Q for the liquid and the supercooled liquid, shown in Figure 3 are, however, different6 IG. 7. The particle positions at t = 0 .
58 in the system with the N ( t ) shown in the previousfigure. The local environment with the critical nucleus is shown in the next figure. Color as inFigure 4. although they have the same mean Gaussian-like shapes. An analysis of the distribution ofparticles with relative high bond-order reveal this fact.From the locations of particles with different values of Q ( i ) for the particles i in thesupercooled liquid one can see, that the distribution is non-uniform, and that there existsbig subdomains with relative high values of Q ( i ) for all the particles in the domain. Thenext figure shows this. Figure 4 is a side view of the positions at a certain time of particles i in the supercooled liquid. The white transparent spheres are particles with Q ( i ) < . . < Q ( i ) < .
35, blue spheres have 0 . < Q ( i ) < .
40, and the redspheres are particles in the supercooled liquid with a value of the bond-order 0 . < Q ( i ).The blue and red particles are particles with a lattice order, which is sufficient for crystalnucleation. (The values of Q are obtained as time averages over thousand time s.jpg, butthe heterogeneous distribution is also obtained from shorter and longer time intervals.) Thepositions of the particles is not uniformly distributed, but contains large areas with relative7igh bond-order.The cluster distribution of particles in the system with a certain quality, e.g. a high Q value, can be obtained directly during the MD simulation and without a significant increasein computer time for the big system by using the nearest neighbour list [12]. The clustersof N particles with 0 . < Q ( i ) and the mean number n N of clusters with N particleswas determined directly during a MD simulation. The number n N (N), of a cluster with Nparticles with a bond-order 0 . < Q within a time interval δt = 0 . Q ( r i ( t )) of the bond-order for a particle i at position r i ( t ) in the supercooledliquid state fluctuates with time, and the values of n N in the figure are the mean for 200independent distributions of clusters with Q ( i ) > .
25, and where Q ( i ) is obtained as themean bond-order of a particle i within the time interval δt = 0 . ≈ one mean vibrationwithin the shell of nearest neighbours). The figure shows the mean number n N of clustersof N particles with relative high bond-order Q ( i ) in the liquid state ( T, ρ )=(4.25,1.10) (redcurve) and the corresponding distribution of clusters in the supercooled state (blue curve).The two distributions are essentially different. The figure shows log ( n N ) as a function of log (N) and the inset gives log ( n N ) as a function of N. The distribution in the liquid state isexponentially declining (red ≈ straight line in the inset), whereas the distribution for biggerclusters in the supercooled liquid is algebraic (blue ≈ linear function and black straightline in the figure). The black straight line is a line obtained from a fit of a × N b to the log ( n N ( log N)) distribution in the interval n N ∈ [20 , n N shows alsothat there is regions of many particles, as can be seen in Figure 4, with relative high liquidbond-order Q > . B. The onset of crystallization
The crystallization in a supercooled liquid appears when an ensemble of particles withlattice order gain free energy by increasing its size. In the CNT by, that the gain in freeenergy of the crystal phase exceeds the cost in surface free energy by the increased surface8f the crystal. In the MD ensemble simulation one primarily observes the crystallization,and the method gives not a direct information about the free energy. For this reason it isnot possible to determine the critical nucleus precisely. But one can locate the successfulnucleus and its environment at the onset of nucleation.The onset of crystallization is determined from the growth of the biggest cluster withbond-order 0 . < Q . In the supercooled liquid the number N c of particles in the biggestcluster with this crystal-like bond-order fluctuates with N c ≤
100 (Figure 6), but from theonset of crystallization the biggest cluster grows very fast, as shown in the figure. The x-,y- and z- positions of the center of mass are shown in the inset. The time record of thesepositions identify the position of the successful nucleus even before it reaches the criticalnucleus size. The rather constant location of the center of mass of the biggest cluster attimes t ≥
56, and the fluctuating locations before t = 56 show, that different crystal likenuclei appear in the system for t <
56, but from then on the crystallization is initiated bythis nucleus. The size of the successful nucleus at t = 56 is, however only N c = 17 andclearly much smaller than the critical nucleus. But after only 2000 time s.jpg the centerconsists of 68 particles with 0 . < Q . The green line N = 70 is an estimate of the criticalnucleus size, based on inspection of the twenty-five N V E simulation. As can be seen fromFigure 6 some nuclei at times t <
56 before the crystallization sometime grow to this size.This behaviour before the onset of nucleation is a typically for all twenty five systems in thesupercooled states.The positions of the particles at t = 58 is shown in Figure 7 and with with the samecolor as in Figure 4. The critical nucleus is enlarged in the next figure. Figure 8 shows thepositions of the 574 particles within a sphere with the center at the center of mass of thecritical nucleus and with the radius 5. The critical nucleus consists of 59 (blue) particleswith 0 . < Q < .
40 and 9 (red) particles with 0 . < Q < .
45 (no particle have a0 . < Q ). There is 223 particles (green) with 0 . < Q < .
35 and the particles with Q < .
25 are white transparent. The mean bond-order of the 68 particles in the crystal-likecritical nucleus is < Q > =0.38. The green particles are mainly located in a shell aroundthe critical nucleus and with a tendency to fit into the lattice planes of the critical nucleus.Even some of the white particles with lower bond-order are located in the critical crystalsplanes.The crystal bond-order < Q > =0.38 in the critical nucleus is, however significant lower9 IG. 8. The local environment of 574 particles within the sphere with radius 5 and with the centerat the center of mass of the of the critical crystal nucleus shown in Figure 7. It consists of N c = 68particles with 59 (blue) particles with 0 . < Q ( i ) < .
40 and 9 (red) particles with 0 . < Q ( i ).The 225 small green particles have 0 . < Q ( i ) < .
35 and the small white transparent partparticles have Q ( i ) < . than the bond-order in a bulk fcc crystal, and the twenty N V E and the twenty
N V T simulations at the higher temperature and density (
T, ρ ) = (2 . , .
30) show the sametendency. All the critical crystal nuclei have a significant lower crystal Q value than aperfect fcc crystals.The means for the twenty-five crystal critical nucleus areMean bond-order in the critical nucleus < Q ( i ) > = 0 . ± . cr = 73 ± IG. 9. Potential energy per particle, u(t), as a function of time at the onset of crystallizationfor the (
T, ρ ) = (2 . , .
13) (right blue point in Figure 1). The inset shows u(t) at the end of thesimulations.
C. The crystal states
Some crystallizations were completed within a short time-interval of 50 to 100 time units,but it took much longer times for many of the crystallizations as shown in Figure 2. More-over, the systems did not end up in the same crystal state with energies close to the energy ofa perfect fcc crystal. Some of the systems ended up in states with significant higher energiesand pressures. For the twenty-five systems at the relative low density two to three of thesimulations ended up in states with significant higher energies than the other systems. Thistendency is more pronounced for the twenty
N V E and the twenty
N V T simulations at thehigher temperature and density state (
T, ρ ) = (2 . , . N V T simulations after thesupercooling, and with the energies at the end of the simulations in the inset of the figure.The energies exhibit two energy bands, with fourteen of the energies a little above the energyfor a perfect fcc crystal (black line), whereas there are six simulations that have significanthigher energies. The same result was obtained for the twenty
N V E simulations at thesame state point. From inspection of the particle positions and from the radial distributionfunctions it is is clear, that the fourteen systems in the lower energy band are fcc latticeswith some defects, but the crystal states in the systems with energies in the upper energyband is more complex.The structures of the systems with energies in the high energy band are investigatedin order to determine their lattice structure. Particles interacting with simple spherical11
IG. 10. Side view of the particles in the system with the highest energy after the crystallization( upper red curve in Figure 9). The particles i are colored as in Figure 4. With: low order with Q ( i ) < . < Q ( i ) < .
30; blue: 0 . < Q ( i ) < .
45; red: 0 . < Q ( i ). symmetrical pair-potentials can exist in many crystal arrangements [5]. But for a LJ systemat the two state point investigated here the fcc crystal structure has the lowest free energy[5, 6]. The positions of particles for the system with the highest energy in Figure 9 (redupper energy function in the inset of Figure 9 is shown in Figure 10. The particles are coloredin accordance with their bond-order, and the colored positions show a complex crystallinestructure with many areas with relative low bond-order and even without crystalline order( Q < . g ( r ). The extremes in the radial distributionfunctions at the high temperature T = 2 .
80 are, however, not sharp due the particlesvibrations at their lattice positions. This thermal “noise” can be removed by cooling thesystem down to a low temperature. Figure 11 shows the radial distribution functions for thesystem with lowest energy (blue) and highest energy (red) in Figure 9, and the distribution g fcc ( r ) for a perfect fcc lattice (green), and after the systems were cooled down to T =0.1.The figure confirms, that both systems mainly consist of particles in a fcc arrangement.12 IG. 11. The radial distribution functions for the system with the lowest energy (blue) and withthe highest energy (red) in the twenty
N V T simulations. After the crystallizations the systemswere cooled down to T = 0 . t = 40 afterthe onset of nucleation. The color of the spheres are the same as in Figure 10, but particles with0 . < Q are (white) transparent. The crystal cluster contains 5516 particles and the periodicboundary plane in front goes through the crystal nucleus. g fcc ( r ) especially atthe radial distance r ≈ .
9. The inset, where the two functions are compared with g hcp ( r )(green) for a hcp crystal reveals, that the crystal (Figure 10) contains a small number ofparticles with hcp structure. The system with the highest energy was simulated over a verylong time (47000 time units) in order to test whether it is stable, and it continued to be inthe polycrystalline state.The polycrystalline state is ensured in the beginning of the crystallization, and it is not afinal size effect of the periodical boundaries. The next figure shows the particle positions 40time units after the onset of crystallization. The crystal nucleus consists of 5516 particlesand the periodic plane in front goes through the crystalline nucleus and reveal, that thenucleus consists of several small crystals with different orientations of the crystal planes.The green particles in the supercooled liquid with relative high bond-order percolate thesystem and a new crystal center (upper left) is created with green particles included in thelattice arrangement. The figure shows that the spontaneous crystallization is preformed alsoby growing coherent order in the supercooled liquid (green particles) , caused by the criticalnucleus. This behaviour of the spontaneous crystallization is found in all the crystallizations. III. CONCLUSION
A supercooled liquid exhibits spatially heterogeneous dynamics [17], which influences thedynamic behaviour of the supercooled and viscous liquid. Here it is determined, that theparticles with relative high bond-order in the supercooled Lennard-Jones system is hetero-geneous distributed, and that the supercooled liquid contains subdomains with significanthigher bond-order than the mean bond-order in the supercooled liquid. Furthermore, thecrystal nucleation is initiated from such region with relative high bond-order Q for a su-percooled liquid. The domains with excess bond-order changes in extension and locationswith time, and this is ”consistent” with the dynamic heterogenities for supercooled liquid,but we have not directly established the connection between the domains of bond-order andthe dynamics of the supercooled liquid, e.g. by determining the mobility of the particles inthe subdomains.The nucleus which initiates the crystallization is a relative compact object with sixfoldsymmetry as assumed in the classical nucleation theory (CNT), but all the N ≈
70 particles14n the nucleus have, however, a significant smaller bond-order than the bond-order in the fcccrystal. And in addition, many of the surrounding particles in the subdomain are alignedwith the particles in the initiating crystal nucleus (Figure 8).The growth of the critical nucleus is fast, and the crystal has percolate the big systemof 80000 LJ particles within a crystallization time of ≈
50 to 100 time units. The systemsdid, however, not always end up in a homogeneous fcc state, but quite often they ended upin a polycrystalline state. The polycrystalline state with with traces of hcp crystallites is,however, established in the crystal nucleus shortly after the onset of crystallization (Figure12). The polymorphism, where the particles crystallizes into different structures have alreadybeen observed in LJ systems [28, 29]. Also here the areas with relative high bond-order inthe supercooled part of the system plays a role for the fast spontaneous crystallization.The heterogeneous distributed areas grow fast and percolate the system long times beforecrystallization.The systems have been simulated by Molecular Dynamics
N V E and
N V T and with thesame qualitative results. This is perh.jpg not surprising, because there is only a marginaldifference between the two MD methods for the big system. This fact is due to, that theparticles in the
N V T dynamics are constrained to a thermostat temperature T th by theexcess of kinetic energy of the hole system . The systems mean temperature T in N V T aresmoothly constrained over longer times to the
N V T (cid:48) s T th value, which removes the latentheat during the spontaneous crystallization. . IV. ACKNOWLEDGMENT
Ulf R. Pedersen, Trond S. Ingebrigtsen and Jeppe C. Dyre is gratefully acknowledged.This work was supported by the VILLUM Foundation’s Matter project, grant No. 16515.
Author contribution statement
The articles idea, programs used in the simulations, the simulations and the writing aredone by the author.
V. APPENDIX
1. The MD system and the computational details N = 80000 Lennard-Jones (LJ) particles in a cubic box withperiodic boundaries, and the crystallizations are obtained by Molecular Dynamics N V E simulations with Newton’s central difference algorithm [30] (“Leap-frog”). The time andlength are given by the length unit l ∗ = σ and energy unit u ∗ = .jpgilon/k B in the Lennard-Jones potential for particles with the mass m . The unit of time is t ∗ = σ (cid:112) m/.jpgilon .The LJ forces are truncated and shifted at the interparticle distance r c = 2 . δt = 0 . r c [33]. Mostsimulations of LJ systems including this one are for r c = 2 .
5, for which the triple point den-sities are ( ρ l , ρ s , T )=(0.8290, 0.9333, 0.618) [34]. The present simulations of crystallizationsare for supercooled liquids at the state points ( T, ρ )=(1.25, 1.10) and (
T, ρ )=(2.80, 1.30) .The crystallizations can be characterized as crystallizations in supercooled condensed liquidstates.The crystallizations are performed by cooling from the liquid states at (
T, ρ )=(4.25,1.10)and (
T, ρ )=(5.25,1.30) down to the supercooled states (
T, ρ )=(1.25, 1.10) and (
T, ρ )=(2.80,1.30), respectively. The cooling and
N V T simulations are performed by a standard
N V T thermostat. The
N V E simulations are performed by cooling the high temperature systemdown in teen thousand time s.jpg by the thermostat. The
N V E supercooled state wasaccepted, if the temperature in the succeeding ten thousand time s.jpg without a thermostatwas within the temperature interval T ∈ . ± . N V E systems were,however, so supercooled, that they ended up in crystal states (T, ρ ) ≤ (1.63,1.10) with totalcrystallization . A LJ fcc crystal at ρ =1.1 melts at T=1.68 [35].
2. Identifying crystal structure
The crystal order is determined by a modified bond-orientation order analysis [16, 25, 36].A complex order parameter q lm ( i, t ) ≡ /n b Σ n b j =1 Y lm ( r αβ ( t )) (1)is calculated for each particle i at time t . The sum over the n b particles j in the localenvironment runs over all neighbors α of particle i plus the particle i itself [36]. A potential16eighbour β to α is defined as a particle j within the first coordination shell of particle α ,given by the first minimum in the radial distribution function. r αβ ( t ) is the vector betweena particle α in n b and a nearest neighbour β . The summation is further restricted to thesum over no more than the twelve nearest neighbours [25] (there are occasionally more thantwelve nearest neighbours in the first coordination shell in the supercooled liquid states).The Y lm are the spherical harmonics, and the Steinhardt order parameter is defined as Q l ( i, t ) = (cid:114) π l + 1 Σ lm = − l | q lm ( i, t ) | (2)In [36] the authors compared Q and Q for different crystal structures. The present LJsystem crystallizes into a fcc crystal, and the best separation between the bond-order in thesupercooled liquid states and in the crystal states is obtained for Q (Figure 3). Finally,a temporarily stable crystal order at particle i is determined by averaging over a shorttime-interval of one time unit Q ( i ) = < Q ( i, t ) > . (3)The clusters of particles with crystal order are obtained directly during the simulations[12]. The threshold value for crystal order is Q ( i ) ≈ .
35 accordingly to Figure 3. A crystalnucleus is determined by, that all particles in the nucleus have an order Q ( i ) > .
35, andall particles in the nucleus have at least one nearest neighbour particle j with Q ( j ) > . [1] B. J. Alder, T. E. Wainwright, J. Chem. Phys. , 1208 (1957)[2] J. P. Hansen, L. Verlet, Phys. Rev. , 151 (1969)[3] S. Auer, D. Frenkel, Nature , 1020 (2001)[4] S. Pronk, D. Frenkel, J. Chem. Phys. , 4589 (1999)[5] J. C. P´amies, A. Cacciuto, D. Frenkel, J. Chem. Phys. , 044514 (2009)[6] J. R. Errington, P. G. Debenedetti, S. Torquato, J. Chem. Phys. , 2256 (2003)[7] M. Volmer, A. Weber, Z. Phys. Chem. , 277 (1926)
8] R. Becker, W. D¨oring, Ann. Phys. (Leipzig) , 719 (1935)[9] For a resent reviewi of CNT see V. Kalikmanov, Nucleation Theory (Springer 2013)[10] A. Kuhnhold, H. Meyer, G. Amati, P. Pelagejcev, T. Schilling, Phys. Rev. E, , 052140(2019)[11] W. Ouyang, B. Sun, Z. Sun, S. Xu, J. Chem. Phys. , 054903 (2020)[12] S. Toxvaerd, J. Chem. Phys. , 154705 (2015)[13] S. Toxvaerd, J. Chem. Phys. , 164502 (2016)[14] S. Jungblut, C. Dellargo, Eur. Phys. J. E, , 77 (2016)[15] H.Tanaka, Eur. Phys. J. E, , 113 (2012)[16] P. J. Steinhardt, D. R. Nelson, M. Ronchetti, Phys. Rev. B , 784 (1983)[17] M. D. Ediger, Annu. Rev. Phys. Chem. , 99 (2000)[18] A. Widmer-Cooper, P. Harrowell, H. Fynewever, Phys. Rev. Lett. , 135701 (2004)[19] A. Widmer-Cooper, P. Harrowell, Phys. Rev. Lett. , 185701 (2006)[20] B. O’Malley, I. Snook, Phys. Rev. Lett. , 085702 (2003)[21] N. C. Karayiannis, R. Malshe, M. Kr¨oger, J. J. de Pablo, M. Laso, Soft Matter , 844 (2012)[22] P. R. ten Wolde, M. J. Ruiz-Montero, D. Frenkel, Phys. Rev. Lett. , 2714 (1995)[23] A. V. Anikeenko, N. N. Medvedev, J. Struct. Chem. , 267 (2006)[24] C. Desgranges, J. Delhommelle, Phys. Rev. Lett. , 235502 (2007)[25] J. Russo, H. Tanaka, Soft Matter , 4206 (2012)[26] J.T. Berryman, M. Anwar, S. Dorosz, T. Schilling, J. Chem. Phys. , 211901 (2016)[27] S. Dorosz, T. Schilling, J. Chem. Phys. , 124508 (2013)[28] C. Desgranges, J. Delhommelle, J. Am. Chem. Soc. , 10368 (2006)[29] C. Desgrange, J. Delhommelle, J. Phys. Chem. B , 1465 (2007)[30] S. Toxvaerd, Eur. Phys. J. Plus , 267 (2020)[31] S. Toxvaerd, J. C. Dyre, J. Chem. Phys. , 081102 (2011)[32] S. Toxvaerd, O. J. Heilmann, J. C. Dyre, J. Chem. Phys. , 224106 (2012)[33] S. Toxvaerd, Condensed Matter Physics , 13002 (2015)[34] S. Toxvaerd, J. Phys. Chem. C , 15620 (2007)[35] M. A. Barroso, A. L. Ferreira, J. Chem Phys. , 7145 (2002)[36] W. Lechner, C. Dellago, J. Chem. Phys. , 114707 (2008), 114707 (2008)