Stretching of cis-formic acid: warm-up and cool-down as molecular work-out

Abstract

A new technique to rotationally simplify and Raman-probe conformationally and vibrationally excited small molecules is applied to the cis–trans isomerism of formic acid. It quintuples the previously available gas phase vibrational data base on this excited form of a strongly anharmonic planar molecule despite its limited spectral resolution. The newly determined cis-formic acid fundamentals allow for a balanced vibrational benchmark on both rotamers of formic acid. Assuming the adequacy of vibrational perturbation theory, it reveals weaknesses of standard methods for these systems like B3LYP-D3(BJ)/aVQZ VPT2 or PBE0-D3(BJ)/aVQZ VPT2. The functionals ωB97-XD and M06-2X additionally suffer from severe integration grid size and symmetry dependencies. The vibrational benchmark suggests B2PLYP-D3(BJ)/aVQZ VPT2 and MP2/aVQZ VPT2 as partially competitive and in any case efficient alternatives to computationally demanding coupled cluster vibrational configuration interaction calculations. Whether this is due to fortuitous compensation between electronic structure and vibrational perturbation error remains to be explored.

---

1 Introduction

Vibrational spectra of small molecules effectively probe the quality of potential energy hypersurface (PES) predictions, when combined with accurate anharmonic calculations.¹ Typically, an intense interplay between theory and experiment initially converges the performance for a set of low quantum number states around the global minimum. To explore the globality of a PES, it is then desirable to add experimental data on a secondary minimum structure. Its quantum states start locally, but evolve into mixed structure states at higher excitation, probing the transition state region as well. For three atoms, HCN/HNC is the paradigmatic example.² For four atoms, the simultaneous description of the nearly isoenergetic formaldehyde molecule and H₂–CO complex is challenging.^3,4 For five atoms, the cis–trans isomerism of formic acid is arguably one of the most interesting systems, calling for suitable experimental reference data for the higher-energy cis-form. These have been surprisingly scarce until very recently, with a single exception.⁵

As the smallest carboxylic acid, the formic acid monomer has been addressed by a plethora of theoretical^6–15 and experimental studies.^5,16–55 The isomerisation from the ground state trans-form to the higher-energy cis-form has been of particular interest,^{5,14,22,41,42} even when looking at processes in the interstellar medium.⁵⁶ When it comes to the vibrations of cis-formic acid, matrix isolation has been the method of choice thus far, because the possibility of long irradiation times allows for a significant formation via laser excitation of the trans-form.^{35,41,42,49,50,57} Since the matrix environment shifts the band positions compared to the gas phase, a direct comparison with predicted band positions of modern quantum chemical methods would require a very challenging description of the environment. Accurate theoretical predictions for the isolated cis-isomer thus suffer from a lack of gas phase experimental reference data. Two recent studies where this applies are by Tew and Mizukami¹⁴ from 2016 and by Richter and Carbonnière¹⁵ from 2018.

Due to the fairly large energy difference of 1365 ± 30 cm⁻¹ between both rotamers of formic acid,²² vibrational gas phase data on the cis-form are rare. The first gas phase band position of cis-formic acid has been published in 2006 by Baskakov and co-workers, who studied the out-of-plane bending vibration with high resolution FTIR spectroscopy.⁵ Only very recently, it was complemented by a second example obtained as a side effect when studying excitonic C=O stretch coupling in jet-cooled carboxylic acid dimers.⁵⁴ This observation has triggered the present work, which represents a systematic study of all four valence stretching modes of cis-formic acid. It is based on a powerful new Raman scattering approach of thermally populated and rapidly re-cooled molecules. Instead of conserving the conformational excitation by cryogenic matrix trapping,⁵⁸ the spectra are rotationally simplified by supersonic expansion. Vibrational and high barrier conformational excitation is largely trapped and can be probed without environmental distortion as a function of initial gas temperature. Back-tunnelling to the trans-form is also not an issue on the time scale of the jet expansion, making it an “easy” experiment.⁵⁹ By a 400% increase of perturbation-free cis-formic acid vibrational frequencies after a decade of stagnation, we provide the first systematic access to the performance of quantum chemical methods towards this model system. This decreases the likelihood of accidental error compensation between electronic structure, vibrational treatment, and matrix shifts for cis-formic acid by orders of magnitude.

The structure of this paper is as follows: we briefly illustrate the general approach of how the spectra of cis-formic acid were recorded, followed by a detailed analysis of the spectra and a first benchmark of vibrational perturbation theory and literature variational data against the experimental data. It is hoped that this progress will trigger further growth in the experimental data base and its use in benchmarking the global PES of formic acid and pentatomic vibrational treatments.

2 Methods

2.1 Experimental

A detailed description of the experimental set-up can be found in previous publications.^54,60 Formic acid (Acros Organics, 98+%) was seeded at 0.2% in helium and expanded through a vertical slit nozzle at 1.0 bar into the evacuated jet chamber (background pressure < 0.1 mbar). The expansion was probed by a 532 nm, 25 W continuous-wave laser from below. The scattered radiation was detected perpendicularly with respect to both the expanding flow and the laser with a monochromator equipped with a charge-coupled device camera. cis-Formic acid was formed in small quantities from the trans-rotamer by heating the nozzle and its feed line to temperatures between 100–190 °C.

2.2 Computations

The quantum chemical calculations shown in this work were performed with the Gaussian 09 program package (revision E.01)⁶¹ using a pruned ultra fine integration grid (99 590). This grid is finer than the default of Gaussian 09 (fine grid, (75 302)).⁶¹ The employed methods are B3LYP,^62,63 B2PLYP,⁶⁴ MP2,^65,66 M06-2X,⁶⁷ ωB97-XD,⁶⁸ PBE0,⁶⁹ and PM3.^70,71 Grimme's pairwise dispersion correction (D3) in combination with Becke–Johnson (BJ) damping has been added for B3LYP, B2PLYP, and PBE0.⁷² For all methods, an augmented quadruple-zeta (aVQZ) basis set has been chosen. Additional augmented double-(aVDZ) and triple-zeta (aVTZ) calculations were carried out for MP2 and B2PLYP-D3(BJ). Moreover, the production calculations have been performed without the use of symmetry, employing opt = tight convergence.

The assignment of cis-formic acid fundamentals has been supported by scaled, harmonic frequency calculations at the B3LYP-D3(BJ)/aVTZ level, which have proven to yield sufficient agreement in a previous study.⁵⁴ For the vibrational benchmark in Section 3.2, anharmonic frequency calculations were performed at all levels listed above using vibrational perturbation theory (VPT2),⁷³ as implemented in Gaussian 09.⁶¹ VPT2 was used under the default settings where resonances identified in a pre-screening are removed and treated variationally.

Additionally, exploratory VPT2 calculations utilising the C_s symmetry as well as a finer integration grid (pruned super fine integration grid (150 974)⁶¹) were carried out in selected cases to probe their impact on the results. A brief discussion can be found in Section 3.3.

3 Results and discussion

3.1 The stretching vibrations of cis-formic acid

To choose suitable spectral regions for cis-rotamer detection, the band positions and Raman scattering cross-sections have been predicted using B3LYP-D3(BJ)/aVTZ alongside those of the trans-form. The results are displayed in Fig. 1. The vibrations have been labelled according to the Herzberg nomenclature. The cis-formic acid vibrations with the largest scattering cross-sections are ν₁, ν₂, ν₃, and ν₆, namely the O–H, the C–H, the C=O, and the C–O stretching vibration. In fact, ν₆ is the only stretching vibration with a distinctly larger scattering cross-section compared to trans-formic acid.

Fig. 1 Harmonically predicted band positions and Raman scattering cross-sections σ (ref. 74) of all nine fundamentals of cis- (violet) and trans-formic acid (black), calculated at the B3LYP-D3(BJ)/aVTZ level.

The experimental spectra of these four vibrational modes of both rotamers (cF, F) can be found in Fig. 2 alongside harmonic, individually F-scaled B3LYP-D3(BJ)/aVTZ calculations below the spectra. For each spectral region, four spectra with increasing nozzle temperature have been recorded. These temperature series have been intensity-scaled to the trans-monomer band of lowest intensity amongst the four. Consequently, any hot band, i.e., cis-formic acid or a non-isomeric hot band originating from thermally populated low-lying energy levels of trans-formic acid, should increase in intensity with nozzle temperature, whereas any formic acid cluster band will decrease due to thermal dissociation.

Fig. 2 Raman jet spectra of formic acid (≈0.2%) in helium recorded at a reservoir pressure of 1.0 bar with increasing nozzle temperatures of 100–190 °C. Four spectral regions are covered at four temperatures with data acquisition times between 28 and 63 minutes. Within a spectral region, the spectra have been intensity-scaled to the trans-formic acid monomer band F with the lowest concentration amongst the four temperatures. Below each temperature series, calculated, v_i(F)-scaled harmonic band positions (B3LYP-D3(BJ)/aVTZ level) of cF and F are shown.

The spectra in the O–H stretching region show one band that increases in intensity with temperature at 3637 cm⁻¹. The band position is in good agreement with the harmonically calculated, ν₁(F)-scaled band position of cF with a deviation of only 5 cm⁻¹. Either the anharmonicity of F and cF is similar or there is error compensation with the density functional used. Another way of validating this assignment is to compare the intensity ratio of the cF (3637 cm⁻¹) and F (3570 cm⁻¹) bands with the energy difference of both forms. The harmonically calculated energy difference of 15.9 kJ mol⁻¹ (B3LYP-D3(BJ)/aVTZ with zero point energy correction) is just below the error bounds of the only experimental value of 1365 ± 30 cm⁻¹ by W. Hocking.²² Neglecting differences in the partition function of the two complexes, this corresponds to a population of 1–2% of cis-formic acid at 190 °C. After correction by the theoretical cross-section ratio, the ratio of the experimental band integrals provides a cis-abundance of 2%, thus reaffirming the cF assignment. The additional bands downshifted compared to the O–H stretching vibration of trans-formic acid at 3560 cm⁻¹ and 3566 cm⁻¹ are most likely trans-formic acid combination bands of ν₂ with the lowest frequency vibrations ν₇ (3560 cm⁻¹) and ν₉ (3566 cm⁻¹), which benefit from the large Raman scattering cross-section of the C–H stretching vibration. The former is in good agreement with the predicted values of Tew and Mizukami (3566 cm⁻¹)¹⁴ as well as Richter and Carbonnière (3558 cm⁻¹)¹⁵ and the latter with a prediction of Tew and Mizukami,¹⁴ who reported (ν₂ + ν₉) in Fermi resonance with (ν₃ + 3ν₉) at 3571 cm⁻¹ and 3579 cm⁻¹.

The ν₂ region is spectrally more congested due to its low sensitivity to hydrogen bonding. In the spectral windows 2970–2945 cm⁻¹ and 2935–2925 cm⁻¹, there are several bands that decrease in intensity with temperature, i.e., cluster bands. The broad underlying signal is due to rovibrational O and S branches of ν₂. As opposed to the O–H stretching region, there are two distinct bands increasing in intensity with temperature at 2925 cm⁻¹ and 2873 cm⁻¹. The latter deviates from the predicted band position of cis-formic acid by 14 cm⁻¹. The amount of cis-formic acid at 190 °C deduced from the integrated intensities of the bands amounts to 1%, which fits the energy difference, as detailed above. Therefore, the band at 2873 cm⁻¹ can be assigned to cF. The second hot band at 2925 cm⁻¹ is shifted by −17 cm⁻¹ from the fundamental of F (2942 cm⁻¹). For an assignment to F, two things need to be considered: firstly, the shift directly yields the off-diagonal anharmonicity constant x_2i between ν₂ and a low-lying energy level v_i that is thermally populated. Secondly, the intensity ratio is dependent on the Boltzmann population of that level and yields the excitation energy of the latter. Hence, the assignment can be checked by comparing the experimentally determined anharmonicity constant and intensity with the calculated values for the lowest-lying energy levels of trans-formic acid. From the anharmonicity matrix elements in Table 1 it is apparent, that the hot band originating from ν₇ will most likely overlap with the fundamental, whereas the hot band originating from ν₉ (and ν₆) could overlap with a cluster band at 2938 cm⁻¹ causing the highest nozzle temperature spectrum (red) and the lowest nozzle temperature spectrum (black) to have similar intensities. However, due to the spectral congestion in this area, reliable assignments are not feasible. Additional depolarisation measurements to subtract the O and S branches from the sharp Q peak are currently ongoing and will be addressed in detail in a subsequent publication. Here we focus on the straightforward assignment of the 2925 cm⁻¹ band. Its observed shift of −17 cm⁻¹ perfectly matches the calculated anharmonicity constant x₂₈. The expected intensity ratio at 190 °C of around 4% approaches the observed ratio of 3%, so that it can be assigned to ν₂ + ν₈ − ν₈.

Table 1 Calculated anharmonic (VPT2) band positions (in cm⁻¹) of trans-formic acid alongside calculated diagonal (italics) and off-diagonal anharmonicity constants x_2i, x_3i, and x_6i (in cm⁻¹) of ν₂, ν₃, and ν₆ with all nine fundamentals

v i	Band position	B3LYP-D3(BJ)/aVTZ VPT2
		x 2i	x 3i	x 6i
1	3533	−3.5	−1.1	−0.9
2	2892	−61.9	−13.4	−2.3
3	1779	−13.4	−9.3	−4.1
4	1374	−20.4	−0.3	−6.2
5	1219	−6.1	+2.8	−14.3
6	1089	−2.3	−4.1	−6.2
7	622	+0.4	−6.2	−3.7
8	1031	−17.1	−5.4	−3.7
9	643	−2.8	−0.4	−5.6

---

A first analysis of the ν₃ spectral region at nozzle temperatures of 23 °C, 110 °C, 140 °C, and 170 °C can be found in a previous publication.⁵⁴ Omitting the room temperature measurement in this work is intentional, as it is heavily congested with cluster bands due to the fairly high concentrations and reservoir pressures chosen. These, however, are essential to obtain high monomer signals of F, and especially cF, after thermal dissociation of the clusters. Due to an improved signal-to-noise ratio compared to the previous measurements as well as the somewhat higher upper nozzle temperature available in the present work, we have reanalysed the ν₃ region. Briefly, the cis-formic acid band can be seen at 1818 cm⁻¹, which deviates from the calculated, ν₃(F)-scaled band position by 3 cm⁻¹. The hot band downshifted by 7 cm⁻¹ from F can most likely be attributed to ν₃ + ν₇ − ν₇, with a negligible (1 cm⁻¹) deviation of the calculated anharmonicity constant x₃₇ compared to the experimentally observed value and a reasonable Boltzmann population match (14% from the level energy and 10% from the Raman spectrum). The hot band intensity qualitatively rules out major contributions from higher energy levels such as ν₈. There are two weaker potential hot bands shifted from F by +6 cm⁻¹ and −13 cm⁻¹ with intensities of around 1–2% compared to F. An assignment is not possible since the shifts do not match the predicted anharmonicity constants (cf. x₃₉, x₃₈, x₃₆, and x₃₅ in Table 1). As previously seen, the intensity ratio gives only a rough estimate of the energy level and as such, cannot serve as a stand-alone assignment criterion. Overall, this highlights the importance of the ν₂ region with its much higher monomer signal due to the large Raman scattering cross-section of the C–H compared to the C=O stretching vibration.

In the C–O stretching region, three hot bands can be seen downshifted from the fundamental ν₆ of trans-formic acid at 1105 cm⁻¹. The shifts amount to −4 cm⁻¹, −7 cm⁻¹, and −11 cm⁻¹ with intensities of around 7%, 3%, and 7% of ν₆ at 190 °C. To assign ν₆ of cis-formic acid, the shifts are compared to the calculated anharmonicity constants x_6i in Table 1. The predicted anharmonicity constants x₆₇ and x₆₈ agree (−3.7 cm⁻¹). In addition, x₆₉ and x₆₆ are very similar (−5.6 cm⁻¹ and −6.2 cm⁻¹). Therefore, it seems likely that the bands at 1101 cm⁻¹ and 1097 cm⁻¹ are a result of overlapping hot bands. The slightly higher intensity of the former is a result of the greater overlap with the fundamental and the lower energy of ν₇ and ν₈ compared to ν₉ and ν₆. The next higher energy level is ν₅ with a predicted band position of 1219 cm⁻¹. A hot band originating from ν₅ is expected to be shifted by −14.3 cm⁻¹ from ν₆ of trans-formic acid (cf.Table 1), which is close to the experimentally observed shift of −11 cm⁻¹ of the third hot band. However, the intensity of that band is with 7% of ν₆ distinctly larger than the expected 2% from thermal population at 190 °C, especially considering that the observed intensities of all other hot bands are smaller than or equal to the predicted values. Hot bands from higher energy levels can therefore also be ruled out as these should have even lower intensities. The predicted band position of cis-formic acid deviates by −5 cm⁻¹ from the band at 1093 cm⁻¹, which falls within the accuracy of the ν₆(F)-scaled harmonic B3LYP-D3(BJ)/aVTZ calculations. Additionally, the observed intensity matches the calculated energy difference between both rotamers, considering the four times larger predicted scattering cross-section of ν₆ of cis-formic acid compared to the trans-form (cf.Fig. 2). Consequently, the band at 1093 cm⁻¹ can be predominantly assigned to the C–O stretching vibration of cis-formic acid.

The band positions of all stretching vibrations of cis-formic acid as well as that of the out-of-plane O–H bending vibration (ν₉) determined from high resolution FTIR measurements⁵ are summarised in Table 2 in comparison to the values obtained in an argon matrix by Maçôas and co-workers.⁴² The argon matrix shifts range from +27 to −20 cm⁻¹ or +23 to −21 cm⁻¹, dependent on the matrix site. This scatter is of a similar order of magnitude as the cis–trans spectral differences themselves, which are also listed in Table 2. It is therefore evident that band positions in a perturbation-free environment are crucial for a direct comparison with theory values such as those of Tew and Mizukami¹⁴ and Richter and Carbonnière.¹⁵

Table 2 Band positions of the fundamentals v_i of cis-formic acid (in cm⁻¹) obtained in a supersonic jet expansion probed with Raman spectroscopy (this work) compared to literature values

v i	Jet (this work)	Gas phase^5	Shift^a (cF–F)	Ar matrix^b^,^42	Matrix shift^c
a gas phase band position shifts from the trans- to the cis-rotamer. b site 1/site 2 are two dominant trapping sites. c band position shift of the values of both matrix sites compared to the gas phase.
1	3637		+67	3617.2/3615.9	−21/−21
2	2873		−69	2899.5/2896.3	+27/+23
3	1818		+41	1808.0/1806.9	−10/−11
4				/1391.8
5				1243.4/1248.9
6	1093		−11	1107.3/1104.6	+14/+12
7				662.3/662.3
8
9		493.420509(7)	−147.30467(1)	502.9/505.3	−9.5/−11.9

---

3.2 Vibrational benchmark

So far, the band positions of cis-formic acid have been compared to v_i(F)-scaled harmonic band positions calculated at the B3LYP-D3(BJ)/aVTZ level, which has shown to be quite valuable in supporting the assignment. The small size of the formic acid monomer and its structural rigidity enable anharmonic vibrational perturbation theory calculations (VPT2),⁷³ which have proven to be robust for the trans-formic acid monomer at various levels of theory in a study of the trans-formic and -acetic acid monomers and their nitrogen clusters.⁵³ This is less the case for the trans-acetic acid monomer, where the presence of the floppy methyl groups resulted in instabilities such as a wavenumber increase of the lowest frequency vibration compared to the harmonic case. The newly determined band positions of cis-formic acid thus enable a significantly extended VPT2 benchmark involving both rotamers, which should not suffer from such methyl torsion instabilities.

The experimental values that will be employed in the benchmark are the five band positions of cis-formic acid as well as the band position difference between the cis- and trans-fundamentals. The methods tested are the same as in ref. 53, namely B3LYP-D3(BJ), B2PLYP-D3(BJ), MP2, PBE0-D3(BJ), ωB97-XD, M06-2X, and in addition also PM3, all as implemented in Gaussian 09.⁶¹ For all methods, an augmented quadruple-zeta (aVQZ) basis set has been used. Additional augmented double-zeta (aVDZ) and triple-zeta (aVTZ) calculations have been performed for MP2 and B2PLYP-D3(BJ). The benchmarking plots can be found in Fig. 3. The accuracy of the band positions with the Raman set-up used in this work is about 1 cm⁻¹ (ref. 60) and the full width at half maximum of all bands is around 2 cm⁻¹, leading to a conservative error estimate of ±2 cm⁻¹ for the band positions and twice the amount for the shift. The green ellipsis in each plot thus shows the area that is in acceptable agreement with experiment. Since ν₉ has been measured with high resolution FTIR spectroscopy with a precision on the order of ∼10⁻⁵ cm⁻¹ and a somewhat lower accuracy,^5,44 a green arrow points towards the exact band position and shift in the bottom panel of Fig. 3. Additionally, the results of Tew and Mizukami (T & M)¹⁴ and Richter and Carbonnière (R & C)¹⁵ have been included for all vibrations. Briefly, Tew and Mizukami have fitted a potential energy surface based on 17 076 energies calculated at the CCSD(T)(F12*)/cc-pVTZ-F12 level. The vibrational levels were obtained by using vibrational configuration interaction (VCI) with an internal coordinate path Hamiltonian for the isomerisation path connecting both rotamers.¹⁴ Richter and Carbonnière have constructed a valence coordinate potential energy surface at the CCSD(T)-F12a/aug-cc-pVTZ level and carried out the vibrational energy calculations with the improved relaxation multi-configuration time-dependent Hartree (MCTDH) method.¹⁵

Fig. 3 Anharmonic (VPT2), theoretical prediction of the shift between the band positions of the cis- and trans-rotamers of formic acid for five fundamentals (ν₁, ν₂, ν₃, ν₆, and ν₉) plotted against the absolute band positions of the cis-form. The basis set size is encoded in the filling of the symbols. All black, filled symbols have been calculated with an aVQZ basis set. For the half-empty and the empty symbols, smaller basis sets have been used, namely aVTZ and aVDZ. The red triangles represent the band positions and shifts obtained from VCI calculations of Tew and Mizukami¹⁴ (T & M) and MCTDH calculations of Richter and Carbonnière¹⁵ (R & C). The green ellipsis shows agreement with experiment, whereby the size indicates the conservatively estimated experimental accuracy of ±2 cm⁻¹ for [small nu, Greek, tilde]_i(cF) and ±4 cm⁻¹ for Δ[small nu, Greek, tilde]_i(cF–F). The experimental band positions of the stretching vibrations (ν₁, ν₂, ν₃, ν₆) have been obtained in this work and the values of the torsional modes ν₉ of both rotamers (tip of the arrow) have been taken from ref. 5 and 44. Further details can be found in the text.

Firstly, the performance of VPT2 calculations at various levels of theory will be discussed before these will be compared with the VCI and MCTDH calculations. One should note that for the O–H stretching vibration this comparison can solely be made with the MCTDH calculations, as the theoretical band position of the cis O–H stretching vibration has not been reported by Tew and Mizukami. Due to error compensation, a better agreement between experiment and the tested methods is typically achieved for the shift between the cis- and trans-rotamers. The absolute band position is predicted correctly in two cases, namely with MP2/aVTZ and MP2/aVQZ for the C–O stretching vibration (ν₆). All other methods fail to predict the cis-formic acid band positions correctly despite generous experimental error bars for the stretching vibrations. An accurate prediction of ν₉ (and the respective shift) is evidently unrealistic due to the high accuracy of the high resolution measurements. The lower resolution Raman spectra are seen to be fully adequate to challenge theory on an absolute wavenumber scale. The vibrations where the shift is predicted within the experimental error for most methods are the C=O and C–O stretching vibrations, whereas the largest divergence is observed for the C–H stretching vibration. This is not surprising as the C–H stretching vibration is prone to stretch-bend Fermi resonance, although the VPT2 code employed⁶¹ attempts to include such pronounced resonances. Consequently, part of the discrepancy may be due to a poor vibrational description by VPT2 rather than the electronic structure calculation. The particularly drastic failure for M06-2X is caused by an inversion of the predicted energy sequence for the C–H stretch fundamental and C–H bending overtone of cis-formic acid, which is amplified by Fermi resonance. If the band labels are switched, the agreement increases significantly – the severe underestimation of the band position of −142 cm⁻¹ (Fig. 3) changes to an overestimation of +10 cm⁻¹. The band position shift improves from −186 cm⁻¹ (Fig. 3) to −34 cm⁻¹, compared to the experimental value of −69 cm⁻¹.

A comparison of the vibrationally averaged, calculated rotational constants for all methods with the experimental values for cis-formic acid obtained by Winnewisser and co-workers³⁹ is shown in Table 3. Small individual deviations on the order of ±0.5% are observed for B2PLYP-D3(BJ), B3LYP-D3(BJ), and MP2, larger deviations of up to 1–2% for M06-2X, ωB97-XD, and PBE0-D3(BJ), and very large deviations for PM3. The average deviation over all three rotational constants (last row in Table 3) supports the overall agreement with the experimental structure. The B3LYP-D3(BJ) structure shows the best agreement with a divergence of −0.1%, directly followed by B2PLYP-D3(BJ) (−0.3%). For MP2, the divergence is slightly larger as all rotational constants are underestimated and thus, do not compensate each other. The same is valid for M06-2X, ωB97D, and PBE0-D3(BJ), where all constants are overestimated.

Table 3 Relative deviations (in %) of VPT2 (aVQZ) rotational constants of cis-formic acid from the experimental values of Winnewisser and co-workers³⁹

	B2PLYP-D3(BJ)	B3LYP-D3(BJ)	MP2	M06-2X	ωB97-XD	PBE0-D3(BJ)	PM3
ΔA0/A0	+0.1	+0.4	−0.3	+1.4	+1.2	+1.1	−18.9
ΔB0/B0	−0.5	−0.3	−0.5	+0.6	+0.6	+0.6	+9.8
ΔC0/C0	−0.4	−0.3	−0.5	+0.7	+0.6	+0.6	+5.3
	−0.3	−0.1	−0.4	+0.9	+0.8	+0.8	−1.3

---

A comparison of the individual performance of the methods for the determination of [small nu, Greek, tilde]_i(cF) and Δ[small nu, Greek, tilde]_i(cF–F) clearly illustrates that there are few reliable methods. In case of PM3, this is not surprising. It is the only method that fails to predict the energetic order of the vibrations correctly with ν₄ and ν₆ switched. Other methods with particularly severe deviations from experiment are ωB97-XD (cf. ν₁ and ν₉) and M06-2X (ν₁, ν₂, and ν₉). The large underestimation of [small nu, Greek, tilde]₂(cF) and Δ[small nu, Greek, tilde]₂(cF–F) of M06-2X is enhanced by a level switch between resonance partners, as discussed above. All other methods predict the correct sequence of fundamental and overtone. Another numerical or fundamental deficiency of M06-2X/aVQZ VPT2 is the incorrect sign of the total anharmonicity of ν₉ of cis-formic acid, which gives rise to a large overestimation of the anharmonic band position (+163 cm⁻¹). In combination with an overestimation of the negative anharmonicity of ν₉ of trans-formic acid, this results in a severe overestimation of the shift (+308 cm⁻¹) between both rotamers. As such, this data point has been omitted from Fig. 3. The PBE0-D3 calculations match the experimental shifts in two cases (ν₃ and ν₆), whereas B3LYP-D3(BJ) predicts the shifts correctly in three of the five cases (ν₁, ν₃ and ν₆). Both exhibit similar deviations with respect to the band positions. Since the rotational constant prediction of B3LYP-D3(BJ) is also more accurate, it is the overall better choice. MP2 is particularly good for the description of the lower frequency modes ν₆ and ν₉ and overshoots for ν₁ and ν₂. For ν₃, an agreement with the shift is reached with the largest basis set aVQZ. It is generally rewarding that basis set sensitive methods tend to move towards the experimental region with increasing basis set size (cf.Fig. 3). Another reliable method is B2PLYP-D3(BJ), which predicts the shifts correctly in three cases (ν₁, ν₃ and ν₆) and shows only small deviations for the other two. The band positions are slightly, but consistently underestimated, apart from ν₂, where a small overestimation occurs for the larger basis sets, and ν₉, which is slightly overestimated for all basis set sizes.

The band positions and shifts obtained from the VCI calculations of Tew and Mizukami¹⁴ show good agreement with experiment. For all stretching vibrations, the band positions are overestimated and ν₉ differs solely by −1 cm⁻¹. With regard to the shifts, only one is predicted within the experimental uncertainty (ν₃), but the shift of ν₉ differs solely by about 1 cm⁻¹. Deviations are generally small and on the same order as for B2PLYP-D3(BJ)/aVQZ VPT2 or MP2/aVQZ VPT2. The agreement of the MCTDH calculations of Richter and Carbonnière¹⁵ with experiment is even slightly better. The band position shifts between both rotamers are predicted accurately for all stretching vibrations apart from ν₃, where the value is with 36 cm⁻¹ just outside the experimental confidence interval (41 ± 4 cm⁻¹). The band position of the C–H stretching vibration is predicted within the experimental accuracy and the ν₉ prediction deviates only by 2 cm⁻¹. However, the latter gas phase value was the only band position of formic acid known in the gas phase before ref. 14 and 15 were published, whereas the other cis-values were true predictions for the isolated molecule.

Another way of visualising the agreement of the theoretical predictions of Tew and Mizukami, Richter and Carbonnière, and results obtained with vibrational perturbation theory (B2PLYP-D3(BJ)/aVQZ VPT2) with experiment is shown in Fig. 4. In these three diagrams, the eight accessible deviations from experiment are plotted in units of experimental confidence interval for all four stretching vibrations in the form of octagons. Each axis connecting two vertices of the octagons corresponds to one of the four vibrations. The two directions of each axis display the two experimental observables for each vibration, namely the cis-formic acid band position (c_i) and the band position shift between cis and trans (Δ_i). The size of the deviation from experiment is encoded in the octagon size. A point on a node with the smallest octagon translates into theoretical agreement within the experimental error bars (±2 cm⁻¹ for the band position and ±4 cm⁻¹ for the shifts). Correspondingly, a point on a node with the nth octagon implies a deviation of that value by up to n confidence intervals from experiment. The predicted band position for the C=O stretching vibration of cis-formic acid by Tew and Mizukami (1824 cm⁻¹)¹⁴ deviates by +6 cm⁻¹ from the Raman jet value of 1818 cm⁻¹. Considering the experimental confidence interval of ±2 cm⁻¹, the prediction for c₃ lies on the third octagon, or in other words, three nodes away from the origin on the c₃ axis. Note that the origin in these diagrams cannot be met due to the experimental uncertainty. The sign of the deviation is illustrated by the colour shade of the symbol, whereby a dark colour shows over- and a light colour underestimation. The intermediate shade represents an indeterminate sign of the deviation, which can be seen for the shift of the C=O stretching vibration Δ₃ of Tew and Mizukami. The aforementioned consistent overestimation of the VCI calculations of Tew and Mizukami (T & M) (apart from Δ₃) can thus be directly seen by the otherwise dark-coloured symbols. The MCTDH method of Richter and Carbonnière falls closer to the origin and varies more in sign. Therefore, it shows superior agreement with experiment compared to the results of Tew and Mizukami. The tendency of the B2PLYP-D3(BJ)/aVQZ VPT2 calculations to underestimate the band position c_i as well as its ability to predict most shifts within the experimental accuracy (smallest octagon) is illustrated. Altogether, Fig. 4 nicely sums up that the MCTDH method utilised by Richter and Carbonnière offers a slightly better description of the vibrations scrutinised here. The VPT2 calculations at the B2PLYP-D3(BJ)/aVQZ level are seen to provide a less expensive alternative. This good performance of the double hybrid functional has recently been illustrated for pyruvic acid by Barone et al.⁷⁵ For formic acid, there are some interesting systematic errors, which have consequences when looking at matrix isolation spectroscopy. Superficially and surprisingly, the comparison of VPT2 anharmonic data for trans-formic acid only improves slightly when moving from a matrix to the gas phase.¹² This is largely due to substantial downshifts of polar (O–H, C=O) stretching vibrations in an Ar matrix, which mimic the underestimation of these vibrations by the B2PLYP functional in the gas phase. Such good agreements for the wrong reason must be avoided in proper benchmarking. Only the gas phase comparison can provide a realistic picture of the electronic structure performance.

Fig. 4 Agreement of the predicted anharmonic band positions of the stretching vibrations of cis-formic acid (c_1,2,3,6) as well as band position shifts between cis- and trans-formic acid (Δ_1,2,3,6) with experiment. The octagon size displays the deviation from experiment in units of experimental confidence interval, i.e., the smallest octagon represents agreement within the error bars of ±2 cm⁻¹ for c_i and ±4 cm⁻¹ for Δ_i (green ellipses in Fig. 3) and the nth octagon agreement within n experimental confidence intervals. The methods tested are the VCI calculations of Tew and Mizukami (T & M),¹⁴ the MCTDH calculations of Richter and Carbonnière (R & C),¹⁵ as well as VPT2 calculations at the B2PLYP-D3(BJ)/aVQZ level. The symbols used are the same as in Fig. 3. The color shade shows whether the experimental observable is overestimated (+), underestimated (−), or met.

With regard to the previous assignment of hot bands of trans-formic acid, the coupling constants to levels with significant thermal population at 190 °C predicted with B3LYP-D3(BJ)/aVTZ (see Table 1) are in good agreement with those at the B2PLYP-D3(BJ)/aVQZ VPT2 level. The largest discrepancy amounts to 0.6 cm⁻¹ (x₃₆), which is below the spectral resolution of the Raman experiment.

3.3 Instabilities of DFT functionals

As previously mentioned in Section 2.2, all production calculations have been carried out without the use of symmetry using the pruned ultra fine integration grid of Gaussian 09.⁶¹ To explore the influence of symmetry and grid size, additional calculations have been performed exploiting the C_s symmetry and a finer integration grid (super fine integration grid, (150 974)).⁶¹ For the following analysis, the five vibrations discussed in this work have been considered for both rotamers, i.e., 10 values.

All density functional theory methods show deviations for anharmonic frequency (VPT2) calculations with and without the use of symmetry when the integration grid size is kept constant, whereas the results obtained with PM3 and MP2 have a negligible (≤0.2 cm⁻¹) dependence on symmetry. The size of the deviation depends largely on the density functional theory method used as well as on the vibration. The most sensitive vibrations of the fundamentals discussed in this work are the O–H stretching (ν₁) and out-of-plane bending vibration (ν₉), while the smallest deviations are observed for the C=O (ν₃) and C–O stretching vibrations (ν₆). For B3LYP-D3(BJ), B2PLYP-D3(BJ), and PBE0-D3(BJ), these deviations are below ±10 cm⁻¹, with mean absolute deviations of 2.5 cm⁻¹, 1.7 cm⁻¹, and 2.1 cm⁻¹ for the ultra fine integration grid, respectively. Particularly severe divergence is observed for ωB97-XD and M06-2X with discrepancies of up to −96 cm⁻¹ and 133 cm⁻¹, respectively. The mean absolute deviations for these methods are as large as 30.2 cm⁻¹ (ωB97-XD) and 59.5 cm⁻¹ (M06-2X). These can be reduced by using the finer integration grid (super fine integration grid). This is illustrated in Fig. 5, where the mean absolute deviation of the band positions using C₁ and C_s symmetry is plotted for both grid sizes (black and blue squares). This decrease in divergence, however, occurs at the expense of distinctly higher computational costs. In case of ωB97-XD and M06-2X, this leads to an mean absolute deviation of 2.8 cm⁻¹ and 48.0 cm⁻¹. The large value for M06-2X is caused by outliers where the deviation between calculations with C_s and C₁ symmetry is enhanced by using the finer grid (ν₉(cF, F) and ν₆(cF)).

Fig. 5 Mean absolute deviations (MAD, in cm⁻¹) of anharmonically (VPT2) calculated band positions of the stretching vibrations (ν₁, ν₂, ν₃, and ν₆) and the O–H out-of-plane bending vibration (ν₉) of cis- and trans-formic acid resulting from the usage of symmetry (C_s) compared to no symmetry (C₁) or the increase of the DFT integration grid size (super fine integration grid compared to the ultra fine integration grid), as implemented in Gaussian 09.⁶¹

When just the integration grid is varied and the symmetry is kept fixed (either C₁ or C_s), the band positions vary on average between 1–2 cm⁻¹ for B3LYP-D3(BJ), B2PLYP-D3(BJ), and PBE0-D3(BJ). This is on the same order of magnitude as the symmetry effects discussed above. Again, a huge impact of the integration grid size is seen for ωB97-XD and M06-2X, where mean absolute deviations of 28.6 cm⁻¹ and up to 63.5 cm⁻¹ are observed (cf. orange and green points in Fig. 5). In both cases, the deviations are larger for the C_s symmetry, whereas for the other methods, it is the other way around.

Altogether, these symmetry and integration grid size dependent variations in anharmonic band positions of the fundamentals of cis- and trans-formic acid are on the order of magnitude of the experimental error bars for B3LYP-D3(BJ), B2PLYP-D3(BJ), and PBE0-D3(BJ). Nonetheless, one should keep in mind that individual outliers are slightly larger. Anharmonic frequency calculations with ωB97-XD and M06-2X on the other hand, show substantial differences with regard to the symmetry and integration grid chosen, so that these results must be viewed with caution, as has been discussed before.^76,77 For most methods, the best agreement with experiment is achieved with the C_s symmetry and the finer integration grid. Since the improvement of the accuracy is below the experimental confidence interval for the more reliable DFT methods, if present at all, Fig. 3 and 4 would only change slightly.

4 Conclusions

Overall, thermal excitation combined with rapid jet quenching and Raman probing as reported in this work provides access to the four stretching vibrations of cis-formic acid in a perturbation-free environment. These reference data points are essential for the validation and comparison of modern quantum chemical methods towards a more global description of this model system. Recent examples are VCI calculations of Tew and Mizukami¹⁴ and MCTDH calculations of Richter and Carbonnière.¹⁵ However, it was also shown that vibrational perturbation theory can be a good compromise between accuracy and computational costs for a reasonably rigid molecule like formic acid, if combined with an adequate method for the electronic structure calculation. In this case, the double hybrid method B2PLYP-D3(BJ)/aVQZ and MP2/aVQZ offer a good compromise between accuracy and cost efficiency, in particular for differences between corresponding cis- and trans-vibrations. A benchmark examining various levels of theory revealed the failure of methods like M06-2X/aVQZ VPT2 or ωB97-XD/aVQZ VPT2 to give consistent results, partly due to numerical grid size and symmetry sensitivity. With the single gas phase value from 2006 (ref. 5) available up to a year ago, these conclusions could not have been drawn. A side effect of the thermal population of cis-formic acid is the significant enhancement of hot bands of trans-formic acid compared to room temperature spectra. The anharmonicity constants that can be deduced from these can help to validate combination band assignments, which are in some cases still under debate.¹⁵ Finally, further experiments such as deuteration or depolarisation experiments will help to shed more light on various debates surrounding trans-formic acid. A prominent example is the assignment of ν₅ and the overtone 2ν₉, where calculations of Tew and Mizukami¹⁴ and Richter and Carbonnière¹⁵ disagree with the experimental, infrared spectroscopic assignments of Freytes and co-workers³⁷ as well as Raman spectra of Bertie and Michaelian.²⁷ Additional Raman data recorded in the fashion shown here, i.e., in combination with thermal excitation, show a distinctly higher intensity for the band previously assigned to the overtone of ν₉ (1305 cm⁻¹ (ref. 37)) compared to ν₅ (1223 cm⁻¹ (ref. 37)), making this a fascinating disagreement of IR and Raman intensity patterns to be resolved.^8,45 Indeed, a very recent IR investigation⁵⁵ points into the same direction.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors would like to thank A. Nejad for very valuable discussions and T. Forsting for the construction of the heatable feed line. This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – 389479699/GRK2455.

Notes and references

1. J. M. Bowman, T. Carrington and H.-D. Meyer, Mol. Phys., 2008, 106, 2145–2182.
2. Z. Bačić and J. C. Light, J. Chem. Phys., 1987, 86, 3065–3077.
3. G. E. Scuseria and H. F. Schaefer, J. Chem. Phys., 1989, 90, 3629–3636.
4. J. Altnöder, K. Krüger, D. Borodin, L. Reuter, D. Rohleder, F. Hecker, R. A. Schulz, X. T. Nguyen, H. Preiß, M. Eckhoff, M. Levien and M. A. Suhm, Chem. Rec., 2014, 14, 1116–1133.
5. O. I. Baskakov, V.-M. Horneman, J. Lohilahti and S. Alanko, J. Mol. Struct., 2006, 795, 49–53.
6. C. W. Bock, M. Trachtman and P. George, J. Mol. Spectrosc., 1980, 80, 131–144.
7. I. Yokoyama, Y. Miwa and K. Machida, J. Phys. Chem., 1991, 95, 9740–9746.
8. J. Demaison, M. Herman and J. Liévin, J. Chem. Phys., 2007, 126, 164305.
9. Y. Scribano and D. M. Benoit, J. Chem. Phys., 2007, 127, 164118.
10. S. Maeda, Y. Watanabe and K. Ohno, J. Chem. Phys., 2008, 128, 144111.
11. L. O. Paulson, J. Kaminský, D. T. Anderson, P. Bouř and J. Kubelka, J. Chem. Theory Comput., 2010, 6, 817–827.
12. V. Barone, M. Biczysko and J. Bloino, Phys. Chem. Chem. Phys., 2014, 16, 1759–1787.
13. M. Tsuge and L. Khriachtchev, J. Phys. Chem. A, 2015, 119, 2628–2635.
14. D. P. Tew and W. Mizukami, J. Phys. Chem. A, 2016, 120, 9815–9828.
15. F. Richter and P. Carbonnière, J. Chem. Phys., 2018, 148, 064303.
16. J. K. Wilmshurst, J. Chem. Phys., 1956, 25, 478–480.
17. R. C. Millikan and K. S. Pitzer, J. Chem. Phys., 1957, 27, 1305–1308.
18. R. G. Lerner, B. P. Dailey and J. P. Friend, J. Chem. Phys., 1957, 26, 680–683.
19. K. Hirota and Y. Nakai, Bull. Chem. Soc. Jpn., 1959, 32, 769–771.
20. B. Zuckerman, J. A. Ball and C. A. Gottlieb, Astrophys. J., 1971, 163, L41.
21. I. C. Hisatsune and J. Heicklen, Can. J. Spectrosc., 1973, 18, 135–142.
22. W. H. Hocking, Z. Naturforsch., A, 1976, 31, 1113–1121.
23. R. L. Redington, J. Mol. Spectrosc., 1977, 65, 171–189.
24. E. Bjarnov and W. H. Hocking, Z. Naturforsch., A, 1978, 33, 610–618.
25. E. Willemot, D. Dangoisse and J. Bellet, J. Mol. Spectrosc., 1978, 73, 96–119.
26. J.-C. Deroche, J. Kauppinen and E. Kyrö, J. Mol. Spectrosc., 1979, 78, 379–394.
27. J. E. Bertie and K. H. Michaelian, J. Chem. Phys., 1982, 76, 886.
28. J. E. Bertie, K. H. Michaelian, H. H. Eysel and D. Hager, J. Chem. Phys., 1986, 85, 4779–4789.
29. L. Nemes, A. R. W. McKellar and J. W. C. Johns, J. Opt. Soc. Am. B, 1987, 4, 1165.
30. Y. Maréchal, J. Chem. Phys., 1987, 87, 6344.
31. W. H. Weber, P. D. Maker, J. Johns and E. Weinberger, J. Mol. Spectrosc., 1987, 121, 243–260.
32. R. E. Bumgarner, J.-I. Choe, S. G. Kukolich and R. J. Butcher, J. Mol. Spectrosc., 1988, 132, 261–276.
33. J. D. Goddard, Y. Yamaguchi and H. F. Schaefer, J. Chem. Phys., 1992, 96, 1158–1166.
34. G. Luiz, A. Scalabrin and D. Pereira, Infrared Phys. Technol., 1997, 38, 45–49.
35. M. Pettersson, J. Lundell, L. Khriachtchev and M. Räsänen, J. Am. Chem. Soc., 1997, 119, 11715–11716.
36. D. Hurtmans, F. Herregodts, M. Herman, J. Liévin, A. Campargue, A. Garnache and A. A. Kachanov, J. Chem. Phys., 2000, 113, 1535–1545.
37. M. Freytes, D. Hurtmans, S. Kassi, J. Liévin, J. Vander Auwera, A. Campargue and M. Herman, Chem. Phys., 2002, 283, 47–61.
38. F. Madeja, P. Markwick, M. Havenith, K. Nauta and R. E. Miller, J. Chem. Phys., 2002, 116, 2870–2878.
39. M. Winnewisser, B. P. Winnewisser, M. Stein, M. Birk, G. Wagner, G. Winnewisser, K. M. Yamada, S. P. Belov and O. I. Baskakov, J. Mol. Spectrosc., 2002, 216, 259–265.
40. A. Perrin, J.-M. Flaud, B. Bakri, J. Demaison, O. Baskakov, S. Sirota, M. Herman and J. Auwera, J. Mol. Spectrosc., 2002, 216, 203–213.
41. M. Pettersson, E. M. S. Maçôas, L. Khriachtchev, J. Lundell, R. Fausto and M. Räsänen, J. Chem. Phys., 2002, 117, 9095–9098.
42. E. M. Maçôas, J. Lundell, M. Pettersson, L. Khriachtchev, R. Fausto and M. Räsänen, J. Mol. Spectrosc., 2003, 219, 70–80.
43. F. Madeja, A. Hecker, S. Ebbinghaus and M. Havenith, Spectrochim. Acta, Part A, 2003, 59, 1773–1782.
44. O. I. Baskakov, E. A. Alekseev, R. A. Motiyenko, J. Lohilahti, V.-M. Horneman, S. Alanko, B. P. Winnewisser, I. R. Medvedev and F. C. de Lucia, J. Mol. Spectrosc., 2006, 240, 188–201.
45. O. I. Baskakov, I. A. Markov, E. A. Alekseev, R. A. Motiyenko, J. Lohilahti, V.-M. Horneman, B. P. Winnewisser, I. R. Medvedev and F. C. de Lucia, J. Mol. Struct., 2006, 795, 54–77.
46. O. I. Baskakov, B. P. Winnewisser, I. R. Medvedev and F. C. de Lucia, J. Mol. Struct., 2006, 795, 42–48.
47. A. Perrin, J. Vander Auwera and Z. Zelinger, J. Quant. Spectrosc. Radiat. Transfer, 2009, 110, 743–755.
48. A. Olbert-Majkut, J. Ahokas, J. Lundell and M. Pettersson, Chem. Phys. Lett., 2009, 468, 176–183.
49. A. Domanskaya, K. Marushkevich, L. Khriachtchev and M. Räsänen, J. Chem. Phys., 2009, 130, 154509.
50. K. Marushkevich, L. Khriachtchev, J. Lundell, A. V. Domanskaya and M. Räsänen, J. Mol. Spectrosc., 2010, 259, 105–110.
51. M. Hakala, K. Marushkevich, L. Khriachtchev, K. Hämäläinen and M. Räsänen, J. Chem. Phys., 2011, 134, 054506.
52. L. O. Paulson, D. T. Anderson, J. Lundell, K. Marushkevich, M. Melavuori and L. Khriachtchev, J. Phys. Chem. A, 2011, 115, 13346–13355.
53. S. Oswald, E. Meyer and M. A. Suhm, J. Phys. Chem. A, 2018, 122, 2933–2946.
54. K. A. E. Meyer and M. A. Suhm, J. Chem. Phys., 2018, 149, 104307.
55. K. Hull, T. Wells, B. E. Billinghurst, H. Bunn and P. L. Raston, AIP Adv., 2019, 9, 015021.
56. S. Cuadrado, J. R. Goicoechea, O. Roncero, A. Aguado, B. Tercero and J. Cernicharo, Astron. Astrophys., 2016, 596, L1.
57. L. Khriachtchev, M. Pettersson and M. Rasanen, J. Am. Chem. Soc., 2002, 124, 10994–10995.
58. S. Ataka, H. Takeuchi and M. Tasumi, J. Mol. Struct., 1984, 113, 147–160.
59. L. Khriachtchev, M. Pettersson and R. B. Gerber, J. Phys. Chem. A, 2015, 119, 2187–2190.
60. T. Forsting, H. C. Gottschalk, B. Hartwig, M. Mons and M. A. Suhm, Phys. Chem. Chem. Phys., 2017, 19, 10727–10737.
61. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery Jr, J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, Ö. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski and D. J. Fox, Gaussian 09 Revision E.01, 2009.
62. A. D. Becke, J. Chem. Phys., 1993, 98, 5648.
63. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 785–789.
64. S. Grimme, J. Chem. Phys., 2006, 124, 034108.
65. C. Møller and M. S. Plesset, Phys. Rev., 1934, 46, 618–622.
66. M. Head-Gordon, J. A. Pople and M. J. Frisch, Chem. Phys. Lett., 1988, 153, 503–506.
67. Y. Zhao and D. G. Truhlar, Theor. Chem. Acc., 2008, 120, 215–241.
68. J.-D. Chai and M. Head-Gordon, Phys. Chem. Chem. Phys., 2008, 10, 6615–6620.
69. C. Adamo and V. Barone, J. Chem. Phys., 1999, 110, 6158–6170.
70. J. J. P. Stewart, J. Comput. Chem., 1989, 10, 209–220.
71. J. J. P. Stewart, J. Comput. Chem., 1989, 10, 221–264.
72. S. Grimme, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2011, 1, 211–228.
73. J. Bloino and V. Barone, J. Chem. Phys., 2012, 136, 124108.
74. N. O. B. Lüttschwager and M. A. Suhm, Soft Matter, 2014, 10, 4885–4901.
75. V. Barone, M. Biczysko, J. Bloino, P. Cimino, E. Penocchio and C. Puzzarini, J. Chem. Theory Comput., 2015, 11, 4342–4363.
76. S. E. Wheeler and K. N. Houk, J. Chem. Theory Comput., 2010, 6, 395–404.
77. N. Mardirossian and M. Head-Gordon, J. Chem. Theory Comput., 2016, 12, 4303–4325.

---