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

Abstract

Analytic second nuclear derivatives for excited electronic state energies have been implemented for the resolution-of-the-identity accelerated CC2, CIS(D_∞) and ADC(2) models. Our efficient implementation with [scr O, script letter O]([scr N, script letter 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⁻¹ on average and we reassign a number of experimental bands on the basis of the ab initio predictions.

---

1 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–7

Coupled cluster methods are among the most accurate ab initio electronic structure methods.⁸ 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–11 Benchmark studies,^12–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⁻¹.

CC2 was designed¹⁶ as approximation to CCSD with an [scr O, script letter O]([scr N, script letter N]⁵) scaling of the computational costs with system size [scr N, script letter N] instead of [scr O, script letter O]([scr N, script letter N]⁶). 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 [scr O, script letter O]([scr N, script letter N]²) and reduces disc storage demands to [scr O, script letter O]([scr N, script letter N]³). These computational savings make it possible to study large molecules such as chlorophylls,¹⁹ which previously could only be studied with CIS or TDDFT.

The same technique can be used to accelerate other second-order methods for excited states, e.g. configuration interaction singles with a perturbative doubles correction,²⁰ CIS(D), and its iterative variant²¹ CIS(D_∞), and the algebraic diagrammatic construction through second-order,²² ADC(2), with similarly pivotal efficiency savings. Compared to the non-iterative perturbative CIS(D) approximation, the iterative methods CIS(D_∞), 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 [scr O, script letter O]([scr N, script letter N]²) 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–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²⁸ 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_∞) 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.

2 Theory

2.1 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²⁶

The Lagrangian is composed of the ground state energy E_CC, the vector functions Ω_μ for the ground state cluster amplitudes t_μ, with Lagrange multipliers t̄^f_μ, 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 ∈ μ₀. 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.²⁶

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:

〈Ĵ^xy〉^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_ν) and right eigenvectors R^f,y_ν 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

C^x = C⁰U^x (3)

where U^x contains the derivative of the orbital rotation parameters κ^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) is the expectation value of an effective second-order Hamiltonian which collects several contributions to the second geometric derivatives of the Hamiltonian in the MO basis: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.³² The evaluation of the expectation value 〈Ĵ^xy〉^ex with [scr O, script letter O]([scr N, script letter N]²) 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_μ, and right eigenvectors, R^f,y_ν. 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_ν and eigenvectors R^f,y_ν 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 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:where F^eff,y_pq is an effective Fock matrix (cf.ref. 26 and 27) and Δ^Q,y_αβ and γ^y_PQ are, respectively, a three-index two-particle density in the AO basis (indices α, β) and a two-index two-particle density:We highlight that Δ^Q,y_αβ 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 and . For undifferentiated eigenvectors and amplitudes the RI approximation is sufficient to implement these contractions efficiently with only [scr O, script letter O]([scr N, script letter N]²) 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 denominatorswhere θ_m and w_m are the integration points and weights, respectively. With this approximation the eigenvector elements can be written aswith

K_Qm,ai = B̂_Q,aie^{(ε_i−ε_a)θ_m}, (13)

K̄_Qm,ai = B̄_Q,ai[R^f]e^{(ε_i−ε_a+ωf)θ_m}, (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:

D^eff,ex is the effective orbital-relaxed excited state and D^SCF the SCF one-particle density. h^[xy]_pq and A^CPHF,[xy]_pqrs are, respectively, the core Hamiltonian and CPHF matrices computed from the second derivatives of the AO integrals. Δ^Q_αβ and Δ^Q,[x]_αβ are auxiliary three-index and γ_PQ and γ^[x]_PQ two-index two-particle densities defined in the ESI.†

All other contributions to 〈Ĵ^xy〉^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 two-electron integrals in the MO basis and with this any [scr O, script letter O]([scr N, script letter N]⁵) scaling steps that depend on two perturbations.

The second and third contribution to Ĵ^xy are combined by introducing a modified first-order Hamiltonian,²⁷, to:

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 have been given in ref. 26 and 27. The last two contributions to Ĵ^xy are combined by introducing the second derivative of the overlap matrix S^xy = S^[xy] + P̂^xy(U^y,S^[x]) to

(S^[xy],Ĥ) + P̂^xy((U^y,S^[x]),Ĥ) = (S^xy,Ĥ) (18)

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.³³ All contributions are evaluated using integral-direct algorithms in the AO basis with [scr O, script letter O]([scr N, script letter N]²) memory demands and strictly avoiding any operations that scale as [scr O, script letter O]([scr N, script letter N]⁵) 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 [scr O, script letter O]([scr N, script letter N]⁶) 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.

2.2 CIS(D_∞) and ADC(2)

With only few modifications the implementation is easily adapted to CIS(D_∞) 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⁽⁰⁾ 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_∞) 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_∞). Here the secular matrix is symmetrised:^22,24

This only requires some small modifications in the singles-singles block of the Jacobian to symmetrise the contributions from 〈μ₁|[[H,T⁽¹⁾₂],τ_ν1]|HF〉 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.

2.3 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 thatT_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 ε_min = 2(ε_LUMO − ε_HOMO) − ω_f and ε_max = 2(ε_vir,max − ε_occ,max), where ω_f is the excitation energy, by minimizing F_Lap with respect to θ_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

Δ^MAD_Hess,CC2 ≤ 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⁻¹⁰ and 11 grid points for the numerical Laplace integration, which corresponds to T_Lap between 2 × 10⁻⁵ and 2 × 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.

2.4 Calculation of anharmonic corrections and frequencies

The analytical implementation of the hessians enables a semi-numerical 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 five-point formulas are used for cubic and quartic derivatives to reduce the error to [scr O, script letter O](δ⁴) in the displacement δ.³⁴

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/2FM^−1/2 and M is the diagonal matrix of the atomic masses, and 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:³⁵

H_vib = H⁰_vib + H¹_vib + H²_vib (24)

The zeroth order term H⁰_vib is a quantum harmonic oscillator in reduced normal coordinates. The first-order correction H¹_vib consists of the cubic force constants ϕ_rst. The second-order correction H²_vib contains the quartic force constants ϕ_rstu and a Coriolos term that depends on the equilibrium rotational constants B^e_α and the vibrational angular momenta j_α.

The anharmonic frequencies are calculated using second order vibrational perturbation theory (VPT2), as it is implemented in the program DYNAMOL.³⁶ Therein, the VPT2 equations and treatment of resonance effects follow the formulation described in the paper of Amos et al.³⁷

3 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–40 which provides a more unique assignment than the Mulliken⁴¹ or Wilson⁴² nomenclature. For modes localised on the substituents, the common spectroscopic notation for stretches (ν), bends (δ) and torsions (τ) 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.

3.1 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–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¹B_2u using CC2 theory.

Spectroscopic studies reveal that para-difluorobenzene exhibits D_2h symmetry in both its ground and first excited state,⁴⁵ 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.⁴⁶

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.

Table 1 CC2 excited state vibrational wavenumbers for para-difluorobenzene

^p D i	Sym.	Exp.	Harm.			VPT2		Exp.
		S0	TZ	QZ	5Z	TZ	QZ	S1	S1
20	30(b3u)	158	120	121	126	120	121	120^47	120
14	8(au)	422	148	147	161	169	168	175^47	175
19	17(b2g)	374	271	276	280	307	305	274^47	274
30	22(b2u)	348	352	353	352	349	350	352^47	352
29	27(b3g)	446	388	388	392	380	381	403^47	403
13	9(b1g)	800	404	406	411	453	476	588^46	*
11	6(ag)	450	411	412	423	407	408	410^47	410
18	29(b3u)	505	472	471	473	466	466	438^47	438
12	7(au)	945	486	490	503	612	676	583^47	*
17	16(b2g)	692	527	546	539	651	647	528^47	*
28	26(b3g)	635	558	558	558	541	542	558^47	558
15	15(b2g)	838	610	616	623	700	658	670^47	670
16	28(b3u)	928	614	621	624	636	813	619^47	*
10	14(b1u)	740	716	716	713	702	703	666^46	666
9	5(ag)	859	824	825	819	812	811	818^47	818
8	13(b1u)	1014	960	962	961	949	951	937^47	937
27	21(b2u)	1085	1025	1025	1022	1014	1006	1100^48
7	4(ag)	1140	1119	1117	1116	1100	1102	1116^46	1116
6	12(b1u)	1228	1228	1222	1218	1203	1195	1015^46
26	25(b3g)	1285	1254	1257	1256	1241	1244	933^46
5	3(ag)	1257	1273	1265	1260	1241	1234	1251^47	1251
25	20(b2u)	1306	1317	1319	1315	1297	1298	1591^48
4	11(b1u)	1514	1449	1444	1439	1412	1408	1335^46	1335
23	24(b3g)	1595	1510	1509	1497	1475	1475	1516^48	1516
3	2(ag)	1615	1557	1552	1540	1507	1503
24	19(b2u)	1437	1680	1677	1664	1635	1632		1591
22	23(b3g)	3085	3253	3252	3246	3118	3116
2	10(b1u)	3073	3257	3255	3249	3118	3116
21	18(b2u)	3091	3267	3265	3258	3107	3103
1	1(ag)	3088	3271	3268	3261	3125	3123

---

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⁻¹ 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₆ and D₂₇ were tentatively assigned to a b_1u band at 1015 cm⁻¹ and a b_2u band at 1100 cm⁻¹, respectively, in gas phase two-photon spectroscopy measurements.⁴⁸ 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₂₆ fundamental can also be discarded. There is a somewhat larger than expected discrepancy between the predicted and observed D₄ fundamental, which we cannot explain. Finally, we make a reassignment of the b_2u band at 1591 cm⁻¹ to D₂₄. The problems encountered with the low frequency out of plane modes prevent meaningful comment of the assignments of modes D₁₂, D₁₃, D₁₆ and D₁₇, 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⁻¹.

3.2 Toluene

The excited state vibrational frequencies of toluene have been studied using dispersed fluorescence spectroscopy,⁵³ UV-IR double resonance spectroscopy,⁵⁴ one-colour resonance-enhanced multiphoton ionization and two-colour zero kinetic energy spectroscopy.^55,56 21 of the 39 fundaments have been assigned, as have the methyl group internal rotor bands, five further overtone bands and several combination bands in the range 200–1500 cm⁻¹.

The CC2/cc-pVTZ optimised structures have C_s symmetry in both the ground and first (¹B₂) 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,⁵³ 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³⁸ for the ring modes, and given pseudo C_2v symmetry labels, where the methyl group is treated as a single pseudo-atom.

Table 2 CC2 excited state vibrational wavenumbers for toluene

M i	Sym.	Harm.			VPT2	Exp.
		TZ	QZ	5Z	TZ	S1
τ		86	118	122	28
20	b 1	141	148	151	140	145
14	a 2	234	238	238	236	226
19	b 1	319	325	325	331	314
30	b 2	323	328	329	325	331
18	b 1	432	441	445	460
11	a 1	449	450	451	442	457
29	b 2	524	524	526	512	532
13	a 2	539	544	560	563
17	b 1	556	560	566	552
16	b 1	639	651	669	734	697
12	a 2	658	667	686	736
15	b 1	717	733	747	839
10	a 1	751	749	749	737	753
28	b 2	932	932	932	915
8	a 1	943	939	938	928	934
9	a 1	969	969	972	957	966
δ as ^−		1003	1003	1002	972
δ as ^−		1048	1047	1047	1024	1021
27	b 2	1157	1155	1156	1137
7	a 1	1162	1160	1160	1143
6	a 1	1220	1214	1214	1196	1193
26	b 2	1302	1304	1305	1281
24	b 2	1386	1385	1386	1364
δ s ^+		1393	1394	1395	1358
5	a 1	1426	1423	1424	1394
δ as ^+		1474	1472	1472	1432
δ as ^+		1485	1484	1483	1443
23	b 2	1534	1526	1524	1482
4	a 1	1557	1547	1545	1505
25	b 2	1685	1668	1661	1635
ν s		3024	3015	3013	2921	2893
ν as		3098	3091	3088	2956	2956
ν as		3142	3135	3133	3013	2988
22	b 2	3209	3203	3202	3076	3048
3	a 1	3209	3204	3203	3077	3063
2	a 1	3218	3212	3210	3079	3077
21	b 2	3231	3224	3222	3087	3087
1	a 1	3240	3233	3231	3098	3097
Overtones
20	a 1	281	296	302	281	290
14	a 1	469	476	476	471	452
19	a 1	638	650	650	666	629
8	a 1	1885	1878	1876	1854	1868
9	a 1	1939	1939	1945	1914	1929
Combination bands
14 + 20	b 2	375	386	389	374	371
19 + 20	a 1	460	473	476	473	462
14 + 19	b 2	553	563	563	566	539
18 + 19	a 1	751	766	770	796	734
29 + 30	a 1	847	852	855	837	864
11 + 29	b 2	973	974	977	956	988
12 + 14	a 1	892	905	924	975	916
10 + 29	b 2	1275	1273	1275	1248	1263
8 + 11	a 1	1392	1389	1389	1370	1390
8 + 29	b 2	1467	1463	1464	1439	1463
9 + 11	a 1	1418	1419	1423	1399	1426
9 + 29	b 2	1493	1493	1498	1468	1494
8 + 9	a 1	1912	1908	1910	1885	1900
6 + 29	b 2	1743	1738	1740	1706	1727
23 + 25	a 1	3219	3194	3185	3121	3101
4 + 25	b 2	3241	3215	3206	3140

---

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⁻¹ with the cc-pVTZ basis, which is well below the intrinsic error bar of 30 cm⁻¹ 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⁻¹ between cc-pVTZ and cc-pV5Z values for the torsion and for modes M₁₂, M₁₅ and M₁₆.

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₁₆, 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₁₂ and M₁₅ 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 examples of such theories. The experimental assignment of the five aromatic C–H stretching frequencies M₁, M₂, M₃, M₂₁ and M₂₂, is complicated by the presence of Fermi resonances, which have not yet been untangled. Our calculations predict that M₂₁ is in resonance with M₄M₂₅ and that both M₁ and M₂ are in resonance with M₂₃M₂₅. The computed M₄M₂₅ and M₂₃M₂₅ combination bands are listed in 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, 54 and 56. The overtone and combination bands agree with experiment to within the error bar expected from the agreement found for the fundamentals. Some Fermi resonances observed in the spectra are not reproduced by our VPT2 calculations. Specifically, we do not find a resonance between M₁₁ and M₁₉M₂₀ and we find no indication that mode M₆ is in resonance with other states. Others, however, are identified by our calculations. The resonance between M₁₀ and M₁₈M₁₉ is found in both theory and experiment, but our calculations predict the fundamental M₁₀ to be lower in energy than M₁₈M₁₉, which if reliable would interchanging Garderner et al.'s assignment of the resonant pair, which has M₁₈M₁₉ lower in energy than M₁₀. Our calculations also verify the resonance between M₁₀M₂₉ and M₁₈M₁₉M₂₉.

Overall, the RMSD between observed and predicted band centres using VPT2 theory with a CC2/cc-pVTZ force field is 16 cm⁻¹ for the fundamental transitions in the excited electronic state.

3.3 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.⁶⁰ 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.⁶² 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.⁶²

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₀ state can be transferred to guide assignment of the S₁ state.^62,63 Our calculations indicate that this assumption was flawed and we report a completely fresh assignment for the observed S₁ 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–67 We use the ^oD_i labeling for ortho-di-substituted benzene rings of ref. 40 and, following the convention of earlier works,⁶⁴ use ν₁, δ₁, τ₁ to refer to the motions of the hydrogen donating OH group, and ν₂, δ₂, τ₂ to refer to the motions of the hydrogen accepting OH group.

Table 3 CC2/cc-pVTZ vibrational wavenumbers for catechol in the S₀ and S₁ states

^o D i	C 2v	S0			S1
		Harm.	VPT2	Exp.	Harm.	VPT2	Exp.
30	a 2	178	220	199	148	145
τ 2	a 2	209	254		475	425	461
29	b 1	293	295	299	130	127	113
21	a 1	303	300	320	295	293	299
τ 1	b 1	432	433		585	576	607
20	b 2	438	434	449	399	390	395
28	b 1	456	454	456	369	348	317
19	b 2	552	545	542	517	486	472
27	a 2	561	590	582	342	347
18	a 1	578	569	564	544	495	488
26	a 2	674	802	721	444	420
25	b 1	743	744	741	495	508	502
17	a 1	772	760	768	732	714	735
24	a 2	820	846	851	595	560	588
16	b 2	854	840	859	829	812	840
23	b 1	897	924	916	761	732	746
22	a 2	930	988	963	837	836	863
15	a 1	1043	1028	1030	968	950	950
14	b 2	1098	1079	1092	1044	1028	1061
13/δ2	a 1	1163	1143	1151	1149	1122
13/δ2	a 1	1168	1150	1151	1166	1159
12	b 2	1210	1190	1195	1122	1104
11	b 2	1268	1240	1251	1254	1220
10	a 1	1309	1278	1263	1286	1252
δ 1	b 2	1367	1338	1365	1343	1312
9	a 1	1458	1419	1324	1402	1368
8	b 2	1494	1455	1479	1450	1415
7	a 1	1540	1502	1504	1635	1581
6	b 2	1646	1601	1607	1511	1466
5	a 1	1651	1609	1616	1566	1517
4	b 2	3186	3048	3051	3177	3046
3	a 1	3211	3092	3051	3250	3102
2	b 2	3224	3101	3060	3190	3065
1	a 1	3235	3105	3081	3204	3065
ν 1	b 1	3750	3553	3605	3612	3390
ν 2	a 1	3813	3626	3663	3542	3252

---

The procedure adopted for assigning the vibrational modes of the S₀ state was as follows: first the normal modes were analysed and categorised according to the ^oD_i nomenclature, resulting in approximate C_2v symmetry labels a₁, b₂, b₁, a₂ with corresponding type A, B, C line shapes for the IR active bands a₁, b₂, b₁, respectively; these were then assigned in frequency order in effective C_2v symmetry blocks to the bands observed by Wilson for vapour phase catechol;⁶⁴ the only exception to this are the four α(CH) bending modes, where the a₁ and b₂ labels are swapped compared to Wilson's and where recent work casts doubt on the original assignment.⁴⁰ 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₂ modes are assigned in frequency order using measurements from solid state catechol; the bands assigned to the OH bending (δ), stretching (ν) and torsional (τ) 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⁻¹. 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₂₆. The only other significant outlier is mode D₉, which is the mode distorting from a aromatic ring to three localised double bonds.

Turning now to the S₁ 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⁻¹ was observed in resonant two-photon ionisation spectra of catechol and was assigned to OH torsional overtones.⁶² 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⁻¹ to the low frequency D₂₉ 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 low-frequency 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⁻¹.

4 Conclusion

Excited state hessians have been implemented for CC2, with a focus on keeping the scaling of the main memory demands to at most [scr O, script letter O]([scr N, script letter N]²). 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_∞) 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⁻¹ 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₁ excited state than in the planar S₀ 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.

Acknowledgements

Open Access funding provided by the Max Planck Society.

References

1. R. N. Dixon, T. A. A. Oliver and M. N. R. Ashfold, J. Chem. Phys., 2011, 134, 194303.
2. T. N. V. Karsili, A. M. Wenge, B. Marchetti and M. N. R. Ashfold, Phys. Chem. Chem. Phys., 2014, 16, 588–598.
3. A. G. Sage, T. A. A. Oliver, G. A. King, D. Murdock, J. N. Harvey and M. N. R. Ashfold, J. Chem. Phys., 2013, 138, 164318.
4. A. L. Sobolewski, W. Domcke, C. Dedonder-Lardeux and C. Jouvet, Phys. Chem. Chem. Phys., 2002, 4, 1093.
5. T. N. V. Karsili, A. M. Wenge, D. Murdock, S. J. Harris, J. N. Harvey, R. N. Dixon and M. N. R. Ashfold, Chem. Sci., 2013, 4, 2434–2446.
6. T. A. A. Oliver, G. A. King, D. P. Tew, R. N. Dixon and M. N. R. Ashfold, J. Phys. Chem. A, 2012, 116, 12444–12459.
7. M. N. R. Ashfold, G. A. King, D. Murdock, M. G. D. Nix, T. A. A. Oliver and A. G. Sage, Phys. Chem. Chem. Phys., 2010, 12, 1218–1238.
8. G. D. Purvis and R. J. Bartlett, J. Chem. Phys., 1982, 76, 1910–1918.
9. J. F. Stanton and J. Gauss, J. Chem. Phys., 1995, 103, 8931–8943.
10. J. F. Stanton and J. Gauss, Theor. Chim. Acta, 1995, 91, 267–289.
11. J. F. Stanton and J. Gauss, J. Chem. Phys., 1996, 104, 9859–9869.
12. O. Christiansen, J. F. Stanton and J. Gauss, J. Chem. Phys., 1998, 108, 3987–4001.
13. K. W. Sattelmeyer, J. F. Stanton, J. Olsen and J. Gauss, Chem. Phys. Lett., 2001, 347, 499–504.
14. H. Larsen, J. Olsen, P. Jørgensen and O. Christiansen, J. Chem. Phys., 2000, 113, 6677–6686.
15. M. Musiał, K. Kowalska and R. J. Bartlett, THEOCHEM, 2006, 768, 103–109.
16. O. Christiansen, H. Koch and P. Jørgensen, Chem. Phys. Lett., 1995, 243, 409–418.
17. O. Vahtras, J. Almlöf and M. Feyereisen, Chem. Phys. Lett., 1993, 213, 514–518.
18. C. Hättig and F. Weigend, J. Chem. Phys., 2000, 113, 5154–5161.
19. T. S. Balaban, P. Braun, C. Hättig, A. Hellweg, J. Kern, W. Saenger and A. Zouni, Biochim. Biophys. Acta, Bioenerg., 2009, 1787, 1254–1265.
20. M. Head-Gordon, R. J. Rico, M. Oumi and T. J. Lee, Chem. Phys. Lett., 1994, 219, 21–29.
21. M. Head-Gordon, M. Oumi and D. Maurice, Mol. Phys., 1999, 96, 593–602.
22. J. Schirmer, Phys. Rev. A: At., Mol., Opt. Phys., 1982, 26, 2395–2416.
23. A. Köhn and C. Hättig, J. Chem. Phys., 2003, 119, 5021–5036.
24. C. Hättig, Adv. Quantum Chem., 2005, 50, 37–60.
25. D. H. Friese, N. O. C. Winter, P. Balzerowski, R. Schwan and C. Hättig, J. Chem. Phys., 2012, 136, 174106.
26. N. K. Graf, D. H. Friese, N. O. C. Winter and C. Hättig, J. Chem. Phys., 2015, 143, 244108.
27. D. H. Friese, C. Hättig and J. Koßmann, J. Chem. Theory Comput., 2013, 9, 1469–1480.
28. J. Almlöf, Chem. Phys. Lett., 1991, 181, 319–320.
29. F. Egidi, D. B. Williams-Young, A. Baiardi, J. Bloino, G. Scalmani, M. J. Frisch, X. Li and V. Barone, J. Chem. Theory Comput., 2017, 13, 2789–2803.
30. J. Liu and W. Liang, J. Chem. Phys., 2011, 135, 184111.
31. T. U. Helgaker, in Geometrical Derivatives of Energy Surfaces and Molecular Properties, ed. P. Jørgensen and J. Simons, Reidel, Dordrecht, Netherlands, 1986, p. 1.
32. T. Helgaker and P. Jørgensen, Adv. Quantum Chem., 1988, 19, 183–245.
33. TURBOMOLE V7.0 2015, a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989–2007, TURBOMOLE GmbH, since 2007, available from http://www.turbomole.com.
34. B. Fornberg, Math. Comput., 1988, 51, 699–706.
35. I. M. Mills, in Molecular Spectroscopy: Modern Research, ed. K. Narahari Rao and C. Weldon Mathews, Academic Press, Inc., New York and London, 1972, p. 115.
36. W. Mizukami and D. P. Tew, J. Chem. Phys., 2013, 139, 194108.
37. R. D. Amos, N. C. Handy, W. H. Green, D. Jayatilaka, A. Willetts and P. Palmieri, J. Chem. Phys., 1991, 95, 8323–8336.
38. A. M. Gardner and T. G. Wright, J. Chem. Phys., 2011, 135, 114305.
39. A. Andrejeva, A. M. Gardner, W. D. Tuttle and T. G. Wright, J. Mol. Spectrosc., 2016, 321, 28–49.
40. W. D. Tuttle, A. M. Gardner, A. Andrejeva, D. J. Kemp, J. C. Wakefield and T. G. Wright, J. Mol. Spectrosc., 2018, 344, 46–60.
41. R. S. Mulliken, J. Chem. Phys., 1955, 23, 1997–2011.
42. E. B. Wilson, Phys. Rev., 1934, 45, 706–714.
43. T. Den, H.-M. Frey, P. M. Felker and S. Leutwyler, J. Chem. Phys., 2015, 143, 144306.
44. R. Zimmerman and T. Dunn, J. Mol. Spectrosc., 1985, 110, 312–325.
45. T. Cvitaš and J. Hollas, Mol. Phys., 1970, 18, 793–800.
46. A. E. W. Knight and S. H. Kable, J. Chem. Phys., 1988, 89, 7139–7160.
47. H. J. Elston, E. R. Davidson, F. G. Todd and C. S. Parmenter, J. Phys. Chem., 1993, 97, 5506–5518.
48. M. Robey and E. Schlag, Chem. Phys., 1978, 30, 9–17.
49. K. L. Reid, T. A. Field, M. Towrie and P. Matousek, J. Chem. Phys., 1999, 111, 1438–1445.
50. J. Long, C. Qin, Y. Liu, S. Zhang and B. Zhang, Phys. Rev. A: At., Mol., Opt. Phys., 2011, 84, 063409.
51. J. M. Martin, P. R. Taylor and T. J. Lee, Chem. Phys. Lett., 1997, 275, 414–422.
52. J. M. L. Martin, T. J. Lee and P. R. Taylor, J. Chem. Phys., 1998, 108, 676–691.
53. C. G. Hickman, J. R. Gascooke and W. D. Lawrance, J. Chem. Phys., 1996, 104, 4887–4901.
54. C. Minejima, T. Ebata and N. Mikami, Phys. Chem. Chem. Phys., 2002, 4, 1537–1541.
55. A. M. Gardner, A. M. Green, V. M. Tamé-Reyes, V. H. K. Wilton and T. G. Wright, J. Chem. Phys., 2013, 138, 134303.
56. A. M. Gardner, A. M. Green, V. M. Tamé-Reyes, K. L. Reid, J. A. Davies, V. H. K. Parkes and T. G. Wright, J. Chem. Phys., 2014, 140, 114308.
57. D. P. Tew, N. C. Handy and S. Carter, Phys. Chem. Chem. Phys., 2001, 3, 1958–1964.
58. D. P. Tew, N. C. Handy and S. Carter, Mol. Phys., 2004, 102, 2217–2226.
59. D. P. Tew and W. Mizukami, J. Phys. Chem. A, 2016, 120, 9815–9828.
60. T. Dunn, R. Tembreull and D. Lubman, Chem. Phys. Lett., 1985, 121, 453–457.
61. M. Gerhards, W. Perl, S. Schumm, U. Henrichs, C. Jacoby and K. Kleinermanns, J. Chem. Phys., 1996, 104, 9362–9375.
62. T. Bürgi and S. Leutwyler, J. Chem. Phys., 1994, 101, 8418–8429.
63. F. Ramírez and J. L. Navarrete, Vib. Spectrosc., 1993, 4, 321–334.
64. H. W. Wilson, Spectrochim. Acta, Part A, 1974, 30, 2141.
65. S. J. Greaves and W. P. Griffith, Spectrochim. Acta, Part A, 1991, 47, 133–140.
66. T. Y. Koh, S. J. Greaves and W. P. Griffith, Spectrochim. Acta, Part A, 1994, 50, 857–873.
67. H. G. Kjaergaard, D. L. Howard, D. P. Schofield, T. W. Robinson, S.-I. Ishiuchi and M. Fujii, J. Phys. Chem. A, 2002, 106, 258–266.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c8cp06952b

---