Symmetry adapted ro-vibrational basis functions for variational nuclear motion calculations: TROVE approach
aa r X i v : . [ phy s i c s . a t m - c l u s ] A ug Symmetry adapted ro-vibrational basis functionsfor variational nuclear motion calculations:TROVE approach
Sergei N. Yurchenko, ∗ , † Andrey Yachmenev, ‡ and Roman I. Ovsyannikov ¶ † Department of Physics and Astronomy, University College London, London, WC1E 6BT,UK ‡ Center for Free-Electron Laser Science (CFEL), DESY, Notkestrasse 85, 22607 Hamburg,Germany ¶ Institute of Applied Physics, Russian Academy of Sciences, Ulyanov Street 46, NizhnyNovgorod, Russia 603950.
E-mail: [email protected]
August 25, 2017
Symmetry adaptation without symmetry operations A A E E i n E i n Y = i n Y H i )( ˆ C ( M ) v Abstract
We present a general, numerically motivated approach to the construction of sym-metry adapted basis functions for solving ro-vibrational Schr¨odinger equations. The pproach is based on the property of the Hamiltonian operator to commute with thecomplete set of symmetry operators and hence to reflect the symmetry of the system.The symmetry adapted ro-vibrational basis set is constructed numerically by solving aset of reduced vibrational eigenvalue problems. In order to assign the irreducible repre-sentations associated with these eigenfunctions, their symmetry properties are probedon a grid of molecular geometries with the corresponding symmetry operations. Thetransformation matrices are re-constructed by solving over-determined systems of linearequations related to the transformation properties of the corresponding wavefunctionson the grid. Our method is implemented in the variational approach TROVE andhas been successfully applied to a number of problems covering the most importantmolecular symmetry groups. Several examples are used to illustrate the procedure,which can be easily applied to different types of coordinates, basis sets, and molecularsystems. Symmetry plays an important role in computing ro-vibrational spectra of polyatomic molecules,particularly in variational solutions of the Schr¨odinger equation. Using a symmetry adaptedbasis set can considerably reduce the size of the Hamiltonian matrix depending on the sym-metry group. For example, in low C s symmetry (with inversion being the only non-trivialsymmetry operation), the use of symmetric and antisymmetric basis functions reduces thematrix by a factor of 2. In higher T d symmetry, the Hamiltonian matrix is split into 10independent blocks, of which only 5 are needed to determine the unique energies and wave-functions of the molecular system (see Fig. 1). For methane, a five-atomic molecule, this isa huge advantage considering the complexity and size of the ro-vibrational computations. If calculating only the energy levels of a molecule, a symmetry adapted basis set is notessential and any sensible basis should lead to a physically meaningful solution. However,knowledge of the symmetry properties of the eigenvectors is vital for generating spectra,2ainly due to the selection rules imposed by the nuclear spin statistics associated withdifferent irreducible representations. Nuclear spin statistical weights give the degeneracy ofthe ro-vibrational states and contribute to the intensity of a transition. Importantly, someenergy levels have zero weights and do not exist in nature. Without knowledge of howthe eigenvectors transform under the symmetry operations, it is impossible to describe themolecular spectrum correctly. From a practical perspective, intensity calculations are alsomuch more efficient in a symmetry adapted representation.The most common symmetry adapted representation is the Wang basis functions, whichare simply symmetric and asymmetric combinations of primitive basis functions. Such com-binations are sufficient for building symmetrized basis sets for Abelian groups, which consistof one-dimensional irreducible representations only, and this is routinely done in most ro-vibrational applications. It is, however, more challenging to symmetrize the basis set fornon-Abelian groups, where the result of the group transformations involve linear combina-tions of basis functions and cannot be described by simple permutations. There exist only ahandful of ro-vibrational methods in the literature capable of dealing with multidimensionalsymmetry group representations. Some examples of the variational approaches include worksby Nikitin et al. , ˇCejchan and Spirko , Boudon et al. , Yurchenko et al. , , Pavlyuchkoet al. , Cassam-Chenai et al. , F´abri et al. .TROVE (Theoretical ROVibrational Energies) is a general method and an associ-ated Fortran 2003 program for computing ro-vibrational spectra and properties of small tomedium-size polyatomic molecules of arbitrary structure. It has been applied to a largenumber of polyatomic species, most of which are characterized by a high degree ofsymmetry ( C , D , D and T d symmetry groups). TROVE has proven very efficient forsimulating hot spectra of polyatomic molecules and is one of the main tools of the ExoMolproject. The most recent updates of TROVE have been reported in Ref. 28,29. Becauseof the importance of symmetry in intensity calculations, TROVE uses an automatic ap-proach for building the symmetry adapted basis set. In this paper we layout the TROVE3ymmetrization approach, which is a variation of the matrix symmetrization method.The matrix symmetrization can be traced back to the original works by Gabriel , Moro-zova and Morozov , , Moccia and was later extensively developed in a series of papers by,for example, Dellepiane et al. , Chung and Goodman , Bouman and Goodman , Jordanovand Orville-Thomas , Chen et al. . The main idea of these studies is to use a diagonalizionof matrices representing specially constructed symmetry operators. Using this technique asymmetry adaptation can be obtain without the use of symmetry operations . For example,Moccia used the nuclear attraction matrix to build symmetry adapted molecular orbitals,or a Wilson G + F − matrix in symmetrized force constants calculations; Dellepiane et al. used a kinetic G -matrix to obtain symmetry adapted representations of vibrational molecu-lar modes; Chung and Goodman used an overlap matrix of atomic orbitals to symmetrizethem. Chen et al. proposed an ‘eigenfunction’ method based on eigenfunctions of a linearcombination of symmetry operations from the so-called complete set of commuting oper-ators (CSCO), which was then extensively employed for constructing symmetry adaptedrepresentations of coordinates and basis functions for (ro-)vibrational calculations .Here we apply the idea of the matrix symmetrization to numerical construction of sym-metry adapted ro-vibrational representations of a ro-vibrational Hamiltonian ˆ H for a generalpolyatomic molecule. In our version of this method, the symmetry adapted basis functionsare generated as eigenvectors of some reduced rovibrational Hamiltonians. These operatorsˆ H (red) are derived from ˆ H such that (i) they represent different vibrational or rotationalmodes and (ii) they are symmetrically invariant to ˆ H . According to the matrix symmetriza-tion method, the eigenvectors of ˆ H (red) necessarily transform according to irreducible repre-sentations (irreps) of the symmetry group.Not only does this allow us to construct the symmetry adapted basis functions, but alsoto improve and contract the basis set via standard diagonalization/truncation procedures.The relative simplicity of this procedure means it can be straightforwardly implemented inmany existing nuclear motion programs. It may also be interesting to apply the method in4uantum chemical approaches, where the initial set of symmetry adapted atomic orbitalscan, for example, be constructed by diagonalizing the bare nuclear Hamiltonian.The explanation of our method will be given in the form of practical illustrative examples,rather than using rigorous group-theoretical formalism. The paper is structured as follows:The main idea of the TROVE symmetrization approach is described in Section 2. Sections 3.1and 3.2 present illustrative examples for XY and XY type molecules. Readers interestedin implementation of the method should read Section 4, where the sampling technique forreconstructing the symmetry transformation properties of vibrational wavefunctions is intro-duced, and Section 5, which details the TROVE reduction method based on the projectionoperator approach. A non-rigid, ammonia-type molecule XY of D (M) molecular symme-try is used as an example to illustrate this part of the method implementation. As a veryspecial case, the degenerate mutidimensional isotropic Harmonic oscillator basis functionsare considered in Section 5.3 with an example shown in the Appendix. Symmetrization ofthe rotational and total ro-vibrational basis functions is realized using standard reductiontechniques and this is discussed in Sections 5.4 and 5.5. In order to introduce the TROVE symmetrization approach, we consider a general multidi-mensional ro-vibrational Schr¨odinger equationˆ H Ψ rv = E Ψ rv , (1)which is to be solved variationally using the ro-vibrational basis set in a product-form:Φ Jk,ν ( θ, φ, χ, q , q , . . . , q N ) = | J, k, m i φ n ( q ) φ n ( q ) . . . φ n N ( q N ) , (2)5 A E E F F F F F F Figure 1: The block-diagonal structure of a Hamiltonian matrix in the T d irreducible repre-sentation. The empty (white) cells indicate vanishing matrix elements. Only 5 blocks of theirreducible representations A , A , E a , F a , and F a are needed as the other five matrices( E b , F b , F c , F b , F c ) contain degenerate solutions.where φ n i ( q i ) is a one-dimensional (1D) vibrational function, n i is a vibrational quantumnumber, q i is a generalized vibrational coordinate, N is the number of vibrational degrees offreedom, | J, k, m i are the rigid-rotor wavefunctions, k = − J . . . J and m = − J . . . J are therotational quantum numbers (projections of the total angular momentum onto the molecule-fixed z and laboratory-fixed Z axes, respectively), ν = { n , n , . . . n N } is a generalizedvibrational multi-index. The primitive basis functions φ n i ( q i ) ≡ | n i i are any vibrational 1Dfunctions from a orthonormal set (e.g. Harmonic oscillator wavefunctions). In the absenceof external fields m does not play any role and can be omitted.Let us assume that the molecule belongs to a molecular symmetry group G consisting of g elements (group operations) R . We aim to construct symmetry adapted basis set functionsΨ J, Γ s µ which transform according to irreducible representations Γ s ( s = 1 . . . r ) of G . Here µ is a counting number and Γ s will be referred to as a ‘symmetry’ or an ‘irrep’ of G . For an l s -fold degenerate irrep, and when we will need to refer to specific degenerate components6f Ψ J, Γ s µ , an additional subscript n = 1 , . . . , l s will be used as, e.g. Ψ J, Γ s µ,n . For example, forthe two-fold degenerate E symmetry, n = 1 and 2 corresponds to the E a and E b symmetrycomponents, while in case of the three-fold degenerate F symmetry, these are F a , F b and F c .Additionally, we will require that the transformation properties of multi-fold irreps (e.g. E or F representations) are known.We now assume that the symmetry adapted basis functions Ψ J, Γ s µ,n can be represented bylinear combinations of the sum-of-product primitive functions from Eq. (2) byΨ J, Γ s µ,n = X k,ν T µ,J, Γ s k,ν,n Φ Jk,ν , (3)where T µ,J, Γ s k,ν,n are symmetrization coefficients. The important advantage of the symmetryadapted basis set is that the corresponding Hamiltonian matrix has a block-diagonal form(see Fig. 1): h Ψ J, Γ s µ,n | H rv | Ψ J, Γ t µ ′ ,n ′ i = H µ,µ ′ δ s,t δ n,n ′ . (4)In practice, this means that each ( J, Γ s , n )-block can be diagonalized independently with J and (Γ s , n ) as good quantum ‘numbers’ (i.e. constants of motion). The main goal of thiswork is to present a general numerical algorithm for constructing symmetrization coefficients T µ,J, Γ s k,ν,n for a molecule of general structure and symmetry.According with the matrix symmetrization method (see, for example, Jordanov andOrville-Thomas ), symmetry adapted set of wavefunctions can be constructed by diago-nalizing matrices representing some operators ˆ A . These operators are chosen to be invariantto the symmetry operations R ∈ G . Our approach is based on the realization that in prin-ciple ˆ H itself would be an ideal choice for ˆ A , as it has the right property to commute withany R from G [ ˆ H, R ] = 0 . (5) Here we assume that there exists isomorphism between the elements R of G and the correspondingrepresentations, and use the same symbol R in both cases. H are also eigenfunctions of R (up to a linear combination ofdegenerate states) and hence transform as one of the irreps of the system (see, for example,the textbook by Hamermesh ). Obviously, it makes no sense to use the ro-vibrationalHamiltonian operator ˆ H for this purpose. Instead, we define a set of reduced Hamiltonianoperators ˆ H ( i ) derived from ˆ H as follows. (i) All ro-vibrational degrees of freedom are dividedinto L symmetrically independent subspaces, which form subgroups of G . (ii) For each i thsubspace ( i = 1 . . . L ) a reduced Hamiltonian operator ˆ H ( i ) is constructed by neglecting orintegrating over all other degrees of freedom. (iii) The symmetry adapted wavefunctions foreach i th subspace are obtained by diagonalizing the corresponding ˆ H ( i ) . (iv) The total basisset is built as a direct product of the subspace bases and then transformed to irreps usingstandard reduction approaches.Symmetrically independent subspaces of coordinates are selected such that each subspacecontains only the coordinates related by symmetry operations of the group. For example, thevibrational motion of a molecule XY spanning the molecular symmetry group C (M) can bedescribed by two stretching and one bending mode, which transform independently and canthus be separated into two subspaces. More specifically, the bond lengths r (X–Y ) and r (X–Y ) are two stretching vibrational modes connected through symmetry transformationsof the group C (M), which form the subspace 1, while the interbond angle α (Y –X–Y )belongs to the subspace 2, with the transformation properties shown in Table 1.To explore Eq. (5) for constructing a symmetry adapted basis, we define and solve aset of eigenvalue problems for reduced Hamiltonian operators ˆ H ( i ) . For each subspace i ( i = 1 . . . L ) a reduced eigenvalue problem is given byˆ H ( i ) ( Q ( i ) )Ψ ( i ) λ i ( Q ( i ) ) = E λ i Ψ ( i ) λ i ( Q ( i ) ) , (6)where Q ( i ) is a set of coordinates { q k , q l , . . . } from a given subspace i , E λ i is an eigenvalueassociated with the eigenfunction Ψ ( i ) λ i and λ i counts all the solutions from the subspace i .8he resulting solutions Ψ ( i ) λ i should transform according with an irrep Γ s of G and one ofits degenerate components n (holds for l s > ( i ) , Γ s λ i will be used, or even Ψ ( i ) , Γ s λ i ,n to further specifyits degenerate components.The reduced Hamiltonian ˆ H ( i ) is constructed by averaging the total vibrational ( J = 0)Hamiltonian ˆ H on the ‘ground state’ primitive vibrational basis functions φ n s ( q s ) = | n s i from other subspaces ( { s } 6∈ { i } ) as given byˆ H ( i ) ( Q ( i ) ) = h p |h q | . . . h r | ˆ H | r i . . . | q i| p i , (7)where | s i is a primitive basis function φ n s ( q s ) with n s = 0 and { p, q, r } are coordinates fromother subspaces, i.e. { p, q, r } 6∈ { i } .For example, in the case of an XY molecule the two reduced Hamiltonian operators canbe formed as ˆ H (1) ( r , r ) = h | ˆ H | i , (8)ˆ H (2) ( α ) = h |h | ˆ H | i| i , (9)where Q (1) = { r , r } and Q (2) = { α } define the partitioning of the three coordinates intotwo subspaces i = 1 and 2.Equation (6) represents the main idea of the method, which will be referred to as TROVEsymmetrization: since ˆ H ( i ) commutes with any R ∈ G , the eigenfunctions Ψ ( i ) λ i ( Q ( i ) ) mustnecessarily span one of the irreducible representations Γ s of the group G . By solving Eq. (6),not only do we get a more compact basis set representation which can be efficiently contractedfollowing the diagonalization/truncation approach, it is also automatically symmetrized. Thetotal vibrational basis set is then constructed as a direct product of L symmetrically adaptedbasis sets followed by a reduction to irreducible representations using standard projectionoperator techniques (see, for example, Ref. 43). The major advantage of this symmetrization9pproach is that it can be formulated as a purely numerical procedure, which is particularlyvaluable for handling the algebra of symmetry transformations to describe high vibrationalexcitations. The required components (Hamiltonian matrices and eigensolvers) are readilyavailable in any variational program and thus the implementation of the present approachinto a variational ro-vibrational calculation should be relatively straightforward.There are however two major problems to overcome: (i) Even though we know thatΨ ( i ) λ i from Eq. (6) should transform as an irrep Γ s , we do not automatically know whichone, except for trivial one-dimensional subspaces; (ii) the degenerate solutions (e.g. forΓ s = E, F, G ) are usually represented by arbitrary mixtures of the degenerate componentsand do not necessarily transform according to standard irreducible transformation matrices(see also examples below). The latter is a common problem of degenerate solutions, sinceany linear combination of degenerate eigenfunctions is also an eigenfunction. For example,when TROVE solves Eq. (6) by a direct diagonalization using one of the numerical linearalgebra libraries (e.g. DSYEV from LAPACK), the degenerate eigenfunctions come out asunspecified mixtures of degenerate components. We will demonstrate in Section 5 that ageneral reduction scheme can be used to recast the degenerate mixtures such that they followthe standard transformation properties upon the group operations. It should be noted thatthe eigenvalue symmetrization method by Lemus , can in principle be used to resolve thedegenerate components by constructing a proper CSCO.Table 1: Transformation properties of the internal coordinates r , r , and α of an XY -typemolecule and the characters of the irreps of the C (M) group.Coordinate E (12) E ∗ (12) ∗ r r r r r r r r r r α α α α α Irrep Γ Characters χA A B B Examples -type molecules In order to demonstrate how TROVE symmetrization based on Eq. (6) works, we againconsider an XY triatomic molecule. It spans the Abelian group C (M) with well-knownsymmetry adapted combinations of vibrational basis functions given by (compare to Eq. (3)):Φ A n ,n ,n = 1 √ φ n ( r ) φ n ( r ) + φ n ( r ) φ n ( r )] φ n ( α ) , n = n (10)Φ B n ,n ,n = 1 √ φ n ( r ) φ n ( r ) − φ n ( r ) φ n ( r )] φ n ( α ) , n = n (11)Φ A n,n,n = φ n ( r ) φ n ( r ) φ n ( α ) , n = n ≡ n, (12)where A and B are two irreducible representations of C (M) (see Table 1). The ‘irre-ducible’ functions Φ A n ,n ,n and Φ B n ,n ,n are also eigenfunctions of the group operators R = { E, (12) , E ∗ , (12) ∗ } , e.g. (12) Φ A ν = Φ A ν , (13)(12) Φ B ν = − Φ B ν , (14)where ν stands for { n , n , n } . The transformation of the ‘reducible’ primitive functions | n i| n i| n i = φ n ( r ) φ n ( r ) φ n ( α ) (for n = n ), that are not eigenfunctions of R = (12),involves two different states: (12) | n i| n i| n i = | n i| n i| n i , (15)(12) | n i| n i| n i = | n i| n i| n i . (16)Now we derive irreducible combinations of | n i| n i| n i using the numerical approach ofEq. (6). As an example here we use the vibrational wavefunctions of the H S molecule11btained variationally with TROVE based on the potential energy surface from Ref. 45.It should be noted however that any computational approach using the same coordinateswould essentially give equivalent expansions. We construct the matrix representations ofthe reduced Hamiltonians in Eqs. (8) and (9) in the basis of 1D functions φ n ( r ), φ n ( r ),and φ n ( α ) determined using the Numerov-Cooley approach as described in Ref. 7. Forsimplicity we use a small basis set limited by the polyad number P max = 2 as given by P = n + n + n ≤ P max . After solving the reduced eigenvalue problem for ˆ H (1) (Eq. (8)), the following variationalwavefunctions were obtained:Ψ (1)1 ( r , r ) = 0 . | , i + 0 . | , i + | , i ) + . . . , (17)Ψ (1)2 ( r , r ) = 0 . | , i − . | , i + | , i ) + . . . , (18)Ψ (1)3 ( r , r ) = − . | , i − | , i ) + . . . , (19)where we have used a shorthand notation | n , n i = | n i| n i . When compared to Eqs. (10,11),the eigenfunctions Ψ (1) λ have the expected symmetrized form and are classified according tothe A and B irreps, i.e. as Ψ (1) ,A , Ψ (1) ,A , and Ψ (1) ,B . Thus the expansion coefficients T µ,J, Γ s k,ν,n = T µ, Γ s { n ,n } in Eq. (3) are obtained numerically without any assumption on the sym-metries. Here J, k = 0 (rotational indices) and n = 1 (degenerate component) are omittedfor simplicity and ν = { n , n } . The numerical error of the symmetrization can be assessedby comparing Eqs. (17–19) to Eqs. (10,11). For example, the differences between T ,A { , } and T ,A { , } , T ,A { , } and T ,A { , } , T ,B { , } and − T ,B { , } are found to be within 10 − .Increasing the size of the basis set (using larger P max ) will lead to analogous expan-sions involving symmetrized contributions from higher excitations | n , n i . The new reducedwavefunctions Ψ (1) , Γ s λ ( r , r ) together with Ψ (2) , Γ ′ s λ ( α ) (eigenfunctions of H (2) in Eq. (9)) canbe utilized to build the new contracted and symmetrized basis set, which is then used to12iagonalize the complete Hamiltonian ˆ H . In this simple example the symmetry proper-ties of the expansion coefficients, as well as of the corresponding wavefunctions, are trivial.However, our goal is to develop a general numerical symmetrization algorithm applicableto arbitrary basis sets, coordinates, symmetries or molecules, which is also in line with theTROVE ideology of a general, black-box like program. As will be demonstrated below, theadvantage of our automatic symmetry classification method becomes more pronounced forlarger molecules with more complicated symmetry, especially for ones containing degeneraterepresentations. -type, C -symmetry Here we present another example of a rigid pyramidal tetratomic molecule XY , charac-terized by the C (M) molecular symmetry group. We choose six internal coordinates as∆ r , ∆ r , ∆ r (bond length displacements) and ∆ α , ∆ α , ∆ α (the interbond angle dis-placements). The associated permutation symmetry operations and characters of C (M)are collected in Table 2. These coordinates, as well as the corresponding 1D primitive basisfunctions | n i i ( i = 1 . . . { ∆ r , ∆ r , ∆ r } and subspace 2 is { ∆ α , ∆ α , ∆ α } . We assume that | n i , | n i , and | n i are the 1D stretching basis functions of ∆ r , ∆ r , and ∆ r , respectively, and | n i , | n i , and | n i are the 1D bending functions of ∆ α , ∆ α , and ∆ α , respectively. The two reducedSchr¨odinger equations for subspaces 1 and 2 are given by:ˆ H (1)str Ψ (1) λ = E λ Ψ (1) λ , ˆ H (2)bnd Ψ (2) λ = E λ Ψ (2) λ , (20)where the reduced 3D Hamiltonian operators ˆ H (1)str and ˆ H (2)bnd are obtained by vibrationallyaveraging the total vibrational Hamiltonian ˆ H over the ground state basis functions from13ubspace 2 and 1, respectively:ˆ H (1)str (∆ r , ∆ r , ∆ r ) = h |h |h | ˆ H | i| i| i , ˆ H (2)bnd (∆ α , ∆ α , ∆ α ) = h |h |h | ˆ H | i| i| i . (21)The C (M) group spans the Γ s = A , A , and E representations, where E is two-fold. Fol-lowing the discussion above, we expect the resulting wavefunctions Ψ ( i ) λ i to be eigenfunctionsof all six symmetry operators R of C (M) from Table 2, i.e. to transform according to A , A , or E .Table 2: The character table and transformation of the internal coordinates of an XY -typemolecule by the symmetry operations of the molecular symmetry group C (M). Variables E (123) (321) (23) ∗ (13) ∗ (12) ∗ r r r r r r r r r r r r r r r r r r r r r α α α α α α α α α α α α α α α α α α α α α CharactersIrrep Γ E (123) (23) ∗ A A E The illustration below is based on the TROVE program again, however it should betransferable, at least in principle, to any method (i.e. any basis set, kinetic energy operatoror potential energy function), as long as a similar choice of vibrational coordinates and aproduct basis of 1D wavefunctions are used. We choose the PH molecule and construct thesymmetrized basis set in TROVE using a polyad number cutoff given by P = 2( n + n + n ) + n + n + n ≤ P max = 10in conjunction with PES of Sousa-Silva et al. . The 1D basis set functions are generated14sing the Numerov-Cooley procedure as described in Ref. 15, where the details on the kineticenergy expansion and the vibrational coordinates can also be found. In practice, we usuallychoose the maximal polyad number P max in the order of 14–20 (see, for example, Refs.15,48–50).The basis set for subspace 1 (stretching) in this case contains only functions with n + n + n ≤ n = n = n = 0, while subspace 2 (bending) basis functions are constructedfrom the contributions n + n + n ≤
10 and n = n = n = 0. The first four variationaleigenfunctions of ˆ H (1)str read (where the shorthand notation | n , n , n i ≡ | n i| n i| n i is used)Ψ (1)1 = 0 . | , , i − .
128 ( | , , i + | , , i + | , , i ) + . . . Ψ (1)2 = 0 . | , , i − . | , , i + | , , i + | , , i ) + . . . Ψ (1)3 , = 0 . | , , i − . | , , i + 0 . | , , i + . . . Ψ (1)3 , = 0 . | , , i + 0 . | , , i − . | , , i + . . . (22)corresponding to the energy term values 0 .
0, 2317.86, 2328.28 and 2328.28 cm − , respec-tively, relative to the zero-point-energy (ZPE) of 5222.59 cm − (the actual coefficients arecalculated in double precision). In fact we find that high numerical accuracy is importantfor numerically reconstructing the symmetries of the eigensolutions. The first two wavefunc-tions Ψ (1)1 and Ψ (1)2 in Eq. (22) exhibit fully symmetric forms (to ∼ − ) and thus belongto A . The symmetry of the solutions Ψ (1)3 , and Ψ (1)3 , cannot be immediately recognized fromtheir expansion coefficients, but the coinciding energy levels (within the numerical error of ∼ − ) indicate a degenerate solution, which for the case of C (M) can only mean the E symmetry. For these two states we use the subscript notation λ, n to refer to two degeneratecomponents ( n = 1 ,
2) of the degenerate state λ = 3.Thus, all four wavefunctions (as well as other solutions not shown here) appear readilysymmetrized. However, guessing the degenerate E symmetries based on the degeneracy ofenergies, which worked here, is not always reliable in actual numerical calculations, espe-15ially for high excitations and arbitrary potential functions with accidentally close energies(accidental resonances). In fact, some reduced Hamiltonian operators can lead to degener-ate solutions of arbitrarily order, such as, for example, Harmonic oscillator or Rigid rotorHamiltonians. The A states as single energy solutions can be also easily mixed up with A .Finally, the degenerate components are usually come out of diagonalizations as arbitrarydegenerate mixtures. For example, the eigenfunctions Ψ (1)3 , and Ψ (1)3 , from Eq. (22) wereobtained using the eigen-solver DSYEV (LAPACK). As a result, they do not necessarilyobey the standard E -symmetry irreducible transformation rules: for example the D [(123)]transformation (which represents the C rotation about the axis of symmetry) does nottransform the two E -symmetry components (Ψ (1)3 , , Ψ (1)3 , ) according to the transformation D [(123)] = cos ǫ sin ǫ − sin ǫ cos ǫ with ǫ = 2 π/
3, as one would expect. In principle, such functions with randomly mixedcomponents would still lead to a block-diagonal representation of the Hamiltonian matrixas in Eq. (4), and thus does not seem to be a problem. However having the functions Ψ ( i ) λ i ,n which transform according with standard transformation rules is very useful for reducing thedirect products of the corresponding irreducible representations. Therefore these randomlymixed degenerate components have to be further recast by a proper linear transformationto the standard form, which will be referred to as standard representations.To conclude this section, the matrix symmetrization method based on reduced Hamil-tonian operators can be efficiently used to produce a symmetry adapted basis set in fullynumerical fashion. However, the method does not tell which irreps these functions belongto and, consequently, which symmetry properties they have; besides, the degenerate compo-nents are mixed by an arbitrary orthogonal transformation which makes it difficult to usein subsequent calculations. This is where the second step of our symmetrization procedure,namely the symmetry sampling, comes in. 16 Symmetry sampling of the eigenfunctions
In this section we show how to reconstruct the symmetries Γ s of the eigenfunctions Ψ ( i ) λ i ,n from Eq. (6) by analyzing their transformation properties and also how to bring their de-generate components n into the ‘standard’ form. Towards this end, we select a grid of N ( i )grid instantaneous sampling geometries and use them to probe the values of the eigenfunctionsΨ ( i ) λ i ,n as well as of their symmetry related images. That is, for a given subspace i and se-lected instantaneous geometries Q ( i ) k ( k = 1 . . . N ( i )grid ), all symmetry related images R Q ( i ) k are generated. These are used to reconstruct the values of the wavefunctions Ψ ( i ) λ i ,n ( R Q ( i ) k )at the new geometries, and to establish the transformation matrices D [ R ] for each operation R of the group G . This is different from the more common practice of directly exploring thepermutational properties of the wavefunctions. At this point it might appear that permutingthe wavefunctions would be easier, at least for the case of | n i| n i| n i in our example of therigid XY molecule. However, as will be shown below, applying the group transformationsto the coordinates instead of the basis functions provides a more general numerical approachapplicable to complex cases when the permutation symmetry properties of the wavefunctionsare not obvious.Let us consider an l λ -fold degenerate eigenstate λ with l λ eigenfunctions Ψ ( i ) λ,n ( n =1 , . . . , l λ ) from a subspace i , and define a grid of randomly selected geometries Q ( i ) k ( k =1 . . . N ( i )grid ). We assume that the transformation properties of the coordinates from a givensubspace with respect to R are known at any specific point k . This can be expressed as: R Q ( i ) k = Q ′ ( i ) k (23)with each subspace being independent from the others by definition. Under the assumptionthat the eigenfunctions Ψ ( i ) λ,n ( Q ( i ) ) can be evaluated at any grid point k , i.e. at any instan-taneous geometry Q ( i ) k including their symmetry related images Q ′ ( i ) k (which is true for the17ROVE program), we can evaluateΨ ( i ) λ,n ( k ) ≡ Ψ ( i ) λ,n ( Q ( i ) k ) , (24)Ψ ′ ( i ) λ,n ( k ) ≡ R Ψ ( i ) λ,n ( Q ( i ) k ) = Ψ ( i ) λ,n ( Q ′ ( i ) k ) , where Q ′ and Ψ ′ are the transformed coordinates and functions, respectively. The eigen-functions Ψ ( i ) λ,n ( Q ( i ) k ) and their symmetric images R Ψ ( i ) λ,n ( Q ( i ) k ) are also related via the trans-formation matrices as given by: R Ψ ( i ) λ,m ( Q ( i ) k ) = l λ X n =1 D [ R ] mn Ψ ( i ) λ,n ( Q ( i ) k ) . (25)It should be noted that we are using the convention by Bunker and Jensen to define theoperations R on the nuclear coordinates and functions. This convension is also referred toas passive (see, e.g., a detailed discussion by Alvarez-Bajo et al. ). For instance, for the E -symmetry wavefunctions from Eq. (22), this expression reads Ψ ′ (1)3 , ( k )Ψ ′ (1)3 , ( k ) = R Ψ (1)3 , ( k )Ψ (1)3 , ( k ) = D [ R ] D [ R ] D [ R ] D [ R ] · Ψ (1)3 , ( k )Ψ (1)3 , ( k ) . (26)It should be noted that the linear system in Eq. (25) does not impose the condition ofunitariness of the solution. As a result the matrices D [ R ] mn can be non-orthogonal andmust be orthogonalizied, for which the Gramm-Schmidt approach is employed.Now by combining Eq. (24) and Eq. (25) we obtain l λ X n =1 D [ R ] mn Ψ ( i ) λ,n ( k ) = Ψ ′ ( i ) λ,m ( k ) . (27)The ( l λ ) elements D [ R ] mn can be determined by solving Eq. (27) as a system of N ( i )grid k = 1 . . . N ( i )grid ) of the type X n A kn x ( m ) n = b ( m ) k . (28)Here the matrix elements A kn = Ψ ( i ) λ,n ( k ) and the vector-coefficients b ( m ) k = Ψ ′ ( i ) λ,m ( k ) areknown, while x ( m ) n = D [ R ] nm are the unknown quantities. Once the system of N ( i )grid linearEquations (28) is solved for each R and all the g transformation matrices D [ R ] are found ( g is the group order), we apply the standard projection operator approach to generate theirreducible representations (see Section 5).The number of degenerate reducible states l λ is simply taken as the number of states withthe same energies. For non-degenerate wavefunctions ( l λ = 1), the sampling procedure willalways produce D [ R ] , = χ [ R ] = ±
1. Accidental degeneracies (e.g. A / A with identicalenergies) are processed as if they were normal degenerate components. In this case theresulting matrices are diagonal, D [ R ] ii = ± D [ R ] i,j = 0 ( i = j ).At least N ( i )grid = ( l λ ) grid points are required to define the linear system (or even lessdue to the unitary property of the transformation matrices). In practice, it is difficultto find a proper set of geometries with all values of Ψ ′ ( i ) λ ( k ) and Ψ ( i ) λ ( k ) large enough tomake the solution of the linear system numerically stable (i.e. with non-vanishing deter-minant). We therefore tend to select more points ( N ( i )grid ≫ ( l λ ) ) and thus solve an over-determined linear system using the singular value decomposition method implemented in theLAPACK/DGELSS numerical procedure. We usually take N ( i )grid = 10–200 geometries Q k randomly distributed within the defined coordinate ranges around the equilibrium geometryof the molecule.This symmetrization procedure can be applied to any primitive functions provided theirvalues can be calculated at any instantaneous geometry. For example, the commonly usedbasis functions in TROVE are 1D eigensolutions of a reduced 1D Hamiltonian determinedusing the Numerov-Cooley procedure and defined on an equidistant grid of geometries, typ-19cally of about 1000 points. In this case the values of the primitive functions φ n k ( q k ) inEq. (2) are obtained by interpolation using the polint procedure. Other popular basissets in TROVE are Harmonic oscillator and Rigid rotor wavefunctions, for which the detailsof the symmetrization procedure are presented below.
Due to the accidental degeneracies and even more so due to the intrinsic degeneracies imposedby some reduced Hamiltonians (e.g. Hamiltonian of isotropic Harmonic oscillators), it iscommon to deal with degenerate solutions of Eq. (6) of high order, which can be much higherthan that of the corresponding irreducible representations. For example, the Hamiltonian ofthe 2D Harmonic isotropic oscillator ˆ H = 12 (cid:2) P a + P b + λ ( Q a + Q b ) (cid:3) has the eigenvalues E v a ,v b = ¯ h √ λ [ v a + v b + 1] , which are ( v a + v b +1)-fold degenerate. As it was discussed above, our numerical symmetriza-tion approach often leads to arbitrarily mixed degenerate representations, which need to befurther transformed to the standard orthogonal form. In the following we show how to usethe standard projection technique to symmetrize such general cases in a fully numericalfashion.In order to reduce a representation Γ red to its irreducible componentsΓ red = a Γ ⊕ a Γ ⊕ a Γ ⊕ · · · ⊕ a h Γ h , (29)the first step is to use the characters χ [ R ] of the reducible representation as traces of the20ransformation matrices D [ R ] mn : χ [ R ] = X n D [ R ] nn and find the number of irreducible representations a s (reduction coefficients) for each irrepsΓ ∈ G as given by: a s = 1 g X R χ [ R ] ∗ χ Γ s [ R ] . (30)Remember that g is the order of the group, R runs over all the elements of the group and χ Γ [ R ] are the group characters. To ensure the numerical stability of the symmetrizationwe usually check if these reduction coefficients are (i) integral and (ii) satisfy the reductionrelations χ [ R ] = X s a s χ Γ s [ R ] , and X R (cid:12)(cid:12) χ Γ s [ R ] (cid:12)(cid:12) = g. (31)If these conditions are not fulfilled (within some numerical thresholds, typically 10 − ), thegrid points are re-selected and the transformation matrices are re-built.In principle a projection onto a non-degenerate irrep Γ s can be generated by the opera-tor: P Γ s = 1 g X R χ Γ s [ R ] ∗ R. (32)However, a non-degenerate function Ψ ( i ) λ obtained using Eq. (6) should already transformirreducibly as Γ s , and therefore this projector P Γ s (Γ s ∈ G ) is not needed. For example,applying the reduction from Eq. (30) to the first two non-degenerate wavefunctions (Ψ (1)1 and Ψ (1)2 ) in Eq. (22) will give a A = 1, a A = a E = 0, which unambiguously defines theirsymmetries.Degenerate solutions require special care. For the sake of generality let us assume thatdegeneracy of the reducible solution l λ can be higher than that of the irreducible representa-tions l s . The degenerate wavefunctions (both accidentally and intrinsically) can be selectedsimply based on the coincidence of energies within a specified threshold (usually 0.001 cm − ).21he corresponding transformation l λ × l λ matrices D [ R ] are constructed using the samplingprocedure of Eq. (25) and then symmetrized with Eq. (30) giving the reduction coefficientsof irreps Γ s .In cases of multiple degenerate states ( l λ > P Γ s mn = l s g X R D Γ s [ R ] ∗ mn R, (33)where D Γ s [ R ] is an irreducible orthogonal transformation matrix of Γ s for an operation R ,and l s is the dimension (degeneracy) of Γ s . Following the symmetrization protocol andusing the reducible D Γ [ R ] mn from Eq. (25) as a representation of R , the m th component ofthe irreducible wavefunction Ψ Γ s λ,m is obtained by acting P Γ s mm (diagonal element) on Ψ Γ red λ,n .Here we distinguish the reducible and irreducible representations by using the superscriptsΓ red and and Γ s , respectively. The off-diagonal elements of the transfer operator P Γ s m,n arethen used to recover other components of Ψ Γ s λ,n .It should be noted that degenerate solutions Ψ Γ red λ,n in general can span more than onerepresentation. Besides the projected vectors are not automatically orthogonal. Thereforethe symmetry classification procedure must include the following steps: (i) The projector P Γ s mm ( m = 1 . . . l λ ) is applied to Ψ Γ red λ,p to form a trial irreducible solution Ψ Γ s λ,m , which isthen (ii) checked against already found functions Ψ Γ s λ,n ( n < m ). The trial function is theneither rejected (if it is already present in the set) or (iii) orthogonalized to this set usingthe Gramm-Scmidt orthogonalization technique. This procedure is repeated until all a s irreducible solutions are found. -type, C -symmetry (Continued) Let us now return to the example above. Choosing 40 points and applying our samplingapproach to the degenerate state Ψ (1)3 ,n in Eq. (22) for all six C (M) group operations listed22n Table 2, the following transformation matrices were determined: D [ E ] = . . . . , D [(123)] = − . . − . − . , (34) D [(321)] = − . − . . − . , D [(23) ∗ ] = − . − . − . . , (35) D [(13) ∗ ] = . − . − . − . , D [(12) ∗ ] = − . . . . . (36)In principle only three matrices are unique, but TROVE currently computes matrices forall representations and does not take the advantage of generators. The characters χ Γ [ R ] ofthese transformations are 2 . − .
0, and 0 . ± − ), which in conjunction with Eq. (30)leads to the following reduction coefficients a E = 1 and a A = a A = 0 ( ± − ) as expectedfor a doubly-degenerate solution.Using the transformation matrices D [ R ] together with Eq. (33), we build a projectionoperator P E and apply it to the degenerate components Ψ a = Ψ (1)3 , and Ψ b = Ψ (1)3 , to obtaina trial vector: ˜Ψ a = 0 . a − . b , which after normalization becomes˜Ψ a = 0 . a − . b . The second component ˜Ψ b is found by applying the transfer operator in Eq. (33):˜Ψ b = 26 X R D Γ s [ R ] ∗ ˜Ψ a which when normalized reads ˜Ψ b = 0 . a + 0 . b .23inally, by applying the transformation vectors to the original (reducible) representation { Ψ a , Ψ b } from Eq. (22) we obtain˜Ψ (1)3 , = ˜Ψ (1) E a = − √ | i| i| i + | i| i| i − | i| i| i ) + . . . (37)˜Ψ (1)3 , = ˜Ψ (1) E b = 1 √ | i| i| i − | i| i| i ) + . . . (38)which is the well-known form that transforms according to the standard E -symmetry rep-resentations of C (see Ref. 4 for example). The expansion coefficients in Eqs. (37,38) aredefined within a numerical error of 10 − .We can check if the new vectors transform correctly, as in this case, i.e. according to thestandard irreducible matrices D Γ s [ R ] as follows R Ψ ( i ) , Γ s λ,m = X n D Γ s [ R ] m,n Ψ ( i ) , Γ s λ,n . As mentioned above, if the projection operator P Γ s mm does not lead to a correct or independentcombination, we would try a different component of P Γ s mm until the correct solution is found(which is guaranteed).With this procedure, symmetries of all eigenstates can be easily reconstructed. For thebasis set P ≤ P max = 10 in this example (see Section 3.2), we computed 38 stretching Ψ (1) λ and 192 bending Ψ (2) λ eigenfunctions, with the symmetries and energies of the first 3 fromeach subspace shown in Table 3.Once all symmetry adapted eigenfunctions for each subspace i = 1 , Γ , Γ λ ,λ = Ψ (1) , Γ λ ⊗ Ψ (2) , Γ λ , which is not irreducible and has to be further symmetrized. We use the same projec-tion/transfer operator approach described above (and even the same numerical subroutine)24able 3: Term energies ˜ E λ (in cm − ), symmetries Γ, degeneracies l λ , TROVE assignments( v i ), and normal-mode assignments ( ν i ) of the eigensolutions for each subspace i =1, 2.Subspace i λ i ˜ E λ i l λ i Γ v v v v v v State1 1 0.0 1 A A ν E ν A A ν E ν by Eqs. (32, 33). The required transformation matrices are obtained as products of thestandard irreducible transformation matrices D Γ , Γ [ R ] = D Γ [ R ] D Γ [ R ]which are well known and also programmed in TROVE for most symmetry groups. Usingstandard transformation matrices is numerically more stable compared to the procedurebased on the matrices D [ R ] evaluated directly as solutions of Eq. (27). This is exhibitedin significantly smaller errors in the computed coefficients a i , which are very close to beingintegral.To illustrate this point, it is informative to look at the product of two degenerate functionsΨ (1) ,E ⊗ Ψ (2) ,E as an example (see Table 3). The four components of the product Ψ (1) ,E ,n Ψ (2) ,E ,m ( n, m =1,2) transform as a direct product of two E -representation matrices D Γ s [ R ] = D E [ R ] ⊗ D E [ R ] . The characters are defined by χ Γ s [ R ] = ( χ E ) and give 4, 1, and 0 for E , (123), and (12) ∗ , respectively, which is the standard textbookexample of a reduction of the E ⊗ E product (see, for instance, Ref. 43). The reduction25oefficients as obtained from Eq. (30) are 1, 1, 1 for A , A , and E , i.e. E ⊗ E = A ⊕ A ⊕ E. The irreducible representations determined using the numerical approach described aboveare Ψ A , = 1 √ h Ψ ( E a ,E a )3 , + Ψ ( E b ,E b )3 , i , (39)Ψ A , = 1 √ h Ψ ( E a ,E b )3 , − Ψ ( E b ,E a )3 , i , (40)Ψ E a , = 1 √ h Ψ ( E a ,E a )3 , + Ψ ( E b ,E b )3 , i , (41)Ψ E b , = 1 √ h Ψ ( E a ,E b )3 , − Ψ ( E b ,E a )3 , i , (42)where the corresponding expansion coefficients ± / √ example as well as the description of the TROVE numerical sym-metrization procedure. The approach is very robust and is applicable to any product-typebasis sets constructed from 1D functions provided the transformation rules for the coordi-nates are known. The most time-consuming part of our numerical implementation is thesampling procedure which relies on the random selection of points and can occasionally leadto poor solutions of Eq. (27) for the transformation matrices. Usually the calculations arequick (seconds) but sometimes they can take hours (remember this is a basis set initializationpart which has to be done only once). 26 .2 An XY molecule of D symmetry: Degenerate and redun-dant coordinates Let us consider a more complicated example of coordinate choice, where some of the coor-dinates transform according to two-fold irreducible representations. Such coordinates arecommonly used to describe the vibrations of non-rigid molecules. For example, the nuclearcoordinates of ammonia can be defined as: q = ∆ r (43) q = ∆ r (44) q = ∆ r (45) q = 1 √ α − ∆ α − ∆ α ] (46) q = 1 √ α − ∆ α ] (47) q = τ. (48)Here, r , r , r are the bond lengths, α , α , and α are the interbond angles and τ is theinversion ‘umbrella’ coordinate measuring the angle between a bond and the trisector (seeRef. 52 for example).In this case the vibrational modes span three subspaces, stretching { q , q , q } , bending { q , q } , and inversion { q } , which transform independently. The symmetry properties of thetwo bending modes special compared to those of the stretching and inversion modes, wherethe effect of the symmetry operations on the latter is just a permutation Rq i = q j , i, j = 1 , , , or a change of sign, Rq = ± q . q and q (which are based on threeredundant coordinates α , α , and α ) are mixed by the degenerate E -symmetry trans-formations: R q q = D E [ R ] q + D E [ R ] q D E [ R ] q + D E [ R ] q . The product-type primitive basis set for NH ( J = 0) is φ ν ( Q ) = φ n ( q ) φ n ( q ) φ n ( q ) φ n ( q ) φ n ( q ) φ n ( q ) , (49)where φ n k ( q k ) ≡ | n k i ( k = 1 . . .
6) are 1D primitive basis functions. Due to the 2D characterof the transformations of q and q , the primitive bending functions φ n ( q ) and φ n ( q ) donot follow simple permutation symmetric properties. For example, by applying the (123)permutation to the product φ n ( q ) φ n ( q ) we get:(123) φ n ( q ) φ n ( q ) = φ n ( − q + √ q ) φ n ( − √ q − q ) , which cannot be expressed in terms of products of φ n ( q ) and φ n ( q ) only. Strictly speaking,an infinite primitive basis set expansion in terms of φ n ( q ) φ n ( q ) is required to represent Rφ n ( q ) φ n ( q ) exactly, except for the special case of Harmonic oscillator functions (seeSection 5.3). In practice, we use expansions large enough to converge the symmetrizationerror below the defined threshold of 10 − . Unlike the two above examples of rigid molecules,the lack of the permutation character of the product-type basis set φ n ,...,n ( Q ) in Eq. (49) alsoprevents its symmetrization using the transformation properties of the functions. However,our approach is based on the transformation properties of the coordinates Q , not functions,which allows a symmetry adapted representation to be constructed even in this case.The first step is to build three reduced Hamiltonian operators for each i = 1 , , H (1) ( q , q , q ) = h |h |h | ˆ H | i| i| i , (50)ˆ H (2) ( q , q ) = h |h |h |h | ˆ H | i| i| i| i , (51)ˆ H (3) ( q ) = h |h |h |h |h | ˆ H | i| i| i| i| i (52)and solve the corresponding eigenvalue problemsˆ H ( i ) ( Q ( i ) )Ψ ( i ) λ i ( Q ( i ) ) = E λ i Ψ ( i ) λ i ( Q ( i ) ) . (53)As discussed above, we expect all eigenvectors of Eq. (53) to transform according to the irre-ducible representations Γ ( i ) s of D (M) despite the non-permutative character of the bendingprimitive functions. It should be noted that in practical calculations, employing a finite basisset affects the accuracy with which the irreducible character of the eigenfunctions can bedetermined, which is particularly true for high vibrational excitations n k .To illustrate this, let us consider a generic variational calculation of several lower eigen-states for ammonia. Here we use the PES from Ref. 12 and the primitive basis set definedby a polyad number P of P = 2( n + n + n ) + n + n + n / ≤ P max = 28 . (54)The primitive basis functions φ n k ( q k ) ( k = 1 ..
6) are obtained as eigensolutions of the cor-responding 1D reduced Hamiltonians using the Numerov-Cooley technique with a com-putational setup as described in Ref. 12. The solution of the reduced stretching problem inEq. (50) is equivalent to the example of the rigid XY example detailed above (see Eqs. (22))and is not discussed further. The first three solutions of the bending reduced problem in29q. (50) are given by (where | n , n i ≡ | n i| n i )Ψ (2)1 = 0 . | , i − . | , i + | , i ) + . . . (55)Ψ (2)2 , = 0 . | , i + 0 . | , i + . . . (56)Ψ (2)2 , = − . | , i + 0 . | , i + . . . . (57)with the energy term values of 0.0, 1679.6324 and 1679.6324 cm − relative to the ZPE =1953.7381 cm − .The wavefunctions Ψ (2)2 , and Ψ (2)2 , are recognized as degenerate ( λ = 2) due to their verysimilar energies (we use a threshold of 10 − cm − ) and should be processed together at thesymmetrization step. The transformation matrices D [ R ] are determined by sampling theeigenfunctions on a grid of 40 points to give the reduction coefficients a i = 0, 0, 1, 0, 0 and0 ( ± − ) for A ′ , A ′ , E ′ , A ′′ , A ′′ and E ′′ , respectively. The projection operator procedureleads to the symmetrized combinations given byΨ (2)2 , = − . | , i − . | , i − . | , i + . . . (58)Ψ (2)2 , = − . | , i + 0 . | , i + 0 . | , i + . . . (59)Reducing the basis set to P max = 2, i.e. taking only n , n ≤
1, leads to similar solutionsbut with larger errors of about 10 − for a i , which is still rather small in this case. However,the wavefunctions corresponding to higher excitations will introduce larger errors and willrequire more basis functions for accurate symmetrization. We use a threshold of 10 − –10 − for reduction coefficients a i to control the symmetrization procedure: the program will acceptsolutions if a i differ from an integer by less than this value.As a final and conclusive test, TROVE also checks the matrix elements of the totalHamiltonian ˆ H between different symmetries, which should be vanishingly small to allow ablock-diagonal form of the Hamiltonian matrix. TROVE uses an acceptance threshold of300 − –10 − cm − to control the quality of the symmetrization procedure. Failure to pass thistest (usually small errors) indicates that the basis set is not large enough for an accuratesymmetrization. Critical failure (huge errors) usually means problems with the model (e.g. inthe potential energy function, coordinate transformation relations, kinetic energy operator,definition of the molecular equilibrium structure etc).This example is a good illustration of how the redundant coordinates can be incorporatedinto a product-type basis of 1D functions. The redundant vibrational coordinates are verycommon, for example they appear as part of multi-dimensional symmetrically adapted coor-dinates. The typical example are the bending modes used to represent vibrational modes ofammonia, Eqs. (46,47) or methane (see, for example, Ref. 53). The TROVE symmetrizationcan still handle this situation even at a cost of a larger basis set. As it will be shown in thenext section, the Harmonic oscillator basis functions have the property of their products toform symmetrized combinations from a finite size basis of functions, which holds also for thecase of the redundant coordinates. Our most common choice of the primitive basis set is based on the Numerov-Cooley approach,where 1D functions are generated numerically on a large grid of 1000–5000 equidistantlyplaced points by solving a set of 1D reduced Hamiltonian problems for each mode. Thisprovides a compact basis set optimized for a specific problem. However, as was discussedin the previous section, some types of degenerate coordinates require large expansions interms of products of 1D functions for accurate symmetrization. A very simple work-aroundof this problem is to use 1D Harmonic oscillators as a basis set. The (degenerate) Harmonicoscillators have a unique property: one can always build an isotropic Harmonic oscillatorwith proper symmetric properties as a finite sum of products of 1D Harmonic oscillators φ HO n i ( q i ) (see, for example, Ref. 43). This is also valid for higher milti-fold degeneracies.As an illustration, in the Appendix we show how to construct a 2D symmetrized basis31et using 1D Harmonic oscillator functions to represent the asymmetric bending modes ofthe ammonia molecule using our symmetrization procedure. In fact this illustration canbe reproduced without the TROVE program as it is solely based on the properties of theHarmonic wavefunctions. This makes up a good toy example to try our symmetrizationapproach without having to deal with TROVE implementation.It should be noted that the eigenfunction methods for many-particle harmonic oscillatorwavefunctions was also explored by Novoselsky and Katriel . TROVE uses the rigid rotor wavefunctions (Wigner D -functions) as the rotational basis set.In principle, for most of the groups (such as C n v , C n h , D n h , or D n d ) the symmetry propertiesof the rigid rotor wavefunctions | J, k, m i are trivial and can be reconstructed based on the k value only. This is possible because all symmetry operations from these groups can beassociated with some equivalent rotations about the body-fixed axes x , y , and z only (see, forexample, the discussion in Ref. 55). Furthermore, symmetrized combinations of the rigid-rotor wavefunctions are trivial and can be given by the so-called Wang wavefunctions. Forexample, TROVE uses the following symmetrization scheme: | J, , τ rot i = | J, k, m i , (60) | J, K, τ rot i = i σ √ (cid:2) | J, K, m i + ( − J + K + τ rot | J, − K, m i (cid:3) , (61)where K = | k | , τ rot is the value associated with the parity of | J, K, τ rot i , σ = K mod 3 for τ rot = 1, σ = 0 for τ rot = 0, and m is omitted on the left-hand side for simplicity’s sake.The symmetry properties of | J, K, τ rot i can be derived from the properties of | J, k, m i underthe associated rotations and depend on J, K and τ rot only. Therefore a more sophisticatedsymmetrization approach like the one presented above is not required in such cases. As anexample, Table 4 lists the symmetries of | J, K, τ rot i for a rigid XY -type molecule ( C (M))32escribed above.Table 4: C (M) symmetries of the rigid-rotor wavefunctions | J, K, τ rot i ( K ≥
0) for the caseof a rigid XY molecule. K = 0 is the special case with τ rot = 0 (even J ) and τ rot = 1 (odd J ). Γ K τ rot A n A n E a n ± E b n ± x , y , and z , such as T d and O h . Consider a rigid XY molecule spanning the T d (M) symmetry group. The permutation (124) is associated with the equivalent rotation R (1 , , π/ x, y, z ) = (1 , , In this case the symmetrized basis can only be formed from a linearcombination of | J, k, m i spanning a range of k values, as was also shown by Alvarez-Bajoet al. . This is where we use the TROVE symmetrization approach to build symmetryadapted rotational basis functions | J, Γ i (see also Refs. 2,56, where this approach was appliedfor J up to 45). The formulation of the symmetrization scheme is given in the Appendix.Once the D Wang [ R ] matrices are known, the numerically-adapted reduction scheme de-scribed above is used to build the symmetrized representation for any J . The rotationalquantum number K cannot be used for classification of these symmetrized rigid-rotor com-binations anymore. Instead we label them as | J, Γ , n i , where n is a counting index. Following the subspace-based approach introduced for symmetrization of the vibrationalpart, the rotational modes are also treated as part of an independent, rotational subspace,which is referred in TROVE to as subspace 0. The symmetry adapted ro-vibrational basisset is then constructed as a direct product of the symmetrized components from different33ubspaces as Ψ (0) , Γ λ ⊗ Ψ (1) , Γ λ ⊗ Ψ (2) , Γ λ · · · ⊗ Ψ ( L ) , Γ L λ L , where L is the number of vibrationalsubspaces. The product of irreducible representations must be further reduced, which ismuch easier when each component is transformed as one of the irreps of the group with stan-dard transformation properties. In this case the same projection operator symmetrizationtechnique is used without further sampling of the symmetric properties of the correspondingcomponents.An efficient alternative to the vibrational basis set as a direct product of subspaces is the J = 0 contraction scheme. According to this scheme the eigenfunctions of the vibrational( J = 0) problem are used as contracted vibrational basis functions for J >
0. The J = 0eigenfunctions represent an even more compact basis set and can be further contracted(referred to as the J = 0 contraction). The symmetry adapted ro-vibrational basis setis then constructed exactly as described above (using the same numerical symmetrizationsubroutines) as a direct product of Φ (0) , Γ λ ⊗ Φ (1) , Γ λ , where the subspace-index i in Φ ( i ) , Γ i λ i runsover 0 and 1 only, and the J = 0 basis functions are combined into subspace 1. A new method for constructing symmetry adapted basis sets for ro-vibrational calculationshas been presented. The method is a variation of the matrix (or eigenfunction) approachesand is based on solving eigenfunction problems for a set of reduced Hamiltonian operatorswithout resorting to rigorous group-theoretical algebra. The advantage of using reducedHamiltonians in the matrix symmetrization is that it also improves the properties of thebasis sets by making them more compact and adjusted to the physics of the problem, thusallowing for efficient contraction. However, it lacks the automatic classification of the basisfunctions by the irreps, which is a useful feature of the CSCO-based eigenfunction approachby Chen et al. . To make up for this, the TROVE symmetrization procedure has to becomplemented by a sampling technique accompanied by a projection-based reduction.34ur symmetrization approach has been implemented in the TROVE program suite andhas been extensively used for a variety of tri-, tetra-, and penta-atomics covering the C s (M), C (M), C (M), D (M), D (M), and T d (M) groups. TROVE symmetrization is general inthat it can be applied to any molecule with arbitrary selection of coordinates provided thesymmetry properties of the latter are known. We are now implementing a general numericaltechnique for C n h (M), C n v (M), and D n h (M) representations, where n is an arbitrarily largeinteger value. Although TROVE symmetrization was developed and used for building ro-vibrational basis sets, we believe it can be useful for many other applications in physicsand chemistry. The symmetrization subroutines (Fortran 95) are written to be as generalas possible and in principle can be interfaced with other variational codes, if there will beinterest from the community. The illustration of the symmetrization approach applied tothe Harmonic oscillator wavefunctions (see Appendix) is an example where using TROVEis not necessary and thus could be a good place to start. Acknowledgements
This work was supported by the ERC Advanced Investigator Project 267219. We also thankthe support of the COST action MOLIM No. CM1405, UCL for use of the Legion HighPerformance Computer and DiRAC@Darwin HPC cluster. DiRAC is the UK HPC facilityfor particle physics, astrophysics, and cosmology and is supported by STFC and BIS. SNYand RIO thank Per Jensen for very helpful discussions and inspiration. The work of RIOwas supported by RFBR No. 15-02-07473, 15-02-07887, and 16-32-00668. We also thankAlec Owens for proofreading the manuscript. His valuable comments and suggestions led toa significant improvement of this work and are greatly appreciated.35
Symmetrized 2D Harmonic oscillator basis sets
Here we illustrate how to build a D (M) symmetrized vibrational basis set for ammonia-type molecules to represent the two asymmetric bending modes q and q from subspace 2(see Eqs.(46,47)) for the example from Section 5.2. The basis set will be formed from theproducts of the degenerate Harmonic oscillator basis set functions as given by: φ HO n ,n ( q , q ) = C n ,n H n ( q ) e − αq H n ( q ) e − αq . where H n is a Hermite polynomial, α is a parameter (the same for all degenerate compo-nents), and C n ,n is a normalization constant. These functions represent solutions of the 2Ddegenerate Harmonic oscillator H HO φ HO n ,n ( q , q ) = ˜ ω ( n + n + 1) φ HO n ,n ( q , q ) , and can be combined to express a solution of the 2D isotropic Harmonic oscillator (IHO):Ψ IHO
N,l = F N,l ( ρ ) e ilφ , where ρ = q q + q , φ = arctan q q ,N = n + n , l = N, N − , ... − N + 2 , − N. Here Ψ
IHO
N,l is an eigenfunction of the corresponding 2D IHO problem: H HO Ψ IHO
N,l = ˜ ω ( N + 1)Ψ IHO
N,l , which transforms as A /A (for l = 0 , , , . . . ) and E (otherwise). That is, there exists alinear transformation that connects Ψ N,l and φ n ,n ( q , q ) subject to N = n + n .36n order to find such a transformation and thus build the symmetry adapted functionsΨ N,l , we apply the TROVE numerical symmetrization procedure. For example, for a givenpolyad number N = 3, we need to combine the following four products φ n ( q ) φ n ( q ) satis-fying n + n = 3: φ φ , φ φ , φ φ , φ φ . (62)These four wavefunctions are degenerate and share the same Harmonic oscillator energy ˜ E HO n ,n = ˜ ω ( n + n + 1) , with ˜ ω = 1679 .
380 cm − and α = 0 . − . We use a sampling grid of 40 geometriesranging between − . ≤ q , q ≤ . × D [ R ] for eachoperation R in D (M).Applying the group operations E , (123), (23), E ∗ , (132) ∗ , and (23) ∗ to the four selecteddegenerate wavefunctions φ = | , i , φ = | , i , φ = | , i , φ = | , i , the following matri-ces D [ R ] are obtained: D [ E ] = D [ E ∗ ] = , D [(23)] = D [(23) ∗ ] − − , D [(132)] = D [(132) ∗ ] = − − − √ − √ √ − √ − √ − √
34 58 − √ − √
38 38 − , where the matrix elements are given to 10 − . With the help of Eq. (31) we obtain the37eduction coefficients a i = 1 , , , , , − ) ( i = A ′ , A ′ , E ′ , A ′′ , A ′′ , E ′′ ), i.e.only A ′ , A ′ , E ′ combinations can be formed. For the 1D representations A and A theprojection operators obtained using Eq. (32) are given by P A = − √ − √ , (63) P A = − √
00 0 0 0 − √
00 0 0 0 . (64)These matrices contain a total of eight vectors that we can choose from to build the ir-reducible combinations of φ , φ , φ and φ , four of which are trivial and only two pairsare linearly independent. Choosing the second column vector from the P A matrix, afternormalization we obtain Ψ A ( q , q ) = √ | i| i − | i| i , where the index 1 indicates the counting number of this state. The only non-trivial andlinearly independent choice is given by (after normalization):Ψ A ( q , q ) = 12 | i| i − √ | i| i . E -symmetry component leads to the following matrices: P E = √ √ , (65) P E = √
00 0 0 0 √
00 0 0 0 . (66)From this space of eight possible vectors we select the following trial vector, which is linearindependent from Ψ A and Ψ A : φ trial = (0 , , , √
34 ) T . This vector is then orthogonalized to Ψ A and Ψ A using the Grand-Schmidt procedure toget Ψ E a , ( q , q ) = 12 | i| i + √ | i| i . The second component of the latter vector, E b , is obtained by applying the correspondingtransfer operator in Eq. (33) to Ψ E a , , which after the Grand-Schmidt orthogonalization reads:Ψ E b , ( q , q ) = − √ | i| i − | i| i . In all these equations above, the coefficients 1 /
2, 3 / √ / √ / − .This completes the example of the reduction of degenerate wavefunctions in the basis of2D isotropic Harmonic oscillator functions. It is also a good illustration of how our method is39pplied to degeneracies of arbitrary order and dimensions. The only limitation is the memoryand time required for sampling wavefunctions and inverting the transformation matrices viaEq. (25). For example in Ref. 1 the highest degeneracy order used was 120 together withthe 2D and 3D symmetries E , F and F .It should be noted that as an alternative to the orthogonalization with Grand-Schmidt,one can impose the orthogonality conditions on the elements of matrix D [ R ] during thesolution of Eq. (27), for example by utilizing the exponential ansatz: D [ R ] = exp( − κ [ R ]) , κ T = − κ , (67)with only one independent element κ in case of a doubly-degenerate irrep, three indepen-dent elements κ , κ , and κ in case of a triply-degenerate irrep, etc. to be determined.Using this representation, the system of equations in Eq. (27) becomes nonlinear and canbe easily solved using the iterative approach described in Ref. 11 for solution of the Eckartequations. We are planning to explore this approach in the future. B Symmetrization of the rotational rigid rotor wave-functions
The symmetry transformation properties (i.e. transformation matrices D [ R ]) of the rigid-rotor wavefunctions required for our symmetrization scheme can be obtained directly fromthe Wigner D -functions, associated with the corresponding Euler angles of the particularequivalent rotation R ( α, β, γ ). These are given by (see also ): R ( α, β, γ ) | J, k, m i = J X k ′ = − J D ( J ) ∗ k ′ ,k ( α, β, γ ) | J, k ′ , m i . molecule ( T d (M)) are given in Table 4of Ref. 41. The transformation properties of the Wang functions | J, K, τ rot i in Eqs. (60,61)can then be deduced using the unitary transformation D Wang [ R ] = U + D ( J ) ∗ ( α, β, γ ) U , where the (2 J + 1) × (2 J + 1) matrix U i,j is given by U , = , even J, − ( i ) J , odd J and U n,n = 1 √ , (68) U n,n +1 = − i ( − σ √ , (69) U n +1 ,n = ( − J + K √ , (70) U n +1 ,n +1 = i ( − J + K + σ √ , (71) U n,n ′ = 0 , for | n − n ′ | > , (72)where n = 2 K , K = 1 . . . J , σ = K mod 3 for τ rot = 1, and σ = 0 for τ rot = 0. Oncethe transformation matrices are known, the standard projection technique described aboveis applied to obtain the symmetry adapted rigid-rotor combinations used as rotational basisfunctions . 41 eferences (1) Yurchenko, S. N.; Tennyson, J.; Barber, R. J.; Thiel, W. Vibrational transition momentsof CH from first principles. J. Mol. Spectrosc. , , 69–76.(2) Yurchenko, S. N.; Tennyson, J. ExoMol line lists - IV. The rotation-vibration spectrumof methane up to 1500 K. Mon. Not. R. Astron. Soc. , , 1649–1661.(3) Nikitin, A. V.; Rey, M.; Tyuterev, V. G. An efficient method for energy levels calculationusing full symmetry and exact kinetic energy operator: Tetrahedral molecules. J. Chem.Phys. , , 094118.(4) ˇCejchan, A.; Spirko, V. Transforming from internal coordinates to Cartesian displace-ments in the Eckart frame: a Taylor series expansion approach. J. Mol. Spectrosc. , , 142–145.(5) Boudon, V.; Champion, J.-P.; Gabard, T.; Lote, M.; Michelot, F.; Pierre, G.; Rot-ger, M.; Wenger, C.; Rey, M. Symmetry-adapted tensorial formalism to model rovibra-tional and rovibronic spectra of molecules pertaining to various point groups. J. Mol.Spectrosc. , , 620 – 634.(6) Yurchenko, S. N.; Carvajal, M.; Jensen, P.; Lin, H.; Zheng, J. J.; Thiel, W. Rotation-vibration motion of pyramidal XY molecules described in the Eckart frame: Theoryand application to NH . Mol. Phys. , , 359–378.(7) Yurchenko, S. N.; Thiel, W.; Jensen, P. Theoretical ROVibrational Energies (TROVE):A robust numerical approach to the calculation of rovibrational energies for polyatomicmolecules. J. Mol. Spectrosc. , , 126–140.(8) Pavlyuchko, A. I.; Yurchenko, S. N.; Tennyson, J. A hybrid variational-perturbationcalculation of the ro-vibrational spectrum of nitric acid. J. Chem. Phys. , , 14.429) Cassam-Chenai, P.; Rousseau, G.; Ilmane, A.; Bouret, Y.; Rey, M. Application of quasi-degenerate perturbation theory to the calculation of rotational energy levels of methanevibrational polyads. J. Chem. Phys. , , 034107.(10) F´abri, C.; Quack, M.; Cs´asz´ar, A. G. On the use of nonrigid-molecular symmetry innuclear motion computations employing a discrete variable representation: a case studyof the bending energy levels of CH +5 . J. Chem. Phys. ,(11) Yachmenev, A.; Yurchenko, S. N. Automatic differentiation method for numerical con-struction of the rotational-vibrational Hamiltonian as a power series in the curvilinearinternal coordinates using the Eckart frame.
J. Chem. Phys. , , 014105.(12) Yurchenko, S. N.; Barber, R. J.; Yachmenev, A.; Thiel, W.; Jensen, P.; Tennyson, J.A variationally computed T =300 K line list for NH . J. Phys. Chem. A , ,11845–11855.(13) Yachmenev, A.; Yurchenko, S. N.; Jensen, P.; Thiel, W. A new ”spectroscopic” potentialenergy surface for formaldehyde in its ground electronic state. J. Chem. Phys. , , 11.(14) Sousa-Silva, C.; Hesketh, N.; Yurchenko, S. N.; Hill, C.; Tennyson, J. High temperaturepartition functions and thermodynamic data for ammonia and phosphine. J. Quant.Spectrosc. Radiat. Transf. , , 66–74.(15) Sousa-Silva, C.; Al-Refaie, A. F.; Tennyson, J.; Yurchenko, S. N. ExoMol line lists VIII:A Hot Line List for Phosphine. Mon. Not. R. Astron. Soc. , , 2337–2347.(16) Underwood, D. S.; Yurchenko, S. N.; Tennyson, J.; Jensen, P. Rotational spectrum ofSO and theoretical evidence for the formation of sixfold rotational energy-level clustersin its vibrational ground state. J. Chem. Phys. , , 244316.4317) Al-Refaie, A. F.; Yurchenko, S. N.; Yachmenev, A.; Tennyson, J. ExoMol line lists IX:A variationally computed line-list for hot formaldehyde. Mon. Not. R. Astron. Soc. ,(18) Al-Refaie, A. F.; Ovsyannikov, R. I.; Polyansky, O. L.; Yurchenko, S. N.; Tennyson, J.A variationally calculated room temperature line-list for H O . J. Mol. Spectrosc. , , 84–90.(19) Owens, A.; Yurchenko, S. N.; Yachmenev, A.; Tennyson, J.; Thiel, W. Accurate abinitio vibrational energies of methyl chloride. J. Chem. Phys. , , 244306.(20) Owens, A.; Yurchenko, S. N.; Yachmenev, A.; Thiel, W. A global potential energysurface and dipole moment surface for silane. J. Chem. Phys. , .(21) Adam, A. Y.; Yachmenev, A.; Yurchenko, S. N.; Jensen, P. Ro-vibrational averaging ofthe isotropic hyperfine coupling constant for the methyl radical. J. Chem. Phys. , .(22) Owens, A.; Yurchenko, S. N.; Thiel, W.; ˇSpirko, V. Accurate prediction of the ammoniaprobes of a variable proton-to-electron mass ratio. Mon. Not. R. Astron. Soc. , , 3191–3200.(23) Al-Refaie, A. F.; Polyansky, O. L.; Ovsyannikov, R. I.; Tennyson, J.; Yurchenko, S. N.ExoMol line lists XV. A new hot line list for hydrogen peroxide. Mon. Not. R. Astron.Soc. , , 1012–1022.(24) Underwood, D. S.; Tennyson, J.; Yurchenko, S. N.; Clausen, S.; Fateev, A. ExoMol linelists XVII: A line list for hot SO . Mon. Not. R. Astron. Soc. , , 4300–4313.(25) Owens, A.; Yurchenko, S. N.; Yachmenev, A.; Tennyson, J.; Thiel, W. A global abinitio dipole moment surface for methyl chloride. J. Quant. Spectrosc. Radiat. Transf. , , 100 – 110. 4426) Owens, A.; Yurchenko, S. N.; Yachmenev, A.; Tennyson, J.; Thiel, W. A highly accurateab initio potential energy surface for methane. The Journal of Chemical Physics , , 104305.(27) Tennyson, J.; Yurchenko, S. N. ExoMol: molecular line lists for exoplanet and otheratmospheres. Mon. Not. R. Astron. Soc. , , 21–33.(28) Tennyson, J.; Yurchenko, S. N. The ExoMol project: Software for computing largemolecular line lists. Intern. J. Quantum Chem. , , 92–103.(29) Al-Refaie, A. F.; Yurchenko, S. N.; Tennyson, J. GPU Accelerated INtensities MPI(GAIN-MPI): A new method of computing Einstein-A coefficients. Comput. Phys. Com-mun. , , 216–224.(30) Gabriel, J. New Methods for Reduction of Group Representations Using an Extensionof Schur’s Lemma. J. Math. Phys. , , 494–504.(31) Morozova, N. K.; Morozov, B. P. Dokl. Akad. Nauk USSR (in Russian) , ,817–820.(32) Morozova, N. K.; Morozov, B. P. Dokl. Akad. Nauk USSR (in Russian) , ,538–541.(33) Moccia, R. A numerical method to obtain a symmetry-adapted basis from the hamil-tonian or a similar matrix. Theor. Chem. Acc. , , 85–88.(34) Dellepiane, G.; Gussoni, M.; Zerbi, G. Symmetry Properties of the Kinetic EnergyMatrix and Their Applications to Problems of Molecular Vibrations. J. Chem. Phys. , , 3450–3452.(35) Chung, A. L. H.; Goodman, G. L. Computer Generation of Molecular Symmetry Or-bitals. J. Chem. Phys. , , 4125–4137.4536) Bouman, T. D.; Goodman, G. L. Computer Algorithms for Symmetry Adaptation: AGeneral Method for Molecular Point Groups. J. Chem. Phys. , , 2478–2479.(37) Jordanov, B.; Orville-Thomas, W. A computational method for symmetry factorizationwithout using symmetry operators. J. Molec. Struct. (THEOCHEM) , , 323 –327.(38) Chen, J.-Q.; Gao, M.-J.; Ma, G.-Q. The representation group and its application tospace groups. Rev. Mod. Phys. , , 211–278.(39) Katriel, J. Some useful results concerning the representation theory of the symmetricgroup. J. Phys. A: Math. Gen. , , 5227.(40) Lemus, R. A general method to obtain vibrational symmetry adapted bases in a localscheme. Mol. Phys. , , 2511–2528.(41) Alvarez-Bajo, O.; Lemus, R.; Carvajal, M.; Perez-Bernal, F. Equivalent rotations as-sociated with the permutation inversion group revisited: symmetry projection of therovibrational functions of methane. Mol. Phys. , , 797–812.(42) Lemus, R. Quantum Numbers and the Eigenfunction Approach to Obtain SymmetryAdapted Functions for Discrete Symmetries. Symmetry , , 667.(43) Bunker, P. R.; Jensen, P. Molecular Symmetry and Spectroscopy , 2nd ed.; NRC ResearchPress: Ottawa, 1998.(44) Hamermesh, M.
Group theory and its application to physical problems ; Dover: NewYork, NY, 1989.(45) Azzam, A. A. A.; Tennyson, J.; Yurchenko, S. N.; Naumenko, O. V. ExoMol molecularline lists XVI. The rotationvibration spectrum of hot H S. Mon. Not. R. Astron. Soc. , , 4063–4074. 4646) Noumeroff, B. Trudy Glavnoi Rossiiskoi Astrofizicheskoj Observatorii ; Moscow, Gosu-darsvennoe Izdatel’stvo, 1923; Vol. 2; pp 188–259.(47) Cooley, J. W. An Improved eigenvalue corrector formula for solving the Schr¨odingerequation for central fields.
Math. Comp. , , 363–374.(48) Yurchenko, S. N.; Carvajal, M.; Thiel, W.; Jensen, P. Ab initio dipole moment and the-oretical rovibrational intensities in the electronic ground state of PH . J. Mol. Spectrosc. , , 71–87.(49) Ovsyannikov, R. I.; Thiel, W.; Yurchenko, S. N.; Carvajal, M.; Jensen, P. Vibrationalenergies of PH calculated variationally at the complete basis set limit. J. Chem. Phys. , , 044309.(50) Sousa-Silva, C.; Yurchenko, S. N.; Tennyson, J. A computed room temperature line listfor phosphine. J. Mol. Spectrosc. , , 28–37.(51) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes –The Art of Scientific Computing , 3rd ed.; Cambridge University Press, 2007.(52) L´eonard, C.; Handy, N. C.; Carter, S.; Bowman, J. M. The vibrational levels of ammo-nia.
Spectrochimica Acta A , , 825 – 838.(53) Halonen, L. Internal coordinate Hamiltonian model for Fermi resonances and localmodes in methane. J. Chem. Phys. , , 831–845.(54) Novoselsky, A.; Katriel, J. Harmonic Oscillator SU States with Arbitrary Permuta-tional Symmetry.
Ann. Phys. , , 55 – 75.(55) Bunker, P. R.; Jensen, P. Spherical top molecules and the molecular symmetry group. Mol. Phys. , , 255–264. 4756) Yurchenko, S. N.; Tennyson, J.; Bailey, J.; Hollis, M. D. J.; Tinetti, G. Spectrum ofhot methane in astronomical objects using a comprehensive computed line list. Proc.Natl. Acad. Sci. U. S. A. ,111