The furan microsolvation blind challenge for quantum chemical methods: First steps

Error cancellation contributes to the success of quantum chemistry in describing experiment, but it can slow down the progress of computational method development, depending on whether it is systematic or fortuitous. The study of a broad set of similar systems with systematic variations is at the heart of chemistry and helps to differentiate between the two effects. Microsolvation,¹ the stepwise addition of solvent molecules to a solute, is a typical area where small errors in individual interactions might either add up or cancel each other with a growing solvation shell. It therefore calls for benchmarking efforts on a molecular scale, to show which quantum chemical models perform well for a good reason.

Furans were previously identified as suitable benchmark systems^2–5 because they offer two docking choices—similar in energy but different in nature—to a protic solvent molecule. One is the polar oxygen atom (O), with an attenuated attractivity due to aromatic delocalization of its electron density. The other is the aromatic carbon skeleton (C), which gains attractivity from the delocalized π cloud. Structurally, the two isomers of a furan/solvent complex can be easily distinguished since their rotational constants and/or dipole components differ significantly. Vibrational spectroscopy can also sense the difference, in particular, when probing the solvating hydride stretching fundamental vibration under cryogenic conditions. The preferred way to generate such a cold environment for neutral molecules makes use of a supersonic jet, i.e., collisional cooling and aggregation in an adiabatic carrier gas expansion. In principle, several O- and C-docking isomers could be formed, but the barriers connecting them are usually rather shallow and narrow such that coexpansion of the solvent and solute typically produces at most one 1:1 complex of each docking type in significant abundance. We are not aware of cases where two or more C-docking isomers on a single aromatic ring were explicitly observed and assigned for a methanol molecule. For O-docking in the particular case of furan, one can imagine that the in-plane (sp²) electron density is sufficiently differentiated from the aromatic (p) cloud density at the oxygen to stabilize two different conformations, but again we are not aware of the simultaneous detection of such isomers in a jet expansion. In contrast, the case of substituted anisoles has shown that the narrow barriers between C- and O-docking can be sufficiently pronounced to stabilize two spectroscopically identifiable species.⁶ If the limiting interconversion barriers between the two sites are shallow and narrow enough, one may approximate the ensemble of complexes in a jet expansion by a low (20–100 K) effective conformational temperature at which the interconversion freezes, together with always lower (≈10 K) rotational temperatures and uncertain as well as non-uniform (50–200 K) vibrational temperatures of individual normal modes. The lower the conformational freezing temperature is, the more sensitive the experimental switch between O- and C-docking will react to the energy difference of the competing sites. Therefore, such low barrier isomerism may serve as an experimental docking energy balance,⁶ which helps to decide which of the competing conformers is lower in energy, in favorable cases down to a sub-kJ/mol energy scale.⁷ 

These experimental C/O-docking ratios (at a very low but still unknown conformational temperature) may then be compared to quantum chemical predictions of the energy difference and help to judge the performance of the employed computational protocol. Several experimental and theoretical questions have to be addressed to render this relative energy comparison useful for benchmarking purposes.

On the experimental side:

• Can one be sure that there is no major barrier between docking conformations which raises the effective temperature at which interconversion freezes?

• Are statistical effects due to different capture areas in the docking process sufficiently relaxed when the complexes are probed?

• Is the probed area of the clustering expansion sufficiently uniform to allow for the determination of a representative average population?

• Does the result strongly depend on expansion conditions such as carrier gas and stagnation pressure, and if so, can one find reasonably stable parameter spaces for a meaningful comparison?

• Could there be band overlap between different conformations, even differing in C- or O-docking, thus potentially concealing isomers?

• Is the assignment to 1:1 complexes unambiguous and free of overlap with n:0 and 0:m as well as 2:1 and 1:2 features in the spectrum?

• Is the assignment of neighboring vibrational bands to C- or O-docking unambiguous?

• Are there sufficiently robust computed IR cross sections to convert signal strength in the linear spectra to relative abundance?

On the theoretical side:

• Is it acceptable to combine different levels of approximation for structure optimization, best electronic energy estimate, zero-point energy correction, and infrared (IR) cross section?

• Is a harmonic treatment of the zero-point energy correction sufficient?

• Are the isomer nuclear wavefunctions sufficiently localized in terms of the docking coordinates and the hydride stretching coordinate, or can one use robust anharmonic corrections to cope with partial delocalization?

• Was the conformational space sampled sufficiently well to identify all low-lying structures and their interconversion barriers?

• Is it meaningful to include thermal and entropic corrections given the non-equilibrium low internal energy nature of the experiment?

• Is it necessary to include symmetry corrections for complexes of different symmetry?

All these call for a concerted effort by theory and experiment and for many different furan and solvent derivatives to make sure or to disprove that such an intermolecular balance experiment can contribute to a rigorous (down to sub-kJ/mol) assessment of computational treatments for non-covalent interactions.

The present contribution forms the nucleation point for such a concerted effort, based on only three different furans [furan (Fu), 2-methylfuran (MFu), and 2,5-dimethylfuran (DMFu), Fig. 1] and one solvent [methanol (MeOH)]. For maximum independence of theory and experiment,⁸ a double-blind character of the study was implemented. The theoretical predictions were submitted before the experiments were carried out. Participants were given approximately three and a half months to carry out the necessary calculations. The data collection and treatment were managed by the Mata group. The experimental measurements were carried out after the theoretical submissions took place, but without any sharing of data. The latter was only compared once the experiments were completed. The result of this comparison is presented in this manuscript, together with a first critical analysis of the theoretical data, which will be seen to trigger future multi-experimental efforts.

A very recent independent matrix isolation study of methanol-furan complexes⁹ adds little to the scope of this benchmark study. The employed Ar matrix perturbs the weakly bound complexes and will certainly impact their relative abundance. The calculations used for the theory-experiment comparison in the aforementioned work should also be considered with caution. Density functional theory was applied without dispersion corrections. The latter are overall significant but particularly so for the description of the OH-π conformers. Unsurprisingly, the authors found a clear preference towards OH–O binding.

As for any experimental benchmark, it would be advantageous to have complementary experiments available to check for technical artefacts of the present FTIR jet characterization which might mislead the theoretical comparison. Powerful UV/IR techniques are difficult to apply due to ultrafast dynamics in the electronically excited states.¹⁰ Otherwise, they might provide information about the docking site from the shift of the electronic spectrum upon complexation. Microwave spectroscopy would be very attractive as it can provide direct structural evidence for the docking site and prove the coexistence of C-/O-docking in a supersonic jet,^11–15 but a quantification of the relative abundance is more challenging. Raman spectroscopy in supersonic jets could provide a higher relative sensitivity for the C-bonded complex and better spatial resolution of the concentration evolution,¹⁶ but the computed scattering cross sections may be somewhat less reliable than IR absorption cross sections.

This manuscript is structured as follows. In Secs. II and III, we describe the experimental setup and report on the spectra measured for the three systems. From this information, we derive an experimental estimate for the relative energetic ordering of the dimers. In Sec. IV, we introduce the theory submissions to the challenge, discussing the latter results in Secs. V–VII. In Sec. VIII, we compare the submitted theoretical fundamental OH band predictions to the experimental data and build an assignment. This, in turn, allows for a stricter comparison of the submitted theoretical results which we again review. In Sec. IX, we draw some initial observations on the challenge/study performed.

All spectra were obtained by probing a 600-mm long supersonic jet expansion through a 0.2-mm wide slit with FTIR spectroscopy in direct absorption.¹⁷ The mildly focused IR beam samples a >1 cm² cross section of the two-dimensional expansion and thus averages over early and later stages of the supersonic gas flow. The gas consists largely of He (99.996%, Linde) with small (≈0.1%) and variable admixtures of furan (Alfa Aesar, 99%, 250 ppm BHT), 2-methylfuran (Roth, ≥99%), or 2,5-dimethylfuran (Acros Organics, 99%), methanol (Sigma Aldrich, ≥99.8%) and/or methanol-d1 (eurisotop, 99% D, HDO + D₂O < 0.1%), and sometimes a relaxation promoter such as Ar (99.999%, Linde). Close to the nozzle exit, atomic collisions cool the molecules and additional molecule-molecule collisions lead to cluster formation. This aggregation process is usually limited to dimers and few trimers by sufficient dilution. Sufficient dilution is indeed the key issue in the present study as methanol prefers to cluster with itself rather than with an unsubstituted furan molecule such that rather low methanol concentrations are needed for the preparation of 1:1 complexes. Self-aggregation of the furans also has to be avoided because aggregates containing a second furan unit could influence the position and intensity of the 1:1 complex absorption in subtle ways. Furthermore, in the presence of major amounts of larger clusters, it cannot be ruled out completely that the 1:1 isomers show different propensities for further aggregation, thus distorting their observed abundance ratio. Therefore, most of the present work uses very high dilution, thus sacrificing precision for accuracy in terms of spectral intensity. It should however be remarked that extreme dilution may also change the aggregation and cooling mechanisms, reducing solvent-exchange mechanisms.^18,19 An important downside of high dilution for OH stretching bands is the distortion of weak band intensities by nearby, insufficiently compensated water rovibrational lines from the dried air in air gaps (1 bar) and in the vacuum spectrometer and detector chamber (<0.005 bar) along the IR beam path. This increases the error bar on the intensity ratio and renders OD stretching bands somewhat more reliable. An important experimental tool was the investigation of isotopic mixtures. These turned out to be useful for a direct and simultaneous comparison of relative intensities for isomers in the OH and OD stretching range under identical expansion conditions and to confirm the presence of only one methanol unit in the complex by comparison to mono-isotopic spectra. If two methanol units are contained in the cluster, the OH (OD) stretching fundamental wavenumber depends slightly on the isotopic composition of the other one, leading to small band shifts.²⁰ By using a 1:2 mixture, the reduced intensity of OD stretching bands was qualitatively compensated. A general disadvantage of isotopically mixed expansions is the statistically reduced intensity in the OH and OD spectral window.

Some experiments involved an admixture of 1%-5% Ar to the expansion with the aim of enhancing the most stable isomer by more effective barrier-crossing collisions. 5% Ar admixture (used for MFu) led to Ar condensation on the complexes, thus distorting the band positions and intensities. 1% Ar admixture (used for Fu) led to no significant relative intensity changes. This indicates in a qualitative way that the observed 1:1 complexes are of very similar energy, without providing information on which is the more stable one. Details of the experimental setup and the first measurement campaign are provided in the supplementary material (E1). In a second measurement campaign, a better compromise between precision (achieved at high concentration at the expense of warmer expansions and cluster interference) and accuracy of the spectral intensities (achieved at high dilution) was found for DMFu (supplementary material E2). The results were added to the data of the first campaign and are clearly superior to the previously published spectra.⁵ 

Figure 2 gives an overview of the key spectra obtained for the three compounds in the full OH, CH, and OD stretching range. The bottom spectrum involves a 2:1 mixture of methanol-d1:methanol and furan (Fu), both at very high dilution. It is dominated by the CH stretching modes of furan and methanol. There are hardly any methanol dimer signals, due to the very high dilution. Whereas the OH stretching region only shows methanol monomer and very weak 1:1 complex features, the OD stretching region is partly intruded by furan modes. For example, the band marked Fu at 2668 cm⁻¹ is a combination band of two CH scissoring modes of furan,²¹ but it does not overlap with two features of 1:1 complexes. The higher wavenumber band has a significantly higher intensity than the lower wavenumber band. To measure the OH signature of these 1:1 complexes with acceptable signal-to-noise ratio, the second-lowest trace shows a deuterium-free spectrum at similar total methanol concentration. It also shows a pattern of two 1:1 complex bands with the stronger one at a higher wavenumber. A qualitatively similar finding is observed for MFu (central four traces), where spectra for regular methanol (including 5% argon in the green trace), deuterated methanol, and a mixture are shown. Here, a significantly larger methanol concentration is chosen, as evidenced by stronger methanol dimer transitions. This gives rise to minor extra peaks partly due to mixed trimers, but the invariance of the peaks assigned to 1:1 complexes with respect to partial deuteration apart from residual water artefacts indicates that the intensities are trustworthy. The upper three spectra compare the previously analyzed isotope effect for DMFu at the highest concentration to a new spectrum based on a moderately diluted expansion with mixed isotope composition. A spectrum at very high dilution is also presented in the supplementary material (E1) and included in the evaluation. The previous assignment⁵ of mixed 1:1 complex bands is confirmed, as is the deuteration trend, but the relative intensity of the lower wavenumber complex in the OH stretching range varies with concentration. Two explanations (or a combination of both) are conceivable: either the higher concentration peak is affected by larger cluster absorptions or the lower concentration peak is affected by overlapping water rovibrational line residuals. The OD stretching range is obviously not affected by the latter problem.

The core OH and OD regions for the complexes of all three furans are repeated in Fig. 3 in magnified form. To analyze subtle anharmonic effects, the OD region is arranged on top of the OH region, after alignment of the monomer band centers and stretching of the OD stretching wavenumber scale by a factor of $2$⁠, the ideal harmonic isotope effect for an infinite oxygen mass. In a one-dimensional approximation, hydride stretching bands with a larger anharmonicity than the isolated methanol (such as methanol dimer²²) will thus be shifted to the right upon deuteration (positive slope of the dashed connecting lines), whereas bands with a smaller anharmonicity may even shift to the left, giving rise to a negative slope of their connecting lines.²³ A complementary interpretation assigns a higher librational stiffness to the OH group upon O-docking, compared to C-docking.²⁴ The resulting zero point vibrational energy (ZPVE) weakens the hydrogen bond, but less so upon deuteration, therefore increasing the downshift compared to C-docking. Other effects such as mode mixing could also play a role, in particular, for the OD region. However, it is striking that the stronger, higher wavenumber bands due to 1:1 complexes consistently show a positive slope, whereas the weaker, lower wavenumber components show a negative slope. Together with a very similar intensity pattern, no additional or missing bands upon deuteration, and a systematic trend in the respective band position with increasing methylation of the furan ring, this strongly suggests a uniform assignment of the band doublet across all three furans. Such a uniform assignment is thus assumed in the following arguments.

Without prejudice to any specific assignment, Table I summarizes all relevant experimental data. Here and in the following, the index h refers to the stronger, higher wavenumber band and l refers to the weaker, lower wavenumber band. One can see that the error bars for the OH intensity ratios are generally larger than those for the OD intensity ratios at high dilution. This is partly due to the water rovibrational line overlap discussed above and makes it difficult to judge whether the intensity ratio increases upon deuteration (hinting at an O-docking assignment for the h band), or decreases (hinting at a C-docking assignment for the h band), or stays the same (hinting at a common C- or O-docking assignment of both bands). However, there is a weakly significant decrease of this intensity ratio with increasing methylation of the furan ring, which indicates that the two isomers slightly approach each other in energy. The h band belongs to the more stable 1:1 complex. Both conclusions assume a sufficiently small interconversion barrier and sufficiently similar IR cross sections for all species.

In principle, the two bands could be due to two C- or O-docking complexes. The lack of significant deuteration effects on the intensity at very high dilution would be in agreement with this. However, in particular, the observation of two C-docking complexes would be inconsistent with previous findings for similar complexes.^4–6,14 Also, it would be difficult to explain why methylation of one or two C atoms in the ring would leave the pattern of bands so unaffected. Furthermore, Ar relaxation experiments without Ar coating carried out for Fu should have resulted in further relaxation. Lastly, it would be difficult to rationalize why the two C-docking isomers differ so much in OH stretching anharmonicity that one gives rise to a negative slope and the other to a positive slope in Fig. 3. Even their spectral splitting of up to 24 cm⁻¹ would be highly unusual for two C-docking methanol units.

Two O-docking complexes are somewhat more difficult to rule out, as the methanol can either approach the furan in its ring plane or perpendicular to it. Ring methylation could then have smoother effects, consistent with observation. Different anharmonicities are conceivable, but one would expect exactly the opposite—the more strongly downshifted band should correspond to a directed, strongly anharmonic hydrogen bond, whereas the compromise hydrogen bond with added dispersion interactions would be shifted less.⁵ This is inconsistent with the deuteration slopes in Fig. 3. Furthermore, it is not evident why Ar relaxation should have no effect on relaxing the more weakly bound complex with directed hydrogen bond to furan into a more compact structure. Finally, why should there be no C-docking isomers given that the furans offer a broad and attractive caption basin with four times as many C-atoms than O-atoms in the ring? Unless the theoretical calculations consistently predict two O-docking variants to outperform all C-docking complexes in energy, this is also an unlikely interpretation.

Assuming the stabilization of one C-docking isomer and one O-docking isomer in different abundance is more in line with findings for related systems. Unfortunately, deuteration does not give a clear indication about the position of the O-docking isomer via intensity changes. Regularly, this isomer gains in intensity upon deuteration because it is destabilized more by ZPVE. The only case with a significant deuteration gain of the stronger band is DMFu, whereas MFu and Fu remain inconclusive within error bars. Therefore, other evidence has to be collected. Assigning the more abundant, less shifted band to C-docking is rather unlikely if one looks at the slopes of the deuteration effects in Fig. 3. C-docking induces less diagonal and off-diagonal anharmonicity increase than O-docking^24,25 and should therefore behave more similar to methanol monomer than to methanol dimer. This is also experimentally observed for the more related benzofuran case;⁴ see Fig. 7 in Ref. 5. Therefore, the most plausible assignment for the presently investigated furans attributes the high wavenumber band to O-docking and the weaker satellite to C-docking. This assignment is fully consistent with the anharmonicity trends and with the previous assignment for the DMFu complex with methanol,⁵ but because of the inconclusive intensity effects upon deuteration, it is not beyond any doubt for MFu and Fu.

Unconventional very indirect evidence for the assignment of the minor isomer may come from the effect of substantial (≈5%) argon addition to the expansion in the case of MFu. It leaves the higher wavenumber band unaffected, but unexpectedly shifts part of the intensity of the lower wavenumber band from 3623 to 3627 cm⁻¹ (see Fig. 2, green trace). Although isomerization along the ring should not be ruled out, this could also be due to further complexation of the 1:1 complex with one or more Ar atoms.¹² Subject to computational confirmation, one may speculate that the preferred Ar complexation site in the C-docking isomer is opposite to the methanol coordination at the furan ring, where it competes for electron density, thus explaining a spectral upshift of the OH stretching band. However, this new kind of assignment aid requires future computational support to become useful for structural assignment.

If the uniform assignment of the major band h to O-docking and of the minor band l to C-docking is accepted for the moment, one may already hint without theoretical intensities that the O-docking isomer could be somewhat more stable than the C-docking isomer for all three furans, but not by a large margin. This margin is roughly bound by 1.5 kJ/mol, if a similar infrared cross section for both species and a conformational temperature of less than 100 K are assumed. The invariance of the spectra towards Ar admixture at least for Fu indicates that the true energy difference may be smaller. The question of symmetry number for symmetric furans (where C-docking is likely to lead to a chiral complex and therefore to enantiomeric degeneracy, whereas O-docking may keep a symmetry plane) does not change this estimate significantly.

The observed intensity trends suggest no major shift in the energetic preference between the two docking sites with increasing methylation. O-docking is found to be systematically favored by a small, perhaps slightly decreasing amount across the series. This is at variance with the simple theoretical predictions in Ref. 5, which imply a systematic shift from C- to O-docking with increasing methylation. Such disparate findings are in line with earlier experience—predicted OH wavenumber trends for hydrogen-bonded clusters at density functional theory (DFT) level are often more robust than delicate energy trends. However, it has to be seen whether these exploratory comparisons persist at the higher levels of theory triggered in this double-blind benchmarking study.

The key results of the experimental study are as follows:

• Clear evidence for two 1:1 complexes between methanol and each of the three furans.

• Rather uniform dominance of the less downshifted OH stretching band by about a factor of 2-3 at a low conformational temperature.

• Likely, but not unambiguous assignment of the higher wavenumber band to an O-docking complex (3654, 3645, and 3636 cm⁻¹ with increasing methylation).

• Likely, but not unambiguous assignment of the weaker and more downshifted band to a C-docking complex (3636, 3623, and 3612 cm⁻¹ with increasing methylation).

• Concentration-dependent and thus inconclusive deuteration effects on intensity ratios and more systematic deuteration effects on relative band positions.

• Small argon relaxation effects, further hinting at small energy differences between the docking sites.

In order to compare to theory, one needs comparative values from the experimental measurements. This is a rather delicate decision, particularly since one needs to define a theory-free criterion, while there is no absolute certainty about the identity of the isomers and no information about the relative intensity of the bands. Nevertheless, we believe there is enough evidence to follow the assignment provided above. It is estimated that the conformational temperature will be between 20 and 100 K, so one will take a 60 ± 40 K temperature value. The integrated bands will be proportional to the populations and the absorption cross sections. The energy difference between the two isomers is provided simply by

whereby $σlσh$ is the ratio in the absorption cross sections of the isomers, including also any potential chirality degeneracy. For now, a value of 1 will be taken, meaning that both isomers would have equal oscillator strengths. However, one could expect variations as large as 50% in the latter ratio. In this energy range, that would mean an added uncertainty of 0.2 kJ/mol. The remaining error bar is determined by the error in the integration of the measured bands and the uncertainty in the temperature, with the effect calculated approximately as

From these considerations, we define the final terms for the evaluation of the theoretical results. The energy difference between the two dimers should be 0.6 ± 0.8, 0.5 ± 0.8, and 0.4 ± 0.5 kJ/mol for Fu/MeOH, MFu/MeOH, and DMFu/MeOH, respectively (OH–O dimer slightly favored). This means that although the OH–O dimer is very likely the most stable, the uncertainties in the temperature and the integration of the band ratios do not allow one to completely exclude the possibility of an OH-π global minimum (that would be within the error bar). Note that this uncertainty could possibly be reduced by 0.2 kJ/mol as soon as an estimate for $σlσh$ is available. In the following, the lower wavenumber peak assigned to the OH-π complex will be set to zero energy. The experimental prediction is then that an OH–O dimer will be found in the energy range of −0.6 ± 0.8, −0.5 ± 0.8, and −0.4 ± 0.5 kJ/mol, with increasing methylation.

There were 12 submissions for the challenge in total. The groups involved as well as the theoretical methods applied are listed in Table II. The theoretical approaches are split into three fields: optimization, zero-point vibrational correction, and energy (the electronic energy used for the final estimate). Quite fortunately, there is a healthy mix of different approaches, from hybrid DFT up to wave-function-based explicitly correlated methods, while keeping some overlap. This overlap will facilitate some of the analysis we will be carrying out on top of the submitted results. The Baptista group provided three different estimates, based on CBS-QB3,^26,27 ωB97xD/6-311G(d,p), and ωB97xD/6-311G(df,pd).²⁸ Given the similarity between the last two methods and corresponding results, we picked the submission with the largest basis set. The values for the 6-311G(d,p) basis set will not be discussed.

From Table II, one can observe a predominance of hybrid dispersion-corrected DFT²⁹ for the optimization of the structures. The only exceptions are entries F (SCS(1.1;0.66)-MP2),³⁰ I (spin-component scaled second-order coupled-cluster method SCS-CC2),³¹ J (coupled-cluster single and doubles with perturbative triple excitations CCSD(T)), and the entries K and L whereby a double hybrid was applied. There is only one entry featuring a non-dispersion-corrected optimization procedure, B3LYP³² in the case of entry A. This is part of the CBS-QB3 procedure, which was established for benchmarks with relatively small molecules. For such cases, dispersion effects have little impact on the accuracy.

Only three of the entries include some form of anharmonic corrections to the ZPVE – C, F, and L. Entries C and F provide an anharmonic estimate considering all vibrations and L corrects for a selected dimension. Otherwise, all results are based on the rigid-rotor harmonic oscillator standard model. Finally, there is a limited variance in the electronic structure methods used for the energy refinement, with a majority of groups attempting to reach the CCSD(T) basis set limit, either through the use of large basis sets, extrapolation procedures, or explicitly correlated variants of the method. Several entries made use of the local pair natural orbitals (DLPNO) approximation for the coupled-cluster calculations.^33–36 The impact of the latter [relative to full canonical CCSD(T) calculations] is not assessed in this work. The interested reader should refer to Refs. 35 and 37. Only two entries (G and H) feature symmetry-adapted perturbation theory (SAPT) calculations,^38–40 while submission B features a final energy at the DFT level.

Making use of the criteria defined in Sec. III, there are several entries which fit into the estimated energy gap between OH–O and OH-π coordination. These are marked in bold in Table II. There are six entries which correctly predict the value for all three systems, entries C, D, I, J, K, and L. It should be noted, however, that several other entries are only a few tenths of a kJ/mol away from the mark.

These considerations constitute a first rough comparison between theory and experiment. Note that the relative energies stem from comparing one OH–O and one OH-π structure, as optimized by each participating group. So far, the structure of the latter has not been discussed. It remains open if these correspond to the minima observed in the IR spectra. In order to further extend our analysis, we need to consider the optimized geometries, judging the sampling made in each entry, the different energy contributions which lead to the results in Table II and finally a comparison of the IR spectral signals with theoretical results in order to provide an assignment. The theoretical data will also be used later to reduce the experimental error bar.

First of all, one should identify which minima were in fact sampled in each submission. We examined the individual geometries in order to identify common patterns. Most submissions provided three structures per system, two minima with an OH–O binding motif and one OH-π structure. There were a few exceptions, but generally all other structures reported were significantly higher in energy, and we will leave the latter out of this discussion. Several submissions were found to report on different enantiomers. This fact is irrelevant to the present study, and in order to facilitate the visual comparison among all entries, we generated the same enantiomer for every minimum by reflecting the atoms relative to an approximate ring plane.

The two types of OH–O binding dimers could be classified as “on-top” (OH–O^t) and “in-plane” (OH–O^p). Most OH-π dimers were rather similar and require no separate classification. In Fig. 4, we present an overlap of the structures from all submissions which could fit into the three mentioned binding patterns. In all cases, a root-mean-square deviation (RMSD) fit of the O, C2, and C5 atoms of the furan ring was used to overlap the geometries.

As one can observe from Fig. 4, the OH–O^t minima are stabilized by both the OH–O hydrogen-bond and the contact (mostly through the methyl group) to the ring plane. In the case of the OH–O^p minima, the OH–O contact is dominating and the OH-π dimers have the methanol lying over the plane, favoring again several interactions with the π-cloud of the furan derivative.

Overall, there is a strong overlap among all structures with only a few exceptions. The A entry structures for the OH–O^t dimers of MFu and DMFu are somewhat displaced relative to the other submissions. The Fu complex is also missing. The latter is due to convergence problems in the B3LYP optimization. As mentioned before, the CBS-QB3 composite scheme makes use of a DFT optimization without dispersion corrections. This in turn disfavors the contact between the methyl and the furan plane and leads to structures with an alignment closer to the OH–O^p conformer. Another difference is found in the MFu OH-π structures from G and H. While all other groups found a minimum with the methanol pointing towards C3, the G and H structures were directed towards C4. In the case of Fu and DMFu, the two carbons are symmetric and there can only be one conformer. Finally, in the case of Fu, the OH–O^p minima show a large variance in the angle to the furan O. Relaxed surface scans keeping the OH–O binding in the same plane while varying the O–H–O angle reveal an extremely shallow surface (see supplementary material T2). A transition state is observed when the methanol aligns with the $σv′$ symmetry plane of the furan, but with a very small barrier. It is then expectable that the different methods used in the optimization of the dimer will find slightly differing positions for methanol. The submitted geometries G and H correspond, in fact, to such a transition state. These were nonetheless kept in the case of Fu, given the very small energy difference.

We carried out one further test to the geometries. The electronic energies were recomputed at the CCSD(T)-F12/cc-pVTZ-F12 level,⁴¹ with the same options and program version (Molpro2012.1).⁴² This is only a rough assessment but does allow us to compare all the structures on the same footing. It also helps to understand the impact that the geometry optimization step had on the final reported estimates. For each dimer type (OH–O^t, OH–O^p, and OH-π), we find the lowest lying structure and plot the relative energies to this reference for all other submissions. The results are shown in Fig. 5. Ideally, the errors within each entry should be similar, meaning a constant error relative to the structures closest to the CCSD(T)-F12/cc-pVTZ-F12 minimum.

We start by discussing the values for Fu/MeOH. As one can observe, most of the submitted structures are very close in energy. The only exceptions are the B structures, 1.2-1.4 kJ/mol above the lowest value. However, this is a rather constant error and should have little impact in the relative stability of the dimers. A similar trend can be observed for the other submissions with hybrid DFT optimization (entries D, G, and H). The geometries show small variations relative to more costly (and in principle accurate) theoretical methods. The largest relative deviation is found for the OH–O^p dimer of entry C. The relative error (compared to the other dimers) is as large as 0.6-1.0 kJ/mol.

Considering the full set of results, the relative errors of entry A confirms the observations made from Fig. 4. The dimers which benefit most from dispersion interactions (OH–O^t and OH-π) reveal large errors. Otherwise there are no major outliers. The lowest lying structures are provided by the SCS(1.1;0.66)-MP2/aug-cc-pVQZ (entry F), CCSD(T)/aug-cc-pVTZ (entry J), and B2PLYP-D3/6-311++G(d,p) (entry K). Several submissions present constant deviations such that the impact of the geometry optimization method on the energy ordering should be rather small.

With all the minima grouped into the three different binding motifs presented in Sec. V, it is possible to compare the submitted results on a better footing. Each submission had a different presentation of the results, some providing a dimerization energy, others providing the raw data. The best way to compare, and one which was possible in all cases, was to take a reference geometry and plot relative energy differences. For all systems, we chose the OH-π dimer as the reference. The results (including all corrections as described for each individual submission in the supplementary material T1) are shown in Fig. 6.

We start by discussing the furan (Fu/MeOH) results. Although this is the smallest system, and theoretically the easiest to handle, the shallow potentials surrounding the different minima as well as the small difference between the stationary points leads to a rather diverse picture. There are no results available for A since only one structure was available. Submissions B, D, F, G, and H predict that the OH-π dimer will be the most stable structure. C and K predict OH–O^p to be the global minimum. There were no results available for this dimer in E, J, and L so that the energy gap between OH–O and the OH-π refers to the “on-top” structure. One already observes large disagreement between the sets which actually fulfilled all criteria (C, D, I, J, K, and L) when comparing to experiment. Although the gap was within experimental error bars, the values referred to different minima.

In the case of MFu/MeOH, there is a much more consistent picture, just as with DMFu/MeOH. With exception of entries A, C, G, and H, all predict the OH–O^t dimer to be the most stable structure. Entry C quite consistently gives OH–O^p as the lowest minimum across the series. This already tells us that the simple evaluation which was carried out before is not enough. There is a disagreement about which of the OH–O dimers is the most stable and to which the gap should be calculated. We investigate this further in Secs. VII and VIII.

From experiment, a very small energy difference is predicted between the OH–O and the OH-π dimer. In fact, the energy is so small that the vibrational energy corrections can easily tip the scales in one or the other direction. The electronic energy contributions to the relative energies are detailed in Table III, together with the contributions from the ZPVE. Overall, some consistency can be observed across the different entries for the relative electronic energies between the OH–O^t and OH-π dimers (exception made to entries A, G, and H). The picture becomes more complicated for the OH–O^p conformer, where even closely related approaches such as those from C and E give differences as large as 2 kJ/mol for E_el. This may in part be attributed to the difficulties found in optimizing the latter minimum, as visible in the larger spread in Fig. 4.

Observing now the E_vib terms, it becomes more difficult to see any type of pattern. There are a few entries which agree quite well over the whole range (e.g., E, J, K, and L), a positive sign of reproducibility. But considering the full picture, there are several deviations which are hard to explain. This is particularly true for entries C and F which include anharmonic E_vib corrections. C predicts a much lower ZPVE in the OH–O^p dimers than in the other conformers. All other entries show only small variations between the two OH–O motifs.

In order to evaluate in further detail the impact of the ZPVE corrections, we have plotted in Fig. 7 the electronic energy differences, together with the ZPVE corrected values from the harmonic calculation plus the anharmonic corrected results for entries C, F, and L. In L, a hindered rotor approach was used, considering the rotation of the methanol with respect to the ring. The latter builds a path between the two minima on top of the ring (which are separated by low barriers) on the potential energy surface. In entry C, all vibrations were included in the anharmonic treatment, but aside the O–H stretch, only diagonal corrections were included. In entry F, the ratio of the vibrational second-order perturbation theory (VPT2) anharmonic and harmonic ZPVEs at the SVWN level of theory was computed and used to scale the SCS(1.1; 0.66)-MP2 harmonic ZPVE corrections.

There is a very stark contrast between the different sets. The electronic energies are only in relatively close agreement for MFu and DMFu, with the OH–O^t dimer being preferred by all sets. The largest difference among the entries is in the result of the anharmonic corrections. In all cases, entry C concludes that the OH–O^p dimer is strongly favored relative to the OH–O^t. According to the authors of the study, the differences in the ZPVE between the OH–O^t and the OH–O^p dimers are mainly due to changes in the potential of the methanol methyl rotation. Nevertheless, the impact of the correction is somewhat surprising. Entries F and L show much smaller variations. This is determining in the final energetic ordering.

Only two entries provided anharmonic values for the OH-stretch, C and F. The computed fundamental frequencies are provided in Table IV. In the case of F, a one-dimensional calculation is performed so that only diagonal anharmonicity effects are included.

Comparing the computed anharmonic values of the OH-π dimers to the tentatively assigned experimental frequencies, the agreement is striking. All computed values are within 10 cm⁻¹ of experiment, and thus the methylation trend is reproduced. The assignment of the higher wavenumber to the OH–O dimer as well as the lower band to the OH-π dimer is confirmed. Less clear is the assignment of the actual OH–O conformer. Results from entry C estimate that the two possible dimers would overlap in the case of Fu/MeOH. In the other two cases, the OH–O^p dimer would better fit the data, with deviations again within 10 cm⁻¹. In the case of entry F, the comparison is difficult, except for the MFu case where again the OH–O^p dimer agrees best with the measured value.

With the available theoretical data, in particular the absorption intensities, one can revisit the energy criteria which were initially proposed from experiment. The problem remains that there is no certain assignment of the OH–O bands across the furan series. To avoid any bias, one can consider the computed $σOH−πσOH–O$ for both options and build compound estimates. We took the anharmonic absorption intensities provided by entry C, removed the 0.2 kJ/mol associated uncertainty, and combined the error estimates from experiment with the span provided by including these quantities in Eqs. (1) and (2). This leads to no change in the case of MFu but slightly smaller energy spans in the other cases. The revised values are −0.6 ± 0.6, −0.5 ± 0.8, and −0.5 ± 0.3 kJ/mol, for Fu, MFu, and DMFu, respectively. The new criteria do not change our analysis of the theoretical results (Table II) greatly, but the only entries consistent for all three systems are reduced to C, I, J, K, and L.

The benchmark could be further refined by experimental confirmation of the OH–O dimer identity, for example, through a microwave spectroscopy study of the three systems. If the OH–O^p conformation is found to be the lowest energy structure, this would certainly highlight the importance of anharmonic ZPVE corrections. It would also stand in disagreement with the original harmonic B3LYP-D3-based assignment from Ref. 5. If the opposite were to be found, meaning that the OH–O^t conformer is the most stable, submissions I, J, and L would be the most consistent set of results available. Methods like D, E, and F basically would agree with the latter entries but differ significantly in their ZPVE corrections. It is also possible that the preference changes across the methylation series, but in that case, the discussion would be much more intricate.

One further criterion for evaluating the different submissions which has not been discussed is the question of where the second OH–O dimer should be found on the energy scale. Only two frequencies were observed in the experiments, the latter had been correctly assigned to two different binding motifs. In the case of Fu, and observing Table IV, it is possible that the OH–O^t and the OH–O^p dimer bands overlap. However, for the methylated furans, the second OH–O dimer is not observed. One could argue that a correct set of results should place the latter higher in energy than the remaining two structures. However, it is also possible that relaxation processes are responsible for the disappearance of the band since the barrier between OH–O^p and OH–O^t should be relatively small. For these reasons, we refrain from using this as a disqualification criterion, as long as no reliable barrier calculations are available.

There are several conclusions to take out of the challenge. We start by the most obvious lessons. Some of the composite methods suggested in the past need to be revised, with the underlying optimization procedures brought up to date. The CBS-QB3 structures showed large deviations from the remaining sets, more than likely due to the lack of dispersion corrections in the optimization step. This brings us to the second, rather clear observation, that the dispersion-corrected density functionals show a remarkable performance in structure optimizations. These include B3LYP-D3, ωB97xD, and PBE0-D3 which have much lower computational costs than other alternatives but only a marginal effect on the relative energies (see Fig. 5). This, of course, only holds if the final energies are computed at a suitable theory level. Most submissions applied CCSD(T) or approximations/variations thereof. Whereas there is a satisfactory convergence for relative electronic energies using such CCSD(T)-based methods and also for harmonic ZPVE corrections, attempts to treat the ZPVE beyond the harmonic approximation lead to surprisingly large deviations. This either points at important and challenging anharmonic effects for the investigated systems or else may indicate a limited robustness of such attempts.

Vibrational spectroscopy alone was not enough to characterize the OH–O dimer. Considering the uncertainty in the conformational temperature and the integration of the measured absorption bands ratio, it was not even absolutely clear whether an OH–O dimer would be more stable than an OH-π structure. However, it did provide from the start a rather strict test. Across the series, the energy difference between OH–O and OH-π dimers should be lower than 1.5 kJ/mol. Most of the entries combining at least hybrid functionals with dispersion corrections for the optimization and a CCSD(T) energy estimate were able to meet this criterion. Even when crossing anharmonic vibrational frequencies with the measured data, the structure assignment was not totally clear. However, it was possible to further reduce the energy window. The criteria could be further tightened if the correct OH–O conformer were to be identified, e.g., by microwave spectroscopy.¹³ Experimentally, it should also be possible (given enough time) to further reduce the error bar in the measured intensities.

One last remark should be made. Although the challenge was double-blind in character, the fact that results were already published for DMFu might have influenced some of the submissions. The three binding motifs which can be identified in the structures were already featured in the first paper. Only in the case of entries G and H did one observe a larger variety in the dimer structures, even if several of them were too high in energy to be considered. For future challenges, it would be of interest to reduce even further the reference points, increasing the significance of structure sampling. In view of the specific properties of jet expansions, this sampling should also include interconversion barriers.

See supplementary material for a detailed description of the original experimental part (E1), of a DMFu remeasurement (E2), and of the laboratory journal (E3) as open data repositories. Computational details are also available, detailing the methods and programs used in each submission (T1). Furthermore, an additional surface scan of the Fu/MeOH system is available (T2) as well as the surface scans used for the hindered rotor calculations of entry L (T3).

Several of the participating groups would like to acknowledge funding support from the SPP 1807 “Control of London dispersion interaction in molecular chemistry.” R.D., H.M., M.M.A.M., and M.H. would like to extend their sincere appreciation to the Deanship of Scientific Research at King Saud University for funding through the Research Group Project No. RGP-333 as well as the support of COST ACTION CM1405 MOLIM. L.B. and M.N.P. would like to acknowledge the financial support of FAPERJ and the infrastructure support of Centro Nacional de Processamento de Alto Desempenho em São Paulo (CENAPAD-SP). K.L. and W.A.K. acknowledge the Cluster of Excellence “Tailor-Made Fuels from Biomass” (EXC 236), which is funded by the Excellence Initiative by the German federal and state governments to promote science and research at German universities. L.C.K. thanks the Deutsche Forschungsgemeinschaft for support within the SFB 985 “Functional microgels and microgel systems.” D.F. is thankful for financial support from the Deutsche Forschungsgemeinschaft (German Research Association) through Grant No. GSC 111. D.M.B. and I.S.U. acknowledge the Viper High Performance Computing facility of the University of Hull and its support team.