Identiﬁcation of Local Structure in 2-D and 3-D Atomic Systems through Crystallographic Analysis

: In the present work, we revise and extend the Characteristic Crystallographic Element (CCE) norm, an algorithm used to simultaneously detect radial and orientational similarity of computer-generated structures with respect to speciﬁc reference crystals and local symmetries. Based on the identiﬁcation of point group symmetry elements, the CCE descriptor is able to gauge local structure with high precision and ﬁnely distinguish between competing morphologies. As test cases we use computer-generated monomeric and polymer systems of spherical particles interacting with the hard-sphere and square-well attractive potentials. We demonstrate that the CCE norm is able to detect and di ﬀ erentiate, between others, among: hexagonal close packed (HCP), face centered cubic (FCC), hexagonal (HEX) and body centered cubic (BCC) crystals as well as non-crystallographic ﬁvefold (FIV) local symmetry in bulk 3-D systems; triangular (TRI), square (SQU) and honeycomb (HON) crystals, as well as pentagonal (PEN) local symmetry in thin ﬁlms of one-layer thickness (2-D systems). The descriptor is general and can be applied to identify the symmetry elements of any point group for arbitrary atomic or particulate system in two or three dimensions, in the bulk or under conﬁnement.


Introduction
In recent decades, computer simulations at the molecular level [1,2] have steadily proved to be an invaluable tool in understanding the connection between atomic structure and behavior and macroscopic properties of materials [3,4]. Significant advances have been made in numerical methods and algorithms enhancing simulation aspects such as speed, accuracy, general applicability, user friendliness and robustness. In parallel, advanced software packages are now available that allow efficient atomic/molecular visualization and/or analysis of the computer-generated structures including animation of the corresponding trajectories [5][6][7][8][9]. For crystallization and phase transitions, visual inspection is vital as it allows a better understanding of structural changes and provides a preliminary identification of the established ordered morphologies. However, visual inspection has two disadvantages: possible observer's bias and lack of systematic quantification. It is thus critical to employ objective and quantitative methods in the analysis of local structure in computer-generated atomic and particulate systems.
The simplest possible form of information on structure can be extracted through the pair radial distribution function, g(r), which calculates the normalized probability of finding a pair of atoms lying at a distance r [10]. In computer simulations the calculation of g(r) is of particular importance as it is related, through a Fourier transform, to the static structure factor S(k) which can be measured experimentally through X-ray scattering [1]. The g(r) data, especially of the first two peaks, can be further used as input for more refined analysis, such as the calculation of normalized abundancies of pairs, as in the work of Honeycutt and Andersen [11]. A simplified variant has been presented by Ganesh and Widom [12] to gauge the structure of liquid and supercooled liquid copper. Common neighbor analysis (CNA) [13] proceeds by decomposing the radial distribution function based on the local environment of pairs. While simple in implementation, it allows for distinction between different crystals in particular of face centered cubic (FCC), hexagonal close packed (HCP), body centered cubic (BCC) and icosahedral symmetry. Very recently, the pair-angle distribution function (PADF) has been introduced with the corresponding three and four body real space distributions constituting a generalization of the pair function [14]. The bond orientation order parameter, introduced in [15], is based on the quantum mechanical concept of rotationally invariant combinations of spherical harmonics and is a widely used tool to study orientational order in general liquids. Polyhedral template matching (PTM) has been proposed in [16], which classifies atomic structure according to the topology of the local environment. In parallel, geometric approaches focus on extracting the corresponding structural information of the local environment through analysis of the Voronoi polyhedron (cell) [17,18], first employed by Finney [19] to study the random close packing (RCP) of simple liquids. Further examples include the topological analysis of Voronoi faces [20] or superposition of simplexes and calculation of the corresponding Procrustean distance [21]. A heuristic approach to gauge local structure is also available based on angular and radial histograms as explained in [22].
Different local structural parameters have been developed over the years. Some examples include the free-volume approach [23], the definition of internal stresses and the symmetry coefficients at the atomic level [24], the concept of flexibility volume [25], the topological cluster classification algorithm [26], the neighbor-distance analysis [27] and the community inference [28].
In an effort to gauge local environment and identify possible phase transitions in computer simulations of general atomic and particulate systems we introduced the Characteristic Crystallographic Element (CCE) norm [29]. The CCE norm is a descriptor able to quantify simultaneously the orientational and radial similarity (or lack of it) of a local configuration with respect to a reference crystal. Furthermore, it can distinguish and differentiate between competing ordered structures. Through its application, we studied crystallization in dense packings of hard sphere (HS) monomers and polymers investigating the effect of volume fraction, chain length and bond gaps [30][31][32][33][34]. Additionally, the structural competition between the formation of close packed crystals composed of face centered cubic (FCC) and hexagonal close packed (HCP) segments and of the fivefold (FIV) local symmetry has been studied for both monomeric and polymeric systems [35][36][37]. The CCE algorithm has also been successfully implemented in the simulation studies of vitrification and crystallization of flexible and semi-flexible polymer systems based on a hybrid bead-spring model [38,39]. A CCE variant, denoted as short-range order symmetry parameter (SSP), has also been adopted to explore short-range order and distortion in bulk metallic glasses [40]. More recently the CCE norm has been applied to study cluster, crystal formation and vitrification, under dilute conditions, of polymers whose sites interact through the square well (SW) potential [41].
The initial implementation of the CCE norm included analysis only with respect to the HCP, FCC and FIV morphologies, a choice justified by the athermal nature and the high packing density of the initial HS systems. Here, we revise and generalize the original CCE norm formalism by expanding the original concept to identify and distinguish among additional reference crystals, including those belonging to the same crystallographic system. In parallel, the new CCE descriptor is readily extended to 2-D systems such as ordinary 2-D assemblies (i.e., of disks) or thin-films of spherical particles under extreme, plate-like confinement where film thickness corresponds to a single layer.
The paper is organized as follows: Section 2 presents the CCE algorithm, its fundamentals and its practical implementation in two and three dimensions, including a short description of the simulation method to create the samples to be used as test cases. Section 3 hosts results from application on reference 2-D and 3-D crystals and on computer-generated, hard sphere and square well configurations. Finally, Section 4 summarizes the main results and lists potential applications of the CCE descriptor.

Characteristic Crystallographic Element Norm
Crystals can be classified and effectively distinguished by means of their point group, loosely speaking, the set of symmetry operations which leave at least one point fixed, under which the crystal is invariant. These symmetry operations or elements of the point group can be intuitively and conveniently organized as being generated by geometric elements of symmetry, such as inversion center, rotation and rotoinversion axes [42][43][44][45].
The CCE algorithm is based on the univocal correspondence between a given crystal and its point group. The identification of the symmetry elements of the arrangement of atoms around a given one permits the assignment of a specific point group to this atom. If the point group, as well as the number of closest neighbors (coordination number), are found to match that of a given crystal, such as FCC, HCP, etc., the atom is declared to have FCC-, HCP-character, etc.
In a strict sense, the search for symmetry elements around a given atom should include all elements of the point group in question; in practice however, it is very often enough to identify a subgroup of the point group in order to discriminate between two competing crystal types.
Previously we have demonstrated [26] that, in practice, identifying a properly selected subset only, does indeed provide the desired accuracy and power of discrimination, while being simpler in algorithmic implementation and faster in computational time.
The CCE norm is always defined with respect to a specific reference crystal X (X-CCE norm). In the present analysis we use as reference templates the following 3-D crystals: hexagonal close packed (HCP), face centered cubic (FCC), hexagonal (HEX) and body centered cubic (BCC), as well as the non-crystallographic fivefold (FIV) local symmetry. In 2-D, the corresponding reference crystals are: triangular (TRI), square (SQU), honeycomb (HON) as well as the pentagonal (PEN) local symmetry. Figures 1 and 2 present fragments of the reference crystals in 3-D and 2-D, respectively, studied here for the demonstration of the CCE descriptor. One aspect that we should consider by comparing the reference 3-D crystals of Figure 1 is the difficulty to identify and distinguish, through visual inspection, each one of them especially in a traditional 2-D image. This further highlights the importance to establish a metric which would allow for a robust characterization of ordered structures. In parallel, visual inspection is enhanced when 3-D images of the computer-generated structure are used. Such images are provided in the Supplementary Material and in the interactive, 3-D version of the present manuscript. red, green, purple and cyan is used throughout the present manuscript. Image created with the Visual Molecular Dynamics (VMD) software [5]. Each figure panel is also available as stand-alone image in 3-D, interactive format in the Supplementary Material. In the continuation, we demonstrate the fundamentals of the CCE norm when applied on a system of N at spherical particles or atoms, where each particle i (i ∈ [1, . . . , N at ]) is defined by its position vector, r i .
First, for a given site i we calculate the number of its nearest neighbors, N(i), and their identities (numerical labels) j (j ∈ [1, N(i)]) through Voronoi tessellation, as implemented in the voro++ software [46] in 3-D and the qhull software [47]    If N(i) < N coord , i.e., if the number of Voronoi neighbors is smaller than the coordination number of the reference crystal X, a penalty function is introduced in the form: where w 0 is a constant whose value is equal to 0.07. This value is empirically determined so that the CCE can differentiate crystals that belong to the same point group but have different number of neighbors and Voronoi polyhedra (for example BCC vs. FCC), as will be explained in the Results section. In computer-generated systems the local environment does not correspond to the Voronoi polyhedron of a perfect crystal due to fluctuations in atom positions. Such typical Voronoi cell has a higher number of vertices and faces as some of the second-nearest neighbors lie closer to the central, reference atom than those of a perfect crystal, i.e., N(i) > N coord (X). In such cases, the penalty function is zero and the CCE proceeds by sorting the neighbors based on their distance from the reference atom and selecting only the N coord (X) closest ones. The following step is also the central one in the CCE execution: given the reference X crystal, the characteristic geometric symmetry elements are identified along with the corresponding symmetry actions for each one of them. As stated earlier, it is not necessary to search for the complete set of geometry symmetry elements for each reference crystal X. It is possible to clearly differentiate between competing and structurally similar crystals by using a subset only. In the current implementation of the CCE norm we consider two geometric symmetry elements: axis (or set of axes) and inversion center (wherever applicable). For example, to distinguish between HCP and HEX-like environments, it is enough to take a single sixfold axis 6 and the five distinct actions 6 6 c = E needs not be taken into account) for the HCP crystal; and a single sixfold rotation axis (6) and its 5 actions: 6 1 c , 6 2 c , 6 3 c , 6 4 c , 6 5 c , plus an inversion center at the reference site i, for the HEX crystal. The nearest neighbors (N coord (X)), the geometric symmetry elements, N el (X), and the corresponding actions for each element, N act (k,X) (k ∈ [1, . . . , N el (X)]) for each X crystal in 3-D and 2-D, are presented in Tables 1 and 2, respectively. The analogous combinations for the fivefold (FIV) local symmetry in 3-D and the pentagonal (PEN) local symmetry in 2-D are also included. Table 1. Number of nearest neighbors, N coord (X), number of distinct geometric symmetry elements N el (X), along with their type, number of actions per element, N act (k,X) and their type for the 3-D crystals (X = HCP, FCC, HEX and BCC) as well as the fivefold (FIV) local symmetry.

Reference
Symmetry Actions of Geometric Element k The geometric symmetry elements and their actions allow us to quantify similarity of the atomic environment with respect to the reference crystal X. This is done by comparing, through the CCE operations, the "real" local environment, formed by the nearest neighbors with coordinates r j , considered relative to the ones of the reference atom, against the "ideal" coordinates R X j ( j ∈ [1, . . . N coord (X)]) of the sites that constitute the coordination polyhedron of the perfectly ordered structure X. In the most general form the CCE norm can be defined as: where ε X,0 i is the penalty function based on the coordination number as defined in Equation (1) and σ is the characteristic length of the system to render the CCE norm dimensionless. In past studies as well as in here, σ is taken to be equal to the diameter of the spherical particles. Summation in Equation (2) runs over all N coord (X) sites, geometric elements N el (X) and N act (k,X) symmetry actions corresponding to the reference crystal X. S X k,m is the orthogonal matrix that performs the m th action of the k th element. For example, the corresponding matrix for the point inversion for the FCC, HEX and BCC crystals is simply: The orthogonal transformation is applied to the position vectors of all neighbors of the reference site and for as many actions and symmetry elements as required by the reference crystal. After its application the new and transformed coordinates of atom j are compared against the positions corresponding to the perfect ordered structure, i.e., the term r j − S X k,m R X j 2 is calculated, selecting the pairs (transformed and ideal) of coordinates that produce the minimum possible discrepancy. Through proper selection of the S X k,m operator the "real" and "ideal" set of coordinates are expressed in the same reference frame: When the geometric symmetry element corresponds to an inversion center, whether it is 2-D or 3-D, its position is fully defined as it coincides with the reference site, i. Similarly, in 2-D the rotation axis is perpendicular to the plane where the atoms or particles lie.
In the general, 3-D case the orientation of the axis(axes) is not known a priori and it has to be determined by a minimization process where the spherical domain is scanned in small steps dϕand dθleading to pairs of polar, θ, (∈ [0, π/2]) and azimuthal angles, ϕ, (∈ [0, 2π]) corresponding to the axis orientation. The minimization scheme is performed for each axis belonging to the characteristic symmetry fingerprint.
For example, for the FCC and BCC crystals the procedure, in practice, is the following: the first symmetry axis S 1 is found by scanning the whole spherical domain (θ 1 , ϕ 1 ); given the S 1 orientation, the second axis is built with an angle of θ 2 = 70.53 • with respect to the first, and a search is now carried out in the azimuthal space: ϕ 2 (∈ [0, 2π/3]). Given S 1 and S 2 the third (S 3 ) and fourth (S 4 ) symmetry axes are fully defined and no further search must be carried out. Each time the geometric element is defined, the analogous symmetry actions (roto-inversions for FCC and BCC) are performed over all coordination neighbors according to Equation (2); as a result, each possible axis orientation, or their combination in case of more than one, provides a different norm value. For reference site, i, and given ordered structure X, the final CCE norm, ε X i , is set to be the minimum of this iteratively refined set of values. By construction, the closer the norm to zero the higher the similarity of the local environment to the ideal X crystal. Moreover, by definition the X-CCE norm of the X reference crystal should be zero, subject to the statistical error imposed by the size of the scanning grid imposed on the spherical domain (see below). Furthermore, by utilizing the characteristic elements and actions as the fingerprint of X-crystal the CCE algorithm is highly discriminating: for two different crystals X and Y, if ε X i → 0 then ε Y i > 0 and vice versa. Once the CCE norm is calculated it is compared against an empirically determined threshold value, ε thres which is the same for all reference crystals. If the CCE norm is lower than this threshold (i.e., ε X i < ε thres ) the reference i site has similarity to that crystal and is thus labelled as X-like. Past studies on bulk, 3-D systems consisting of non-overlapping [29][30][31][34][35][36] or fused [38,39] spheres adopted thresholds of 0.245 and 0.200, respectively. The threshold value is identified through parity (ε X i vs. ε Y i ) plots over all particles in the system as the ones to be presented below. As it will be demonstrated in the results section the same threshold (ε thres = 0.245) can be adopted in the case of hard spheres in 2-D films.
To illustrate specific examples, given the elements and actions, as reported in Table 1, the CCE norm adopts the following forms for HEX (Equation (5)) and FCC (Equation (6)): where S HEX 2,m and S FCC 5,m correspond to the inversion orthogonal matrix (Equation (3)). S HEX 1,m and S FCC k,m transform coordinates according to six-fold rotations and three-fold roto-inversions, respectively, and the minimizations are carried out over axis/axes orientations. The CCE analysis, as described above, is applied to all atoms/particles of the system and for all required reference crystals. For the ordered structures demonstrated here (HCP, FCC, HEX and BCC) due to the increased number of geometry elements (four roto-inversion axes) the CCE descriptor is significantly more expensive computationally for the crystals of the cubic class (FCC and BCC). Among them, the FCC analysis is slower than the one for BCC because of the higher coordination number (12 vs. 8). The required computational time as well as the precision (numerical error) in the calculation of the norm are both significantly affected by how fine the scanning grid is in the minimization in the multi-parameter space of all possible polar and azimuthal angles. Obviously, finer scanning grids result in higher accuracy but also increased processing time.
The CCE analysis described above is repeated over each atom/particle of the system, i (i ∈ [1, . . . , N at ]) and with respect to any desired reference crystal, X. After this information is accumulated over all atoms, the probability distribution function of norms, P(ε X ), is computed for each recorded system configuration (frame). Given that once the CCE norm value is below a certain threshold (ε X i < ε thres ) the site i is labeled as having X similarity, an order parameter with respect to the reference crystal X, S X , can be calculated: with S X ∈ [0, 1]. The value of the order parameter (Equation (7)) is specific to the given crystal or local structure X. If we focus on the reference crystals, excluding pentagonal symmetry, in a system of dimensionality d (d = 2 or 3) N s,d corresponds to the number of such crystal templates. In the present demonstration of the CCE norm N s,3 = 4 (HCP, FCC, HEX and BCC) and N s,2 = 3 (TRI, SQU, HON). Additionally, a degree of ordering or crystallinity, τ c , can be calculated as the sum of all CCE order parameters: where index k runs over all different reference crystal structures. For the atomic systems to be used for the demonstration of the CCE descriptor, Equation (8) becomes: where crystallinity in Equations (9) and (10) corresponds to 3-D and 2-D systems, respectively.

Molecular Simulations
All systems have been generated and equilibrated through extensive Monte Carlo (MC) simulations. The protocol behind these simulations is a generalization of the one presented in Refs. [41,48,49] containing chain-connectivity-altering moves, cluster identification/displacement algorithms, and compression/volume fluctuation processes in the bulk or under various conditions of confinement. In the present work, two different models have been used to describe interactions between atoms. For a pair of atoms i and j, whose centers lie in a distance r ij , according to the hard sphere (HS) model the energy is either zero or infinite: where σ 1 is the collision diameter, also taken as the characteristic length of the system. The second model is the square well (SW), which further includes an attractive range (if r ij ∈ [σ 1 , σ 2 ]) with a coresponding intensity, ε SW : Initial simulations are based on the HS potential on linear chains of tangent hard spheres of uniform size under constant volume in the bulk [48] or under confinement. The latter is realized through the presence of parallel, flat and impenetrable walls in one dimension [49]. Extreme conditions correspond to the case where the distance between the walls, d wall , is equal (within a tolerance of 10 −4 ) to the sphere diameter, σ 1 . Once this situation is met the system is practically converted into a 2-D, thin film whose thickness corresponds to a single layer (d wall → σ 1 ). Two different sets of simulations were carried out containing 100 chains of average length N av = 12 and 50 chains of N av = 24. Due to the application of chain-connectivity-altering moves dispersity is introduced in chain lengths, which fluctuate uniformly in the intervals N ich ∈ [6,18] and [12,36] for N av = 12 and 24, respectively.
In all 3-D cases, initial system configurations are borrowed from fully equilibrated athermal polymer samples at very low packing density (φ = 0.05). Then, two patterns are followed: (i) a SW potential is activated on the polymer configurations at dilute (3-D) or dense (2-D) conditions under constant volume and temperature (T = 1/k B , k B being the Boltzmann constant); (ii) all bonds are eliminated so that the system consists now of monomers and the SW potential is further activated. For the polymer-based systems, the equilibration algorithm and the corresponding MC moves are similar to the ones used for HS [30,33,34,[48][49][50][51][52][53] and SW [41] chains. For the monomer-based systems two moves are employed: simple sphere displacements whose amplitude is auto-configured based on the acceptance rate and cluster-based displacements following the cluster detection/translation algorithm presented in Ref. [41]. Simulations on 3-D SW polymers and monomers are conducted under conditions of constant volume (NVT ensemble) and pressure (NPT ensemble), respectively.
2-D, thin-film configurations are generated either from fully confined 3-D analogs according to the wall-wrapping algorithm described in [49] or through adsorption onto the surface with simultaneous repulsion from the opposite side until full surface coverage is achieved. Two different realizations have been used: 48-chain N av = 100 system at φ = 0.40 under full confinement and 100-chain N av = 12 system at φ = 0.48 with confinement imposed on the short dimension and periodic boundary conditions enforced on the long ones. In both cases, the wall thickness is approximately equal to the monomer diameter, effectively corresponding to 2-D thin films. Successive steps include compaction of the athermal HS chain systems to even higher volume fractions or activation of SW attractive interactions under constant volume (chains) or constant pressure (monomers).
Once simulation trajectories are collected, the CCE analysis is performed on every frame or using regular intervals (i.e., every 100 frames) depending on the system under study. First, based on the atomic coordinates and the presence of periodic boundary conditions or impenetrable walls, a Voronoi tessellation is performed through the voro++ [46] and qhull [47] programs for 3-D and 2-D systems, respectively. Then, the CCE descriptor is employed, in a trivially and massively parallel manner, with all trajectory frames being analyzed simultaneously. Technically, the present implementation of the CCE algorithm allows the inclusion/exclusion of specific geometric symmetry elements, for example the points of inversion as presented in Tables 1 and 2.
The output of the CCE analysis consists of a file with a list of the atom index, i, coordinates, r i , the (minimum) CCE norms with respect to all tested reference crystals, ε X i , and optionally the orientation(s) of the corresponding symmetry axi(e)s. Additionally, a pair of pdb/psf files is provided to be used for successive visualization through appropriate software (for example VMD [5], as used in the present work).

CCE Norm Application to Perfect Crystals
According to the definition of the CCE descriptor, when applied to local structures corresponding to a perfect crystal X, the corresponding norm should be equal to zero (ε X i = 0). However, in practice a discretization error is associated with the CCE application in the case of 3-D systems where the orientation of the axis (or axes) is not known a priori. As explained in the methods section, the spherical coordinate system has to be divided into a discrete mesh. The finer the mesh, the lower the associated error but the higher the required computational time. Except otherwise stated, all results to be presented in the continuation have been obtained with a discretization step of 0.1 rad for all polar and azimuthal angles needed for the identification of the symmetry axi(e)s that minimize the CCE norm.
As a first test of validity we apply the revised and extended CCE descriptor on selected reference crystals: HCP, FCC, HEX and FIV in 3-D, as seen in Figure 1 (segments of the infinite crystal) and Figure 3 (local environment), and TRI, SQU and HON in 2-D, as seen in Figure 2 (infinite crystal) and Figure 4 (local environment). Tables 3 and 4 host the results from the CCE application to perfect reference crystals in 3-D and 2-D, respectively. In the CCE analysis we have further included the non-crystallographic fivefold (X = FIV) and pentagonal (X = PEN) local symmetry in 3-D and 2-D, correspondingly. Inspection of the data in both tables reveals that the X-CCE produces, as expected, a norm equal to zero when applied on the X reference crystal. The only two exceptions to this rule are the FCC and BCC crystals where the corresponding norms are very low (ε FCC = 2.44 × 10 −4 and ε BCC = 1.90 × 10 −4 ) but non-zero. This is because both crystals belong to the cubic system and the fingerprint corresponds to a set of four threefold roto-inversion axes. The determination of the orientation of the first axis follows the same procedure as for the other crystals but another search (plus discretization) in space is required to establish the orientation of the second axis, as explained in the methods section. Accordingly, for a given mesh size (0.1 rad) it is more difficult to find the optimal set of axes for the BCC and FCC cases than the other crystals. Unavoidably, these crystals will have higher uncertainty in their detection compared to the other ordered motifs. However, we should note that the established values are already 10 3 times lower than the identification threshold, ε FCC = 0.245.
The dependence of the value of the CCE norm on mesh discretization, which is the same for any polar or azimuthal angle required to identify the orientation of the first and second axes for the FCC crystal is presented in Figure 5. The increase in accuracy in the detection, as quantified by the reduction in the corresponding CCE norm is not monotonic for specific mesh intervals. This is because the axis orientations may be closer to specific values of the polar and azimuthal angles even if these are produced by a coarser grid of the 3-D space. For a discretization step of 0.1 rad, the required time for the CCE analysis per atom is 0.5 s. The current algorithmic implementation of the CCE norm allows for a trivially parallel execution so that the computational time of the analysis of the simulation trajectory is only marginally longer than the time required per frame.  Another point to notice is the ability of the CCE descriptor to distinguish between competing crystals. The X-CCE norm when applied on the Y reference crystal (Y X) in 2-D provides a norm which, in the vast majority of the cases, is significantly higher than the threshold value. For 3-D there are X-Y combinations where the value is lower than the threshold. In this case, which is only very seldom encountered in computer-generated configurations, the site is labelled according to the lowest norm. Furthermore, we should note that when the HCP-CCE descriptor is applied on the perfect HEX structure the corresponding value (0.28) is a result of the ε HCP,0 i contribution in Equation (2) as calculated by Equation (1). This is because the Voronoi polyhedron of the local environment that corresponds to the perfect HEX structure, as seen in the third panel of Figure 3, is a hexagonal prism with eight faces and thus eight closest neighbors. When the actions of the 6-fold roto-inversion axis of the HCP are applied on the HEX Voronoi cell they produce zero HCP-CCE norm according to Equation (5). However, as the HEX Voronoi polyhedron has 8 neighbors, which are fewer than the 12 ones used for the HCP-based operations, according to Equation (1) we have ε HCP,0 i = 0.07 (12 − 8) =0.28. So that the HCP-CCE norm when applied on the reference HEX crystal has this high value (0.28) as seen in Table 3. Again, this is a particular case of an "ideal" (reference) structure that is never encountered in practice in computer-generated configurations. In the latter the number of closest neighbors, as calculated by the Voronoi tessellation, is either equal or higher than the one that corresponds to the perfect HEX (as will be demonstrated in the continuation).

CCE Norm Application to Computer-Generated, 3-D Bulk Systems
Our test cases consist of computer-generated packings of monomers or polymers whose sites interact through the hard-sphere (HS) or the square-well (SW) potential. In all systems studied here the collision diameter is taken as the characteristic length of the system (σ = σ 1 in Equation (2)). The aim of the present work is not to present new physics on 2-D (thin film) or 3-D (bulk) crystallization of athermal or attractive chains but rather to demonstrate the efficiency of the CCE algorithm to accurately identify local structure and to detect and quantify crystal nucleation and growth in general particulate and atomic systems. Initial configurations correspond to bulk, 3-D HS chains at low packing density (φ = 0.05), as can be seen in Figure A1 of the Appendix A.
Once the SW potential is activated, chains or monomers form clusters, whose stability depends strongly on the values of the intensity and range of attraction [41]. Under specific conditions the clusters grow in size and aggregate, eventually leading to a single cluster which contains all sites. Crystal nucleation and growth inside such clusters, as well as the type and morphology of the ordered structure, also depend strongly on attraction intensity and range. Typical configurations for SW polymers and monomers at the end of the NVT (chains) and NPT (monomers) simulations, after application of the CCE algorithm, are presented in Figures 6 and 7.
Sites are color-coded according to their similarity to an X reference crystal, as identified by the value of the CCE norm. As mentioned in the method description the condition for a local environment to be labeled as X-like is that the value of the corresponding CCE norm should be lower than the threshold: ε X < ε thres = 0.245. Throughout all images and snapshots to be presented in the continuation the coloring scheme is the same: red, blue, green, purple and cyan colors correspond to HCP, FCC, FIV, HEX and BCC crystals, respectively. All sites that show no similarity for none of the reference crystals are labeled as "amorphous", or more precisely as "none of the above", and are represented in yellow. In parallel, to render the depth of the snapshot visually accessible, crystal-like and amorphous sites are shown with reduced dimensions in a 3:5 and 1:5 scale, respectively. In the specific examples, the established ordered morphologies, as identified by the CCE algorithm, are almost perfect HCP and FCC crystals for the SW chain systems at dilute conditions ( Figure 6) while HEX and BCC are the dominant ones for the SW monomeric clusters in Figure 7. A detailed analysis on the phase behavior of chains and monomers in the whole spectrum of the intensity and range parameters as well as an explanation on why certain ordered morphologies are favored for specific simulation conditions will be presented in a future work. A significant difference between the low-density chain crystals, as obtained from the NVT simulations, and the high-density ones for monomers, as obtained from the NPT simulations, is the presence in the former of an external surface of unidentified character that surrounds the core of the cluster as seen in the snapshots of Figure 6. Atoms on the surface of the cluster cannot be assigned to any reference crystal because they lack a complete set of first neighbors, and hence have incomplete Voronoi polyhedra [41].  Focusing on the level of atoms or particles dozens of final configurations of MC trajectories were scanned in an effort to discover local structures with an environment as close as possible to the reference tested crystals. Such X-like structures should be characterized by very low values of the corresponding X-CCE norm. Figure 8 hosts some of the "best" computer-generated local structures with respect to the four reference crystals tested. The corresponding CCE values are all very low with the minimum being registered for an FCC-like environment (ε FCC = 0.0223), followed by HCP-like (ε HCP = 0.0340), HEX-like (ε HEX = 0.0466) and BCC-like (ε BCC = 0.0536), the latter being approximately 2.4 times higher than the minimum recorded for FCC. All these CCE norm values are very small, close to zero, and significantly lower than the threshold value of 0.245. Visual inspection of the computer-generated local environment as seen in Figure 8, especially when compared against the reference perfect structures of Figure 3, verifies the CCE results: the similarity is striking not only at the level of nearest neighbors but also with respect to the shape of the Voronoi cell. We should remind here that out of all neighbors that form the Voronoi polyhedron only the nearest N coord (X) ones are taken into account for the calculation of the X-CCE norm.
How the CCE descriptor quantifies increasingly poor similarity to a specific reference structure can be seen in Figure 9. The four snapshots correspond to different local structures whose HCP-CCE norms are approximately ε HCP = 0.050 (left-most panel), 0.10, 0.20 and 0.30 (right-most panel). The lowest HCP-CCE norm is linked to an almost perfect HCP ordered structure. As the norm increases small structural defects appear, including minor radial and/or orientational deviations. This is most evident in the rightmost snapshot with ε HCP = 0.30. This value is actually higher than the threshold (0.245) used to assign HCP character. Accordingly, based on the HCP-CCE norms in Figure 9, all structures are identified as HCP-like except the one on the right-most panel. Increasing (weakening) the criterion for the detection to a value higher than 0.30 would also include this structure in the HCP population. However, caution should be exercised as too large an increase in ε thres may result in sites having dual or multiple similarity eliminating the ability of the CCE method to discriminate between competing, structurally similar crystals. It is tempting to attempt a comparison between the characteristics of the environment of the perfect reference crystals against the one of some of the "best" computer-generated structures, as for example the ones depicted in Figure 8. Parameters related to the topology of the enclosing Voronoi polyhedron for each reference crystal, X, include: coordination number, N coord (X), volume, V VP (X), surface area, A VP (X), number of faces, F(X), number of edges, E(X) and number of vertices, V(X). Furthermore, local packing density can be calculated as: which for spherical objects, of radius r = σ 1 /2, is connected to packing density through: ϕ(X) = π 6 ρ n (X) = π 6V VP (X) (14) In parallel, for 3-D objects, the standard isoperimetric quotient, Q(X), is defined as [54]: For convex polyhedra faces, edges and vertices are connected through the Euler-Poincaré equation: F(X) + V(X) = E(X) + 2. The space filling tessellations of the reference crystals studied here correspond to trapezo-rhombic dodecahedron (HCP), rhombic dodecahedron (FCC), hexagonal prism (HEX) and truncated octahedron (BCC) with the corresponding parameters presented in Table 5. It is interesting to check whether the absolute values but also the well-established trends for the reference crystals-φ(HCP) = φ(FCC) > φ(BCC) > φ(HEX) and Q(BCC) > Q(HCP) = Q(FCC) > Q(HEX)-will be also valid for the computer-generated, defect-ridden counterparts. As discussed earlier, the enclosing Voronoi polyhedron of a computer-generated structure, even a highly ordered one, is more complex than the one of the reference structures because the second-nearest neighbors are brought to closer distances than the ones who adopt regular positions in crystal lattices. This can be seen for example in Table 6 where the number of faces, edges and vertices is systematically higher in the computer-generated structures with high crystal similarity compared to the reference templates (Table 5). If we focus our attention on the volume of the Voronoi cell and accordingly the packing density, all computer-generated structures are less dense than the perfect crystals. In the best comparison which corresponds to the FCC-like crystal, the difference in volume is approximately 3.4% with a similar but higher value (6.6%) being registered for the HCP-like; in the worst cases for HEX-and BCC-like it increases to 19.9 and 19.7%, respectively. The relative percentage differences for V VP (X), A VP (X) and Q(X) between the computer-generated structures of high crystal similarity and the reference ones, are summarized in Table 7. A remarkable trend, which must be studied in more detail through an extended sample of configurations, is that the generated HEX and BCC crystals are significantly more spacious than the perfect analogs. In parallel, the standard isoperimetric quotient seems to be in very good (3% for HEX) to an almost perfect (0.01% for FCC) agreement between the generated structures and the analogous perfect templates. The computer-generated BCC structure is still, with a slightly reduced margin compared to the perfect templates, the local structure with the highest isoperimetric quotient value, even larger than the ones of the close packed HCP and FCC crystals (Q(BCC) > Q(HCP) ≈ Q(FCC) > Q(HEX)). Table 5. Statistics of the Voronoi polyhedron of perfect, reference 3-D crystal X: coordination number N coord (X), number of faces F(X), vertices V(X) and edges E(X), volume V VP (X), surface area A VP (X), local number density ρ n (X), packing density φ(X) and standard isoperimetric quotient Q(X).

X N coord (X) F(X) V(X) E(X) V VP (X)
A VP (X) ρ n (X) φ(X) Q(X)  Table 6. Same as in Table 5 but for computer-generated 3-D structures with high crystal similarity (or equivalently low CCE norm), as visualized in Figure 8. Coordination number is deduced from the number of closest neighbors so as to match the one of the reference crystal.

X-like N coord (X) F(X) V(X) E(X) V VP (X)
A VP (X) ρ n (X) φ(X) Q(X)  Table 7. Relative percentage difference between the computer-generated local structures with very high crystal similarity (as visualized in Figure 8) and the reference, perfect 3-D crystals with respect to volume, surface area and standard isoperimetric quotient of the enclosing Voronoi polyhedron. Also reported in the last column is the corresponding X-CCE norm for the computer-generated structure. The ability of the CCE descriptor to discriminate between similar and competing crystal structures can be demonstrated by the parity plots where the Y-CCE norms are contrasted against the X-CCE norms (X Y), in all possible combinations of crystal pairs for all atoms present in the system. Such parity plots can be seen in Figure 10 as obtained from the last configuration from NPT MC simulations on 1200 monomers interacting through the SW potential (ε SW = 2.1 and σ 2 = 1.2). Dashed horizontal and vertical lines denote the threshold (ε thres = 0.245) below which a site is identified as Y-like and X-like, respectively. The discrimination ability of the CCE algorithm is demonstrated by the empty square whose borders are defined by the points (0, 0), (0, ε thres ), (ε thres , 0) and (ε thres , ε thres ). In practice, absence of data points in this area implies that no site possesses dual ordered character by having simultaneously low CCE norm with respect to the two different X and Y reference crystals. Vertical and horizontal dashed lines denote the CCE norm threshold below which a site is labeled as Xand Y-like, respectively. Blue, red, green, purple, cyan and yellow colors correspond to HCP-, FCC-, FIV-, HEX-, BCC-like and amorphous (or unidentified) sites, respectively.

CCE Application to Computer-Generated, 2-D Thin-Film Systems
In the present section, we will demonstrate the application of the CCE descriptor on 2-D systems. As explained earlier these systems could be inherent 2-D systems consisting for example of disks, or extremely confined thin-films of 3-D objects (like spheres), with a film thickness of a single layer (d wall → σ 1 ). Local packing density, ρ * n (X), is defined here as the inverse of the surface area, A VP (X), of the Voronoi polygon of crystal X: Thus, surface coverage, ϕ * (X), is the ratio of the occupied area divided by the total area of the Voronoi polygon of crystal X: A thin film of one-layer thickness, and thus any sub-volume of it, is characterized by a packing density as defined by Equation (14). Combining Equations (14) and Equation (17) for the extremely confined thin film (d wall → σ 1 ) we get: If A VP (X) and P VP (X) are the surface area and the perimeter of the Voronoi polygon for the 2-D crystal X the standard isoperimetric ratio, q(X), is defined as [54]: Table 8 summarizes the properties of the Voronoi polygons for the reference 2-D crystals studies here: TRI (regular hexagon), SQU (square) and HON (equilateral triangle). Table 8. Statistics of the Voronoi polygon of perfect, reference 2-D crystal X: coordination number N coord (X), number of vertices V(X) and edges E(X), surface area A VP (X), perimeter P VP (X), local packing density ρ * n (X), surface coverage ϕ * (X) and standard isoperimetric ratio q(X).

X N coord (X) V(X) E(X) A VP (X) P VP (X)
ρ * n (X) ϕ * (X) q(X)  Figures 11 and 12 present snapshots of MC simulations on thin-film systems of linear chains whose sites interact with the HS and SW potential, respectively. In all cases, atoms are colored according to the X-crystal similarity as quantified by the 2-D norm: blue, red, green and cyan colors correspond to TRI, SQU, PEN and HON symmetries, respectively, while amorphous (or unidentified atoms) are shown in yellow. Given that in all cases film thickness is approximately equal to the diameter of monomers packing density can be related to surface coverage according to Equation (18). The 48-chain N av = 100 system under full confinement (left panel in Figure 11) is characterized by low density (φ = 0.40) and equivalently low surface coverage (ϕ * = 0.60). As a consequence, only a very small fraction of sites shows ordered structure. As density increases (right panel in Figure 11 corresponding to 100-chains of N av = 12 at φ = 0.48 or equivalently ϕ * = 0.472) the population of sites with crystal local structure increases appreciably, especially that of the TRI character.
The situation is quite different as the SW potential is turned on and simulations are conducted under constant volume; this is as demonstrated by the snapshots of Figure 12. The configuration of the left panel (ε = 0.5 and σ 2 = 1.3) shows a clear tendency for close packing, as most sites adopt a triangular structure, which possesses the highest density among all 2-D crystals. The snapshot of the right panel (ε = 0.5 and σ 2 = 1.6) also presents higher degree of ordering than the original athermal packing (right panel of Figure 11) and different mixture of crystal sites. Here, the population of square-like sites is comparable to that of triangular ones. In both cases, no honeycomb structures are detected, an expected trend given that the honeycomb crystal shows significantly lower surface coverage than the square and triangular ones. Additionally, in the configuration of mixed character (TRI and SQU) there is a small, but non-zero population of sites with non-crystal, pentagonal local symmetry. In contrast, in the TRI-dominated system, pentagonal sites are completely absent. A systematic analysis of the phase behavior as a function of attraction intensity and range in thin films of SW chains as well as an interpretation of the observed tendencies will be presented in a future work.   Figure 13 hosts local structures belonging to the athermal N av = 12 system (which can be seen in Figure 11, right panel) with progressively reduced similarity to triangular crystal as quantified by the TRI-CCE norm. From left to right the values are approximately 0.05, 0.10, 0.20 and 0.30. The last structure with ε TRI = 0.30 is not recognized as TRI-like, as the norm is higher than the threshold value of ε thres = 0.245. Visual inspection qualitatively confirms the diminishing similarity as established quantitatively by the CCE descriptor. This is further demonstrated by the data on the statistics of the Voronoi polygon, as reported in Table 9, where they are further compared against the corresponding one of the perfect, reference triangular crystal. It can be deduced that the computer-generated local structure with the closest TRI similarity differs by approximately 5.8, 2.9 and 0.2% in surface area, perimeter and standard isoperimetric ratio with respect to the reference crystal. For the least similar structure the corresponding percentages (of the ones reported in Table 9) increase to 15.5, 9.0 and 2.9%, for A VP (TRI), P VP (TRI) and q(TRI), respectively. In both cases and as expected, the computer-generated polygons are larger than the reference ones and thus less dense.  Table 9. Statistics of the Voronoi polygon of computer-generated local structures, visualized in Figure 13, with progressively poor triangular singularity as quantified by the corresponding CCE norm, ε TRI : coordination number N coord (X), number of vertices V(X) and edges E(X), surface area A VP (X), perimeter P VP (X), local packing density ρ * n (X), surface coverage ϕ * (X) and standard isoperimetric ratio q(X). Also shown for comparison in the first row are the results that correspond to the perfect, reference 2-D triangular crystal. The discriminating ability of the CCE algorithm is demonstrated by the parity plots for all possible reference pairs as shown in Figure 14. As in the case of 3-D analogs, for every site (atom or particle) in the system we plot the Y-CCE norm versus the X-CCE one (Y X) in all possible combinations. An empty box near the origin, whose borders are marked by the threshold value (ε thres = 0.245), means that no site exists with dual crystal character. Out of all possible comparisons of crystal pairs (= N at × N s,2 ! = 1200 × 6 = 7200) there is only one site which is detected with dual character (ε TRI , ε HON < ε thres ) as seen in the left panel of Figure 14, with the HON similarity (ε HON ≈ 0.24) being the weakest one and very close to the detection threshold. In the very rare case of a site having two CCE norms lower than the threshold, it adopts the one that has the lower value. Vertical and horizontal dashed lines denote the CCE norm threshold below which a site is labeled as Xand Y-like, respectively. Blue, red, green, cyan and yellow colors correspond to TRI-, SQU-, PEN-, HON-, and amorphous (or unidentified) sites, respectively.

Conclusions
We have presented a revised an extended version of the characteristic crystallographic element (CCE) descriptor, used to characterize local structure in general atomic and particulate systems in 2 or 3 dimensions, in the bulk or under confinement. The CCE descriptor is able to quantify the structural similarity of any given environment with respect to reference crystals in both radial and orientational terms. Its definition is based on the characteristic set of geometric symmetry elements and corresponding actions which are uniquely related to a reference crystal. The CCE norm is effectively able to quantify phase transition and to distinguish between different competing crystals even if they are structurally similar and belong to the same point group.
We have demonstrated the efficiency of the CCE norm by applying it on atomic systems of polymers or monomers interacting with the hard sphere and square well potentials. We have tested structural similarity with respect to hexagonal close packed, face centered cubic, fivefold, hexagonal and body centered cubic in 3-D and triangular, square, pentagonal and honeycomb in 2-D.
The CCE descriptor can be easily extended to more crystal templates. Moreover, in 3-D, its precision and speed, or better their balance, can be tuned by the fineness of the grid used to scan the spherical domain for the detection of the optimal symmetry elements (axes).
Finally, the CCE descriptor can be also applied to 3-D systems showing heterogeneous crystallization for example due to the presence of walls. This can be achieved by treating differently the atom layers touching the walls (using 2-D crystal templates) and the bulk system (using 3-D reference crystals). It can be also used to study the connection between local structure and relaxation in out-of-equilibrium systems as the ones reported in Refs. [55,56].
The analogous initial configurations for the simulations on 2-D, thin films are shown in Figure A2. On the left panel a 48-chain N av = 100 system at φ = 0.05 is shown under confinement in all dimensions. The MC protocol to produce such polymer systems under extreme confinement is described in [49,57,58]. The right panel hosts a snapshot of a 100-chain N av = 12 system at φ = 0.48. Here, periodic boundary conditions are applied on the long dimensions of the cell while flat parallel, impenetrable walls exist in the short dimension. In both systems of Figure A2 the thickness in the short (x) dimension corresponds to one layer, i.e., (d wall → σ 1 ).