Analytic Cubic and Quartic Force Fields Using Density-functional Theory

We present the first analytic implementation of cubic and quartic force constants at the level of Kohn–Sham density-functional theory. The implementation is based on an open-ended formalism for the evaluation of energy derivatives in an atomic-orbital basis. The implementation relies on the availability of open-ended codes for evaluation of one-and two-electron integrals differentiated with respect to nuclear displacements as well as automatic differentiation of the exchange–correlation kernels. We use generalized second-order vibrational perturbation theory to calculate the fundamental frequencies of methane, ethane, benzene, and aniline, comparing B3LYP, BLYP, and Hartree–Fock results. The Hartree–Fock anharmonic corrections agree well with the B3LYP corrections when calculated at the B3LYP geometry and from B3LYP normal coordinates, suggesting that the inclusion of electron correlation is not essential for the reliable calculation of cubic and quartic force constants. © 2014 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution 3.0 Unported License.


I. INTRODUCTION
Vibrational spectroscopy provides a rich and diverse source of information about molecular structure and functionality.For this reason, methods for calculating molecular vibrational spectra were developed already in the early days of quantum-chemical computations, to help elucidate molecular structure and to provide insight into experimental observations.Already in 1958, Bratoz recognized the benefits of calculating the forces acting on the nuclei of a molecule in an analytic manner. 1 However, the breakthrough in terms of computing molecular gradients and force fields came with the efficient implementation of molecular gradients by Pulay in 1969. 2 Molecular gradients are the first-order derivatives of the molecular energy with respect to nuclear displacements and can be determined from the unperturbed electron density and from differentiated one-and two-electron integrals.The complexity increases significantly when going to the molecular force fields or molecular Hessians (corresponding to the second-order derivatives of the molecular energy with respect to nuclear displacements) as also the perturbed electron density matrix is needed in this case.Thomson and Swanstrøm presented the first implementation of molecular Hessians in 1973. 3In the early 1980s, implementations of molecular Hessians for unrestricted and restricted open-shell Hartree-Fock (HF) wave functions were published by Yamaguchi, Schaefer, and co-workers. 4,5 ][21] In parallel with the development of force-field calculations for correlated wave functions, Schaefer, Handy, and co-workers extended the evaluation of geometrical derivatives of the HF energy to third 22,23 and fourth 24 orders.Despite the importance of cubic and quartic force fields for determining, for instance, anharmonic corrections to vibrational frequencies, Fermi resonances, 25 rotation-vibration constants, 26 vibrationally averaged geometries, [27][28][29] and ldoubling constants, 25,30,31 these implementations have seen little use in the literature.Part of the reason for this may be that anharmonic corrections are important only in highaccuracy studies, [32][33][34] where the HF approximation may not be sufficient.Also, the implementation of efficient schemes for obtaining anharmonic force fields by numerical differentiation of forces and force constants at highly correlated levels of theory has proven feasible, even for relatively large molecules. 35,36  some extent, the low cost of Kohn-Sham DFT has changed this picture.In a series of studies, Barone and coworkers have shown that high-quality harmonic force fields in combination with DFT anharmonic corrections provide reliable estimates of anharmonic force fields 37,38 and of anharmonic corrections to intensities, 39 thus demonstrating that anharmonic effects can now be studied straightforwardly for large and complex molecular systems.In their studies, geometrical derivatives beyond second order were determined by finite differences. 32,40,41 Te high cost of such finitedifference schemes and their numerical instability provide a strong motivation for developing analytic methods for thirdand fourth-order geometrical derivatives at the DFT level.Also, fifth-and sixth-order geometrical derivatives are needed in fourth-order vibrational perturbation theory (VPT4); recently, we have calculated analytically quintic and sextic force fields at the HF level of theory. 42n this work, we describe the analytic calculation of cubic and quartic force constants at the level of Kohn-Sham DFT, using generalized-gradient-approximation (GGA) and hybrid functionals.This work builds on several developments in our groups in recent years.In particular, we use the framework of an atomic-orbital (AO) based, open-ended quasienergy response-theory formalism described by Thorvaldsen et al., 43 which for force fields reduces to regular energyderivative theory in the AO basis, 44 extended to fourth-order derivatives.The geometrical derivatives of the one-electron integrals arising from the geometry dependence of the AOs are evaluated using the one-electron integral framework of Gao, Thorvaldsen, and Ruud. 45The evaluation of geometrical derivatives of the two-electron repulsion integrals follows the approach of Reine, Tellgren, and Helgaker, 46 expanding solid-harmonic Gaussians directly in Hermite Gaussians.We furthermore extend automatic differentiation of exchangecorrelation kernels 47 to include corrections arising from the dependence of the AOs on the nuclear positions.Finally, we demonstrate the usefulness of the code by evaluating the cubic and quartic force constants and anharmonic force-field corrections for selected molecules.
The bulk of this paper is organized as follows.In Sec.II, we give a brief account of the AO-based energyderivative framework used by us and give the expressions for the cubic and quartic force constants.A brief description of the evaluation of the exchange-correlation contribution is also given.Section III contains computational details of the calculations, while Sec.IV presents and discusses the results.Finally, in Sec.V, we give some concluding remarks and perspectives for the analytic calculation of higher-order properties that involve geometrical distortions.

II. THEORY
We here present the theory behind our AO-based implementation of DFT cubic and quartic force fields, the AOformulation ensuring that the approach is suitable for linearscaling methodology. 48The approach builds on the general AO-based framework for time-and perturbation-dependent basis sets by Thorvaldsen et al., 43 here applied to timeindependent perturbations.We note that, even though explicit equations are given for the evaluation of the cubic and quartic force fields, our implementation uses a recursive scheme, for which explicit expressions for energy derivatives are not needed. 49Compared with the earlier implementations by Schaefer, Handy, and co-workers, [22][23][24] our implementation is extended to DFT and formulated entirely in the AO basis, although at present a molecular-orbital-based response solver is used to determine the perturbed density matrices.
In Sec.II A, we review the basic theory of the analytic calculation of geometrical derivatives at the DFT level, providing explicit expressions for the cubic and quartic force constants; next, in Sec.II B, we describe the evaluation of exchange-correlation contributions to the cubic and quartic force constants, combining the perturbation dependence of the overlap distributions with the use of automatic differentiation to evaluate the higher-order exchange-correlation kernel derivatives. 47

A. AO-based energy derivative theory
We follow the notation of the AO-based response theory for self-consistent-field (SCF) methods with time-and perturbation-dependent basis sets by Thorvaldsen et al., 43 specialized to static perturbations.In Kohn-Sham DFT, the energy can in the AO basis be written as In this expression, "Tr" indicates that the trace of the matrix products on the right-hand side of the equation is taken, E xc [ρ (D)] is the exchange-correlation energy, which is a functional of the generalized density vector ρ (see Sec. II B), h nuc is the nuclear repulsion energy, and D is the AO density matrix.The density matrix fulfills the idempotency relation We have in Eqs. ( 1) and ( 2) also introduced the one-electron integral matrix, h, the two-electron integral matrix constructed from D with γ fractional exchange, G γ (D), and the overlap matrix, S, whose elements are given by The χ μ are spherical-harmonic Gaussian AOs and the summation in Eq. ( 3) is over atomic nuclei at R K and with charge Z K .The two-electron integrals are defined in the conventional manner as with integration over all spin and spatial coordinates.
The optimized Kohn-Sham density fulfills the SCF condition where the Kohn-Sham matrix F is defined as with Following Refs.43, 49, and 50, we take as our starting point for the generation of higher-order derivatives the energy-gradient Lagrangian L a defined as (10)   where the superscript a indicates differentiation with respect to an applied perturbation of strength ε a .In Eq. (10) we have introduced the energy-weighted density matrix as well as matrices that represent the constraints on the unperturbed reference state-in particular, the idempotency and SCF-state matrices, respectively, For an optimized SCF state, Eqs. ( 2) and ( 7) can be written compactly as Z = 0 and Y = 0. Additionally, we have introduced the Lagrange multipliers λ a and ζ a , respectively, for these constraints.In Ref. 43, it is shown that Eq. ( 10) is variational in D if the zeroth-order multipliers are defined as The subscript a on the multipliers does not indicate differentiation, merely the relation to L a .Since we require Eq. ( 10) to be variational in the density D and the multipliers λ a and ζ a , we can take advantage of the 2n + 1 rule for the density and the 2n + 2 rule for the multipliers 50,51 when differentiating the energy gradient Lagrangian.

Molecular gradient
Let us first consider the first-order geometrical derivatives of the molecular energy.In this case, Eq. ( 10) simplifies to Pulay's expression from 1969 where no derivatives of the Lagrange multipliers are required because of the 2n + 2 rule. 50We note that the molecular gradient can be determined from a knowledge of the unperturbed density alone, in accordance with the 2n + 1 rule.

Molecular Hessian
Differentiating Eq. ( 10) with respect to b, keeping only terms that fulfill the 2n + 1 and 2n + 2 rules, we obtain for the molecular Hessian the expression As for the gradient, no zeroth-order multipliers are required.However, the first-order perturbed density matrix is needed to calculate the molecular Hessian; we return to the evaluation of the perturbed densities in Sec.II A 5. Henceforth, we adopt the subscript notation of Ref. 43, writing the Hessian as where superscripts denote total derivatives.The subscripts (k, n) specify the maximum order of the perturbed density matrix D: to order k for collections of perturbations involving perturbation a, and to order n for collections of perturbations not involving perturbation a.The notation n W for the SW term specifies in a similar manner the maximum order of differentiation of W in this term, as dictated by the value of n.
The Hessian expression in Eq. ( 17) is not explicitly symmetric in a and b (the numerical values, of course, are).As shown by Sellers, an explicitly symmetric formula can be advantageous from a numerical point of view. 52

Cubic force constants
For the calculation of cubic force constants, the thirdorder energy derivative is needed.Proceeding as above, we differentiate the gradient Lagrangian twice, keeping only terms that fulfill the 2n + 1 and 2n + 2 rules where the subscript 1 on the right-hand side indicates that only terms with density matrices up to first order are to be included.For example, in the expressions there are no terms containing D bc .In Ref. 43 it is shown how to rewrite the above expression in a more symmetric form Here, a prime on the subscript means that the subscript refers to the maximum order of differentiation of S, D, and F (rather than the order of D), for example, Even though the two first terms E abc 1,1 and (SW) abc 1 W in Eq. ( 22) are explicitly symmetric in abc, the remaining terms are not.It is possible to symmetrize the expression but this does not give any computational benefits.

Quartic force constants
Following the same procedure as for the cubic force constants, it is possible to derive the following expression for the fourth-order energy derivative: 43 As expected, we need up to second-order perturbed density matrices and first-order multipliers.By retaining the secondorder density matrices involving the perturbations bcd, it is possible to eliminate the first-order multipliers 43 (26) In this notation, the latter expression appears more compact, but if one expands the different terms, one finds that Eqs. ( 25) and ( 26) are of similar complexity.Finally, we note that neither expression is explicitly symmetric in the perturbation labels.

Perturbed density matrices
To evaluate the expressions for the energy derivatives, we need first-and second-order perturbed density matrices.Since the unperturbed density matrix satisfies the idempotency condition and the SCF equations, the perturbed density can be obtained by differentiating Eqs. ( 2) and (7).In other words, the matrix D b is the solution to the simultaneous equations We rewrite Eq. ( 27) by collecting only terms containing D b on one side, yielding which has a solution of the general form where D b has been partitioned into a particular part D p and a homogeneous part D h .In the homogeneous equation, N = 0 and the equation is automatically satisfied by the ansatz To determine X, we use the differentiated SCF equation.Inserting Eq. ( 30) into Eq.( 28) and collecting all terms containing X on the left, we arrive at the coupled-perturbed Kohn-Sham equations (or more generally the linear response equations) E [2] X = M, (33)   where we can identify the electronic Hessian E [2] in the AO basis and the right-hand side In the same way, the second-order perturbed density matrix D bc can be determined from the equations Z bc = 0 and Y bc = 0.The resulting equations have the same structure as in Eqs. ( 30) and ( 33), the matrices N and M now being To summarize, we must solve one set of linear response equations for each perturbed density matrix, where the right-hand side depends on (perturbed) density matrices of lower orders.We refer to Ref. 43 for further details.

B. Evaluation of exchange-correlation contributions
We employ an exchange-correlation energy E xc defined as the integral over a local function xc (r) that depends on the density n(r) and its Cartesian gradient ∇n(r): where ∇n(r) To simplify notation, we henceforth collect the density variables n(r) and ∇n(r) in a generalized density vector ρ(r).This notation also simplifies a generalization of our implementation to other density variables such as the kinetic-energy density in meta-GGA functionals and density variables in spin DFT and current DFT.The exchange-correlation energy and the exchangecorrelation potential matrix are integrated on a numerical grid defined by a set of suitably chosen grid points r i and grid weights w i , according to When differentiating the exchange-correlation energy and potential matrix we ignore the contribution from the gridweight derivatives.The importance of grid-weight derivatives in the evaluation of geometrical derivatives at the DFT level has been discussed by Baker et al. 53 and Johnson and Frisch. 17The extension to higher order is not straightforward, and for this reason we use very large grids in order to minimize the errors arising from the lack of grid-weight derivative contributions, and the quality of the results has been verified by test calculations against numerically calculated derivatives.
For the implementation of DFT analytic cubic and quartic force constants, we need up to fourth-order geometrical derivatives of the exchange-correlation energy density ( a xc , ab xc , abc xc , and abcd xc ) and up to second-order geometric derivatives of the exchange-correlation potential matrix contributions (v a xc and v ab xc ).The exchange-correlation energy density derivatives are evaluated using the following expressions: where the arguments of densities, functional derivatives, and overlap distributions have been omitted for notational clarity.
In our code, the contractions of the functional derivative vectors with the perturbed generalized density vectors are not explicitly programmed.Instead, we obtain the perturbed exchange-correlation energy densities xc directly from the XCFUN program 47,54 by forming a generalized density Taylor series expansion (ρ, ρ a , ρ b , ρ ab , . . .), which is internally contracted with the density functional Taylor expansion.This approach significantly reduces the complexity of the exchange-correlation integrator.
If the total energy had been variational with respect to the density ρ, then, according to the 2n + 1 rule, we would only need the first-order (second-order) perturbed densities for the cubic (quartic) force field.In our case, the energy is variational with respect to the AO density matrix D, which means that we still need the second-and third-order (thirdand fourth-order) perturbed densities.
This dependence is given as Depending on the perturbation order, many of the above terms are omitted when applying the 2n + 1 rule to the density matrix.Note that the omitted terms are not zero by themselves, but only in combination with non-exchangecorrelation terms containing the same density matrices.Here we have used into the generalized overlap distribution vector ( p,q ρ ) μν .Having discussed exchange-correlation energy density contributions, we now turn to the exchange-correlation potential matrix contributions.The perturbations can either act on the generalized overlap distribution or on the functional derivative term, giving where we have used Finally, we note that an efficient implementation of the density evaluation and matrix distribution routines is essential, bearing in mind the large number of terms that need to be evaluated.We evaluate both the densities and the matrix elements in a blocked manner, allowing mathematical matrixmatrix multiplication libraries to be used in conjunction with efficient pre-screening techniques.

III. COMPUTATIONAL DETAILS
To calculate the cubic and quartic force constants, the recursive implementation 49 of the open-ended response-theory framework by Thorvaldsen et al. 43 has been used, as provided by the OPENRSP program package.We use the DAL-TON program package 55 as a backend for the calculation of undifferentiated integrals and the unperturbed and perturbed density matrices, which are obtained with the linear response solver of Jørgensen et al. 56 The calculation of properties associated with one-electron integrals was carried out using the GEN1INT library, 57 building on the flexible integral evaluation scheme of Gao and co-workers. 45The differentiated twoelectron integrals were mainly calculated using Thorvaldsen's CGTO-DIFF-ERI code, 58 which uses the scheme of Reine et al. for the evaluation of differentiated two-electron integrals using solid-harmonic Gaussians, 46 but some of the lowerorder contributions were calculated using DALTON.The differentiated exchange-correlation energy and potential contributions needed for the cubic and quartic force constants were calculated using the XCFUN library, 54 which uses automatic differentiation for evaluating the derivatives of the exchangecorrelation energy. 47We have used an in-house integrator to perform the integration of the exchange-correlation contributions.
Cubic and quartic force constants in the Cartesian basis were calculated at the HF and DFT levels of theory for methane, ethane, benzene, and aniline.For methane, we performed a basis-set convergence study using the 6-31G 59 and the correlation-consistent basis sets 60 of double-, triple-, and quadruple-zeta quality (cc-pVDZ, cc-pVTZ, and cc-pVQZ).For the other molecules, we have used the cc-pVTZ basis set for ethane and the cc-pVDZ basis set for benzene and aniline.In order to explore the sensitivity of the results to the choice of exchange-correlation functional, both the BLYP [61][62][63][64] and the B3LYP 65 functionals have been used.
For the HF and B3LYP calculations, the geometry was optimized and the molecular Hessian was calculated at the DFT (B3LYP) level of theory with the DALTON program package, 55 using the same basis as in the anharmonic force field calculations.The B3LYP Hessian was used in the vibrational analysis-both for evaluating the harmonic vibrational frequencies and for transforming the anharmonic force constants to a reduced normal coordinate basis 26 before evalu-ating the fundamental frequencies (vide infra).Although not consistent, this approach circumvents the well-known deficiencies of the HF method for harmonic frequencies and allows us to get a better impression of the quality of the HF cubic and quartic force constants.For the calculations involving the BLYP functional, the geometry optimization, the vibrational analysis, and the cubic and quartic force constants were calculated using this functional, allowing us to compare directly the results obtained using the BLYP and B3LYP functionals.
In the calculations, we have converged the coupledperturbed Kohn-Sham equations to a relative norm of 10 −6 , observing no problems with convergence of the response equations.To reduce the errors arising from the lack of gridweight derivative contributions, we have used an ultrafine grid with a radial quadrature accuracy of 2 × 10 −15 and with an angular expansion order of 64.
From the cubic and quartic force constants, anharmonic frequency corrections were calculated using the generalized vibrational second-order perturbation (GVPT2) model, 66,67 in which terms that are too large because of Fermi resonances are excluded from the perturbational treatment 68 and treated variationally. 41The threshold criteria for the identification of Fermi resonances are those used by Bloino and Barone 69 except for ethane, where the threshold for the Martin parameters was increased to 1.5 cm −1 from the default value of 1 cm −1 , to avoid a splitting of degenerate modes due to unevenly distributed interactions between these modes and a different set of two degenerate modes, which would otherwise lead to an unsymmetric identification of Fermi resonances.We refer to Refs.41 and 69 for more details about the GVPT2 model and the treatment of Fermi resonances.All rotational effects, as described by the rotational constants and the Coriolis coupling constants in the GVPT2 scheme, are disregarded in the present work.

IV. RESULTS AND DISCUSSION
In Tables I-IV, we have listed the calculated fundamental frequencies of methane, ethane, benzene, and aniline, respectively.Regarding the basis-set dependence of the anharmonic corrections, the methane results in Table I indicate it is rather weak, with small differences between 6-31G and cc-pVQZ anharmonic corrections, the largest difference between the HF/6-31G and HF/cc-pVQZ results being 6 cm −1 (4%).
Regarding the differences between the various levels of electronic-structure theory, we see from Table I that the difference between the HF and the B3LYP corrections (both calculated using B3LYP geometries, harmonic frequencies, and normal coordinates) are not large for methane, the absolute value of the B3LYP corrections being on average smaller than 10%.From the results for the larger molecules in Tables II-IV, this behaviour appears to be a general trend, but with some discrepancies being slightly larger than 10%.Also, in a few cases, the B3LYP anharmonic corrections are larger than the corresponding HF corrections-for modes 8 and 12 in benzene and for mode 13 in aniline, for instance, the B3LYP correction is substantially more negative than the HF correction.These differences are mostly the result of differences between the HF and B3LYP values for the associated diagonal (iiii) quartic force constants; for mode 12 in benzene, differences in the semidiagonal (iijj) quartic force constants are also important.Overall, the BLYP anharmonic corrections are in good agreement with the B3LYP corrections, but with some exceptions.For mode 2 in benzene, for example, the anharmonic correction is substantially less negative with the BLYP exchange-correlation functional, because of different identifications of Fermi resonances at different levels of theory.In general, however, the BLYP anharmonic corrections are slightly more negative than the B3LYP corrections.
Among the different levels of theory applied here, the B3LYP results are in best agreement with the experimental fundamental frequencies.The listed HF results are of comparable quality but have been obtained at the B3LYP geometry and are based on the B3LYP harmonic vibrational analysis; they would have been considerably worse had they been based on HF quantities alone.In any case, the calculation of HF anharmonic corrections based on DFT geometries and harmonic frequencies should be a viable approach in many cases.Likewise, we expect harmonic frequencies calculated at higher levels of theory-for instance, at the coupled-cluster level of theory-to perform well in combination with SCF anharmonic corrections; indeed, such an approach has been used in earlier works. 38,70,71 T BLYP results consistently show the poorest agreement with experiment.However, much of the discrepancy arises from inaccurate harmonic frequencies.For all systems, the BLYP corrections are mostly close to the HF and B3LYP corrections-however, for methane, the BLYP correspondence to experiment for the harmonic frequencies is clearly inferior to the B3LYP correspondence.This is further accentuated when we note that, for the high-frequency modes in methane, the derived experimental anharmonic corrections are in general smaller than the calculated ones, suggesting that the experimental harmonic frequencies are underestimated.In general, the agreement between the B3LYP and BLYP anharmonic corrections is good, supporting the notion that it is the harmonic frequencies that are poorly described by the BLYP functional.
In order to get a better global understanding of the performance of the different computational levels, we have in Tables II-IV also collected the mean absolute errors (MAEs) for the different computational levels compared to experiment.We note that the harmonic B3LYP frequencies in general are about 2.5% off the experimental frequencies, in line with the recommended scaling factors often used for B3LYP calculations of 0.9679. 72We note that for aniline, the MAE is somewhat larger than for the other molecules, but this is largely due to mode 30 which displays a MAE of almost 25%.Excluding this mode, the MAE for the B3LYP harmonic frequencies is 2.4%.Including anharmonic corrections to the B3LYP harmonic frequencies, either at the HF or B3LYP levels of theory, brings the MAE down to about 1% with the HF error being slightly larger and the B3LYP error slightly smaller than 1%, the differences in general being small.However, once again mode 30 in aniline is an interesting case, clearly showing the superiority of the B3LYP method over the HF method for difficult cases, as the HF fundamental frequency for this mode is off by 23% from the experimental value, B3LYP instead being only 0.8% off the experimental data.
The BLYP model in contrast provides a much better MAE than the B3LYP model for harmonic frequencies.However, this agreement is fortuitous and when adding the anharmonic corrections to the BLYP data, the MAE actually increases, from about 1% to 3%-4%.If we instead add the BLYP corrections to the B3LYP harmonic frequencies, the MAE becomes comparable to that obtained with HF (about 1%).
The differences between the anharmonic corrections obtained at various levels of theory are relatively small.However, for certain spectroscopic processes that have recently received increased attention-for example, the doubly vi-brationally enhanced four-wave-mixing using two incident infrared lasers discussed in Ref. 73-the principal twodimensional spectroscopic features may consist of closely spaced peaks separated by a distance related to the anharmonic coupling between the modes involved, in addition to lower-order contributions.The shape of such features can be very sensitive to the values of the anharmonic corrections, putting higher demands on the accuracy in the calculated anharmonic corrections for two-dimensional than using the HF model: the HF anharmonic corrections are in good agreement with the B3LYP corrections when the HF calculations are based on the B3LYP structure and harmonic force field.
In future work, we intend to apply the calculated cubic and quartic force constants to obtain fundamental vibrational frequencies for use in various spectroscopic designs and to obtain anharmonic corrections to spectroscopic intensitiesin particular with an eye towards multidimensional vibrational spectroscopies, where anharmonic effects are important.

TABLE I .
Harmonic fundamental vibrational frequencies ω, corrected fundamental frequencies ν, and anharmonic vibrational corrections δ for methane.All values are in cm −1 .
a Experimental data taken from Ref.74and ordered by decreasing frequency.

TABLE II .
Harmonic fundamental vibrational frequencies ω, corrected fundamental frequencies ν, and anharmonic vibrational corrections δ for ethane using the cc-pVTZ basis set.All values are in cm −1 .Experimental data taken from Ref.75and ordered by decreasing frequency. a

TABLE III .
Harmonic fundamental vibrational frequencies ω, corrected fundamental frequencies ν, and anharmonic vibrational corrections δ for benzene using the cc-pVDZ basis set.All values are in cm −1 .
a Modes 7 and 8 and modes 13 and 14 were switched after an analysis of the normal coordinates to agree with the B3LYP ordering.bExperimentaldata taken from Ref.76and ordered by decreasing frequency.
This work has received support from the Research Council of Norway through a Centre of Excellence Grant (Grant No. 179568/V30) and a research grant (Grant No. 191251), the European Research Council through a ERC Starting Grant (Grant No. 279619), and from the Norwegian Supercomputer Program through a grant of computer time.T.H. and U.E. acknowledge support from the European Research Council under the European Union Seventh Framework Program through the Advanced Grant ABACUS, ERC Grant Agreement No. 267683.The work was also supported by COST (Action CoDECS: COnvergent Distributed Environment for Computational Spectroscopy).We are grateful to Julien Bloino and Vincenzo Barone for helpful discussions.