Choosing a solvation model for simulating reactions in KOH(KOBut)/DMSO superbasic media

The ability of protocols based on the B2PLYPD/6-311+G**//B3LYP/6-31+G* method with various schemes for accounting for nonspecific solvation to reproduce C-H, N-H, O-H and S-H acidity in a dimethyl sulfoxide medium is considered. For a selected set of 20 compounds, typical reagents for reactions in superbasic media, the IEFPCM scheme with UFF cavity and α = 1.35 multiplier yields better results than the popular SMD model.


Introduction
Modern advances in the chemistry of acetylene are largely associated with the use of media with increased basicity (superbasic media) of the type "hydroxide or alkoxide of an alkali metaldimethyl sulfoxide" [1]. The use of such media made it possible to significantly expand the possibilities of classical acetylene reactions such as ethynylation, vinylation, and acetylene-allene rearrangement [2], and also led to the discovery of a number of new reactions. In particular, the recently discovered reaction of ketone vinylation formed the basis for an impressive variety of cascade transformations ensuring assemblies of complex heterocyclic systems [3].
Modeling such reactions by means of quantum chemistry contributes to a better understanding of the laws governing their course. Adequate models, methods, and protocols are required to reliably describe the mechanisms of chemical transformations in a superbasic environment.
We have recently proposed several models of the reaction center of the superbasic system [4]. The most complete of them considers an alkali metal cation surrounded by up to five explicitly included solvent molecules. The remaining part of the bulk solvent is taken into account within the polarizable continuum model. At the same time, it was shown that simpler models, namely the "monosolvate" one, which includes explicitly only one solvent molecule, and even the so-called "anionic", which does not include the metal cation at all and takes into account only the effects of nonspecific solvation, are capable of providing thermodynamic and kinetic characteristics of the reactions under study that are close to those obtained using a more complete model.
When choosing the level of theory, we settled on the B2PLYP functional [5], which provides good performance with relatively low resource intensity. It is noteworthy that this functional successfully copes with the problems arising in the description of the reactions of acetylenes and ketones in the presence of bases, which are traditionally difficult for most popular functionals, as well as for the ab initio MP2 approach, such as acetylene-allene rearrangement [6] or the aldol reaction of acetone [7].
The subject of this publication is the selection of a design scheme for describing the effect of a solvent. Over the years, we have used the polarizable continuum model (PCM) using the integral equation formalism variant (IEFPCM) proposed by Tomasi et al. [8], limiting our consideration to electrostatic contributions only. At the same time, a different version of IEFPCM proposed by Truhlar et al. [9] and known as SMD (Solvation Model based on the quantum mechanical charge Density), has become widespread. In addition to the electrostatic term, the SMD scheme includes cavitation, dispersion and repulsion energies. In the current version of Gaussian program suite, this scheme is the recommended choice for computing ΔG of solvation.
The SMD model has presented a good performance for challenging neutral molecules in water solvent. The advantages of this model for describing anions in non-aqueous solvents are not so obvious. Thus, Pliego and Riveros [10] found that for six nucleophilic substitution reactions involving different anions, the SMD model is not quite suitable for dipolar aprotic solvents, giving a root mean square error (RMS) of 5.6 kcal/mol, while the PCM method with atomic cavities, suggested by Pliego and Riveros [11], gives an RMS of only 3.2 kcal/mol.
A good test for the method's ability to describe both neutral and anionic systems in solution is the calculation of C-H, N-H, O-H and S-H acidities. Here we present the acidity estimates of compounds that are typical reagents, intermediates and products of the superbase-promoted reactions obtained using the classical PCM and SMD approaches.

Computational details
According to Pliego et al. [12], the calculation of the pKa value in DMSO solution can be undertaken through the use of the following proton-transfer reaction: HA + OH -→ A -+ H2O (1) This approach is more adequate than using the direct ionization of the HA acid, because it does not require the value of the experimental solvation free energy of the H + ion [12]. Considering the general chemical equilibrium relationship one can obtain pKa(HA) = pKa(H2O) + ΔG/2.303RT , (2) where ΔG is given by The Gibbs free energy in the solution is renormalized to a concentration of 1 mol/L, taking into account the entropy losses during solvation.
and that the entropy in the DMSO solution Ssol was obtained from the entropy for ideal gas in the harmonic approximation Sharm as follows [6]: We have also used the experimental pKa value of 31.4 for water reported by Bordwell [13]. Two approaches were used to estimate the free energies in the gas phase. One of them, considered as a reference, was the composite CBS-Q//B3 precision method [14,15] implemented in the Gaussian 09 program suite [16].
Another approach, less resource-intensive and, therefore, more practical for studying the mechanisms of complex reactions, is based on the double-hybrid B2PLYP-D functional, including the dispersion correction [17]. In this approach, the structural parameters of the studied molecules and their anions were optimized in the gas phase using density functional theory (DFT) at the B3LYP [18,19] level of theory with the 6-31+G* basis set. The vibrational corrections to enthalpies and Gibbs free energies were calculated at the same level of theory (B3LYP/6-31+G*) at a standard temperature 298.15 K. Furthermore, the energies at the stationary points were refined by using the B2PLYP-D functional in combination with the extended 6-311+G** basis set.
Calculations of solvation free energies within the SMD model were carried out at the B3LYP/6-311+G** level of the theory with radii and non-electrostatic terms for Truhlar and coworkers' SMD

Results
The following types of acidity are of interest for describing the reactions in superbasic media: i) O-H acidities, primarily O-H acidities of alcohols (their assessment is extremely important for describing the comparative strength of superbases based on alkali metal hydroxides and alkoxides, as well as for describing the classical reactions of vinylation of alcohols with acetylene and its derivatives) and oximes, intermediates of cascade transformations with the participation of hydroxylamine and its derivatives; ii) C-H acidities of acetylene and its derivatives (associated with ethynylation reactions), ketones (vinylation reaction of ketones and further transformations of the resulting vinyl-and allyl ketones) and ketimines (of particular interest in connection with studies of the possibility of their vinylation with acetylenes); N-H acidities of amines, pyrroles (primarily due to the unique reaction of the formation of vinylpyrroles) and imines (unfortunately, for the latter, there are very few experimental data available) and iv) S-H acidities of sulfides (starting materials on the way to divinylsulfides).

Gas phase acidities
Calculated within the framework of two approaches, CBS-Q//B3 and B2PLYP-D/6-311+G**//B3LYP/6-31+G, the Gibbs free energies of proton detachment are collected in table 1. The Gibbs free energy of the proton in gas phase was assumed to be -6.29 kcal/mol [21]. Comparison of theoretical estimates with experimental ones shows that the CBS-Q//B3 method provides, as expected, "chemical" accuracy, and the mean absolute deviation (MAD) for it is close to 1 kcal/mol. The B2PLYP-D/6-311+G**//B3LYP/6-31+G* approach yields a slightly larger MAD, 1.44 kcal/mol. It can be seen, however, that the largest error exceeding 3 kcal/mol is contributed by CF3CH2OH and 2indanone. Elimination of these two compounds reduces the MAD to 1.12 kcal/mol. Unexpectedly, the CBS-Q//B3 precision approach predicts a large enough error for the abstraction of a proton from a water molecule. Such an error can cause a systematic underestimation of the pKa values calculated by Equations 2-5 by 1.2 units.

CBS-Q//B3 based acidities
To calculate DMSO acidities, we used several schemes for accounting for the effects of solvation: i) the B3LYP/6-311+G** SMD model (SMD); ii) the B3LYP/6-311+G** PCM with non-electrostatic terms included and scaling factor α = 1.35 (PCM-LB); iii) same as PCM-LB, but taking into account electrostatic term only, as recommended in Ref. 12 (PCM-ES). The results of these calculations collected in table 2 together with the experimental data show that the SMD scheme gives significantly worse results (MAD = 3.25 pKa units) than the classical PCM approach (MAD = 2.26 and 2.18 pKa units for PCM-LB and PCM-ES, correspondingly). The largest error, reaching 7.7 and even 7.8 pKa units in the calculation with the SMD method, is associated with the C-H acidities of phenylacetylene and cyclopentadiene. After excluding these points, the MAD decreases to 2.6 units, but still remains greater than the error of the PCM-LB and PCM-ES schemes. Table 2. CBS-Q//B3 based acidities obtained using various schemes (see text for notation) compared to experimental [29] ones. The values in parentheses were obtained using the experimental value of the water deprotonation energy in the gas phase.  The values of the mean error (ME) show that all methods are characterized by an underestimation of the pKa., which is largely due to the above error in the gas phase acidity of water. The use of the experimental value for the water deprotonation energy brings a significant improvement in the results for the PCM-LB and PCM-ES schemes, while the error of the SMD scheme is still large, and even after elimination of the problem values for phenylacetylene and cyclopentadiene the MAE of the SMD scheme is 0.4 units higher than that of the PCM circuit. Thus, the PCM scheme with a scaling factor α = 1.35 seems to be more preferable for describing both neutral molecules and their anions in dimethyl sulfoxide.
The MAD value obtained for the PCM-LB scheme turns out to be twice as large as that found within the framework of a similar calculation scheme (CBS-Q//B3 for the gas phase and B3LYP/6-311+G** IEFPCM for the energy of solvation) by Khursan and Ovchinnikov [30] in their article,  [30] used different scale factors α for various acids, and, in addition, the error of the pKa calculation was decreased by correlation allowances for each kind of C-H, O-H, N-H and S-H acids. In contrast, our goal is not obtaining the accurate pKa values, but a comparison of different schemes for calculating the solvation energies of neutral and anionic forms in a uniform manner.

B2PLYP based acidities
In this series, the Gibbs free energy in the gas phase was calculated using the B2PLYP-D/6-311+G**// B3LYP/6-31+G* method. To calculate the solvation energies, in addition to the SMD, PCM-LB and PCM-ES schemes described above, we used i) the PCM-SB scheme, similar to PCM-LB, but with a reduced basis set 6-31+G*, and ii) the PCM-D scheme, also using the 6-31+G* basis set, the default in G09 value α = 1.10, and taking into account only the electrostatic component. Table 3. B2PLYP/6-311+G**//B3LYP/6-31+G* based acidities obtained using various schemes (see text for notation) compared to experimental [29]  The results presented in table 3 show that all three PCM schemes using a factor of α = 1.35 give similar acidities for the test compounds, which are in reasonable agreement with the experimental data. An unfortunate exception is 1,1-diphenylmethanimine, for which the calculated pKa value is 6 units higher than the experimental one. In the case of PCM-D, which uses a default multiplier, this error reaches 9 units, while the SMD scheme with the same cavity shape has no difficulty with this molecule. The exclusion of 1,1-diphenylmethanimine from the test set reduces the MAD value significantly, for example, for the PCM-ES the MAD decreases to 1.7 pKa units.

Discussion
The above results allow us to recommend the B2PLYP-D/6-311+G**//B3LYP/6-31+G* method in combination with the calculation of solvation energies in the framework of the continuous model (PCM)-B3LYP/6-311+G** with a scaling factor α = 1.35 and taking into account only the electrostatic component (PCM-ES) to simulate the reactions in KOH(KOBu t )/DMSO superbasic media. In conclusion, it is appropriate to make some remarks about its applicability to describe different classes of compounds exhibiting different types of acidity.