Anharmonic excited state frequencies of para-difluorobenzene, toluene and catechol using analytic RI-CC2 second derivatives

Analytic second nuclear derivatives for excited electronic state energies have been implemented for the resolution-of-the-identity accelerated CC2, CIS(DN) and ADC(2) models. Our efficient implementation with O(N) memory demands enables the treatment of medium sized molecules with large basis sets and high numerical precision and thereby paves the way for semi-numerical evaluation of the higher-order derivatives required for anharmonic corrections to excited state vibrational frequencies. We compare CC2 harmonic and anharmonic excited state frequencies with experimental values for para-difluorobenzene, toluene and catechol. Basis set problems occur for out-of-plane bending vibrations due to intramolecular basis set superposition error. For non-planar molecules and in plane modes of planar molecules, the agreement between theory and experiment is better than 30 cm 1 on average and we reassign a number of experimental bands on the basis of the ab initio predictions.


Introduction
The characterisation of molecules in excited electronic states remains a challenge, both for experimental and theoretical chemistry. Electronic excitation is often accompanied by significant structural change and complex intramolecular vibrational energy redistribution processes, 1-3 resulting in rich spectra that are sometimes difficult to interpret. Better theoretical treatments of electronic excited states are key to understanding photochemical phenomena and ultimately to harnessing photochemistry as a route to controlling molecular bond fission processes. [4][5][6][7] Coupled cluster methods are among the most accurate ab initio electronic structure methods. 8 Coupled cluster methods for excited state properties have been developed extensively by Stanton and Gauss in the framework of equation of motion coupled cluster (EOM-CC) theory. [9][10][11] Benchmark studies, [12][13][14][15] found that EOM-CCSD ground and excited state harmonic frequencies agree with values derived from experiment with a root mean squared deviation (RMSD) of 20-30 cm À1 .
CC2 was designed 16 as approximation to CCSD with an O(N 5 ) scaling of the computational costs with system size N instead of O(N 6 ). CC2 conserves the order in the fluctuation potential through which single excitation dominated transitions are described correctly in CCSD theory. In the original implementations with exact four-index electron repulsion integrals the high prefactor of the computational costs and the high memory demands for transforming and storing the two-electron integrals severely limited system and basis set sizes for which calculations could be performed on commonly available computer hardware. This bottleneck is removed by combining CC2 with the resolution of the identity (RI) approximation for the two electron integrals. 17,18 The RI approximation accelerates the calculations by one to two orders of magnitude, depending on the basis set, reduces memory demands to O(N 2 ) and reduces disc storage demands to O(N 3 ). These computational savings make it possible to study large molecules such as chlorophylls, 19 which previously could only be studied with CIS or TDDFT.
The same technique can be used to accelerate other secondorder methods for excited states, e.g. configuration interaction singles with a perturbative doubles correction, 20 CIS(D), and its iterative variant 21 CIS(D N ), and the algebraic diagrammatic construction through second-order, 22 ADC(2), with similarly pivotal efficiency savings. Compared to the non-iterative perturbative CIS(D) approximation, the iterative methods CIS(D N ), ADC (2), and CC2 have the advantage that they provide a consistent description of excited state potential energy surfaces (PES), even in the region of avoided crossings, and can thus be used more straightforwardly for searching and characterizing stationary points on the excited state PES.
Implementations of analytic excited state gradients of these approximate second order methods have been reported in ref. 23 and 24, where it is demonstrated that memory demands for first derivatives scale at most as O(N 2 ) when using the RI approximation. Recently, second-order electronic response properties for ground and excited states and ground state nuclear Hessians as analytic second derivatives of the energy have been implemented. [25][26][27] To preserve both the high computational efficiency and the low storage demands for second derivatives the RI approximation is combined with a numerical Laplace transformation of orbital energy denominators 28 for the contribution of double excitation amplitudes to first derivatives of the density matrices.
In the current work, we extend the theory and implementation to analytic geometrical second derivatives for CC2, CIS(D N ) and ADC(2) excited state energies. Analytic second derivatives are particularly important for obtaining the high numerical accuracy required to calculate semi-numerical third and fourth derivatives for polyatomic molecules. Only recently has a similar route has been pursued for TDDFT. 29,30 With the implementation of analytic Hessians for CC2 and ADC(2) it becomes possible to compute anharmonic vibrational spectra of polyatomic molecules with a correlated ab initio wavefunction method. We demonstrate the applicability of our implementation by computing harmonic and anharmonic excited state frequencies for medium sized molecules, which we compare to experimentally observed band centres.

Excited state Hessian for RI-CC2
The theory and implementation of orbital-relaxed electric second-order response properties for excited states at the RI-CC2 level has been presented in ref. 26. In the current work we focus on the additional theory required for geometric second derivatives. We use identical notation to that of ref. 26 and 27 and, rather than repeat them here, refer the reader to that work for the full definition of all of the terms and intermediates.
In coupled cluster theory, properties of an excited state f can be obtained as derivatives of the excited state quasienergy Lagrangian 26 The Lagrangian is composed of the ground state energy E CC , the vector functions O m for the ground state cluster amplitudes t m , with Lagrange multipliers % t f m , and the excitation energies The fourth term is the orbital-rotation constraint, which imposes vanishing Fock matrix elements F pq for the relevant orbital pairs pq A m 0 . The last and second last terms determine the phase and normalisation of the left and right eigenvectors of the Jacobi matrix A, by coupling them to the eigenvectors L f,(0) and R f,(0) of the unperturbed limit. 26 The hessian of the excited state is obtained by differentiating the Lagrangian twice with respect to the nuclear positions. The result can be obtained in the same way as for excited state polarisabilities, presented in ref. 26 hĴ xy i ex is the expectation value of an effective second order Hamiltonian (vide infra) evaluated with the density for the excited state. The second and third terms depend on first derivatives of integrals and first derivatives of the cluster amplitudes (t y n ) and right eigenvectors R f, y n contained in the derivative one-and two-particle densities. The final two terms collect contributions that are bilinear in first derivatives of the cluster amplitudes and eigenvectors.
The only difference to polarizabilities are that for derivatives with respect to nuclear coordinates the AO basis and consequently all AO one-and two-electron integrals depend on the perturbation. This also gives rise to additional contributions to the derivatives of the MO coefficients related to the changes in the AO overlap. We write the first derivative of the molecular orbital (MO) coefficients with respect to nuclear displacements as where U x contains the derivative of the orbital rotation parameters j x and the derivative of the overlap matrix S [x] in the basis of the unmodified molecular orbitals (UMOs).
Here and throughout, a square bracket indicates a derivative of the UMO integrals, which includes the derivative of the AO integrals, but not the derivative of the MO coefficients. 31,32 The first term in eqn (2) For electric response properties only the second and third terms are present and H [ y] has no two-electron part. The round brackets indicate a one-index transformation of the integrals. 32 The evaluation of the expectation value hĴ xy i ex with O(N 2 ) memory demands will be described below.
The last two terms in eqn (2) depend only indirectly on the derivatives of the one-and two-electron integrals, through the derivatives of the cluster amplitudes, t x m , and right eigenvectors, R f,y n . These are evaluated as described in ref. 26 for electronic derivatives, with the difference that now the first derivatives of the AO two-electron integrals have to be included in the calculation of the amplitude and eigenvector derivatives. Detailed expressions for the derivatives of the cluster amplitudes are given in ref. 27 and the changes to the expressions for the derivatives of the eigenvectors are obtained in an analogous way.
The second and third terms in eqn (2) combine (via the densities) first derivatives of the amplitudes t y n and eigenvectors R f,y n for a coordinate y with first derivatives of the MO Fock matrix F x pq and two-electron integrals (pq|rs) x for a coordinate x. Here, as for all two-electron integrals in the correlation treatment, we employ the RI approximation with three index intermediates B Q;pq ¼ P P ðpqjPÞ V À 1 2 ! PQ composed of two-index V PQ = (P|Q) and three-index ( pq|P) electron repulsion integrals. The indices P and Q denote functions from an auxiliary basis set. For an efficient evaluation we rewrite the two-particle term as: X pqrsd nsep;ex;y pqrs where F eff,y pq is an effective Fock matrix (cf. ref. 26 and 27) and D Q,y ab and g y PQ are, respectively, a three-index two-particle density in the AO basis (indices a, b) and a two-index two-particle density: We highlight that D Q,y ab is evaluated without ever building the four-index two-particle density d nsep,y pqrs . Full working expressions for the effective Fock matrix and for the two-and three-index densities are provided in the ESI. † The one-electron densities D F,ex,y require contractions of the doubles parts of the left eigenvectors and the derivatives of the right eigenvectors of the form P abk L f jakb R f;y iakb and P ijc L f iajc R f;y ibjc . For undifferentiated eigenvectors and amplitudes the RI approximation is sufficient to implement these contractions efficiently with only O(N 2 ) memory demands and without storing doubles vectors as four-index quantities on disc. For differentiated eigenvectors, however, this is not the case. There are additional terms that are not simply dressed two-electron integrals, but come from contractions of the undifferentiated eigenvectors with derivatives of the Fock matrix: For full definitions of the intermediates see ref. 26. To convert the last three terms into an expression that is separable in the index pairs ai and bj we apply a numerical Laplace transformation of the orbital energy denominators where y m and w m are the integration points and weights, respectively. With this approximation the eigenvector elements can be written as with (14) and the transformations with the derivatives of the Fock matrix can be performed on the three-index intermediates. This makes it possible to evaluate the contributions from R f to the right hand side of R f,y within the same loops as the contributions from the dressed integrals, just with an additional summation over the Laplace grid points. We now turn to the evaluation of the expectation valueĴ xy of the second order Hamiltonian for an excited state. The contributions toĴ xy are grouped into three terms. The first term contains the second derivatives of the AO integrals for the Hamiltonian. It is evaluated in the same manner in the AO basis as the respective contribution from H [x] to the gradient described in ref. 23: ab are auxiliary three-index and g PQ and g [x] PQ two-index two-particle densities defined in the ESI. † All other contributions to hĴ xy i ex are rewritten as contractions of effective Fock matrices with derivatives of the overlap matrices and U x . This is done to avoid the evaluations of twoelectron integrals in the MO basis and with this any O(N 5 ) scaling steps that depend on two perturbations.
The second and third contribution toĴ xy are combined by introducing a modified first-order Hamiltonian, 27 Such expectation values of one-index transformed Hamiltonians can be evaluated as contraction of the transformation matrices U x with effective Fock matrices: The definition and calculation of the effective Fock matrices for the modified first-order Hamiltonian F eff and is evaluated as: where F eff,ex pq is the (unperturbed) excited state effective Fock matrix, also known as energy weighted density matrix. Its implementation for RI-CC2 has been described in ref. 23.
The above formulas have been implemented in the development version of the ricc2 program of the TURBOMOLE package. 33 All contributions are evaluated using integraldirect algorithms in the AO basis with O(N 2 ) memory demands and strictly avoiding any operations that scale as O(N 5 ) with the basis set size and at the same time quadratically with the number of perturbations, so that the computation of the full Hessian still scales only as O(N 6 ) with the system size. The time-determining steps are the solution of the linear equations for the first derivatives of the amplitudes and eigenvectors and the calculation of the first-order density matrices that have to be done for each perturbation.
The implementation of the excited state hessians and the excited state polarizabilities enables also the calculation of the derivatives of the (excited state) dipole moment as mixed derivatives by differentiating once with respect to the strength of an electric field and once with respect to the nuclear coordinates.

CIS(D N ) and ADC(2)
With only few modifications the implementation is easily adapted to CIS(D N ) and ADC (2). 21,22 In both methods the coupled cluster ground state amplitudes are replaced by the amplitudes from first-order Møller-Plesset perturbation theory (MPPT). We assume that the Brillouin condition is fulfilled and thus the singles ground state amplitudes t (0) and their derivatives t x vanish in first order MPPT, as do the singles Lagrange multipliers % t f,(0) . Discarding all contributions from these singles parameters turns the CC2 code into a CIS(D N ) code. As a side effect, the equations for the derivatives of the cluster amplitudes can be inverted directly, in the canonical implementation, bypassing the iterative solution procedure.
The implementation for ADC (2) is closely related to the implementation CIS(D N ). Here the secular matrix is symmetrised: 22,24 This only requires some small modifications in the singlessingles block of the Jacobian to symmetrise the contributions from hm 1 |[[H,T (1) 2 ],t n 1 ]|HFi and the corresponding contributions to the right hand side of the equation for R f,y and to the derivatives of the Jacobian. Due to the symmetric secular matrix, the left eigenvectors are identical with the right eigenvectors and only one set of eigenvalue equations has to be solved.

Thresholds and numerical accuracy
The numerical accuracy of the excited state Hessians, and also the computational costs, depend mainly on two thresholds: T Lap and T LRE . T Lap controls the numerical Laplace transformation such that T LRE is the threshold for the residual of the equations for the first derivatives of the cluster amplitudes and eigenvectors. The grid points for the Laplace transformation were optimised for the interval from e min = 2(e LUMO À e HOMO ) À o f and e max = 2(e vir,max À e occ,max ), where o f is the excitation energy, by minimizing F Lap with respect to y m and w m . We investigated the dependence of the results on the thresholds for a test set composed of the lowest excited states of glyoxal, methanethial, propinal, benzene and naphthalene at their equilibrium structures, the lowest two excited states of water and formaldehyde at the ground and the excited state equilibrium structure, and the third excited state of formaldehyde and the lowest two excited states of thiophene again at their equilibrium structures. For all 16 cases we computed the Hessian for the cc-pVDZ and the aug-cc-pVDZ basis sets with different values for both thresholds, T Lap and T LRE , and evaluated the deviations in the elements of the Hessian from a reference calculation with very tight thresholds. For the aug-cc-pVDZ result the mean absolute deviations of the Hessian elements follow roughly the relation D MAD Hess,CC2 r 5T Lap + 10T LRE (22) with the maximum deviations 20-50 times larger. For the cc-pVDZ basis set the deviations are between one and two orders of magnitude smaller. In general, we observe that, not unexpectedly, the thresholds have to be set about two orders of magnitude tighter than for ground state calculations and that tighter thresholds are needed if other states are close by. In this work we used a tight threshold for the response equations T LRE = 10 À10 and 11 grid points for the numerical Laplace integration, which corresponds to T Lap between 2 Â 10 À5 and 2 Â 10 À10 . These tight values were required because the third and fourth derivatives for VPT2 theory are obtained by numerical differentiation of the hessian. We kept the number of Laplace points the same for all hessians used in a finite difference formula to prevent any numerical noise from changes of the grid points.

Calculation of anharmonic corrections and frequencies
The analytical implementation of the hessians enables a seminumerical evaluation of third and fourth derivatives of the energy, which are sufficiently accurate to be used as cubic and quartic force constants for the calculation of anharmonic frequencies. The third and fourth derivatives are calculated by central finite differences of hessians. Since only semi-diagonal quartic derivatives are required for VPT2 theory, the total number of hessians required is N acc Á(3N À 6) + 1, where N acc is either two or four. N acc = 4 is required when four-and fivepoint formulas are used for cubic and quartic derivatives to reduce the error to O(d 4 ) in the displacement d. 34 We did not use the normal coordinates l m of the equilibrium structure for the displacement vectors since they are normalized in the mass weighted coordinate system, which leads to unbalanced step sizes. Instead, we used rescaled coordinates which ensures that the displacements become nearly independent of the atomic masses. Here, the vector l m denotes the mth eigenvector of the mass weighted hessian M À1/2 FM À1/2 and M is the diagonal matrix of the atomic masses, and m m 1 is the reduced mass of the normal mode m.
The anharmonic force constants are used for the calculation of anharmonic frequencies, by setting up the vibrational Hamiltonian as: 35 The zeroth order term H 0 vib is a quantum harmonic oscillator in reduced normal coordinates. The first-order correction H 1 vib consists of the cubic force constants f rst . The second-order correction H 2 vib contains the quartic force constants f rstu and a Coriolos term that depends on the equilibrium rotational constants B e a and the vibrational angular momenta j a . The anharmonic frequencies are calculated using second order vibrational perturbation theory (VPT2), as it is implemented in the program DYNAMOL. 36 Therein, the VPT2 equations and treatment of resonance effects follow the formulation described in the paper of Amos et al. 37

Results
We have used our new efficient CC2 excited state hessian implementation to construct harmonic and quartic force fields of toluene, para-difluorobenzene and catechol in their first excited electronic state, from which we obtain the band centres of vibrational transitions using VPT2. These medium sized molecules have been studied experimentally and we compare the computed and experimental wavenumbers, assessing the accuracy of the CC2 predictions and correcting assignments where appropriate. The nomenclature of the normal modes is adopted from the systematic studies of the Wright group for substituted benzenes, [38][39][40] which provides a more unique assignment than the Mulliken 41 or Wilson 42 nomenclature. For modes localised on the substituents, the common spectroscopic notation for stretches (n), bends (d) and torsions (t) is used. Outputs from all of our VPT2 calculations are provided in the ESI. † These contain a full list of fundamentals, overtones, combination bands, all Fermi and Darling-Dennyson resonances, and all effective Hamiltonians for the polyads used to treat these resonances variationally.

para-Difluorobenzene
The molecule para-difluorobenzene has been the focus of a considerable number of spectroscopic explorations to determine its ground and excited state structure and dynamics [43][44][45][46][47][48] including the intramolecular vibrational redistribution pathways in the excited state, which are mediated by rotational coupling and through Fermi resonance. 49,50 Here we examine the fundamental vibrational transitions of para-difluorobenzene in its first excited state 1 1 B 2u using CC2 theory.
Spectroscopic studies reveal that para-difluorobenzene exhibits D 2h symmetry in both its ground and first excited state, 45 and the CC2/cc-pVTZ optimised geometries indeed retain D 2h symmetry. The optimised structures are listed in the ESI. † The computed 0-0 transition is 4.67 eV, which is close to the experimentally determined energy of 4.57 eV. 46 Table 1 reports harmonic and anharmonic wavenumbers for the fundamental vibrational transitions of para-difluorobenzene in its first excited state, computed using CC2/cc-pVTZ and CC2/ cc-pVQZ levels of theory using the optimised CC2/cc-pVTZ excited state geometry. The harmonic wavenumbers computed using the cc-pV5Z basis set are also listed. Using the normal mode vectors, each transition has been assigned using the D i labeling convention of ref. 39 for para-di-substituted benzenes.
For ease of reference, the irreducible representation, Mulliken labeling and ground state experimental wavenumbers have also been included.
Examining the basis set convergence of the harmonic frequencies, we find that while the in plane modes converge rapidly with the basis set, the out of plane modes (a u , b 1g , b 2g and b 3u ) converge slowly, with deviations of more than 10 cm À1 between cc-pVTZ and cc-pV5Z values for the low frequency vibrations. As highlighted by Martin, Taylor and Lee for benzene and acetylene, 51,52 modes that break planarity suffer from basis set inconsistency errors that artificially lower the frequency and overestimate anharmonic terms in the quartic force field.
This artificial enhancement of the anharmonic couplings to out-of-plane vibrations complicates ab initio assignment the experimental bands. These difficulties notwithstanding, our calculations confirm the majority of the assignments of the 24 fundamental bands observed experimentally. The rightmost two columns of Table 1 list the previous assignment as collated in ref. 47 and our own ab initio assignment, respectively.
The fundamentals D 6 and D 27 were tentatively assigned to a b 1u band at 1015 cm À1 and a b 2u band at 1100 cm À1 , respectively, in gas phase two-photon spectroscopy measurements. 48 These assignments can be confidently discarded on the basis of our calculations and it is likely that these bands do not correspond to fundamental transitions. Similarly, Knight and Kable's tentative assignment of the b 3g band at 933 to the D 26 fundamental can also be discarded. There is a somewhat larger than expected discrepancy between the predicted and observed D 4 fundamental, which we cannot explain. Finally, we make a reassignment of the b 2u band at 1591 cm À1 to D 24 . The problems encountered with the low frequency out of plane modes prevent meaningful comment of the assignments of modes D 12 , D 13 , D 16 and D 17 , which is indicated by an asterix in the table. For the in plane modes, the overall agreement between anharmonic CC2/cc-pVQZ and experimentally observed transitions is very good, with an RMSD of 26 cm À1 .

Toluene
The excited state vibrational frequencies of toluene have been studied using dispersed fluorescence spectroscopy, 53 UV-IR double resonance spectroscopy, 54 one-colour resonance-enhanced multiphoton ionization and two-colour zero kinetic energy spectroscopy. 55 The CC2/cc-pVTZ optimised structures have C s symmetry in both the ground and first ( 1 B 2 ) excited states. In the ground state the carbon atoms fall in a plane, one hydrogen of the methyl group orientated perpendicular to the plane. In the excited state, the methyl group moves slightly out of plane and has a dihedral angle of 4 degrees with the benzene ring. The 0-0 transition energy for the first excited state is at 4.65 eV, 53 which is closely reproduced by CC2/cc-pVTZ theory, which yields 4.86 eV. Table 2 summarises the calculated and experimental frequencies of toluene in its first excited state. We list harmonic frequencies in cc-pVTZ, cc-pVQZ and cc-pV5Z basis sets, computed at the optimised cc-pVTZ structure, and VPT2 anharmonic frequencies using the cc-pVTZ basis set. The normal mode vectors were analysed and classified using the M i nomenclature of Gardner and Wright 38 for the ring modes, and given pseudo C 2v symmetry labels, where the methyl group is treated as a single pseudo-atom.
Concerning the basis set convergence of the harmonic wavenumbers, we find a similar pattern as for para-difluorobenzene. The in plane modes are converged to within 10 cm À1 with the cc-pVTZ basis, which is well below the intrinsic error bar of 30 cm À1 commonly ascribed to CC2 theory (for ground state frequencies) due to missing higher order correlation effects. As for para-difluorobenzene, the out-of-plane modes display a much slower basis set convergence, with differences of more than 30 cm À1 between cc-pVTZ and cc-pV5Z values for the torsion and for modes M 12 , M 15 and M 16 .
Overall, the CC2/cc-pVTZ anharmonic frequencies agree very well with the experimental band centres for the fundamental transitions. The only outlier is mode M 16 , which appears to have artificially enhanced positive anharmonic corrections due to the basis set incompleteness errors. By the same token, we expect the predicted fundamentals for modes M 12 and M 15 to also lie above the experimental bands, should they be measured in the future. The methyl internal rotation is not expected to be well described through VPT2 theory and requires a more advanced vibrational treatment using a potential energy surface that exhibits the three equivalent minima, see ref. 57-59 for  Table 2 and lie at higher frequency than the principle C-H stretches. Since we lack predicted intensities or experimental band shape information, we pragmatically assign the experimental bands in order of increasing frequency, which results in excellent agreement between computed and experimental values.
In Table 2 we also list selected overtones and combination bands for comparison with those determined experimentally, collected from ref. 53 Overall, the RMSD between observed and predicted band centres using VPT2 theory with a CC2/cc-pVTZ force field is 16 cm À1 for the fundamental transitions in the excited electronic state.

Catechol
Catechol (l,2-dihydroxybenzene) is biochemically important since the catecholamines adrenaline, noradrenaline, and dopamine are active in neurotransmission. It has been the subject of extensive spectroscopic studies, many of which have focused on the structure and dynamics of the low energy rotamers formed through changing the relative orientations of the hydroxyl groups, and the resulting differing levels of intra-and intermolecular hydrogen bonding. 60 The vibrational frequencies of the excited state have been probed using resonant two-photon ionization, fluorescence emission techniques, and molecular-beam hole-burning experiments. 61,62 Catechol is planar in the lowest energy isomer of the ground electronic state and the hydroxyl groups form an intermolecular hydrogen bond, but the structure of the excited state is assumed to be slightly distorted out of plane. 62 Our structural investigations, using CC2/cc-pVTZ theory to optimise geometries of the ground and excited states, confirm that the ground state is planar with an intermolecular hydrogen bond. We find that the excited state retains the intermolecular hydrogen bond, but is significantly distorted from planarity, puckering at the carbons with the hydroxyl groups. The optimised structures are reported in the ESI. † CC2/cc-pVTZ predicts the 0-0 transition at 4.53 eV which is in line with the experimentally determined 0-0 excitation energy of 4.42 eV. 62 The calculated and experimental frequencies of the ground and excited states are compiled in Table 3. Only 16 fundamentals of the excited state have been reported. The previous assignment was predicated on the assumption that the selection rules and band shapes for the planar S 0 state can be transferred to guide assignment of the S 1 state. 62,63 Our calculations indicate that this assumption was flawed and we report a completely fresh assignment for the observed S 1 vibrational bands. We therefore present frequencies also for the ground state, comparing our purely ab initio assignment procedure to the more comprehensive experimental assignments available for this state. [64][65][66][67] We use the o D i labeling for ortho-di-substituted benzene rings of ref. 40 and, following the convention of earlier works, 64 use n 1 , d 1 , t 1 to refer to the motions of the hydrogen donating OH group, and n 2 , d 2 , t 2 to refer to the motions of the hydrogen accepting OH group.
The procedure adopted for assigning the vibrational modes of the S 0 state was as follows: first the normal modes were analysed and categorised according to the o D i nomenclature, resulting in approximate C 2v symmetry labels a 1 , b 2 , b 1 , a 2 with corresponding type A, B, C line shapes for the IR active bands a 1 , b 2 , b 1 , respectively; these were then assigned in frequency order in effective C 2v symmetry blocks to the bands observed by Wilson for vapour phase catechol; 64 the only exception to this are the four a(CH) bending modes, where the a 1 and b 2 labels are swapped compared to Wilson's and where recent work casts doubt on the original assignment. 40 Where additional or more precise measurements are available, Wilson's band centres are replaced or supplemented by the modern values and the Raman active a 2 modes are assigned in frequency order using measurements from solid state catechol; the bands assigned to the OH bending (d), stretching (n) and torsional (t) modes are identified from examination of the normal mode coordinates.
The overall agreement with experiment and CC2/cc-pVTZ theory is excellent, with an RMSD of 29 cm À1 . The difficulties associated with out-of-plane vibrations are much less pronounced here and there is only one outlier of this type, the ring puckering mode D 26 . The only other significant outlier is mode D 9 , which is the mode distorting from a aromatic ring to three localised double bonds.
Turning now to the S 1 state, comparison with experiment is problematic. Only 16 frequencies have been assigned and the assignment appears flawed. A vibronic progression with a spacing of 113 cm À1 was observed in resonant two-photon ionisation spectra of catechol and was assigned to OH torsional overtones. 62 Our calculations do not predict low frequency torsional modes, but in fact predict a significant redshift of the torsional frequencies upon electronic excitation. The more tetrahedral arrangement at the oxygen centres in the excited state structure leads to a stronger intermolecular hydrogen bond, higher torsional frequencies and lower OH stretching frequencies. Instead, our calculations suggest assigning the progression of 113 cm À1 to the low frequency D 29 ring bend. Note that the selection rules based on spatial symmetry do not rigorously apply since the excited state structure is significantly distorted out of plane. Having discounted the presence of lowfrequency torsional modes, many of the spectral features observed in the experimental works must be reassigned. Proceeding to match computed and experimental values, accounting for band shape information where available, we report a fresh ab initio assignment in Table 3. The RMSD between VPT2 fundamentals using a CC2/cc-pVTZ force field and the observed bands with our fresh assignment is 22 cm À1 .

Conclusion
Excited state hessians have been implemented for CC2, with a focus on keeping the scaling of the main memory demands to at most O(N 2 ). This has been realized by exploiting the RI approximation for the two electron integrals and by choosing a Laplace decomposition of orbital energy denominators in the calculation of the first-order density matrices. The implementation is an extension of that for excited state polarizabilities and was extended straightforwardly to excited state hessians for the CIS(D N ) and ADC(2) methods. The code is parallelized with OpenMP and MPI to make use of modern computer hardware. The analytic implementation enables the semi-numerical calculation of third and fourth derivatives for anharmonic corrections.
We applied VPT2 theory based on CC2 quartic force fields to para-difluorobenzene, toluene and catechol and compared the computed frequencies to experimentally observed vibrational bands of the first excited states. In contrast to previous benchmark studies, we find that the accuracy of CC2 frequencies for excited electronic states is comparable to MP2 or CC2 frequencies for the ground electronic state, with typical average deviations between theory and experiment of less than 30 cm À1 for fundamental transitions when using a cc-pVTZ basis. However, we find that out-of-plane modes carry a much larger uncertainty due to internal basis set superposition errors, which leads to a strong basis set dependence of the force field terms, and unphysically large anharmonic corrections to typically underestimated harmonic frequencies.
In addition to assessing the accuracy of CC2 theory, our calculations have revealed some anomalies in the assignment of some of the experimentally observed bands. Our calculations discount some of the more tentative assignments in the spectrum of para-difluorobenzene. More significantly, our calculations indicate that the assignment of the lowest frequency features in the vibrational bands of catechol to OH torsions is incorrect, and that these lie much higher in energy due to the stronger hydrogen bond in the distorted S 1 excited state than in the planar S 0 state. Instead, we assign these features to ring modes, which become symmetry allowed transitions due to the low symmetry of the relaxed excited state structure. Due to this re-interpretation of the experimental bands it was necessary to perform a completely fresh assignment, and our new assignment can be considered a pure ab initio assignment.

Conflicts of interest
There are no conflicts to declare.