Boron-substituted aromatics - part two: thermochemical calculations

In my last post I described how I calculated vibrationally resolved UV/vis absorption spectra for a couple of boronaromatics and their acetonitrile (ACN) adducts and how this helped to resolve the mystery of their odd solvatochromism (solvent dependent absorption spectra). The one question that remained open was: Do DBI and DBA exclusively form 1:1 monoadducts with ACN or also the 1:2 diadducts? Since the latter are predicted to have essentially no absorption in the relevant spectral range, i.e., between 300 and 700 nm, I could neither confirm nor exclude their formation.

To answer this question, I started to conduct what was planned as a "brief" thermochemical analysis of formation of the respective mono- and diadducts. However, since the first results were inconclusive the I ended up employing a hierarchy of methods of increasing sophistication (and computational demand), which only eventually (at the very highest level of theory) cumulated in a good enough agreement with the experimental observations. I think the convergence of the results with respect to level of theory and employed basis-sets is quite instructive (without having much experience in quantum-thermochemistry), but at the same time too long and theoretical to be included in the article, and hence provides a nice topic for this blog. Lets get to it!

The initial plan was to just collect and investigate the results from the calculations I had already conducted to model the vibrationally resolved spectra. These include all optimizations and normal-mode analyses for naked DBI, DBA and DBP as well as their monoadducts and trans-diadducts, such that solely the respective calculations for ACN were missing. Having completed the latter and putting together all the numbers (computed at the same level of theory as the vibronic spectra, i.e., B3LYP-D3BJ/SVP) afforded the following picture:

Thermochemical data for the formation of the ACN mono- and diadduct of DBI, DBA and DBP (only trans-diadduct of DBI is shown) calculated at the B3LYP-D3BJ/SVP level of theory (using the COSMO solvation model to mimic the influence of the solvent ACN here and in all of the following results). The leftmost value termed "electronic" is computed from the electronic energies only. "ZPE" additionally includes the energy changes due to zero-point vibrations, which make the formation of additional B-N lewis-bonds more expensive compared to electronic energies only (since some energy goes into the zero-point vibrations of these new bonds). "RT" adds thermal corrections for room-temperature on top of that, which is the energy is is stored in the low-frequency (torsional and collective) modes that are thermally excited already at room temperature. ZPE and RT corrections are computed from the result of the normal-mode analyses of the molecules.

The good news is that at least the ordering suits the experimental obervsations with the DBI monoadduct attaining the lowest energy, followed by DBA and eventually DBP. The bad news is that all values are about 20 kJ/mol too low to explain the experiment. The negative values of the free energies of formation ("+RT" values) for all mono-adducts suggest a quantitative formation for all of them, even if only traces of ACN are present. Remembering the experimental observations (no changes in the spectrum of DBP even in pure ACN, weak influence for DBA at high ACN concentrations, major changes only for DBI), the values are apparently systematically too low.
Considering that thermochemical calculations are (in contrast to the excited-state calculations) quite sensitive to basis-set size and that I have employed the small SVP basis, this should not really be a surprise. The magnitude of the effect, however, did surprise me. The underlying problem is that small incomplete basis sets are susceptible to the so-called basis-set superposition error (BSSE). In a nutshell, BSSE is the result of an artificially constrained electronic wavefunction. Consequently, any calculation for a larger system (e.g. adducts, dimers) attain energies that are systematically lowered compared to respective calculation for the respective subs systems (their isolated constituents, mononers), since in the supermolecular calculation the wavefunction of a fragment A can also use the basis-functions on fragment B and vice versa.
Although there are corrections for the BSSE, the straightforward approach is typically to use a larger (augmented) basis set, which is what I did. To save time in the calculations with the larger (def2-TZVP) basis set, I skipped the optimization of the the geometries and used the ones from before, which is a quite common thing to do since geometries are not as sensitive as energies to basis-set size. To save even more time I used another common trick: I did not redo the normal-mode analyses for the ZPE and RT corrections but used the ones obtained with SVP, since they are also much less basis-set dependent. Such a combined approach is abbreviated B3LYP-D3BJ/def2-TZVP//SVP (behind  the double-slash follows the level used for the geometries/corrections).

Electronic energies obtained at the B3LYP-D3BJ/def2-TZVP//SVP level of theory.

With the larger basis set, the results are already in much better agreement with the experimental observations. With the main difference being a systematic shift to lower energies of formation by about 25 kJ/mol, also the energy of the DBP adduct has moved up slightly relative to the DBI and DBA adducts. Although the values do not yet fully agree with the experiment, we can at this point most certainly answer the initial question about the formation of diadduct: Judging from these results, their formation can most certainly be excluded. The values for the DBA and DBP diaddcuts (not shown) are very similar to those of DBI with all free energies of formation well above +30 kJ/mole.

However, I wasn't quite satisfied with the agreement and more importantly just curious how ab-initio methods would compare to the DFT results. So I eventually conducted additional SCS-MP2 and CEPA/1 (a coupled-cluster variant) calculations. To also improve on the structures, I re-optimized them at the SCS-MP2/SVP level of theory and later refined them at the SCS-MP2/def2-TZVP level of theory (this was JUST possible for the DBI diadduct). For the final energies at the SCS-MP2 and CEPA/1 levels of theory, I even conducted complete basis-set (CBS) extrapolations using the def2-SVP and def2-TZVP sets. This technique estimates the energy that would be obtained with a hypothetical, complete set of basis functions from the differences between the energies of two limited sets.  The free energies of formation for the monoadducts of DBI DBA and DBP are summarized in this final figure:

Summary of the free energies of formation of the monoadducts at (from left to right) increasing levels of theory.

SCS-MP2 with the small SVP basis yields an even stronger overbinding than B3LYP, most certainly for the very same reasons (mainly BSSE). This is not too surprising since correlated methods are said to be in general more prone to basis-set incompleteness errors than DFT.
The systematic overbinding is improved but not eliminated at the mixed approach with the larger basis (def2-TZVP//SVP), and becomes slightly worse at the fully consistent SCS-MP2/def2-TZVP level of theory. Note that the difference between the mixed SVP//def2-TZVP and fully consistent def2-TZVP approaches is quite small, whereas the mixed approach is MUCH cheaper.
The CBS extrapolation, which I would have expected to correct the systematic errors, actually worsens the agreement with the experiment. All energies of formation are again systematically reduced, such that even the DBP monoadduct is suggested to be stable at this level. I had heard that MP2 tends to overbind organic molecules, but this is a larger error than I expected. In particular since I have used the spin-component scaled variant, which should be less prone to these problems.
Only at the CEPA/1/CBS level of theory do the calculated free energies of formation ultimately agree with the observed behaviour: DBI showing quantitative monoadduct formation is weakly bound, DBA showing weak but significant monoadduct formation at high MeCN concentrations is energy-neutral, and DBP, which does not show any signs of adduct formation, is strongly endothermic. To make the coupled-cluster calculations possible for these already quite large systems, I employed the very handy domain-localized pair natural orbital (DLPNO) approximation implemented in Orca, and still burned a lot of computer time.

At the time, I found it quite instructive (and ultimately satisfying) to see how these numbers eventually converge to agree with the experimentally observed behaviour at the highest level of theory, and how all of the methods show the issues I had heard about, but never observed before. I hope you has a similar experience and can take something home. If you have any questions please leave a comment.

The next post will be about the investigation of the mechanism underlying the different fluorescence quantum yields of DBI, DBA and DBP.

So long!
Jan