University of Birmingham Generating symmetry-adapted bases for non-Abelian point groups to be used in vibronic coupling Hamiltonians

The vibronic coupling Hamiltonian is a standard model used to describe the potential energy surfaces of systems in which non-adiabatic coupling is a key feature. This includes Jahn–Teller and Renner–Teller systems. The model approximates diabatic potential energy functions as polynomials expanded about a point of high symmetry. One must ensure the model Hamiltonian belongs to the totally symmetric irreducible representation of this point group. Here, a simple approach is presented to generate functions that form a basis for totally symmetric irreducible representations of non-Abelian groups and apply it to D 1 h (2D) and O (3D). For the O group, the use of a well known basis-generating operator is also required. The functions generated for D 1 h are then used to construct a ten state, four coordinate model of acetylene. The calculated absorption spectrum is compared to the experimental spectrum to serve as a validation of the approach.


Introduction
For decades, vibronic coupling models [1][2][3][4] have served as bridges connecting nuclear dynamics studies with the static studies of electronic structure calculations [5]. The vibronic coupling model is a simple polynomial expansion of diabatic potential energy surfaces and couplings. The expansion coefficients are chosen so that the eigenvalues of the potential operator map on to the adiabatic potential surfaces. This diabatisation by ansatz circumvents many of the problems of describing non-adiabatic systems. It is also the inspiration for a diabatisation scheme that is used in modern, direct-dynamic methods that include non-adiabatic effects [6]. For a model Hamiltonian to correctly approximate the eigenvectors of the true Hamiltonian it has to span the totally symmetric irreducible representation (IrRep) of the point groups the molecule belongs to, at the appropriate symmetric geometries [7]. In recent times, many articles have demonstrated the advantages of using symmetry when constructing analytic model potentials [8][9][10][11][12], most often in the context of permutation-inversion groups [13].
Vibronic coupling Hamiltonians predominantly use a Taylor expansion of the nuclear coordinates that suitably represent the (quasi)-diabatic electronic potential operator elements around the point of high symmetry. Quasi-diabatic states conserve the property of having slowly varying potential energy [14], thereby requiring few, low-order polynomial terms to converge its Taylor series in some pertinent region of interest.
Generally, electronic excited states at point-symmetric nuclear configurations will form a bases for an IrRep of the point group the molecule belongs to. By choosing coordinates that also form a basis for IrReps, symmetry allows us to determine whether a given monomial term is allowed in some element or whether monomials across different elements share coefficients. A textbook example are the linear terms in the E e Jahn-Teller diabatic model describing E degenerate states with a branching space along e degenerate modes. The symmetry of this system dictates that the linear coupling and gradient should share the same coefficient, correctly resulting in the well-known ''mexican hat" adiabatic potential. When constructing such models, one must always ensure such relationships are maintained, since these ensure the symmetry of the system is kept. It is therefore most convenient to work with a symmetry-adapted basis of matrices which obey the desired symmetry considerations. This is especially relevant in the case of non-Abelian groups, since it is there we find situations like the one just described.
In this paper, we present approaches for generating such diabatic matrices that form bases for totally symmetric IrReps, starting from functions representing electronic states and nuclear coordinates which transform as known IrReps of the group for which the matrix representations for the operations are known. We generate matrices for the O and D 1h point groups and provide all the polynomial expansions representing nuclear coordinates (up to fourth and sixth order respectively) and states for many of the IrReps of each group.
Molecules that belong to either of these groups will exhibit Jahn-Teller and Renner-Teller effects respectively and are archetypal examples of the breakdown of the Born-Oppenheimer approximation. The Jahn-Teller (JT) effect occurs when a molecule in a high-symmetry configuration has degenerate electronic states. Vibronic coupling between these states along degenerate nuclear coordinates form a conical intersection at the high symmetry point, stabilising the system away from it and thus lowering the symmetry of the system. The Renner-Teller (RT) effect is a special case that happens in linear molecules. Due to symmetry considerations a glancing intersection rather than a conical intersection is formed [4].
Cederbaum et al. were the first to study the form of vibronic Hamiltonian of linear systems for interacting degenerate and non-degenerate electronic states to different orders and their effects on resulting spectra [1,15]. To generate the D 1h Hamiltonian, we followed a similar approach to Viel and Eisfeld [16], who generated matrices up to sixth order for E e Jahn-Teller Hamiltonians. They used them to fit to the 2 E 0 anharmonic surfaces along twofold e 0 stretches of NO 3 obtained from ab initio MR-SDCI calculations. They compared results of several diabatic potential models, one of which includes extra ad hoc functions that do not obey the C 3 requirements. They warn over the use of such ad hoc functions, showing they can result in wavefunctions with different expectation values, population transfer and autocorrelation values to those models obeying the C 3 requirements. That work was followed by the generation of bases properly describing pseudo Jahn-Teller (pJT) coupling between E degenerate states to nondegenerate A states [11]. They tested them on the pJT coupling between the ground state 2 A 00 2 and exited 2 E 0 states along the twofold e bending mode surfaces of cation NH 3 + obtained from ab initio MRCI calculations. A large decrease in fitting error was obtained as higher order terms were included. For the O group, although the form of the 3D Jahn-Teller Hamiltonian had been known since the first half of last century [17] and later for all IrReps of the group [18], it has been only recently that Opalka and Domcke used invariance theory to generate higherorder expansions as bases for the E and T 2 IrReps [12,19]. Here we used a basis generating operator, devised by Wigner, Lowdin and Shapiro, amongst others [20] to generate Hamiltonians for degenerate states and coordinates of A 1 ; E and T 1 symmetry. This is a well known method for generating symmetry-adapted functions and only briefly presented here.

Basis generating operator
A given function from a set of n v orthonormal basis functions classes, under a given transformation operation O R corresponding to the symmetry operator R of the group, must satisfy: owing to the great orthogonality theorem: where g is the order of the group. In other words, operating with P R D l ij ðRÞ Á O R on a basis function f v q will generate another basis function f v i of the same IrRep or else annihilate the function (if vl).
A second operator can be constructed by choosing i ¼ j and summing over j: This operator P l has the property of annihilating any function that does not belong to the lth IrRep space, or else project out any function which does. With both these operators it is therefore possible, starting with an arbitrary function within the lth space (called a generator of expansion), to generate the complete set of orthonormal functions belonging to this IrRep. To generate the totally symmetric IrReps of the group, it is therefore only necessary to utilise the latter operator.
To generate all polynomials described here we used the opensource mathematics software SAGE [21].

D 1h Renner-Teller symmetry-adapted basis
There are an infinite number of possible gerade/ungerade E n IrReps indexed by n where a ¼ p n ; a being the angle of the C a rotation required to interchange any two real basis functions forming a twofold degenerate representation of this group.
The dominant contribution of the first few singlet excited states of organic molecules tend to be from functions with low orbital angular momentum. For linear organic molecules we therefore expect that the low singlet states, being formed from p and r functions built of l ¼ 2p functions would result in electronic states forming a basis for low n representations. Similarly, it is rare to find a full orthonormal set of coordinates with high n index number for small molecules (for example, normal coordinates do not exceed P symmetry [7]).
Given these considerations, we will restrict the construction of symmetry-adapted bases of matrices by solely using functions that form a basis for R þ=À ; P u=g and D u=g IrReps (A n and E n for n ¼ 1; 2) to represent diabatic states and only P u=g functions to represent nuclear coordinates.
The resulting matrices were used to construct a 10-state 4dimensional diabatic model in the subspace of P u=g coordinates exhibiting Renner-Teller coupling. Theoretical absorption spectra were calculated which compare well to experimental one, validating the use of the bases generated (Section 3). Herzig and Altmann have published the most comprehensive book of point group tables to date [22]. One can use spherical harmonic functions as a basis to form a representation of any point group. They provide the matrix representations of the P g=u and D u=g IrReps for all the operations of the group, in the complex symmetrised spherical harmonic basis: where Y m l are the ortho-normalised spherical harmonics using the Condon-Shortley phase convention [22] (the D 1h tables can be found in Herzig and Altmann and are reprinted in Appendix A). The effect of operating on singly degenerate IrReps (R þ=À ) can be obtained from character tables. The corresponding real Cartesian basis functions forming the same representation are related by (ignoring common factors): where P x=y u=g and D x=y g=u are a shorthand notation for the real, Cartesian basis functions. For a real R À representation, one must multiply by a complex factor i. One advantage of working in the complex representation is immediately apparent from the tables: every operation reduces to a permutation operation multiplied by a factor (termed generalised permutation). Another is that the basis functions for each IrRep are complex conjugates of one another (ungerade IrReps are also multiplied by a factor of À1), from which we can deduce that the on-diagonal Hamiltonian elements are equal and real and that some elements in the upper-triangular elements are complex conjugates of one another. For example: The above relation implies that the Taylor coefficients of such Hamiltonian elements are related. If we know one term to form part of an invariant expression, we can infer the other term must also form part of that invariant expression. Fig. 1 shows the relationship between elements of states of symmetries pertinent to this discussion.
Concretely, for a term such as to form a basis for some IrRep of the group it must be part of the combination since there are some operations in the group that will permute the degenerate basis functions. Operating with any arbitrary operation of the point group:Ô where K R are the product of factors obtained from the matrix elements of operation R for each IrRep making up the first and second terms in the expression above. If each term in the above expression results in K R ¼ þ1 under all operations, we can form a totally symmetric term by summing both terms. Additionally, for operations whose representation gives an anti-diagonal matrix (r v and C 0 2 ), both terms leading to K R ¼ -1 can also imply a totally symmetric term, so long as the expression in Eq. (9) is a difference rather than sum of terms permuted by the symmetry operation. Notice that the on/offdiagonal elements of the representation of the operations given in the D 1h tables are complex conjugates. Since in general the factors K R are complex conjugates of each other, it cannot be that one term gives +1 while the other does not. It follows that it is only necessary to test either K R in order to test for invariance. A similar approach was devised by Eisfeld and Viel [11,16] for Jahn-Teller systems.
Using the projection operator formalism would have been another possible approach (and is used in the following sub-section). It would, however, lead to the generation of many unimportant terms, as is show in Section 2.2 for the O case. It is also unneccessary as starting with one term in the expression, it is straightforward to see how the operator would construct the totally symmetric bases by operating with r v ð/Þ and adding both terms (times the appropriate eigenvalue factor).
For the example in Eq. (9), for n ¼ 0; m ¼ 1, operating with S þ 1 ðwÞ on the first term we would indeed get: as well as for any other operation in the group, making the above term, when combined with the second term, a totally-symmetric IrRep. Note that the symmetry of the electronic states is being taken into account when doing this and so strictly speaking the Hamiltonian element and not the polynomial is totally-symmetric. Finally, to obtain real polynomials and states for the model diabatic potential one is required to transform to a real basis: So the above example would yield: This particular example shows the form linear coupling takes between states of P and D symmetry, termed pseudo Jahn-Teller coupling. First and second order polynomial elements generated in this way are provided in Appendix A and Supplementary information lists them up to sixth order. One can easily generate polynomials of higher order.

O (3D) Jahn-Teller symmetry-adapted basis
The strategy described above breaks down for the case of 3dimensional non-Abelian groups, where the operations of the group become harder to disentangle. In the available representation [22], the operations are also generalised permutation (monomial) matrices. But because of the threefold degeneracy of T 1 , the possible number of terms that make up the invariant expression increase (unlike D 1h with a predictable number of 2 terms). Here we take recourse to the projection operator. Although it should be mentioned that other, more general techniques have been devised to generate T 2 related polynomials [12].
The complex matrix representations for the O group was obtained from Herzig and Altmann. They use the following representation in the symmetrised spherical harmonic basis (for matrices see reference [22] or for their real representation in Appendix B): which relate to the cartesian tensor basis (ignoring common factors) as: so that the transformation matrices to the real representation are given by: where E a=b and T y=x=z are shorthands for the appropriate Cartesian tensors given in Eq. (16). The projection operator is expensive to construct since it involves summing all polynomials generated by operating with all 24 operators on a single generating function. For IrReps T 1 ; E and A 1 (6 functions) testing all possible terms such as jA 1 iQ i b Á Q j c Á Q k d hEj would mean that we could construct up to, say fourth order $14,500 combinations. It is therefore desirable to find a more efficient way of discriminating terms without immediate recourse to the projection operator. We decrease this number drastically if we can work in a representation which reduces some operation to a multiplicative factor and test the trial terms for invariance with respect to that operation.
For the E IrRep, the representation available gives us a similar situation to the D 1h group, where most operations are in a diagonal or anti-diagonal representation. For T 1 , we could achieve a similar result by switching to the representation of the eigenfunctions of some operator of the group whose eigenvalues are preferably distinct. Under the chosen representation, we can immediately evaluate trial terms for invariance to that operation, in a similar manner to that described for the D 1h group. To clarify, consider T 1 under the representation of the eigenvectors of C þ 31 , which belongs to the class of operators with eigenvalues Àe Under this new representation we have basis vectors which are complex conjugates of each other, and so we again find that many elements share relationships like those given in Fig. 1 for IrReps of D 1h but for T 1 ; E and A 1 . The effect of operating with C þ 31 on any a trial term, such as will simply be to multiply out the product of eigenvalues of the different terms forming the element (note that E is already diagonal in the available representation for C þ 31 ). Concretely, the first term in the combination above, when acted upon with C þ 31 indeed gives À Á Ã ¼ 1, and so must the second term. This way, the number of trial terms (also termed generator of the expansion) that need to be acted on by the projection operator can be trimmed into a constellation that satisfy these more stringent conditions. This representation allowed us to cut down by more than a factor of eight the number of times we needed to construct P A 1 , from 14,280 to 1710 trial terms. Having collected those trial generators that are invariant with respect to the C þ 31 operation, we rotate them back to the real basis representation having used Eq. (18) for T 1 .
The real representation of matrix operations O R for T 1 are also generalised permutation matrices, facilitating the construction of P A 1 (provided in Appendix B). We then operated on the trial generators with P A 1 ¼ P R O R to project out the full A 1 expression. Finally, one can rotate to the real representation for E, using the second equation in (17). First order polynomial elements generated in this way are provided in Appendix B for states and coordinates of {A 1 ; E; T 1 } symmetry. In contrast to the T 2 Â t 2 Jahn-Teller Hamiltonian [12,18], degeneracy is not lifted to first order in t 1 , but instead occurs along modes of e symmetry. Supplementary information lists all matrices up to fourth order. For practical reasons, only three dimensional complex polynomials were tested to generate them. This guarantees all three dimensional terms possible.

Test case: vibronic model using symmetry-adapted bases for D 'h ; the absorption spectrum of acetylene
We demonstrate how these polynomials can be used to fit a diabatic potential model to EOM-CCSD/aug-cc-pVTZ energies for the first ten singlet states of acetylene (Molpro [23] was used to obtain the ab initio energies). Details for the construction of vibronic coupling Hamiltonian and their use in dynamics calculations and calculation of spectra can be found elsewhere [3,4]. The states considered are coupled by Renner-Teller and pseudo Jahn-Teller coupling along P g and P u curvilinear coordinates. To correctly describe the first three singlet states (R þ u and D u ) inside this 4D subspace, it was necessary to include coupling terms to the next three higher lying doubly degenerate states (P u ; P g and D g ). Since the resulting model has the direct product expansion form required for the MCTDH algorithm [24], we subsequently used it to calculate the 4D absorption spectrum. Comparison to the experimental spectrum served as a validation of the basis functions presented in Section 2.1. We have also expanded this model to include all internal degrees of freedom to study the vibrationally mediated photodissociation of acetylene [25].

Model system
Mass-scaling the cartesian coordinates relative to the reference geometry: ðx i À x i0 Þ ! 1 ffiffiffiffi m i px i and diagonalising the Hessian gives us seven internal normal modes, four of which belong to degenerate P g and P u IrReps and which describe trans and cis displacements. By rotating the same-plane components of P g and P u normal modes we obtain a set of coordinates approximately describing the cartesian displacements of the hydrogen atoms away from their equilibrium geometry (termed here as 'quasi-cartesian' coordinates): where Q ay is the quasi-cartesian y-coordinate of H-atom 'a' (similarly for Q by ) q cy and q ty refer to the cis and trans normal coordinates along the yz-plane. The xz-plane coordinates, and the symmetric and antisymmetric stretches (z-axis) are treated in a similar manner. We termed them 'quasi' since the positions of other atoms change slightly so as to fix the centre of mass and to maintain the linear space of rotations and translations outside the coordinate space.
Since we are using mass-scaled coordinates, the KE operator is invariant with respect to rotations. From there we performed a spherical polar transformation of each quasi-cartesian axis describing the position of a Hydrogen atom to obtain two sets of quasi-spherical polar coordinates: and correspondingly for H-atom 'b'. The ground state minimum is linear and thus lies at h a % p 2 ; / a % p; R a % Q z 0 and similarly for Hatom 'b'. These coordinates form representations for the C 2v group. The kinetic energy expression for these coordinates result in a pair of radial and orbital angular momentum operators,L, in (quasi) spherical coordinates for the H atoms a/b: where a radial linear term R i has been factored into the radial part of the wavefunction, U R ðR i Þ R i Á UðR i Þ, to simplify the expression [26]. Finally, one can rotate coordinates / a and / b in the inverse manner to Eq. (19): where H x t and H x c correspond to trans and cis angular motion along the zx-plane. With the same treatment for the zy-plane coordinates, fH x t ; H y t g and fH x c ; H y c g form a basis for P u and P g IrReps respectively. Notice they resemble the original normal coordinates, spanning the same D 1h IrReps, but with an angular character at large displacements. Fig. 2 displays a schematically the steps described above. Although fH x t ; H y t g and fH x c ; H y c g enabled us to use the invariant matrices generated by the method described in Section 2.1, we subsequently rotated the resulting fitted model (following section) to fh; /g coordinates, since they leave the KE operator in the simple form given in Eq. (21). We note that as the model is based on normal mode coordinates the Coriolis coupling between rotations and vibrations is not accounted for.

Fitting model parameters to adiabatic energies
Given the high symmetry of the coordinates there are many cuts one can generate with resulting identical potential energy.
For example, H x t þ H y c results in the same energy as H y t þ H x c ; related by a C 2 rotation they are physically indistinguishable. We explored the surface along all possible 2 dimensional diagonal and anti-diagonal vector cuts for the four D 1h coordinates (Fig. 2). We then rotated to the C 2v symmetry curvilinear coordinates fh; /g (obtained from transformation (20) and before applying Eq. (22)) and similarly explored all possible 2 dimensional cuts. We found ten distinct cuts for which the first 10 singlet state EOM-CCSD/aug-cc-pVTZ energies are shown in Fig. 3; they are labelled in the fh; /g, C 2v coordinates since these are used in the final model (following section). A sparse four dimensional grid with 6 energy points per coordinate was also calculate (these ignore the potential energy symmetry considerations just mentioned).
Diabatic model parameters for the aforementioned energy points were obtained by fitting 78 invariant matrices in this subspace of Renner-Teller coordinate, 10-state landscape with strong coupling terms. A genetic algorithm tailored for the optimisation of these problems was used and will be described elsewhere [27]. This includes up to four dimensional, fourth order polynomials to fit the surfaces shown in Fig. 3. Where appropriate, any fourth order terms with negative coefficients where omitted in order to ensure a bounded potential. All states and modes are coupled by inclusion of all first and second order terms, and the values of these coefficients are provided in Table 1. Only a few third and fourth order terms were included, as the increasing number of parameters pushed the limits of the optimisation algorithm. Fewer terms were used to fit S 0 -no coupling was assumed to occur between S 0 and any excited state (fit not shown in Fig. 3). Against the suggestion that the model might be over-parametrised, an 8 state model with terms up to the same order was attempted but yielded no satisfactory fitting.
The S 1 -S 2 crossing mentioned in the literature [28] can be seen in the top row, third sub-figure in Fig. 3 at an out-of-plane geom- Fig. 3. Adiabatic potential surfaces from a 10-state model that uses curvilinear coordinates compared to EOM-CCSD/aug-cc-pVTZ calculations (points). Energy in eV, coordinates in radians. S 0 is not shown in plots. Each slide shows a physically distinct vector generated by exploring diagonal and anti-diagonal vectors of D 1h and C 2v coordinates. From left to right: first row:

Table 1
Coefficients of first and second order polynomials of the invariant matrices fitted to EOM-CCSD/aug-cc-pVTZ ab initio energies in the 10-state, 4D Renner-Teller coordinates space Pg) and similarly for zy-plane (/) coordinates. k are linear inter-state coupling constants, c and l the quadratic intra and inter-state coupling constants respectively. The sub-script denotes the symmetries of the coordinates involved in that polynomial term. The sup-script refers to the matrix basis used as indexed and provided in Appendix A. Note that some matrices like 2.5 and 2.6 span inter and intra-state elements. Where no index has been provided, it refers to quadratic terms that are identical for degenerate coordinates.
The last two bottom-right sub-figures show the S 1 and S 2 potential energy around the dihedral angles that takes the molecule between the cis and trans minima (in this 4D subspace) of the S 1 and S 2 surfaces respectively. The model therefore allows for dihedral rotation to occur properly. Its worth noting that the degeneracy of D states is only broken to fourth order, in contrast to second order for P states.
According to literature [29], there should be no stabilisation away from the linear geometry arising from intra-state, diagonal diabatic elements. Instead, it has been proven that any contributions to the energy that stabilise away from the linear geometry should arise solely via pJT coupling. Fig. 4 shows the model diabatic surfaces of the intra-state elements along the P g and P u curvilinear coordinates. Inter-state couplings are also shown in the smaller sub-figures. The diagonal elements exhibit some stabilisation away from the linear geometry; it is therefore clear that the model does not account for all diabatic interactions to higher lying states and cannot be said to be a separable Hilbert subspace. But given the strong pJT effect present in this system, it is an approximation that, as we show below, yields good qualitative agreements. It is worth mentioning that a wavefunction relaxation on the S 1 state retained a D 1h symmetry; the position expectation value results in a linear molecule. That the matrices indeed obey D 1h operations is supported by this result, with no region in the potential being asymmetrically favoured.

Absorption spectrum
The S 1 minimum has been found experimentally [30] and theoretically [28] to be at a trans bent geometry with an isomerisation barrier to its cis conformer via a torsional motion. Consequently, the absorption spectrum is dominated by a long progression arising from S 1 trans-bending mode as well as contributions arising from CC and possibly CH (totally symmetric) stretch modes. The S 0 ! S 1 maximum intensity is not discernible as it lies higher than the onset of the S 0 ! S 2 transition (6.71 eV). The S 0 ! S 1 maximum is believed to lie 1.5 eV above the band origin [31].
The absorption spectrum [32] in the range provided in experiments must be generated by the first three singlet excited states. Theoretical spectrum from ab initio studies were first calculated by Peric et al. [33] and more recently Köppel and co-workers have built vibronic models to characterise the absorption spectrum. The latter used the Liu et al. [34] Jacobi coordinate formulation with two angles that approximately describe CCH angles, a CC bond and a torsional angle. They calculate the PES for S 1 using CASPT2 and later for S 1 and S 2 using MR-CISD+Q levels of theory (Ref. [31,35] respectively). The last model of the series includes non-adiabatic transfer between S 1 and S 2 along the torsional motion using a regularised diabatic representation [36]. Their spectra are in very good qualitative agreement with experiment. Compared to this earlier work, the coordinates used here should be more flexible and able to describe the photo-dissociation of the C-H bond, while including non-adiabatic effects.
To reproduce the spectrum with this model, a nuclear wavepacket was propagated for 200 fs on these surfaces using the MCTDH algorithm [24]. Mode combination of coordinates fh a ; / a g (and similarly for b) were used, since these are coupled via the kinetic energy operator (Eq. (21)). At least ten two-dimensional single particle functions (SPF) were used to represent the time dependent basis for each state, with fifteen for the first three excited states. The number of primitive grid points for h and / coordinates were 69 and 85 respectively, having used the 2D Legendre DVR basis for pairs h a ; / a and h b ; / b . The relaxed ground state wavefunction is initially operated on by the transition dipole from the hS 0 j l j S 1 i elements of this coordinate subspace using a linear approximation. By multiplying the autocorrelation of the resulting propagated wavefunction by a trigonometric damping function and Fourier transforming this product, we obtain a theoretical estimate for the experimental absorption spectrum (Fig. 5).
The resulting 4D spectrum (top left in Fig. 5) is in qualitative agreement with experiment around the energy range corresponding to S 1 and S 2/3 cis and trans progressions, but is clearly missing progressions at energies 6.75-7.5. It is likely that the missing progressions at higher energies arise from vibrational states arising from the CC and CH (totally symmetric) stretch modes, with higher ground state harmonic frequencies of x 0 % 0:125 and % 0:22 eV respectively. To corroborate this, we calculated the spectrum for the CC stretch mode, shown in the top right plot of Fig. 5, labelled 1D. Superimposed in red, the spectrum for the 4D model is shown in the same plot; one can appreciate the relative intensity contributions to the final spectrum arising from these different coordinates. The lower left panel of Fig. 5 also shows the spectrum obtained from the convolution of the spectra from the four dimensional Renner-Teller subspace and CC stretch mode models. This is effectively a 5D model with no coupling terms between the CC stretch mode and 4D Renner-Teller coordinates. The resulting spectrum is in qualitative agreement with the experimental spectrum, also provided in Fig. 5, as well as with work done by Köppel et al. [31].

Conclusions
We have presented two strategies to generate symmetryadapted bases of matrices of real polynomials representing nuclear coordinates over electronic diabatic elements that are invariant with respect to non-Abelian point group operations. The two cases for which we generated these bases are the D 1h and O, but this approach should work for any group. Working in the simplest possible representation, i.e. a complex symmetrised bases representation which leaves the matrix representation of the IrReps of the groups as generalised permutation matrices, and the use of the projection operator were two indispensable tools that made the generation of these polynomials tractable. The D 1h polynomials generated were subsequently used to construct a four dimensional, ten state diabatic model of acetylene which adequately fits the ab initio EOM-CCSD energies. This includes a coupled manifold of nine excited states and an uncoupled ground-state. The absorption spectrum was calculated to validate the model thus constructed. Owing to the strong pTJ coupling exhibited by S 1 to high lying states, without recourse to the D 1h polynomial basis presented here it would not have been possible to construct such a model.

Conflict of interest
No conflict of interest.

Acknowledgements
We would like to thank Dr. S. Neville and Prof. P. Herzig for valuable discussions and advice. Christopher Robertson received funding for this work from The University of Birmingham.  Matrix representations of the O group in the representation given in Eq. (15), but with T1 in its real representation given by Eq. (17). g ¼ expð2ip=3Þ. Modified from [22].