Analytic Calculations of Anharmonic Infrared and Raman Vibrational Spectra

a Using a recently developed recursive scheme for the calculation of high-order geometric derivatives of frequency-dependent molecular properties [Ringholm et al.,, we present the first analytic calculations of anharmonic infrared (IR) and Raman spectra including anharmonicity both in the vibrational frequencies and in the IR and Raman intensities. In the case of anharmonic corrections to the Raman intensities, this involves the calculation of fifth-order energy derivatives—that is, the third-order geometric derivatives of the frequency-dependent polarizability. The approach is applicable to both Hartree–Fock and Kohn–Sham density functional theory. Using generalized vibrational perturbation theory to second order, we have calculated the anharmonic infrared and Raman spectra of the non-and partially deuterated isotopomers of nitromethane, where the inclusion of anharmonic effects introduces combination and overtone bands that are observed in the experimental spectra. For the major features of the spectra, the inclusion of anharmonicities in the calculation of the vibrational frequencies is more important than anharmonic effects in the calculated infrared and Raman intensities. Using methanimine as a trial system, we demonstrate that the analytic approach avoids errors in the calculated spectra that may arise if numerical differentiation schemes are used.


Introduction
The calculation of spectroscopic phenomena involving molecular vibrations is an example of a successful application of theoretical chemistry to aid the interpretation of experimental observations. The calculation of molecular structure and vibrational spectra was made possible by the pioneering work of Bratoz, Pulay, Pople, Schaefer and others in developing methods for calculating analytical geometric derivatives of the molecular energy. [1][2][3][4][5] Second-order geometric derivatives of the energy have since been derived and implemented for a wide range of correlated wave functions, 3,[5][6][7][8][9][10][11][12][13] as well as at the level of density functional theory (DFT). [14][15][16] For a detailed historical account we refer to recent reviews of molecular properties in general and molecular force fields in particular. [17][18][19] A commonly used approximation in the study of vibrational spectroscopies is the double-harmonic approximation, 20 where molecular vibrations are described as harmonic oscillators and where the fundamental properties describing the spectroscopic intensities are determined by the first-order geometric derivative of the polarization property governing the spectroscopic phenomenon under study. This corresponds to a description which obeys the well-known selection rules for e.g. infrared (IR) and Raman spectroscopies, 21 where it is the magnitude of the first-order geometric derivative of the molecular electric dipole moment and polarizability, respectively, that determines the spectroscopic intensity associated with the excitation to a singly excited state of a particular normal mode. Although coupledcluster theory can provide vibrational frequencies of high accuracy, [22][23][24][25][26][27][28] its computational cost prevents its routine use for larger molecules. For this reason, density-functional theory (DFT) has been gaining increasing popularity in recent years and has been used for calculations of Raman spectra of molecules as large as buckminsterfullerene. 29 The calculations have often been done in combination with scaling of the frequencies in order to account for anharmonicities and errors inherent in the exchange-correlation functional used. 30,31 A typical computational protocol has been to determine the frequencies of the (harmonic) normal modes by DFT using the B3LYP functional coupled with intensity calculations performed at the HF or DFT level, [32][33][34] the choice of level of theory for the intensity calculations depending on the computational tools available for the calculation of the necessary geometric derivative of the pertinent polarization property.
By leaving the double harmonic approximation, it is possible to obtain a more accurate description of both the vibrational frequencies and the spectroscopic intensities. 23,[35][36][37][38] For the latter, the introduction of anharmonic effects will enable both an improved description of the intensities associated with a single excitation of a particular normal mode as well as introduce the leading-order contributions to intensities associated with transitions corresponding to the simultaneous excitation of two or more vibrational quanta, either involving only a single normal mode or several of them, often referred to as overtone and combination bands, respectively.
Calculations of anharmonic contributions for the purpose of correcting vibrational frequencies have regularly been carried out, 17,23,[39][40][41] requiring at least third-order geometric derivatives of the molecular energy, from here on referred to as the cubic force constants. We will refer to the corresponding fourth-order derivatives as the quartic force constants. Calculations of the cubic and quartic force constants have previously almost without exception been done using numerical differentiation. 23,[41][42][43] The only exception is the analytic calculation of cubic and quartic force constants at the HF level reported by Handy and coworkers. 44,45 Recently, we presented an analytic implementation of cubic and quartic force constants at the DFT level 46 by the use of a newly developed recursive code 47 for the calculation of molecular properties by response theory. 48 For the IR and (regular) Raman spectroscopies, programs that allow for the analytic calculation of the required first-order geometric derivatives of the dipole moment and polarizability, respectively, have been available for some time. [49][50][51] The calculation of anharmonic corrections to the intensities in these spectroscopies requires both the development of the necessary vibrational perturbation theory 38 to obtain expressions for these corrections and the possibility of calculating second-and third-order geometric polarization property derivatives, as well as the cubic and quartic force constants, that enter into these expressions. Programs that would allow for the analytic calculation of some of these properties are available, but such calculations have mainly been restricted to the HF level of theory, and for some of the properties (and more so if a DFT description is desired), the researcher has had to resort to numerical differentiation. Analytic calculation offers several advantages over numerical methods such as higher attainable accuracy and ease of computation, 51 as numerical derivatives are sensitive to the finite perturbation/geometry displacements employed, and this can have significant effects on the results if not managed carefully. [52][53][54] For these reasons, analytic methods are preferred over numerical ones.
In this work, we present the first application of our recursive approach for the analytic calculation of the anharmonic vibrational frequencies and infrared and Raman intensities of methanimine as well as nitromethane and its mono-and di-deuterated isotopomers. Methanimine has been shown to be very sensitive to the numerical differentiation parameters 52 and thus provides a good illustration of the advantages of the analytic approach. The nitromethane isotopomers have been selected because experimental spectra display a large number of combination and overtone bands, for which calculation calls for the use of an anharmonic treatment. We remark that anharmonic effects have also been found to contribute appreciably to the spectroscopic intensities for several other molecules. 38 The rest of the paper is organized as follows: in Section 2, we outline the theoretical foundation for the analytic calculation of anharmonic corrections to vibrational frequencies and IR and Raman intensities. In Section 3, we provide details about the computational setup used for the calculations on our chosen systems. We present and discuss the results of our calculations in Section 4, and make some concluding remarks in Section 5.

Theory
We will begin in Section 2.1 by outlining how the high-order molecular properties used in this work can be calculated analytically through the use of our recently developed recursive response code and then in Section 2.2 proceed to show how these properties can be used to determine anharmonic corrections to vibrational frequencies and IR and Raman intensities. Although the general framework has been described previously, [46][47][48] this work is the first report of fifth-order analytic derivatives involving geometrical distortions.

Analytic calculation of response properties
A detailed presentation of the response theory, which in our approach is fundamental for the analytic calculation of the cubic and quartic force fields and the high-order geometric derivatives of the dipole moment and polarizability that are needed in this work, is too long to show here, and we will therefore restrict ourselves to the most salient features. We refer to the original work 48 for a more thorough treatment, and to our recent work 47 for a description of the recursive implementation used in the present work.
Our analytic scheme uses as a starting point that linear response functions described by perturbations a and b can be formulated as perturbation strength e i (i = a, b,. . .) derivatives of a quasienergy Lagrangian gradient, expressed in a densitymatrix (D) formulation as 48 where the derivative is evaluated at zero perturbation strength and where higher-order response functions can be found by further differentiation of eqn (1). A tilde is used to represent a quantity considered at an arbitrary perturbation strength, and the absence of a tilde denotes evaluation at zero perturbation strength. The quasienergy Lagrangian L a is given bỹ where we have introduced the atomic orbital (AO) overlap matrix S as wherew is an atomic orbital, and where the energy-and frequency-weighted Fock matrix W is defined as where the generalized Kohn-Sham Fock matrixF is given bỹ We also introduced the generalized energyẼ as In eqn (5)- (7), we introduced the half-time-differentiated overlap matrix T, the one-electron matrix h, the external field operator Ṽ t and the two-electron matrix G g with g-fractional exchange asT and also the exchange-correlation contributions F xc and Ẽ xc [r(D)] in addition to a nuclear potential operator h nuc . Here and throughout the paper, atomic units are used unless otherwise stated. Molecular properties characterized by a perturbation tuple abc. . . can therefore be formulated as derivatives of the quasienergy Lagrangian gradient as where we have introduced a short-hand notation for differentiation and tracing by and respectively. This theory is sufficient to define any response function using the so-called n + 1 rule formulation, 55 where the calculation of a response property of order n + 1 requires the calculation of the density matrix perturbed to order n. However, other formulations placing other conditions on which perturbed density (and Fock) matrices must be calculated are possible. 55 Let us represent the idempotency of the density matrix and the time-dependent self-consistent field (TDSCF) conditions as the matrices Y and Z, respectively, so that and where the notation and has been introduced, and where adjungation is defined to happen before time differentiation. It can be shown that the ansatzl for the multiplierl a for Y leads to the definition of the multiplier z a for Z as It is then possible to make a general expression for the quasienergy Lagrangian for the calculation of response properties as where the values of k and n in the various terms denote, with minor variations, to which orders perturbed Fock and density matrices must be calculated in order to evaluate this expression: the value of k determines to which order the perturbed matrices must be calculated for perturbation tuples involving perturbation a, whereas the value of n determines the same for perturbation tuples not involving perturbation a. We have that k + n = N À 1, where N is the order of the property considered, and k must be chosen as an integer in the interval k A [0,(N À 1)/2], where (N À 1)/2 is rounded down for even N. In this work, we do not discuss how the necessary perturbed Fock and density matrices can be calculated, as it is described in detail in ref. 48. We remark, however, that since the calculation of high-order properties requires solving linear response equation systems and since this part of the calculation is computationally expensive, a judicious (k,n) rule choice may give a significant reduction in the number of such systems to be solved, both compared to other rule choices and to numerical differentiation schemes. For instance, for the calculation of cubic force constants, the choice (k,n) = (1,1) makes it necessary to solve M systems, where M is the number of geometrical coordinates, whereas (k,n) = (0,2) results in M 2 such systems. Similarly, a scheme where an analytically calculated molecular Hessian is differentiated numerically by nuclear displacements results in the number of such systems being of the order of M 2 . Similar savings can be achieved for other properties.
With the recursive program developed by our group, 47 it is possible to evaluate eqn (23) for any response property, including the calculation of the required perturbed Fock and density matrices, as long as external routines are available that can provide the necessary (un)perturbed one-and two-electron integral contributions, 56-58 exchange-correlation contributions 59,60 and perturbed nuclear potential contributions, and solve the response equations 61,62 that arise during the evaluation of perturbed Fock and density matrices. More information about the external modules used in this work is given in Section 3. All such modules used in the present work have been parallelized; see e.g. ref. 63.

Anharmonic corrections to vibrational frequencies and spectroscopic intensities
Having determined the harmonic vibrational frequencies and normal modes of vibration from the well-established eigenanalysis of the molecular Hessian, 20 it is possible to make anharmonic corrections to fundamental vibrational frequencies and frequencies corresponding to combination or overtone excitations of the normal modes by a second-order perturbational approach, where the resulting expressions involve the cubic and quartic force constants and Coriolis vibration-rotation coupling constants. In the VPT2 approach, 23,35,36 the corrected fundamental vibrational frequencies n i 1, first overtone frequencies n i 2 and first combination frequencies n i 1 j 1 are given as, respectively where O ijk is defined as In the above expressions, o i denotes a harmonic fundamental frequency, f ijk and f ijkl are cubic and quartic force constants, respectively, B a is the rotational constant for axis a, and z a ij is a Coriolis coupling constant.
The method chosen in the present work is the so-called generalized vibrational second-order perturbation (GVPT2) model. 38,41 In this method, a VPT2 treatment of the molecular vibrations is used, except for the cases where Fermi resonances are considered to have occurred. In these cases, the terms in the VPT2 treatment that are affected by the Fermi resonance are not included, 23 and the affected frequencies are instead resolved in a variational approach.
Expressions for corrections to spectroscopic intensities can also be identified by a perturbation-theory approach. In a recent work by Bloino and Barone, 38 GVPT2 expressions for IR and Raman intensities have been derived. The expression for the IR intensities is (30) and in classical Raman spectroscopic measurements, the unpolarized (as well as polarized) scattering intensity at a temperature T, related to the Raman cross section, is given by where and where n i = o i in the harmonic approximation and is given by eqn (24)- (26) in the anharmonic GVPT2 treatment, n 0 is the frequency of the incident laser in the Raman experiment, and h i 0i represents the transition moment of the relevant polarization property from the vibrational ground state to the ith vibrational excited state. In the double-harmonic treatment, these transition moments are determined by first-order geometric derivatives of the polarization property (hPi 0i 1 = qP/qq i ), whereas the anharmonic expressions also involve the secondand third-order geometric derivatives of the polarization property and the cubic and quartic force constants. The resulting expressions in the anharmonic case are large and we refer to the work of Bloino and Barone 38 where the complete expressions are reported. Altogether, the expressions used in the complete VPT2 treatment involve the first-, second-, and third-order geometric derivatives of the molecular electric dipole moment and polarizability in the IR and Raman case, respectively, in addition to the cubic and quartic force constants, meaning that the highest-order property that must be calculated, i.e. the cubic force constants of the frequency-dependent polarizability, is a -order energy derivative. The contributions to this property can be identified from eqn (23) and are shown here in order to demonstrate the complexity involved in the analytic calculations performed in this work and to justify the use of a recursive approach.
The third-order geometric derivative of the polarizability can be defined from a perturbation tuple (a, b, c, d, e), where perturbations a, b and c correspond to differentiation with respect to geometrical displacements, and d and e to differentiation with respect to a frequency-dependent electric dipole perturbation. Denoting a geometric perturbation as g and the two electric dipole perturbations as f o and f Ào , where, respectively, each perturbation is associated with a positive or negative frequency o, eqn (23) takes the form where the rule choice (k,n) = (2,2) is used because this will give the lowest computational cost. Omitting terms that must be zero straightforwardly or because, in the differentiation carried out, there was lack of dependence on the perturbation operators and, for the sake of brevity, writing contributions that are permutations of identical operators only once, the terms in eqn (34) can be written as and where, for example, W gfofÀo 2 0 from eqn (37) is and where the other differentiated W, Y, and Z terms are of a similar complexity. We consider the length of these expressions, in particular eqn (40), and the corresponding complexity in treating them, as strongly supporting the use of a recursive approach for calculations of the high-order properties required for the GVPT2 treatment, and in a similar manner, automated approaches based on automatic differentiation are needed in order to evaluate the differentiated exchange-correlation energy and kernel E gggfofÀo xc . 59

Computational details
To compute the cubic and quartic force constants and the first-, second-and third-order geometric derivative tensors of the electric dipole moment and of the electric dipole polarizability, the recursive implementation 47 of the open-ended response theory framework of Thorvaldsen et al. 48 has been used. This formalism has been implemented in a development version of the Dalton2013 program package. 64,65 The linear response solver of Jørgensen et al. 61 has been used for the solution of the response equations. Differentiated one-and two-electron integrals were computed using the Gen1Int 56,57 and Cgto-Diff-Eri 58,66 programs, respectively, except for some of the lower-order two-electron integral geometric derivatives which were computed using existing functionality in Dalton. The differentiated exchange-correlation (XC) energy and potential contributions up to fifth order needed in the DFT calculations were computed using the XCFun library, 59,60 where the integrator XCInt has been used for the integration of the XC contributions. The calculation of the Coriolis coupling constants is not done in a response theory framework, but have been calculated in the manner outlined in ref. 67. All calculations have been performed at the DFT level of theory using the B3LYP hybrid functional. [68][69][70] This functional has already been shown to give good results for the calculation of higher-order properties in earlier work. 46, 71 Dunning's correlation-consistent polarized triple-z (cc-pVTZ) basis set 72 has been used. The study was conducted for methanimine (CH 2 NH), and nitromethane (CH 3 NO 2 ) and its mono-(CH 2 DNO 2 ) and di-deuterated (CHD 2 NO 2 ) isotopomers. Two conformations (eclipsed and staggered) have been considered for the nondeuterated isotopomer and four (H-eclipsed, D-eclipsed, H-staggered and D-staggered) for each deuterated isotopomer (cf. Fig. 1).
For each system, the geometry was optimized and the molecular Hessian and the rotational constants were computed using the Dalton2013 program package. 64,65 The other relevant molecular properties were computed at the optimized geometry using the recursive response property implementation, and the Coriolis coupling constants have been implemented in a development version of Dalton2013. The molecular Hessian was then used in a vibrational analysis to find the harmonic vibrational frequencies and to transform the geometric differentiation in the property tensors from a Cartesian basis to a reduced normal coordinate basis 73 to calculate anharmonic frequencies and spectral intensities.
Anharmonic corrections to the fundamental frequencies, as well as first overtones and combination band frequencies were calculated from the cubic and quartic force constants, the rotational constants and the Coriolis coupling constants using a scheme based on vibrational second-order perturbation theory 35,36 as described in Section 2.2, where terms found to be affected by Fermi resonances are taken out of the perturbational treatment 23 and resolved variationally 41 using the GVPT2 model. 38 First-order geometric derivatives of the electric dipole and electric dipole polarizability in the reduced normal coordinate basis were used for the evaluation of the harmonic IR intensities and Raman scattering cross-sections, respectively. Anharmonic corrections to the spectral intensities were calculated by further considering the second and third geometric derivatives of the corresponding properties and the cubic and quartic force constants, in a reduced normal coordinate basis, using the GVPT2 model, resulting in features associated with corrections to the fundamental bands and the appearance of the first overtone and combination bands.
For methanimine, the cubic and quartic force fields have also been evaluated by numerical differentiation from the molecular Hessians calculated for Cartesian displacements dx of 10 À2 and 10 À3 Å with Dalton2013 using the expressions where E xixj , E xixjxk and E xixjxkxl represent, respectively, the second-, third-and fourth-order derivatives of the energy with respect to the Cartesian components in superscript, and using convergence thresholds of 10 À8 for both the molecular orbital (MO) coefficients and relative to the norm of the perturbed MO coefficients when solving the response equations. The same convergence criteria have been applied to all fully analytic calculations. We remark that the errors in the calculated properties resulting from these strict thresholds are negligible. The first, second and third geometry derivatives of the electric dipole moment and polarizability have also been evaluated this way and with the same convergence thresholds, but using the expressions where P denotes either the electric dipole moment or the electric polarizability, and P xi , P xixj and P xixjxk represent respectively the first, second and third derivatives with respect to geometry distortions. The spectral bands have been modeled using Lorentzian functions for the band shape with a 10 cm À1 full width at half maximum. A 1 cm À1 resolution was used to plot all spectra. Raman spectra have been evaluated considering an incident laser wavelength of 514 nm, corresponding to an Ar + laser at 298.15 K.

Results and discussion
4.1 Reliability of the approach: methanimine In this section, we will illustrate the need for analytic differentiation techniques by calculating the infrared and Raman spectra of methanimine (CH 2 NH), comparing the analytic approach to the results obtained by numerical differentiation using different step lengths. The sensitivity of methanimine to numerical differentiation parameters 52 makes it a suitable system for illustrating the advantages of using an analytic approach.
The theoretical vibrational frequencies obtained using the different approaches are compiled in Table 1. Experimental values 74,75 are also presented for comparison.
In the case of numerical differentiation using a step length of dx = 10 À3 Å, the difference in Hessian values between some of the displaced systems was smaller than the numerical precision, thus illustrating one of the problems of this approach. This can be illustrated by the anharmonic correction to the vibrational frequency of the 9 1 mode which is positive, whereas anharmonic corrections are generally expected to be negative, as is obtained in the analytic approach and when a step length of dx = 10 À2 Å is used in the numerical differentiation approach. The anharmonic corrections to the vibrational frequencies of the high-frequency modes appear less sensitive to this problem. The analytic approach does not depend on the energy difference between slightly displaced systems and is therefore free from this source of numerical error.
Using a step length of dx = 10 À2 Å for the numerical differentiation, numerical noise is largely avoided and the anharmonic frequencies are in better agreement with experimental fundamental frequencies. This is also observed for the anharmonic frequencies obtained by analytic differentiation. Nevertheless, the numeric anharmonic corrections are still on average in error by about 10% compared to the analytic corrections, the latter being always larger than the former. Fig. 2 and 3 show, respectively, the calculated infrared and Raman spectra of methanimine for the analytical and numerical approaches. In the calculated infrared spectrum, using a step length of dx = 10 À3 Å in the numerical differentiation, not only are the anharmonic corrections to the frequencies in poor agreement with the analytic ones, but so are also the corrections to the intensities, most strikingly so for the lowfrequency peaks. In this case, for the IR spectrum, the dx = 10 À3 Å numerical differentiation reproduces the analytic anharmonic spectral intensities almost perfectly for the peaks of frequency above 2900 cm À1 but overestimates (in absolute value) drastically the intensity for the other peaks, the lower the frequency the larger the overestimation.
Numerical (dx = 10 À2 Å) and analytic anharmonic corrections to the spectral intensities both go in the same direction for each individual peak, but the magnitude of the corrections differs. The difference in the intensity of the anharmonic corrections to the infrared intensities between the numerical and analytic values varies from 10 to 230% of the analytic correction depending on the peak considered, with the majority   of the corrections being in error by 35-85%, the only exceptions being the low-energy modes 8 1 and 9 1 . However, there is no trend as to whether the numerical corrections under-or overestimate the analytic results. As the anharmonic corrections to the total intensity of the peaks is small, these differences are not easily visible from the spectra plotted in Fig. 2.
For the Raman spectra, numerical noise does not affect the derivatives of the electronic polarizability when using a step length of dx = 10 À3 Å. The numerical anharmonic corrections to the spectral intensities are thus in better agreement with the analytic ones than in the infrared case, but the corrections to the vibrational frequencies remain wrong. Considering the intensities, the numerical spectrum obtained using a step length of dx = 10 À3 Å shows differences of less than 10% compared to the analytic spectrum, and is thus in better agreement than the spectrum obtained using a step length of dx = 10 À2 Å, where these differences may be as large as 20%. The only exception is the 8 1 mode, for which both step lengths give corrections that are far from the analytic one. As for the IR spectra, the calculated corrections can be both larger and smaller than the analytic result and whether the corrections are over-or underestimated also depends on the step length. It should also be noted that, depending on the step length used, the ordering of the intensity of the peaks can differ. For example, in the case of dx = 10 À2 Å, 1 1 is slightly more intense than 2 1 , whereas with dx = 10 À3 Å, the 2 1 peak is more intense than 1 1 , in agreement with the analytic differentiation results.
This example illustrates that even if the use of numerical differentiation can lead to qualitatively sound results, it still depends strongly on the step length used. While methanimine is still a rather small molecule, it could still be expected that these difficulties will be present in larger systems. On this note, we now turn our attention to using the analytic approach to calculate anharmonic vibrational spectra and compare these with available experimental observations.

Comparison with experiment: nitromethane
In Section 4.2.1, we will present and discuss the computed vibrational frequencies, before we in Sections 4.2.2 and 4.2.3 turn to a discussion of the theoretical IR and Raman spectra, respectively, comparing our theoretical results to available experimental data.
All experimental and theoretical studies 76-79 on the geometry of nitromethane agree that the barrier (DE = 9.6mE h 76 ) for the rotation of the methyl group around the CN axis is very small, with the staggered conformation being slightly more stable. Our results reproduce quantitatively the barrier height (DE B3LYP = 10mE h ). Such a low barrier makes it necessary to consider several rotamers when modeling the theoretical spectra, and for this reason all the geometries corresponding to the extrema of the energy along the rotation of the methyl group are considered in this study (cf. Fig. 1). A Boltzmann averaging at room temperature of these rotamers would give a quasi-equal weight for each of the conformers, and for this reason all rotamers will thus be considered of equal weight in the averaging of the spectra from the different rotamers. We note that such a treatment for the low-frequency internal rotation of the methyl group has to be considered approximate, and that this vibration mode probably should be treated by a non-local representation going beyond the normal-mode approximation. For this reason, we will in the following not include this mode in the anharmonic treatment.

Vibrational frequencies.
For all rotamers of each isotopomer, the frequency corresponding to the rotation of the methyl group is found to be quite small at the harmonic level and negative at the anharmonic level, which is consistent with what can be expected for a quasi-free rotating methyl group. 76,77 The observed spectra should therefore come from the average over all the rotamers. For this study, only the extremum rotamers (staggered and eclipsed) have been considered (cf. Fig. 1), and the system has been treated as having only 14 normal modes (instead of 3N À 6 = 15) by not considering the derivatives with respect to the methyl rotation mode in the anharmonic calculations.
Using partially deuterated isotopomers lowers the symmetry of the system, thus allowing new rotamers to be spectroscopically active and giving rise to band splittings. 77,78 Calculated (harmonic and anharmonic) frequencies for the fundamentals of the non-, mono-and di-deuterated isotopomers of nitromethane are compiled in Tables 2, 3 and 4, respectively. Experimental frequencies 77,[79][80][81] are also given for comparison.
The computed harmonic fundamental frequencies are, in line with previous findings, 79 found to be overestimated compared to experiment. Even with anharmonic corrections, the frequencies are in many cases overestimated compared to the experimental data, but lead to a significantly better agreement with experiment. The differences in the calculated vibrational frequencies for the different rotamers are in general very small. Indeed, very similar vibration frequencies are found for the two rotamers of CH 3 NO 2 at both the harmonic and anharmonic level of calculation, the largest difference being 9 cm À1 . The calculated vibrational frequencies are also in very good agreement with the experimental assignments of the modes. [80][81][82][83]  Experimentally, two peaks are assigned in the infrared spectrum to the stretching mode of the C-D bond in CH 2 DNO 2 : A strong mode at 2266 cm À1 corresponding to the stretching perpendicular to the plane of the nitro group, and a weak one at 2276 cm À1 corresponding to the stretching parallel to the same plane. 79 We find that the maximum frequency for the n(CD) mode is found for the D-eclipsed geometry in both the harmonic and anharmonic treatment, and the frequency decreases the further away the deuterium atom is from the NO 2 plane.
A similar behaviour is also observed for the stretching of the CH bond in CHD 2 NO 2 , with a strong band at 3029 cm À1 corresponding to stretching perpendicular to the plane of the nitro group, and a weak band at 3014 cm À1 corresponding to stretching parallel to the same plane. 79 The maximum frequency for the n(CH) mode is found for the H-eclipsed geometry in both the harmonic and anharmonic treatment, and the frequency then decreases the further away the hydrogen atom is from the NO 2 plane, in analogy to the observations for CH 2 DNO 2 .

Infrared spectra.
The infrared spectrum is calculated by summing the calculated spectra of the three rotamers. The harmonic and anharmonic calculated infrared spectra for non-, mono-and di-deuterated isotopomers are shown in Fig. 4-6, respectively. Tables 5-7 show the calculated infrared spectral intensities (before Lorentzian normalization) for the normal modes of the non-, mono-and di-deuterated isotopomers.
For the three isotopomers considered in this study, the anharmonic corrections do not substantially change the relative intensities of the different fundamental bands below 2000 cm À1 . In this region, the main improvements arising from the anharmonic treatment is in the calculated vibrational frequencies, as discussed in Table 4 Calculated (harmonic and anharmonic) and experimental normal vibrational frequencies of CHD 2 NO 3 (in cm À1 )   Table 3 Calculated (harmonic and anharmonic) and experimental normal vibrational frequencies of CH 2 DNO 2 (in cm À1 )  the previous section. This observation is in agreement with the findings of Bloino and Barone 38 using the GVPT2 approach with numerical calculation of the anharmonic IR spectra for a series of molecules.
For all three isotopomers, a weak feature appears in the anharmonic spectrum around 2940 cm À1 arising mainly from the combination of the symmetric (mode 4 for all isotopomers) and asymmetric (mode 7 for CH 3 NO 2 , mode 6 for CH 2 DNO 2 and mode 5 for CHD 2 NO 2 ) stretching modes of the NO 2 fragment.
The experimental gas-phase spectrum of CH 3 NO 2 from ref. 84 is reproduced in Fig. 4 for comparison. As already noted, a low-intensity feature around 3000 cm À1 , also appearing in the experimental spectrum, is introduced with the anharmonic treatment. The main peak of this feature, at 2955 cm À1 , arises mainly from the 4 1 7 1 combination band and a minor contribution from the 3 1 fundamental band. Apart from the 1 1 and 2 1 fundamental bands (3038 and 3005 cm À1 , respectively), another low-intensity combination band, 4 1 8 1 at 2933 cm À1 , appears from the anharmonic treatment. Other low-intensity peaks, also present in the experimental spectrum, appear due to the 4 1 11 1 combination band at 2471 cm À1 and the 7 2 overtone band at 2768 cm À1 .
In the CH 2 DNO 2 spectrum, a small shoulder arising from the 10 1 12 1 combination band of the H-eclipsed conformer and the 9 1 13 1 combination band from the D-staggered conformer appears on the most intense peak. These two bands are combinations of an angular vibration of the NO 2 fragment and an angular motion of the CD fragment. A low-intensity band appears from the anharmonic treatment around 3000 cm À1 . The main peak of this band, around 2933 cm À1 , corresponds to the 4 1 6 1 combination band from all four conformers. The rest of the features of this band arise from the 1 1 and 2 1 fundamental bands of the four conformers. Other low-intensity peaks appear in the anharmonic spectrum around 2738 cm À1 due to the 6 2 overtone band of the four conformers and around 2462 cm À1 due to the 4 1 11 1 combination band from the H-eclipsed and D-staggered conformers.
In the CHD 2 NO 2 spectrum, the anharmonic treatment gives a modification in the shape of the (broad) band between 1230 and 1280 cm À1 , mainly due to the 6 1 and 7 1 fundamental bands   and from the 12 2 overtone band. A weak band appears around 1505 cm À1 , due to the 10 1 12 1 and 11 1 12 1 combination bands from the four conformers. A low-intensity peak, corresponding to the 4 1 5 1 combination band, appears around 2930 cm À1 from the anharmonic treatment. In addition to the peak at 3020 cm À1 , corresponding to the 1 1 fundamental band of the H-eclipsed conformer, another low-intensity peak appears in the anharmonic spectrum around 2740 cm À1 and is due to the 5 2 overtone band.

Raman spectra.
As done for the IR spectra, the Raman spectra were obtained as the sums of Raman spectra for the individual rotamers. The calculated harmonic and anharmonic Raman spectra for non-, mono-and di-deuterated isotopomers are shown in Fig. 7-9, respectively. Tables 8-10 show the calculated Raman spectral intensities (before Lorentzian normalization)   for the normal modes of the non-, mono-and di-deuterated isotopomers.
The relative intensities of the bands corresponding to CH or CD vibrations compare rather well with the experimental data, as well as the relative intensities of the bands corresponding to the CN and NO 2 motions. However, the agreement between theory and experiment for the relative intensities of these two different vibrations is poor. The anharmonic treatment gives slightly better agreement with experiment, though the differences are very small.
The major correction arising from the anharmonic treatment occurs for the frequencies, as also noted for the IR spectra. Anharmonic corrections do not modify significantly the relative intensities of the fundamental bands, except for the band at 1380 cm À1 corresponding to the n s (NO 2 ) mode that is weakened by the anharmonic treatment for all three isotopomers relative to the harmonic model. However, for both of the deuterated isotopomers, a slight improvement in the intensities of the band corresponding to the 3 1 mode for the mono-deuterated and the 2 1 mode for the di-deuterated isotopomer (both corresponding to a stretching of the CD bonds) seems to appear from the anharmonic treatment, where they are reduced slightly more than the other bands that correspond to hydrogen-related motions. This modification of the relative intensities leads to somewhat better agreement with experiment.
We note that the liquid-phase IR spectrum for CH 3 NO 2 (as can be seen in Fig. 1(a) of ref. 85, not reproduced here) presents strikingly different spectral intensities to the gas-phase spectrum. Therefore, a corresponding effect can be expected also for the Raman spectra. This could explain at least some of the differences between the calculated (gas-phase) and experimental (liquidphase) spectra presented here, and solvent effects will be investigated in future work.
The experimental liquid-phase spectrum of CH 3 NO 2 from ref. 85 is reproduced in Fig. 7 for comparison. A weak peak appears around 1280 cm À1 corresponding to the 12 2 overtone. The small shoulder appearing in the experimental spectrum on the low-frequency side of the peak around 1400 cm À1 could be related to this overtone band. Other very weak peaks appear between 2700 cm À1 and 2850 cm À1 and are due to several different overtones and combination bands. These features could be related to the weak peak appearing in the experimental spectrum around the same frequency. A large number of features in the theoretical spectrum appear between 1600 and 2900 cm À1 and arise from the contributions of different combination and overtone bands.
The experimental liquid-phase spectrum of CH 2 DNO 2 from ref. 81 is reproduced in Fig. 8 for comparison. A major improvement in comparison to the experimental spectrum is observed when including anharmonic corrections for the splitting of the 3 1 band with respect to the position of the deuterium atom, around 2195 cm À1 if the CD bond stretches almost parallel to the NO 2 plane and around 2295 cm À1 if it stretches almost perpendicular to it, as seen in the experimental spectrum. A weak peak appears in the anharmonic spectrum around 2250 cm À1 , corresponding to the 6 1 11 1 combination band, mainly from the D-eclipsed rotamer. Another weak feature appears around 2310 cm À1 , corresponding to the 8 1 9 1 combination band present in all the four rotamers.
The experimental liquid-phase spectrum of CHD 2 NO 2 from ref. 81 is reproduced in Fig. 9 for comparison. The band between 2250 and 2290 cm À1 , arising from the 2 1 mode in the harmonic treatment, is broadened in the anharmonic spectrum,  mainly due to the 6 1 9 1 combination of the four rotamers on the low-frequency side and to the 6 1 8 1 combination bands on the high-frequency side of the D-eclipsed and H-eclipsed rotamers. Two weak peaks appear around 1920 and 2060 cm À1 due, respectively, to the 8 2 and 9 2 overtones, primarily from the D-staggered rotamer. A very weak band also appears around 2530 cm À1 , corresponding to the 6 2 and 7 2 overtones, mainly from the D-staggered rotamer.

Summary and concluding remarks
We have presented the first analytic calculations of anharmonic corrections to both the vibrational frequencies and intensities in infrared and Raman spectroscopies. This has been made possible by our recent development of a recursive scheme for the calculation of high-order molecular properties, including properties involving frequency-dependent perturbations and perturbation dependence in the basis set. 47 The approach is applicable to single-determinant self-consistent field models such as Hartree-Fock theory and Kohn-Sham DFT, and being matrix-based, it can also be extended to linear-scaling approaches, 48,62 as well as to the relativistic four-component level of theory. 86 We have previously applied our approach to the calculation of anharmonic corrections to vibrational frequencies using density functional theory 46 and to the analytic calculation of Raman optical activity 87 and hyper-Raman scattering, 71 and in this work, we have used it to calculate anharmonic infrared and Raman spectra of nitromethane and its partially deuterated isotopomers. We find that anharmonic effects lead to an improvement in the quality of the computed IR and Raman spectra. The major improvements arising from the anharmonic treatment occur for the vibrational frequencies, while the effects of the anharmonic corrections on the infrared and Raman intensities are smaller and show only minor influences on the relative intensities of the fundamental bands, which is in agreement with earlier observations. 38 Nevertheless, the anharmonic corrections are important in order to capture the overtone and combination bands. The anharmonic corrections are found to be somewhat more important for the Raman spectra, even if very small, than for the infrared spectra. Overall, the anharmonic spectra are in better agreement with experiment than the corresponding harmonic spectra. We have also shown that evaluating the energy and property derivatives by numerical differentiation is prone to numerical instabilities, as also noted elsewhere, 51 so that obtaining reliable numerical derivatives can prove difficult for general molecular systems, and we have seen that the errors thus introduced can significantly affect the calculated results, whereas analytic approaches would be free of these sources of error. Because of this, we believe that analytical derivatives of high order is an important step in making the inclusion of anharmonic corrections in calculated infrared and Raman spectra routine, leading to an improved understanding of the importance and occurrence of anharmonic effects in vibrational spectroscopies.
Finally, solvent effects are known to affect vibrational spectroscopies. 88 It is therefore important that the scheme presented here is extended to include solvent effects, either in the form of polarizable continuum models or polarizable embedding approaches, 89 and work in this direction is in progress. 90