QSym2: A Quantum Symbolic Symmetry Analysis Program for Electronic Structure

Symmetry provides a powerful machinery to classify, interpret, and understand quantum-mechanical theories and results. However, most contemporary quantum chemistry packages lack the ability to handle degeneracy and symmetry breaking effects, especially in non-Abelian groups, and they are not able to characterize symmetry in the presence of external magnetic or electric fields. In this article, a program written in Rust entitled QSym2 that makes use of group and representation theories to provide symmetry analysis for a wide range of quantum-chemical calculations is introduced. With its ability to generate character tables symbolically on-the-fly and by making use of a generic symmetry-orbit-based representation analysis method formulated in this work, QSym2 is able to address all of these shortcomings. To illustrate these capabilities of QSym2, four sets of case studies are examined in detail in this article: (i) high-symmetry C84H64, C60, and B9– to demonstrate the analysis of degenerate molecular orbitals (MOs); (ii) octahedral Fe(CN)63– to demonstrate the analysis of symmetry-broken determinants and MOs; (iii) linear hydrogen fluoride in a magnetic field to demonstrate the analysis of magnetic symmetry; and (iv) equilateral H3+ to demonstrate the analysis of density symmetries.


Introduction
Symmetry provides a systematic framework to categorize and classify various mathematical quantities that are of interest to quantum chemists, such as electronic wavefunctions and densities, via the lenses of group and representation theories.The ability to examine these quantities based on their symmetry enhances one's arsenal of analysis tools that facilitate the assignment of such quantities calculated from approximate numerical methods to true eigenfunctions of the electronic Hamiltonian of the system.In such studies, having a robust method to unambiguously identify and label the symmetries of the quantities being investigated ensures that their properties can be correctly tracked and assigned to known or expected ground and excited electronic states of the system.This is especially true when the underlying equations that govern such quantities yield multiple solutions with differing degrees of physical relevance, thus making the task of understanding them much more challenging.][3][4][5][6] In both SCF HF and KS theories, spinorbitals are one-electron wavefunctions that form the cornerstones upon which relevant quantities of interest, e.g., single-determinantal wavefunctions in HF 7 and electron densities in KS, 8 are constructed.The spin-orbitals themselves have long been deemed to be of great importance, for they provide chemists with a useful means to interpret the underlying multielectron quantities which are often too complicated to examine directly.In fact, in HF theory, the spin-orbitals that result from the variational optimization of a single-determinantal ansatz form the starting point for many families of post-HF correlated methods such as configuration interaction (CI), 9 coupled cluster (CC), 10 and complete active space (CAS). 11,12On the other hand, in KS theory, there have long been discussions that the KS spin-orbitals are just as useful as their HF counterparts in chemical theories based on molecular-orbital (MO) models [see Refs.13,14 and also contributions (2.2.4)-(2.2.7) in Ref. 15].In either case, it is imperative that the shape and symmetry properties of spin-orbitals be identified so that they can be used effectively in the qualitative investigations of chemical phenomena 14 and the quantitative calculations of physical properties such as correlation energies (via post-HF correlated treatments), ionization potentials, 16 and vertical excitation energies. 17owever, it is not only the symmetry of spinorbitals that is important, since spin-orbitals are only one-electron functions and hence do not fully represent electronic states in any multi-electron system.In fact, in wavefunction theories, one often needs to obtain a good understanding of symmetry properties of multielectron wavefunctions before one can confidently attribute them to actual electronic states of the system, especially when one is interested in more than just the ground state, such as in the computation of electronic spectra. 18,19A few studies in which symmetry is used to assist the interpretation of ground-and excitedstate correlated wavefunctions can be found in Refs.20-22.In addition, the multiple, generally non-orthogonal, SCF solutions that arise from the HF equations may interact with each other in a CI expansion, if their symmetries are compatible, to give improved multideterminantal wavefunctions describing certain electronic states with definitive symmetries.Some examples of this include the examinations of low-lying HF solutions in the NO 2 radical, 23 in various classes of hydrocarbons, 24 in octahedral transition-metal complexes, 25 and in avoided crossings in LiF. 26 Furthermore, a thorough insight into the symmetry properties of wavefunctions and densities proves necessary to ensure formal correctness in the fundamental development of DFT and the interpretation of DFT calculation results.This is particularly important in degenerate systems where care must be taken to handle any symmetry breaking in the densities correctly to avoid the well-known symmetry dilemma that often arises in the KS formalism where the KS effective potential has a different symmetry from that of the physical external potential, as discussed in great detail by many authors including Görling, 27, 28 Savin, 29 and Chowdhury and Perdew. 30n addition, since chemistry is hardly ever static, it is often of great interest to follow electronic states as the symmetry of the system is varied.6][37] As the symmetry of the system changes, degeneracies might be lifted and broken symmetry (i.e., when a function and its symmetry partners span multiple irreducible representations of the full symmetry group of the system) might be restored. 25A knowledge of wavefunction and density symmetry allows one to correlate electronic states from low-symmetry configurations to high-symmetry configurations, thus gaining additional insight into their behaviors and properties.
Unfortunately, to the best of our knowledge, despite the importance of symmetry in quantum-chemical theory and computation, there does not yet exist any implementation for a general analysis of symmetry properties of electronic wavefunctions, densities, and potentials.In fact, many existing generalpurpose quantum chemistry packages such as Q-Chem, 38 Orca, 39 PySCF, 40 Dalton, 41 [Open]Molcas, 42 Psi4, 43 CFOUR, 44 and TURBOMOLE 45 come with features to carry out symmetry analysis to some extent, but most (with Q-Chem and TURBOMOLE being exceptions) opt to work in D 2h or one of its subgroups, all of which are Abelian groups whose irreducible representations are real and onedimensional, and hence are unable to take into account any spatial degeneracy in wavefunctions properly.Moreover, none of these packages is able to cope with symmetry breaking, nor are they programmed to examine symmetry properties of quantities other than wavefunctions, and as far as we are aware, no existing software provides options to analyze symmetry in the presence of external fields.
In this article, a framework for general symmetry analysis is introduced.This framework is implemented in a Rust 46,47 program named QSym 2 , which stands for Quantum Symbolic Symmetry, and which seeks to address some of the needs for symmetry in electronic-structure theory and computation that are currently not fulfilled by existing quantum chemistry packages.In particular, QSym 2 is designed to work with all possible finite point groups, Abelian or not, for which necessary character tables are automatically and symbolically generated on-the-fly, so that degeneracies and symmetry breaking can be represented correctly.In addition, this framework is sufficiently general to be applicable to any linear-space quantities and not just wavefunctions or densities.Furthermore, QSym 2 is capable of performing symmetry analysis in the presence of external fields, particularly ensuring that complex irreducible representations, which occur frequently when a magnetic field is present, are handled explicitly.In addition, QSym 2 is able to provide transformation matrices that enable the generation of symmetry-equivalent partners of any linear-space quantities, as long as they can be expanded in terms of atomic-orbital (AO) basis functions or products thereof.All of this is possible thanks to one governing design principle that QSym 2 undertakes, which insists that all of its computational elements (e.g., symmetry operations and irreducible representation characters) are treated symbolically as much as possible, so that defining properties of groups such as closure and existence of inverses are respected and utilized to guarantee accuracy and efficiency.
The article is organized as follows.In Section 2, the theoretical foundation for the symmetry analysis framework implemented in QSym 2 is laid out.In particular, the various aspects of group and representation theories involved in the determination of molecular symmetry groups, the management of symmetry operations, and the in situ generation of character tables are explained.This is followed by the formulation of a general method for representation symmetry analysis applicable to any linear space.Then, Section 3 presents several case studies to illustrate the usefulness of symmetry analysis via QSym 2 in interpreting and understanding electronic-structure calculations.Finally, Section 4 concludes the article with a few remarks on the capabilities and limitations of the symmetry analysis framework implemented in QSym 2 , and also charts possible directions for QSym 2 to be extended in the future.

Unitary symmetry of the electronic Hamiltonian
For a molecular system with N e electrons and N n nuclei in a uniform external electric field E and magnetic field B = ∇ × A(r), where A(r) denotes the magnetic vector potential, the electronic Hamiltonian is given by In atomic units, the first contribution has the form and is the zero-field Hamiltonian which has an explicit dependence on the multiplicative exter-nal potential v ext whose form is governed by the geometric arrangement of the nuclei, In Equation (3), r i denotes the position vector of the i th electron and R A that of the A th nucleus.The second contribution, describes the interaction between the electrons and the external electric field, 48 and the third contribution, where pi is the linear momentum operator for the i th electron, ŝi the spin angular momentum operator for the i th electron, and g s the electron spin g-factor, gives the interaction of the electrons with the external magnetic field. 49,50The unitary symmetry group G of the system consists of all unitary transformations û that leave Ĥ invariant: Clearly, G is the intersection of the unitary symmetry groups of Ĥ0 , Ĥelec , and Ĥmag , which we shall denote G 0 , G elec , and G mag , respectively.We further restrict the elements in these groups to be point transformations acting on the configuration space where physical systems such as atoms, molecules, and fields are described. 51hen, G 0 is also commonly known as the point group of the molecular system.A robust algorithm to determine the name and elements of G 0 for any molecular system has already been described by Beruski and Vidal. 52s shown formally in Appendix A of Ref. 34, the group G mag consists of orthogonal transformations in three dimensions [i.e., elements of the group O(3)] that would map the uniform magnetic field B onto itself and is commonly known as C ∞h , 32,53 which is an infinite Abelian group with principal axis parallel to B. A similar approach can be used to show formally that G elec consists of three-dimensional orthogonal transformations that would leave the uniform electric field E unchanged and is commonly recognised as C ∞v , 32 which is an infinite, but not Abelian, group with principal axis parallel to E. Hence, a naïve procedure to locate all elements of G is to first identify all elements of G 0 , and then to filter out only those elements that would leave E and/or B invariant.However, this procedure is unnecessarily wasteful as it requires additional efforts to be spent on finding a large number of elements of G 0 that would eventually be discarded, since the presence of external fields almost always leads to a reduction of unitary symmetry.In fact, for highly symmetric molecular systems where G 0 is large, these additional efforts can be non-trivial.

Including external fields: method of fictitious special atoms
It is desirable to make use of the algorithm by Beruski and Vidal 52 as much as possible to locate all elements of G directly without having to go through the intermediary of G 0 in the presence of external fields.To this end, we propose that fictitious special atoms be introduced to represent the external fields such that the combination of the molecule and fictitious atoms has the same unitary symmetry group G as the combination of the molecule and the external fields.Each fictitious special atom is characterized by a pair of parameters (t, R t ), where t encodes its type and R t denotes its position.A uniform electric field E is represented by one fictitious atom of type t = e placed at R e = R com +kE, where R com is the position vector of the center of mass of the molecule and k a scalar factor chosen to ensure that this fictitious atom does not coincide with any actual atom in the molecule, and that the subsequent unitary symmetry group determination is numerically stable.The vector R e − R com is therefore parallel to E, and as E is a polar vector, 54 it is im-posed that fictitious atoms of type e transform under all operations in the group O(3) just as any ordinary atom does.It is easily seen that the combination of the molecule and the fictitious atom has the same unitary symmetry group as the molecule in the external E field [Figure 1(a)].
On the other hand, a uniform magnetic field B is represented by two fictitious atoms, one of type b+ and the other of type b−, placed at R b± = R com ± kB, where the +/− signs in the type names signify the polarities of the fictitious atoms.The vector R b+ − R b− is parallel to B, and these two fictitious atoms transform under all operations in the group O(3) almost like any ordinary atom, but since B is an axial vector, 54 it is additionally required that the polarities of the fictitious atoms be reversed under improper transformations.This ensures that the combination of the molecule and the fictitious atoms has the same unitary symmetry group as the molecule in the external B field [Figure 1(b)].
With the introduction of fictitious special atoms, external fields are no longer required to be treated separately in the unitary symmetry group determination.In fact, fictitious atoms can be incorporated directly into the Beruski-Vidal algorithm, 52 provided that the following modifications are taken into account: (i) Fictitious atoms must be included in the calculation of the principal moments of inertia of the system and the subsequent classification into four main rotational symmetry types: spherical top, symmetric top, asymmetric top, and linear.For this purpose, a mass of 100.0 u is chosen for the fictitious atoms: there is no physical significance to this value; it simply has been found to ensure numerical stability in all test cases.
(ii) Fictitious atoms must be included in the determination of distance-based symmetrically-equivalent-atom (SEA) groups.This means that b+ and b− can be in the same SEA group if they both have the same distance signature to all other atoms in the molecule, despite their different polarities.
(iii) The possibility that polyhedral SEAs be arranged in a spherical top must also be taken into account.This additional possibility was not originally considered in Ref. 52 as it can only arise when a spherical top molecule is placed in an external magnetic field.Figure 2(a) shows an example where a magnetic field is applied along one of the C 3 axes of tetrahedral adamantane: this molecule-field combined system is now a symmetric top with the unique axis along the field direction, but the six carbon atoms highlighted in orange constitute a group of SEAs that are arranged in a regular octahedron.
(iv) The symmetric top rotational symmetry may also result in the C s group.This additional possibility was not considered in Ref. 52 either as it can only occur when an external field is applied to a spherical top in a manner that eliminates all symmetry elements of the system apart from a single mirror plane.Figure 2(b) illustrates an example of this when either a magnetic or an electric field is applied to tetrahedral CH 4 such that it is simultaneously parallel to one of the molecular mirror planes and perpendicular to another.

Magnetic symmetry of the electronic Hamiltonian
When antiunitary operations are taken into account, the unitary symmetry group G might no longer be the largest symmetry group of the electronic Hamiltonian Ĥ .In fact, in many studies involving magnetic phenomena and magnetic materials, [55][56][57][58][59] it is necessary to consider a supergroup of G that also contains antiunitary symmetry operations that leave Ĥ invariant.Such a group is called the magnetic symmetry group M of the system, and it can easily be seen 56 that M must admit G as a normal subgroup of index 2, so that we can write where â0 can be any of the antiunitary elements in M but must be fixed once chosen.The left  Let us now consider the time-reversal operation θ, which is an archetype of antiunitary operations (see Chapter 26 of Ref. 60 for an in-depth discussion of the time-reversal operation in quantum mechanics).It turns out that, with respect to θ, magnetic symmetry groups can be classified into just two kinds. 56,61,62The first kind are those that contain θ, in which case one can choose â0 = θ so that These are called magnetic grey groups.The second kind are those that do not contain θ; however, one can always find a unitary operation û0 not in the group such that the product θû 0 is an antiunitary operation that occurs in the group.This then enables one to write where â0 has been chosen to be θû 0 .Such groups are called magnetic black-and-white groups.It is then clear that, in the absence of an external magnetic field, θ is a symmetry operation of the system.However, this ceases to be the case when an external magnetic field is applied: the magnetic field vector B is timeodd 54,56 and thus gives rise to terms in the electronic Hamiltonian [Equation (5)] that do not commute with θ (see Appendix A of Ref. 34 for a detailed explanation).Therefore, the following general rules can be deduced: (i) in the absence of an external magnetic field, the system always has a magnetic symmetry group which must be one of the magnetic grey groups; (ii) in the presence of an external magnetic field, if the system exhibits any antiunitary symmetry, then it has a magnetic symmetry group that must be one of the magnetic black-and-white groups, but if the system exhibits no antiunitary symmetry, then it only has a unitary symmetry group.
For both kinds of magnetic groups, it is often useful to consider a unitary group M ′ that is isomorphic to M. In cases where M ′ is identifiable with a subgroup of the full rotationinversion group in three dimensions O(3) and can thus be given a Schönflies symbol, the magnetic group M can be written as M ′ (G). 59,62hen this is not possible, however, the antiunitary coset form with respect to the unitary symmetry group G and a representative antiunitary operation â0 [Equations ( 7)-( 9)] can always be employed to uniquely denote M because it is always possible to assign a Schönflies symbol to G, which is guaranteed to be a subgroup of the molecular point group G 0 (cf.Section 2.1.1). Figure 1 depicts two examples of how M is typically denoted.
To determine M and all of its elements given a molecular system in a uniform external field, the Beruski-Vidal algorithm 52 can once again be exploited with an additional modification that any unitary transformation considered in the algorithm can also be accompanied by the antiunitary action of time reversal.For all ordinary atoms and fictitious atoms of type e representing an applied electric field, time reversal has no effects.However, for fictitious atoms of types b+ and b− representing an applied magnetic field, their polarities must be reversed under time reversal due to the time-odd nature of the magnetic field vector B. 54,56 2.2 Abstract group construction

Computational representation of symmetry operations
In QSym 2 , symmetry operations located using the method described in the previous Section are stored as instances of the SymOp structure.It shall henceforth be written "SymOp(ĝ)" to denote an instance of the SymOp structure that represents the actual ĝ symmetry operation computationally.In order for this representation to be efficient and to respect discretegroup-theoretic properties, most notably compositability and closure, of the underlying symmetry operations, it is imposed that the SymOp structure fulfill the following traits: (i) equality comparisons that are equivalence relations -reflexivity, transitivity, and symmetry must be satisfied for the "=" relation between SymOp instances, which must take into account the 2π-periodicity of spatial rotations; (ii) hashability -each SymOp(ĝ) instance must be able to produce an integer hash value hash[SymOp(ĝ)] that allows itself to be looked up from a hash table with an average constant time O(1), and that must be compatible with equality comparisons: ) where " * " denotes the composition operation between two SymOp instances.
The design of the SymOp structure in QSym 2 is detailed in Section S1 of the Supporting Information to illustrate how the above traits are satisfied.

Unitary conjugacy class structure
Prior to generating the character table of the symmetry group, its conjugacy class structure must first be determined.The conjugacy class structure of a group in turn depends on how the conjugacy equivalence relation between group elements is defined.In this article, only the familiar unitary conjugacy equivalence relation is considered: which holds when all elements in the group are represented as mathematical unitary operators on linear spaces, even if some of them are actually physical antiunitary operators.A different conjugacy equivalence relation called magnetic conjugacy equivalence relation holds if some of the elements in the group are represented on linear spaces as mathematical antiunitary operators, 63 which leads to a different conjugacy class structure. 63,64Although magnetic conjugacy classes have also been implemented in QSym 2 , their uses in magnetic symmetry via corepresentation theory 60 will be examined in a future study.The classification of elements of finite molecular symmetry groups in QSym 2 is carried out via the Cayley table C of the group: which encodes the group's multiplicative structure in a two-dimensional array of integers.

Generation of character tables of irreducible representations
Once an abstract group structure has been obtained for the underlying symmetry group, its character table then needs to be computed to allow for subsequent symmetry analysis.This can indeed be performed on-the-fly in QSym 2 .Algorithms for the automatic generation of symbolic character tables [65][66][67] are well-known and have been implemented before, most notably in GAP. 68However, no such implementations exist for molecular symmetry applications in quantum chemistry.These algorithms are thus re-implemented in QSym 2 with additional functionalities to ensure that the generated character tables respect conventions that are familiar to most chemists, e.g., the labeling of irreducible representations using Mulliken symbols. 69The details of these algorithms are recapitulated in Section S2 of the Supporting Information.
The ability to generate character tables automatically and symbolically enables QSym 2 to completely circumvent the need for hard-coded character tables which would limit the number and types of groups available for symmetry analysis.In particular, as demonstrated in Section 3.1, QSym 2 is able to handle degeneracy in non-Abelian symmetry groups which are encountered in highly symmetric molecular structures such as disk-like boron clusters (D nh ), 70,71 hydro-clusters of group-14 elements in quantum dots (T d ), 72,73 and buckminsterfullerenes (I h ).QSym 2 is also capable of tackling complexvalued representations that frequently arise in the presence of an external magnetic field, e.g., eight out of the twelve one-dimensional irreducible representations in C 6h , which is the unitary symmetry group of benzene in the presence of a uniform perpendicular magnetic field, are complex.
In all cases, there is no requirement for the molecule and/or external fields to be in any predefined standard orientation for the charactertable generation algorithm to work.In fact, as long as the symmetry operations of the system can be computationally represented and composited as described in Section 2.2.1 and Section S1 of the Supporting Information, the group structure can be abstracted away from these concrete representations of symmetry operations to allow for the character table to be computed entirely algebraically without recourse to any other knowledge exterior to the group structure.Only when labels of computed irreducible representations are to be deduced is information about molecular structures and symmetry operation orientations required in order to satisfy Mulliken's conventions.This ensures that molecules and external fields can be placed in whatever orientation is the most sensible or convenient for chemical computation while still being able to benefit from the symmetry analysis offered by QSym 2 .

Representation analysis
An initial formulation of the method for representation analysis has been discussed by one of the authors in a previous article (Appendices B and C of Ref. 25).However, this formulation only focuses on wavefunctions in Hilbert spaces and therefore leaves out other quantumchemical quantities that are not wavefunctions but that still have symmetry properties.Examples of such quantities include electron densities, vibrational coordinates, and magnetically induced ring currents.In this Section, a more general formulation of this method will be presented in which all linear-space quantities are covered.It will also be pointed out how the computational availability of group multiplicative structures via Cayley tables leads to a reduction in the representation analysis time complexity from O(|G| 3 ) to O(|G|) by taking advantage of group closure.

Formulation of linear-space representation analysis
Characters and representation matrices.
Let V be a linear space and w an element in V whose symmetry under a prevailing group G is to be computationally determined.To this end, the linear subspace W ⊆ V that is spanned by the orbit where ĝi denotes the action of g i in V is first determined.The symmetry of w in G is then given by the decomposition of W into known irreducible representation spaces on V of the group G. Technically, ĝi and g i are two very different quantities: the former is an operator acting on V and thus a member of GL(V ) (i.e., the group of all general linear operators acting on V ), whereas the latter is a member of an abstract group G.However, this distinction is unnecessarily pedantic for the purpose of this article and will therefore be ignored: we will use the hatted forms almost exclusively, refer to them as members of the group G, and make no attempt to distinguish operators that represent actions of the same abstract element g i but on different linear spaces.
To characterize W , we seek its character function χ W whose value for each element ĝ in the group is given by where D W (ĝ) is the representation matrix of ĝ in some finite basis chosen for W .Let B = {e m : 1 ≤ m ≤ dim W } be such a basis.The elements of D W (ĝ) satisfy the set of equations one for each element e m in the basis.
Representation matrix determination.Equation ( 13) now needs to be solved in a suitably chosen basis to determine the diagonal D W nn (ĝ) elements so that the character value χ W (ĝ) can be computed.In principle, these equations can be viewed as a set of simultaneous equations that can be solved algebraically to give the required matrix elements.However, such an approach would be tedious and it is much more common for equations of this type to be solved using a projection operator, which requires the existence of an inner product.
Thus far, no reference has been made to any inner products on V , because inner products are not required in the definition of representation symmetry.In fact, symmetry is a linearspace property rather than an inner-productspace property.This realization has an important implication: we are at liberty to define any inner product that is the most convenient to compute for a given linear space V in order to construct a projection operator solely for the purpose of inverting Equation ( 13); the value of the character χ W (ĝ) must be independent of this choice of inner product, even if V itself does not possess an intrinsic inner product.
Let us now endow V with an inner product ⟨•|•⟩ with which the overlap matrix S between elements in the orbit G • w is defined: It would then be ideal to use the orbit G • w as a basis B for W with which Equation ( 13) can be solved to give the character values.How-ever, the elements in G • w are not necessarily linearly independent, and the matrix S is thus not necessarily of full rank.In this case, a tall rectangular matrix X can be constructed: where λ m is a non-zero eigenvalue of S and U im the i th component of the corresponding eigenvector.The matrix X allows a linearly independent basis for W to be defined: such that the overlap matrix in this basis, is of full rank and Equation ( 13) becomes where primed subscripts have been used for later convenience.If S is already of full rank, then we simply set B = G • w, and so S = S.
In either case, the square matrix S is invertible, with which a non-orthogonal projection operator Pm can be constructed: where S−1 mn = ( S−1 ) mn .This projection operator satisfies Applying Pm to both sides of Equation ( 17) and making use of Equation (19) gives Multiplying both sides by ⟨ wm | and using the definition of Pm in Equation ( 18) then yields or equivalently, by canceling out the ⟨ wm | wm ⟩ term on both sides, By reintroducing the original terms in the orbit G • w using Equation ( 16), we obtain S−1 mn X ♢ in ⟨ĝ i w|ĝ|ĝ j w⟩ X jm ′ , where The above result can be conveniently written in a matrix form: where which gives a closed-form expression for the representation matrix D W (ĝ) to be computed from elements in the orbit G • w. and In both cases, as long as the overlaps between the elements in the orbit G •w and the orbit origin w have been evaluated, which costs O(|G|) time, all matrix elements of S and T(ĝ) can be deduced without any further involvement of the expensive overlap computation.However, this optimization is only possible if one can make the identifications ĝk = ĝ−1 j ĝi [Equation (22)] and ĝl = ĝ−1 j ĝ−1 ĝi [Equation (23)] which require knowledge of the multiplicative structure of the group G.The implementation of SymOp in QSym 2 that enables the construction of the Cayley table [Equation (11)] achieves exactly this and allows QSym 2 to perform extremely efficient representation analysis of linear-space quantities.

Examples of linear-space representation analysis
Single-determinantal wavefunctions and spin-orbitals.Let us consider an N e -electron single-determinantal wavefunction: (24) In the above expression, Â = 1 N e !P ∈Sym(Ne) (−1) π( P ) P is the antisymmetrizer with P an element of Sym(N e ), the symmetric group of degree N e , and π( P ) the parity of P .The antisymmetrizer acts on the electronic spin-spatial coordinates x i in terms of which the spin-orbitals χ i are written.In these cases, the linear space V (cf.paragraph "Characters and representation matrices" under Section 2.4.1) is chosen to be an N e -particle Hilbert space denoted H Ne .It should be noted that the spin-orbitals χ i are special cases of Ψ det with N e = 1.Being a Hilbert space, H Ne comes equipped with the familiar inner product x (x 1 , . . ., x Ne ) dx 1 . . .dx Ne (25)   which can be leveraged to compute the orbit overlap matrix S [Equation ( 14)] for Ψ det under the action of the prevailing group G.In particular, for single determinants, the inner product in Equation ( 25) has a particularly simple form: 75 where χ w,i denotes the i th occupied spin-orbital of the Ψ det w determinant.Each spin-orbital can be expanded in terms of the AO basis functions according to where φ µ (x) is an AO spin-spatial basis function and µ a composite spin-spatial index.The required spin-orbital overlaps can then be written as where the two-center overlap integrals can be easily obtained from many available integral packages for Gaussian AO basis functions (e.g., Libint 76 and Libcint 77 ) or London AO basis functions (e.g., QUEST, 78 London, 79 BAGEL, 80,81 and ChronusQ 82 ).QSym 2 also implements its own generic n-center overlap integral routine based on the recursive algorithm by Honda et al. 83 that is capable of handling both Gaussian and London AO basis functions. 84he calculation of overlaps between single determinants [Equation ( 26 where the composite spin-spatial coordinate x 1 has been relabeled and separated into a spin coordinate s and a spatial coordinate r in the integrand.In an AO basis, ρ(r) can be expanded as where ϕ γ (r) and ϕ δ (r) are spatial AO basis functions, γ and δ spatial indices, and P δγ elements of the corresponding density matrix P in this basis.The containing linear space for ρ(r) is well known to be the Banach space X = L 3 (R 3 ) ∩ L 1 (R 3 ). 85Being a Banach space, X does not have an intrinsic inner product, but, as explained in Section 2.4.1, it is possible to endow X with an inner product for the purpose of representation analysis.The simplest such inner product can be defined as follows: It is straightforward to show that the above definition for ⟨•|•⟩ satisfies all required properties of an inner product: • conjugate symmetry: • linearity in the second argument: • positive-definiteness: where ⟨ρ|ρ⟩ = 0 if and only if ρ(r) = 0 identically, otherwise there would exist regions in R 3 where |ρ(r)| 2 < 0, which is not possible.
Using the basis-expanded form of the electron density in Equation ( 28) in the inner product definition in Equation (29) gives where are four-center overlap integrals computable using the generic n-center overlap integral routine as described earlier.
The calculation of overlaps between electron densities [Equation ( 29)] is available in QSym 2 , thus enabling the symmetry analysis of electron densities obtained from a wide range of electronic-structure methods from single-and multi-determinantal wavefunctions to DFT.

Results and Discussion
In this Section, several case studies showcasing the capabilities and utility of QSym 2 are presented.Each case study is based around a distinct computational chemical problem whose results can be better understood by a detailed analysis of symmetry provided by QSym 2 .

Degeneracy in non-Abelian groups
The first set of case studies consists of three molecules of various sizes and symmetries: tetrahedral C 84 H 64 quantum dot (Figure 3a), icosahedral C 60 (Figure 3b), and octagonal B 9 -(Figure 3c).The non-Abelian symmetry of these three molecules allows for degeneracy to occur in the SCF MOs that arise from either a HF or a KS-DFT description of the ground state of the system.By examining the symmetry of such degenerate MOs in the ground SCF solutions of these three molecules, we seek to demonstrate the capability of QSym 2 to determine degenerate irreducible representation la-bels accurately, irrespective of the size of the molecules or the complexity of the MOs.

Computational details
For each molecule, a ground-state KS-DFT calculation using an appropriate exchangecorrelation functional and basis set was performed in Q-Chem 6.1.0.In all calculations, suitable symmetry thresholds were chosen to ensure that Q-Chem produced symmetry assignments for the computed MOs in the highest possible group.Afterwards, all geometry information, basis set information, and MO coefficients from these calculations were passed to QSym 2 where the unitary symmetry group G of the system was deduced, following which the representations of G spanned by the MOs were identified and analyzed according to the formulation given in Section 2.4.The MO symmetry assignments from QSym 2 were then compared with those from Q-Chem.Q-Chem was chosen as the benchmarking program for these case studies because of its ability to assign degenerate symmetry labels in certain non-Abelian groups.As stated in Section 1, most other quantum-chemistry programs are only able to perform symmetry analysis in Abelian groups and are thus not suitable for this purpose.

Degenerate symmetry benchmarks
We begin with tetrahedral C 84 H 64 , which proves to be a straightforward case.Using the geometry reported by Karttunen et al., 72 a groundstate unrestricted CAM-B3LYP/6-31+G* calculation was performed and a set of KS MOs were obtained.Table 1 shows the symmetry assignments that have been produced by both Q-Chem and QSym 2 for the frontier MOs.These particular MOs have been chosen because they have been identified in Ref. 73 to be responsible for the most intense transition in the electronic absorption spectrum of this molecule, and are therefore the most interesting to examine from a symmetry perspective.For this system, Q-Chem is able to identify its symmetry group as T d , as is QSym 2 .Both programs are also able to agree on their symmetry assignments of the frontier MOs, be they degenerate or not, thus confirming that the representation analysis formulation implemented in QSym 2 (Section 2.4) is valid.We consider next icosahedral C 60 , which has a higher symmetry than tetrahedral C 84 H 64 , and for which an unrestricted CAM-B3LYP/6-31+G* calculation was also carried out to yield a set of KS MOs.It turns out that, even though Q-Chem is able to identify the symmetry group of this system as I h when the symmetry tolerance value is set at 1 × 10 −4 , it seeks recourse to C 5 (a subgroup of I h ) to perform symmetry analysis, but then fails to produce any symmetry assignments for almost all MOs except those that are non-degenerate.Tightening the symmetry tolerance value to 1 × 10 −5 forces Q-Chem to identify the symmetry group as C 2h instead, but this then allows for a successful assignment of symmetry labels to MOs, albeit only under C 2h .On the other hand, QSym 2 is able to both identify the symmetry group correctly as I h and classify the symmetry of MOs using the irreducible representations of I h .Table 2 shows these symmetry assignments for the highest three degenerate sets of occupied MOs.
An inspection of Table 2 4).In C 2h , the dipole moment operator transforms as A u ⊕ 2B u , and so the integrand transforms as  therefore all optically forbidden in I h .This example shows that, such MO symmetry considerations, when done accurately in the full symmetry group of the system, can improve the efficiency of algorithms that make use of symmetry to screen molecular integral evaluations [87][88][89][90] or sharpen the interpretation of spectroscopic results.2).This MO has A g (I h ) and A g (C 2h ) symmetry.
Finally, we consider an octagonal boron disc, B 9 -, which has D 8h symmetry and should, in principle, be simpler than the I h symmetry of C 60 .Unfortunately, no matter the symmetry tolerance value, Q-Chem is only able to determine the symmetry group of this anion as D 4h and classify the MOs calculated at the B3LYP/def2-TZVP level of theory using the irreducible representations of this group.QSym 2 , on the other hand, is capable of identifying the symmetry group as either D 4h or D 8h , depending on the choice of the distance threshold (Section S1.4 of the Supporting In-formation), and then assigning MO symmetry labels using the irreducible representations of the corresponding groups.These are summarized in Table 3 for the valence canonical MOs of B 9 -as reported by Ðorđević et al. 71 It can be seen that, in D 4h , the symmetry assignments from both Q-Chem and QSym 2 are in agreement, while those in D 8h computable only by QSym 2 provide a full symmetry classification of the MOs that is consistent with their nodal structures.

Symmetry breaking in SCF solutions and orbitals
The next case study seeks to demonstrate the ability of QSym 2 to detect and quantify symmetry breaking in SCF solutions and orbitals.For this purpose, an octahedral complex Fe(CN) 6 3has been chosen.Having an unpaired electron due to the low-spin d 5 configuration on the Fe 3+ metal center surrounded by an octahedral ligand field, this complex is expected to give rise to multiple low-lying SCF UHF solutions, some of which break spatial symmetry. 25

Computational details
In all calculations, the structure of Fe(CN) 6

-
was held fixed at an O h geometry with Fe−C = 2.025 602 3 Å and C−N = 1.157 074 6 Å.Multiple SCF solutions at the UHF/def2-TZVP level of theory were found for this geometry in Q-Chem 6.1.0using metadynamics 91 combined with the Direct Inversion in the Iterative Sub-Table 1: Comparison of symmetry assignments of frontier canonical MOs in C 84 H 64 calculated at the CAM-B3LYP/6-31+G* level of theory using the geometry reported by Karttunen et al. 72 For each MO χ(r), the isosurface is plotted at |χ(r)| = 0.008.

MO Isosurface
space (DIIS) algorithm. 92SCF convergence was set at a DIIS error value of 1 × 10 −10 as implemented in Q-Chem.For each converged solution, its geometry information, basis set information, and Pipek-Mezey-localized 93,94 MO coefficients were read in by QSym 2 from which symmetry assignments for the individual MOs, as well as for the overall wavefunction, were determined.Unfortunately, no benchmarking symmetry assignments were available because no existing programs were able to handle symmetry-broken quantities.
For each quantity that is symmetry-analyzed in QSym 2 , a threshold λ thresh S must be chosen to determine which eigenvalues of the orbit overlap matrix S [Equation ( 14)] are nonzero so that the transformation matrix X can be constructed [Equation (15)].Whether or not a particular choice of threshold is reasonable depends on the gap between the eigenvalue of S that is immediately above the threshold, λ > S , and the eigenvalue of S that is immediately below the threshold, λ < S .In all cases for Fe(CN) 6  3 -, the threshold was chosen such that log 10 λ > S − log 10 λ < S ≳ 4.

-
alongside their symmetry assignments from QSym 2 determined at the linear independence threshold λ thresh S = 1 × 10 −7 .All four solutions are found to be symmetry-broken (i.e., they each span a reducible representation space of O h ), and interestingly, the lowest two solutions, A and B, do not contain the T 2g term naïvely expected of the ground state for a d 5 configuration in an octahedral strong-field environment, 95 whereas the other two solutions, C and D, do.In addition, the A and B solutions have a different inversion symmetry (ungerade) compared to that of the C and D solutions (gerade).This strongly suggests that there is a qualitative difference between these two pairs of solutions.
For illustrative purposes, we focus only on the lower-energy solution in each of the two pairs, namely the A and C solutions.To acquire a Table 2: Comparison of symmetry assignments of frontier occupied canonical MOs in C 60 calculated at the CAM-B3LYP/6-31+G* level of theory using an I h -symmetrized geometry.For each MO χ(r), the isosurface is plotted at |χ(r)| = 0.008.The four-and five-dimensional irreducible representation labels follow the convention specified in Ref. 86.
Table 3: Comparison of symmetry assignments of valence canonical MOs in B 9 -calculated at the B3LYP/def2-TZVP level of theory using the geometry optimized at the same level by Ðorđević et al. 71 For each MO χ(r), the isosurface is plotted at |χ(r)| = 0.04.In QSym 2 , the distance thresholds yielding D 4h and D 8h are 10 −5 and 10 −4 , respectively (see Section S1.4 of the Supporting Information for an explanation of this threshold).
Table 4: Symmetry-broken SCF solutions of Fe(CN) 6 3calculated at the UHF/def2-TZVP level of theory using an O h geometry with Fe−C = 2.025 602 3 Å and C−N = 1.157 074 6 Å.Each solution has M S = +1/2.All symmetries were determined using QSym 2 with a linear independence cut-off λ thresh S = 1 × 10 −7 .See Section 3.2.1 in the main text for the description of λ > S and λ < S .
(a) Symmetries of four lowest-lying M S = +1/2 UHF solutions located in Q-Chem using SCF metadynamics and DIIS with a convergence threshold of 1 × 10 −10 .Solutions are labeled alphabetically in ascending order of energy.crude understanding of the origin of this qualitative difference, we turn to the Pipek-Mezeylocalized MOs 93,94 obtainable from the canonical MOs of these two solutions, mainly because localized MOs have been known to provide a useful link between detailed quantumchemical calculations and classical chemical concepts such as non-bonding orbitals, lone pairs, and multiple bonds with which most chemists have gained great familiarity and intuition. 96,97In the particular case of Fe(CN) 6 3 -, localized orbitals help quantify the number of delectrons on the iron center and allow for a discussion on the nature of the UHF solutions in terms of the metal d n -electronic-configuration and oxidation-state descriptors that are common in coordination chemistry.

SCF
In Tables 4b and 4c, the Pipek-Mezeylocalized d-MOs for the A and C solutions are listed, respectively.Each of these MOs has a dshell Mulliken population of at least 0.980 and can therefore be regarded as being predominantly contributed to by a d-electron on the iron center.Clearly, the iron center in solution A admits a d 6 configuration, whereas that in solution C admits a d 5 configuration.This is also confirmed by a LOBA oxidation-state analysis formulated by Thom et al.: 97 the iron center in solution A has an oxidation state of +2, whereas that in solution C has an oxidation state of +3.Noting that all d-orbitals on the iron center must have gerade inversion symmetry, and also that all p-orbitals on the iron center have been confirmed to be paired, we conclude that the ungerade inversion symmetry in solution A must arise from an unmatched ungerade ligand orbital between the two spin spaces.The seemingly innocent ungerade inversion symmetry found in solution A turns out to be a manifestation of a ligand-to-metal charge transfer process.
The fact that the four UHF solutions A-D are symmetry-broken means that none of them is able to provide a physical description of the ground state of the system. 98However, it has been demonstrated elsewhere 25 that post-HF methods such as NOCI can yield multi-determinantal wavefunctions that conserve symmetry and thus provide more ap-propriate approximations of the ground state.For this to be viable, either a basis B spanning a complete representation space W [Equation (16)] or a full symmetry-equivalent orbit spanning the same space [Equation (12)] must be provided as the basis for NOCI -both of which can be generated by QSym 2 .
We conclude this case study with a remark that, as expected, the symmetry breaking of the overall determinants shown in Table 4a can be traced back to the symmetry breaking of the constituting orbitals.This is in fact demonstrated by the symmetry assignments for the Pipek-Mezey-localized d-MOs of the A and C solutions in Tables 4b and 4c, respectively.It should be noted, however, that symmetry breaking effects can sometimes be subtle and difficult to discern from a mere visual inspection of isosurface plots.A detailed analysis based on the formulation in Section 2.4 should therefore be preferred to obtain reliable symmetry information.For instance, consider the MO χ β 34 of solution A whose isosurface at 0.014 is shown in Table 4b.At first glance, this orbital appears just like a typical d yz orbital (with some distortions due to interactions with the ligands) and should just have T 2g symmetry.However, a close inspection of this isosurface, with the aid of the contour plot in the yz-plane shown in Figure 5, reveals that the interactions with the ligands on the y-axis are not equivalent to those on the z-axis, thus causing the relation Ĉx 4 χ β 34 = −χ β 34 to fail to hold.In other words, χ β 34 and Ĉx 4 χ β 34 are linearly independent, which gives rise to the T 1g ⊕ T 2g symmetry breaking as observed.The large gap between the boundary orbit overlap eigenvalues λ > S and λ < S (ca. 10 orders of magnitude) indicates that this symmetry breaking is in fact not just a numerical artifact of the analysis, but rather an intrinsic feature of this MO.

Symmetry in external fields
Thus far, we have demonstrated the symmetry analysis capability of QSym 2 for real orbitals and determinants in the absence of any external fields.In this next case study, we illustrate how QSym 2 can be used to understand quantum- chemical behaviors when external fields are introduced.In particular, we shall show that, for a hydrogen fluoride molecule in a uniform magnetic field, a knowledge of the symmetries of the complex-valued MOs helps rationalize the reversal of the electric dipole moment along the inter-nuclear axis observed by Irons et al. 34 at strong fields perpendicular to the molecule but not, curiously, at parallel fields.

Computational details
We followed Irons et al. 34 and performed current-DFT calculations with the cTPSS functional in the uncontracted aug-cc-pVQZ basis set 99,100 employing the resolution-of-theidentity approximation with the AutoAux auxiliary basis 101 in QUEST. 78The obtained KS MOs were then passed to QSym 2 for symmetry analysis in the appropriate unitary symmetry group G of the molecule-plus-field system (cf.Section 2.1.3).In all cases, the gauge origin of the magnetic field and the center of mass of the molecule were placed at the origin of the Cartesian coordinate system so that the gauge origin would always be left invariant by the applications of all symmetry operations during the orbit generation [Equation ( 12)]. 102r the parallel and perpendicular field orientations, complex MO isosurfaces were also plotted in VMD 103 using the method described by Al-Saadon et al. 104 In the cases where G is an infinite group (i.e., C ∞v at zero field or C ∞ at parallel-field orientations), a finite integer order n ≥ 2 is chosen for the infinite-order rotation axis C ∞ so that G is restricted to a suitable finite subgroup G n (i.e., C nv or C n , respectively) in which the representation analysis of Section 2.4 is carried out by QSym 2 .The actual representations in G can then be deduced by the representations in G n produced by QSym 2 according to the following subduction rules: where Γ k and Γk in C ∞ and C n are complexconjugate one-dimensional irreducible represen-tations with character functions For   Figure 6) at various strengths and orientations of the external magnetic field together with the values of the electric dipole moment component along the inter-nuclear axis, p z .The H -F bond length is kept fixed at its zero-field equilibrium value, 0.928 97 Å, for all field strengths and orientations.It can be seen that, in the region where p z becomes less negative and approaches zero (|B| ≥ 0.5B 0 ≈ 1.18 × 10 5 T and ϕ ≈ 90 • ), the frontier MO landscapes display significant curvature.This suggests that these MOs interact with one another strongly in this region, and these interactions might be responsible for the observed electric dipole reversal.However, in parallel fields, even at very high field strengths (|B| ≥ 0.7B 0 ≈ 1.65 × 10 5 T), p z remains at ca. −0.7 a.u.which is approximately the same value as that at zero field.The energy landscapes of the frontier MOs also show very little curvature in the vicinity of ϕ = 0 • or 180 • , thus implying a lack of interaction between these MOs and further strengthening the conjecture that these MOs must interact in some way to result in a reversal of the electric dipole moment.
In Figure 7b, cross-sections through the frontier MO energy landscapes and p z plots at ϕ = 0 • and 90 • are shown together with MO symmetry assignments from QSym 2 and complex isosurface plots as described earlier.It becomes immediately obvious that, at ϕ = 0 • , the three frontier MOs have different symmetries at all values of |B| and are therefore unable to interact via the KS operator.In fact, the Σ and Γ 1 energy curves are able to cross because of their different symmetries.Even though their energies vary quite significantly as |B| increases, this variation is due only to the interactions of these MOs with the applied field.Although these interactions do lead to qualitative changes in the shapes of the MOs, most notably the disappearance of nodal planes in the two Π MOs at zero field that become Γ 1 / Γ1 MOs at |B| > 0, these changes do not actually affect the distribution of electrons along the inter-nuclear axis in any significant way.This is indeed confirmed by the near equality of the ⟨χ|μ z |χ⟩ values of these MOs at |B| = 0 and 0.74B 0 .Consequently, the electric dipole moment along the inter-nuclear axis remains almost unchanged.
The situation is markedly different at ϕ = 90 • .The three frontier MOs now have the same symmetry in C s and are thus permitted to interact via the KS operator.Indeed they do, as  is evident from the distortions in their energy curves for |B| ≥ 0.5B 0 and also in the shapes of their isosurfaces.Most significantly, the highest occupied MO (HOMO) shows the most drastic change from a 2p orbital localized entirely on the fluorine atom to a laterally delocalized MO with a pronounced lobe on the hydrogen atom.Associated with this change is the large increase in the value of ⟨χ|μ z |χ⟩ for this MO from −1.680 a.u. at zero field to −1.047 a.u. at |B| = 0.74B 0 , which more than outweighs the decreases in the values of ⟨χ|μ z |χ⟩ for the other two frontier MOs.There is thus a partial charge transfer from the fluorine atom to the hydrogen atom in the HOMO induced by the perpendicular magnetic field, which is responsible for the observed dipole reversal.

Symmetry of electron densities
In the final set of case studies, we demonstrate the ability of QSym 2 to perform symmetry analysis for electron densities, as formulated in Section 2.4.2.In particular, we show how the symmetry of the electron density is intimately related to that of the underlying electronic wavefunction, both at zero field and in the presence of external electric and magnetic fields.

Computational details
For the above purpose, we chose the equilateral geometry of H 3 + that has been found in Ref. 37 to be the optimal geometry for the lowest M S = −1 electronic state when a uniform magnetic field of strength |B| = 1.0B 0 is applied perpendicular to the plane of the molecule.For this geometry, the lowest M S = −1 wavefunctions and densities were computed in QUEST 78 at the UHF/6-311++(2+,2+)G** level of theory in three cases: at zero field, in the presence of a perpendicular uniform electric field with strength |E| = 0.1 a.u., and in the presence of a perpendicular uniform magnetic field with strength |B| = 1.0B 0 .The symmetry assignments for the resulting determinantal wavefunctions and the corresponding electron densities were then determined by QSym 2 .In all cases, the H 3 + structure was placed in the yz-plane so that any external field applied perpendicular to the molecule would be along the x-direction.

Density symmetries in H 3 +
Table 5 shows the wavefunction and density symmetries of the lowest M S = −1 UHF wavefunction in the three cases described above.We examine first the perpendicular magnetic field case (labeled B x in Table 5) where the unitary symmetry group of the molecule-plusfield system is C 3h .The lowest M S = −1 UHF wavefunction has already been reported in Ref. 37 to have Γ ′ (C 3h ) symmetry, which is a one-dimensional irreducible representation in C 3h whose character function satisfies χ Γ ′ ( Ĉ3 ) = exp(2iπ/3) and χ Γ ′ (σ h ) = 1.As this is a non-degenerate wavefunction, the corresponding density must be totally symmetric in C 3h , which is indeed the case verified by the density symmetry assignment and also by the density isosurface and contours in the yz-plane.
Here, the electron cloud can be seen to be equidistributed over the three symmetry-equivalent hydrogen nuclei.
The situation is rather different in the other two cases.In the absence of any external fields (labeled 0 in Table 5), for which the unitary symmetry group is D 3h , the lowest M S = −1 UHF wavefunction at the same geometry turns out to exhibit symmetry breaking due to its A ′ 1 ⊕ E ′ (D 3h ) symmetry.The symmetry analysis of the corresponding total density shows that the density also exhibits an A ′ 1 ⊕ E ′ (D 3h ) broken symmetry.Similarly, in the presence of a perpendicular electric field (labeled E x in Table 5), for which the unitary symmetry group is reduced to C 3v , the UHF wavefunction now has A 1 ⊕ E(C 3v ) symmetry, as does the corresponding density.The symmetry breaking of the total density in these two cases can be visualized most easily via the density isosurfaces and yz-contours: the electron cloud is not equally distributed over the three symmetry-equivalent hydrogen nuclei.
The external fields can also be applied parallel to the plane of the molecular frame of H 3 + .The  wavefunction and density symmetries and density isosurfaces resulting from the UHF calculations in these field orientations are summarized in Table S1 in the Supporting Information.In all cases, it can be observed that, whenever the wavefunction is non-degenerate, the density is totally symmetric in the prevailing symmetry group, and whenever the wavefunction exhibits degeneracy, be it because of symmetry breaking or not, the corresponding density becomes non-totally-symmetric.
We note that it is also possible to obtain density symmetries without having access to any symmetries of the underlying wavefunction, meaning that the same approach can be readily applied to densities from correlated wavefunction methods as well as HF or KS-DFT, with no additional implementation.In the context of KS-DFT, this means that symmetry analysis can be carried out directly on the density, rather than the KS orbitals and non-interacting wavefunction as a proxy for the physical wavefunction.
Examples of symmetry analyses for KS densities obtained with the r 2 SCAN0 functional are also shown in Table S1 in the Supporting Information for various external field orientations.By considering the symmetries of electron densities, we have a qualitative way to compare HF and KS-DFT calculations: if HF and KS densities have the same symmetry, then there is a likelihood that both calculations describe the same electronic state of the system, but if HF and KS densities differ in their symmetries, then it must be concluded that they are qualitatively different, perhaps because they describe different electronic states, which can occur especially if there are multiple SCF solutions that occur close to one another (see Refs. 3-5,25,91 and also Section 3.2 above).In fact, for the cases listed in Table S1, based on the symmetries of densities, both HF and KS-DFT calculations in each field orientation can be said to approximate the same electronic state.

Conclusion
A new program for quantum symbolic symmetry analysis, QSym 2 , has been presented in this work.A key feature of the program is its capability to generate character tables symbolically on-the-fly, which endows it with the ability to perform symmetry analysis for general groups automatically.This flexibility means that QSym 2 can yield reliable symmetry assignments for systems exhibiting degeneracy and symmetry breaking effects, where standard implementations cannot be applied.In addition, QSym 2 handle reduced symmetries arising in electric and/or magnetic fields, thus providing a valuable tool for analysis and insight into systems under less chemically intuitive conditions.
The ability of QSym 2 to perform analysis of high-symmetry systems was demonstrated for C 84 H 64 , C 60 , and B 9 -, where in each case the full molecular symmetry group could be correctly identified and carried through to classify the symmetry of the resulting MOs, including their associated degeneracies.The octahedral transition-metal complex Fe(CN) 6 3was then used to demonstrate how QSym 2 is able to deduce representation spaces spanned by symmetry-broken determinants and MOs, giving a way to classify and understand symmetry breaking effects.Furthermore, the capability of QSym 2 to analyze symmetry in external magnetic fields was demonstrated for the hydrogen fluoride molecule, where the symmetry of the MOs under a magnetic field was shown to provide a rationalization of the behavior of the molecular electric dipole moment as a function of field strength and orientation.
An important benefit of the generic symmetry-orbit-based representation analysis framework formulated in this article and used in QSym 2 is the ability to analyze symmetry of quantities other than wavefunctions and MOs that arise in quantum-chemical calculations.As a simple example, the changes in the electron density of the equilateral H 3 + as a function of electric and magnetic field were analyzed, with the density symmetry analysis revealing the symmetry breaking or conservation of the un-derlying wavefunction.This approach can be applied on an equal footing to densities arising from SCF calculations or more elaborate post-HF correlated calculations without the need to explicitly perform symmetry analysis on the correlated wavefunctions.
The generality in the code design of QSym 2 opens up many avenues for future research.In particular, the applicability of the symmetryorbit-based representation analysis to members of general linear spaces makes it possible to directly consider the symmetry of many important quantities in quantum chemistry.One group of such quantities includes normal coordinates that describe translational, rotational, and vibrational modes of molecules.Another interesting class of such quantities includes functions of electron density and/or density matrix, of which the Fukui function, 105 which encapsulates chemical reactivity information, and the magnetically induced current density, 106,107 which provides an interpretation for observations in magnetic spectroscopic methods, are prime examples.Moreover, QSym 2 can already provide a more complete analysis of magnetic symmetry using corepresentation theory 56,60,62 and is not limited to uniform external fieldsthese developments will be reported in future publications.
Finally, we emphasize that the Rust implementation of QSym 2 is also flexible, such that the program can operate either as a tool to be applied subsequent to a quantum-chemical calculation in a stand-alone manner (as used in this work with Q-Chem), or as a library readily integrated into existing programs (as used in this and earlier 37,108 work with QUEST 78 ).Since the program is available open-source, we hope that it will become a useful tool for application in a wide range of chemical simulations.(84) Strictly speaking, the integrals obtained from these packages only give the spatial part of the spin-spatial integral in Equation (27).The spin part is typically handled separately and implicitly, especially for the commonly used spin functions α and β whose orthogonality is known.

Figure 1 :
Figure 1: Equivalence between systems in external fields and systems with fictitious special atoms.(a) A single fictitious special atom of type e is placed at R com + kE to represent a uniform electric field.(b) Two fictitious special atoms, one of type b+ and the other of type b−, are placed at R com ± kB to represent a uniform magnetic field.

Figure 2 :
Figure 2: Two special cases involving a uniform external field where the original Beruski-Vidal algorithm 52 needs to be modified.(a) A tetrahedral adamantane molecule placed in a uniform external magnetic field oriented along one of its C 3 axes.This illustrates a possible scenario in which a polyhedral SEA group (the six carbon atoms highlighted in orange) is arranged in a spherical top fashion (a regular octahedron).(b) A tetrahedral methane molecule placed in a uniform external magnetic or electric field oriented simultaneously parallel to one of the molecular mirror planes and perpendicular to another.This illustrates a possible scenario of the C s unitary group arising from the symmetric top rotational symmetry.

Figure 5 :
Figure 5: Contours of Pipek-Mezey-localized MO χ β 34 of solution A in the yz-plane.The inequivalence between the interactions of the (CN) -π orbitals in the y-and z-directions with the Fe-based d yz orbital accounts for the T 1g ⊕ T 2g symmetry breaking of this MO.

Figure 6 :
Figure 6: A simplistic depiction of the MOs in hydrogen fluoride at zero field.

3. 3 . 2
Figure7ashows the landscapes of the m s = +1/2 frontier MOs in hydrogen fluoride (cf.Figure6) at various strengths and orientations of the external magnetic field together with the values of the electric dipole moment component along the inter-nuclear axis, p z .The H -F bond

Figure 7 :
Figure 7: (a) Energy landscapes of the frontier m s = +1/2 MOs in hydrogen fluoride at various magnetic field strengths and angles.(b) The cross-sections through these landscapes at parallel (top) and perpendicular (bottom) field orientations.Annotated on these cross-sections are the MO symmetries and isosurfaces plotted at |χ(r)| = 0.100.The color at each point r on an isosurface indicates the value of arg χ(r) at that point according to the accompanying color wheel.The numerical value next to each isosurface gives the value of the orbital electronic dipole moment ⟨χ|μ z |χ⟩ for the associated MO.In C ∞ , the one-dimensional irreducible representation Γ 1 has character function χ Γ 1 [ Ĉ∞ (ϕ)] = exp(iϕ), and the corresponding complex-conjugate one-dimensional irreducible representation Γ1 has character function χ Γ1 [ Ĉ∞ (ϕ)] = exp(−iϕ), where Ĉ∞ (ϕ) denotes an anticlockwise rotation through an angle ϕ as viewed down the z-axis.
The compositions ĝi ĝj are effected computationally through the corresponding compositions SymOp(ĝ i ) * SymOp(ĝ j ) of the SymOp structure.Once the Cayley table C has been computed and stored, any operations that call for the multiplicative structure of the group, such as the determination of the conjugacy class structure or the construction of the group's character table (Section S2 of the Supporting Information), only need to make cheap queries to C without having to repeatedly recalculate group element compositions.
Optimization by group closure.It is clear from Equation(20)that the computation speed of D W (ĝ) is limited by the computation speed of the orbit overlap matrix S and the matrices T(ĝ), for all ĝ ∈ G. Naïvely, explicit constructions of S based on Equation (14) and of T(ĝ) based on Equation (21) would incur time complexities of O(|G| 2 ) and O(|G| 3 ), respectively.However, closure of G under composition allows all matrix elements of S and T(ĝ) to be identifiable with only |G| unique values:

Table 2 (
Figure each choice of integer n ≥ 2, irreducible representations up to and including E ⌈n/2⌉−1 in G = C ∞v and Γ ⌈n/2⌉−1 and Γ⌈n/2⌉−1 in G = C ∞ remain irreducible in the respective subgroups G n = C nv and C n .These irreducible representations in G can therefore be deduced unambiguously from those in G n .In the current case study of hydrogen fluoride, it is known from basic MO theory (Figure6) that MOs of up to Π symmetry at zero field are occupied in the ground state, and so we require n ≥ 3 so that C nv and C n have enough irreducible representations to describe Π(C ∞v ), Γ 1 (C ∞ ), and Γ1 (C ∞ ) symmetries unequivocally.In fact, for good measure, we chose n = 8 in all infinitegroup symmetry analyses for hydrogen fluoride in QSym 2 .

Table 5 :
Electronic wavefunction and total density symmetry of the lowest M S = −1 state of H 3 + in the presence of external electric and magnetic fields.Calculations were performed at the UHF/6-311++(2+,2+)G** level of theory.The magnitude of the applied electric field is 0.1 a.u. and that of the applied magnetic field is 1.0B 0 .Isosurfaces of total densities are plotted at |ρ(r)| = 0.050.
Table of Contents Graphic.