{'fraction_non_alphanumeric': 0.03848219080007159, 'fraction_numerical': 0.01449794165025953, 'mean_word_length': 4.935552407932011, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fingerprint region of the formic acid dimer: variational vibrational computations in curvilinear coordinates

Curvilinear kinetic energy models are developed for variational nuclear motion computations including the inter- and the low-frequency intra-molecular degrees of freedom of the formic acid dimer. The coupling of the inter- and intra-molecular modes is studied by solving the vibrational Schr\"odinger equation for a series of vibrational models, from two up to ten active vibrational degrees of freedom by selecting various combinations of active modes and constrained coordinate values. Vibrational states, nodal assignment, and infrared vibrational intensity information is computed using the the full-dimensional potential energy surface (PES) and electric dipole moment surface developed by Qu and Bowman [Phys. Chem. Chem. Phys. 18, 24835 (2016); J. Chem. Phys. 148, 241713 (2018)]. For further progress in comparison with jet-cooled vibrational spectroscopy experiments, improvement of the potential energy surface is necessary. As to the fingerprint region, the multi-dimensional couplings representation of the PES requires further improvement.


I. INTRODUCTION
In a recent article, Nejad and Suhm reviewed the spectroscopy of the formic acid dimer (FAD) [1] with a focus on the intermolecular vibrational range. The formic acid is the simplest carboxylic acid and its dimer is a prototype for a cyclically arranged pair of hydrogen bonds. The formic acid dimer has been studied in spectroscopy experiments for decades, but the rotational and temperature effects made the detection of the precise vibrational band positions challenging [1].
Coupling infrared and Raman spectrometers to jet-cooled helium beams seeded with a small amount of molecules in the gas phase made it possible to gain precise information on the vibrational dynamics of molecular complexes and clusters [2][3][4][5].
In particular, all intermolecular vibrational fundamentals and several overtone and combination bands of FAD [6][7][8][9][10] have been assigned over the past 16 years with an experimental uncertainty on the order of 1 cm −1 . The amount and the quality of the experimental data call for detailed and high-level quantum dynamics computations [1].
A detailed quantum dynamics computation rests on two, or in fact, three pillars: a) electronic structure methodology that provides good approximate solutions to the electronic Schrödinger equation at a series of nuclear configurations over the coordinate range relevant for the nuclear motion; b) high-dimensional fitting or interpolation methods that build a potential energy surface (function) from the electronic energies available at points; and c) rovibrational methodology that provides solution to the (ro)vibrational Schrödinger equation. A full-dimensional ab initio potential energy surface for the formic acid dimer is already available in the literature (pillars a & b) [11].
In present work, we aim to contribute with a multi-dimensional, internal coordinate, (quasi-)variational (ro)vibrational framework (pillar c) and present computational results specifically targeting the intermolecular range. For this work, we take the already available potential energy surface computed and fitted by Qu and Bowman [11], henceforth labelled with QB16-PES. The QB16-PES was obtained as a least-squares fit of a permutationally invariant potential energy function for the 10-atomic formic acid dimer to 13475 CCCD(T)-F12a/haTZ electronic energies. Qu and Bowman used a maximum of fourth order polynomials in the fitting and they report an 11 cm −1 'energy-weighted' root-mean-squared deviation and an absolute error of about 14 cm −1 for their points below 4400 cm −1 . The global minimum structure of the PES has C 2h point-group symmetry and the corresponding harmonic frequencies are in reasonable agreement with earlier theoretical and experimental work.
Qu and Bowman [12] performed full(24)-dimensional vibrational configuration interaction (VCI) computations with the MULTIMODE computer program using the normal-coordinate representation of the kinetic energy operator (KEO) and a 4mode-representation of QB16-PES. They write about the low-frequency vibrational energies (below 1000 cm −1 ) obtained in the VCI computations that the energies are 'slightly' up-shifted most likely due to the use of rectilinear normal coordinates, which are usually not well suited for describing floppy degrees of freedom, and possible deficiencies in the fitted PES. Otherwise, they estimate the VCI energies to be converged within about 10 cm −1 [12].
In the next section (Sec. II), we define curvilinear internal coordinate vibrational models for FAD with focusing on the intermolecular dynamics. Then (Sec. III), the computed vibrational states are analysed and compared with the experimental data. Section IV is about an assessment of the kinetic and the potential energy representations. The article ends (Sec. V) with a summary of the results, conclusions, and outlook for possible future work.

II. THEORETICAL AND COMPUTATIONAL DETAILS
In the present work, the quantum dynamical computations were carried out using the GENIUSH [13,14] computer program. The general rovibrational Hamiltonian implemented [15][16][17][18][19][20] in this program iŝ whereĴ a (a = 1(x), 2(y), 3(z)) are the body-fixed total angular momentum operators andp k = −i∂/∂q k with the q k (k = 1, 2, . . . , D) internal coordinates. The G kl = (g −1 ) kl coefficients andg = det(g) are determined from the g matrix, defined as follows, and r i are the body-fixed Cartesian coordinates for the i-th atom and e a represent the body-fixed unit vectors.
D ≤ 3N − 6 is the number of the active vibrational dimensions in the system. If D < 3N − 6, then this definition of the kinetic energy operator (KEO) corresponds to imposing rigorous geometrical constraints for the fixed part of the system, and the results depend only on the constrained geometry, but they are independent on the actual coordinate representation of the constrained moiety [13]. This procedure is sometimes referred to as 'reduction in the g matrix' that can be contrasted with 'the reduction in the G matrix' (that is also common and) for which the results would depend on the actual coordinate representation of the constrained fragments.
Regarding the present system, the formic acid dimer has N = 10 atoms and 3N −6 = 24 vibrational degrees of freedom. A fully coupled, variational computation is currently out of reach for such a high number of degrees of freedom, except for highly efficient normal-coordinate based computations that have been performed for FAD [11,[21][22][23] and other systems of similar size [24,25]. In this work, a series of reduced-dimensional curvilinear internal coordinate models are defined and we solve the corresponding vibrational Schrödinger equation.
For a specific coordinate choice, the GENIUSH program automatically computes the KEO coefficients over a grid, and uses the matrix representation of the Hamiltonian (Eq. 1) constructed with discrete variable representation (DVR) [26] for the active vibrational degrees of freedom. The lowest-energy eigenvalues and corresponding eigenvectors are computed using an iterative Lanczos algorithm [27,28].

A. Curvilinear internal coordinates and body-fixed frame
The GENIUSH program uses the t-vector representation to construct the g, and then the G matrices appearing in the KEO, Eqs. (2)-(4). To compute the vibrational t-vectors, Eq. (3), it is necessary to define the Cartesian coordinates in the body-fixed frame, r i , with respect to the internal coordinates, q k . Based on this definition (added to the program as a subroutine), the coordinate derivatives, the KEO coefficients, and the Hamiltonian matrix terms are constructed in an automated fashion. Any coordinates can be set to constrained, to a fixed value provided by the user, or active, for which an appropriate coordinate range and DVR must be defined.
First, we define the monomer structures using a Z-matrix-type notation. These coordinates belong to the intramolecular modes that are set as active or fixed simultaneously in both monomers to respect the compositional symmetry of the system.
Next, we define the coordinates that describe the relative position and orientation of the two monomers similarly to Refs. [29][30][31].
The monomer coordinate definition is summarized in Figure 1. The coordinate axes and the molecular plane are attached to the OCO fragment in both monomers and the following algorithm is implemented to define the KEO: 1. To define the intra-monomer coordinates, the Cartesian coordinates of both monomers are defined as a function of the interatomic distances, r ij , angles, ϑ ijk , and torsion angles, τ ijkl according to the following expressions (see also   Definition of the intermolecular coordinates, (R, θ, φ, α, β, γ), shown for the equilibrium geometry of the dimer.
Up to this point, both monomers have identical positions.
3. To define the inter-molecular coordinates, we rotate both monomers from their original orientation using the rotation matrices O 1 (0, θ, φ) and O 2 (α, β, γ) parameterized with five Euler angles, (θ, φ, α, β, γ), and we shift monomer B by R in the positive direction along the Z axis ( Figure 2): 4. Shift to the dimer center of mass: 5. (optional) Change of the body fixed frame of the dimer: the overall dimer can be rotated to a new body-fixed frame, for example, to the Eckart frame.
The curvilinear nature of the coordinates results in singularities in the KEO.
We have defined the coordinates so that the singularities (most importantly, at cos θ = ±1 and cos β = ±1) of the KEO are possibly far from the equilibrium structure and the dynamically important coordinate range. For example, it would be a more natural choice to align the C-H bond of both monomers along the z axis, but then the equilibrium structure of the dimer would correspond to θ = 0 (cos θ = 1) where the KEO has a singularity. For this reason, we align the C=O bond in both monomers along the z axis ( Figure 2) and perform the O 1 (0, θ, φ) and O 2 (α, β, γ) rotation from this initial orientation. As a result the values of the angles at the global minimum (Table I)

B. Vibrational models and matrix representation
For a start, the constrained coordinates were fixed at their equilibrium value of the dimer's global minimum. Regarding the active coordinates, their initial range was determined from inspection of 1-dimensional cuts of the QB16-PES ( Figure 3).
Based on these 1D cuts, we may expect that the fingerprint region can be well described. For every coordinate, we performed a 1D vibrational computation using the 1D (unrelaxed) potential energy cut over a broad, physically meaningful interval (highlighted in orange in the figure) that is not affected by an unphysical behavior of the PES. For this 1D model, we employed a large number of a discrete variable representation (DVR) points [26] scaled to the selected interval. These computations were used to define potential-optimized DVRs (PO-DVR) [33,34].
This construction (and the semi-rigid properties of the system) allowed us to   retain only a small fraction of grid points for the multi-dimensional vibrational computations. Table I summarizes the coordinate parameters (minimum value and interval), the DVR grid type and the number of points, as well as the number of PO-DVR points that were found to be sufficient to converge the multi-dimensional variational computations presented in this work within 0.01 cm −1 . This convergence threshold is orders of magnitude better than the fitting error reported for the PES, but we used this threshold to make sure that all states appearing in our energy list are true (converged) states corresponding to the vibrational model.  For this reason, we present (Table III) the lowest-energy vibrational energies computed with GENIUSH using the 6D(I) intermolecular, the 8D(Ib) intermolecularbending, the 8D(It) intermolecular-torsional, and the 10D(Itb) intermolecularbending-torsional models and the full QB16-PES. For exploratory reasons, we have constructed other reduced dimensionality models as well, e.g., by including the anharmonic O-H stretching modes, but they did not seem to qualitatively change the intermolecular vibrational energy pattern (we note that the O-H tunneling effects may be accounted for using the same methodology, but it is left for future work).
We have also studied the effect of the precise value of the constrained structure by using the equilibrium, r eq , or vibrationally averaged, r 0 , parameters (Table II).
Since we observed only minor shifts (a few cm −1 ), we continue the discussion for the results obtained with using r eq values for the constraints.
Adding the HOCO torsional modes to the active degrees of freedom appears to be qualitatively important in the present case. At the same time, the OCO bending modes have only a small effect on the energy pattern. We do not provide a detailed assignment for all 6D, 8D, and 10D vibrational computations, but we list (Table III) the computed vibrational energies together with the infrared intensities, to facilitate comparison of the various computed and experimental vibrational band origins (VBOs). In the jet-cooled experiments, we may assume that initially only the vibrational ground state is populated, so we compute only the transition intensity for excitation to the final state ('f') only from the ground state ('0'). The µ α (α = x, y, z)  For the 8D(It) vibrational model, we provide a detailed assignment in Table IV for the computed vibrational states below 350 cm −1 (the first 18 vibrational states above the zero-point vibration). The computed vibrational states were assigned based on their nodal structure. Example wave function plots are shown in Figure 4, in which a clean nodal structure can be observed for the case of the fundamental, the first, the second, and the third overtone of the ν 16 intermolecular twist vibration.
Further wave plots used for the assignment are deposited in the Supplementary

Information.
We see a less clear nodal pattern for the case of the close-lying states with 191 and 208 cm −1 vibrational energies, nodal features along the R intermolecular stretching appear in both states. The latter band appears to have a more pronounced stretching character, although making a decision between the in-plane bending (ν 9 ) and intermolecular stretching (ν 8 ) would be ambiguous based on these results. The two modes belong to the same A g irreducible representation (irrep) of the C 2h point group, and thus their mixing is allowed (and a mixing effect has been observed already in lower-level electronic structure computations [1]). In comparison with the experimental results, we observe a strong and most likely erroneous blueshift.
The erroneous behavior is indicated by the strong, positive 'anharmonicity correction' (deviation of the variational energy from the harmonic oscillator energy) that indicates a problem in the theoretical description.
We can identify all experimentally observed fundamental vibrations, overtone, and combination bands in this range, except for 319 cm −1 peak tentatively assigned to the 2ν 9 (A g ) band in Ref. [9]. This VBO may be missing from our energy list shown up to 347 cm −1 due to the erroneous blueshift of the ν 9 fundamental vibration in the computations. Of course, the computed energy list contains a few more combination and overtone bands that have not been observed in experiment yet.
So, apart from the erroneous ν 9 /ν 8 system, the 8D(It) computational results appear to be in a reasonable, although not spectacularly good, agreement with experiment.
Switching on the OCO bending modes, resulting in the 10D(Itb) model, does not change much the computed results. On the contrary, by freezing both the bending and the torsional degrees of freedom at their equilibrium values, resulting in the 6D(I) model, we obtain a qualitatively different energy pattern and the problematic blue-shifts appear to be even more pronounced (Table III). We find this result surprising, we would have expected that already the 6D(I) intermolecular model is qualitatively correct and it can be further improved by adding the lowfrequency intramolecular modes. To better understand the origin of this behavior, we performed a couple of test calculations using different KEOs that are reported and analyzed in the Section IV.
In summary, 8D(It) appears to be the simplest vibrational model for which we get meaningful results with the present QB16-PES, but the results are far from perfect.
For this reason we constrained the discussion to the first 18 states of the parent