How Accurate is the SMD Model for Predicting Free Energy Barriers for Nucleophilic Substitution Reactions in Polar Protic and Dipolar Aprotic Solvents ?

In the past seven years, the SMD (Solvent Model Density) method has been widely used by computational chemists. Thus, assessment on the reliability of this model for modeling chemical process in solution is worthwhile. In this report, it was investigated six anion-molecule nucleophilic substitution reactions in methanol and dipolar aprotic solvents. Geometry optimizations have been done at SMD/X3LYP level and single point energy calculations at CCSD(T)/TZVPP + diff level. Our results have indicated that the SMD model is not adequate for dipolar aprotic solvents, with a root of mean squared (RMS) error of 5.6 kcal mol, while the Polarizable Continuum Model (PCM) method with the Pliego and Riveros atomic cavities has led to RMS error of only 3.2 kcal mol. For methanol solvent, the SMD method has a RMS error of 3.0 kcal mol. The classical protic to dipolar aprotic solvent rate acceleration effect was not predicted by the SMD model for the tested systems.


Introduction
In the middle of the twentieth century, several experimental studies have led to foundation of empirical rules for qualitative prediction of solvent effects on chemical reactions, authoritatively reviewed by Parker 1 in a classical paper.1][12][13][14][15][16][17] On the other hand, gas phase dynamic of chemical reactions can open new mechanistic possibilities. 18,19mong the different solvents used to conduct ionic chemical reactions, polar protic and dipolar aprotic solvents are paramount due to their property of solubilize ionic species. 20Methanol is a polar protic solvent and its ability to solvate ions is similar to water. 21,22In addition, it has the advantage of solubilize many organic molecules not soluble in water.Some usual dipolar aprotic solvents are dimethyl sulfoxide (DMSO), dimethylformamide (DMF) and acetonitrile.Because anions are less solvated in these solvents, anion-molecule reactions are accelerated when transferred from protic to dipolar aprotic solvents.Our ability to describe theoretically this effect is an important goal and a challenge for continuum solvation models once requires these models be able to describe specific solutesolvent interactions. 23Considering that the Born formula for solvation of a spherical ion of charge q and radius R into a solvent with dielectric constant ε is given by: (1)   it could be noted that the use of the same "molecular cavity" for an ionic solute into a protic and dipolar aprotic solvent with high dielectric constant would lead to very similar solvation free energy of ions.A possibility to overcome this problem would be to use different molecular cavity for protic and dipolar aprotic solvents.This approach has reached some successful ten years ago using the polarizable continuum model (PCM) in the study of anion-molecule reactions. 24n the present time, continuum solvation models like SMD (solvation model density) of Truhlar and co-workers 25 are widely used to study chemical reactions.In particular, this model has the advantage of being parametrized for more than 100 solvents.It is very important to know how reliable the SMD method is, using different tests to assessment the reliability of the model.7][28][29][30][31][32][33] The SMD model has presented a good performance for challenging neutral molecules in water solvent. 34Other successful application of the SMD model was the prediction of phase separation. 35Nevertheless, we should emphasize that these tests were restrict to neutral molecules.][38] However, even more important than predicting absolute solvation free energy values is to predict correct variation of the solvation free energy between transition states and reactants as well as between products and reactants.This ability would lead to reliable prediction of solvent effects on chemical reactions.
The aim of this work is to evaluate the performance of the SMD continuum model for anion-molecule reactions in protic and dipolar aprotic solvents.The SMD model has been proposed seven years ago and has become a very popular model to investigate solvent effects on organic reactions.Thus, an evaluation of its performance for ion-molecule reactions in solvents like methanol, DMSO and DMF is worthwhile.We have chosen six ion-molecule nucleophilic substitution reactions that are presented in the Table 1.

Theoretical methods
The potential energy surfaces for the studied reactions were investigated using density functional theory, including the solvent effect through a continuum solvation model.Thus, potential of mean force of the reaction in solution were explored.The search for minimum and transition state structures were done using the X3LYP functional 39,40 with the Dunning DZV basis set augmented with d polarization and sp diffuse functions on the heavy atoms as implemented in GAMESS. 41,42The continuum solvation model SMD 25 was used in the geometry optimization and in the harmonic frequency calculations.Considering that the solvent effect on geometry and harmonic frequency in polar solvents are close, methanol was used as solvent for all the calculations of geometry and frequency.
In order to improve the gas phase electronic energy, single point energy calculation was done with the functionals X3LYP and M11 43 using more extended TZVPP basis set of Ahlrichs, 44 augmented with sp diffuse functions on the O, N, S, Cl and Br atoms.For the M11 method it was used the JANS = 2 option in GAMESS to use a finer grid.Accurate wave function CCSD(T) method was also used with this TZVPP + diff basis set.In some cases involving large molecules (reactions 4 to 6), the CCSD(T)/TZVPP + diff energies were estimated using the additivity approximation through calculations at MP2 and CCSD(T) levels with the SVP + diff basis set and MP2 calculations with the TZVPP + diff basis set.
The solvation free energy of the reactants and transition states were calculated at SMD/X3LYP/DZV + P(d) + diff level for each solvent investigated (methanol, DMSO and DMF).Additional calculations were done with the PCM/X3LYP/DZV + P(d) + diff method and using the Pliego and Riveros 45 atomic cavities for DMSO solvent.These parameters were used for both DMSO and DMF solvents, once previous studies in DMF also suggest similar atomic cavities. 46he activation free energy values were calculated through the equations: (2) (3) In equation 2, the first term in the right side is equivalent to gas phase free energy contribution (ΔG ‡ mol ).However, because the geometries and frequencies were obtained in solution, it is more correct to name this term as "molecular" contribution to the free energy.The second term is the solvation free energy contribution calculated by the SMD method (ΔΔG ‡ solv ).In equation 3, the first term in the right side is the electronic energy contribution obtained from single point energy calculations (ΔE ‡ el ), and the second term is the vibrational, rotational and translational energy contributions obtained from the harmonic frequency calculations in solution (ΔG ‡ vrt ).Correction for standard state of 1 mol L -1 was done for all the structures.
The calculations using the X3LYP and M11 methods, as well as the SMD model, were done with the GAMESS program. 41,42The MP2 and CCSD(T) calculations were done with the ORCA program. 49

Results and Discussion
The geometries of the transition states determined in this work are presented in Figure 1.It can be observed the usual backside attack mechanism of S N 2 reactions, with exception of TS6, which corresponds to S N Ar reaction.
It was used three levels of theory to obtain the electronic energy contribution.Considering that density functional theory is widely used in the investigation of organic reactions, it is interesting to know how reliable the functionals are.In this work, it was tested the X3LYP and the M11 functionals.The X3LYP method is slightly more accurate than the very popular B3LYP one, whereas the M11 functional of Peverati and Truhlar 43 has been claimed to have a good performance on a wide range of properties.These electronic energies were compared with the gold standard CCSD(T) wave function approach.The results are presented in Table 2 and Figure 2. It can be observed that the M11 functional is slightly more accurate than the X3LYP one.Based on previous report, 43 an even better performance of the M11 functional was expected.However, an interesting point to notice is that the X3LYP method systematically underestimate the activation barrier for S N 2 reactions, while it overestimates in the case of S N Ar reaction.On the other hand, the M11 functional seems to be more regular.
The next aspect to be investigated is how the SMD model predicts the solvent effect in the free energy of activation in both protic and dipolar aprotic solvents.) and the squares the theoretical solution phase barriers (ΔG sol ‡ ), both compared to experimental solution phase barriers.This comparison is useful to point out how the theoretical SMD model correct the "gas phase" barriers.It can be noted that the ΔG mol † values are highly dispersed.For example, there is a very high deviation of TS6 in methanol, and a small deviation of TS5 in DMSO.This behavior is somewhat expected.Previous study has already indicated that anionmolecule S N Ar reactions have higher solvent effect than S N 2 reactions. 50On the other hand, the PhS -nucleophile has a considerable charge dispersion, resulting in small solvent effect in TS4 and TS5.
An inclusion of the solvent effect (ΔG sol ‡ ), even by a simple continuum model, leads to considerable improvement in the theoretical activation barriers (Figure 3).In fact, the barriers become less dispersed and close to the y = x curve.The overall root of mean squared (RMS) error is 4.3 kcal mol -1 .When analyzing each solvent separately, the RMS error becomes 3.0 kcal mol -1 for methanol and 5.6 kcal mol -1 for dipolar aprotic solvents.These results are unexpected considering that it is well documented that continuum models work better for anions in aprotic solvents than in protic solvents. 51,52In fact, several theoretical studies of anion-molecule S N 2 reactions in DMSO using the PCM method and the Pliego and Riveros 45 atomic cavities have led to reliable activation barriers. 24,53,54In part, the reason for this high deviation can be attributed to reaction of the PhS -nucleophile in DMF, with a deviation of 8.3 kcal mol -1 for TS4.These results suggest that the radius of sulfur atom of the SMD model in aprotic solvents needs to be improved.
Analysis of the classical rate acceleration effect on going from protic to dipolar aprotic solvent indicates a disappointing result.Observing Table 2, we can notice that the barriers in methanol and DMF (or DMSO) are very close.Similarly, Figure 3 points out the squares are displaced to the right in relation to the triangles.To reproduce the rate acceleration effect, the displacement should be parallel to the y = x curve.These observations point out the SMD model is not able to predict this effect.The reason for this fail can be attributed to the definition of the cavities in the model.Indeed, continuum solvation models have two contributions to the solvation free energy (ΔG * solv ).The first term is the electrostatic contribution (ΔG el ) and is related to the dielectric constant of the solvent while the noneletrostatic term (ΔG nel ) is related to cavity formation and dispersion contribution as indicated below: (4) Considering that many protic and dipolar aprotic solvents has high dielectric constant, the resulting continuum electrostatic contribution to the solvation free energy are very similar, as observed in equation 1.Thus, one possibility is introducing a variable solute cavity for each solvent.Such approach has been used partially in the SMD model, because solely the oxygen atom has a variable radius.Consequently, the barriers calculated in this article are not sensible to the solvent.
In previous test of the performance of the SMD model for pKa prediction in methanol, it was used a proton exchange reaction. 55That study has indicated that even an isodesmic reaction leads to reasonable deviation of the experimental data, suggesting a fail of the SMD model for ions in methanol.The modest error observed in this report can be due the fact the studied ions have high charge dispersion, with exception of the methoxide ion in TS6.This barrier has presented the highest deviation in methanol (4.6 kcal mol -1 ).
Considering that previous studies on anion-molecule S N 2 reactions in DMSO using the PCM method with the Pliego and Riveros 45 cavities has worked fine, this approach was also tested.The results are presented in Table 2, with the values in parenthesis.There is a considerable improvement in the results and the RMS error becomes 3.2 kcal mol -1 .The highest deviations occur for TS4 and TS5, which involves sulfur nucleophiles.These findings reinforce the idea that the SMD model needs to be improved for dipolar aprotic solvents.
It would be worthwhile to know the performance of diverse continuum solvation models for studying ionic reactions in solution.For example, the COSMO-RS model 56,57 has been tested for proton exchange and recombination reactions with reasonable accuracy. 58Other promising approaches are the SMVLE 59 and CMIRS [60][61][62] methods.Essentially, these methods include a term to describe the effect of extremum electric field generated by charged solutes.This high electric field is not properly accounted for the dielectric continuum solvation.In fact, the reported studies have indicated this field-extremum correction term is able to make a substantial improvement on the ionic solvation thermodynamics.][68] The use of full explicit solvent molecules coupled with free energy perturbation was pioneered by Jorgensen, 69 and the rate acceleration effect was predicted for a classical S N 2 reaction.Acevedo and Jorgensen 70 have also investigated an anion-molecule S N Ar reaction in different solvents using explicit solvent (quantum mechanics/molecular mechanics methods, QM/MM) and the PCM method.While explicit solvent has reproduced the methanol to DMSO rate acceleration effect, the PCM method has failed; recent reviews on applications of QM/MM methods can be found. 71,72he issue of determining the best solvent for a specific chemical reaction has been recently addressed by combining quantum chemistry calculations with computeraided molecular design. 73Thus, accurate prediction on solvent effects remain an important goal and can be very useful in the future studies of the design of best solvent for diverse chemical reactions of practical interest. 74

Conclusions
Analysis of the performance of the SMD model to predict free energy of activation of S N 2 and S N Ar reactions in methanol, DMF and DMSO solvents has shown that this model is able to capture the most important solvent effect.Thus, the theoretical barriers are in reasonable agreement with the experimental data.When finer details such as the protic-dipolar aprotic solvents rate acceleration effect was analyzed, it becomes evident the SMD method fails to predict this effect in the investigated reactions.In addition, the performance of the model is worse for aprotic solvents, while a simpler PCM model with the Pliego and Riveros 45 atomic cavities has a better performance.Although part of the error in the predicted barriers could be attributed to the gas phase contribution (error in the CCSD(T)/TZVPP + diff method), our results suggest the SMD model needs to be improved to describe ionic reactions in different solvents.

Figure 2 .
Figure 2. Activation electronic energy barriers.Performance of the X3LYP and M11 functionals in relation to the CCSD(T) method.The RMS error is 3.04 and 2.55 kcal mol -1 for X3LYP and M11, respectively.All the calculations with the TZVPP + diff basis set and SMD/X3LYP/ DZV + P(d) + diff geometries.

Figure 3
Figure3presents the results.The circles correspond to the theoretical "gas phase" barriers (ΔG mol ‡

Table 2 .
45eoretical and experimental activation properties aStandard state of 1 mol L -1 ; b electronic energy contribution.All calculations with the TZVPP + diff basis set and geometries obtained at SMD/X3LYP/ DZV + P(d) + diff; c molecular (gas phase) contribution to the free energy of activation, using CCSD(T) electronic energies; d solution phase free energy of activation using the SMD model; e reactions 1 to 4 in dimethylformamide (DMF) and 5 in dimethyl sulfoxide (DMSO); f values in parenthesis correspond to the PCM/X3LYP/DZV + P(d) + diff method using the Pliego and Riveros45atomic cavities.Theoretical free energy of activation versus experimental value.Circles correspond to theoretical "gas phase" barrier versus experimental solution phase values.Squares and triangles correspond to theoretical barriers versus experimental solution phase values.Units in kcal mol -1 .Vol. 27, No. 11, 2016 a