Calculating the symmetry number of flexible sphere clusters
Emilio Zappa Miranda Holmes-Cerfon





1]Emilio Zappa
2]Miranda Holmes-Cerfon

[1]Mathematics Department, Fordham University, NY
[2]Courant Institute of Mathematical Sciences, New York University, NY.


Calculating the symmetry number of flexible sphere clusters
    [
    
===========================================================





We present a theoretical and computational framework to compute the symmetry number of a flexible sphere cluster in ℝ^3, using a definition of symmetry that arises naturally when calculating the equilibrium probability of a cluster of spheres in the sticky-sphere limit. 
We define the sticky symmetry group of the cluster as the set of permutations and inversions of the spheres which preserve adjacency and can be realized by continuous deformations of the cluster that do not change the set of contacts or cause particles to overlap. The symmetry number is the size of the sticky symmetry group. We introduce a numerical algorithm to compute the sticky symmetry group and symmetry number, and show it works well on several test cases. 
Furthermore we show that once the sticky symmetry group has been calculated for indistinguishable spheres, the symmetry group for partially distinguishable spheres (those with non-identical interactions) can be efficiently obtained without repeating the laborious parts of the computations. 
We use our algorithm to calculate the partition functions of every possible connected cluster of 6 identical sticky spheres, generating data that may be used to design interactions between spheres so they self-assemble into a desired structure. 





§ INTRODUCTION


Symmetry plays an important role in the study of molecules and clusters of more general particles.  A group of particles, like atoms or colloids, can assemble into a variety of different clusters, and to determine which of these forms requires evaluating the number of symmetries of each cluster  <cit.>. 
The definition of a symmetry varies depending on the  
particular system, quantities of interest, and method of calculating the partition function from statistical mechanics, but for a great many clusters – those that are close to rigid, i.e. they don't move far from a reference configuration – the definition coincides with the geometrical symmetries of the configuration, such as a symmetry upon reflection or upon rotation by certain angles around certain axes. Much effort has gone into developing methods to evaluate the geometrical symmetry group of molecules, an effort which is ongoing since this computation is a challenge for large molecules <cit.>. 



When a cluster is flexible, i.e. it can deform by some non-negligible amount along its internal degrees of freedom without a significant change in energy, the geometrical symmetry of the cluster is less meaningful, because the cluster almost always adopts a configuration with no or few true geometrical symmetries. However, the concept of symmetry is still meaningful, as long as one extends symmetry elements to include operations that continuously transform the cluster to one with a similar energy, called  “feasible” transformations in the seminal paper by Longuet-Higgins <cit.>. 
This is a useful definition of symmetry because one can show that the number of such symmetries equals the number of times one has overcounted the cluster for the most common methods of evaluating its partition function. 
Calculating this kind of symmetry group requires more sophisticated methods since it requires understanding the kinetic pathways the cluster can follow, and not just its static geometrical symmetries. 
Some methods to compute such symmetry elements find pathways by searching a database of local minima and transition states on a cluster's energy landscapes, a process which can work well for atomic clusters, whose energy changes smoothly with the configuration <cit.>.
Other methods are based on identifying a set of generators for the symmetry group, and then building the full group using group-theory software <cit.>. 
These methods often start with a particularly symmetric reference configuration with many geometical symmetries, or else identify symmetries by eye <cit.>. 
Sometimes it is even possible to prove more general symmetry results about a cluster, and to extend these proofs to the more restrictive symmetry groups obtained when an energy function is present, for example as in <cit.>.
Because identifying generators by eye or by first finding a particularly symmetric configuration, or proving statements about individual clusters, requires an external observer to provide input, 
these latter methods cannot be used to automatically evaluate the symmetries of a large collection of clusters. 

We are interested in particles with diameters of nano- to micrometres (colloids), which are much larger than atoms, and form the building blocks for a wide range of materials <cit.>. 
Such particles interact attractively over scales typically much smaller than their diameters, so it is effective to model them in the sticky limit, where the particles are treated as hard, classical rigid bodies that can't overlap, such as spheres, and the interaction potential is a delta function at the point of contact between a pair of bodies <cit.>. Therefore, the energy of a cluster of particles is proportional to the number of pairs of particles that are exactly in contact. In contrast to atomic systems, this energy function changes abruptly at discrete locations in configuration space, so concepts developed for smoother energy landscapes, such as local minima and transition states, are no longer as meaningful. Therefore, algorithms for calculating the symmetry number of molecules with smooth interaction potentials will not work directly for these clusters. 

Our goal is to clarify the concept of symmetry for sticky-sphere clusters, possibly flexible, and to provide an algorithm that can evaluate the symmetry number of a sticky-sphere cluster automatically. 
We are interested both in indistinguishable particles and particles which can be partially distinguished. The theory developed here is a natural application of the theory of molecular symmetries from chemical physics, which is concerned with the types of symmetries present in a system's Hamiltonian. However,  this theory is usually presented in the context of quantum mechanics, and is not always straightforward to adapt to a purely classical setup <cit.>.  
 Our goal therefore is to present the theory of symmetry for sticky-sphere clusters in a mathematical framework, highlighting those aspects of the theory that depend only on geometry or topology and the overall connection to group theory. We hope this presentation will be accessible to those without a background in physics or chemistry, and will make it easy to adapt the theory to other situations where a modified definition of symmetry is needed. That the concept of symmetry for classical colloidal clusters needs clarification is evident in light of the numerous recent papers that attempt to explain it and the link to entropy in the statistical mechanics of such systems (e.g. 
<cit.>.)

Our theory starts with an equivalence relation between clusters, which says that two clusters are the same if one cluster can be continuously deformed into another without breaking any contacts or causing spheres to overlap. 
We define the sticky symmetry group of a sphere cluster to be the set of permutations or reflections of a cluster that can be achieved by such a continuous deformation, and the symmetry number to be the size of this group.  
This is a natural extension of the definition for rigid clusters, and we show how the sticky symmetry group is related to other symmetry groups commonly studied, the point group of the cluster and the automorphism group of its corresponding adjacency graph. 
Among the symmetry groups for flexible molecules studied in the literature, our approach is most closely related to the topological group studied by Flapan <cit.>. The difference is that the symmetry group she considers arises from purely topological properties of the cluster's contact graph, while our symmetry group contains some geometry, since we require spheres to not overlap. 

We introduce a numerical algorithm to automatically compute the sticky symmetry group of a sphere cluster. 
The main component of the algorithm is a numerical method to find a continuous deformation linking one cluster to another, based on the steepest descent method in optimization <cit.> and simulated annealing <cit.>. 
 We apply this algorithm to several test cases and show it works well for small clusters. In addition we show that once the sticky symmetry group has been calculated for indistinguishable spheres, it can be efficiently computed when some groups of spheres are distinguishable from others, for any possible partitioning of the spheres into distinguishable groups. 
 
Ultimately, we wish to use the algorithm in an exhaustive computation of the partition functions of small clusters of sticky spheres (which depend on the symmetry numbers), and then ask how to design interactions between spheres so they self-assemble into a desired structure. Toward this aim, we compute the symmetry numbers and partition functions of every connected cluster of N=6 spheres, and comment on some of the physical insight this data gives us. 

The outline of the paper is as follows. In section <ref> we give an overview of how the symmetry number arises when studying the statistical mechanics of sticky-sphere clusters. In Section <ref> we develop the mathematical theory required to define a symmetry number for a flexible sticky-sphere cluster. In Section <ref> we introduce a numerical algorithm to compute the symmetry number for a flexible cluster, and in Section <ref> we apply our theory and algorithm to several examples, including the exhaustive calculation for N=6 spheres.  Section <ref> concludes and discusses further applications in which this algorithm may be used. 




  
Mathematical setup. A cluster of N spheres is a pair (,), where = (_1, …, _N) ∈^3N is the vector of sphere centers, with _i ∈^3  the center of i-th sphere, and = (r_1,…,r_N) is the vector of sphere radii, with r_i>0 the radius of the ith sphere. 
We suppose that m spheres are in contact, given by the set E = { (i_1, j_1), …, (i_m,j_m)}. 
When two spheres are in contact, their centers are related by the equation

    |_i-_j|^2 = (r_i+r_j)^2 ,    (i,j) ∈ E .


We assume that all non-contacting pairs of spheres do not overlap, so their centers must satisfy the inequality 

    |_i -_j|> r_i+r_j,    (i,j) ∉ E .

It is sometimes convenient to represent the set of contacts by the adjacency matrix A,  an N × N matrix whose entries are given by

    A_ij = {[            1 if (i,j) ∈ E;            0 if (i,j) ∉ E ].
     .

Notice that the pair (,) is sufficient to characterize the cluster, since from it we can determine the set of contacts E and hence the adjacency matrix A. 

We define M_A^() to be the set of all configurations  with adjacency matrix A and radii , i.e. the set of points in ^3N which satisfy (<ref>), (<ref>):

    M_A^() = {∈^3N : q_ij() = (r_i+r_j)^2,  if  A_ij = 1,   q_ij() > (r_i+r_j)^2,  if  A_ij = 0 } ,

where q_ij() =  |_i-_j|^2. 
In the following, we will usually drop the dependence on the radii , and write M_A for M_A^() and call a cluster , since  doesn't change for a given problem. 
 If the gradients {∇ q_ij() }_(i,j)∈ E are linearly independent for every ∈ M_A, then M_A is a manifold of dimension p = 3n-m.  
In general, M_A is not connected, as we illustrate with examples in Section <ref>. We denote by M_A, the connected component of M_A to which cluster  belongs.  

The cluster  is rigid if every point  in the connected component M_A, can be obtained as a rotation or translation of , otherwise it is non-rigid or flexible <cit.>. In the latter case, the cluster has internal degrees of freedom, other than translations and rotations: it can be continuously deformed without breaking contacts.  







§ OVERVIEW OF THE SYMMETRY NUMBER AND PARTITION FUNCTION FOR A STICKY-SPHERE CLUSTER




We now give a brief overview of the partition function for a sphere cluster in the sticky-sphere limit, and  explain why the symmetry number enters. Readers interested in only in the mathematical definitions may skip to section <ref>. 

We start with a cluster  of identical spheres with adjacency matrix A with sticky interactions, and we wish to know the probability it will occur in equilibrium. We consider a system where we cannot tell  apart from a  cluster obtained by performing any of the following transformations: (i) rotating it, (ii) translating it,  (iii) deforming it along internal degrees of freedom without changing its adjacency matrix, (iv) reflecting it about an axis, or (v) permuting the particles. 
 In a classical system, such indistinguishability could arise because we cannot experimentally distinguish such clusters, or it could be imposed artificially even if we can distinguish them because we are only interested in properties of a group of clusters, and not individual clusters within the group. (In a quantum system, particles may be truly indistinguishable in a sense not possible in classical mechanics, but we do not consider quantum effects here.)
The equilibrium probability of finding  or any of its transformed versions is proportional to the partition function Z_, namely the integral of the Boltzmann distribution over the region in configuration space corresponding to all the possible transformations above of .[
Specifically, given the partition function for a collection of distinct (not related by any of the 5 transformations) clusters 1,2,…,k, the equilibrium probability to observe a cluster i (or a transformation of it), given that at least one of the set above was observed, is Z_i / ∑_j=1^kZ_j. 
] 
For sticky-sphere interactions where all pairwise interactions are identical, and when the constraints {∇ q_ij() }_(i,j)∈ E are linearly independent[
When these constraints are not linearly independent, a similar expression may still be available but requires an additional parameter beyond κ, see <cit.>.
] for every ∈ M_A,, the partition function can be written as <cit.>

    Z_ = C  n_ κ^m ∫_M_A,  f(') dσ(') .

Here κ>0 is a system-dependent constant depending on the properties of the interaction potential and temperature, 
f() is a function depending only on the relative locations of the spheres (physically it is related to the vibrational entropy), 
 dσ() is the natural surface measure on M_A, induced by restriction from the ambient Euclidian space, which also only depends on the relative locations of the spheres, and C is a constant that is the same for all clusters with the same number of spheres N (provided the space they live in is large enough that excluded volume effects do not matter.) 
We usually put additional restrictions on the cluster to make the integral (<ref>) finite, for example by assuming the cluster is contained in a large box, or by fixing its center of mass. 

The remaining quantity in (<ref>) is n_, which is the number of disconnected copies of M_A, that we obtain by permuting or reflecting particles. 
Usually this factor is written as 2N!/σ, where σ is the so-called symmetry number of the cluster. Because the definition of symmetry is not fundamental, but depends on the method of calculating the partition function, we prefer to call n_ the counting number, to highlight that it counts the number of geometrically isomorphic copies of M_A, obtained by such transformations. 


To see how the expression in (<ref>) accounts for the five transformations listed above, notice that 
 M_A, contains all possible transformations of  available through rotations, translations, and deformations along the internal degrees of freedom of a cluster. Indeed, these transformations don't alter the constraints in (<ref>) defining the manifold, and they are continuous so they take  to another point in M_A,. We only need to explicitly account for reflections and permutations, which can sometimes put  on a different connected component of M_A, or a different manifold M_A' for A≠ A'. This accounting is done through the counting number. 

To see how, it is easiest to first consider a rigid cluster. 
An upper bound on the counting number is obtained by counting the total number of permutations and reflections: there are N! ways of permuting the particles, and for each permutation we can reflect or not the cluster through some plane, so there are up to 2N! copies of the original cluster. 
This set of permutations and reflections are called permutation-inversion (PI) operations <cit.>. 
However, not all of these 2N! copies are distinct: some of these copies might correspond to a rotation or a reflection of the cluster, which we have already accounted for via the integral over M_A,.  
The number of such repeated copies is defined to be the symmetry number σ of the cluster, because each repeated copy corresponds to a geometrical symmetry element in the point group of the cluster <cit.>. 
The number of distinct copies of M_A, available by permutations and reflections is then n_ = 2N!/σ.

For clusters with the same number of contacts, the symmetry number can play an important role in determining the relative probabilities, as (<ref>) shows that any ratio of symmetry numbers between clusters affects their relative probabilities by an equal ratio. 
Indeed, experiments have shown that in an ensemble of 6 particles interacting with a short-range depletion interaction, the least symmetric cluster (the polytetrahedron, with symmetry σ=4, see Figure <ref>) occurs 24 times more frequently than the most symmetric cluster (the octahedron, with σ=48, see Figure <ref>), and most of this discrepancy is due to the ratio of symmetry numbers <cit.>. 

For a flexible cluster, the counting number must account for the fact that the cluster has internal degrees of freedom. Therefore, even if none of the 2N! permutation-inversion operations correspond to a perfect geometrical symmetry of , they may still transform the cluster so it lies on the same connected component M_A,, and therefore has already been accounted for in the integral over M_A,. 
For example, consider a chain of four identical spherical particles, labelled 1–4, as in Figure <ref>. This particular embedding  of the chain is purposely chosen so the chain has no geometrical symmetries; its point group is trivial. However, the permutation (14)(23), which reverses the ordering of the spheres, results in a cluster that can be obtained from the original one by a continuous deformation, for example by straightening the chain, rotating it by 180^∘, then crinkling it up again. The permuted cluster therefore lies on the same connected component M_A,.
In this case, to integrate over the correct space of transformations if we start with the integral expression in (<ref>),  we need to define the symmetry number σ to be the number of permutation-inversion operations that can also be realized by a combination of transformations (i–iii), namely translating, rotating, and internal deformations; none of these need correspond to an actual geometrical symmetry. 
The counting number is still the same, n_=2N!/σ. 

These ideas may be easily extended to the case of partially distinguishable spheres, i.e. where certain groups of spheres are distinguishable. 
Physically, this situation may arise when the particles have different kinds of interactions with each other, or are physically different in some way such as in the materials they are made with. Such situations are of interest in self-assembly problems where one wants to design interactions between particles to so they assemble into a desired structure <cit.>. 
For example, we may consider spheres in some set C_1 = {1,2} to be indistinguishable from each other, and those in set C_2={3,4,…,N} to be indistinguishable from each other, but we can determine whether a sphere belongs to set C_1 or C_2. 
It is convenient to think of the spheres as having different colors, so 
 we may call spheres in C_1 “blue” and spheres in C_2 “red.” In general there could be anywhere from 1 to N colors. 

When there are two or more colors, what changes in the partition function (<ref>) are the following: (i) the factor κ^m is replaced by ∏_(i,j)∈ Eκ_ij, a product of sticky parameters κ_ij for each pair (i,j) in contact in the cluster; and (ii) the counting number n_ must be recomputed to account for the number of different copies of M_A, one obtains under the more restrictive set of permutations which preserve the colors. Notably, the integral over M_A,, which is laborious to compute <cit.>,  does not change. An important contribution of this paper will be to show that the counting number for colored spheres can be efficiently obtained from the counting number for indistinguishable spheres. Therefore, using this framework one can easily compute the partition function after changing the interactions (sticky parameters) and colorings of the spheres, once the integral and counting number have been computed for indistinguishable spheres. 






§ COUNTING NUMBER, SYMMETRY NUMBER, AND STICKY SYMMETRY GROUP


This section more precisely elaborates on the ideas introduced in section <ref>. 
Our aim is to define the counting number and symmetry number of a cluster and show how they are related to the number of isometric copies of manifolds M_A,. 






We start by fixing the radii = (r_1, …, r_N) and defining the set Y  of all the physical realizations of sticky-sphere clusters with that set of radii:

    Y= ^3N∖{∈^3N : |_i-_j| < r_i+r_j,  ∃ i, j,   i≠ j }.

Within set Y the adjacency matrix of a cluster  A()=A(,) is well-defined. 


We introduce in Y an equivalence relation ∼ between clusters, that tells us which clusters are assumed to be “the same,”  written ∼. We define

    ∼  ⇔   A() = A(),  and ∃φ:[0,1]→ M_A, φ cts,  s.t. φ(0)=, φ(1)=.
 
It is trivial to verify that ∼ is an equivalence relation. 

In words, ∼ if these clusters have the same adjacency matrix, and there is a continuous deformation, formed from some combination of rotations, translations, and motion along internal degrees of freedom, from  to . 

By construction, if ∼ then they belong to the same connected component of M_A(), i.e. M_A(), = M_A(),. 

The set of clusters which can be obtained by a continuous deformation of  is the equivalence class of , written []. A trivial but important relation is that

    []=M_A, .

Two clusters are distinct if they belong to different equivalence classes: they either have different adjacency matrices, or there is no continuous deformation from one to the other. 

We denote by X the quotient set of Y modulo ∼:

    X = Y/∼.
 
X is the set of all distinct equivalence classes of clusters. 

We point out that each of Y, ∼, X depends on , but we suppress this dependence in the notation for brevity. 


We also need to consider the set of permutations and reflections of a cluster .
A reflection of  is simply -; any reflection about any other plane is then obtained by rotating a particular reflection. 
To construct a permutation of  we start with a permutation matrix P, an N × N matrix with entries P_ij=1 if i→ j after the permutation is applied, and P_ij=0 otherwise.  The permuted cluster is then (P ⊗ I_3), where I_3 is the 3 × 3 identity matrix, and ⊗ is the Kronecker product of matrices  <cit.>.  
For brevity we write P=P ⊗ I_3 so that P is the permuted version of . 

Here are two useful facts about a permutation matrix P: one,  P^-1 = P^T, since P is an orthogonal matrix, and two, if  the adjacency matrix for  is A, then the adjacency matrix for P is PAP^T. 

To make the link to symmetry groups, we need to consider how the group of permutation-inversion operations acts on []. 
It is simplest to first consider spheres with identical radii that are indistinguishable; we do this in Section <ref>. Then, we consider extensions to partially distinguishable spheres (Section <ref>) and spheres with different radii (Section <ref>). 




 §.§ Indistinguishable spheres with identical radii


Throughout this section we assume  is a cluster of N  indistinguishable spheres with identical radii. 



  
Basic results from group theory.
Consider the product group 

    P(N) × C_2 = { (P,δ) : P∈ P(N),  δ∈ C_2 },

where P(N) is the group of N × N permutation matrices, and C_2 = {± 1}. An element in P(N) × C_2 is called a permutation inversion (PI) operation <cit.>.  

We define an action of P(N) × C_2 on the quotient set X given in (<ref>) as 

    (P,δ) · []   =  
     [δP]  .

That is, the action permutes the spheres in  and possibly reflects it. 
The proof that · is a well-defined group action is given in the Appendix. 


We define the counting number n_ of a cluster ∈^3N to be the size of the orbit of [] with respect to the action ·, written ([]):

    n_ = |([])| =  |{ [δP ] : P ∈ P(N), δ∈ C_2 }|.

In words, the counting number equals
the number of distinct copies of [] that we obtain by permuting spheres in  or reflecting it. 

Some permutations or inversions of  leave its equivalence class unchanged. 
The set of such operations forms the stabilizer of [], ([]): 

    ([]) = { (P,δ) ∈ P(N) × C_2 : (P,δ) · [] = []}.

It is a fact of group theory that ([]) is a subgroup P(N)× C_2 <cit.>. 

The orbit-stabilizer theorem relates the sizes of the orbit and stabilizer of [] <cit.>: 

    |([])| = |P(N) × C_2|/|([])| .





  
Sticky symmetry group and relation to manifolds M_A(),.

We will now show how this theorem is related to the number of distinct copies of M_A(),, and consider a more explicit way to define the stabilizer. 

Because of (<ref>), 
the group action (<ref>) can also be thought of as acting on manifolds as

    (P,δ) · M_A(), = M_PA()P^T, δP .


Therefore, an element (P,δ) belongs to the stabilizer if and only if PA()P^T = A(),  and both  and δP belong to the same connected component of the manifold M_A().





This relation gives another way to define the stabilizer that will prove useful in computations. 

Recall that  the automorphism group  of an adjacency matrix A is the set of all permutation matrices P that preserve adjacency:

    = Aut(A) = { P ∈ P(N) : PA = AP }.

An element P∈ corresponds to a permutation that, when applied to spheres in a cluster, doesn't change who each sphere is in contact with. 
Notice that the automorphism group is independent of the embedding of  ∈^3N and the radii ; it is a property only of the graph associated with the adjacency matrix. 

Using this definition we can write the stabilizer more explicitly as:

    _ := ([]) = { (P, δ) ∈× C_2  : ∃φ : [0,1] → M_A, φ cts, : φ(0) = ,  φ(1) = δP} .




We call _ the  sticky symmetry group of the cluster. 
From this definition it is clear that _ is a subgroup of × C_2, a group we call the automorphism-inversion group.
In particular, _ is the set of the elements of the automorphism-inversion group  that can be obtained as either a continuous deformation of the cluster, or a deformation combined with a reflection. Notice that _ depends on the class [], in contrast to × C_2, which doesn't. 

We define the  symmetry number σ_ of the cluster  to be the cardinality of _:

    σ_ = |_|.


Combined with the orbit-stabilizer theorem (<ref>), we obtain the relationship 

    n_ = 2N!/σ_.


Because of (<ref>) and (<ref>),  n_ equals the number of disconnected manifolds one obtains by applying all the permutation-inversion operations to M_A(),. 
Why is it reasonable for n_ to appear in the partition function (<ref>)?
We show in the Appendix that  the mapping →δP is an isometry, and therefore 
∫_M_A,  f() dσ() = ∫_M_PA()P^T, δPf() dσ(). 
Therefore, this factor accounts for the integral over the parts of configuration space that we wish to include in the partition function, but that are identical to the factor already computed in the integral over M_A(),. 








We point out that _ provides information on the connectivity of the manifold M_A. In particular, if _⊊× C_2, then M_A is disconnected, with certain components that are not related by reflections.[
If a cluster is chiral, i.e. there is no continuous transformation between  and -, then M_A will be disconnected simply because M_A(),≠ M_A(),-. 
]
 To see why, suppose (P,δ) ∈× C_2 ∖_ is an element in the automorphism-inversion group but not in the sticky symmetry group, and let = δP. 

 Since P ∈,  and  belong to the same manifold M_A. However, since (P,δ) ∉_, there exists no path in M_A connecting  with , and therefore  and  belong to different connected components of M_A. 
We will provide a concrete example where this happens in Section <ref>. 


The converse is not true in general: M_A could be disconnected even if |_|=|| (note we do not say _= because each element of _ is associated with an inversion while the elements in  are not. However, in most cases only one of +P,-P is in _.) We will provide a counterexample in section <ref>. 





  
Relationship to the point group.

It is useful to relate the sticky symmetry group to another group widely studied in physics and chemistry, the point group of a cluster. The point group _ of the cluster  is the set of all the elements of the automorphism-inversion group that can be realized as a rotation of  <cit.>:

    _ = { (P,δ) ∈× C_2 : ∃ R ∈ SO(3): δP = (R ⊗ I_N) } .



There is a useful relation between the point group, the sticky symmetry group, and the automorphism-inversion group of a cluster:

    _⊆_⊆× C_2.

That is, the sticky symmetry group _ is a subgroup of the automorphism-inversion group, and it contains the point group _. 

To see that _⊆_, let (P,δ) ∈_. By the definition of point group  (<ref>), there exists an orthogonal matrix R ∈ SO(3) such that (R ⊗ I_N) = δP. We now show we can achieve this transformation continuously.  Since SO(3) is connected, there exists a continuous path R(t) : [0,1] → SO(3) such that R(0) = I_3 and R(1) = R. We define the path σ(t) = (R(t) ⊗ I_N). For every t ∈ [0,1], σ(t) ∈ M_A since 

    q_ij(σ(t)) = |R(t)_i-R(t)_j|^2 = |R(t)(_i-_j)|^2 = |_i-_j|^2 = q_ij(),

and therefore q_ij(σ(t)) = (r_i+r_j)^2 if A_ij = 1, and q_ij(σ(t))> (r_i+r_j)^2 if A_ij = 0. Moreover, σ(t) is continuous and σ(0) =, σ(1) = (R ⊗ I_N) = δP. Therefore, (P,δ) belongs to the sticky symmetry group _. 


Relationship (<ref>) is useful for computing the sticky symmetry group of clusters that are not too large, because for such clusters the automorphism and point groups may be calculated on reasonable timescales.  
The automorphism group may be computed using algorithms from graph theory <cit.>, and 
the automorphism-inversion group is obtained as a direct product. 
Note that calculating the automorphism group becomes a challenge
 for clusters that are not small, as its size can grow extremely rapidly with N (e.g. see examples in <cit.>.) 

If one can compute the automorphism-inversion group, then one can obtain the 
 point group _ by checking which automorphisms also preserve the set of pairwise distances <cit.>. Let D_ be the N × N matrix whose entries measure the squared distance between the centers of the spheres:

    (D_)_ij = |_i - _j|^2 .


Then  (P,δ) ∈_ for some δ if and only if PD_P^T = D_; in other words, all the pairwise distances are the same after applying the permutation. (Of course, to compute the point group one also needs to determine δ, which can be done straightforwardly). 
Although this result is widely used in physics and chemistry to compute the point group of molecules, we have not found a rigorous proof of it in the literature (we found the standard reference <cit.> to be incomplete), and therefore include one in the Appendix for completeness.



Remark. In physics and chemistry, the counting number n of a rigid molecule with N atoms is often computed as

    n = ξ N!/σ,


where ξ is 2 if the molecule is chiral and 1 otherwise, and σ is the number of rotations which are equivalent to a permutation of the atoms. 
Recall that a cluster  is chiral if there is no rotation-translation operation which maps the cluster to its reflection -. 
This is consistent with our formula (<ref>): if the cluster is chiral, then ξ = 2 and our symmetry number is the same as the one above, σ_ = σ, since no reflections belong to the point group in this case. If the cluster is achiral, then ξ = 1 and our symmetry number is σ_=2σ, since half of the elements of the point group are rotations. 



Remark.
In many applications one wishes to distinguish reflections of a cluster. This framework can be adapted to such a situation, by removing the outer product with C_2 in all the groups under consideration: the sticky symmetry group 𝒯_^0 would be the set of permutations that can be achieved by a continuous deformation, and the counting number would be |𝒢| / |𝒯^0_| = N! / |𝒯^0_|. 










 §.§ Colored particles


The theory developed so far can easily be adapted to colored (partially distinguishable) particles. 
Let  be a cluster with adjacency matrix A, and let V = {1,2,…,N} be labels, one for each sphere. 
We still assume the spheres have identical radii. 
Suppose V is partitioned into k disjoint subsets = { C_i }_i=1^k, where k is the number of colors, so that V = ⋃_i=1^k C_i. 
We denote by ^ the cluster colored according to the partition . 

A permutation P ∈ P(N) acts on the partition as 

    P · C_i = {π(j_1^(i)), …, π(j_p_i^(i)) },      i = 1, …, k,


where π∈ S_N is the permutation associated with P, and  C_i = { j^(i)_1, …, j_p_i^(i)}. 
To count the number of distinct copies of M_A, that we obtain by permuting the particles, we now consider only those permutations which preserve the partition .  
In particular, we define the group 

    P(N)^ = { P ∈ P(N) : P · C_i = C_i,  ∀ i=1, …, k},


which is a subgroup of P(N), and consider the product group P(N)^C × C_2.

The sticky symmetry group _^ of the colored cluster ^ is
    _^ = { (P,δ) ∈_ : P · C_i = C_i,  ∀ i = 1, …, k }.

That is, it is those elements of  the sticky symmetry group of the cluster when the spheres are indistinguishable, _ (see (<ref>)), that preserve the partition. 

  The counting number n_^ for the colored cluster is then given by 

    n_^ = 2|P(N)^|/σ_^,


where σ_^ = |_^|. 

It is a useful fact that once the sticky symmetry group has been computed for indistinguishable particles, with k=1, we may easily obtain the sticky symmetry group and hence symmetry number for any coloring, simply by checking which elements of _ preserve the particular color labels. 





 §.§ Particles with different radii


The framework above may also be adapted to spheres with different radii. 
Specifically,  let (, ) be a cluster in ^3N with radii = (r_1, …, r_N), of which there are n_r distinct radii. 
There is a natural partition ℛ = { R_i }_i=1^n_r of the spheres V={1,2,…,N} obtained by grouping together spheres with the same radius.  
We may then proceed as in the case of colored particles, substituting ℛ for the partition . 
The only significant change occurs when computing the sticky symmetry group of particles with different radii, since the manifold of configurations M_A^() depends on the radii (see (<ref>)). 
Contrary to the case of colored particles, changing the radii  and partition ℛ implies a new computation of the sticky symmetry group _^ℛ. 

Finally, if a cluster has particles with different colors and radii, we first determine the sticky symmetry group _^ℛ according to the partition ℛ, then check which elements of _^ℛ preserve the partition into colors 𝒞 as in (<ref>). We have to make sure the colors respect the radii partition, i.e. two spheres with different radii cannot have the same color. 






§ NUMERICAL ALGORITHM



In this section we provide a numerical algorithm to compute the sticky symmetry group of a cluster , and hence its symmetry number. The algorithm first computes the automorphism and point groups of , and then it checks each automorphism P which is not in the point group to see if either of (P,+1) or (P,-1) are in the sticky symmetry group. To do this, we search for a continuous path in M_A connecting  with P and with -P. If we find a path for either of these cases, we add (P,δ) to the sticky symmetry group. 
Although there is no guarantee that we will find a path if there is one, we show that the algorithm performs very well on several test cases. 

Here is a summary of our algorithm. 



  * We compute the automorphism group  of the cluster (see (<ref>)). 
This can be done using graph-theoretical algorithms such as Nauty <cit.>. 
For very small cases one can simply check each of the N! permutations to see if it preserves adjacency. 

  * We compute the point group _: we calculate the distance matrix D_ of  as in (<ref>), and check which elements in  preserve the distance matrix (see (<ref>)). We then determine δ. In our implementation we found δ using the path algorithm specified in Section <ref> below, simply because it was convenient, but there are simpler methods. 

  * We compute the sticky symmetry group _. 
Since it contains _ (see (<ref>)), we must only check each element in × C_2 ∖_.  For every element (P,δ) in this set, we check if there exists a continuous path in the manifold M_A connecting  and δP, using the algorithm specified in Section <ref> below. If such a path exists, we add (P,δ) to the sticky symmetry group _. 

  * If the particles are colored with partition , we compute the sticky symmetry group _^𝒞 as in (<ref>), by checking which elements of _ preserve the partition. 




We remark that we can use the fact that _ is a group to check for errors in our algorithm. If for some (P,δ) we fail to find a continuous path from  to δP, but one exists, we can sometimes find it afterwards by checking to see if our numerically-computed _ is a group, for example using group theoretical algorithms such as GAP <cit.>. If the elements we have found do not form a group, the algorithm computes the smallest group containing these found elements. 
We found this check particularly useful when developing our path-finding algorithm, where we did sometimes fail to find certain path, but in the path-finding algorithm's current state we haven't found examples where it fails. 
We expect this check to be useful if the error rate is small, as then the check can fill in the rare missing elements. 

We further remark that in computing the sticky symmetry group in step 3  above, it is actually sufficient to find only the generators of the group _. Specifically, one could just find the elements which, together with _, generate (as a group) the sticky symmetry group _, and then generate the entire group using a group theoretical algorithm. We are not aware, though, of a method to test if an exhaustive set of generators of _ have been found. 



The major part of the algorithm is looking for paths in the manifold M_A. We describe a numerical algorithm in the next paragraph. This algorithm may be applied more generally to manifolds defined by equality and inequality constraints, as we will discuss. 





 §.§ Finding continuous deformation paths


We consider a slightly more general setup than that in section <ref>. 
Let M be a set in ^d implicitly defined by 

    M = {∈^d : q_i() = 0,   i = 1, …, m,   h_j() > 0,   j = 1, …, l },

where q_i, h_j : ^d → are smooth functions. We assume the gradients {∇ q_i() }_i=1^m are linearly independent for every ∈ M, which implies that M is a differentiable manifold of dimension p= d-m. Let _0,_1 ∈ M be two distinct points in M. Our aim is to determine if _0 and _1 belong to the same connected component of M, i.e. if there exists a continuous path φ : [0,1] → M such that φ(0) = _0 and φ(1) = _1. 

We start with the following observation. Let's forget about being on a manifold, and suppose we want to minimize in ^d the function 

    U() = | - _1|^2

using the method of steepest descent starting from point _0 <cit.>. We do this by computing a sequence of points _0,  _1,… with _0 = _0 and 

    _k+1 = _k - t_k ∇ U(_k),    k = 0, 1, … .

The best choice of  t_k is the one that minimizes the function 

    ϕ_k(t) = U(_k-t∇ U(_k)) = (1-2t)^2|_k-_1|^2.

The minimum occurs at t_k=1/2, 
which implies that the second point in the sequence is 
_1 = _0 - 1/2· 2 (_0 -_1) = _1. 
The steepest descent in this case is trivial, leading us in one step to the true minimum. 

Suppose now that we want the sequence (<ref>) to approximate a continuous path and to be constrained to the manifold M. To approximate a continuous path we may put an upper bound on the step size, say . To be constrained to the manifold, we project the steepest descent direction _1-_k to the tangent space to the manifold, take a step in this direction, and then project back to the manifold. We choose a projection such that the correction step is perpendicular to the tangent space at _k, so that if there were no limit on the step size, the complete step including tangent step plus projection step would take us directly to _1.  

Specifically,  suppose _k ∈ M is the point generated at step k, and let T_kM denote the tangent space of M at _k and T_k^⊥ M its complement, the normal space. Let P_k : ^d → T_kM be the matrix which orthogonally projects a vector to the tangent space. 
We generate a sequence of points as 

    _k+1 = _k + Δ s _k  + _k,      k = 0,1, …,


where Δ s >0 is the step-size, _k∈ T_kM is a unit vector in the tangent space at y_k, and 𝐰_k ∈ T_k^⊥ M is a vector in the normal space at y_k. 
The direction in the tangent space is chosen by projecting the steepest descent direction as 

    _k = P_k(_1-_k)/|P_k(_1-_k)|.

The vector _k in the normal space is chosen so that _k+1∈ M. 
Specifically, since {∇ q_ij(_k)}_(i,j)∈ E spans T^⊥_k, we let 
= ∑_j=1^m a_j ∇ q_j(_k) for some unknown coefficients {a_j}_j=1^m, and then 
solve the system of equations q_i(+∑_j=1^m a_j ∇ q_j(_k)) = 0, i=1, …, m using Newton's method. Generically the solutions for _k will be isolated. 

The step size is chosen as

    Δ s = min{,  |P_k(_1-_k)|}.

This choice is motivated by the following observation: if Δ s = |P_k(_1-_k)|, then we would find a solution for the normal step is _k = P_k^⊥ (_1-_k), where P_k^⊥ is the orthogonal projection matrix to T_k^⊥ M, so _k+1 = _1: the algorithm would bring us to the optimal point in one step. 
The step in the optimal direction is simply broken up into a step in the tangent space, and a step in the normal space. We impose an upper bound  on the step size to approximate a continuous path, but take a smaller step if it would be better. 


If ever  |_k-_1|< the algorithm would take us to _1 in one step, so we stop generating points, and say we have found a path from _0 to _1. 

When the manifold has a boundary ∂ M, because there are one or more inequalities in (<ref>), the sequence _k may fall outside of the boundary. Therefore, we check at each step if _k satisfies the inequalities, i.e. if h_j(_k) > 0, for every j = 1, …, l. If not, we stop the steepest descent and switch to a different, random method to generate points. 

We switch to random steps if any of the following occur: 


  * A point _k lies outside the boundary. 

  * The projection onto the manifold fails, i.e. we fail to find a solution _k to the nonlinear equations. 

  * The point _k is stuck in a local mininum, i.e. |∇ q_ij·_k| <, for some i,j = 1, …, N, and  is a numerical tolerance. 



If we decide to switch to random steps, we generate some number  of random steps, and then resume the steepest descent. 
We tried two different methods for generating random points. 



  * Random steps. Given a point _k ∈ M, we generate a random vector ∈ T_kM according to the isotropic Gaussian density

    p() = 1/(2π)^p/2σ^pexp(-||^2/2σ^2),
 

where σ >0 is a parameter and p is the dimension of the manifold M. We then project the point  _k+ back to the manifold using the same method as before and generate the point
_k+1 =  _k +  + _k.
In our implementation, we typically chose σ comparable to .  This could occasionally produce steps that are larger than , so for a strict tolerance on the step size one might wish to truncate the larger steps. 



  * Sampling. We sample points from the density

    ρ() = 1/Zexp(-β U()),     where   
    U() = {[ |-_1|^2   x ∈ M;       ∞   x ∉ M ].


where β∈ is a real parameter and Z>0 is a normalization constant. We sample on the manifold from the density ρ using the Markov Chain Monte Carlo algorithm described in <cit.> and a step size parameter σ, again usually comparable to . The parameter β is called the “inverse temperature" in simulated annealing or other sampling techniques <cit.>, and it controls how peaked the density is near the minimum of U(): large β means the density is strongly peaked near the minimum so sampling pushes points toward _1, β≈ 0  means the density is relatively flat so the manifold is sampled nearly uniformly (similar to the first method but with a Metropolis step to ensure we sample the correct density), and β <0 means the density is lowest at the optimum, so sampling should push us away from _1 in general. We found that β <0  helped to overcome boundary obstacles, as we describe later in our numerical experiments. 



We remark that since steepest descent is like sampling with β→∞, our method of switching between steepest descent and random sampling is basically a form of simulated annealing with temperature cycling. 

Note that another option would be to use an active set method <cit.>, and add in additional equations when one hits a boundary to navigate along the boundary while still performing steepest descent. We leave this option for future work. One may still have to resort to a random method for an arbitrary manifold with an arbitrary boundary. 



Our algorithm terminates when one of two conditions is met. One, if ever |_k-_1|<, we stop, and declare that we have found a path. Two, if ever the total number of points exceeds some maximum , we stop, and declare there is no path. 



To summarize, our path-finding algorithm consists of the following steps: 


  * Generate points according to steepest-descent algorithm, as in (<ref>). 

  * If one of conditions (a-c) is met, switch to random mode, and generate  random steps as in (<ref>) or (<ref>). 

  * Stop if either (a) |_k-_1|< (declare path found), or total number of points exceeds   (declare no path.) 




The algorithm could fail, either by failing to find a path that exists (false negative), or by finding a path that doesn't actually exist (false positive.) More commonly is to fail to find a path when one exists.  Similar to the problem of finding low-energy paths between points on an energy landscape <cit.>, no known algorithm can guarantee to find paths between points on arbitrary nonlinear manifolds, with arbitrarily complicated boundaries to navigate; certainly the longer one looks, the more likely one is to find a path, but we don't expect to ever be able to provide any guarantees. 
It is helpful that we can check a posteriori if _ is a group, as then we can fill in missing elements if the number of such elements is small. 

The algorithm can also find a path that doesn't exist in the continuum problem, for example if the numerical path jumps between disconnected components of the manifold, a possibility if the components approach closer than the continuity parameter . We never found an example where this happened, but it is possible with inappropriate parameter values or a particularly complicated manifold geometry. 










  
A toy example in ^2. 

We apply our algorithm to a toy example to illustrate how it works. We consider the set 

    D = { (x,y) ∈^2 : |x|<3,   0<y<x^2+5 }.


This set is visualized in Figure <ref>. We fix the points = (-2.9, 12.5) and = (2.9, 12.5), and apply our algorithm to find a path in D connecting  with . The path has to navigate around a boundary, moving in the opposite direction to the steepest descent direction - for quite some time, to reach . 
We tested both methods for generating random points, and the resulting paths are shown in Figure <ref>. Notice that the sampling with negative parameter β results in a path with significantly fewer step points than the method with random steps. This is because the points are “pushed away" from , and hence from the boundary ∂ M, resulting in the sequence hitting the boundary fewer times. 






§ EXAMPLES


In this section we apply the theory developed so far to several examples, both to show that it works  and to analyze how well, as well as to point out pedagogical examples that illustrate properties of the sticky symmetry group. In the following, unless specified, we assume that the particles are spheres with diameter d = 1. 
In addition, we remove the three translational degrees of freedom from each cluster by fixing the center of mass to the origin, adding the additional three constraints ∑_i=1^N _i = 0 to the manifold (<ref>). This doesn't change any of the theory regarding the symmetry number or how to compute it. 
In our implementation we usually set = 10^4-10^5.





  
Loops and chains. We consider loops (L) and chains (C) of N identical particles (see Figure <ref>). A loop of N particles has N contacts, while a chain has N-1 contacts. These are flexible clusters. In particular, if we denote by M_L,N and M_C,N the manifold of configurations of a loop and a chain of N particles, respectively, we have M_L,N = 2N-3 and M_C,N = 2N-2. 

We compute the automorphism group, the point group and the sticky symmetry group of loops and chains for  N=4-10,15,20. In order to test that our algorithm works, we randomly choose a point  in the manifold of configurations using the sampling algorithm in <cit.>, and compute _ and _. In this way,  will have no “a priori" symmetry axes.  For a loop of N particles, the automorphism group _L,N is isomorphic to the planar dihedral group D_N, the symmetry group of a regular N-gon in the plane. The point group of a random configuration _L,N is trivial (it consists of only the identity (E, 1)).
Using our algorithm, with σ= = 0.1, and = 20 random steps, we find that _L,N = _L,N× C_2, i.e. it is the entire automorphism-inversion group. One can visualize why this must be: starting from a configuration ∈ M_L,N, the loop can deform continuously until it reaches its most symmetrical configuration, a polygon that lies in a plane, and then it rotates according to the permutation P ∈_L,N in the automorphism group, and finally it deforms again until it reaches the new configuration ±P. The symmetry number of a loop of N particles is then σ_L,N = |_L,N| = 4N, and the counting number is  n_L,N = 2N!/(4N) = (N-1)!/2.  


For a chain of N particles, the automorphism group _C,N consists of two elements: the identity E and the permutation 

    π = {[ (1   N)(2   N-1) (3 N-2) …( N/2 N/2+1 )                            if N is even;  (1   N)(2   N-1) (3   N-2) …( N-1/2+2)                             if N is odd ].


which corresponds to a reflection along the axis passing through the center of mass of the chain, when the chain assumes its most symmetrical configuration, with the spheres lying in a line. The point group of a random configuration is trivial, and again applying our algorithm we find that _C,N = _C,N× C_2. Therefore, the symmetry number is σ_C,N = 4 and the counting number is n_C,N = 2N!/4 = N!/2.  






In Figure <ref> we plot the average number of steps to find a path for loops and chains of different sizes, where the average is over all paths that were successfully found, i.e. all elements in the corresponding sticky symmetry groups (we didn't consider variations with the initial condition or with different realizations of the noise.) The average number of steps increases roughly linearly with N, though possibly faster than linearly for chains. 

 We point out that, in the context of molecular symmetry, the loop of N = 6 particles corresponds to the benzene molecule.  Our algorithm finds the symmetry number σ_L,6 = 24 for the 6 loop, in agreement with previous results <cit.>, where the symmetry number of benzene is computed using GAP[Actually, the symmetry number for the benzene in <cit.> was 12, since the authors only consider permutations, and not permutation-inversion operations.].








  
Cluster of particles with different radii. 
We consider again a loop of N = 6 particles, but this time we set the radius of the particles 2,4 and 6 to be r = 0.6, and the radius of particles 1,3 and 5 to be r' = 0.4 (see Figure <ref>). In this case, we first of all check which element of the automorphism group _L,6 of the loop (with identical particles) preserves the partition R of the vertices according to the radii (compare with (<ref>)). We find that the automorphism group ^ℛ of the cluster is  

    ^ℛ = { (135)(246), (153)(264), (26)(35), (13)(46), (24)(15), E }.


In particular, ^ℛ is generated by the rotation (135)(246) and the reflection (13)(46), and it is isomorphic to the dihedral group D_3, the symmetry group of an equilateral triangle.  
Using our numerical algorithm (applied to the manifold of configurations with different radii as in (<ref>)), we find that the sticky symmetry group _^ℛ corresponds to ^ℛ× C_2. 





  
Cluster of N = 6 spheres with two bonds broken. Consider clusters with N = 6 identical spheres. There are two rigid clusters with m = 12 contacts, the octahedron and the polytetrahedron <cit.>. We consider all the clusters with m=10 bonds, formed by deleting two bonds from a rigid cluster (and keeping only clusters with nonisomorphic adjacency matrices.)
The dimension of each manifold is 18-10-3 = 5. 
For each cluster we compute its symmetry and counting number. The results are given in Table <ref>, where each cluster is referred to as a “mode",  adopting terminology from <cit.>. In the computations, we used random steps to explore the boundary, and parameters σ =  = 0.1. These results agree with the computations done in <cit.> for the same clusters, where the symmetry number of flexible clusters is computed using a combinatorial argument[The counting number in <cit.> is actually half the counting number reported here, since the reflections were not taken into consideration.]. 







  
A cluster of N = 6 spheres with _⊊_⊊× C_2. 

We consider in detail a cluster  with N = 6 spheres and m=10 bonds among the ones analyzed above (mode 14 in <cit.>.) It is plotted in Figure <ref>, where we show our labeling convention of the spheres. 
This is an example of a cluster where _⊊_⊊× C_2, and the inclusions are strict. 

Indeed, 
the automorphism group  consists of 16 elements (we write them as permutations of the labeled particles):

    =    { E, (14)(23)(56), (14)(23), (1234)(56), (1234), (13)(24)(56), (13)(24), (13)(56), 
        (13), (1432)(56), (1432), (12)(34)(56), (12)(34), (24)(56), (24), (56) },


where E denotes the identity element. 

We compute the point group _ of the cluster by checking which of the permutations in  preserve the distance matrix D_ (compare with (<ref>) and (<ref>)), and additionally determining the associated inversion. This consists of four permutation-inversions: 

    _ = {E, (12)(34)^*, (56)^*, (12)(34)(56)},


where the * indicates that the permutation is combined with a reflection (δ = -1). 
This group is generated by one reflection, and one rotation along the axis passing through the particles 5 and 6 and the plane which contains the particles 1,2,3,4. 

We compute the sticky symmetry group by checking which elements in × C_2 ∖_ belong to the sticky symmetry group. For any P ∈× C_2∖_, we use our numerical procedure to check if there is a continuous path in the manifold of configurations M_A connecting  with ±P, where A is the adjacency matrix of . We set = σ = 0.1, and we adopt random sampling with β = -0.1. We find the sticky symmetry group _ consists of eight elements

    _ ={E, (12)(34)^*, (56)^*, (12)(34)(56), (14)(23)(56)^*, (14)(23), (13)(24)(56), (13)(24)^*}.
 


Therefore _⊊_⊊× C_2. Because there exists an element P ∈ such that (P,δ) ∉_, for any δ∈ C_2, the manifold M_A is not connected and the disconnected components are not related by a reflection. For example, the permutation (24) which switches particles 2 and 4 preserves adjacency, however it is not realizable as a continuous deformation or rotation or reflection of the cluster. 



  
A colored cluster.
Next, we consider the same example as above but distinguish the particles using three different colors (Figure <ref>, right). We form the partition ={C_1,C_2,C_3} with  C_1 = {1,2}, C_2 = { 3,4} and C_3 = {5,6}. We compute the sticky symmetry group _^ of the colored cluster by checking which elements in _ preserve the partition = {C_i : i =1,2,3} (compare with (<ref>) and (<ref>).) We have that _^ coincides with the point group _ of the original cluster , since all the elements of _ preserve the partition into colors. 

We point out that, in the case of colored particles, we only need to compute the sticky symmetry group _ for identical particles once, then check which elements of it preserve the coloring. In other words, we can change the partition without having to check again the connectivity of the original manifold of configurations M_A. This is not true if we change the radii, because then the manifolds themselves will change, so the paths need to be recomputed. 










  
Two clusters of N = 11 spheres where |_|= || but M_A is disconnected. 

We consider in detail two rigid clusters _1 and _2  (discovered in <cit.>) with N = 11 particles, which have the same adjacency matrix A, however there is no continuous transformation from one to the other or its reflection (see Figure <ref>). This implies that they belong to different connected components of the manifold of configurations M_A, i.e. M_A,_1≠ M_A, _2.  

We use these clusters to build two flexible clusters with the same adjacency matrix, such that the sticky symmetry group has the same size as the automorphism group[
Note that we don't say they are the same because the sticky symmetry group includes inversion operations whereas the automorphism group does not; see (<ref>).
], but that live on disconnected components of M_A. To do this we break one bond (bond 6-11) in each cluster _1 and _2. The resulting two flexible clusters, which we denote by _1 and _2, have the same adjacency matrix A' so they belong to the same manifold M_A'. However, we do not find a continuous path in M_A' connecting _1 with either ±_2, which implies that they belong to two different connected components of M_A' that are not related by a reflection. The automorphism group  of _1 and _2 is given by 

    = {E, (23)(45)(78)(9   10)}.


Using our algorithm, we find that __1=__2, and every element in  induces a single permutation-inversion operation in the sticky symmetry groups __1,__2. 

Therefore, the manifold M_A' is not connected, even if |__1|=|__2|=.   
Note that, in contrast to the previous example, we expect the two distinct components to be nonisomorphic, since clusters on each component are not related even by a permutation-inversion operation; we verified this heuristically by building physical models of the clusters with balls and magnetic sticks. 




  
Partition functions of every connected cluster of N=6 identical spheres 







As a final example we compute the partition function (<ref>) of every connected cluster of N=6 indistinguishable sticky spheres, at every temperature. With these partition functions in hand, one can determine essentially any equilibrium quantity one wants at any temperature for a system of 6 spheres (conditional on the spheres being connected), and, with further similar computations, one could begin to ask about non-connected systems, non-identical spheres, and spheres with non-identical interactions. Performing such an exhaustive calculation of the free energy landscape for sticky particles is an accomplishment in itself. Such an exhaustive calculation is not possible for non-sticky particles, i.e. those with smoother interaction potentials, where free energy landscapes are typically characterized by a collection of local minima and saddle points <cit.>,  which vary with temperature and cannot be guaranteed to be found by existing numerical algorithms. 

We start with the two rigid clusters, which are known to be the only clusters with m=12 contacts <cit.>. On each cluster, we break, in turn, each bond, each pair of bonds, each triple of bonds, etc. For each graph α so obtained we check that it is connected, and that it is not isomorphic to a graph we have already seen.[For larger systems one would also need to check that each isomorphic graph is also on the same connected component of M_A, (after applying the permutation), but we don't anticipate disconnected copies to be a problem in such a small system.] If it passes both tests, then we compute the integral in (<ref>) using the method described in <cit.> (with center of mass of the cluster fixed), which is essentially a form of thermodynamic integration on a manifold, and call the resulting quantity I_α. Then we calculate the symmetry number σ_α using the method described in this paper. The so-called “geometrical” partition function for graph α (i.e. the partition function without the factors of κ) is, up to a constant that is the same for all clusters under consideration, 
z^g_α = I_α/σ_α. 

The table in Figure <ref> reports the total number of clusters found for each dimension d, where d is, equivalently, the number of bonds broken from a rigid cluster, the number of internal degrees of freedom of the cluster, or the dimension of the corresponding manifold M_A, in (<ref>) (minus 3, to account for rotational degrees of freedom.) 
The number of clusters increases up to d=4, and then decreases to d=7, for a total of 99 clusters; for d>7 there are no  connected clusters. 
This table also reports the total geometrical partition function Z_d for clusters with dimension d, formed by summing the geometrical partition functions z^g_α over all graphs α with dimension d. 
The ratios Z_1/Z_0 = 7.1, and Z_2/Z_1 = 4.6, agree within numerical error with those reported in <cit.> (7.1, 4.5 respectively), partially verifying that our computations are correct; it is worthwhile to note that the methods used to compute Z_0,Z_1,Z_2 in <cit.> cannot be extended to higher-dimensional clusters. 

Figure <ref> also shows the histograms of symmetry numbers for clusters of each dimension. The most symmetric cluster, at d=7, has a symmetry number σ=120. This is the cluster formed when 5 spheres touch a central sphere with a single contact each. The elements in the symmetry group are all of the 5! permutations of the outer spheres. 

The actual partition function for any cluster with graph α is obtained from the geometrical partition function as z_α = κ^m z^g_α, where m=3N-d is the number of contacts in the cluster, and κ is the sticky parameter, which measures how strong the bonds are between spheres: large κ means spheres like to spend more time sticking together; κ small means they come apart easily. For an experimental system where the pair interactions are not perfect delta functions, κ would be a function of temperature, and the width and depth of the actual attractive interaction potential between the particles <cit.>. It is large when temperature is small, and/or when the pair potential is deep and wide. Experimentally, clusters tend to form and rearrange when κ≈ O(10)-O(100); any smaller, and the system evaporates, and any larger, and it takes a very long time for clusters to rearrange <cit.>. 

From the partition functions one can calculate the equilibrium probability to find the system in any given configuration. For example, 
Figure <ref> (left) shows the equilibrium probabilities of finding the system in a cluster of a given dimension, at different values of κ (recall this probability is conditional on the cluster being connected), equal to  κ^18-dZ_d / Z, where Z=∑_d=0^7 κ^18-dZ_d is the total partition function for the system. For large κ, the system is most likely to be in the lowest-dimensional configurations (those with the most contacts), while for small κ, it is likely to be in the highest-dimensional ones. For κ≈ 1-7, the probability has a maximum at an intermediate value of d; a range of κ that coincides with the range of ratios Z_d+1/Z_d (from smallest to largest d: 7.1, 4.6, 4.5, 3.8, 2.2, 1.6, 1.0.) 

Figure <ref> (right) shows the probability of finding the system in a floppy configuration (d>0) as a function of κ. This probability is very close to 1 for small κ, and very close to 0 for large κ, however it has a notably wide transition region, κ≈ 4-70, where the probability is not close to either endpoint (between 0.1-0.9); it crosses 0.5 at κ≈ 11.5. In this transition region one would expect interesting dynamics, with the system forming clusters but rearranging substantially on observable timescales; hence, this is the range of κ values one should aim for experimentally. 





§ CONCLUSION


We developed a theoretical and computational framework to compute the symmetry number of flexible sticky-sphere clusters, i.e. hard spheres interacting with a delta-function attractive interaction potential. 
We started from an equivalence relation which says that two clusters are the same if they are related by any combination of rotations, translations, or deformations that don't change the sphere-sphere contacts, 
 transformations which are all continuous and preserve the energy of a sticky-sphere cluster. 
 We showed how to count the number of distinct equivalence classes that one obtains by considering all permutations and reflections of a given cluster, a number we called the counting number. The counting number is related to a cluster's symmetry number, which is in turn obtained from its sticky symmetry group. 
We analyzed the sticky symmetry group and showed how it is related to two other groups commonly used to study molecular symmetries, the point group and the automorphism group of the graph describing the cluster's pairwise contacts. 

We introduced a numerical algorithm to compute the sticky symmetry group of a cluster. The key part of the algorithm is finding continuous paths in the manifold of configurations of the cluster. For this, we provided a numerical procedure alternating between steepest descent and random sampling, a form of simulated annealing. The algorithm comes with no guarantees, and could produce both false positives or false negatives, however we found with the right choice of parameters it worked extremely well for small clusters. In addition, we can use the fact that the collection of paths forms a group, to catch rare false negatives, i.e. paths that aren't found by the algorithm. An interesting question for future research, would be for what false negative rate is it possible to compute the entire group with high probability. 

We tested our algorithm on small clusters, some with up to N=20 spheres. 
Our algorithm is efficient when it is possible to compute the automorphism group , and gives the symmetry number many times faster than any hand computation would. 
For larger clusters or those with very high symmetry, several problems can occur. First, the automorphism group  can be prohibitively large to compute; for example <cit.> gives an example of a cluster of 33 particles whose automorphism group has more than 10^35 elements. Nevertheless, it may be possible to overcome this challenge by adapting algorithms for rigid molecules that cleverly avoid  computing the entire automorphism group <cit.>. Secondly, as a cluster becomes larger, the geometry of its manifold can become more complicated, so it can require longer to find paths along it. 
Finally, for large clusters, there are a large number of constraints defining the manifold (<ref>), so computing the tangent space, and projecting back to the manifold after taking a step in the tangent space, are more computationally expensive. 


Our method can be used as an essential component of an algorithm that computes the entire free energy landscape of small clusters. As a step in this direction, we computed the partition functions for every possible connected cluster of N=6 identical spheres. One could perform a similar calculation for any N≤ 8, and for spheres with different radii, however for identical spheres larger N would require dealing with contact constraints that are not linearly independent, leading to singularities in the manifolds of configurations that are challenging (though not impossible) to deal with numerically <cit.>. Of course, for large enough N one cannot exhaustively enumerate every configuration, but nevertheless symmetry or related topological considerations can be important in detecting different regimes of behaviour in systems of hard particles <cit.>.



With such data, one can then ask how the free energy landscape varies as one varies the pairwise interactions between spheres. 
This question is important in materials science, where one might want to design particles, such as colloids, to assemble into a particular cluster, by varying the strengths and specificity of the interactions <cit.>, control that is possible by coating colloids with strands of sticky DNA <cit.>.  In the sticky-sphere limit, such variation in interactions would correspond to changing the colorings of the particles and the strengths of the interactions between colors, represented by the coloring partition  and the sticky parameters κ_ij where i,j are colors (see (<ref>), and discussion thereafter.) 
Our method leads to an efficient way to solve this problem, since 
once we have calculated the integrals in (<ref>) and the sticky symmetry groups _ for each connected component M_A, of the landscape, we may 
obtain the partition functions for all clusters, for any coloring of the particles and any interaction strengths between colors, without repeating these arduous computations. 
Therefore, we expect the framework outlined in this paper to be an important component of an algorithm which aims to efficiently compute cluster probabilities over a wide range of interactions strengths and specificities, and to solve the inverse problem of asking how to make a particular cluster with high probability, given experimental constraints. 



    



  
Acknowledgments. 
We thank Louis Theran for useful discussions.  
This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research under award DE-SC0012296.
M.H.-C. acknowledges support from the Alfred P. Sloan Foundation.




§ APPENDIX: PROOFS


In this Appendix we provide the proofs of some statements of Section <ref>. 

First of all, while it seems like an obvious fact, we will need to know that if we apply a permutation-inversion (P,δ) to a cluster , we obtain a manifold M_PAP^T,δP which is an isometry of M_A,, i.e. it is geometrically indistinguishable from it. We prove this statement in the following lemma, which in addition will imply that continuous paths in M_A, get mapped to continuous paths in M_PAP^T,δP after applying the permutation-inversion operation (P,δ) to each element along the path. 


Let ∈^3N be a cluster of N identical spheres. Let ' = δP, with P ∈ P(N) and δ∈ C_2. Let A' denote the adjacency matrix of '. The map 

    ϕ_P, δ : M_A,    ⟶ M_A','
       ⟼δP

is an isometry.
 


The function ϕ_P,δ is clearly smooth and bijective, with inverse ϕ_P, δ^-1 = ϕ_P^-1, δ^-1. Let T_M_A and T_'M_A' denote the tangent spaces at  and ' of the manifolds M_A and M_A', respectively. 
Because ϕ_P,δ is a linear function, 
the tangent map Tϕ_P,δ : T_M_A ⟶ T_'M_A' may be easily calculated to be Tϕ_P,δ() = δP. 
We have, for every , ∈  T_M_A, 

    ⟨ Tϕ_P,δ(), Tϕ_P,δ()⟩_' = ⟨δP, δP⟩ = δ^2 ^T P^T P = ^T P^TP = ⟨, ⟩_,


where ⟨, ⟩_ (respectively ') is the standard Euclidean metric in ^3N restricted to in M_A (respectively M_A'.) This proves our claim. 




This lemma allows us to prove that the function given in (<ref>) is in fact a well-defined action of P(N) × C_2 on the quotient set X given in (<ref>). 


The function · : (P(N) × C_2) × X ⟶ X given by 

    (P,δ) · [] = [δP] = M_PA()P^T, δP


is a well-defined action of P(N) × C_2 on the quotient set X.
 

We need to show two things: one, that the action is well-defined on equivalence classes, and two, that it respects the property of multiplication. 

For the first, 
suppose ∼ are two equivalent clusters. We want to show that δP∼δP or, equivalently, that [δP] = [δP], for each (P,δ) ∈ P(N)× C_2. We have 

    [δP] = M_PA()P^T, δP,      [δP] = M_PA()P^T, δP.


Since ∼, we have that M_A(),  = M_A(),. In particular, this implies A() = A(), and the existence of a continuous path φ : [0,1] ⟶ M_A() such that φ(0) = and φ(1) =. We now consider the isometry ϕ_P,δ as in (<ref>) and construct the path 

    ψ(t) = (ϕ_P,δ∘φ)(t) : [0,1] ⟶ M_PA()P^T ,


which is a continuous path in M_PA()P^T connecting ψ(0) = δP with ψ(1) = δP. Therefore, δP and δP belong to the same connected component of M_PA()P^T, implying M_PA()P^T, δP = M_PA()P^T, δP, which proves our claim. 

Finally, · is an action since, for every (P,δ), (Q,μ) ∈ P(N)× C_2,

    (P,δ) ·((Q,μ) · []) = (P,δ) · [μQ] = [δμPQ] = (PQ, δμ) · [] = ((P,δ) (Q, μ)) · [],


and clearly (I_N,1) · [] = [], where I_N is the N × N identity matrix. 


Finally, we prove the connection between point group and distance matrix of a cluster, a result that is widely used but for which we have found no accessible proof in the literature, so we prove it here for completeness.  
A canonical (though incomplete) reference for this result is  <cit.>, and <cit.> provides a clear explanation of distance matrices and some common manipulations with them. 


Let ∈^3N be a cluster of N spheres, and let _ be the point group of  as in (<ref>). Then 

    (P,δ) ∈_ for some δ∈ C_2   ⟺   PD_ = D_P.





Let 

B = (_1| … |_N)^T


 be the N × 3 matrix whose rows are _i, for i = 1, …, N. Let G_ be the N × N Gram matrix of , computed from B as
 
    G_ = BB^T .
 
 We first prove that preserving the distance matrix is equivalent to preserving the Gram matrix, i.e. 

    PD_ P^T = D_  ⟺   PG_P^T = G_,


for every P ∈. 
To see this, first notice that the matrices G_ and D_ are related by the formula 

    D_ = diag(G_)1^T -2G_+1diag(G_)^T,

where 1 is a N × 1 column vector of all ones <cit.>. 
Next, notice that P(diag(G_)1^T)P^T = diag(G_)1^T. In fact, the i-th row of the matrix K = diag(G_)1^T is a 1 × N vector with the same entries, i.e. of the form (c_i, …, c_i), for some c_i ∈. Therefore, the permutation matrix P permutes the rows of K, which are then permuted again using the inverse permutation P^T, resulting in the original matrix K.  
We use this and (<ref>) to compute

    PD_P^T 
       = diag(G_)1^T-2(PG_P^T) + 1diag(G_)^T .

Subtracting (<ref>) from (<ref>) gives

    PD_P^T - D_ = -2(PG_P^T -G_ )

from which implication (<ref>) is clear. 



We now shift the attention to the Gram matrix G_. Specifically, we want to prove that, given P ∈, then (P,δ) ∈ P_ for some δ∈ C_2, if and only if PG_P^T = G_. 



Showing that a permutation in the point group preserves the Gram matrix follows by direct calculation. 
Suppose (P,δ) ∈ P_. Then, by definition of _ (see (<ref>)), there exists R ∈ SO(3) such that P = δ (R ⊗ I_N). Let = P, and C = (_1| … |_N)^T. Then C = PB. The Gram matrix G_ of  is 
G_ = CC^T = (PB)(PB)^T = PG_P^T. 
On the other hand, since = δ(R ⊗ I_N), then C = δ BR. This implies 
G_ = (δ BR)(δ BR)^T = BB^T = G_. 

Therefore PG_P^T = G_.


For the other direction, suppose P is such that PG_P^T = G_. Then  (PB)(PB)^T = BB^T,  so the clusters  (formed from the rows of C=PB) and  (formed from the rows of B) have the same Gram matrix. 

It remains to show that if two clusters have the same Gram matrix, then they are related by an orthonormal transformation. 
This is a result in linear algebra that we reproduce here. 
Let 

    B = U_1Σ_1 V_1^T,      C = U_2Σ_2 V_2^T

be the singular value decompositions of B,C∈^N× 3, 
where U_i∈^N× N, V_i∈^3× 3, Σ_i∈^N× 3, for i=1,2. 
We are given that G_=BB^T=CC^T and so U_1Σ_1Σ_1^TU_1^T = U_2Σ_2Σ_2^TU_2^T. But this is an eigenvalue decomposition of the symmetric matrix G_, which is unique up to reordering of eigenvalues and up to eigenvalues that are the same. We may order the eigenvalues in order of decreasing absolute value, and therefore the diagonal elements of the diagonal matrices Σ_1Σ_1^T, Σ_2Σ_2^T may be chosen to be the same, which implies, since the diagonal elements of Σ_1,Σ_2 are nonnegative, that Σ_1=Σ_2. 
For eigenvalues that are the same, we may choose any orthogonal basis for the corresponding columns of U_i among the available eigenvectors, and therefore we may choose the bases such that U_1=U_2. Therefore, we may write C=U_2Σ_2V_2^T=U_1Σ_1 V_2^T = U_1Σ_1V_1^TV_1 V_2^T =BV_1V_2^T = BQ where Q=V_1V_2^T∈^3× 3∈ O(3), and the result is proven. 








unsrt



