Asymmetric Assembly of Lennard-Jones Janus Dimers
Sina Safaei, Caleb Todd, Jack Yarndley, Shaun Hendy, Geoff R. Willmott
AAsymmetric Assembly of Lennard-Jones Janus Dimers
Sina Safaei, , Caleb Todd, Jack Yarndley, Shaun Hendy, , , and Geoff R. Willmott , , The MacDiarmid Institute for Advanced Materials and Nanotechnology, New Zealand Department of Physics, University of Auckland, New Zealand Te P¯unaha Matatini, Department of Physics, University of Auckland, New Zealand School of Chemical Sciences, University of Auckland, New Zealand
Self-assembly of Janus (or ‘patchy’) particles is dependent on the precise interaction between neigh-bouring particles. Here, the orientations of two amphiphilic Janus spheres within a dimer in anexplicit fluid are studied with high geometric resolution. Molecular dynamics simulations and first-principles energy calculations are used with hard- and soft-sphere Lennard-Jones potentials, andtemperature and hydrophobicity are varied. The most probable centre-centre-pole angles are in therange 40 ◦ to 55 ◦ , with pole-to-pole alignment not observed due to orientational entropy. Anglesnear 90 ◦ are energetically unfavoured due to solvent exclusion, and we unexpectedly found that therelative azimuthal angle between the spheres is affected by solvent ordering. Janus particles (JPs) are micro- or nanoparticles thathave at least two sides with distinct physical or chem-ical properties [1–5]. JPs can have asymmetry in (forexample) optical, electrical, or magnetic properties; theycan be synthesised with different shapes (e.g. spheroids,rods, platelets); and properties can be altered over par-ticular areas, resulting in ‘patchy’ particles. Amphiphilicmicrospheres, which have different wettability on eachhemisphere, are a commonly-studied type of JP [6–11].Directional interactions between JPs can lead to self-assembly of complex structures [3, 4, 12], and there hasbeen extensive interest in related design rules for Janusand patchy particles [13–16]. In 3D, experiments haveproduced aggregates ranging from small clusters up tonon-equilibrium helices with variable chirality [8]. Sim-ulations have predicted phases including small clusters[17], micellar or wire-like structures [17–21], and sheets[11, 18]. In 2D, close-packed particles [22, 23] can pro-duce tiled patterns [24–26], while less dense arrangementsinclude clusters [26], chains [27] and open lattices [28].Phases depend on parameters such as solvent properties,temperature, and patch geometry [21, 29].An analogy can be made between particle self-assemblyand formation of chemical bonds, albeit with differentgeometric rules [14, 30, 31] that are yet to be fully un-derstood. For example, the number and directions of‘bonds’ can be controlled via the positions and sizes[13, 14, 16, 26, 28, 30] of patches. Orientational en-tropy determines the flexibility (or ‘floppiness’) of bonds[8, 14, 24, 28, 30], and can determine the stability of aphase [28, 30–34].The precise nature of interactions between particlesultimately determines the possible aggregates. Simula-tions are an important predictive tool in this field [14],and a variety of anisotropic directionally-dependent pointpotentials have been used to study JP assembly [22].These potentials typically consist of an isotropic repul-sion such as an infinitely large ‘hard-sphere’ potentialbarrier [11, 18, 35, 36], as well as an anisotropic in-teraction representing the JP asymmetry. Kern and Frenkel’s approach [36], which has been widely used[11, 18, 22, 24, 25, 27, 28], applies a square-well potentialwhen attractive surfaces are in contact on the line joiningthe centres of two particles, and zero otherwise. This issuitable for modelling short-range forces, and producesan energetic step function with respect to JP orientationat the boundary of a hemisphere (or patch). As an al-ternative, the energy can be minimised when attractivehemispheres are in pole-to-pole contact, and increase ac-cording to vector dot products as the JPs are rotated[17, 19, 22, 35]. These potentials have been simplified forthe purpose of many-body calculations. There are lim-ited examples of more detailed geometric calculations ofthe energy between two JPs, in which the interaction hasbeen based on electrostatics (i.e. DLVO theory) [37, 38]and critical Casimir forces [39, 40].In this Letter, we study the 3D orientation dynamics ofamphiphilic spherical JPs in self-assembled dimers. Weuse a Lennard-Jones potential and two geometrically rig-orous calculation methods; none of these features haveappeared in previous studies of the interaction betweentwo JPs. Results obtained using molecular dynamics(MD) simulations of many-atom JPs are presented along-side calculations made using a numerical integration (NI)method we have developed. In the latter method, JPhemispheres are modelled as continuous surfaces, andpotentials are integrated to calculate the first-principlesinteraction energy for any particular configuration. Ananalytic equivalent of our NI method has been used tomodel other Lennard-Jones nanostructures [41] includinguniform spheres [42], but to our knowledge this has notbeen applied to Lennard-Jones JPs.Put together, our results produce a detailed and nu-anced description of particle interactions for short-rangedpotentials, and specifically the Lennard-Jones potentialwhich approximates van der Waals forces. The MD andNI approaches both directly include solvent, have non-zero separation of JP and solvent surfaces, and can ac-count for orientational entropy. In addition, the MDsupports study of solvent ordering, changing surface sep- a r X i v : . [ c ond - m a t . s o f t ] M a r aration lengths, and other dynamics. Solvent orderingeffects are of particular interest, as they introduce newtemperature dependent asymmetries in the dimer config-urations. Otherwise, results from MD and NI are broadlyconsistent at different temperatures and hydrophobici-ties, enabling discussion of the geometric and entropiccontributions to the preferred orientations of JPs in adimer. Aside from improving fundamental understand-ing in this way, the results can inform the developmentof design rules in self-assembled structures using ‘floppy’bonds. Dimers of JPs are particularly important becausethey are the first kinetic step in any self-assembly [8],they are possible stand-alone modular building blocks,and they appear as a motif in 2D tiling [24–26]. Molecular dynamics - MD simulations were carried outusing the large-scale atomic molecular massively parallelsimulator (LAMMPS), and the open visualisation tool(OVITO) was employed for visualisation and analysis ofthe results [43, 44]. All simulations were performed inLennard-Jones reduced units, where σ is the unit of dis-tance, (cid:15) the unit of energy, m the unit of mass and τ theunit of time ( (cid:112) mσ/(cid:15) ), without loss of generality. Theequations of motion were integrated using the standardvelocity-Verlet algorithm with a timestep of ∆ t = 0 . τ .All simulations were carried out in a canonical ensem-ble (NVT) in which the temperature was set at 1.1 (cid:15)/k B and controlled by a Langevin thermostat with a damp-ing parameter of 10 τ (for methodological details, see theSupplemental Material [45]).MD studies of JPs can typically be categorised intothose which model Janus particles as single ‘atoms’ us-ing point potentials [17, 20, 22, 46] or those which rep-resent spherical JPs as collections of atoms distributedover the sphere’s surface (Fig. 1(a)) [47–51]. A many-atom JP is important for the study of hydrodynamic ef-fects such as the effect of slip on particle rotation andtranslation [51–54]. Here we use the latter approach,with Janus spheres modelled as in our previous studies[50, 51], where the atoms were distributed over the sur-face of the spheres with a surface density of 1 . σ − . Ineach simulation, the positions of atoms on the JP surfacewere varied based on a randomised process [51], avoidingartefacts associated with regular structures (see the Sup-plemental Material [45]). Anisotropy of JPs emerges bydefining different potentials for atoms on different partsof the surface. Hence, the spheres were divided geomet-rically into two hemispheres referred to as and , andthey interacted with atoms of an explicit solvent referredto as . To design the spheres as similarly as possibleto experimentally fabricated JPs, the boundary betweenthe hemispheres is not sharp. This eliminates any re-lated effects [51] and means that the number of atoms ondifferent hemispheres may not be equal.The potential between atoms separated by distance r was a 12-6 Lennard-Jones potential, FIG. 1. (a) A dimer of Janus spheres (blue: hydrophobic,red: hydrophilic) showing the polar angles θ i and θ j and theazimuthal angle φ i . The angle φ j = 90 ◦ is omitted for clarity. C i and C j are the centres of the spheres, and the vectorsconnect the centres to hydrophobic poles. Insets illustratesphere orientations at specified points on plot (b). (b) Theoccurrence probability for θ i and θ j in MD simulations at T = 1 . (cid:15)/k B , A = 0 . (cid:15) , and R = 10 σ . U ij ( r ) = 4 (cid:15) (cid:20)(cid:16) σr (cid:17) − (cid:16) σr (cid:17) (cid:21) , (1)with the interaction strengths between different types ofatoms defined as: (cid:15) (1 − = A o (cid:15) (0 − (cid:15) (2 − = A i (cid:15) (0 − . ≤ A o ≤ A i ≤ . . Here A (the ‘interaction strength’) is a parameter allow-ing the potential to be adjusted for the hydroph o bic ( A o )and hydroph i lic ( A i ) hemispheres in a Janus sphere,which are depicted in all figures as blue and red, re-spectively. As the interaction strength increases, the hy-drophobicity of the surface decreases. The interactions (cid:15) (1 − , (cid:15) (2 − , and (cid:15) (1 − were calculated using the Berth-elot mixing rule. Parameters m and σ were the same forall atom types.To study the orientation of individual JPs in a dimer,two spherical coordinate systems ( r i,j , θ i,j , φ i,j ) wereemployed (Fig. 1(a)), where the origins are at the cen-tres of each of the spheres C i,j , and the line connectingthe centres ( −−−→ C i C j ) indicates the direction of the z -axis.The orientation of each sphere is defined by a vector con-necting its centre to its hydrophobic pole. The axes aredefined so that φ j = 90 ◦ and the relative orientation ofthe two Janus spheres can be fully specified using θ i , θ j ,and ∆ φ (where ∆ φ = φ j − φ i ).Fig. 1(b) gives an example of a simulated occurrenceprobability distribution as a function of the polar angles θ i and θ j . Here, the observed angles are most commonlyin the range from 30 ◦ to 60 ◦ , and orientations with an-gles higher than 75 ◦ are not observed. Orientations withangles less than 30 ◦ are not favoured due to the loweravailable surface area at those angles - this is the effect oforientational entropy. The sudden decrease above 60 ◦ isrelated to screening of solvent atoms from the hydropho-bic hemispheres, as discussed below. Numerical integration - In the numerical integration(NI) approach, we calculated the interaction between twoJanus spheres over a range of configurations. Each sphereconsisted of two continuous hemispherical surfaces sepa-rated by a uniform, equatorial boundary. The energy atany given configuration was calculated by integrating theLennard-Jones interaction over the different geometries,resulting in a total of 8 integrals (4 surface-surface and4 surface-solvent). In addition to the hard-sphere (12-6)potential in Eq. 1, we carried out calculations using asoft-sphere (9-6) Lennard-Jones potential. This was ofinterest due to frequent use of soft directional-dependentpotentials in previous studies of Janus spheres [19, 55–61].The centres of the two spheres were positioned at(0 , ,
0) and (0 , , R + d s ) in Cartesian coordinates(Fig. 1(a)), where R is the sphere radius and d s is thedistance between the spheres’ surfaces. A fixed distancebetween the solid surfaces and the solvent ( d f ) was alsodefined. In MD, d s and d f can continuously adjust inresponse to the configuration, and d f is (on average) dif- FIG. 2. The occurrence probability of orientations with differ-ent θ angles for Janus spheres in a dimer. The inset shows theenergy landscape of the spheres configurations as a functionof the polar angles. Lines are drawn to guide the eye. ferent for the hydrophilic and hydrophobic sides. Defaultvalues of d s = 0 . σ and d f = 0 . σ were chosen by ex-tracting averages from MD, but calculations were alsomade using several different values to show that theseparameters have a minimal effect on orientation proba-bility distributions (see Supplemental Material). Other-wise, NI calculation parameters (e.g. surface and fluiddensities) matched the values used in MD simulations.Given a configuration energy E ( θ j , θ i , φ i ), the probabil-ity density can be determined using p ∝ sin( θ j ) sin( θ i ) exp (cid:20) − Ek B T (cid:21) , (2)where k B is Boltzmann’s constant and the number ofways a configuration of θ i , θ j can be obtained is propor-tional to sin( θ j ) sin( θ i ). The Supplemental Material con-tains further detail regarding the NI method [45]. Potential effect - Fig. 2 compares the occurrence prob-ability of polar angles calculated using MD and NI meth-ods. The curves are broadly similar, and in particular themaximum occurrence probability occurs at very similarangles for the MD simulations and the NI calculationwith a matching hard-sphere potential. The likely rea-sons for discrepancies between these plots are the knownminor differences between MD and NI models, i.e. thefixed values of d s and d f and the sharp boundary be-tween hemispheres in the NI model. For the soft-spherepotential, the probability distribution shifts slightly to-ward smaller angles.To further understand these probability distributions,consider the calculated energy for the interaction betweenthe two spheres (inset to Fig. 2) which is approximatelyconstant for all configurations with θ i,j < ◦ . In this re-gion, the reason for variation in occurrence probability isthe lower configurational entropy (due to smaller surfacearea) near the poles. This is characteristic of potentialswhich have little or no dependence on orientation, rep-resenting short-range interactions [25]. In contrast, ap-proaches which use geometrically rigorous long-range in-teractions (electrostatic or DLVO) produce continuouslyvarying potential landscapes [37, 38].At θ > ◦ , the interaction energy gradually increases,producing a sharp decline in the probability distributionfor 60 ◦ < θ i,j < ◦ . This decline is not present whenusing a Kern-Frenkel type potential. In those cases, themaximum probability is found adjacent to the step func-tion (at θ i,j ≈ ◦ for hemispheres), where the configu-rational entropy is greatest. Labb´e-Laurent and Dietrich[40] derived a geometrically rigorous potential for a Janusdimer interacting via Casimir forces, which produced asimilar combination of constant and varying potentialsas in Fig. 2 inset, albeit without solvent.The main reason for the increase in energy over therange 60 ◦ < θ i,j < ◦ is that there is no solvent in theregion near the point of closest approach between thespheres. It is energetically favourable for the hydrophobicsurfaces to be shielded from the solvent near this contactpoint. The volume in which there is no fluid (an ‘exclu-sion zone’) is determined by d f . The polar angle at whicha sphere’s hydrophilic side enters this exclusion zone (andthus the hydrophobic side is not optimally hidden fromthe fluid) is θ i,j = 90 ◦ − arccos( R + d s / R + d f ) (see Supplemen-tal Material [45]), or θ i,j = 72 ◦ for our parameters, whichmatches well with Figs. 1(b) and 2. The precise shape ofthese distributions additionally depends on the form andrange of the interaction potential.The shift of the probability distribution to lower polarangles when using a soft-sphere (9-6) potential can beexplained by the increased importance of long range in-teractions. For example, the minimum of the soft poten-tial occurs at a larger radius than for the 12-6 potential.Increasing the range of the potential decreases the polarangle required for the hydrophilic side to be affected bythe exclusion zone, thereby producing a decline in prob-ability at lower polar angles.Temperature was varied between 1 . (cid:15)/k B and 1 . (cid:15)/k B to study the effect on JP orientation. Fig. 3(a) shows po-lar angle probability distributions for two different tem-peratures calculated using both MD (via repeated simu-lations) and NI (using Eq. 2). At higher temperatures,the spheres spend more time in high energy orientations( θ i,j > ◦ ) even up to unfavourable angles between 75 ◦ and 80 ◦ . Figure 3(b) indicates the expected value of θ i,j at different temperatures. The results and trends fromthe MD and NI methods match closely, with a small dif-ference possibly caused by fluid layering in MD (discussed FIG. 3. (a) Occurrence probability as a function of polarangle for Janus spheres in a dimer at different temperatures( R = 10 σ ). Lines are drawn to guide the eye. (b) and (c),expected value of the polar angle as a function of temperatureand hydrophobicity, respectively. further below).To study the effect of hydrophobicity on orientation,the hydrophobic interaction strength was varied from A o = 0 .
05 to 0 . (cid:15) at the default temperature (1 . (cid:15)/k B ).This produces very similar changes in the overall proba-bility distribution to variations in temperature (Fig. 3(a);see Supplemental Material [45]) and similar trends in theexpected value of θ i,j (Fig. 3(c)) albeit for a differentreason. As the hydrophobicity decreases, the energy ofthe fluid-hydrophobic interactions becomes more nega-tive which allows the hydrophobic sides of the spheresto be more exposed to the solvent atoms. Nonetheless,the spheres avoid forming orientations with polar angleshigher than 70 ◦ .While the NI calculations showed little dependence onthe azimuthal angle (see the Supplemental Material [45]),an unexpected variation in probability was observed forMD results. As shown in Fig. 4, smaller values of ∆ φ arefavoured. Data for two different temperatures (Figs. 4(a)and (b)) suggest that this variation is greater at highertemperature. This result can be explained by layeringof solvent atoms around the spheres, which shows lessordering at smaller ∆ φ angles, increasing the entropy of FIG. 4. The occurrence probability for dimer orientations in MD simulations ( A = 0 . (cid:15) , R = 10 σ ) at (a) T = 1 . (cid:15)/k B , and(b) T = 1 . (cid:15)/k B . (c) The spatial distribution of solvent atoms near the hydrophilic sides of JPs at ∆ φ = 5 ◦ and ∆ φ = 175 ◦ .Insets, which represent the orientations from (a) at the corresponding labels, illustrate regions of layered solvent atoms used inthe histogram. (ii) does not cover a full semi-circle to maintain equivalence with (i), which has an overlapping region. the system. The free energy for orientations with smaller∆ φ is therefore lower, and decreases with increasing tem-perature.The layering can be observed in Fig. 4(c), which showsthe spatial distribution of solvent atoms close to hy-drophilic surfaces for ∆ φ = 5 ◦ and ∆ φ = 175 ◦ when55 ◦ < θ i , θ j < ◦ . The first three layers of solventatoms are approximately the same for both cases, whilethe fourth and fifth peaks are somewhat stronger at∆ φ = 175 ◦ . The insets to Fig. 4(c) illustrate the rea-son for this difference. At low ∆ φ , the hydrophilic sidesare close to each other so that regions of adjacent lay-ered solvent overlap. Layers are observed at short range( (cid:46) σ ) but the particles each influence longer range or-dering of solvent adjacent to their neighbour. At higher∆ φ , the hydrophilic faces are further away from eachother and the layering of solvent atoms is not disturbedby the neighbouring sphere. Near the hydrophobic faces,the range of ordered solvent atoms is much shorter andthere is negligible contribution to ordering.To conclude, our calculations using two independentmethods have yielded consistent results for the favouredorientations of a dimer of amphiphilic JPs interactingvia Lennard-Jones potentials in a solvent. The mostfavoured configurations are at polar angles between 40and 55 ◦ . Pole-to-pole configurations are unfavoured dueto orientational entropy and the short-range nature ofthe potential, while orientations near 90 ◦ are energeti-cally unfavourable due to solvent exclusion between theparticles. MD revealed that dynamic changes in the sep-aration between particles and solvent atoms have a smalleffect on orientation probabilities, as do the details of theboundary between hemispheres. MD also showed thatsolvent ordering can break the degeneracy with respectto ∆ φ . The NI method, which directly calculates theinteraction energy, could be extended to ranges of prac-tically important parameters, such as different potentials or non-hemispherical patches. Considering the role of ori-entational entropy, these methods should be particularlyuseful for designing patches which yield desired ‘floppybonds’ geometries. The calculation efficiency of thesemethods suggest that they could be extended to studyclusters, and then used to find simplified yet accuratepotentials for simulating many-particle structures.This research was funded by The MacDiarmid Insti-tute for Advanced Materials and Nanotechnology, andsimulations were run on the NeSI Mahuika and M¯auiClusters, part of the Centre for eResearch hosted by theUniversity of Auckland. [1] L. M. Bergstr¨om, Application of Thermodynamics to Bi-ological and Material Science , 289 (2011).[2] Z. Yang, A. H. Muller, C. Xu, P. S. Doyle, J. M. DeSi-mone, J. Lahann, F. Sciortino, S. Glotzer, L. Hong, D. A.Aarts, et al. , Janus particle synthesis, self-assembly andapplications (Royal Society of Chemistry, 2012).[3] E. Poggi and J.-F. Gohy, Colloid and Polymer Science , 2083 (2017).[4] M. N. Popescu, Langmuir , 6861 (2020).[5] A. Perro, S. Reculusa, S. Ravaine, E. Bourgeat-Lami,and E. Duguet, Journal of Materials Chemistry , 3745(2005).[6] A. Walther and A. H. M¨uller, Soft Matter , 663 (2008).[7] A. Walther and A. H. Muller, Chemical Reviews ,5194 (2013).[8] Q. Chen, J. K. Whitmer, S. Jiang, S. C. Bae, E. Luijten,and S. Granick, Science , 199 (2011).[9] K.-H. Roh, D. C. Martin, and J. Lahann, Nature Mate-rials , 759 (2005).[10] S. Jiang, Q. Chen, M. Tripathy, E. Luijten, K. S.Schweizer, and S. Granick, Advanced Materials , 1060(2010).[11] Z. Preisler, T. Vissers, G. Munao, F. Smallenburg, andF. Sciortino, Soft Matter , 5121 (2014). [12] N. Safaie and R. C. Ferrier Jr, Journal of Applied Physics , 170902 (2020).[13] J. Zhang, B. A. Grzybowski, and S. Granick, Langmuir , 6964 (2017).[14] J. Zhang, E. Luijten, and S. Granick, Annual Review ofPhysical Chemistry , 581 (2015).[15] M. F. Hagan, O. M. Elrad, and R. L. Jack, The Journalof Chemical Physics , 104115 (2011).[16] B. A. Grzybowski, K. Fitzner, J. Paczesny, andS. Granick, Chemical Society Reviews , 5647 (2017).[17] G. Rosenthal, K. E. Gubbins, and S. H. Klapp, TheJournal of Chemical Physics , 174901 (2012).[18] G. Munao, Z. Preisler, T. Vissers, F. Smallenburg, andF. Sciortino, Soft Matter , 2652 (2013).[19] Z.-W. Li, Z.-Y. Lu, Z.-Y. Sun, and L.-J. An, Soft Matter , 6693 (2012).[20] F.-F. Hu, Y.-W. Sun, Y.-L. Zhu, Y.-N. Huang, Z.-W. Li,and Z.-Y. Sun, Nanoscale , 17350 (2019).[21] M. M. Moghani and B. Khomami, Soft Matter , 4815(2013).[22] (cid:32)L. Baran, M. Bor´owko, and W. R˙zysko, The Journal ofPhysical Chemistry C (2020).[23] I. Kretzschmar and J. H. K. Song, Current Opinion inColloid & Interface Science , 84 (2011).[24] Y. Iwashita and Y. Kimura, Soft Matter , 7170 (2014).[25] H. Shin and K. S. Schweizer, Soft Matter , 262 (2014).[26] Q. Chen, J. Yan, J. Zhang, S. C. Bae, and S. Granick,Langmuir , 13555 (2012).[27] Y. Iwashita and Y. Kimura, Soft Matter , 10694 (2013).[28] X. Mao, Q. Chen, and S. Granick, Nature Materials ,217 (2013).[29] W. L. Miller and A. Cacciuto, Physical Review E ,021404 (2009).[30] F. Sciortino, La Rivista del Nuovo Cimento , 511(2019).[31] Y. Wang, Y. Wang, D. R. Breed, V. N. Manoharan,L. Feng, A. D. Hollingsworth, M. Weck, and D. J. Pine,Nature , 51 (2012).[32] Q. Chen, S. C. Bae, and S. Granick, Nature , 381(2011).[33] F. Sciortino, A. Giacometti, and G. Pastore, PhysicalReview Letters , 237801 (2009).[34] F. Smallenburg and F. Sciortino, Nature Physics , 554(2013).[35] L. Hong, A. Cacciuto, E. Luijten, and S. Granick, Lang-muir , 621 (2008).[36] N. Kern and D. Frenkel, The Journal of Chemical Physics , 9882 (2003).[37] L. Hong, A. Cacciuto, E. Luijten, and S. Granick, NanoLetters , 2510 (2006). [38] R. Hieronimus, S. Raschke, and A. Heuer, The Journalof Chemical Physics , 064303 (2016).[39] A. Squarcini, A. Macio(cid:32)lek, E. Eisenriegler, and S. Diet-rich, Journal of Statistical Mechanics: Theory and Ex-periment , 043208 (2020).[40] M. Labb´e-Laurent and S. Dietrich, Soft Matter , 6621(2016).[41] D. Baowan and J. M. Hill, Advances in Mechanical En-gineering , 1687814016677022 (2016).[42] D. Baowan, B. J. Cox, T. A. Hilder, J. M. Hill, andN. Thamwattana, Modelling and mechanics of carbon-based nanostructured materials (William Andrew, 2017).[43] S. Plimpton, Journal of Computational Physics , 1(1995).[44] A. Stukowski, Modelling and Simulation in Materials Sci-ence and Engineering , 015012 (2009).[45] See Supplemental Material at [URL will be inserted bypublisher] for further details on MD and NI methods. [46] T. Vissers, Z. Preisler, F. Smallenburg, M. Dijkstra,and F. Sciortino, The Journal of Chemical Physics ,164505 (2013).[47] Y.-C. Li, N.-B. Zhang, Z. Wei, B.-Y. Li, M.-T. Li, andY. Li, Molecular Simulation , 759 (2019).[48] J. Xu, Y. Wang, and X. He, Soft Matter , 7433 (2015).[49] Y. Kobayashi, N. Arai, and A. Nikoubashman, Soft Mat-ter , 476 (2020).[50] S. Safaei, S. C. Hendy, and G. R. Willmott, Soft Matter , 7116 (2020).[51] S. Safaei, A. Y. Archereau, S. C. Hendy, and G. R.Willmott, Soft Matter , 6742 (2019).[52] G. Willmott, Physical Review E , 055302 (2008).[53] G. R. Willmott, Physical Review E , 066309 (2009).[54] A. Kharazmi and N. V. Priezjev, The Journal of Chemi-cal Physics , 234503 (2015).[55] Z.-W. Li, Z.-Y. Lu, Y.-L. Zhu, Z.-Y. Sun, and L.-J. An,RSC Advances , 813 (2013).[56] Z.-W. Li, Z.-Y. Lu, and Z.-Y. Sun, Soft Matter , 5472(2014).[57] Z.-W. Li, Y.-L. Zhu, Z.-Y. Lu, and Z.-Y. Sun, Soft Mat-ter , 741 (2016).[58] Z.-W. Li, Y.-L. Zhu, Z.-Y. Lu, and Z.-Y. Sun, PhysicalChemistry Chemical Physics , 32534 (2016).[59] Z.-W. Li, Y.-L. Zhu, Z.-Y. Lu, and Z.-Y. Sun, Soft Mat-ter , 7625 (2018).[60] Q.-Z. Zou, Z.-W. Li, Z.-Y. Lu, and Z.-Y. Sun, Nanoscale , 4070 (2016).[61] Q.-Z. Zou, Z.-W. Li, Y.-L. Zhu, and Z.-Y. Sun, SoftMatter15