The determination of point groups from imprecise molecular geometries

We present a new approach for the assignment of a point group to a molecule when the structure conforms only approximately to the symmetry. It proceeds by choosing a coordinate frame that minimises a measure of symmetry breaking that is computed efficiently as a simple function of the molecular coordinates and point group specification.


Introduction
In principle, the discovery of the point-group symmetry of an isolated molecule with known atomic coordinates is a straightforward task that can be achieved with simple algorithms that search for particular symmetry elements and eliminate candidate groups based on the outcome. However, except in very simple cases, there are usually complicating factors associated with the provenance of the molecular structure. If it has been obtained from a crystal structure measurement, or from the numerical minimisation of an energy function, there will inevitably be noise in the structure that results in some or all of the likely actual symmetry elements being formally absent. Further complication can arise from the coordinate frame in which the structure is expressed; the task of finding symmetry elements involves varying the position of the coordinate origin and the orientation of the coordinate axes until a match is found. Again, in simple cases, especially if the underlying symmetry is in fact exact, orienting the molecule so that its axes and planes are in the expected places may be straightforward, but in unfavourable cases care is needed.
The imperfect satisfaction of symmetry relations arising from numerical noise is typically desirable to condone, because the structure is taken to stand for an ideal exact-symmetry structure. However it may also arise that the molecule has a structure that really is symmetry broken, in the sense that it is close to satisfying all the elements of some particular point group, but would remain distorted even if the precision of data were improved. It then becomes important to have clear criteria that allow one to distinguish unambiguously, against some threshold-based measures, these two possibilities.
The work of Avnir and coworkers [1][2][3][4][5][6][7][8][9][10] is largely based on discovering approximate symmetry elements, such as rotational axes, by numerical inspection of approximate motifs, such as regular polygons, in subsets of the atoms. Once an approximate motif has been found, a corresponding exact form is generated by averaging of coordinates; a symmetry measure is then constructed as mean square displacement of the actual structure from the idealised.
Largent, Polik and Schmidt adopt a somewhat simpler approach, again based on the atomic coordinates, but using the actual desired point group operations mapped onto the molecular frame in order to calculate the aggregate deviation from exact symmetry [11]. The coordinate transformation is defined by the principal inertial axes for asymmetric tops, with additional partially-specified criteria for symmetric and spherical tops. An efficient computational implementation is reported. Similar schemes are implemented in other software packages [12][13][14].
Closely related is the analysis of approximate point group symmetry from the viewpoint of fuzzy sets [15,16], and using similarity measures based on electron density [17][18][19][20], and introducing the syntopy concept [21]. These methods have in common the need to already know the electron density, or even the potential energy surface, for a particular electronic state. The latter [15], in particular, is appealing, since it can give information on whether or not significant symmetry breaking is in play against an energy scale defined by temperature. These methods offer important insights and quantification that can be applied to any molecule. However, in the context of the present work, these approaches cannot be used directly, since we seek an analysis of near symmetry based only on the coordinates of the nuclei, typically to be used before an electronic structure computation is carried out.
Oakley et al. [22] describe use of the symmetry measure of Ref. [1] as a continuous function of the relative orientation of the molecule and point group coordinate frames that can be minimised by reorientation. This approach, which could be superior to simple use of inertial axes when symmetry is not exact, is developed further in the present work.
Grimme [23], and later Casanova and Alemany [24], adopt an approach where the focus is on quantum-mechanical wavefunctions. A symmetry measure can be defined based on the overlap of the wavefunction with its image induced by the action of the operator on the electronic coordinates, which deviates from unity as symmetry is broken.
In this paper, we introduce an alternative scheme that extends and combines these previous approaches, and defines measures of broken symmetry by providing clearly-defined algorithms that are invariant to coordinate system, and which support the automatic evaluation and discovery of point-group symmetry. The approximate symmetry is searched for by optimisation of a symmetry measure with respect to the position and orientation of the molecule, with the target symmetry specified entirely by the definition of the particular point group being tested, rather than looking for specific features in the atomic coordinates, leading to a modular and robust algorithm. We also discuss the proper refinement of atomic positions in order to satisfy the elements of a chose point group more precisely. The work is accompanied by a freely-available software implementation.

Symmetry measure
The molecule is defined by the coordinates and charges of a set of atoms, { ⃗ R A , Z A } for A = 1, … , N . We then consider a set of symmetry operators G would normally consist of all of the elements of a particular point group, but need not do so, as none of the closure or other group properties will be used; it could as a special case consist of just a single operator. Each operator acting on each atom gives rise to an image and if the molecule conforms exactly to G, ⃗ R At ∈ G for all A, t. Since symmetry operators effect specific geometric transformations, we must be able to specify the location and orientation of any axes, planes, etc., and we will assume initially that the operators are defined in a fixed global coordinate system with cartesian axes {⃗ e x , ⃗ e y , ⃗ e z } via real unitary representation matrices that act on the components in this coordinate system of any given position vectors: In order to describe inexact conformity with the group, we define for each image its matcher, the closest atom to the image, in which deviations of d At from zero indicate symmetry breaking. The deviations can be used to define an overall continuous symmetry measure, for example [1] (

At
Following instead the signature of symmetry in the electronic wavefunction [23,24] where the symmetry operator in wavefunction space is defined as If T t is an exact symmetry operation, the n-electron wavefunction is an eigenfunction with unit eigenvalue; otherwise, its expectation value is always less than one.
The coordinate-and wavefunction-based approaches can be unified by adopting a simple model for the wavefunction that reflects most of the features of the symmetry deviation of the true wavefunction induced by displacement of the nuclei from exact symmetry positions. We start with the case of a collection of hydrogen atoms-or more generally, a system where in the vicinity of the nuclei, the wavefunction is dominated by a single orbital, with the overall wavefunction an orbital product. Each orbital has the form multiplied by an appropriate damping factor giving zero value near any other nucleus, and near the nucleus, the overall wavefunction is the product of this orbital and n − 1 slowly-varying normalised orbitals. We adopt (7), but drop the n −1 factor; this gives a formulation that is formally extensive (if all atoms show the same symmetry deviation), and locally intensive (if additional atoms that obey the symmetry exactly are added, the measure does not change). We then obtain This agrees with the simple point-charge formula except that near each nucleus, the distance unit is the inverse of the nuclear charge. Additionally, with the closed formula (11), the measure is bounded, 0 ≤ F 0 ≤ N g , supporting interpretation of the value of F 0 . Note that this formula is additive in both atoms (extensivity) and symmetry operations (descent to subgroups), allowing one to look at a chain of symmetry groups and evaluate the effect of lowering symmetry. As an example, water with one bond infinitely stretched gives F 0 = 0 in C s , F 0 = 2.56 in C 2v and F 0 = 8.26 in D 2h because of additional contributions from new operators.
For non-hydrogenic atoms, the model can in principle be extended. A simple way to do this is to recognise additional non-zero occupied orbitals near each nucleus that, because of the Coulomb singularity, will have the local shape of the 2s, 3s, 2p, … orbitals of the 1-electron ion. The overlap integrals are similar in form to that in (10), varying as d 2 At but with a different coefficient [25]. Since the symmetry measure is meant to be only an arbitrary, but well-defined, quantity, we have chosen to ignore the consequences of many-electron wavefunction structure, and use (27) for all atoms without modification. This will tend to underemphasise the importance of positional deviations of heavier elements, but compared to the simple distance-squared criterion, they already have an additional Z 2 A weight.

Optimisation of symmetry operators
In principle, the measure defined above can be calculated for any suspected point group, and if the value is less than some chosen threshold, the molecule can be taken to be belonging to that group. This opens the way to using the standard point-group decision tree found in many textbooks [26] to systematically discover the highest-order compliant group. However, the molecular coordinates may happen to be expressed in any reference frame that is completely unrelated to the coordinate system used to define G, and it is then quite unlikely that a small measure will result. One approach to this challenge is to first ensure that the definition of the group is through operators that are in some standard orientation. For example, for an axial group, an obvious choice is for the unique high-order axis to be aligned along ⃗ e z , and for some of the other axes or reflection plane normals to be chosen to coincide with ⃗ e x or ⃗ e y . Then the molecule can be translated so that its centre of charge is at the origin, and rotated so that the axes coincide with appropriate principal axes-for example, the eigenvectors of the second moment tensor. This works well for molecules with exact symmetries that are asymmetric tops, i.e., the three eigenvalues of the second moment tensor are distinct, but is incomplete when there are degeneracies (symmetric tops and spherical tops). If the molecule has only approximate symmetry, it is not obvious that this origin and axis choice is the best one.
Instead, we adopt an approach that is agnostic with respect to the frames in which the point group and molecule are specified, and effect a relative realignment that is varied until the symmetry measure F is minimum.
We choose to do this by leaving the molecular coordinates untouched, and specifying a rotation and translation of frame holding the symmetry operations.
We need to differentiate F with respect to the parameters that define the coordinate system (origin ⃗ o , axes ⃗ a , = 1, 2, 3 ) in which the operators are defined. They consist of = {o 1 , o 2 , o 3 , q 1 , q 2 , q 3 } defined by where ⃗ are the vectors defining the global coordinate system, and is a function that produces a unitary matrix from three unconstrained parameters. Consideration where we make use of the matrix representation of the symmetry operator in local coordinates, t is fixed by the nature of the symmetry, but T t ⃗ b depends on it through the parameters .
We can now express the image errors in terms of the local symmetry matrices, which differentiate as We can then proceed with and then varying to minimize F using the Broyden-Fletcher-Goldfarb-Shannon (BFGS) algorithm [27]. In practice, in order to increase the convexity of the objective function, we minimize F 1 which uses f 1 (x) = x 2 instead of f 0 (x) specified in (11); if the optimum F is small, the minima of F 0 and F 1 will be close. Computational efficiency can be improved without loss of generality by summing, in (11), over a reduced set of operators that are a generator set, i.e. a minimal set of operators that when combined sufficient times generate the complete point group. This leads to a different symmetry measure, but one which will be zero for exact symmetry, but not otherwise. For any given group, there is typically more than one valid generator set, and each will give rise to a numerically different symmetry measure. For this reason, we use F 1 with a generator set for optimisation, but for comparison and testing of symmetry measures, F 0 with the full group is used.

Choice of coordinate frame: further detail
For some point groups, there are multiple feasible coordinate frame orientations for which a molecule with exact symmetry will give F = 0 . In these cases it is desirable to introduce further criteria that lead to an unambiguous standard orientation. They include Asymmetric tops-axis permutations: We first construct and diagonalise the second moment tensor Asymmetric top molecules give three distinct eigenvalues, and restrict the possible point groups to D 2h and its subgroups. The coordinate axes are defined by the eigenvectors for an exact structure, which form good starting guesses for minimising F otherwise. There are 6 possible coordinate axis permutations, of which 4 are infeasible for those groups with a unique C 2 axis or unique mirror plane. For D 2h and planar C 2v molecules, we follow Recommendations 5b, 5a respectively of Ref. [28], but otherwise any remaining freedom is satisfied by assigning the x axis to the minimum eigensolution. Symmetric tops-choice of perpendicular axes: Symmetric tops are characterised by double degeneracy in the second moment eigensolutions. The unique eigensolution maps to the z axis, and we follow Recommendations 5c and 5d of Ref. [28] where possible. Spherical tops: In spherical tops, all three eigenvalues of the second moment tensor are equal, and the eigenvectors offer no help in finding the orientation that (26) matches the molecule to the point group. Furthermore, the optimisation of F with respect to coordinate axes can be very ill-conditioned; for example, in the icosahedral groups, the lowest rank non-zero multipole moment has angular momentum 6, meaning that the system is very close to spherical. For these systems, we proceed by first determining the maximum-order rotational axis of the point group, and then looking for an approximate regular polygon of that order amongst the atoms. The normal vector of this polygon defines an axis which is then aligned to one of the point group's axes. Before entering BFGS optimisation, a discrete scan of rotations about that axis is performed, in order to find a starting guess with the lowest F . This procedure is relatively costly, but does succeed in aligning even icosahedral molecules such as C 60 .

Purification
Often, the application of this methodology will be to discover the point group, or to determine the extent of deviation of the structure from a given point group. But a further use is to identify a point group, and then refine the geometry so that it conforms as exactly as possible to the group. We will then seek the least-motion distortion of the coordinates that will result in zero F . We can proceed by minimising F with respect to the atomic coordinates using BFGS. The atomic-coordinate gradient of F is evaluated straightforwardly as where we define the inverse image map B B At ,t = A . However, this is an ill-posed problem, since when the symmetry conditions are obeyed exactly, any arbitrary step in the direction of any combination of coordinates that is a basis for a totally symmetric irreducible representation of the group will retain F = 0.
In principle, the arbitrariness can be overcome by always moving orthogonal to the null space, thereby satisfying, at least conceptually, an objective that might be expressed as finding the symmetrised structure that is as close as possible to the original. However, until the minimum is reached, the null space is not pure and exact, because of the slightly broken symmetry, and the steps in the numerical optimisation may introduce unnecessary additional motion in the symmetric directions.
Instead, a simple way to remove the arbitrariness is via a tie-breaking penalty function. We define a measure of displacement of the structure { ⃗ R A } from its unre- 1) , and measures the changes in all interatomic distances; f 0 is used to map the semi-infinite range to [0, 1) and has the effect of de-emphasising the weight of changes of the distances between very distant atoms. We then vary the coordinates to minimise F + P , where is a small chosen parameter.

Software implementation
All of the methodology described above is incorporated in a freely available software library [29]. The library is written in C++ with additional bindings for C (including Fortran-callable functions) and command line. Its principal functions are optimisation of coordinate frame for a given molecule and point group, calculation of the symmetry measure, point group discovery, and structure refinement. Table 1 illustrates the effect of fully optimising the coordinate frame to minimise the symmetry measure. Exact D 2h atomic coordinates for ethene are displaced randomly using a uniform distribution of specified width. The table shows, for a number of values of the noise parameter, the mean D 2h symmetry measures obtained by (a) not adjusting the coordinate frame; (b) adopting centre of charge and inertial axes; (c) full frame optimisation. A sufficiently large sample was taken to ensure the convergence of the means to the three significant figures quoted. It is seen that the use of inertial axes reduces the deviation by a factor of approximately 4, and that a further smaller, but significant, improvement results from full optimisation. Table 1 Mean D 2h symmetry measures for ethene contaminated with random noise in atomic coordinates. Each atomic coordinate is displaced by a random value drawn from a uniform distribution between plus and minus 'Noise'. 'Unoptimised' is the large-sample mean symmetry measure without any readjustment of coordinate frame. 'Inertial' and 'Optimised' give the mean symmetry measures after adjustment to centre of mass and inertial axes, and after full frame optimisation, respectively Noise Unoptimised Inertial Optimised

Conclusion
A new approach for fuzzy assignment of a point group to a molecule has been described. It produces the best match possible by choosing a coordinate frame that minimises a measure computed as a simple function of the molecular coordinates and point group specification. Except for improving algorithm speed and robustness in high symmetry cases, no inspection of the detail of the structure is needed to identify individual symmetry operations, with the consequence that the resources needed scale only linearly with the number of atoms and with the size of the group.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.