From Sticky-Hard-Sphere to Lennard-Jones-Type Clusters
Lukas Trombach, Robert S. Hoy, David J. Wales, Peter Schwerdtfeger
FFrom Sticky-Hard-Sphere to Lennard-Jones-Type Clusters
Lukas Trombach, Robert S. Hoy, David J. Wales, and Peter Schwerdtfeger
1, 4, ∗ Centre for Theoretical Chemistry and Physics, The New Zealand Institute for Advanced Study,Massey University Auckland, Private Bag 102904, 0632 Auckland, New Zealand Department of Physics, University of South Florida, Tampa, Florida 33620, USA University Chemical Laboratories, Lensfield Road, Cambridge CB2 1EW, UK Centre for Advanced Study (CAS) at the Norwegian Academy of Science and Letters, Drammensveien 78, NO-0271 Oslo, Norway (Dated: February 2, 2018)A relation M SHS → LJ between the set of non-isomorphic sticky hard sphere clusters M SHS and the sets oflocal energy minima M LJ of the ( m , n )-Lennard-Jones potential V LJ mn ( r ) = ε n − m [ mr − n − nr − m ] is established.The number of nonisomorphic stable clusters depends strongly and nontrivially on both m and n , and increasesexponentially with increasing cluster size N for N (cid:38)
10. While the map from M SHS → M
SHS → LJ is non-injectiveand non-surjective, the number of Lennard-Jones structures missing from the map is relatively small for clustersizes up to N =
13, and most of the missing structures correspond to energetically unfavourable minima evenfor fairly low ( m , n ). Furthermore, even the softest Lennard-Jones potential predicts that the coordination of13 spheres around a central sphere is problematic (the Gregory-Newton problem). A more realistic extendedLennard-Jones potential chosen from coupled-cluster calculations for a rare gas dimer leads to a substantialincrease in the number of nonisomorphic clusters, even though the potential curve is very similar to a (6,12)-Lennard-Jones potential. I. INTRODUCTION
The nucleation of atoms and molecules in the gas phase, orliquid, to the solid state is still an active research field [1–8].Rowland noted in 1949 that “The gap between theory and theexperimental approaches to nucleation has been too wide” and“the subject [nucleation] is still in the alchemical stage” [9].More than half a century later, despite all the advancementsmade in cluster physics, “there is still a large gap betweenexperiment and theory” as Unwin noted [10].The underlying reason for this rather slow progress is thatcluster formation is a dynamic process, and fully charac-terizing the corresponding high-dimensional potential energylandscape is typically an NP-hard problem, since there are(presumably) exponentially many local minima at any giventemperature and pressure [1, 11–16]. Moreover, phase transi-tions between di ff erent morphologies as a function of size N usually occur where N is too large for an accurate quantum-theoretical treatment [17–21]. For example, Krainyukova ex-perimentally studied the growth of argon clusters [22], andfound that small, initially icosahedral clusters transform intoanti-Mackay clusters for N > N > atoms, in qualita-tive agreement with theoretical predictions using Lennard-Jones (LJ) type potentials [2, 23, 24]. The notorious rare gasproblem was solved only very recently by accurate relativisticquantum methods, correctly predicting a slight preference ofthe fcc over the hcp phase due to phonon dispersion [25].Simple models often have to be used to simulate clustergrowth and nucleation [26–29]. The simplest model potentialsthat can be applied to theoretical studies of atomic cluster for-mation are “HCR-SRA” potentials with isotropic hard-core-like repulsive and short-range-attractive interactions [30]. The ∗ Email: [email protected] simplest HCR-SRA potential is the “sticky hard sphere”(SHS) potential [31] V SHS ( r ) = ∞ , r < r s , − ε, r = r s , , r > r s , (1)where r s and ε can be arbitrarily set to 1 (unit sphere andreduced units, respectively). Eq. 1 can be used as a pertur-bative basis for finite-ranged HCR-SRA potentials [32, 33].Since sticky hard spheres are impenetrable and their energy E = − N c ε is a function only of the number of interpar-ticle contacts N c , SHS cluster structure and energetics canbe uniquely mapped to their adjacency matrices ¯ A , where N c = (cid:80) Ni < j A i j . This mapping allows them to be exactly charac-terized via complete enumeration [34–36]; recent studies haveidentified all mechanically stable SHS clusters for N ≤
14, andputatively complete sets for N ≤
19 [34–39]. Note, however,that di ff erent SHS structures can have the same adjacency ma-trix for N ≥
14 [38], and the mapping is therefore only surjec-tive.From the Gregory-Newton kissing-number argumentproved in 1953 by Sch¨utte and van der Waerden [40], nosphere can be surrounded by more than 12 spheres of equalradius [41]. For small clusters, graph-theoretic arguments dic-tate max( N c ) ≤ N ( N − /
2. Thus a loose bound on the maxi-mum contact number N c ( N ) is N max c ( N ) ≤ min { N ( N − / , f ( N ) } (2)with f ( N ) = N . This upper bound has been tightened severaltimes, most recently by Bezdek and Reid [42] to f ( N ) = N − / π − / N / . (3)In Refs. [37, 38] it was shown that N max c ( N ) = { , , , , , , , , , , , , , , , } for4 ≤ N ≤
19. While determining N max c ( N ) for arbitrary N is a r X i v : . [ phy s i c s . a t m - c l u s ] F e b -1-0.5 0 0.5 1 0 0.5 1 1.5 2 E / ε r/r e n =2 n =6 n =12 n =24 n =48 n =96 FIG. 1: Lennard-Jones potentials for di ff erent exponents( m , n ) with fixed n = m . As the exponents grow larger thewell of attraction becomes narrower and its shape approachesthe SHS potential. The dashed line shows the extended LJpotential for the xenon dimer [48].equivalent to the still-unsolved Erd¨os unit distance problem[43], it is clear that N max c ( N ) = N − + m ( N ), where m ( N )grows slowly from zero to around f ( N ) − (3 N −
6) withincreasing N .While the maximum contact number increases (sub)linearlywith N , the number of non-isomorphic cluster structures |M ( N ) | and transition states is assumed to increase exponen-tially [11, 44, 45] (here we denote M ( N ) as the set of allnon-isomorphic cluster structures of size N , and |M ( N ) | asthe number of structures in M ( N )). Stillinger showed that un-der certain conditions lim N →∞ |M ( N ) | ∝ exp( α N ) [44]. ForSHS clusters, the complete set M SHS ( N , N c ) has been exactlydetermined for N ≤
14 and 3 N − ≤ N c ≤ N max c ( N ) via ex-act enumeration studies employing geometric rejection rules[37, 38]. Unfortunately, such precise calculations are verydi ffi cult for finite-ranged potentials since exhaustive searchesfor energy minima are computationally intensive [46]. Only afew such studies have been performed, e.g. recent studies of N ≤
19 clusters interacting via short-range Morse potentials[13, 15, 47].It remains unclear how the HCR-SRA models commonlyused in cluster physics relate to more physically relevant,softer interaction potentials such as the ( m , n )-Lennard-Jones(LJ) form: V LJ m , n ( r ) = ε n − m (cid:20) m (cid:18) r e r (cid:19) n − n (cid:18) r e r (cid:19) m (cid:21) (with n > m ) . (4)Here ε > r e the equilibriumtwo-body interparticle distance. To simplify the presentation,we (without loss of generality) adopt reduced units ( ε = r e =
1) below. For m , n → ∞ , V LJ m , n ( r ) → V SHS ( r ) (Fig. 1);the energy landscapes of the two potentials converge in thislimit. However, real systems are not in this limit. For exam-ple, for N =
13, there are |M SHS | = ,
221 stable SHS clusters[37, 38], but only |M LJ | = ,
510 stable ( m , n ) = (6 ,
12) LJ clus-ters [49]. This di ff erence is understood qualitatively – energy landscapes are well known to support more local minima asthe range of the interaction potential decreases [50, 51]. Thereare several e ff ects that will cause the set of stable LJ clus-ters to increasingly deviate from the set of stable SHS clustersas interactions become longer ranged. As n and m decrease,second-nearest-neighbor attractions become increasingly im-portant, producing stable structures with r i j ≤
1. Fold catastro-phes [51, 52] progressively eliminate stable SHS clusters, andseveral stable SHS structures may collapse into a single sta-ble LJ cluster. However, detailed quantitative understandingof such e ff ects remains rather limited.In this paper, we quantitatively examine how stable N ≤ m , n )decrease. We focus on both the topography of the energy land-scape (decreasing |M LJ ( N ) | ) and the evolving topologies ofthe stable cluster sets. We examine these changes in furtherdetail for specific N = −
14 clusters discussed by Gregoryand Newton in the 1600s in the context of the kissing numberproblem [40], and also for a more realistic two-body poten-tial that has been shown to accurately model rare-gas clusters[23].
II. COMPUTATIONAL METHODS
The pele program [53] was used to generate putatively com-plete sets of local minima for ( m , n )-Lennard-Jones potentials V LJ mn ( r ) as defined in Eq.(4). This program applies a basin-hopping algorithm that divides the potential energy surfaceinto basins of attraction, e ff ectively mapping each point inconfiguration space to a local minimum structure [54–56].The results confirmed the number of local minima reported inprevious work [57]. Finite computer time limited our searchto clusters of size N ≤ N = N c < N −
6) [38], we carried out geometry optimisations with ( m , n )-Lennard-Jones potentials using the multidimensional func-tion minimiser from the C ++ library dlib [58]. The opti-misation scheme was either the Broyden-Fletcher-Goldfarb-Shanno (BFGS) or the conjugate gradient (CG) algorithm.The optimisations were terminated when the change in en-ergy (in reduced units) over the course of one optimizationcycle was smaller than 10 − . Subsequently, the eigenvaluesof the Hessian were checked for all stationary points. If nega-tive eigenvalues were found, the a ff ected structures were reop-timized following displacements in both directions along thecorresponding eigenvectors to locate true local minima. Thisprocedure assures that the floppy SHS packings are success-fully mapped into LJ minima.As the optimisations often result in many duplicates, espe-cially for small values of n and m where we have |M ( m , n ) − LJ | (cid:28)|M SHS | , the final structures were further analysed and sorted.Nonisomorphic SHS clusters can be distinguished (apart frompermutation of the particles) by their di ff erent adjacency ma-trices for N ≤
13 [38]. This is not the case for soft potentialslike the LJ potential since drawing edges (bonds) between thevertices (atoms) becomes a matter of defining the distance cut-o ff criterion for a bond to be drawn. Therefore, we comparethe interparticle distances { r i j } instead: two clusters are iso-morphic (structurally identical) if they have the same orderedset of inter-particle distances { r i j } . While enantiomers cannotbe separated using this methodology, permutation-inversionisomers are usually lumped together, since the number of dis-tinct minima is analytically related to the order of the corre-sponding point group [52]. To verify the number of distinctstructures we introduced a second ordering scheme using theenergy and moment of inertia tensor eigenvalues.Two sets of structures are obtained from our optimiza-tion procedure: the first set contains all possible LJ minima M LJ from the basin-hopping algorithm, while the second set M SHS → LJ contains the LJ minima obtained using only the M SHS sticky-hard-sphere cluster structures as starting pointsfor the geometry optimization. To compare and identify cor-responding structures between the two sets, the N ( N − / { r i j } were again used as an identifyingfingerprint.Two-body “extended Lennard-Jones” (ELJ) potentials thataccurately model two-body interactions in rare-gas clusterscan be written as expansions of inverse-power-law terms [23]: V ELJ ( r ) = (cid:88) n c n r − n , (5)where in reduced units the condition (cid:80) n c n = − ffi cients for the ELJ potential (in reducedunits): c = − . c = − . c = − . c = + . c = − . c = + . c = − . c = + . c = − . c = + . c = − . c = + . III. RESULTSA. Exploring the limits of Lennard-Jones
To study the convergence behavior of the number of distinct(nonisomorphic) LJ minima in the SHS limit, we performedgeometry optimisations, starting from all nonisomorphic SHSstructures. We will show later that the number of unique min-ima obtained in this procedure |M SHS → LJ | only misses out ona small portion of minima obtained from the more exhaustivebasin-hopping approach, i.e. |M SHS → LJ | ≈ |M LJ | . The resultsfor a constant chosen ratio of LJ exponents n / m = |M SHS → LJ | smoothly converges towards the SHS limit(dashed line, values in Table I) from below, thusdemonstrating that for LJ systems the number of dis-tinct minima does not grow faster than exponentially.The (48,96)-LJ potential has ∆ M ≡ |M LJ | − |M SHS → LJ | =
10 20 30 40 50 60 70 80 90 | M S H S → L J | nN =8 N =9 N =10 N =11 N =12 N =13 N =14
20 30 40 50 60 70 80 90 | M S H S → L J | n FIG. 2: Convergence of the number of distinct LJ localminima |M SHS → LJ | obtained through geometry optimisationsstarting from the nonisomorphic SHS structures withincreasing LJ exponent n . Permutation-inversion isomers andenantiomers are not distinguished. The dashed line gives theexact SHS limit |M SHS | . Top panel: m = n /
2. Bottom panel:fixed m = { , , , , , , } fewer stable minima than theSHS potential. The fractions of missing minima ∆ M / |M SHS | for this potential grow with increasing N and are respec-tively { . , . , . , . , . , . , . } %. Note thatfor N ≥
10 most of these missing minima correspond to highenergy ( N c < N max c ) structures.If the exponent n for the repulsive part of the LJ potentialis increased with m kept constant, the LJ potential becomesequivalent to the SHS potential in the repulsive range but re-mains attractive at long range. Figure 2 (bottom) shows theconvergence of the number of unique structures with respectto n at set m = α SHS ≈ α (6,12)-LJ ≈ l n | M | N SHS(6,12) LJ
FIG. 3: Growth behaviour of |M ( N ) | of SHS and (6,12)-LJclusters and corresponding asymptotic exponential rise rateparameter α for N ≥
12 as defined in Eq.(6). The interceptsln |M ( N = | are − .
19 and − .
94 for the SHS and(6,12)-LJ cases respectively.minima compared to the rigid SHS model.To see if the asymptotic increase in the number of distinctminima |M ( N ) | ∼ e α N is indeed exponential, we use Still-inger’s expression for the asymptotic exponential rise rate pa-rameter [44] α = lim N →∞ (cid:16) N − ln |M ( N ) | (cid:17) . (6)Figure 3 shows the number of distinct minima for SHS clus-ters obtained from the data shown in Table I. The N ≥
12 SHSdata gives α SHS ≈ .
21. Figure 3 also shows the (6,12)-LJresults obtained using basin-hopping; these yield α LJ ≈ . α = . |M SHS | / |M LJ | with N isexplained by the much larger values of α for the SHS com-pared to the LJ clusters.Using the results for N ≥
13 from Figure 2, we can calculatehow α depends on the LJ range parameter n . As shown inFigure 4, a general function of the form α ( n ) = α max + a ( n − n ) p (7)fits the results nicely, allowing the prediction of growth be-haviour for di ff erent LJ potentials. For |M ( n / , n ) − LJ | , α max isequivalent to α SHS = . a = − . n = − .
386 and p = .
473 (Figure 4). Wealso show the ratio α ( |M SHS → ( n / , n ) − LJ | ) /α ( |M SHS → (6 , n ) − LJ | )between the two di ff erent LJ asymptotic exponential rise rateparameters, which shows that larger cluster sizes need to bestudied to correctly describe the asymptotic limit.The distribution of minima as a function of (free) energywas suggested to be Gaussian [60]. Figure 5 shows the en-ergy distribution of minima for di ff erent LJ ( n / , n ) potentialsderived from SHS initial structures. We do not see a Gaus-sian type of distribution; this result does not change if we α n M SHS → ( n /2, n )-LJ M SHS → (6, n )-LJ r a t i o FIG. 4: Convergence behaviour of the asymptoticexponential rise rate parameter α (Eq.(6)) towards the SHSlimit with respect to the LJ exponent n . The inlet shows theratio of the two quantities α ( |M SHS → ( n / , n ) − LJ ( N ) | ) /α ( |M SHS → (6 , n ) − LJ ( N ) | ). r e l a t i v e f r e q u e n c y E / ε n =12 n =18 n =24 n =84SHS 0.0001 0.001 0.01 0.1 1-40 -38 -36 -34 FIG. 5: Histogram of the energies (bin size ∆ ε = .
1) ofminima M SHS → ( n / , n ) − LJ ( N ) for N =
13 and di ff erentexponents n up to the SHS limit. For better visibility, theheight of the bars are set to ∆ |M| / |M| in the interval ∆ ( E /(cid:15) ).The inlet shows the same data in logarithmic scale.take the free energy at finite temperatures. The results in-dicate a “phase transition” in the potential energy landscapeaway from low-energy to high energy minima as n increases.The transition occurs at fairly small n . Results for the (9 , , N c =
34 and N c =
35 SHS clusters, respectively. It is alsoclear that (as expected) the distributions narrow with increas-ing n .It is well known that the global minimum for rare gas clus-ters with 13 atoms is the ideal Mackay icosahedron [61–63].Simple geometric considerations imply that such a symmetricTABLE I: Number of distinct local minima |M SHS | forcluster size N (from Refs. [36–38]) and contact number N c from the exact enumeration, compared to the number ofdi ff erent structures obtained from a geometry optimisationstarting from the set M SHS → LJ ( N , N c ) for a (6,12)-LJpotential. The overall number of unique minima |M SHS → LJ | = (cid:80) N c |M SHS → LJ ( N c ) | − ( |M LJ | ). The di ff erence ∆ M = |M LJ | − |M SHS → LJ | is also listed. N N c |M SHS ( N c ) | |M SHS → LJ ( N c ) | |M SHS → LJ | |M LJ | ∆ M a a a b (94) b
34 707 10135 10537 41036 872992 393937 10280 100238 878 23739 79 4240 4 3 a The largest value for |M SHS | has been taken from Refs. [36–38]. b Estimated. cluster is not possible for sticky hard spheres; all vertices ofa regular icosahedron with unit edge length lie on a circum-scribing sphere with radius r c ≈ . m , n ) at which the icosahedral N =
13 LJ cluster breaks sym-metry to form a rigid cluster. For the n = m case consideredabove, this symmetry-breaking occurs at m (cid:39)
15. We also explored a more realistic extended LJ potential(Eq. 5; Figure 1) for one of the rare gas dimers (xenon) incomparison with other LJ potentials. We see that the repul-sive part agrees nicely with the conventional (6,12)-LJ po-tential, while for r > |M SHS → ELJ | = { , , , , , , } for N = { , , , , , , } . For N =
13 the number of distinct min-ima is 44% larger than it is for the simple (6,12)-LJ potential,which shows that |M ( N ) | is rather sensitive to the potentialchosen. Hence, to correctly describe the topology of real sys-tems, one has to take care of the correct form of the 2-bodycontribution (as well as higher n -body contributions) [25]. B. (6,12)-Lennard-Jones clusters from basin-hopping
Table I shows the number of distinct minima found byour cluster geometry optimisation procedure using the (6,12)-LJ potential compared to results from exact enumeration forSHSs and from basin-hopping for the (6,12)-LJ potential. Asthe SHS clusters for a specific N value can be grouped bytheir contact number N c , the geometry optimisations were car-ried out separately for each group of M SHS ( N c ). Hoy [36, 37]and Holmes-Cerfon [38] reported slightly di ff erent results for N =
11 and N =
13; we find that upon geometry optimisation,their datasets yield the same final clusters |M SHS → LJ ( N c ) | . Asidentical LJ clusters appear in multiple groups with di ff erentcontact numbers, we remove the duplicates to create the set M SHS → LJ of distinct minima, which can be directly comparedto the set of LJ minima M LJ obtained from the basin-hoppingmethod. It should be noted that including the hypostatic clus-ters and the di ff erent |M SHS | for N =
11 and N =
13 fromRef. [38] did not change our results, implying that hypostaticclusters are not an important feature for the LJ energy land-scape.Interestingly, our gradient-based minimisation procedurestarting from the SHS packings does not in general lead to acomplete set of LJ minima; the mapping from SHS minima toLJ minima is non-injective and non-surjective. Clearly, somestructural motifs found in LJ clusters are not found in SHSclusters and vice versa, and the topology of the hypersurfacechanges in a non-trivial fashion from SHS to LJ. However, itis surprising that the fraction of structures that are missed bythis optimisation procedure is so small (see Table III). To gainfurther insight, we analysed the energetics and structure of theunmatched clusters in more detail.Figure 6 shows an analysis of the di ff erence between thelongest to the shortest bond lengths d ∆ = d max − d min obtainedfor the largest clusters in M LJ with N = { , , } [64]. Thehistograms show that the clusters most commonly have a d ∆ of about 0 .
03. In contrast, as shown by the orange bars, theunmatched structures have significantly larger d ∆ values of atleast 0 .
05, with most of them having d ∆ (cid:39) .
06. This is afirst indication of why these structures are not found by start-ing from SHS packings. The latter only form bonds of length f r e q u e n c y d ∆ M SHS → LJ M LJ − M SHS → LJ (a) N = f r e q u e n c y d ∆ M SHS → LJ M LJ − M SHS → LJ (b) N = f r e q u e n c y d ∆ M SHS → LJ M LJ − M SHS → LJ (c) N = FIG. 6: Histograms of the di ff erence between the longest and shortest bond distances d ∆ = d max − d min for the complete set ofdistinct LJ minima M LJ ( N ) for N = { , , } . Orange bars give the number of distinct structures not contained in M LJ asobtained from the basin-hopping algorithm.TABLE II: Range [ E , E max ] of the energy spectrum of all LJminima, position of the second lowest minimum structure E and position of the first unmatched (UM) structure E UM0 relative to the respective global minimum (in reduced unitsand E = N E max E E UM0 one, and a large variation in bond length could imply that aSHS packing similar to the LJ structure does not exist as theSHS boundary conditions are not satisfied. The data in Ta-ble II show that the unmatched (UM) structures for a specific N value have much higher energies compared to the one of theglobal minimum (which is set to zero, i.e. E = d ∆ and the energetic positionof the LJ clusters.Last, we checked the geometries of the missing structuresin more detail. As it turns out, almost all of the missing sta-ble LJ clusters can be created from a smaller set of missingclusters by capping some of their triangular faces. Therefore,these groups of clusters can be referred to as “seeds” [35].The corresponding starting structures of each seed are shownin Figure 7. None of these structures are stable SHS packings.For example, structure (d) can be described as three octahedraconnected via triangular faces sharing one edge. Geometricconsiderations [35, 36] immediately show that this structurecannot be a stable SHS packing; the dihedral angle in an octa-hedron is approximately 109 . ◦ , which means three octahedraonly fill 328 . ◦ of a full circle, leaving a gap between twofaces.Table III shows the number of missing minima belonging TABLE III: Number of missing structures after optimisationbelonging to the same ”seed” (Fig. 7). N = seed N = N = N = N = N =
13a 1 1 - 3 8b - 1 3 4 12 a c - - 1 1 a -d - - 1 1 5e - - - 1 6f - - - 1 1remaining - - - - 2total 1 2 5 11 34% 4.76 3.13 2.94 2.14 2.25 a Some structures do not resemble a perfect capped cluster, but undergo aslight rearrangement. Specifically, two structures belonging to seed (b)and one structure belonging to seed (c) were found to deviate slightlyfrom the perfect arrangement, but minor rearrangements of thesestructures lead to the desired geometry and they can be reasonablyassociated with these seeds. to each seed. Over 60 % of the unmatched structures belongto seeds (a) and (b). From a graph theoretical point of view[34, 35], grouping structures into seeds means that all struc-tures belonging to the same seed contain the graph of the start-ing structures as a subgraph in their respective connectivitymatrix. This approach simplifies the analysis to a great extent,as the feature that prevents the structures from being foundby geometry optimisation is the same for each of the struc-tures arising from a specific seed. The smallest unmatchedstructures that cannot be associated with any of seeds (a)-(f)have N =
13; these could be the starting structures for two newseeds.Finally, we note that the starting SHS minima in our op-timisation procedure are not stationary points on the LJ hy-persurface, and we therefore optimise to most but not all lo-cal and available LJ minima. This observation explains whysome high-energy structures were not found by our optimisa-tion procedure. For a smooth change in the topology of the (a) (b) (c)(d) (e) (f)
FIG. 7: Graphical representations of the structures that arestarting new seeds, but are not contained in M SHS → LJ . SeeTable III and text for more details.potential energy surface from SHS to LJ type clusters one hasto continuously vary the exponents ( n , m ) in real space, whichis computationally too demanding. C. The special case of a Gregory-Newton cluster
We call a cluster “Gregory-Newton” (GN) when it belongsto the set of all clusters consisting of 12 spheres kissing acentral sphere. The canonical Gregory-Newton cluster is theicosahedron, which is perhaps the most common form studiedin cluster chemistry and physics [2, 41, 65, 66]. We thereforeinvestigate this cluster type in more detail here.For monodisperse SHS clusters, the Gregory-Newton argu-ment (as proved by Sch¨utte and van der Waerden [40]) that nomore than 12 equally sized spheres can touch a central sphereof same size holds. We note that the problem of the number ofkissing spheres in k dimensions, or even in three dimensionswith sphere size smaller than that of the central sphere, re-mains largely unsolved [41, 67]. For unequally sized spheres,some simple results related to spherical codes [68] are known;for example, 13 hard spheres of radius r s can touch a centralsphere of unit radius only if r s ≤ . V LJmn ( r ), however,the situation is far more complicated since systems minimizeenergy rather than di ff erences in the distances between neigh-boring particles, and few general results are known. Nonethe-less, this latter problem is important for understanding realsystems such as coordination compounds [69], which have re-cently been shown to possess coordination numbers as high as17 [70] or even 20 [71].Motivated by these recent results, we investigated longerrange potentials by decreasing the LJ exponents ( m , n ), to seewhether the restriction of no more than 12 kissing equal-sizedspheres still holds. As it is impossible to distribute 13 points on a sphere evenly (there is no triangulation of a sphere with13 vertices of degree 5 and 6 [72]), we used the Fibonaccisphere algorithm [73, 74] to find an approximate distributionof points on a sphere and added a center sphere. By optimisingthe coordinates for this N =
14 cluster with di ff erent LJ expo-nents and calculating the distance of every sphere to the cen-ter sphere, we can deduce at which “softness” a 13th sphereis (perhaps) allowed to enter the first coordination shell, i.e. totouch the center sphere.Figure 8 shows the di ff erence between the largest and thesmallest center-to-outer sphere (COS) distances in relation tothe LJ exponents m and n . Interestingly, none of the ( m , n )-LJ potentials lead to equal distances around a central sphere.While this result could be due to the lack of symmetry, onesphere is clearly further away from the central sphere even forthe softest “Kratzer” (1,2)-LJ potential [75]. For this poten-tial the largest and smallest COS distances are r max = . r min = . r = .
845 and r = . r max / r min ratio is 1.097 and much shortercompared to r max / r min = √ r GN14 = .
347 (see discussion below). Hence the 13th sphere“almost” touches the center sphere.Note that all COS distances for the N =
14 (1,2)-LJ clus-ter are significantly shorter than r =
1, due to the N ( N − / V LJ mn ( r ) with n > m >
3, onecan prove [23] that the nearest neighbor distance is r NN ( m , n ) = (cid:16) L n L − m (cid:17) n − m . (8)Here L n is the Lennard-Jones-Ingham lattice coe ffi cient for aspecific lattice determined from 3D lattice sums. Since L n < L m for n > m , we see that r NN <
1, and lim m , n →∞ r NN ( m , n ) = r min ( N ) are: r min (8) = . r min (9) = . r min (10) = . r m in (11) = . r min (12) = . r min (13) = . r min (12) is smaller than r NN (6 , r NN (6 ,
12) values are 0 . . . N =
13 and N =
14 SHS clusters from Ref. [38]. This set contains all non-isomorphic SHS structures that can be considered GN clusters( N =
13) and the N =
14 structures that can be derived fromthem by attaching a 14th sphere. We find a surprisingly largenumber (737) of nonisomorphic N =
13 GN-SHS structures( { , , , } for N c = { , , , } ), that all optimise tothe ideal icosahedral arrangement ( I h symmetry) if a (6,12)-LJpotential is applied. An even larger number of clusters existsfor N =
14 (14529), which is ≈ . |M SHS (14) | . All of these m n FIG. 8: Relation of LJ exponents m and n to the di ff erence oflargest and smallest center-to-outer sphere (COS) distances.A value of zero would imply that all surrounding spheres aretouching the center sphere.structures optimise to just one of two possible (6,12)-LJ min-ima of GN type. The first is the Mackay icosahedron capped atone of its triangular faces, and the second is an elongated pen-tagonal bipyramid (belonging to the class of Johnson solids)with the 14th sphere capping a square face.Most of these N =
14 clusters are minimally rigid ( N c = N − = N c > N −
6) and none are hypostatic ( N c < N − { , , , , } such clusters with N c = { , , , , } and N =
14. The clusters with N c =
40 are hcp and fcccore-shell structures capped at a square face; these arrange-ments maximise N c . Most of the clusters with N c = { , } are deformed versions of the elongated pentagonal bipyra-mid mentioned above, indicating that this arrangement is afavored route to these intermediate-energy structures. How-ever, N c =
39 also contains hcp and fcc structures capped at atriangular face. The first example of a cluster derived from aperfect icosahedral symmetry shows up at lower value N c = N =
14 cluster with the closest central-to-outer sphere (COS) distance r COSmin was not known. Herewe close this gap by determining the COS distance for allGregory-Newton type clusters. We find one single cluster with r COSmin = . r COS cluster ( r COS = . (a) r GN14 = . N c =
39 (b) r GN14 = . N c = r GN14 = √ N c =
40 (d) r GN14 = (cid:113) , N c = FIG. 9: Graphical representations of SHS packings with N =
14, where a center sphere is maximally contacting. Theorange sphere in each cluster is the 14th outer sphere, notable to touch the center sphere (in black). (a) distortedelongated pentagonal bipyramid (Johnson solid); (b)distorted icosahedron; (c) hcp capped on a square; (d) hcpcapped on a triangle. f r e q u e n c y r COS
FIG. 10: Frequency of distances from the cluster center to themost distant sphere for all Gregory-Newton-like clusterscontained in the structures from Ref. [38]. The width of thebars is 0 . r COS values forthe full set of GN clusters is shown in Figure 10. Motifs withlarger r COS are far more prevalent. For example, the peak at r COS = .
41 corresponds to structures where the 14th sphereis touching 4 other spheres that are part of a tetragonal pyra-mid, therefore forming a regular octahedron with a tip-to-tipdistance of √ r COS value (1 . √ / . , .
58 and1 .
55 are derived from the regular trigonal bipyramid and re-sult from breaking its axial bonds. In these structures, themore bonds are broken, or the further the axial spheres areseparated, the shorter the center-to-outer sphere distance be-comes.
IV. CONCLUSIONS
We have characterized the sets of ( m , n )-LJ-potential min-ima obtained using complete sets of nonisomorphic SHSpackings with 8 ≤ N ≤
14 [34–38] as initial states for energyminimization. The number of distinct minima (i.e. excludingpermutation-inversion isomers) is far smaller than the numberof SHS packings for the standard Lennard-Jones exponents( m , n ) = (6 , m , n ) increase. We characterized how the number of distinctminima M ( N ) increases with cluster size N by determiningStillinger’s rise rate parameter α (Eq. 6 [44]). The increase of α from ≈ . ≈ . m , n )-LJ energy landscape towards the SHS energy land-scape as ( m , n ) increase.Using a more realistic extended LJ potential obtained fromcoupled cluster calculations for the xenon dimer [23, 48] leadsto M values close to those obtained for the (6,12)-Lennard-Jones potential, but our results indicate the the topology ofthe energy hypersurface is very sensitive to the model poten-tial applied. For softer potentials, we showed that it is still un-favourable for a 13th outer sphere to touch the center sphere.Indeed, the Gregory-Newton argument still holds true for eventhe softest ( m , n ) = (1 ,
2) potential.Finally, we compared our optimisation results to thepreviously published results for the (6,12)-LJ potential.The mapping from M SHS to M SHS → LJ is non-injective andnon-surjective, however, the number of structures missed bythe optimisation procedure is relatively small. The unmatchedstructures belong to the high energy region of the potentialenergy hypersurface and possess rather large variations intheir bond lengths. An analysis of their geometries revealedthat most of the larger structures can be constructed from asmaller cluster by capping some of the triangular faces. Thisprocedure e ff ectively sorts almost all unmatched structuresinto six seeds for clusters up to N = V. ACKNOWLEDGEMENTS
We acknowledge financial support by the Marsden Fund ofthe Royal Society of New Zealand (MAU1409). DJW grate-fully acknowledges financial support from the EPSRC. PSacknowledges financial support by the Centre for AdvancedStudy at the Norwegian Academy of Science and Letters(Molecules in Extreme Environments Research Program). Wethank Drs. Lukas Wirz and Elke Pahl for useful discussions. [1] F. H. Stillinger and T. A. Weber, Science , 983 (1984).[2] T. Martin, Phys. Rep. , 199 (1996).[3] J. P. K. Doye and D. J. Wales, Science , 484 (1996).[4] E. Vlieg, M. Deij, D. Kaminski, H. Meekes, and W. van Enck-evort, Faraday Discuss. , 57 (2007).[5] G. Meng, N. Arkus, M. P. Brenner, and V. N. Manoharan, Sci-ence , 560 (2010).[6] C. R. A. Catlow, S. T. Bromley, S. Hamad, M. Mora-Fonz, A. A.Sokol, and S. M. Woodley, Phys. Chem. Chem. Phys. , 786(2010).[7] S. Karthika, T. K. Radhakrishnan, and P. Kalaichelvi, CrystalGrowth & Design , 6663 (2016).[8] M. Holmes-Cerfon, Ann. Rev. Cond. Matt. Phys. , 77 (2017).[9] P. R. Rowland, Discuss. Faraday Soc. , 364 (1949).[10] P. R. Unwin, Faraday Discuss. , 409 (2007).[11] A. R. Oganov and C. W. Glass, The Jour-nal of Chemical Physics , 244704 (2006),https: // doi.org / / , 037101 (2007).[13] D. J. Wales, ChemPhysChem , 2491 (2010). [14] A. R. Oganov, A. O. Lyakhov, and M. Valle, Accountsof Chemical Research , 227 (2011), pMID: 21361336,http: // dx.doi.org / / ar1001318.[15] F. Calvo, J. P. K. Doye, and D. J. Wales, Nanoscale , 1085(2012).[16] D. J. Wales, J. Chem. Phys. , 130901 (2015).[17] B. W. van de Waal, J. Chem. Phys. , 3407 (1989).[18] C. L. Cleveland and U. Landman, J. Chem. Phys. , 7376(1991).[19] B. W. van de Waal, Phys. Rev. Lett. , 1083 (1996).[20] J. P. K. Doye, D. J. Wales, and R. S. Berry, J. Chem. Phys. ,4234 (1995).[21] B. W. van de Waal, G. Torchet, and M. F. de Feraudy, Chem.Phys. Lett. , 57 (2000).[22] N. V. Krainyukova, R. E. Boltnev, E. P. Bernard, V. V. Khme-lenko, D. M. Lee, and V. Kiryukhin, Phys. Rev. Lett. ,245505 (2012).[23] P. Schwerdtfeger, N. Gaston, R. P. Krawczyk, R. Tonner, andG. E. Moyano, Phys. Rev. B , 064112 (2006). [24] N. V. Krainyukova, The European Physical Journal D , 45(2007).[25] P. Schwerdtfeger, R. Tonner, G. E. Moyano, and E. Pahl,Angew. Chem. Int. Ed. , 12200 (2016).[26] H. Cox, R. L. Johnston, and J. N. Murrell, J. Sol. St. Chem. , 517 (1999).[27] S. Y., O. K., T. T., and O. M., Scientific Reports , 13534(2015).[28] C. Leitold and C. Dellago, J. Chem. Phys. , 074504 (2016).[29] M. B. Sweatman and L. Lue, J. Chem. Phys. , 171102(2016).[30] R. J. Baxter, J. Chem. Phys. , 2770 (1968).[31] S. B. Yuste and A. Santos, Phys. Rev. E (1993).[32] T. W. Cochran and Y. C. Chiew, J. Chem. Phys. , 224901(2006).[33] M. Holmes-Cerfon, S. J. Gortler, and M. P. Brenner, Proc. Nat.Acad. Sci. , E5 (2013).[34] N. Arkus, V. N. Manoharan, and M. P. Brenner, Physical Re-view Letters , 1 (2009).[35] N. Arkus, V. Manoharan, and M. Brenner, SIAM Journal onDiscrete Mathematics , 1860 (2011).[36] R. S. Hoy, J. Harwayne-Gidansky, and C. S. O’Hern, PhysicalReview E , 051403 (2012).[37] R. S. Hoy, Phys. Rev. E , 012303 (2015).[38] M. C. Holmes-Cerfon, SIAM Rev. , 229 (2016).[39] Y. Kallus and M. Holmes-Cerfon, Phys. Rev. E , 022130(2017).[40] K. Sch¨utte and B. L. van der Waerden, Math. Ann. , 325(1952 / Sphere packings, lattices andgroups , Vol. 290 (Springer Science & Business Media, 2013).[42] K. Bezdek and S. Reid, J. Geom. , 57 (2013).[43] P. Erd¨os, Am. Math. Mon. , 248 (1946).[44] F. H. Stillinger, Phys. Rev. E , 48 (1999).[45] Y. Forman and M. Cameron, J. Stat. Phys. , 408 (2017).[46] S. Heiles and R. L. Johnston, International Journal of QuantumChemistry , 2091 (2013).[47] J. W. R. Morgan, D. Mehta, and D. J. Wales, Phys. Chem.Chem. Phys. , 25498 (2017).[48] P. Jerabek, O. Smits, E. Pahl, and P. Schw-erdtfeger, Molecular Physics , 1 (2018),https: // doi.org / / , 8417 (1999). [50] P. A. Braier, R. S. Berry, and D. J. Wales, J. Chem. Phys. ,8745 (1990).[51] D. J. Wales, Science , 2067 (2001).[52] D. J. Wales, Energy Landscapes: Applications to Clusters,Biomolecules and Glasses (Cambridge Molecular Science,2004).[53] “Pele: Python energy landscape explorer,” https://github.com/pele-python/pele (2017).[54] Z. Li and H. A. Scheraga, Proc. Natl. Acad. Sci. USA , 6611(1987).[55] D. J. Wales and H. A. Scheraga, Science , 1368 (1999).[56] D. J. Wales and J. P. K. Doye, The Journal of Physical Chem-istry A , 5111 (1997).[57] J. P. K. Doye and D. J. Wales, J. Chem. Phys. , 3777 (2002).[58] D. E. King, Journal of Machine Learning Research , 1755(2009).[59] D. C. Wallace, Phys. Rev. E , 4179 (1997).[60] F. Sciortino, W. Kob, and P. Tartaglia, Phys. Rev. Lett. , 3214(1999).[61] M. R. Hoare and P. Pal, Adv. Phys. , 645 (1975).[62] M. R. Hoare and J. McInnes, Faraday Discuss. Chem. Soc. ,12 (1976).[63] M. R. Hoare, Adv. Chem. Phys. , 49 (1979).[64] We define spheres that have a equilibrium distance between0 . − . , 916 (1962).[66] J. Uppenbrink and D. J. Wales, J. Chem. Soc., Faraday Trans. , 215 (1991).[67] O. R. Musin and A. S. Tarasov, Discrete & Computational Ge-ometry , 128 (2012).[68] C. Phillips, A. Jankowski, M. Marval, and S. Glotzer, Phys.Rev. E (2012).[69] A. Hermann, M. Lein, and P. Schwerdtfeger, Angew. Chem.Int. Ed. , 2444 (2007).[70] N. Kaltsoyannis, Angew. Chem. Int. Ed. , 7066 (2017).[71] T. D. Della and C. H. Suresh, Phys. Chem. Chem. Phys. ,14588 (2016).[72] P. Schwerdtfeger, L. N. Wirz, and J. Avery, Wiley Interdis-ciplinary Reviews: Computational Molecular Science , 96(2015).[73] ´A. Gonz´alez, Mathematical Geosciences , 49 (2010).[74] B. Keinert, M. Innmann, M. S¨anger, and M. Stamminger, ACMTrans. Graph. , 193:1 (2015).[75] A. Kratzer, Z. Phys.3