Benchmarking Continuum Solvent Models for KetoâEnol

Article pubs.acs.org/JPCA

Benchmarking Continuum Solvent Models for Keto−Enol Tautomerizations Billy W. McCann, Stuart McFarland, and Orlando Acevedo* Department of Chemistry and Biochemistry, Auburn University, Auburn, Alabama 36849, United States S Supporting Information *

ABSTRACT: Experimental free energies of tautomerization, ΔGT, were used to benchmark the gas-phase predictions of 17 different quantum mechanical methods and eight basis sets for seven keto−enol tautomer pairs dominated by their enolic form. The G4 method and M06/6-31+G(d,p) yielded the most accurate results, with mean absolute errors (MAE’s) of 0.95 and 0.71 kcal/mol, respectively. Using these two theory levels, the solution-phase ΔGT values for 23 unique tautomer pairs composed of aliphatic ketones, β-dicarbonyls, and heterocycles were computed in multiple protic and aprotic solvents. The continuum solvation models, namely, polarizable continuum model (PCM), polarizable conductor calculation model (CPCM), and universal solvation model (SMD), gave relatively similar MAE’s of ∼1.6−1.7 kcal/mol for G4 and ∼1.9−2.0 kcal/mol with M06/6-31+G(d,p). Partitioning the tautomer pairs into their respective molecular types, that is, aliphatic ketones, β-dicarbonyls, and heterocycles, and separating out the aqueous versus nonaqueous results finds G4/PCM utilizing the UA0 cavity to be the overall most accurate combination. Free energies of activation, ΔG‡, for the base-catalyzed keto−enol interconversion of 2-nitrocyclohexanone were also computed using six bases and five solvents. The M06/6-31+G(d,p) reproduced the ΔG‡ with MAE’s of 1.5 and 1.8 kcal/mol using CPCM and SMD, respectively, for all combinations of base and solvent. That specific enolization was previously proposed to proceed via a concerted mechanism in less polar solvents but shift to a stepwise mechanism in more polar solvents. However, the current calculations suggest that the stepwise mechanism operates in all solvents.

■

INTRODUCTION Keto−enol tautomerism entails both a molecular interconversion of a proton and a shifting of π electrons that results in the equilibration of two unique isomers, that is, the keto and enol tautomers (Figure 1).1,2 Thermodynamically driven, the

Table 1. Average Energies (kcal/mol) for Bonds Distinguishing Keto and Enol Tautomersa keto form C−H C−C CO sum a

equilibrium preference for either tautomer may be rationalized energetically. A simple distinction can be made between the compounds via three unique bonds; that is, the keto tautomer possesses C−H, C−C, and CO bonds, whereas the enol tautomer contains O−H, CC, and C−O bonds. A straightforward analysis of their average bond energies (Table 1) finds the keto form to be more stable by ca. 12 kcal/mol.1 This often renders the enolic form undetectable, as the tautomers establish an equilibrium constant KT that is proportional to their free energy difference eq 1.3 [enol] = exp( −(Genol − G keto)/kBT ) [keto] © 2015 American Chemical Society

O−H CC C−O

110−111 146−151 85−91 341−353

Reference 4.

There are cases where the enol may be present to a significant degree.5 A classic example is phenol, which is completely dominated by the enolic form, KT = 4 × 1013,6 as tautomerization to the keto would endure an unfavorable loss of aromaticity (Figure 2).7 Bulky substituents at the α-carbon position can also shift the equilibrium toward the enolic form.8

Figure 1. Keto−enol tautomerism.

KT =

enol form 96−99 83−85 173−181 352−365

Figure 2. Phenol keto−enol tautomerization. Received: April 29, 2015 Revised: July 21, 2015 Published: July 23, 2015

(1) 8724

DOI: 10.1021/acs.jpca.5b04116 J. Phys. Chem. A 2015, 119, 8724−8733

Article

The Journal of Physical Chemistry A

form.28 Solvation of acetylacetone by apolar aprotic solvents finds KT ratios to be similar to gas-phase values, but increasing the solvent polarity by using, for example, N,N-dimethylformamide considerably enhances the formation of the keto tautomer.29,30 The solvent effect can be rationalized by an increased competition for hydrogen bonding to the carbonyl oxygen atoms between the intramolecular proton of the enolic form and the polar solvent molecules. Calculating solvation energy using continuum models has a long history of success31,32 but can be challenging for certain classes of molecules. For example, predicting tautomer ratios for cyclic alkanes using continuum solvent models has been shown to be problematic.33 In this study, the ability of quantum mechanical (QM) methods to reproduce solvent effects on keto−enol ratios was investigated using the integral equation formalism polarizable continuum model (IEFPCM is synonymous with PCM)34 and the polarizable conductor calculation model (CPCM)35 with six different cavities each, namely, BONDI, PAULING, UA0, UAHF, UAKS, and UFF.36 Truhlar’s universal solvent model (SMD) was also examined.37 Twenty-three pairs of compounds were organized into three groups containing similar structures or functional groups that included aliphatic ketones (Figure 5), β-dicarbonyl molecules

For example, steric hindrance derived from bulky aryl groups is better tolerated in the enol tautomer through a 120° separation of the substituents compared to the more crowded 109.5° coordination present in the keto form.9,10 A final example of enol preference may be found in β-dicarbonyl compounds, such as acetylacetone (Figure 3), where the enol tautomer benefits from the simultaneous formation of a conjugated π-system and an intramolecular hydrogen bond.11−13

Figure 3. Keto−enol tautomerization of acetylacetone.

In the current computational study, the first objective was to reproduce experimentally determined free energies of tautomerization, ΔGT, for compounds favoring the enol form in the gas phase. Seven compounds were selected based on the availability of experimental data: 2-nitrocyclohexanone, acetylacetone, ethyl acetoacetate, methyl acetoacetate, acetylacetamide, 2pyridone, and 6-chloropyridone (Figure 4). In the case of the

Figure 5. Tautomerizations for aliphatic ketones.

(Figure 6), and heterocycles (Figure 7). The enolization ratios were computed in multiple solvents including water, chloroform, dichloromethane, carbon tetrachloride, acetonitrile, and cyclohexane. The heterocycles, in particular, provided tautomerization examples beyond the traditional keto/enol, for example, amino/imino tautomerism in cytosine and 3methylcytosine, and lactam/lactim tautomerism in 2-pyridone and 6-chloropyridone. Finally, the free energies of activation, ΔG‡, for the basecatalyzed enolization of 2-nitrocyclohexanone in cyclohexane, carbon tetrachloride, chloroform, dichloromethane, and acetonitrile were computed and compared to experimental values. Six bases were tested in the catalysis of 2NCH that included: pyridine (Pyr), 3-methylypyridine (3MePyr), 3methyloxypyridine (3MeOPyr), 4-methylpyridine (4MePyr), 4-methoxypyridine (4MeOPyr), and triethylamine (Et3N). A detailed investigation using M06/6-31+G(d,p) and the CPCM and SMD solvent models provided ΔG‡ values for all combinations of bases and solvents. All benchmarking results are discussed in detail herein, and recommendations are given based upon the findings.

Figure 4. Tautomerizations computed in the gas-phase. In all cases the enol form is preferred.

pyridones, a lactam (keto) and lactim (enol) tautomerization is occurring. A wide variety of theoretical methods were employed to benchmark ΔGT values between the tautomers that included: Hartree−Fock (HF); eight density functional theory (DFT) methods: B3LYP,14 X3LYP,15 M05, M05-2X, M06, M06-2X,16,17 ωB97X-D,18 and PBEh1PBE;19 five composite thermochemical models: CBS-Q,20 CBS-4M,21 CBS-QB3,21 CBS-APNO,22 and G4;23 and finally, three semiempirical quantum mechanical (SQM) models: AM1,24 PDDG/PM3,25 and PM6.26 Nine different basis sets were tested for each of the HF and DFT methods. Temperature and concentration considerably affect enolization. eq 1 underlines the correlation between temperature and tautomer equilibrium. For example, whereas gas-phase acetylacetone (Figure 3) has a 95% enol abundance at 50 °C,27 an increase of temperature to 155 °C reduces the amount to 78%.12 Solvent can also dictate the stability of the enolic 8725


Article


implemented in Gaussian 09 using the default settings. The SMD model included solute−solvent repulsion interaction, solute cavitation, and solute−solvent dispersion interaction energies. Molecular energies and other properties are dependent upon the choice of atomic radii used to build the cavity into which the solute is placed. Therefore, all cavity types included in Gaussian 09, namely, BONDI, PAULING, UA0, UAHF, UAKS, and UFF, were explored.36 The COULOMB atomic radii was optimized for the SMD model.37 Radii values for atoms in this study are given in Table 2. Table 2. Radii (Å) of Atomic Spheres for Each Cavity Type BONDI PAULING UA0 UAHF UAKS UFF COULOMB

Figure 6. Tautomerizations for β-dicarbonyl molecules.

C

O

N

Cl

S

H

1.70 1.50 2.00 1.80 1.80 1.93 1.85

1.52 1.40 1.80 1.50 1.50 1.75 2.02

1.55 1.50 1.90 1.60 1.60 1.83 1.89

1.75 1.80 1.97 1.98 1.98 1.97 2.38

1.80 2.26 2.22 2.11 2.11 2.02 2.49

1.20 1.20

1.44 1.20

■

RESULTS AND DISCUSSION Gas Phase. The gas-phase free energies of tautomerization, ΔGT, for 2-nitrocyclohexanone, acetylacetone, ethyl acetoacetate, methyl acetoacetate, acetylacetamide, 2-pyridone, and 6chloropyridone (Figure 4) were computed using 17 different QM methods and eight basis sets. Table 3 provides the mean absolute error (MAE) in the ΔGT values for the most accurate QM/basis set combination per molecule set. Figure 8 provides a graphical representation of Table 3, including the maximum and minimum deviation from experiment for each method. Unsigned errors in the calculated gas-phase ΔGT for every QM method and basis set combination are exhaustively provided in the Supporting Information. Of the composite methods tested, CBS-APNO and G4 gave the lowest MAE values in ΔGT with 0.82 and 0.95 kcal/mol, respectively, for the gas-phase tautomerizations. Intriguingly, all composite methods gave a large error for the predicted ΔGT of 2-nitrocyclohexanone (2NCH). Both the axial and equatorial conformations of the 2NCH keto tautomer failed to achieve the kilocalorie accuracy normally associated with these thermochemical QM models. A detailed exploration of the enolization mechanism and solvent effects for 2NCH is discussed later. Visual inspection of Figure 8 will readily reveal the M06/631+G(d,p) theory level to be energetically the most accurate method, more so than the widely used B3LYP and the computationally demanding G4 and CBS-APNO methods. M06 yielded an MAE of 0.71 kcal/mol and had the smallest maximum deviation from experiment for the gas-phase molecule set. This is in contrast to an earlier study by Piacenza and Grimme that reported poor gas-phase energy predictions for six DNA base pairs and related molecules when employing DFT methods, such as HCTH407, PBE, BP86, BLYP, B3LYP, and BHLYP, relative to QCISD(T)-derived reference energies; high-level ab initio methods provided better comparisons, for example, SCS-MP2 gave a root-mean-square deviation (rmsd) of 0.7 kcal/mol.44 Whereas it was hypothesized that the problem may lie in DFT’s inability to accurately model both aromatic and nonaromatic tautomeric forms simultaneously using the same exchange terms,44 a more recent study by Galstyan and Knapp found that another one of Truhlar’s DFT

Figure 7. Tautomerizations for heterocycle molecules.

■

COMPUTATIONAL METHODS Eight basis sets were tested for the HF and DFT methods that included 6-31+G(d,p), 6-311++G(2d,p), 6-311++G(3df,3pd), cc-pVDZ, cc-pVTZ, cc-pVQZ, aug-cc-pVDZ, and aug-cc-pVTZ. All molecules were fully optimized. Frequency calculations were performed to verify that all resultant geometries were minima on the potential energy surface and to provide thermodynamic corrections to the electronic energy in both vacuum and solution. All calculations were performed using the Gaussian 09 program38 on computers located at the Alabama Supercomputer Center. Solution-phase calculations employed the continuum solvation models PCM,34 CPCM,35 and SMD37 as 8726


Article


Table 3. Individual Unsigned Errors and the Mean Absolute Error (MAE) in the Calculated Gas-Phase Free Energies of Tautomerization, ΔGT (kcal/mol), at 25 °C for Seven Molecules moleculea

2NCH

AcAc

EAA

MAA

2Pyr

6ClPyr

AAM

exptl ΔGT method CBS-QB3 CBS-Q CBS-APNO CBS-4M G4 HF/aug-cc-pVTZ B3LYP/6-311++G(2d,p) X3LYP/6-31+G(d,p) M05/cc-pVTZ M05−2x/6-31+G(d,p) M06/6-31+G(d,p) M06-2X/6-311++G(2d,p) wB97XD/6-31+G(d,p) PBEh1PBE/6-311++G(2d,p) AM1 PDDG/PM3 PM6

−3.53b

−2.2c

−0.08c

−0.85d unsigned errors 0.32 1.74 0.04 1.65 0.26 2.98 1.47 1.74 0.39 3.88 0.97 1.78 0.00 2.43 2.46 5.94 1.13

−0.69e

−2.1f

−0.37g

2.83 4.09 2.52 5.73 3.72 6.28 0.68 1.03 2.48 0.52 0.07 0.51 1.62 1.76 2.15 2.13 3.73

0.68 0.94 0.15 2.58 0.07 3.42 2.23 1.69 0.54 2.83 1.29 0.69 0.80 3.47 3.23 8.50 0.70

0.92 3.23 0.77 0.93 0.99 2.29 2.19 2.43 0.80 3.24 0.84 1.78 2.15 3.15 1.48 5.16 2.71

0.59 0.69 0.47 1.60 0.25 0.23 1.16 1.01 2.15 1.24 0.82 1.17 0.97 0.76 0.37 2.71 6.42

1.55 0.25 h 2.38 1.21 1.53 0.27 0.50 0.63 2.65 0.77 2.42 0.63 0.66 1.66 1.11 9.00

0.79 0.77 0.94 2.49 0.16 2.66 0.75 0.86 1.86 1.88 0.22 0.75 0.80 1.71 0.73 3.19 6.10

MAE 1.10 1.67 0.82 2.44 0.95 2.77 1.25 1.32 1.27 2.32 0.71 1.30 1.00 1.99 1.73 4.11 4.26

Molecule abbreviations are defined in Figure 4. Negative experimental ΔGT values indicate a preference for the enol form. b25 °C, ref 39. c25 °C, ref 27. d36 °C, ref 40. e67 °C, ref 41. f120 °C, ref 42. g74 °C, ref 43. hChlorine is not implemented in the CBS-APNO method. a

Figure 8. Mean absolute error (●) and the maximum (□) and minimum (−) unsigned error in the calculated gas-phase free energies of tautomerization, ΔGT (kcal/mol), at 25 °C for each QM combination given in Table 3.

functionals PW6B9545 did yield accurate gas-phase tautomer energy differences for a set of 36 tautomer/isomer pairs of lactams.46 Error ranges were similar to the results reported here for M06. For example, PW6B95 gave an rmsd of 0.3 kcal/mol compared to five experimental values and an rmsd of 0.67 kcal/ mol relative to energies derived using the QCISD(T)(q-ζ) method on 15 compound pairs.46 A detailed comparison of M06 with multiple basis sets is given in Table 4. While the modestly sized 6-31+G(d,p)

produced the lowest MAE, alternative basis sets gave more accurate unsigned errors in ΔGT for some of the molecules on an individual case-by-case basis. For example, 6ClPyr gave an unsigned error of 0.02 kcal/mol when using aug-cc-pVTZ compared to 0.77 kcal/mol for 6-31+G(d,p) (Table 4). While the current gas-phase test set is far from exhaustive with only seven tautomer pairs (Figure 4), it does provide a good starting point for selecting QM methods for a more expansive test in solution. 8727


Article


Table 4. Individual Unsigned Errors and the Mean Absolute Error (MAE) in the Calculated Gas-Phase Free Energies of Tautomerization, ΔGT (kcal/mol), at 25 °C for Seven Molecules Using the M06 Method with Different Basis Setsa 2NCH basis set 6-31+G(d,p) 6-311++G(2d,p) 6-311++G(3df,3pd) cc-pVDZ cc-pVTZ cc-pVQZ aug-cc-pVDZ aug-cc-pVTZ a

AcAc

0.07 0.12 0.35 0.56 0.27 0.66 0.61 0.05

1.29 1.53 1.41 2.14 1.25 0.83 1.89 1.52

EAA

MAA

2Pyr

6ClPyr

AAM

0.84 1.35 1.69 1.47 1.49 1.57 1.76 1.96

unsigned errors 0.97 0.91 0.60 0.90 1.10 0.73 1.63 0.92

0.82 0.85 1.12 0.56 1.16 1.22 0.37 1.17

0.77 0.48 0.06 1.06 0.24 0.18 1.24 0.02

0.22 0.19 0.04 0.48 0.02 0.47 0.61 0.43

MAE 0.71 0.78 0.75 1.02 0.79 0.81 1.16 0.87

Molecule abbreviations are defined in Figure 4.

Table 5. Continuum Solvent Model and Cavity Combination for Predicting Solution-Phase Free Energies of Tautomerization, ΔGT (kcal/mol), at 25 °C for 23 Compoundsa M06/6-31+G(d,p) MAE

G4

MAX

STDDEV

MAE

PCM BONDI PAULING UA0 UAHF UAKS UFF

1.90 2.08 1.98 2.38 3.17 1.99

BONDI PAULING UA0 UAHF UAKS UFF

1.87 2.09 1.99 2.40 3.42 1.99

COULOMB

2.05

5.47 7.04 6.67 6.60 9.17 5.70

5.18 5.34 5.02 3.96 4.15 5.32


1.75 1.96 1.56 2.62 2.86 1.71

5.17 6.15 6.70 6.69 9.94 5.73

4.74 5.21 4.99 3.99 4.24 5.30


1.80 2.05 1.60 2.57 2.96 1.75

5.11

5.15

COULOMB

1.91

CPCM

STDDEV

5.80 5.92 5.86 8.31 8.31 5.10

4.33 4.40 4.24 3.61 3.63 4.61

5.85 5.98 5.89 7.37 8.22 5.13

4.32 4.35 4.16 3.53 3.61 4.57

7.98

4.13

CPCM

SMD a

MAX PCM

SMD

Mean absolute error (MAE), maximum unsigned error (MAX), and standard deviation (STDDEV).

Table 6. Mean Absolute Error (MAE) in Calculated Free Energies of Tautomerization, ΔGT (kcal/mol), at 25 °C Separated by Water and Nonaqueous Solvents M06/6-31+G(d,p) MAE BONDI PAULING UA0 UAHF UAKS UFF BONDI PAULING UA0 UAHF UAKS UFF COULOMB

water PCM 2.18 2.34 2.30 1.78 2.33 2.67 CPCM 2.05 2.28 2.28 1.81 2.46 2.67 SMD 2.42

G4

nonaqueous

MAE

1.69 1.87 1.73 2.81 3.76 1.45


1.73 1.94 1.77 2.86 4.16 1.46


1.78

COULOMB

Condensed Phase. Solvent is known to have a large effect on tautomer preference.47−49 Accurately predicting the energetic difference between the keto−enol tautomers in

water PCM 2.02 2.21 1.94 2.40 2.19 2.47 CPCM 2.07 2.23 1.92 2.42 2.24 2.50 SMD 2.34

nonaqueous 1.53 1.77 1.28 2.77 3.34 1.11 1.59 1.91 1.35 2.68 3.52 1.17 1.58

solution is not trivial, as many factors contribute including the potential of multiple rotameric forms and the specific interaction of solvent molecules with hydrophilic functional 8728


Article


Table 7. Mean Absolute Error (MAE) in Calculated Free Energies of Tautomerization, ΔGT (kcal/mol), at 25 °C Separated by Aliphatic Ketones, β-Dicarbonyl Molecules, and Heterocycles in All Solventsa M06/6-31+G(d,p) MAE

β-carb

ketone

G4 heterocyc

MAE

ketone

PCM BONDI PAULING UA0 UAHF UAKS UFF

4.51 5.63 3.24 1.01 0.83 4.68


4.30 5.39 3.21 1.06 0.85 4.68

COULOMB

4.18

heterocyc

1.55 1.88 1.57 3.36 3.28 1.39

1.70 1.66 1.58 2.22 2.78 1.67

1.63 2.02 1.57 2.99 2.94 1.44

1.73 1.74 1.66 2.38 3.16 1.71

1.38

2.18

PCM 1.08 1.21 1.68 2.82 2.92 1.43

1.86 1.92 1.82 2.43 3.92 1.66


2.55 3.36 1.50 2.39 2.18 2.87

1.07 1.24 1.60 2.51 2.68 1.34

1.95 1.99 1.94 2.59 4.36 1.74


2.56 3.35 1.48 2.46 2.25 2.90

0.96

2.34

COULOMB

1.80

CPCM

CPCM

SMD a

β-carb

SMD

See Figures 5−7.

homogeneous, dielectric continuum ignores explicit hydrogen bonding. Discrete-continuum models that include one or more explicit water molecules hydrogen bonded to the solute to represent first solvation shell interactions along with an implicit solvent model to simulate the bulk phase have been shown to improve energetic agreement with experiment.57,58 Alternatively, calculating the electrostatic energy of solvation and gasphase free energies independent of one another has been shown to yield excellent results in both predicting pKa values in water59 and computing one-electron reduction potentials in protic and aprotic solvents for a set of diverse organic molecules.60 Further partitioning the 23 compounds into three groups containing similar structures or functional groups, that is, aliphatic ketones (Figure 5), β-dicarbonyl molecules (Figure 6), and heterocycles (Figure 7), revealed that MAE predictions vary widely depending on the solvent model and cavity chosen (Table 7). The aliphatic ketones group faired the poorest, particularly when using M06/6-31+G(d,p) with the Gaussian 09 defaulted UFF cavity. However, that set of compounds only had experimental values reported in aqueous solution for MAE comparison, which was shown in Table 6 to give inferior predictions when compared to nonaqueous solvents. Surprisingly, switching to the UAHF and UAKS cavities dramatically improved the MAEs; for example, M06/PCM gave values of 1.01 and 0.83 kcal/mol, respectively (Table 7). Nevertheless, those cavities were generally unreliable as many of the βdicarbonyl compounds spontaneously deprotonated in solution when employing UAKS and UAHF. These cavities give the βdicarbonyl molecules the smallest solute volumes within the united atom class (Table 2). As a test an extra sphere was added to the hydrogen on the α-carbon using the addsphere keyword in Gaussian, but deprotonation still occurred. This effect is not likely related to the pKa value of the solute, as deprotonation did not ensue for trifluoroacetic acid using either the UAKS or UAHF cavities. Furthermore, deprotonation was not observed for any of the aliphatic ketone or heterocycle compounds. Overall, G4 in combination with PCM/UA0 generally gave the lowest solution-phase MAEs in ΔGT across all solvents and compound classes. The method also performed well for even

groups. Given the significant implications of tautomerism in areas such as drug discovery,50 a standard operating methodology for the solution-phase calculation of ΔGT would be beneficial. Implicit solvent models have been shown to be promising in a blind prediction challenge of tautomeric equilibria for 23 compounds.51−55 In that previous challenge, researchers championed their continuum solvent method of choice, for example, IEF-MST52 and SM8,53 with varying success in the predicted accuracy. In the current work, 23 alternate pairs of compounds were tested and organized into three groups containing similar structures or functional groups that included aliphatic ketones (Figure 5), β-dicarbonyl molecules (Figure 6), and heterocycles (Figure 7). Two of the more accurate methods from the gas-phase work, G4 and M06/6-31+G(d,p), were used to optimize the geometries and compute the ΔGT values in solution using the PCM, CPCM, and SMD solvation models. The free energies of tautomerization were computed in multiple solvents including water, chloroform, dichloromethane, carbon tetrachloride, acetonitrile, and cyclohexane (see the Supporting Information for a detailed description of every molecule/ solvent combination). PCM, CPCM, and SMD all gave relatively similar MAEs of ∼1.6−1.7 kcal/mol for G4 and ∼1.9−2.0 kcal/mol with M06/6-31+G(d,p) (Table 5). Changing cavities between UFF and UA0 (the Gaussian 09 and 03 default cavities,38,56 respectively) did not yield significant differences in the MAEs. However, substituting for either the UAHF or UAKS cavity did appreciably increase the MAEs, for example, G4/PCM/UAKS error of 2.86 kcal/mol, which is not surprising as their radii were optimized for the HF/6-31G(d) and PBE1PBE/6-31G(d) levels of theory, respectively. Separating the ΔGT predictions exclusively by solvent type found that continuum models gave improved MAEs in nonaqueous solvents compared to water (Table 6). For example, M06/6-31+G(d,p) with PCM/UFF gave MAEs in water and nonaqueous solvents of 2.67 and 1.45 kcal/mol, respectively; SMD gave reasonably similar values of 2.42 and 1.78 kcal/mol. The poor results in water are not unexpected when one considers that representing the solvent as a 8729


Article


methyloxypyridine (3MeOPyr), 4-methylpyridine (4MePyr), 4methoxypyridine (4MeOPyr), and triethylamine (Et3N). The CPCM/Pauling combination was selected as it yielded the most accurate MAE in ΔGT for 2NCH with unsigned errors as low as 0.14 kcal/mol (see Supporting Information). As the Pauling cavity is an unusual radii choice (despite its pronounced accuracy in this specific case), the systems were also explored using the more standard SMD solvent model as a point of comparison. Note also that calculation convergence proved difficult for some base/solvent combinations and are listed as not determined (nd) in Tables 8 and 9. The solvation model and cavity found to best reproduce the ΔGT for 2NCH, that is, CPCM/PAULING, unsurprisingly yielded the most accurate ΔG‡ results with a MAE and MAX of 1.5 and 3.3 kcal/mol, respectively (Tables 8 and 10). The SMD model gave a comparable MAE and MAX of 1.8 and 5.2 kcal/ mol, respectively (Tables 9 and 10). The Et3N-catalyzed activation barriers were substantially overestimated when using SMD, for example, a 5.2 kcal/mol deviation from experiment in dichloromethane compared to an error of 1.6 kcal/mol when using the CPCM/PAULING. The discrepancy may arise from steric hindrance due to a less favorable geometry of the methyl groups within the transition state for the SMD versus CPCM (Figure 9). Interestingly, the transition-state geometry computed for Et3N using SMD was similar to that of the optimized ground-state geometry of the isolated amine. Removal of the Et3N data lowers the SMD-computed MAE and MAX to 1.2 and 2.1 kcal/mol, respectively, whereas removal of triethylamine results for CPCM/PAULING lowers the MAE by a negligible 0.01 kcal/mol (Table 10). Regarding the mechanism of enolization, the stepwise transition structure was found with nearly identical geometries for all combinations of bases and in polar and apolar solvents. The base-catalyzed proton abstraction was computed to have an angle of ∼100 degrees for the atoms: carbonyl carbon, αcarbon, and proton (OC−Cα−H). This bond angle value is more consistent with a stepwise deprotonation, whereas a more acute angle would be expected for a concerted mechanism. All attempts to locate a concerted transition structure similar to Scheme 1 failed; the methodology was expanded from M06 to the CBS-QB3, MP2, B3LYP, and M06-2X methods with larger basis sets, for example, aug-cc-pVTZ, but every theory level gave a stepwise transition state. A separate mixed quantum and molecular mechanical (QM/MM) study utilizing Monte Carlo sampling and free energy perturbation theory (MC/FEP)63−65 is currently underway to further examine the mechanism in explicit solvent, including unique condensed-phase environments such as room-temperature ionic liquids.

the difficult aliphatic ketones group; for example, G4/PCM/ UA0 gave an MAE of 1.50 kcal/mol compared to 3.24 kcal/mol for M06/6-31+G(d,p)/PCM/UA0 (Table 7). In light of the significant computational resources needed for the G4 method, the M06/6-31+G(d,p)/PCM with either the UA0 or UFF radii still proved to be reasonably accurate for predicting ΔGT in solution. Enolization of 2-Nitrocyclohexanone. Accurate gasphase reproduction of ΔGT for 2NCH proved difficult for the routinely precise composite methods (Table 3), and calculation of the solvation energy for 2NCH placed the compound functionally somewhere in between the aliphatic ketones and the β-dicarbonyl molecules. Highlighting the importance of properly modeling solvent effects, a change in mechanism has been proposed for the base-catalyzed enolization of 2-nitrocyclohexanone when performed in polar versus apolar solvents (Scheme 1).61 In addition, solventScheme 1. Proposed Solvent Dependent Mechanism for the Enolization of 2-Nitrocyclohexanone is (A) Concerted in Less-Polar Solvents and (B) Stepwise in More-Polar Solvents

dependent rate enhancements for the reaction have been reported by Angelini et al. and are postulated to originate from a stabilization of the enolate-like transition state from their study featuring multiple solvents.62 Rate-constant measurements have found a dependency on both the base and solvent, as ΔG‡ for ketonization and enolization significantly decrease on passing from less polar solvents, for example, cyclohexane, to more polar solvents, such as acetonitrile. To assess solvent effects upon the base-catalyzed keto−enol interconversion of 2NCH, the M06/6-31+G(d,p) in combination with CPCM and the Pauling cavity was employed to calculate the free energies of activation, ΔG‡, in cyclohexane, CCl4, CHCl3, CH2Cl2, and CH3CN. The bases utilized in the reaction were pyridine (Pyr), 3-methylypyridine (3MePyr), 3-

Table 8. CPCM/PAULING Calculated Free Energies of Activation, ΔG‡ (kcal/mol), at 25 °C for the Enolization of 2Nitrocyclohexanonea,b c-C6H12

a

CCl4

CHCl3

CH2Cl2

CH3CN

base

calc

expt

calc

expt

calc

expt

calc

expt

calc

expt

Et3N Pyr 4MeOPyr 4MePyr 3MeOPyr 3MePyr

16.4 22.1

18.5 22.4

16.2 21.6

17.0 21.9

13.5 19.1 18.0 19.3 19.5 nd

15.6 20.3 18.5 19.1 20.0 19.5

13.7 18.2 17.1 18.5 18.4 17.3

15.3 20.3 18.5 19.1 20.3 19.5

13.1 16.9 16.5 nd 17.7

15.1 20.2 18.6 19.2 20.4

M06/6-31+G(d,p). nd = not determined. bExperimental values from refs 61 and 62. 8730


Article


Table 9. SMD Calculated Free Energies of Activation, ΔG‡ (kcal/mol), at 25 °C for the Enolization of 2-Nitrocyclohexanonea,b c-C6H12

a

CCl4

CHCl3

CH2Cl2

CH3CN

base

calc

expt

calc

expt

calc

expt

calc

expt

calc

expt

Et3N Pyr 4MeOPyr 4MePyr 3MeOPyr 3MePyr

21.0 23.2

18.5 22.4

21.2 22.9

17.0 21.9

18.9 22.0 20.4 nd nd 21.4

15.6 20.3 18.5 19.1 20.0 19.5

20.5 21.6 20.6 19.5 22.2 21.6

15.3 20.3 18.5 19.1 20.3 19.5

17.2 21.0 19.9 19.0 20.3

15.1 20.2 18.6 19.2 20.4

M06/6-31+G(d,p). nd = not determined. bExperimental values from refs 61 and 62.

whereas the SMD model yielded an MAE of 1.8 kcal/mol. However, SMD gave particularly poor predictions of the activation barriers in solution when catalyzed by triethylamine (Et3N), and removal of the data specific to this base reduces the MAE to 1.2 kcal/mol; removal of Et3N results for CPCM/ PAULING lowers the MAE by a negligible 0.01 kcal/mol.

Table 10. Mean Absolute Error (MAE) and Maximum Unsigned Error (MAX) in the Calculated Solution-Phase Free Energies of Activation, ΔG‡ (kcal/mol), at 25 °C for the Base-Catalyzed Enolization of 2-Nitrocyclohexanone M06/6-31+G(d,p)

MAE

MAX

CPCM/Pauling SMD

1.5 1.8

3.3 5.2

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpca.5b04116. Computed free energies of tautomerization, ΔGT, for all molecules discussed including every combination of QM method, basis set, solvation model, and cavity; complete refs 38 and 56. (XLSX)

■

AUTHOR INFORMATION

Corresponding Author

Figure 9. Calculated transition structures for the enolization of 2nitrocyclohexanone in dichloromethane as calculated using M06/631+G(d,p) with (a) SMD and (b) CPCM/PAULING.

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

■

CONCLUSIONS In summary, the gas-phase benchmarking of seven enol-favored tautomerizations was performed using 17 different QM methods and eight basis sets. The G4 method and M06/631+G(d,p) provided the most accurate free energies of tautomerization, ΔGT, with MAE values of 0.95 and 0.71 kcal/mol, respectively. Consequently, the two methods were used in conjunction with three continuum solvent methods and six cavity models to assess their ability to reproduce an additional 23 keto/enol ratios in solution. The G4/PCM/UA0 combination provided the lowest solution-phase MAE in ΔGT with a value of 1.56 kcal/mol. The method excelled in accuracy across all solvents and compound classes tested, namely, aliphatic ketones, β-dicarbonyls, and heterocycles. In the absence of the computational resources necessary for the G4 method, the M06/6-31+G(d,p)/PCM with either the UA0 or UFF radii yielded reasonably accurate ΔGT values in solution. Correct gas-phase reproduction of ΔGT for 2-nitrocyclohexanone proved difficult for the routinely precise composite methods, for example, CBS-QB3 and G4, and calculation of the solvation energy placed the compound functionally somewhere in between the aliphatic ketones and the β-dicarbonyl molecules. To assess solvent effects upon the base-catalyzed keto−enol interconversion of 2-nitrocyclohexanone, the free energies of activation, ΔG‡, were benchmarked with six different bases in five solvents using the M06/6-31+G(d,p) with the CPCM/PAULING and SMD continuum methods. The CPCM method gave an MAE in ΔG‡ of 1.5 kcal/mol,

ACKNOWLEDGMENTS Gratitude is expressed to the National Science Foundation (CHE-1149604) and the Alabama Supercomputer Center for support of this research.

■

REFERENCES

(1) Smith, M. B.; March, J. In March’s Advanced Organic Chemistry, 6th ed.; John Wiley & Sons: Hoboken, NJ, 2007; pp 98−105. (2) Kresge, A. J. Ingold Lecture. Reactive intermediates: carboxylic acid enols and other unstable species. Chem. Soc. Rev. 1996, 25, 275− 280. (3) Kresge, A. J. Generation and study of enols and other reactive species. Pure Appl. Chem. 1991, 63, 213−221. (4) Lovering, E. G.; Laidler, K. J. A System of Molecular Thermochemistry for Organic Gases and Liquis: Part II. Extension to Compounds Containing Sulpher and Oxygen. Can. J. Chem. 1960, 38, 2367−2372. (5) Kresge, A. J. Flash photolytic generation and study of reactive species: from enols to ynols. Acc. Chem. Res. 1990, 23, 43−48. (6) Guthrie, J. P.; Cullimore, P. A. The enol content of simple carbonyl compounds: a thermochemical approach. Can. J. Chem. 1979, 57, 240−248. (7) Capponi, M.; Gut, I.; Wirz, J. The Phenol 2,4-Cyclohexadienone Equilibrium in Aqueous Solution. Angew. Chem., Int. Ed. Engl. 1986, 25, 344−345. (8) Rappoport, Z.; Biali, S. E. Sterically crowded stable simple enols. Acc. Chem. Res. 1988, 21, 442−449. (9) Kaftory, M.; Nugiel, D. A.; Biali, S. E.; Rappoport, Z. Stable simple enols. 23. Solid-state structures of 1-alkyl-2,2-dimesitylethenols. Application of the principle of structural correlation to ring-flip

8731


Article

The Journal of Physical Chemistry A processes in 1,1-diarylvinyl systems. J. Am. Chem. Soc. 1989, 111, 8181−8191. (10) Biali, S. E.; Rappoport, Z. Stable simple enols. 8. Synthesis and keto.dblharw. enol equilibria of the elusive 2,2-dimesitylethanal and 1,2,2-trimesitylethanone. Conformations of 1,2,2-trimesitylethanone and 1,2,2-trimesitylethanol. J. Am. Chem. Soc. 1985, 107, 1007−1015. (11) Caminati, W.; Grabow, J.-U. The C2v Structure of Enolic Acetylacetone. J. Am. Chem. Soc. 2006, 128, 854−857. (12) Srinivasan, R.; Feenstra, J. S.; Park, S. T.; Xu, S.; Zewail, A. H. Direct Determination of Hydrogen-Bonded Structures in Resonant and Tautomeric Reactions Using Ultrafast Electron Diffraction. J. Am. Chem. Soc. 2004, 126, 2266−2267. (13) Moriyasu, M.; Kato, A.; Hashimoto, Y. Kinetic studies of fast equilibrium by means of high-performance liquid chromatography. Part 11. Keto−enol tautomerism of some β-dicarbonyl compounds. J. Chem. Soc., Perkin Trans. 2 1986, 515−520. (14) Becke, A. D. A new mixing of Hartree−Fock and local densityfunctional theories. J. Chem. Phys. 1993, 98, 1372−1377. (15) Xu, X.; Goddard, W. A., III The X3LYP extended density functional for accurate descriptions of nonbond interactions, spin states, and thermochemical properties. Proc. Natl. Acad. Sci. U. S. A. 2004, 101, 2673−2677. (16) Zhao, Y.; Truhlar, D. G. The M06 Suite of Density Functionals for Main Group Thermochemistry, Thermochemical Kinetics, Noncovalent interactions, Excited States, and Transition Elements: Two New Functionals and Systematic Testing of Four M06 Functionals and Twelve Other Functionals. Theor. Chem. Acc. 2008, 120, 215−241. (17) Zhao, Y.; Truhlar, D. G. Density Functionals with Broad Applicability in Chemistry. Acc. Chem. Res. 2008, 41, 157−167. (18) Chai, J.-D.; Head-Gordon, M. Long-range corrected hybrid density functionals with damped atom−atom dispersion corrections. Phys. Chem. Chem. Phys. 2008, 10, 6615−6620. (19) Ernzerhof, M.; Perdew, J. P. Generalized gradient approximation to the angle- and system-averaged exchange hole. J. Chem. Phys. 1998, 109, 3313−3320. (20) Montgomery, J. A., Jr.; Frisch, M. J.; Ochterski, J. W.; Petersson, G. A. A complete basis set model chemistry. VI. Use of density functional geometries and frequencies. J. Chem. Phys. 1999, 110, 2822−2827. (21) Montgomery, J. A., Jr.; Frisch, M. J.; Ochterski, J. W.; Petersson, G. A. A complete basis set model chemistry. VII. Use of the minimum population localization method. J. Chem. Phys. 2000, 112, 6532−6542. (22) Ochterski, J. W.; Petersson, G. A.; Montgomery, J. A., Jr. A complete basis set model chemistry. V. Extensions to six or more heavy atoms. J. Chem. Phys. 1996, 104, 2598−2619. (23) Curtiss, L. A.; Redfern, P. C.; Raghavachari, K. Gn theory. WIREs Comput. Mol. Sci. 2011, 1, 810−825. (24) Dewar, M. J. S.; Zoebisch, E. G.; Healy, E. F.; Stewart, J. J. P. AM1: A New General Purpose Quantum Mechanical Molecular Model. J. Am. Chem. Soc. 1985, 107, 3902−3907. (25) Repasky, M. P.; Chandrasekhar, J.; Jorgensen, W. L. PDDG/ PM3 and PDDG/MNDO: Improved semiempirical methods. J. Comput. Chem. 2002, 23, 1601−1622. (26) Stewart, J. J. P. Optimization of parameters for semiempirical methods V: Modification of NDDO approximations and application to 70 elements. J. Mol. Model. 2007, 13, 1173−1213. (27) Folkendt, M. M.; Weiss-Lopez, B. E.; Chauvel, J. P., Jr.; True, N. S. Gas-phase proton NMR studies of keto-enol tautomerism of acetylacetone, methyl acetoacetate, and ethyl acetoacetate. J. Phys. Chem. 1985, 89, 3347−3352. (28) Tureček, F.; Vivekananda, S.; Sadílek, M.; Polásě k, M. Lactone Enols Are Stable in the Gas Phase but Highly Unstable in Solution. J. Am. Chem. Soc. 2002, 124, 13282−13289. (29) Reichardt, C.; Welton, T. Solvents and Solvent Effects in Organic Chemistry, 4th ed.; Wiley-VCH: Weinheim, Germany, 2011. (30) Karelson, M.; Maran, U.; Katritzky, A. R. Theoretical study of the keto-enol tautomerism in aqueous solutions. Tetrahedron 1996, 52, 11325−11328.

(31) Tomasi, J.; Persico, M. Molecular Interactions in Solution: An Overview of Methods Based on Continuous Distributions of the Solvent. Chem. Rev. 1994, 94, 2027−2094. (32) Cramer, C. J.; Truhlar, D. G. Implicit Solvation Models: Equilibria, Structure, Spectra, and Dynamics. Chem. Rev. 1999, 99, 2161−2200. (33) Zacharias, M. Continuum Solvent Modeling of Nonpolar Solvation: Improvement by Separating Surface Area Dependent Cavity and Dispersion Contributions. J. Phys. Chem. A 2003, 107, 3000− 3004. (34) Cossi, M.; Scalmani, G.; Rega, N.; Barone, V. New developments in the polarizable continuum model for quantum mechanical and classical calculations on molecules in solution. J. Chem. Phys. 2002, 117, 43−54. (35) Cossi, M.; Rega, N.; Scalmani, G.; Barone, V. Achieving LinearScaling Computational Cost for the Polarizable Continuum Model of Solvation. J. Comput. Chem. 2003, 24, 669−681. (36) Barone, V.; Cossi, M.; Tomasi, J. A new definition of cavities for the computation of solvation free energies by the polarizable continuum model. J. Chem. Phys. 1997, 107, 3210−3221. (37) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Universal Solvation Model Based on Solute Electron Density and on a Continuum Model of the Solvent Defined by the Bulk Dielectric Constant and Atomic Surface Tensions. J. Phys. Chem. B 2009, 113, 6378−6396. (38) Frisch, M. J. Gaussian 09, Gaussian, Inc.: Wallingford, CT, 2013 (39) Angelini, G.; Chiappe, C.; De Maria, P.; Fontana, A.; Gasparrini, F.; Pieraccini, D.; Pierini, M.; Siani, G. Determination of the Polarities of Some Ionic Liquids Using 2-Nitrocyclohexanone as the Probe. J. Org. Chem. 2005, 70, 8193−8196. (40) Belova, N. V.; Oberhammer, H.; Girichev, G. V. Tautomeric and Conformational Properties of Methyl Acetoacetate, CH3OC(O)CH2-C(O)CH3: Electron Diffraction and Quantum Chemical Study. J. Phys. Chem. A 2004, 108, 3593−3597. (41) Nowak, M. J.; Lapinski, L.; Fulara, J.; Les, A.; Adamowicz, L. Matrix isolation IR spectroscopy of tautomeric systems and its theoretical interpretation: 2-hydroxypyridine/2(1H)-pyridinone. J. Phys. Chem. 1992, 96, 1562−1569. (42) Beak, P.; Covington, J. B.; White, J. M. Quantitative Model of Solvent Effects on Hydroxypyridine-Pyridone and MercaptopyridineThiopyridone Equilibria: Correlation with Reaction-Field and Hydrogen-Bonding Effects. J. Org. Chem. 1980, 45, 1347−1353. (43) Belova, N. V.; Girichev, G. V.; Shlykov, S. A.; Oberhammer, H. Tautomeric and Conformational Properties of Acetoacetamide: Electron Diffraction and Quantum Chemical Study. J. Org. Chem. 2006, 71, 5298−5302. (44) Piacenza, M.; Grimme, S. Systematic Quantum Chemical Study of DNA-Base Tautomers. J. Comput. Chem. 2004, 25, 83−98. (45) Zhao, Y.; Truhlar, D. G. Design of Density Functionals That Are Broadly Accurate for Thermochemistry, Thermochemical Kinetics, and Nonbonded Interactions. J. Phys. Chem. A 2005, 109, 5656−5667. (46) Galstyan, G.; Knapp, E. W. Computations of 36 tautomer/ isomer equilibria of different lactams. J. Phys. Chem. A 2012, 116, 6885−6893. (47) Powling, J.; Bernstein, H. J. The Effect of Solvents on Tautomeric Equilibria. J. Am. Chem. Soc. 1951, 73, 4353−4356. (48) Yamabe, S.; Tsuchida, N.; Miyajima, K. Reaction Paths of KetoEnol Tautomerization of β-Diketones. J. Phys. Chem. A 2004, 108, 2750−2757. (49) Earle, M. J.; Engel, B. S.; Seddon, K. R. Keto-Enol Tautomerism as a Polarity Indicator in Ionic Liquids. Aust. J. Chem. 2004, 57, 149− 150. (50) Martin, Y. C. Let’s not forget tautomers. J. Comput.-Aided Mol. Des. 2009, 23, 693−704. (51) Geballe, M. T.; Skillman, A. G.; Nicholls, A.; Guthrie, J. P.; Taylor, P. J. The SAMPL2 blind prediction challenge: introduction and overview. J. Comput.-Aided Mol. Des. 2010, 24, 259−279. (52) Soteras, I.; Orozco, M.; Luque, F. J. Performance of the IEFMST solvation continuum model in the SAMPL2 blind test prediction 8732


Article

The Journal of Physical Chemistry A of hydration and tautomerization free energies. J. Comput.-Aided Mol. Des. 2010, 24, 281−291. (53) Ribeiro, R. F.; Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Prediction of SAMPL2 aqueous solvation free energies and tautomeric ratios using the SM8, SM8AD, and SMD solvation models. J. Comput.Aided Mol. Des. 2010, 24, 317−333. (54) Klamt, A.; Diedenhofen, M. Blind prediction test of free energies of hydration with COSMO-RS. J. Comput.-Aided Mol. Des. 2010, 24, 357−360. (55) Meunier, A.; Truchon, J.-F. Predictions of hydration free energies from continuum solvent with solute polarizable models: the SAMPL2 blind challenge. J. Comput.-Aided Mol. Des. 2010, 24, 361− 372. (56) Frisch, M. J. Gaussian 03, Gaussian, Inc.: Wallingford, CT, 2004. (57) Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. Adding explicit solvent molecules to continuum solvent calculations for the calculation of aqueous acid dissociation constants. J. Phys. Chem. A 2006, 110, 2493−2499. (58) Tunon, I.; Rinaldi, D.; Ruiz-Lopez, M. F.; Rivail, J. L. Hydroxide Ion in Liquid Water: Structure, Energetics, and Proton Transfer Using a Mixed Discrete-Continuum ab Initio Model. J. Phys. Chem. 1995, 99, 3798−3805. (59) Schmidt am Busch, M.; Knapp, E.-W. Accurate pKa Determination for a Heterogeneous Group of Organic Molecules. ChemPhysChem 2004, 5, 1513−1522. (60) Schmidt am Busch, M.; Knapp, E.-W. One-Electron Reduction Potential for Oxygen- and Sulfur-Centered Organic Radicals in Protic and Aprotic Solvents. J. Am. Chem. Soc. 2005, 127, 15730−15737. (61) Siani, G.; Angelini, G.; De Maria, P.; Fontana, A.; Pierini, M. Solvent effects on the rate of the keto−enol interconversion of 2nitrocyclohexanone. Org. Biomol. Chem. 2008, 6, 4236−4241. (62) Angelini, G.; De Maria, P.; Chiappe, C.; Fontana, A.; Gasbarri, C.; Siani, G. The Base-Catalyzed Keto−Enol Interconversion of 2Nitrocyclohexanone in Ionic Liquids. J. Org. Chem. 2009, 74, 6572− 6576. (63) Acevedo, O. Simulating Chemical Reactions in Ionic Liquids Using QM/MM Methodology. J. Phys. Chem. A 2014, 118, 11653− 11666. (64) Acevedo, O.; Jorgensen, W. L. Quantum and Molecular Mechanical (QM/MM) Monte Carlo Techniques for Modeling Condensed-Phase Reactions. WIREs Comput. Mol. Sci. 2014, 4, 422−435. (65) Acevedo, O.; Jorgensen, W. L. Advances in Quantum and Molecular Mechanical (QM/MM) Simulations for Organic and Enzymatic Reactions. Acc. Chem. Res. 2010, 43, 142−151.

8733


Benchmarking Continuum Solvent Models for KetoâEnol

Recommend Documents

Benchmarking Continuum Solvent Models for KetoâEnol