On the vibrations of formic acid predicted from first principles

In this article, we review recent first principles, anharmonic studies on the molecular vibrations of gaseous formic acid in its monomer form. Transitions identified as fundamentals for both cis- and trans form reported in these studies are collected and supported by results from high-resolution experiments. Attention is given to the effect of coordinate coupling on the convergence of the computed vibrational states.


Introduction
Vibrational spectroscopies continue to be a shared interest of both theoretical and experimental approaches in chemistry, with the harmonic approximation 1 as the dominating theoretical description in practice. In contrast, methods that describe vibrations from first principles, obey to the rules of quantum mechanics and go beyond the harmonic approximation are far from routine, despite the apparent anharmonicity of many molecular vibrations. Since a clear hierarchy is not well established for these methods, investigating their performance and characteristics in a joint discussion with regard to the accuracy of the predicted vibrations of a selected molecule may be of interest in addition to the comparison to high-accuracy experimental data. Formic acid, the smallest carboxylic acid, emerges as a good candidate in this respect; it exhibits structural cistrans isomerism and it forms a hydrogen-bonded dimer at room temperature, giving rise to complex nuclear dynamics which are difficult to model accurately. 2 In a broader perspective, understanding the nuclear dynamics of carboxylic acids could also give insight into the conformation and folding mechanisms of proteins. After a brief discussion of selected concepts employed in the anharmonic and variational computation of molecular spectra, we review experimental results and discuss ab initio calculations of the anharmonic molecular vibrations of the monomer forms of formic acid in a chronological order.

Selected concepts of the variational calculation of molecular vibrations
Arguably one of the most difficult aspects of the calculation of molecular vibrations beyond the harmonic approximation is the definition of an appropriate coordinate system and the derivation of the corresponding vibrational Hamiltonian. The separation of the purely vibrational motion of the molecule requires internal coordinates, 1,3 but the conditions which define a unique and an optimal coordinate system which is appropriate for various types of nuclear motion are not known. 4 The Watson Hamiltonian which was derived for non-linear 5 and linear 6 reference configurations and defines the rovibrational motion in terms of rectilinear normal coordinates remains to be one of the most popular choices, since its definition ensures minimal rovibrational coupling and an optimal representation of molecular vibrations with small displacements from the reference geometry. 1,4 It follows, that difficulties are encountered if vibrations corresponding to e.g. torsional motions or significant displacements from the equilibrium structure are computed in this representation. Slow convergence of the computed vibrational eigenstates or inaccurate results in approximate methods are often encountered in these cases. 7 Since the Watson Hamiltonian is invariant under the unitary transformations of the normal modes, attempts for improvement have been made by optimizing [8][9][10] and localizing 11 the normal modes. The choice of the vibrational Hamiltonian also affects the ease of interpretation of the computed results. Whereas vibrations which are well described by the harmonic approximation are found to be mostly localized on the respective normal coordinate, anharmonic vibrations which are encountered e.g. in the high-energy region, can have contributions along multiple normal coordinates, leading to increasingly difficult assignments. 12 Whereas an ill-defined coordinate system can introduce artificial correlation, a well-choosen one may minimize the correlation between the different degrees of freedom. 13 The main focus of this section is not a to give a comprehensive review of available vibrational methods, for which already vast literature is available, 4,14-19 but rather to give an overview of selected concepts and methods necessary for discussing the literature on formic acid, reviewed below. We also shall not discuss the computation of intensities, but emphasize the importance of the development of e.g. accurate electronic dipole surfaces aiding the interpretation of IR spectra.
In order to present a joint discussion of variational methods and the literature reviewed below, we introduce the following notation for the vibrational coordinates q = {q 1 ,. . .,q M }. Once defined, the vibrational Hamiltonian may be separated into a kinetic and potential part, where the form of both the kinetic and potential term depends on the employed coordinate system. The potential term, the so called Potential Energy Surface (PES) V(q) can be calculated either by direct evaluation of ab initio electronic energies along the vibrational coordinates, or it may be pre-fitted to an appropriate analytic form. Most commonly, either a Taylor expansion 20 @VðqÞ @q m @q n @q l q 0 q m q n q l þ . . . ; (1) where q 0 is a reference geometry, or a many-body expansion 21 is employed, where the V m (q m ) are the so-called one-, V mn (q m ,q n ) are the two-, and V mn. . .M (q m ,q n ,. . .,q M ) are the Mbody terms, respectively, each dependent on a subset of coordinates. The sum-of product (SOP) form 22,23 VðqÞ is also commonly employed. Here, the potential is expressed as a sum with expansion coefficients c r of potential terms V r , where the latter is expressed as a product of terms dependent only on coordinate Q k , respectively. The Q k s may combine multiple coordinates of q, therefore reducing the computational effort, i.e. f o M and r = [r 1 ,. . .,r f ]. We note, that the accuracy of any vibrational structure calculation is limited by the accuracy of the PES, and therefore the electronic structure method employed for its construction.
Beyond the well-known harmonic approximation, for which the solutions can be obtained exactly due to approximating the potential as a second order Taylor expansion in normal coordinates, variational solutions 24 to the vibrational Schrödinger equation require an ansatz. In the Vibrational Self-Consistent-Field (VSCF) method, 18,25 the vibrational states are computed in the product form where the functions dependent on a single vibrational coordinate are eigenfunctions of the mean-field operator 26 where we have assumed that the kinetic energy operator can be written as a sum of independent terms, T m (q m ). n = [n 1 ,. . .,n M ] collects the available quantum numbers associated with the single-coordinate functions f nm (q m ) determined self-consistently, and % n denotes all quantum numbers but the m-th. The vibrational energies are then the expectation value of the full Hamiltonian with the VSCF wavefunction (eqn (4)) which can be obtained in a state-specific, for a chosen set of quantum numbers, or in a non state-specific fashion. The VSCF approximation has been improved by the correlation-corrected VSCF approach (CC-VSCF) developed by Gerber and co-workers 18,27,28 also known as second order vibrational Møller-Plesset (VMP2) theory. In this method, the deviation of the exact wave function from the mean-field solution is assumed to be small, such that a second order perturbation correction is applied for a correlation-corrected energy expression. The latter differs in the zeroth order problem from the well-known second order vibrational perturbation (VPT2) theory in which the harmonic approximation is corrected.
Going beyond the product form, a linear combination of configurations with expansion coefficients C n may be used. A full direct product basis can be employed as a basis in eqn (6) (e.g. products of weighted classical orthogonal polynomials, 29 harmonic oscillator basis functions 30 ), however, without further approximations, this leads to an ansatz with an exponential scaling with regard to the size of the system treated and therefore a blowup in the number of terms in the expansion. The truncation of the direct product basis according to some condition is therefore common practice. In Vibrational Configuration Interaction (VCI) 26,31,32 the configurations are usually taken from preceding VSCF calculations. The ansatz can be written as where c i (q) is the reference configuration andt l generates an excited configuration from the quantum numbers i = [i 1 ,. . .,i M ].
If the full excitation space is used, the method is called full VCI (FVCI). Similarly to CI in electronic structure theory, singles, doubles, and higher excitations are generated, and the expansion is truncated to allow for a subset of the available quantum numbers and/or excitations, therefore introducing a bound on the, in

Review PCCP
principle, infinite sum. The expansion coefficients C n of the normalized wavefunction are determined variationally as eigenvectors of the VCI matrix hc n (q)|H(q)|c l (q)i. However, the straightforward construction of the untruncated VCI matrix and computation of its eigenvectors is prohibitively expensive and in practice approximations are introduced. 33,34 By partitioning the modes into a set of active and a set of bath modes Mizukami and Tew 35 have introduced the vibrational active space self-consistent field theory (VASSCF), vibrational active space configuration interaction (VASCI) and vibrational active space second order perturbation theory (VASPT2) methods. These methods employ a product of a CI type ansatz for the active modes and a simple product ansatz for the bath modes. Whereas in their VASSCF approach all CI coefficients and all single-coordinate functions are optimized, only a subset of single-coordinate functions are optimized in their VASCI approach.
In the Multi-configurational time-dependent Hartree 22,36 (MCTDH) method a time-dependent ansatz is employed with expansion coefficients C r where f corresponds to the number of combined coordinates (particles) and the orthonormal single-particle functions f rk (Q k ,t) are determined as linear combination of primitive basis functions with time-dependent expansion coefficients, determined variationally. In practice an upper bound is introduced for the sum by limiting the number of single-particle functions. The eigenstates can be computed via the improved relaxation or block-improved relaxation method of MCTDH. 23 Since MCTDH uses time-dependent configurations, the configurational space can naturally adapt during the evolution of the wave function and can therefore be compared to the vibrational multiconfigurational self-consistent field method. 37 Generally speaking, the underlying coordinate system defines the level of coupling among the coordinates that has to be captured by the Hamiltonian and the wave function of each vibrational state. The convergence of the Hamiltonian with respect to the convergence of the computed energy levels can be tested by e.g. including higher-order terms of a Taylorexpanded PES or higher-order coupling terms in a many-body expanded PES. A systematic improvement of variational wave functions is possible through e.g. the enlargement of the underlying basis set(s) and by employing multiple configurations. Finally, SCF and dynamic methods have the advantage of adapting the single-coordinate functions to a given computed state. Besides obtaining converged states within the respective approach, the accuracy of the latter can be quantified by comparison to high-accuracy experimental frequencies, or a higher level of theory, if available. It follows, that an inaccurate PES computed at a low-level of electronic structure theory or not supportive of the relevant nuclear positions will not yield good results, as will not a single-configurational approach for a coordinate system which is inappropriate for the computed vibration (e.g. rectilinear coordinates for vibrations involving torsional motions). In addition, resonances between states cannot be captured by single-configurational approaches such that any variational or perturbative approach relying only on one configuration will fail in this respect. 35 We note, that although multi-configurational approaches can provide the most accurate results, due to their inherent scaling with system-size, they are difficult to apply to larger systems or highly-coupled motions requiring large bases. In this respect, perturbative approaches represent cost-effective alternatives beyond the harmonic approximation.
Once the vibrational states are determined, the interpretation of the computed levels presents the next challenge. In the harmonic approximation the vibrations are decoupled, such that a single quantum number is sufficient for characterizing a molecular vibration, i.e. n n 1 corresponds to the n-th solution of the harmonic oscillator in the coordinate q 1 . If n = 1 the vibration is called a fundamental and otherwise with n 4 1 an overtone vibration. As mentioned at the beginning of this section, it is advantageous to find a representation which allows each vibration to be localized along its coordinate and therefore to be decoupled from the remaining degrees of freedom. In this case, the product form is already sufficient and an unequivocal assignment in terms of a single set of quantum numbers is possible. Then, the harmonic or VSCF approximation may already present a satisfactory description of the vibrational state, and the corresponding VCI wave function is characterized by a single large CI coefficient. In the opposite case, the absence of a prevalent configuration may be regarded as the evidence of the vibrational counterpart of the concept of static correlation in electronic structure theory. As a dynamic method, in MCTDH the single-particle functions naturally adapt to the computed state, such that in principle a compact wave function can be achieved. In certain cases a singleconfigurational description completely fails, as is the case for resonances: in Fermi resonance, 38 for instance, an energy splitting, ascribed to the interaction of an overtone and a fundamental accidentally close in energy, is observed.

Formic acid
Despite its size, formic acid exhibits a rich structural complexity, which has led to an ongoing theoretical and experimental interest in its nuclear dynamics over the years. 30,35,[39][40][41][42][43][44][45][46][47][48][49][50][51][52][53][54][55][56][57] In its monomer form (Fig. 1), two stable configurations of the C s point group symmetry, cis and trans have been identified via microwave 57-59 and submillimeter spectroscopy 60,61 in the gas phase. With a difference of about 1365 cm À1 in zero point energy, as measured by gas-phase microwave relative intensity measurements, 57 and 1412 cm À1 62 and 1415 cm À1 56 determined ab initio, the trans form is significantly more abundant at room temperature. This relative stability has been assumed to be due to an intramolecular hydrogen bond, 20 but this idea was discarded later on. 20,63 The gas-phase structure of the trans form was determined via electron diffraction. 64 Based on the ground-state rotational constants determined by microwave spectroscopy, approximate equilibrium structures have been derived. 65 The geometries can be connected via an internal rotation of hydrogen of the OH group around the CO axis, 66 where the earliest experiments found a barrier height of 4827 cm À1 . 57 At low temperatures, proton tunneling through the torsional barrier was found to limit the lifetime of cis-HCOOH as measured in Argon matrix by Petterson et al., 42,67 whereas the deuterium tunneling rate was found to be significanty slower. 68 The first IR absorption and Raman spectra of gaseous formic acid were measured in 1938 69 and 1940 70 by Bonner et al., and numerous experimental studies determining the vibrational transitions have been performed ever since at an increasingly high resolution. [39][40][41][42][43][44][45][46][47][50][51][52]54,[71][72][73][74][75][76][77][78][79] The signals arising due to the presence of the cis and trans monomers or the dimer have been differentiated experimentally by signal intensification via population control with thermal 54 or laser 42 excitation. Most variational theoretical methods consider purely the vibrational motion at zero temperature which can be compared to low temperature or high-resolution (e.g. FT absorption) rotationally analysed data. Two popular low temperature spectroscopic methods have emerged in this respect. In matrix isolation spectroscopy 80,81 the analyte is condensed with an inert gas, e.g. a noble gas or nitrogen at low temperature, such that the molecules are ''trapped'' with different local environments in a matrix. This results in a shift of the measured frequencies and a splitting of absorption bands. For highresolution data, the best agreement with theory is therefore achieved if the matrix-analyte interactions are also modeled in the theoretical description. As opposed to matrix isolation spectroscopy, jet-cooled spectroscopy 76 does not suffer from such effects. In this method the analyte is passed through a nozzle with a neutral carrier gas and expanded into an evacuated jet chamber at low temperatures. Population enhancement has been achieved in a jet-setup via thermal control therefore recording and spectra at different temperatures have been compared. 54 Complementary high-resolution IR and Raman spectra of partially and fully deuterated formic acid have been recorded for all isotopologues of the cis and trans monomer (DCOOH, HCOOD, DCOOD). 41,44,45,50,[82][83][84][85][86][87][88] Some spectral regions of formic acid are dense and therefore require a careful rotational analysis, hot bands and overlapping dimer signals further complicate the spectrum with the density of states generally increasing in the higher energy region. 54,82 Amongst the many experimental studies performed for trans HCOOH, 50,55,61,83,88,89 we mention the work of Freytes et al. 82 , who measured gas-phase room-temperature Fourier-transform (FT) and intracavity laser absorption spectra of trans-HCOOH. They reported band origins from 626 cm À1 to 13 284 cm À1 and performed a detailed analysis on the rotational structure of the first CH stretching vibrational overtone n 2 2 and the second OH stretch overtone n 3 1 of formic acid. Their work built on a previous study of Hurtmans et al., who presented a rovibrational analysis in the range of the first and third overtone of the OH stretching vibration n 2 1 and n 4 1 of trans-HCOOH based on FT spectroscopy and intracavity laser absorption spectroscopy measurements. 83 Room temperature Raman measurements were performed by Bertie et al. 40 in the range from 70 cm À1 to 4000 cm À1 for all isotopologues of the dimer, where the monomer signals could also be identified. Recent work of Nejad, Suhm and Meyer 54 provided gas-phase Raman-jet measurements for all four deuterated isotopologues of formic acid monomer of both the cis and trans configuration, where the spectra were recorded with a nozzle temperature of 160 1C. They detected Fermi resonance doublets for the fundamentals n 1 and n 5 of trans-HCOOH, n 3 and n 5 for trans-DCOOH, n 2 and n 5 for trans-HCOOD, and n 2 and n 3 for trans-DCOOD. A summary of available benchmark quality experimental values for the vibrational transitions of trans-HCOOH and deuterated forms up to 4000 cm À1 in the far and mid-IR region with additional Raman-jet measurements has been recently compiled by Nejad and Sibert 70 from gas-phase IR and Raman spectra reported in literature.
Vibrational data for cis-HCOOH and isotopologues is more scarce 48,66,68,[90][91][92][93][94][95] due to the experimental difficulty of populating the cis state. Matrix-isolation spectroscopy measurements in Argon have been performed for cis-HCOOH by narrowband IR pumping of the first OH stretching overtone of trans-HCOOH n 2 1 by Pettersson et al. 42 Maçôas et al. recorded near-and mid-IR spectra for HCOOH, DCOOH and HCOOD in solid argon matrix at 8 K from 400 cm À1 to 7800 cm À1 48 where the cis population was enhanced by pumping a trans to cis transition by narrowband tunable IR radiation. Raman-jet spectra for cis-HCOOH and isotopologues have been recorded recently by Meyer, Nejad and Suhm 54,87,96 via thermal excitation between 100-190 1C prior to the jet expansion. In Table 1 we summarize transition frequencies of gaseous formic acid monomer from highresolution experiments, which have been assigned as fundamentals of the cis and trans monomer, respectively. We note, that many of the fundamentals have been found to be affected by rovibrational 97 and resonance interactions, 20 of which the n 5 /n 2 9 Fermi resonance pair of trans-HCOOH has been the most prominent due to an early experimental misassignment of the more intense overtone band as the fundamental. 82 Throughout this work, the assignment of the fundamentals in terms of nuclear displacements is discussed with the valence coordinate definition of Richter and Carbonniére 56 (Fig. 2), although some of the studies reviewed below employ different coordinates. 4 Ab initio studies on the formic acid monomer

First ab initio studies
One of the first ab initio, anharmonic studies on the formic acid monomer was concerned with the overtone series of the OH stretching vibration n n 1 of trans-HCOOH. 83 Hurtmans et al. recorded rotationally resolved overtone spectra of formic acid with high-resolution gas-phase FT and intracavity laser absorption spectroscopy and performed rovibrational analysis for the first n 2 1 and third n 4 1 overtones in the same work. They constructed an effective one-dimensional Hamiltonian dependent on the normal coordinate of the OH stretching mode for the analysis of the OH stretching vibration by relaxing the remaining geometrical parameters along the coordinate. The resulting effective potential was constructed at the MP2/cc-pVTZ level, and the corresponding one-dimensional variational solutions were computed for n n 1 , n r 4 with an agreement within 50 cm À1 compared to the experimental band origins of the same work. It was found that the excitation of these n n 1 vibrations affects the HOC angle with increasing n. The induced electric dipole moment was found to be initially almost parallel to the OH bond and progressively tilted towards the CQO bond with increasing n. Based on these results Hurtmans et al. proposed a possible proton exchange mechanism promoted by the excitation of the OH overtones.
Subsequently Maçôas et al. 48 performed CC-VSCF calculations to support the assignment of the mid-and near-IR transitions measured by matrix-isolation spectroscopy in the same work. The transitions were determined at two trapping sites in solid argon for cis-and trans-HCOOH. Mid-IR spectra of trans and cis-DCOOH and near-IR spectra of trans-DCOOH and trans-HCOOD were also recorded, where the cis transitions were determined by narrowband IR pumping of the first OH overtone transition n 2 1 of the trans configuration. For cis-HCOOH about half of the expected fundamentals could be recorded and strong matrix effects, such as splitting of the recorded signals, was observed, making the assignments difficult. The CC-VSCF calculations were performed in the normal mode representation without rovibrational coupling terms and the PES was expressed with pairwise coupling terms, constructed directly at the MP2/ 6-311++G(2d,2p) level of theory for each pair of normal coordinate gridpoints. They computed all fundamentals and selected overtone and combination bands in the mid-and near-IR region for both trans-and cis-HCOOH and -DCOOH. We have included in Tables 2 and 3 the computed transitions identified as fundamentals of cis-and trans-HCOOH reported in their work. Anharmonicity was found to be strongest for the highfrequency stretching vibrations of trans and cis-HCOOH, i.e. n 1 and n 2 were found to have a downshift of up to 200 cm À1 compared to their calculated harmonic counterpart, whereas most fundamentals were found to reproduce the matrixisolation experiments reasonably well. Exceptions are fundamentals corresponding to torsional and deformation type motions, that is, n 5 and n 9 were predicted to be below the experimental values, e.g. for n 9 around 40 cm À1 and 70 cm À1 for trans-HCOOH and cis-HCOOH, respectively. In the near-IR region, the first overtones n 2 1 and n 2 2 of trans-HCOOH and cis-HCOOH were found to be reasonably accurate compared to the measured experimental values. On the other hand, most combination and overtone bands involving n 5 and n 6 were inaccurate. In addition to limitations based on e.g. the level of ab initio theory and the order of coupling in the potential, the reason for this was found to be the normal modes associated with the aforementioned fundamentals. We note, that this inaccuracy may be due to the inability of VSCF to reproduce Fermi resonances. The well-known n 5 /n 2 9 Fermi resonance was observed experimentally for trans-HCOOH and no analogous resonance was identified for cis-HCOOH. For DCOOH an n 3 /n 2 8 Fermi resonance was observed for both the trans and cis form with the involvement of the OH stretching vibration n 1 .
We also briefly mention the work of Scribano and Benoit, 101 who four years later computed the OH stretching frequency n 1 of trans-HCOOH with CC-VSCF, VCI and single-to-all (STA) CC-VSCF and STA-VCI. The STA methods were introduced in their work and rely on the idea that the modes could be partitioned into an active and a ''spectator'' set. In these calculations, the rovibrational coupling terms were neglected and the potential was expressed with up to two-body coupling terms in normal coordinates. The direct evaluation of the potential was performed at the MP2 and CCSD(T) level with TZP and polarization augmented SBK bases. The OH stretch fundamental n 1 was computed using standard VCI and STA-VCI with only one active degree of freedom. At the CCSD(T)/TZP level of theory a the STA-CC-VSCF and STA-VCI energies were found to be only 8 cm À1 worse compared to the CC-VSCF and VCI calculations for the frequency of the n 1 fundamental, where only singles and doubles excitations were included in the VCI expansion.

Advanced calculations with analytic potential energy surfaces
A more accurate description of the nuclear dynamics of formic acid has become available due to the development of analytic ab initio PESs, the availability of which leads to a significant speedup in the potential energy evaluation and therefore allows for the use of higher accuracy vibrational structure methods.  The first anharmonic ab initio force-fields of formic acid were developed for trans-HCOOH by Demaison et al. 20 at the MP2/ VTZ, MP2/aVTZ and CCSD(T)/VTZ(ae) level of theory. In addition to analytic harmonic terms, the cubic and semidiagonal quartic force constants were determined via displacements along the normal coordinates with the method of finite differences. With this force field they computed anharmonic fundamentals using vibrational perturbation theory for trans-HCOOH. In Table 2 we included transitions obtained with the CCSD(T)/ VTZ(ae) force-field. Vibration-rotation interaction constants as obtained by standard perturbation theory were also reported. The n 7 , n 9 , and n 3 fundamentals of trans-HCOOH were found to be affected by strong rovibrational interactions. They carried out an analysis of the resonances in terms of the magnitude of the force constants, where the Fermi resonant states were assumed to be connected via the cubic constants of force constants of type k ijj and Darling-Dennison type resonances 102 via the quartic constant k kkll . In contrast to the experimentally derived force constants, several calculated force constants indicating Fermi resonance (e.g. k 499 for n 4 /n 2 9 , k 677 for n 6 /n 2 7 ) differed significantly in magnitude and sign whereas the n 5 /n 2 9 Fermi resonance matched qualitatively. The resonance interactions were further analyzed by constructing a 5 Â 5 symmetric matrix representing interactions involving the n 2 9 , n 2 7 , n 6 , n 5 , n 4 states and computing the corresponding eigenstates.
In 2013, Tew and Mizukami reported their first PES for formic acid, initially fitted for the description of the trans form ref. 35 with 5873 randomly generated configurations within the energy range 0-15 000 cm À1 at the CCSD(T)(F12*)/ccpVDZ-F12 level of theory. With this preliminary surface, they performed VSCF, VMP2, VASCI, and VASPT2 calculations based on the normal-coordinate Watson Hamiltonian for the fundamentals of trans-HCOOH (see Table 2). Rovibrational coupling terms were neglected and up to three body coupling terms were included in the many-body expansion. Up to seven singlemode functions were employed for the state-specific VSCF calculations and for the VASCI calculations an active space of four modes was selected, involving the normal modes associated with the n 1 , n 5 , n 7 , and n 9 fundamentals. The authors note a good performance of VMP2 except for the fundamentals n 1 and n 5 , where the former was overestimated by 22 cm À1 compared to the rovibrational-interaction-corrected experimental value. The largest deviations of the VASPT2 calculations were found for the n 2 and n 4 fundamentals. Since the isomerisation barrier was reached already at six quanta of excitation of the singlecoordinate function along the torsional mode associated with the n 9 fundamental, only a limited excitation space was available. This led to multiple unconverged states due to the involvement of the torsional mode in multiple states. The authors concluded, that for fully converged values adopting a different coordinate system that can treat both the cis and trans minima would be necessary.
Consequently, two semiglobal PESs capable of describing the cis-trans isomerization have been developed. Tew and Mizukami published a PES in 2016 62 constructed using LASSO-based regression 62 with electronic energies obtained at the CCSD(T)(F12*)/cc-pVTZ-F12 level of theory using 17 076 random structures in the energy range 0-15 000 cm À1 relative to the trans minimum. The structures were generated by random displacements from cis and trans equilibria as well as from the transition state of the OH torsional motion and with points generated in the internal coordinate path (ICP) coordinate system which was constructed for formic acid. The latter corresponds to curvilinear path for the torsional motion and 8 normal coordinate displacements. The analytic potential fit consists of a zeroth order surface, which is a sum of Morse functions of atomatom separations, and a correction surface, which is a sum of distributed multivariate Gaussian functions which are a combination of atom-atom separations. In the same work, the authors constructed an internal-coordinate path Hamiltonian (ICPH) using their PES, where a curvilinear path was used for the torsional motion connecting the trans and cis rotamers and the remaining mutually orthogonal rectilinear modes were defined such that the rovibrational coupling constant associated with the motion along the path was minimized. Trigonometric functions were used for the expansion of the single-mode functions along the torsional path coordinate, and harmonic oscillator functions were used for the remaining coordinates. State-specific VSCF and VMP2 calculations were carried out with the ICP Hamiltonian of HCOOH, where the potential included up to four-body coupling terms, whereas in the ICPH-VCI calculations up to five-body coupling terms were included with a maximum overall excitation of up to 10 quanta in the VCI expansion. The change of the zero point energy by including, instead of four-body terms, five-body terms in the Hamiltonian was 2 cm À1 and for the fundamentals n 2 , n 3 , n 5 , n 8 n 9 it was within 1 cm À1 . For the n 1 , n 4 , n 6 , and n 7 fundamentals slower convergence with a change of up to 6 cm À1 was observed with respect to the increase in the coupling. The ICPH-VCI vibrational energies with respect to the trans-minimum were reported up to 4750 cm À1 with respect to the global trans minimum and 14 bands in the region 0-4720 cm À1 were reassigned based on the VCI coefficients compared to the assignment of Freytes et al. 82 The cis-zero point energy (ZPE) and cis states were identified from the quantum number of the torsional mode. The authors note the overestimation of the n 2 6 and n 3 6 states in their ICPH-VCI calculations compared to experiment. They also discuss the ambiguity of the assignment of some states in the high-energy region and of increasing excitation, further complicated by the high density of states in this region. They observed a coupling between the n 4 , n 5 , n 6 , and n 7 cis states and concluded that the coupling could be an artifact from the definition of the coordinate system. Due to this coupling and the associated slow convergence only states with up to two quanta of excitations relative to the cis-ZPE were reported. The highest cis-HCOOH fundamental n 1 was not reported in this work, and the n 2 4 overtone was identified to be part of a Fermi doublet with n 2 . We report only the transitions identified as fundamentals in their work in Tables 2 and 3.
The second semiglobal ab initio PES was developed by Richter and Carbonnière 56 in 2018 who fit 660 single point calculations at the CCSD(T)-F12a/aug-cc-pVTZ level of theory in the energy range 0-6000 cm À1 relative to the trans minimum. The PES was fit in terms of internal valence coordinates 1 (i.e. bond stretches, angles and torsions) with the points generated by the Adaptive Generation of Adiabatic PES (AGAPES) procedure. 103 The coordinate definition used is depicted in Fig. 2. The barrier between the two minima along the torsional coordinate t 2 was found to be 4442 cm À1 , where the barrier height was found to be affected by the bond-length r 1 and the bending angle y 2 . The barrier of this PES is 20 cm À1 higher compared to the value reported by Tew and Mizukami. 56 For the vibrational Hamiltonian, the exact analytical kinetic energy operator was constructed in the SOP form with the TANA program 104 in polyspherical coordinates, and the PES was expanded as a many-body expansion in valence coordinate displacements. The ''dead branching'' 105 option of the AGAPES procedure was used, in which only certain coordinate combinations were selected based on a remoteness measure defined by geometric mean values for the potential coupling terms. Out of the 256 possible coordinate combinations in the four-body expansion, 36, 84 and 126 were kept, for two, three and fourbody terms, respectively. A maximal number of terms was kept for the coupling between the y 3 and t 2 coordinates, associated with the n 5 and n 9 fundamentals. Below 3500 cm À1 with respect to the trans ZPE only half the period of the grid of the torsional coordinate was used (t 2 r p/2), such that no cis-trans delocalization effects were included in the computation of the trans states. Full-period calculations were performed between 3500 cm À1 -3640 cm À1 . The high-lying fundamentals associated with the CH and OH stretch vibrations (n 2 , n 1 ) were computed with the improved relaxation scheme of MCTDH, and all states below 3100 cm À1 for trans and 2700 cm À1 for cis with respect to their respective ground state were computed with the block Davidson scheme of MCTDH. We have included the transitions identified as fundamentals in Tables 2 and 3. The MCTDH reference state was found to have only a minor effect on the convergence of the states. All vibrational states below 3796 cm À1 were converged in either the cis or the trans well, as differentiated by the torsional coordinate, and the first delocalized state was found at 3796 cm À1 above trans ZPE, which was labeled as the n 6 9 overtone of cis-HCOOH. Four states could not be converged to sufficient accuracy, these were the states at (assignments of Richter Carbonnière in brackets) 3086 cm À1 (n 4 + n 3 ), 3116 cm À1 (n 2 8 + n 2 6 ) for trans-HCOOH and at 2357 cm À1 (n 9 + n 7 + n 5 ), 2479 cm À1 (n 2 7 + n 5 ) for cis-HCOOH. For most states which were previously identified in literature to be involved in Fermi resonances, a single-configurational description was found to be sufficient, except for the well known n 5 /n 2 9 Fermi resonance, where a mixing of the reference configurations was observed. The computed states of trans-HCOOH were in good agreement with the values reported by Tew and Mizukami with a root mean square deviation (RMSD) of 11 cm À1 : the RMSD of the fundamentals, two, three and four quanta excited states were 5 cm À1 , 8 cm À1 , 9 cm À1 and 25 cm À1 , respectively with maximal deviations (MAXDs) for n 3 , n 3 2 , n 5 + n 2 7 and n 2 7 + n 2 9 . The RMSDs of the cis states were about 49 cm À1 compared to the levels of Tew and Mizukami, where the RMSDs for fundamentals, two and three quanta excited states were 14 cm À1 , 49 cm À1 and 111 cm À1 , respectively, and the MAXDs were found for the n 8 , n 4 + n 8 and n 8 + n 2 9 states. This discrepancy of the cis results compared to the work of Tew and Mizukami was further investigated by the authors. 56 For cis-HCOOH the authors validated their results by performing the AGAPES construction referenced on the cis, instead of the trans geometry: the levels computed with the global PES were found to be consistent with the latter. VPT2 and MCTDH calculations were performed with both the present PES and the 2016 PES of Tew and Mizukami, for which the latter was re-fitted in the SOP form. The VPT2 fundamentals were found to compare well with each other and with the original MCTDH values such that the authors concluded that the disagreement for the cis levels stems from the ICPH-VCI calculations, rather than from the PES. For the cis levels, a mixing of the n 5 and n 6 levels was identified.

Recent developments
Fundamentals, combination bands and overtone transitions up to approximately 3000 cm À1 with respect to the global trans minimum were computed by Aerts et al. 53 in 2020 for the deuterated forms of cis and trans formic acid using the above mentioned PES of Richter and Carbonnière 56 with the blockimproved relaxation method of MCTDH. The assignment of the states was performed via visualization of the reduced densities of the computed vibrational wavefunctions. The RMSDs with respect to experimental data of the fundamental transitions were found to be 8 cm À1 , 7 cm À1 and 3 cm À1 for trans-DCOOD, trans-HCOOD and trans-DCOOH, respectively. The calculations were performed by defining modes combining r 1 and y 3 , y 1 and y 2 , and t 1 and t 2 , respectively (see Fig. 2). Aerts et al. Vibrational transition energies were computed by Nejad and Sibert 79 up to 4000 cm À1 for HCOOH and isotopologues with sixth order Canonical Van Vleck Perturbation theory (CVPT6) employing both surfaces of Tew and Mizukami and Richter and Carbonnière. 56 In their work, the PESs were represented as a Taylor expansion in terms of stretch, bend, and dihedral angles with the valence coordinate definition of Richter and Carbonnière. 56 Up to four-body coupling terms were used for the PES, and the elements of the exact kinetic energy operator were expanded up to the sixth order in the internal coordinates. The CVPT6 calculations were performed with the Hamiltonian constructed in curvilinear normal coordinates based on only one geometry. As a basis for the matrix representations, a product of harmonic oscillator functions was used where the excitation degree of each oscillator was limited. A slightly faster convergence of the trans-HCOOH states obtained with the 2016 PES of Mizukami and Tew as compared to the 2018 PES of Richter and Carbonnière was observed with regard to the order of perturbation, where fourth and sixth order corrections were found to be especially important for higher lying vibrational transitions. Slower convergence was observed for the cis-HCOOH states on the 2016 compared to the 2018 PES. It was found, that the 2016 surface generally overestimates, whereas the 2018 PES often underestimates the transition energies for all isotopologues as compared to experimental values, where increasing deviations were found for higher lying vibrational states. The CVPT6 fundamentals of trans-HCOOH compared within 2 cm À1 to ICPH-VCI and MCTDH results of the respective surface (see Table 2). For cis-HCOOH it was suggested by the authors that the ICPH-VCI results of Tew and Mizukami were not fully converged. The band assignment was performed based on the leading coefficient in the VCI expansion and in comparison to previous calculations and previously assigned IR and Raman bands reported in literature. 11 new vibrational band centers were assigned for trans-HCOOH and 53 for the deuterated isotopologues by the authors with multiple reassignments in comparison to previous assignments from theory or experiment. For trans formic acid strong coupling among the states which share the polyad quantum number N p = n 5 + n 9 /2 was found, with an involvement of the n 1 fundamental at higher energies.
Most recently, Martín Santa Daría, Avila, and Mátyus performed calculations using the 2016 PES of Tew and Mizukami with a Hamiltonian constructed using the cis-trans torsional coordinate and eight curvilinear normal coordinates defined with respect to an instantaneous reference configuration along the torsional motion. 30 The curvilinear normal coordinates were adapted to the cis-trans isomerization of formic acid by relaxing the internal coordinates such that the potential energy was at minimum along the torsional motion. In principle, the kinetic and potential energy coupling among the coordinates was reduced by employing curvilinear normal coordinates. The kinetic energy operator was constructed numerically with the kinetic and potential energy integrals evaluated over the direct product basis with the GENIUSH (general rovibrational code with numerical, internal-coordinate, user-specified Hamiltonians) code. 106 For the curvilinear normal coordinates harmonic oscillator basis functions and for the torsional degree of freedom a Fourier basis was used. A Smolyak grid was then used for the truncated direct product basis for the evaluation of the multidimensional integrals. The authors computed converged vibrational states up to ca. 4700 cm À1 with respect to the trans-ZPE, slightly above the cis-trans isomerization barrier and therefore the highest-energy cis n 1 fundamental was not computed. The convergence was evaluated by increasing the basis set size. The energy difference between calculations using the same basis set size was performed with rectilinear and relaxed curvilinear normal coordinates and the difference was found to be around 8-10 cm À1 , with an increasing difference in the higher energy region. The assignment of the cis-states was performed in terms of the contribution of the 1D torsional basis function, which in the lower energy ranges was localized in either well. A variational improvement of up to 5-10 cm À1 was reported compared to the CVPT6 results and by 10-40 cm À1 compared to the previously computed ICPH-VCI levels. Beyond 3700 cm À1 , mixed cis-trans states were reported and the ''tunneling splittings'' associated with these mixed states were computed to be below 1-5 cm À1 , almost within the convergence uncertainty of the calculations.
Beyond 3900 cm À1 , non-negligible contributions from delocalized torsional states were observed. We note, that although the results of the authors unequivocally indicate the switching of the n 5 and n 6 cis-HCOOH states as opposed to previous calculations, in Table 3 we list the labeling to be consistent with previous results.
Finally, we discuss the recent work of Aerts et al., who simulated the intramolecular vibrational redistribution (IVR) dynamics of formic acid on the 2018 PES of Richter and Carbonniére with the goal of identifying a possible laserinduced trans-cis isomerization pathway. The dynamical evolution following the excitation of an infrared-active vibration was simulated by the time-evolution of an initial vibrational wave function which was prepared by the excitation of the single particle function of the torsional mode t 2 associated with the isomerization. The probability of finding the molecule in the trans well, as defined by the value of the torsional coordinate, was followed during the evolution, and local energies were ascribed to each mode. By investigating the time-evolution of the fractional energies defined for each overtone of the local torsional vibration of trans-HCOOH, multiple mode-couplings and resonances were identified and confirmed. The n 5 /n 2 9 Fermi resonance was observed as a strong, reversible energy flow between the two fundamentals, also for states sharing the previously mentioned polyad quantum number N p = n 5 + n 9 /2. A coupling between the modes of the n 9 and n 4 fundamentals was identified. The probability to find the wavepacket in the initial trans well was found to diminish at n 6 9 at which part of the wavepacket transfered from trans to cis within 100 fs. The non-periodic probability indicated, that this transition between trans and cis is irreversible.

Summary
We have reviewed first principles, anharmonic studies on the molecular vibrations of formic acid in its monomer form and collected all reported computed transitions identified as fundamentals for both cis-and trans-HCOOH. Formic acid has been a widely explored albeit challenging system to both experiment and theory. On one hand, the spectrum is dense and complex in certain spectral ranges and requires a careful rovibrational analysis 83 and experimental differentiation of the signals arising due to its possible geometries. 107 Undistorted, high-resolution experimental data has still not been recorded for all of the available cis transitions predicted by theory. On the other hand, the higherenergy vibrations of the molecule are expected to explore trans and cis geometries or both, 30 a formidable task for vibrational methods. The minimal energy path connecting the two minima was found to be a function of the t 2 torsional coordinate, 56,57 but possible isomerization pathways involving the excitation of the n n 1 OH overtones 83 or involving n 6 and n 9 were also suggested. 2 Early on, multiple resonances affecting the vibrational levels of formic acid were identified, most notably the n 5 /n 2 9 Fermi resonance of trans-HCOOH. 48 Hurtmans et al. reported the involvment of the HOC angle in the n n 1 excitations, 83 however, the first ab initio studies employed rectilinear normal coordinates which were found to be ill-suited for torsional and deformation vibrations. 48 The first PESs were computed directly along the normal coordinates 48 or as a Taylor 20 or a local many-body expansion, 35 such that vibrations with a local trans character could be predicted. The development of analytic semiglobal PESs capable of reproducing transitions both from the cis and trans minima followed. 56,62 Transitions of cis-and trans-HCOOH were computed with a Hamiltonian constructed with normal coordinates along the curvilinear isomerisation path, and the computed high-energy cis states were found to be coupled significantly. 62 The vibrational transitions were also computed with the dynamical MCTDH method in valence coordinates, 56 and the computed cis-states differed significantly from the previously reported results. A comparison between the PESs was made with transitions computed based on a curvilinear normal coordinate Hamiltonian, where the coordinates were defined based on a single geometry. 79 The transitions of the deuterated species 53,79 and assignments up to 4000 cm À1 with respect to trans-ZPE were also described. 79 A curvilinear normal coordinate Hamiltonian adapted to the cis-trans torsional motion was constructed, 30 and vibrational states slightly above the isomerisation barrier were reported including cis-trans entangled and delocalized states. Small splittings associated with the high torsional isomerisation barrier were reported. The study of Martín Santa Daría et al. reported the lowest energy variational solutions up to now 30 with a good agreement with high-accuracy experimental values for both cis (Table 3) and trans ( Table 2) fundamentals.
Finally, we discuss the computed vibrations of formic acid which have been identified as fundamentals across the studies mentioned in this work, as compiled in Tables 2 and 3. Although such an analysis is not necessarily warranted considering that different PES and Hamiltonians are used across these approaches, it is still worth analyzing the results in terms of the employed methods. The harmonic analysis performed by many of the studies shows the worst agreement to experimental values with an RMSD ranging from approximately 80-125 cm À1 where generally a better agreement is achieved for the lower energy states. Single-configurational approaches follow with an RMSD ranging from approximately 3.4-31 cm À1 . Finally, all multi-configurational methods achieve an RMSD better than 4 cm for the trans-states and better than 8 cm À1 for the cis states, where for the latter conclusions should only be drawn with care due to the few available experimental results.
As a small molecule exhibiting large-amplitude motion, the vibrations of formic acid are comparatively difficult to model, where the computation of converged cis states and highly excited states above the delocalization barrier remains to be challenging due to the strong coupling of the coordinates in this region. We note that these difficulties are illustrative of the challenges faced for larger systems, since the nuclear dynamics of bigger molecules and clusters can become increasingly intricate but multiconfigurational methods become prohibitively expensive due to their scaling with respect to system size.

Conflicts of interest
There are no conflicts to declare.