The first microsolvation step for furans: New experiments and benchmarking strategies

The site-specific first microsolvation step of furan and some of its derivatives with methanol is explored to benchmark the ability of quantum-chemical methods to describe the structure, energetics, and vibrational spectrum at low temperature. Infrared and microwave spectra in supersonic jet expansions are used to quantify the docking preference and some relevant quantum states of the model complexes. Microwave spectroscopy strictly rules out in-plane docking of methanol as opposed to the top coordination of the aromatic ring. Contrasting comparison strategies, which emphasize either the experimental or the theoretical input, are explored. treatment itself is largely justified for the zero-point energy, likely and by design due to the systematic cancellation of important anharmonic contributions between the docking variants. Therefore, discrepancies between experiment and theory for the isomer abundance are tenta-tively assigned to electronic structure deficiencies, but uncertainties remain on the nuclear dynamics side. Attempts to include anharmonic contributions indicate that for systems of this size, a uniform treatment of anharmonicity with systematically improved performance is not yet in sight.


I. INTRODUCTION
Any rigorous comparison of electronic structure theory and experiment in the field of solvation has to cope with zero-point motion and thermal motion of the nuclei, before firm conclusions about the underlying electronic structure theory are possible. Thermal motion can be minimized by preparing solute-solvent clusters in supersonic jet expansions. 1 Zero-point motion effects 2 can often be approximated harmonically, in particular, when different solvation sites for the same species are compared such that most of the anharmonic contributions cancel out to a significant degree. 3 However, this approximation needs careful verification, e.g., by isotope substitution, because the non-covalent binding is soft and of large amplitude.
With the goal of elucidating the role of electronic structure and zero-point motion effects in mind, a blind challenge on the microsolvation of furans by methanol (MeOH) and its O-deuterated isotopolog (MeOD) was launched recently. 4,5 Several theory groups predicted the docking preference of the first methanol molecule on furan 6,7 (Fu), 2-methylfuran 8,9 (MFu), and 2,5-dimethylfuran 9,10 (DMFu), while the corresponding mixed complexes were independently characterized by Fourier transform infrared (FTIR) jet spectroscopy. The initial results were published in Ref. 5, but the comparison between theory and experiment raised several challenging questions on both sides, which the present contribution addresses. Apart from technical improvements for several of the theoretical submissions, which are described in Sec. IV, the original predictions from Ref. 5 are used.
The experimental challenge is twofold. Competing solvation structures need to be spectrally identified and ranked in relative energy. Electronic structure theory can help in the first task and can then be tested with respect to the experimental energy ordering, provided that zero-point energy can be reasonably well described and the isomer interconversion barriers are sufficiently low. Quantitative spectroscopic binding energies are usually restricted to special cases, [11][12][13] but for subtle energy differences between isomers, qualitative and semiquantitative rankings such as stability sequences and trends can already be helpful.
The stepwise solvation of furans with hydrogen bond donors provides a good example for subtle solvation preferences. 14,15 A protic solvent molecule X-H may approach either one of the four C atoms or the O atom of the furan ring with its hydrogen atom positioned on top of the ring (t), thus overlapping with a section of the π cloud. Alternatively, it may dock sideways on the O atom, overlapping mostly with its σ lone electron pair in the furan plane (p), as illustrated in Fig. 1. In each case, the solvent residual X can adjust its position to optimize secondary interactions with the ring or with its substituents R. As the primary binding preference between C-and O-docking is delicate, such secondary R⋯X interactions may easily tip the balance in favor of one or the other docking site. This can be used to investigate subtle secondary interactions between the solvent and the solute. 3,16 Because the barriers connecting the non-equivalent docking sites are usually rather low, X-H will preferentially localize in the lowest potential well upon rapid cooling and aggregation of the two binding partners in a supersonic expansion. 3 Depending on the carrier gas, some secondary population will still survive even behind shallow barriers in metastable solvation sites [17][18][19] because the number of cooling collisions on the microsecond timescale of the supersonic expansion is finite and relaxation from initial docking sites to the most stable one may come to a halt. For a single metastable conformation, this can be expressed as an effective conformational freezing temperature Tc (vide infra) for a Boltzmann-type distribution among the two conformational populations. 3 The higher and the wider the barrier, the higher this temperature will be. If more than one metastable conformation survives, the analysis becomes more complex. In Sec. IR1 of the FIG. 1. Schematic overview illustrating the structural essence of the binding situation between a solvent molecule (X-H) such as methanol (X= =MeO) and furans. The solvent molecule may form a hydrogen bond to the π cloud over one of the carbon atoms (label "C") or to the π or σ electrons of the oxygen atom (label "O"). The small label "t" indicates a hydrogen bond on top of the furan plane, i.e., to the π system of the furan, whereas "p" indicates a hydrogen bond in the furan plane, i.e., to a σ orbital. The positions at which the furan (Fu) is substituted to form methylfuran (MFu) or dimethylfuran (DMFu) are indicated by the labels "R." The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp supplementary material (IR), we show that this is actually the case for the methanol complex of MFu. In the spirit of the benchmarking goal, only clear situations should be included, and this complex is better omitted from the analysis. Different supersonic jet spectroscopies can be used to detect the competing docking structures. Linear infrared (IR) spectroscopy 20 provides separate X-H stretching signals for different docking conformations and helps to quantify the frozen low-temperature population ratio because relative IR absorption band strengths I OH can be predicted quite accurately. Coherent microwave (MW) spectroscopy 21 provides direct structure data 22 and also some semiquantitative information on the relative abundance of the docking species based on the calculated dipole moment components μα and rotational constants A, B, C. A combination of low-temperature gas phase IR and MW spectroscopy thus allows us to compare different R⋯X interactions in energy and structure.
After clarifying the structural preferences of initial Fu and DMFu microsolvation by methanol, the present contribution explores different meeting points between theory and experiment, all resting on a simplified two-state Boltzmann distribution for the O/C docking ratio N O /N C = exp(−ΔE 0 (O − C)/RTc), which neglects differences in the rovibrational partition function due to the low temperature and high similarity of the compared complexes. 3 This cancellation of thermal excitation effects among the isomers may not be perfect, but the residual vibrational and rotational entropy differences are expected to be smaller than or at most comparable to the zero-point energy approximations involved. One can either extract approximate experimental ground state energy differences ΔE 0 from an estimated freezing temperature range, 5 or one can analyze theoretical predictions of ΔE 0 in terms of physical or unphysical (negative or too positive) conformational freezing temperatures Tc. This falls short of reaching what is occasionally termed spectroscopic accuracy for vibrational states because these quantities are based on spectral intensities. Ultimately, only the ratio ΔE 0 /RTc is straightforwardly accessible by the experiment, but already this can be useful to explore limitations of density functional theory (DFT) functionals and approximate wavefunction methods. [23][24][25][26][27][28][29][30][31] The aim of this work is to go beyond the benchmarking solely based on docking preferences and evaluate the different theoretical approaches on a multitude of experimentally derived quantities. Of particular interest is the inclusion of microwave spectroscopy data, providing a structural aspect to the theory-experiment challenge. We explore the following: (a) the total energy difference between different structural isomers (as in the first benchmark work) 5 and derived conformational temperature estimates, (b) the description of zero-point vibrational effects by considering different isotopologs (we separate this into two descriptors, deuteration effects and isotope substitution temperature), (c) correlations between spectroscopic parameters, and (d) structural information from microwave spectroscopy (rotational constants and planar moments).
Differentiated ratings of the theoretical model performance for the two experimentally most robust systems are provided in Fig. 11  in the conclusions (Sec. VII). The overall goal is to find a physically motivated protocol yielding a consistently good comparison of theory and experiment. While the present systematic theory-experiment comparison focuses on methanol complexes of Fu and DMFu, the long-term goal is much broader because success for just these two specific systems and their isotopologs could always be fortuitous. The initial experimental studies 4,33 and the first benchmarking steps 5 have indeed already triggered further experimental furan-XOH dimer conformational preference studies. 34,35,38 They are listed in Table I along with systems that were previously studied. 32,33

A. Theory protocol
Most theoretical predictions were taken directly from the original blind challenge. 5 However, in about 40% of the predictions, there were technical changes such as errata, reduction in the number of modes treated anharmonically, inclusion of intramolecular modes into the zero-point vibrational energy (ZPVE) correction, use of an improved algorithm, and tighter optimization or relaxation of symmetry constraints. These are described in more detail in Sec. IV and the supplementary material (TH). To compare experimental data to different theoretical predictions, a strict protocol had to be followed for the theory could also be provided for a single method to be detailed in the supplementary material (TH). To convert the energy predictions to spectral simulations, unscaled harmonic wavenumbers ω OH and infrared band strengths I OH /km mol −1 in the double-harmonic approximation had to be supplied, whereas their anharmonic counterpartsν OH and I anh OH were optional. The lowest (an)harmonic fundamental wavenumber ω low (ν low ) provides information on the floppiness or transition state character of the isomer. To facilitate the comparison of the IR spectroscopic parameters to the experiment, ω OH , ω low , and I OH were also requested for free methanol, again with the option to provide the corresponding anharmonic values. An important quantity is relative degeneracy, the number of equivalent C1, C2, O1, and O2 structures. If O1 is predicted to have definitive Cs symmetry and C1 clearly C 1 symmetry, then there may be a statistical advantage for C1 by a factor of 2 as it comes in two enantiomeric forms. In this case, Tc determination in the Boltzmann approximation requires to divide the measured intensity of C1 by 2, relative to O1. If both isomers are predicted to have the same symmetry, this correction does not apply. Optionally, electronic conversion barriers from C2 to C1, from O2 to O1, and from C1 to O1 could be provided, as they give valuable information on the expected conformational temperature range. For the purpose of microwave spectroscopy, rotational constants A, B, C and dipole moment components μa, μ b , μc along the principal axes of inertia had to be provided at least in the rigid body approximation (el, in the following labeled by an index "e"), optionally also after anharmonic ZPV averaging (0). Finally, Cartesian coordinates x, y, z, by preference in the inertial axis frame, had to be submitted in a standardized format. The supplementary material (TH) also contains details about the employed methodsincluding relevant references.

B. Microwave spectroscopy protocol
The complexes of methanol with Fu and DMFu were studied using chirped pulse Fourier transform microwave (CP-FTMW) spectroscopy between 2 GHz and 12 GHz with 1.0 × 10 6 and 3.6 × 10 6 averages per experiment, respectively. The details of the setup are described elsewhere. 21,40 Methanol (Aldrich, 99.8%) or methanol-13 C (Aldrich, 99%) was mixed with 0.5% Fu (Aldrich, 99%) or DMFu (Aldrich, 99%) in a carrier gas and expanded through a pulsed pin-hole nozzle with an orifice of 1.1 mm (General valve Series 9) at 2.5 bars backing pressure of the carrier gas (3 bars for DMFu-MeOH). For Fu-MeOH, individual experiments with both helium and neon as carrier gas were performed.
For the analysis of the spectra, the PGOPHER 41 and JB95 42 programs were used. Refined fits (DMFu-MeOH) were obtained with Pickett's SPFIT/SPCAT programs 43 assuming a rigid rotor Hamiltonian with a Watson A reduction in the I r representation, with centrifugal distortion terms included. A further, detailed analysis of the observed tunneling splitting arising from internal rotation of the methanol methyl group in the furan-methanol Ct complex (vide infra) as well as the three-rotor problem in DMFu-MeOH was performed using the XIAM program. 44 XIAM is a least squares fitting program specifically designed for the fitting of up to three internal rotors by employing the combined axis method of Woods 45,46 to account for internal rotation through a potential barrier in C 1 or Cs symmetry, as in the case of DMFu-MeOH.

C. FTIR spectroscopy protocol
The FTIR measurements of mixtures of furan (Alfa Aesar, 99%, 250 ppm BHT) or 2-methylfuran (Roth, ≥99%) with methanol (Sigma Aldrich, ≥99.8%) or methanol-d1 (eurisotop, 99% D, HDO + D 2 O < 0.1%) in varying concentrations of ≈0.1% in helium (Linde, 99.996%) were carried out with the setup already described in the original blind challenge report. 5 A supersonic jet expansion of the gas mixtures through a 600 × 0.2 mm 2 slit nozzle is probed by the mildly focused IR beam of a Bruker IFS 66v/S spectrometer. 20 For a single spectrum, 175-700 scans were averaged. The low analyte concentrations lead to isolated molecules and small clusters (mostly dimers) at low conformational temperatures (vide infra), thus providing a good meeting point with quantum chemical calculations. For each donor-acceptor combination, several spectra at slightly different concentrations were recorded.
As documented in the blind challenge report, 5 the IR measurements were characterized by a high dilution and therefore a poor S/N ratio for Fu and MFu to ensure that there are no interferences from larger clusters. DMFu had already been remeasured at an intermediate concentration, resulting in a smaller statistical error bar. 5 Here, we add the results of similar remeasurements for Fu and MFu in combination with regular and deuterated methanol. Strictly speaking, this reduces the double-blind character of the original study. The previous and new measurements were incorporated into the best estimate for the docking isomer intensity ratio, weighted by their inverse variance, as before. The calculation of the intensity ratios and their error bars is detailed in Sec. IR2 of the supplementary material (IR). The new values overlap with the former results within their combined error bars in all cases, suggesting that our previous error estimates had been realistic. of Chemical Physics

III. EXPERIMENTAL RESULTS AND THEIR IMPLICATIONS FOR THE STRUCTURES AND ABUNDANCES
A. Rotational spectroscopy of methanol/furan The first round of the blind challenge 5 led to an ambiguity between the majority of harmonically evaluated predictions, which favor a π or top coordination of the ether oxygen by methanol (Ot), and a single anharmonic calculation, which favors a more planeparallel and thus more σ hydrogen-bonded coordination (Op). Although the infrared spectra have hinted at Ot because the small OH downshift indicated competition between hydrogen bonding and methyl dispersion interaction, conclusive evidence requires rotational resolution. Microwave spectroscopy can easily distinguish between such isomers, whereas dynamics that is much faster than nanoseconds, such as the question of true vs effective symmetry in low barrier situations, can be more challenging. We come back to that topic toward the end of this section.
Part of the experimental spectrum obtained from a Fu-MeOH mixture is displayed in Fig. 2 (top trace). Its analysis results in the assignment of two dimer species with pronounced a-type transitions, which point to significant μa dipole moment components. Simulations based on the fitted rotational parameters (Table II and  Tables MW2 and MW3 Tables MW17-MW22) for the two Fu-MeOH complexes as well as singly substituted 13 C isotopologs of the Ot complex (vide infra) are given in the lower trace of Fig. 2. By comparison of the experimental rotational constants with the results from the quantum-chemical calculations, these two dimer conformations can be assigned to Ot and Ct, both involving a top coordination of the furan ring (Table II). From an early stage in our analysis, we can rule out the possibility that the observed rotational transitions belong to the Op isomer. For the assigned isomers of Fu-MeOH, the A rotational constant of Op is significantly larger in the prediction compared to the experiment (≫10%). Additionally, in the warmer helium expansion where higher-energy isomers are less likely to relax to lower-energy isomers, the Op isomer was not found. As such, only rigorous comparisons to the Ot and Ct isomers are carried out.
Among the two observed isomers, Ot clearly dominates the spectrum, and Ct forms a minor constituent. Relative abundances are determined using predicted dipole moments and transition intensities (only a-type transitions were used here). Using predicted μa dipole moment components from each of the computational methods, relative abundances can range from 4.8(Ot):1(Ct) to 1.7(Ot):1(Ct). An overview of each methods' performance is given in Table MW1 of the supplementary material (MW). The different prediction methods used in this study give either C 1 or Cs symmetries for Ot(Fu-MeOH), which gives a larger uncertainty in the relative abundance determined by the microwave experiment. Despite this uncertainty, and other small ones also described in the supplementary material (MW), Ot is determined to be preferred over Ct, regardless of the source of the predicted dipole moments. This confirms the infrared results of a more stable Ot isomer. The use of different carrier gases in the expansion can provide useful insights into the energy ordering of molecular isomers or conformers. For Fu-MeOH, we observe no major change in the relative abundance ratio for the two complexes, Ot and Ct, under the two different expansion conditions. However, for the helium measurement, only a smaller frequency region was probed so that only one set of rotational transitions could be compared, which limits the statistics, but certainly allows us to determine an overall trend. A comparison of these spectra can be found in Fig. MW1 of the supplementary material (MW).
Further experimental information on the structures and symmetries of the clusters can be obtained by isotopic substitution. The high dynamic range of the observed rotational spectrum allowed for the detection of Ot 13 C singly substituted isotopologs in natural abundance (1.1%). Using these additional datasets, the carbon backbone of the Ot isomer can be determined. A first hint at the effective symmetry of the complex was found from the intensity of the 13 C isotopologs. They exhibited twice the expected abundance relative to the parent isotopolog, i.e., 1:50 instead of 1:100, which is due to a pair-wise equivalence of the carbon atoms of the furan ring and thus effective Cs symmetry. Additional structure information was obtained with 13 C-labelled methanol. The corresponding rotational constants for both isomers and their 13 C isotopologs, Ot and Ct, are listed in Tables MW2 and MW3 of the supplementary material (MW).
The changes in the moments of inertia of the clusters upon isotopic substitution enable us to determine the experimental atom positions in the principal axis system by using Kraitchman equations (rs) 47 and least-squares fitting routines. The two most common least-squares fitting methods result in (a) the effective ground state structure, r 0 , and (b) the so-called mass-dependent r (1) m structure, which allows us to compare experimental ground state structures with computed equilibrium geometries accounting for vibrational contributions to the observed moments of inertia. [47][48][49][50] Detailed parameters from the structure determination and their respective coordinates can be found in the supplementary material (MW) for the Ot (Tables MW6-MW12) and Ct (Tables MW13-MW16) structures.
Comparing both the Kraitchman and least-squares methods (rs and r Most of the quantum-chemical calculations employed in the blind challenge predict a non-symmetric minimum structure with the carbon atom of methanol being somewhat off the vertical symmetry plane of Fu. This may well be compatible with effective Cs symmetry in the vibrational ground state, as in related cases of quasi-Cs-symmetric molecules such as 2-oxazoline, 51 oxetane, 52 or also the water pentamer. 53 If the barrier is lower than the zero-point energy along the planarization coordinate, as in the former two examples, the microwave spectra can still be qualitatively similar to those of a rigid planar rotor despite a shallow double-minimum potential. Even extremely floppy molecules can display quite regular spectra at low temperatures. 54 Computed barriers to Cs symmetry for the Ot complex were not requested from the theory groups in the double-blind study, even though the predicted barrier would help to understand the structure of the Ot complex. The C 1 symmetric structure of Ct is shown in Fig. MW5 using the single rs coordinate for 13 C methanol as a qualitative guide. Structural comparisons in the Fu-MeOH complexes serve as an example for the needed symmetry considerations when predicting equilibrium structures that will be used to simulate ground state rotational levels.
Both Ot and Ct show internal motion of the methanol methyl group, resulting in characteristic line splitting in the rotational spectra into two tunneling components designated as A and E states. In Table II, we present two fits for the Ct isomer, where fit 1 only considers the A state components, which are typically less affected by internal motion than the E state components and which result in useful effective rotational constants. In fit 2, obtained using the program XIAM, both tunneling components A and E are fit together using additional Hamiltonian terms to describe the torsion-rotation coupling. The determined internal rotation barrier, V 3 , of about 2.7956(29) kJ/mol [233.69(24) cm −1 ] is a typical value for the internal rotation of the methanol methyl group in complexes. For the Ot complex, a global fit accounting for such internal motion could not be achieved, which points to the existence of additional large amplitude motions within the cluster. These motions may also involve the vibration of the methanol unit across a shallow barrier at the symmetry plane of the complex, as it is predicted by several calculations. If this barrier is low enough, the zero-point energy for this motion will exceed the barrier height. The first excited state then still corresponds to vibrational excitation and may not be populated significantly in the jet or have rather different rotational constants. In contrast to a tunneling excitation in the high barrier limit, it could escape detection and assignment. Such vibrational dynamics on the picosecond timescale only has indirect effects on ground state rotational spectra with their nanosecond dynamics and may explain the better standard deviation (one order of magnitude) attained using the Cs structure over the C 1 structure in the r For the Ct isomer, the overall lower line intensities compared to Ot in the spectrum prevented us from observing the singly substituted 13 C isotopologs in natural abundance. Therefore, the structural analysis was based on the isotopically labeled methanol 13 CH 3 OH. Better results were obtained using the r Determining the symmetry of the Ct complex is complicated by furan being such an oblate molecule that its orientation in the Ct cluster is ambiguous, as further discussed in the supplementary material (MW). We cannot provide a clear statement if the Ct complex is of C 1 or of Cs symmetry on the timescale of our microwave spectroscopy experiment. The set of rotational constants (determined with and without considering internal rotation, B. Rotational spectroscopy of methanol/2,5-dimethylfuran The analysis of the DMFu-MeOH broadband rotational spectrum revealed the presence of one strong spectrum with a rich fine structure due to internal dynamics, which we can assign to the Ot complex, as displayed in Fig. 3, based on a comparison with the quantum-chemical calculations. No second isomer, such as Ct, could be identified in the rotational spectra, and the Ot complex can be concluded to be the global minimum. The analysis of the rich fine structure due to the internal rotation of the three methyl groups in the complex is illustrated in more detail in Fig. MW6 of the supplementary material (MW). The results of the global fit including internal dynamics are summarized in Table III.
As in the case of the Ct isomer of Fu-MeOH, the overall intensity of the transitions did not allow us to obtain the 13 C isotopologs in natural abundance so that a detailed structure analysis, as achieved for the Ot isomer of Fu-MeOH, is not possible here. The analysis of the internal dynamics and the magnitudes of the dipole-moment components, however, provides good evidence that we obtain an overall Cs symmetric structure of the Ot isomer of the DMFu-MeOH complex, as discussed in the supplementary material (MW).
In summary, the most impactful result of the microwave investigation for the relative energy benchmark is that methanol preferentially coordinates the O atom of Fu from the top, and not in the ring  plane, and that the same is the case for DMFu. Carbon coordination is less abundant (Fu) or not detectable (DMFu).

C. New benchmark-relevant FTIR results
The previously reported FTIR results concerning band positions and assignments 5 remain unchanged, but the remeasurements for MFu reveal a third binary complex of an unclear docking structure (O2 or C2) in the spectral region of C1. This is elaborated upon in Sec. IR1 of the supplementary material (IR). The unclear docking structure leads to a lower bound for the O/C intensity ratio, when the new peak is included in the C-docking integral. It is thus listed with a ">" sign in Table IV and not considered further in the relative energy benchmark. This table compares intensity ratios including the improved spectra with the preliminary ones reported before. 5 In all cases, the preliminary and new ratios overlap within their error bars.
An interesting consequence of the reduced statistical error bar is that, now, the deuteration effect becomes statistically significant for Fu. Satisfactorily, it has the same direction as in DMFu and further supports the common O-docking assignment of the higher frequency signal due to the more localized hydrogen bond.

Band intensity ratio Fu MFu DMFu
2.4(2) Preliminary 5 3.3(10) 1.9(5) 2.4 (2) Whether or not a symmetry correction factor of 2 has to be applied to the IR intensity ratio when converting it to Boltzmann energy differences has to remain open. If the O-structure is Cs symmetric and the C-structure is not, the ratio entering a Boltzmann analysis must be doubled because the chiral C-docking structure is enantiomerically degenerate and the O-structure is not. While the Ot methanol complex with Fu has been unambiguously shown to have effective Cs symmetry (vide supra), Cs symmetry cannot be strictly ruled out for the corresponding Ct structure, although it is less likely. Furthermore, the large amplitude motion in the Ot structure may increase its partition function and thus qualitatively compensate for the symmetry disadvantage. For DMFu, there is no evidence for the Ct structure from the microwave experiment, again leaving open the symmetry correction factor. While no microwave structure data are available for Fu and DMFu complexes of deuterated methanol, the new IR data support that there is no switch in docking preference by deuteration.
As a cautionary consequence, the Fu and DMFu intensity ratios in Table IV will be converted into Boltzmann energy differences both with and without a symmetry correction factor of 2, thus increasing the experimental indeterminacy.

IV. EXPERIMENT APPROACHES THEORY FOR THE SOLVATION OF SYMMETRIC FURANS
There are different ways to evaluate the performance of theoretical methods in this furan microsolvation benchmark. The simplest is to count how many of the former individual computational results 5 for Fu, MFu, and DMFu predict the correct global minimum structure. Assuming no major IR band strength anomalies and incorporating the microwave proof for an Ot structure in MeOH-Fu and to some extent also MeOH-DMFu, we now have substantial evidence that all 1:1 complexes, including the deuterated variants and likely also the classical structures without ZPVE, prefer the O-docking structure where the methyl group of methanol is bent over the ring. Without ZPVE, this was predicted for 27 out of the 35 original calculations. With addition of ZPVE (anharmonic in 9 cases and harmonic in the other 26), this number dropped to 23. Less than one third of all methods made this prediction consistently for all three complexes 5 (see also Fig. 4). However, this simplified binary analysis falls short of what is available from the experiment, and at the same time, it runs the risk of overrating some approximations.
In the preceding paper, 5 an attempt was made to estimate an experimental energy range based on a plausible (20-100 K) range of conformational temperatures and the calculated IR intensities of a selected anharmonic method (building on contribution C). Only about one quarter of all methods passed this test for all three complexes. As the original 5 anharmonic variant of method C is now demonstrated with the help of microwave spectroscopy to make the wrong structural prediction, this analysis shall not be adopted to the present improved experimental data. However, we will repeat the more conservative analysis presented in Table II of Ref. 5, which was restricted to purely experimental evidence and conservative estimates of missing information on IR band strength ratios (up to ±50% deviation from the value 1.0) and plausible conformational temperatures (Tc = 60 ± 40 K) together with worst case error propagation. This has led to experimental error bars for the docking isomer energy difference, which span both signs, although quite asymmetrically. Because the error bar is dominated by Tc indeterminacy, the improved experimental IR intensities from this work for Fu only lead to small reductions in the total experimental error (0.6 instead of 0.8 kJ/mol). The results are shown in Table V for the two symmetric furans for which firm experimental conclusions are possible.
Among the original theoretical entries A-L, entry J is omitted in this work because it essentially provides an improved structure for method I. Where some kind of anharmonic calculation was made available, the entry is marked with a " * ," whereas the harmonic values are displayed without the " * ." Therefore, methods A, C, D, and F now have two entries each, and method L is equivalent to K * . Where available, results for MeOD are added in the lower part of this table. Wherever there is a qualitative contradiction with the experiment in terms of preferred docking geometry at a given atom type (C or O), the entry for that species is replaced by a ". . .". This is  5 and detailed in the supplementary material (TH)] were used. Entries where the results are qualitatively incompatible with the results of the experiments are marked with ". . .". "HO" means the harmonic oscillator model for the calculation of the ZPVE, whereas a " * " marks different kinds of anharmonic corrections. "VCI" indicates the vibrational configuration interaction method and "VPT2" vibrational second-order perturbation theory. The utilization of an anharmonic lower level model Hamiltonian to obtain a ZPVE scaling factor is indicated by "AS," and hindered rotor calculations are marked by "HR." For each of the donor molecules, two estimated experimental energy difference spans (in kJ/mol) are given-once assuming C 1 symmetric O and C structures and once for a C 1 symmetric C in combination with a Cs symmetric O structure. To avoid rounding artifacts, all energies are provided with two decimal places, irrespective of their significance. the case when theory predicts no C structure (A, A * ) or the Op structure as the more stable O-docking conformation, instead of the microwave-verified Ot structure for furan (A, A * , K). K * is a borderline case because the anharmonic analysis for the Op structure is unavailable. We therefore keep it in the discussion because Ot might become the global minimum structure on the anharmonic level. The pre-selection of methods is also described in Sec. TH1 of the supplementary material (TH). When no value for a certain species was submitted, this is indicated by a blank space. The DMFu value for entry A has changed due to an erratum. Entry C * now replaces the VSCF-based anharmonic correction of C for all modes in the original benchmark study 5 by a mixed harmonic/anharmonic When compared to the experimental values assuming no symmetry degeneracy, only methods D, I, and K * reproduce both MeOH complexes within the error bars and only methods I and K reproduce both MeOD complexes. This is further discussed in the supplementary material (IR). If Ot symmetry degeneracy is assumed, method E can be added for MeOH and MeOD, method K * for MeOD in addition to MeOH, and method F for MeOD only. Thus, only methods E, I, and K * provide fully consistent predictions for all investigated species, but only method I does so quite independently on the symmetry issue. One may argue that Cs symmetry is experimentally proven for the MeOH-Fu complex and to some extent also for MeOH-DMFu, but there are several reasons why we hesitate to impose it systematically in this initial comparison: It has not been experimentally proven for MeOD-Fu or MeOD-DMFu, so far. Furthermore, the symmetry disadvantage of Ot might be partially compensated by relaxation of the competing Op structure into Ot at a late stage of the expansion or by excessive floppiness along the Op-Ot interconversion coordinate. Finally, there is microwave evidence for floppiness along the double-minimum potential coordinate breaking Cs symmetry (vide supra). All this may give Ot more statistical abundance even at low temperature than provided by a simple harmonic Boltzmann prediction.
Note that most successful methods treat ZPVE harmonically, whereas of the anharmonic calculations, only method K * matches the experiment in a systematic way. Furthermore, the proposed anharmonic corrections are quite irregular in size and sign, depending on the employed method and system. In part, this irregularity stems from the difficulty of treating the strongly coupled and largely delocalized large-amplitude motion, which does not affect the isomers in a systematic way. For example, it has been documented long ago 57 that VSCF can perform poorly for low frequency modes of neutral hydrogen bonded complexes and this has recently been confirmed for a well-characterized system. 58,59 Overall, the unsatisfying performance of the anharmonic methods could be interpreted as a lack of mature global anharmonic treatments for systems of this size, or it may be due to fortuitous error compensation in the successful harmonic methods. The rather regular IR spectral patterns across deuteration and methylation point to the absence of major anharmonic effects on the relative docking propensity, but only an extension of the experimental database to further systems will be able to support or disprove this.
For MeOH, where the most complete dataset is available, the results of Table V are graphically represented in Fig. 4. A similar figure for MeOD can be found in Fig. IR2 of the supplementary material (IR). MeOD predictions for method D would be required to decide whether it is comparable to or even superior to methods E, I, and K * .

V. THEORY APPROACHES EXPERIMENT FOR THE SOLVATION OF SYMMETRIC FURANS
The previous analysis is attractive for theory because it puts all the burden on the experiment and just asks for electronic energy differences, complemented with some estimate of the ZPVE difference between the competing docking isomers. It would be even more attractive if ZPVE issues could be separated based on the experiment, which is far from trivial. The MeOD data come closer to the electronic limit and essentially suggest that the best O-docking isomer is electronically 1 ± 1 kJ/mol more stable than the best Cdocking isomer for both furans. This is met by methods C, D, E, F, H, I, and K, whereas methods A, B, and G tend to favor C-docking in at least one of the cases.
The experiment can be more rigorous in ruling out theoretical approaches if all the burden of complex geometry and symmetry, electronic energy, ZPVE, and IR band strength is put on the theoretical side and the predictions from different groups are individually combined with experimental IR intensity ratios [within their error bar, see Sec. IR2 of the supplementary material (IR)] to predict the most indeterminate experimental quantity, conformational freezing temperature Tc. For a consistent high quality prediction, each individual Tc value should be strictly below the starting temperature of the expansion (here 300 K) and above the effective rotational temperature (here ≈10 K) corresponding to barrierless relaxation. Negative Tc values are equivalent to a wrong energy sequence prediction. More realistically, the low, but finite barrier situation in furan-methanol complexes requires high quality theoretical models to predict 20 K < Tc < 100 K, as discussed before. 3,5 If the predicted Tc values for MeOH and MeOD are similar, this indicates a reasonable description of experimental zero-point energy effects by the employed ZPVE model. If theory predicts a switch in docking preference upon deuteration at variance with the experiment, the two isotopologs' Tc values differ in sign. If the predicted Tc values across homologous systems are similar, this implies similar barriers for conformational interconversion. If they differ strongly among homologous compounds, either the barriers vary with substitution or theory does not capture the substitution effect properly.
We should note that the proposed Tc prediction method fails for a number of special cases, which are discussed in Sec. IR4 of the supplementary material (IR) because they are not critical in the present case. Figure 5 plots the obtained Tc values for MeOD vs those for MeOH for different levels of harmonic and anharmonic ( * ) theory and the two furans (Fu squares, DMFu diamonds). Symmetry is used in the prediction whenever it is found in the calculation based on the supplied structures as detailed in Sec. TH14 of the supplementary material (TH). While the resulting Tc values range from about −130 K (unphysical) to +160 K (large, but still physical), in particular, the positive ones fall close to the bisector of isotope-independent temperature. This provides experimental evidence for a reasonably accurate and regular theoretical description of ZPVE effects by the employed methods, indicating that major deviations from the experiment might be blamed on structural and electronic deficiencies.
Methods C, C * , F * , G, and H are found to significantly overestimate the stability of the C-docking isomer for at least one furan species. Among the methods for which MeOD predictions are available, only method E lies inside the narrow experimental window for both furans, with methods I and K * coming close to it, followed by method F. Generally, DMFu calculations provide more physical conformational temperatures than Fu calculations. Anharmonic error cancellation appears to work better for DMFu than for Fu. The performance of theoretical methods for zero-point vibrational energy effects can be more clearly separated from electronic structure performance by plotting the conformational temperature T Δ , which is obtained from the deuteration effect on the theoretical energy difference itself together with the deuteration effect on the experimental concentration ratio, thus removing the electronic offset. If the symmetry of the complexes does not change upon deuteration, this analysis is also independent of symmetry corrections. The resulting T Δ values are the same as those in Fig. 5, if the latter lie exactly on the bisector, otherwise they can deviate quite strongly, and the deviation is sensitive to experimental integration error. T Δ is only a measure of how well theory describes the deuteration effect and not on how well it describes the electronic energy difference between the docking sites. If T Δ is negative or larger than the nozzle temperature, this suggests a ZPVE calculation, which disagrees with the experiment. Therefore, we plot the T Δ results for those submissions that provide predictions for all four species (Fu, DMFu with MeOH, and MeOD) in Fig. 6. They are seen to be quite reasonable in all cases (light gray), and within their experimental error bars, they touch the expected temperature window (white area).

FIG. 6.
Correlation of effective isotope substitution temperatures for MeOH docking on Fu (horizontal) and DMFu (vertical axis) with maximum error bars due to the experimental intensity uncertainty for different harmonic and anharmonic ( * ) computational levels, which are structurally consistent with the experiment. All predictions are centered in the physically reasonable temperature range (light gray), and their error margin touches the expected temperature range (white).
The best performance is found for I, with method F offering the most consistent performance for both furans. Interestingly, both use SCS-MP2 or -CC2 harmonic frequencies, but the DFT-based predictions are comparable within the large error bars. In principle, the effect of deuteration on the energy difference between two docking isomers, which enters this analysis, can also be estimated experimentally by analyzing the spectral shifts between the isomers for all relevant fundamental excitations, most importantly for hydrogen bond modes and OH/OD stretching modes. This would allow for a purely experimental estimate of the conformational temperature 60 but is quite challenging in practice. Even without this experimental option, Fig. 6 supports the view that the included methods describe zero-point energy effects acceptably and that this study can therefore provide performance information on electronic energy as well. Figure 7 finally plots the Tc values for MeOH complexes with DMFu vs those for Fu. It analyses whether the combined prediction of electronic energy difference, zero-point energy difference, symmetry, and infrared band strength ratio for the two isomers is consistent with the observed infrared intensity ratio according to the simplified two-state Boltzmann distribution (see Sec. I). Most Fu predictions yield negative Tc values, indicating that they do not get the energy sequence right. I and K * predict very low positive conformational temperatures, and only E falls in the experimental window. For DMFu, the overall performance is better, with negative Tc values (wrong energy sequence) only for G and H. C, C * , D, E, and F match the narrow experimental window for DMFu. Therefore, E provides the most consistent match with the experiment for MeOH, followed by I, K * , and D.
In summary, the theory-approaches-experiment strategy identifies I, K * , D, and, in particular, E as the most promising methods to describe the structure, energetics, and spectral intensities of of Chemical Physics methanol-furan complexes. They all combine dispersion-corrected DFT or CC2 methods for the structure with CCSD(T) variants for the energy, and only one of them employs anharmonic ZPVE estimates. Extension to further systems (Table I)

VI. THEORY MEETS EXPERIMENT FOR SPECTROSCOPIC PARAMETERS
Besides the fundamental challenge of matching subtle conformational preferences in furan-alcohol complexes, it makes sense to exploit and analyze the available theoretical predictions of spectroscopic observables as assignment tools. This will be attempted in Subsections VI A and VI B for infrared and microwave spectroscopy. While it is very difficult to match experimental differences in OH stretching fundamental wavenumber or OH stretching fundamentals themselves and almost 61 impossible to match rotational constants within experimental accuracy, one can still explore how close the individual theoretical levels come to the experiment. As before, one can compare harmonic/equilibrium values and hope for cancellation or smallness of anharmonic corrections. Alternatively, one can use approximate anharmonic treatments such as vibrational perturbation theory (VPT2) or other approximations proposed by some participants of this blind challenge. We restrict the following discussions on methods, which are in qualitative structural agreement with the experiment (vide supra). The complete set of raw data is available in the supplementary material (TH).

A. OH stretching wavenumbers and downshifts
The most useful and robust theoretical input for a series of homologous complexes involving competing OH docking interactions is the substitution trend in the hydrogen-bonded OH stretching fundamental. The individual hydrogen bond-induced downshifts Δν from the isolated monomer may still be strongly affected by major anharmonic effects, but if the downshifts are similar for different docking situations, one may expect that a balanced theoretical model predicts the difference in downshifts ΔΔν among closely related systems such as Fu and DMFu quite well in the harmonic approximation (ΔΔω). This assumption is tested in Fig. 8. A corresponding figure for deuterated methanol is provided in Sec. IR5 of the supplementary material (IR; Fig. IR3). The underlying spectroscopic data are listed in Table IR1 of the supplementary material (IR), restricted to those theoretical predictions, which agree qualitatively with the experiment in terms of the docking structure. Figure 8 shows how much further the DMFu OH stretching transitions of O and C isomers are shifted to lower wavenumbers than the corresponding Fu transition. Several conclusions can be drawn for the methods used to optimize the structure and calculate the Hessian: • Where anharmonic values ( * ) are available, they correlate less well with the experiment than the associated harmonic ones-the pragmatic spectroscopy approach of using harmonic predictions as relative assignment aids due to the The multi-level strategy to describe structures and energies with different theoretical methods seems to work well for the method combinations in E and I related to OH-stretching and relative energies.
For the prediction of absolute OH stretching fundamental wavenumbers, anharmonic effects have to be captured in some way. Table VI compares such experimental vibrational wavenumbers with anharmonic predictions, where available. Besides the improvement from adding the diagonal anharmonicity correction for the OH oscillator, which can also be achieved by simple scaling, it is difficult to see a systematic relative performance among the different systems and docking sites. Stretching frequencies for O-docking are mostly predicted too high, whereas the predictions for C-docking scatter around the experimental values. It is likely that both offdiagonal anharmonicity contributions and deficiencies in the predicted Hessian are responsible for the deviations, but this needs future systematic investigations.
A room temperature gas phase and matrix isolation study of furan-alcohol complexes 35 qualitatively confirms the importance of jet-cooled spectra for benchmarking (we note that some key entries for the energetics of Fu complexes in Table 2 of Ref. 35 appear to be inconsistent). Shifts in the gas phase are too small due to the thermal excitation of the hydrogen bond, and the bands are too broad to

B. Rotational constants
For benchmarking structures/rotational constants, there are three aspects that complicate a direct comparison of the manifold of experimental results with the calculations: (a) the effective Cs symmetric structures, (b) the methyl group internal rotation, and (c) the fact that the microwave experiments probe the molecules in the ground vibrational state, v = 0, resulting in the rotational constants A 0 , B 0 , and C 0 . The calculations, if not corrected for vibrational ground state averaging, provide equilibrium rotational constants at the non-physical minimum of the potential well. 66 For many molecules, the difference between A 0 , B 0 , C 0 and Ae, Be, Ce is on the order of a few MHz. When comparing the global fit of Fu(Ct) to a fit including only the A state transitions, the differences are almost negligible for B and C and more pronounced for A (with 1.4% difference, see Table II). To take all of these effects into account, we consider a generous error for vibrational averaging effects of −1(±3)% from the equilibrium values to the experiment following the recent suggestion of Oswald and Suhm. 67 Comparison between calculations performed as part of this blind challenge (mostly Ae, Be, and Ce) and our experiment (A 0 , B 0 , and C 0 ) shows that most theoretical approaches predict the rotational constants (i.e., the molecular structures) of the two Fu-MeOH isomers reasonably well with deviations on the order of a few percent, as shown in Fig. MW3 of the supplementary material (MW).
There seems to be a general trend for all calculations that there is better agreement with the experiment for Fu(Ct) than for Fu(Ot). It is a common observation that structures with a rather localized OH⋯O hydrogen bond contain more zero-point vibrational energy than the less localized OH⋯π structures. This trend may also be taken as additional evidence for large amplitude motion in Ot along a shallow symmetric double well, which gives rise to effective Cs symmetry.
There are several different ways to compare experimental and computed rotational constants in molecular complexes such as Fu-MeOH and DMFu-MeOH, and we have to restrict ourselves to a few, referring to the supplementary material (MW) and later work for others.
One is very pragmatic and inspired by the practical usefulness of structural predictions in assisting the assignment of microwave spectra. In Fig. 9, we plot the percent error for the two assigned Fu complexes and the single DMFu complex for the different structure optimization methods listed in Table V in terms of (B + C)/2. This average of the two smaller rotational constants is the most informative in predicting transitions of near-prolate tops with μa dipole moment components (Fu-MeOH), which is why it is chosen for this analysis, even though predictions of the A rotational constant typically have larger errors. From the plot, (Be + Ce)/2 for method C formally appears to best predict (B 0 + C 0 )/2 for all three measured complexes (it is closest to the upper right corner of the gray auxiliary plane), except that it is typical for predicted equilibrium structure rotational constants to be larger than experimental rotational constants, which C fails to achieve for the Fu(Ot) complex (a positive percent error in Fig. 9). of Chemical Physics If one instead allows for a generous uncertainty window for vibrational averaging effects from equilibrium to experimental rotational constants of −1(±3)% for (B + C)/2, methods C, D, E, G, H, and I all fit well for the two Fu complexes and the DMFu complex. Methods F and, in particular, B are outside the vibrational averaging margin for at least one of the complexes. For method F, this somewhat contradicts the very good performance in the highlevel electronic energy test shown in Fig. 5 of Ref. 5. This might be due to differences in performance for stiff and soft degrees of freedom. Overall, the comparison of computed equilibrium structures with experimental ground state rotational constants can profit from anharmonic corrections. 68 In the case of Fu complexes, where anharmonic corrections are available within method F, it is reassuring that they amount to −2(±1)% and thus bring methods D, E, G, H, and I in even more satisfactory agreement with the experiment, when transferred to those methods. Analysis of the A, B, and C rotational constants, instead of (B + C)/2, compared to predicted values from the methods used in Table V is carried out in Table MW4 of the supplementary material (MW).
Another form of analysis is more specific to the investigated systems and builds on the fact that in all three complexes, the aaxis corresponds qualitatively to the dissociation coordinate such that the planar moment Paa describes best how well a computed structure matches the intermolecular distance. Essentially, experimental Paa values give relative comparisons of the intermolecular interactions between the two molecules and can be compared to the computed values to determine how well they predict the overall interaction between the two molecules. This is like the treating of van der Waals complexes in the pseudo-diatomic model to obtain dissociation energies. Planar moments are related to the mass along a given axis by the formula Paa = ∑ i mia 2 i , where ai is the a-coordinate of each atom, i, with the respective mass, mi. Paa can also be obtained experimentally from the measured moments of inertia Paa = 1 2 (I b + Ic − Ia). The side-on view to represent the out of bc-plane mass is shown in Fig. 10. Again, vibrational averaging has to be considered and is likely dominated by anharmonicity such that one expects theoretical equilibrium predictions of Paa to be smaller in comparison to experimental values. This is indeed the case for most predictions, but it is unclear how large the anharmonic effect is. With the assumption that anharmonic ground state vibrational averaging is the largest contributor to the difference between calculated equilibrium and observed ground state planar moments, methods with a smaller, positive percent error in Table VII are considered to be more accurate than those with larger values. Entries D, E, and G are ahead in this analysis, narrowing down on the relatively low discrimination based on (B + C)/2. Entries H and I, which did well in the (B + C)/2 analysis, have slightly larger percent errors in their planar moments, and entry C has the wrong sign.

VII. CONCLUSIONS
The reliable comparison of theory and experiment for the purpose of theory benchmarking requires considerable effort and great care on both sides. Experimental data are susceptible to misinterpretation, 69 and theory has to cope with several layers of approximation. 2 The present work explores different meeting points between theory and experiment for the first microsolvation step for furans, further elaborating on the first results of a recent blind challenge on the topic 5 and pointing the way for future extensions. Blind testing has a longer tradition in the field of crystal structure, 70 where the primary quantities, structure and energy for non-covalent interactions, are well accessible even at low temperature. 71 Supersonic jet techniques give access to such quantities for non-periodic, much smaller model systems and are therefore valuable for theory benchmarking and blind testing purposes, if the non-equilibrium aspects are sufficiently under control. Furan microsolvation is promising in this context, as the results show. Figure 11 provides a color-coded summary of the performance of the methods investigated in the furan microsolvation challenge. Only two methods with a full dataset (methods E and I) show good or reasonable performance for all tested properties, whereas the less expensive method D performs satisfactorily but misses data for isotope substitution effects. Method K * could only be tested for energy and deuteration effect predictions, and all further methods fail for at least two of the employed test protocols. Not surprisingly, the best discrimination of the performance of the methods can be reached if the burden of auxiliary information is put on the side of theory, as the comparison of the benchmarks of the energy predictions via energy difference and conformational temperature illustrates.
Microwave spectroscopy provides valuable constraints on the docking structure and symmetry, but the rigorous benchmarking of rotational constants would require extensive anharmonic modeling, including the effective symmetry issue. Vibrational frequency trends also rule out relatively few methods, when evaluated harmonically.  Table V, which were characterized in Figs. 4-9 and Table VII in terms of their structure, energy, and Hessian performance based on the selected properties. Best performance is marked in green, still reasonable performance in yellow, and poor performance in red. No color means non-coverage due to missing entries or qualitative failure (see Table V for details). Note that this rating may change substantially, once the database is extended beyond the two complexes and their OD isotopologs.

ARTICLE scitation.org/journal/jcp
A long-term challenge is thus the anharmonic modeling of the solvent docking landscape, which includes rather low barriers between in-plane oxygen and on-top oxygen or carbon coordination. Although the present study suggests that a harmonic description is not too far from the reality for energy differences, quantitative and even occasional qualitative consequences of anharmonicity should not be ruled out completely. Already the effective symmetry issues discussed in the microwave Secs. III A and III B point at non-trivial large amplitude effects. Furthermore, success for just two very similar furan-methanol complexes and their isotopologs could always be fortuitous.
A possible strategy to extend the furan microsolvation challenge is to address some of the systems listed in Table I with the most promising computational methods. While some of these systems are substantially larger in size than the present ones and the experimental evidence is frequently limited to a qualitative ranking of different docking structures, they may still serve to discriminate between different computational approaches. The ultimate goal is to find a physically reasonable computational protocol to yield a systematically good performance when comparing to the experiment.

SUPPLEMENTARY MATERIAL
See the supplementary material for additional information on the microwave spectroscopy analysis (MW), the infrared spectroscopy analysis (IR), and the theoretical database (TH).