Isomer Equilibria of Different Lactams

Article pubs.acs.org/JPCA

Computations of 36 Tautomer/Isomer Equilibria of Different Lactams Gegham Galstyan and Ernst-Walter Knapp* Department of Biology, Chemistry and Pharmacy, Institute of Chemistry and Biochemistry, Freie Universität Berlin, Fabeckstr. 36a, D-14195 Berlin, Germany S Supporting Information *

ABSTRACT: Gas-phase energies of 36 tautomer/isomer pairs of 18 six-membered N-heterocyclic compounds were computed quantum chemically. Among the considered B3LYP, BH&HLYP, BH&HLYP(G), and PW6B95 DFT functionals, the latter two provide accurate tautomer/isomer pair energies with root-mean-square deviations (rmsd) relative to experiments of 0.2 and 0.3 kcal/mol, respectively. Since only few (namely five) experimental data are available, 15 tautomer/isomer pair energies were computed with the very precise QCISD(T)(quadruple-ζ) method serving as reference. Relative to this reference the PW6B95 DFT functional is slightly superior to the BH&HLYP(G) functional, yielding an rmsd of 0.7 and 0.8 kcal/mol, respectively. In contrast to BH&HLYP(G), the PW6B95 DFT functional yields also accurate tautomer/isomer pair energies if zwitterionic structures are involved. The tautomer/isomer pair states possess different amounts of aromaticity. This is characterized by nucleus-independent chemical shift (NICS) values. The tautomer/isomer pair reference energies, from which the energies computed with PW6B95 are subtracted, correlate linearly with the corresponding differences in the NICS values. This correlation is used to construct a correction term for the pair energies computed with PW6B95, yielding tautomer/isomer pair energies with rmsd of 0.3 kcal/mol with respect to the more CPU time demanding QCISD(T)(quadruple-ζ) method.

1. INTRODUCTION Due to its important role in organic and medical chemistry as well as in biochemistry, prototropic tautomerism in heterocyclic compounds has been the subject of many studies over the last 60 years.1−4 Experiments in this field are difficult, since it is necessary to work with ultrapure compounds of unequivocal structures. The high mobility of protons and the often very large energy difference between two tautomeric forms impede measuring ratios of tautomeric equilibria of molecular compounds even when using state of the art physicochemical methods. Therefore, high-quality experimental data on such systems are scarce in particular for lactam−lactim tautomeric pairs,5,6 which made the exploration of tautomeric equilibria an attractive target for theoretical studies in gas7−24 and liquid12,22,25−35 phases. A systematic theoretical study was performed by Grimme et al.23 exploring the energetics of seven tautomer/isomer pairs in vacuum (formamide and pyridone as well as the five nucleic acid bases: cytosine, uracil, thymine, guanine, and isocytosine). They used the QCISD(T) method (quadratic configuration interaction with single and double excitations and triple excitations added perturbatively)36 as a reference to compare with results from QCISD (quadratic configuration interaction with single and double excitations),36 SCS-MP2 (spin-component scaled second-order Møller−Plesset perturbation theory),37 MP2 (second-order Møller−Plesset perturbation theory),38,39 HF (Hartree−Fock theory), and six different DFT (density functional theory)40,41 methods, using triple-ζ basis sets. These © 2012 American Chemical Society

studies revealed that SCS-MP2 performs best in reproducing the QCISD(T) reference energies with a root-mean-square deviation (rmsd) of 0.7 kcal/mol in the gas phase.23 All six tested DFT methods23 tend to underestimate the stability of the lactim relative to lactam tautomer. The authors23 speculated that this behavior was caused by the different aromatic character of heterocyclic rings in the two tautomeric forms (aromatic vs nonaromatic), which cannot be properly described by the same exchange terms in the DFT functional. According to ref 23, the DFT functional with the largest contribution of exact exchange (50%) mimicked the reference energies best among the tested functionals (rmsd of 0.9 kcal/mol in gas phase). For formamide, Grimme et al.23 showed that enlarging the basis set from triple-ζ to quadruple-ζ can improve the predicted tautomeric energy differences by 0.2 kcal/mol and further by another 0.1 kcal/mol using pentuple-ζ basis sets. Leszczynski et al.5 performed a systematic experimental and theoretical study of tautomeric equilibria of 2-pyridinone, 4pyrimidinone, and 2-pyrazinone and their benzo-fused derivatives in the gas phase. Accordingly, benzo-fusion leads to significant stabilization of the lactam form, unless the transformation from lactim to lactam tautomer entails loss of aromaticity, which happens for instance with 3-hydroxyisoquinoline (form B of the tautomer pair 3 in Figure 1). This study also Received: March 16, 2012 Revised: May 22, 2012 Published: May 31, 2012 6885

dx.doi.org/10.1021/jp302569g | J. Phys. Chem. A 2012, 116, 6885−6893

The Journal of Physical Chemistry A

Article

0(W) Figure 1. Structures of 36 tautomer/isomer pairs (V/W) are shown, for which the equilibrium energies ΔE(V/W) = E0(V) electronic − Eelectronic are computed in the gas phase. The tautomer/isomer pairs of two molecular species named as V and W are denoted by (V/W). To address the isomer pairs individually, they are enumerated with bold numbers. The label R indicates that a computed reference pair energy [ΔE(V/W)(ref)] based on the QCISD(T)(q-ζ) level is available. Analogous tautomer/isomer pairs are grouped together, yielding nine different families denoted by the letters a−i. The letters A, A, C, or Y denote the tautomer/isomer states with deprotonated oxygens (i.e., oxo-groups) or nitrogens which are not part of the conjugated ring system. The letter Z denotes zwitterionic states of tautomer/isomer pairs. Two types of prototropic tautomerisms are considered: pairs of lactam/lactim (1−6, 8− 10, 12, 14, 16, 18, 20, 22, 24, 28−36) and of amine/imine (21, 23, 25) type. Two types of prototropic isomerisms are considered: positional (7, 11, 13, 15, 17, 19, 27) and rotational (26). The states denoted A (or A) are always the oxo form of the lactam/lactim tautomer pairs, while B (or D) denote the respective hydroxy form. If a compound contains two nitrogens in the ring, positional isomerism is possible; i.e., a hydrogen can be attached to either one of the nitrogens. Such isomer pairs are denoted by (A/C) or (A/Z) (see for instance the pairs 7 and 15). If compound A possesses an amino-substituent at a ring position adjacent to nitrogen in the ring, amine/imine tautomerism is possible and the Y state refers to the imino form of the corresponding amine/imine pair (see for instance 21). The letters with upper bar describe tautomer pairs involving two independent pairs of equivalent functional groups (two ring nitrogens and two oxo/hydroxyl groups) that can undergo lactam/lactim tautomerism (see for instance 29). A̅ represents compounds with two oxo groups, D̅ with two hydroxyl groups, while B̅ and C̅ refer to compounds with one oxo and one hydroxy group (see family i). Within the same family compounds denoted by the same one-letter symbol in the same color are identical compounds. CAS names of all structures are given in Table S1 of the Supporting Information.

The influence of different factorsincluding also π-electron distributionon gas-phase tautomeric equilibria of uracil and related azaphenols was studied theoretically in ref 24 on the basis of AM1 computations42 and HOMED (harmonic oscillator

clearly showed that the presence of a second aromatic ring reduces the aromaticity difference between the two tautomeric forms as compared with their single ring analogues (see for instance tautomer pairs 1 and 4 in Figure 1). 6886



Article

model of electron delocalization) indices.43 According to that study, the preference of a tautomeric state is mainly determined by the stability of the relevant molecular groups carrying the hydrogen. Another important factor affecting tautomeric preference is the presence of aza-nitrogens in the 2,4- or 2,6positions of a phenol ring. Similarly to polarity and intramolecular H-bond-like interactions, aromaticity had only a minor influence on the preferred tautomeric state of the compounds studied in ref 24. It was proposed that the discrepancies between tautomerization energies of DNA bases obtained by quantum chemical DFT and QCISD(T)/TZV(2df,2dp) computations are connected to the problem of DFT functionals describing aromatic compounds with insufficient precision.23 They tend to underestimate the stability of the hydroxy form relative to the oxo form of a tautomer pair, which was shown in ref 23 to originate from the type of exchange terms used in the DFT functionals. In the present study we computed tautomerization energies for 18 different molecular compounds (shown in Figure 1) in the gas phase, which can serve as the basis to compute these energies in aqueous solution. To calculate the electronic energy difference 0(W) ΔE(V/W) = E0(V) electronic − Eelectronic of the tautomer/isomer pair (V/ W) we used the quantum chemical method QCISD(T) with double-ζ and for relatively small molecules also triple-ζ and quadruple-ζ basis sets [QCISD(T)(d/t/q-ζ)], as well as DFT with the functionals B3LYP,44−47 BH&HLYP,45,48 BH&HLYP(G) (modified BH&HLYP from Gaussian,49 as described in the methods part) and PW6B9550 using quadruple-ζ basis set. Furthermore, we compared the calculated energies of tautomer equilibria with the available experimental data for five compounds. For those tautomer/isomer pairs for which no measured energies are available, we used for the compounds with relatively small number of atoms the energies computed with QCISD(T)(q-ζ) as reference, denoted as ΔE(V/W)(ref) for the tautomer/isomer pair (V/W). When no energies with the QCISD(T)(q-ζ) method are available, we used the DFT functional PW6B95 as reference as explained below. To explore how the DFT approach can deal with aromatic compounds, we also computed nucleus-independent chemical shift indexes of the considered compounds, to find correlations between aromaticity and systematic deviations of tautomer/isomer pair energies obtained with the DFT method. These correlations were used to construct an empirical correction to obtain more precise tautomer/isomer pair energies with the PW6B95 functional. The resulting accuracy is 0.3 kcal/mol rmsd relative to the tautomer/isomer energies computed with the precise quantum chemical method QCISD(T)(t-ζ) and where possible QCISD(T)(q-ζ). Tautomer/isomer equilibria are difficult to measure. They have often been computed quantum chemically for the gas phase only. However, for practitioners, tautomer/isomer energies in aqueous solution are relevant. So far, precise results for the gas phase were obtained with quantum chemical approaches using the enormously CPU demanding QCISD(T) methods with triple-ζ basis functions. The present study provides a detailed analysis of five DFT functionals. It is demonstrated in this study that, combined with electrostatic energy computations, appropriate DFT functionals yield very precise tautomer/isomer energies in aqueous solution. This allows one to compute tautomer/isomer energies with a moderate amount of CPU time.

2. METHODS 2.1. Geometry Optimization and Vibrational Energies. The quantum chemical program Jaguar v7.751 was used for all tautomer/isomer states (individual states are denoted by X) of the 18 compounds considered in this study to optimize the geometries and to compute the zero point vibrational energies (X) (E(X) ZPVE) as well as the thermal vibrational free energies (Gvib (T)) 44−47 using the B3LYP DFT functional with Pople’s split-valence 6-31G** basis set.52−54 These geometries and vibrational energies are combined with ground-state electronic energies computed with different quantum chemical methods as explained below. In the present study, we refrained from using frequency rescaling factors to compute the vibrational energies, since their actual influence on the tautomer/isomer pair energies was small. All vibrational energies have been performed at T = 298.15 K, but for the five compounds (1, 4, 8−10, Figure 1), where experimental tautomerization energies in the gas phase are available, these energies were also computed at the temperatures at which the measurements of tautomerization equilibria were performed. 2.2. Ground-State Electronic Energies. The ground state 0(X) electronic energies Eelectronic of tautomeric species X were obtained from single point energy calculations with Jaguar v7.7,51 carried out for the optimized geometries, using B3LYP, 4 4 − 4 7 BH&HLYP, 4 5 , 4 8 BH&HLYP(G), 4 9 and PW6B9550 DFT functionals with the relatively large Dunning’s correlation-consistent basis set cc-pVQZ(-g).55 The BH&HLYP(G) functional, a modified version of BH&HLYP implemented in Gaussian 0949 EBH&HLYP(G) = 0.5EexHF + 0.5EexLSDA + 0.5ΔEexBecke88 + ECLYP

(1)

contains an additional nonlocal exchange term 0.5ΔEBecke88 as ex compared to the original “half and-half” functional proposed by Becke.48 For the computationally more expensive electronic energy computations with the QCISD(T) (quadratic configuration interaction with single and double excitations and triple excitations added perturbatively) method,36 we used Gaussian 09.49 Such computations where performed for all 36 tautomer/ isomer pairs using the relatively small 6-31G** basis set52−54 [short notation QCISD(T)(d-ζ)]. For the tautomer/isomer pairs (1, 5, 8, 10, 11, 14, 15, 20−36) involving the smaller compounds only, we performed also the more CPU time demanding QCISD(T) computations using Dunning’s correlation-consistent basis sets cc-pVTZ and cc-pVQZ,55 abbreviated as QCISD(T)(t-ζ) and QCISD(T)(q-ζ), respectively. For the tautomer pairs 4 and 9 involving two rings, we also performed computations with QCISD(T)(t-ζ), since for those pairs experiments are available. 2.3. Computation of Gas-Phase Energies. To compare quantum chemically computed energies with corresponding experimental values of gas-phase tautomer/isomer equilibria, the gas-phase energies G(X) gas (T) of the individual tautomeric species X = A, B are computed as the sum of ground state electronic energy (X) (E0(X) electronic), zero point vibrational energy (EZPVE), and thermal vibrational free energy (G(X) (T)) vib (X ) (X ) 0(X ) (X ) Ggas (T ) = Eelectronic + EZPVE + Gvib (T )

(2)

(A) (B) yielding the energy difference ΔG(A/B) gas (T) = Ggas (T) − Ggas (T) for the tautomer/isomer pair (A/B).

6887



Article

2.4. Chemical Shift and Aromaticity. The nucleusindependent chemical shift (NICS) tensor can be used to monitor magnetic shielding in the neighborhood of a molecular compound. It was first used in ref 56 to characterize the degree of aromaticity of molecular compounds. Different definitions of scalar NICS values derived from the NICS tensor are in use. The NICS(0) index is for instance the total isotropic nuclear magnetic shielding (with reversed sign) in the center of the ring plane (geometric center of the heavy atoms of the ring). The more general NICS(d) index refers to a position d Å above the center of the ring plane. The recently introduced NICS(0)πzz57,58 index (being superior in conception and performance), measures the magnetic shielding in the ring center for the out-of-plane zzcomponent of the NICS tensor, which is due to the π-molecular orbitals only. In benchmark tests against 75 aromatic stabilization energies, NICS(0)πzz was shown59 to be the method of choice among other variants of NICS to characterize the aromaticity of planar aromatic rings. However, the computation of the NICS(0)πzz index requires localized or canonical molecular orbital (LMO or CMO) dissection provided by natural bond orbital analysis, which can only be performed with special software (for more details see ref 59). The more readily accessible NICS(1)zz aromaticity index, which we used in the current work, is a useful alternative to the NICS(0)πzz index.59 NICS(1)zz considers the out-of-plane zz-component of the isotropic NICS tensor. It is computed by a gauge-independent atomic orbital approach,60 based on wave functions computed with DFT PW91PW91/6-311+G** 61,62 for all tautomer/ isomer pairs listed in Figure 1. For these computations, we employed the Gaussian 09 program.49 It should be noted, that the NICS(1)zz(X) value of a tautomer/isomer species X containing two rings is computed as sum of NICS(1)zz values of the individual rings NICS(1)zz(X ) = NICS(1)zz a(X ) + NICS(1)zz b(X )

Figure 2. Correlation diagram of computed and measured5 free energies of tautomerization in the gas phase (ΔG(A/B) = G(A) − G(B), where A and B refer to the tautomeric species with the hydrogen atom at the nitrogen and oxygen atom, respectively) for five tautomer pairs (indicated by numbers 1, 4, 8−10; see Figure 1). The computed pair energies are based on QCISD(T)(d/t/q-ζ) and the DFT functionals B3LYP, BH&HLYP, BH&HLYP(G), and PW6B95. The numerical values are listed in Table 1. According to ref 5 the experimental error is ±0.12, ±0.36, ±0.24, ±0.24, and ±0.07 kcal/mol for pairs 1, 4, 8, 9, and 10, respectively.

Table 1. Free Energies of Tautomerization ΔG(A/B) gas (T) (kcal/ mol) in the Gas Phasea

(3)

where NICS(1)zza(X) and NICS(1)zzb(X) are the NICS(1)zz values for the main ring (denoted by a) and the fused ring (denoted by b), respectively.

3. RESULTS AND DISCUSSION 3.1. Comparison with Measured Tautomer Pair Energies. Four different quantum chemical methods based on DFT are considered to compute equilibrium energies in the gas phase for five tautomer pairs (1, 4, 8−10, Figure 1), where measured values are available.5 Results are shown in Figure 2 and Table 1. Among the four DFT functionals considered in this study [B3LYP, 44−47 BH&HLYP, 45,48 PW6B95, 50 and BH&HLYP(G)], B3LYP and BH&HLYP poorly reproduce the measured tautomeric energy differences. Results computed with the two DFT functionals PW6B95 and BH&HLYP(G) agree well with experiment, as do results obtained with QCISD(T)(d/t/q-ζ). The latter two very CPU time demanding methods [QCISD(T)(t/q-ζ)] were applied only for tautomer/ isomer pairs involving only a relatively small number of atoms (1, 5, 8, 10, 11, 14, 15, 20−23, 29−32). According to the results for the five compounds where measured tautomer pair energies are available, the tautomer pairs 1 and 8 are more stable in the hydroxy forms (B), while the pairs 4, 9, and 10 are more stable in the oxo forms (A). Root mean square deviations (rmsd) of computed versus measured energies of these five tautomer pairs are 1.0, 0.8, 0.3, 0.2, 0.5, 0.4, and 0.2 kcal/mol for B3LYP, BH&HLYP, PW6B95, BH&HLYP(G),

method

1 (A/B)

4 (A/B)

8 (A/B)

9 (A/B)

10 (A/B)

B3LYP B3LYP-LOCc BH&HLYP PW6B95-corrd PW6B95 BH&HLYP(G) QCISD(T)(d-ζ) QCISD(T)(t-ζ) QCISD(T)(q-ζ) expe Texp (K)e

−0.58 1.08 1.40 0.76 0.31 0.81 1.07 1.23 0.97 0.69 340

−4.47 −2.81 −3.12 −3.29 −3.85 −3.77 −3.14

0.53 2.19 2.66 1.96 1.49 2.07 1.71 2.04 1.82 1.91 360

−3.76 −2.10 −2.48 −2.60 −3.18 −3.11 −2.67

−1.76 −0.10 0.19 −0.44 −0.88 −0.44 −0.54 −0.22 −0.44 −0.57 400

−4.06 450

−3.35 450

rmsdb 1.0 0.8 0.8 0.5 0.3 0.2 0.5 0.4 0.2 0.0

a The energies are given for the five tautomer pairs (1, 4, 8-10, defined in Figure 1), where measured values are available. bThe rmsd values are root-mean-square deviations of computed relative to the five measured energies. cIt should be noted that the parameters we used here for localized orbital corrections of energies of tautomeric species were originally fitted for the cc-pVQZ basis set against the G3 data set.64 An appropriate set of parameters obtained exactly for the ccpVQZ(-g) basis set, used in DFT computations in the current paper, should result in a better match with the measured pair energies. dOwn method using PW6B95 with empirical correction based on eq 5 adjusted relative to the pair energies from QCISD(T)(q-ζ). Considering the larger number of compounds in Table 2, where no measured values are available, the performance of PW6B95-corr is better than that of PW6B95 and BH&HLYP(G) (see Table 2). e Measured energies and temperatures are taken from ref 5.

QCISD(T)(d-ζ), QCISD(T)(t-ζ), and QCISD(T)(q-ζ), respectively (Table 1). All DFT calculations were performed with the Dunning’s correlation-consistent quadruple-ζ valence basis set55 cc-pVQZ(-g) as implemented in Jaguar v7.7.51 The tautomer pair energies computed with the B3LYP (BH&HLYP) functional systematically underestimate (overestimate) the measured pair energies (see Figure 2). A main reason for the weak performance of the DFT functionals B3LYP and 6888



Article

0(W) Figure 3. Diagram showing the deviation of computed electronic energy differences ΔE(V/W) = E0(V) electronic − Eelectronic of tautomer/isomer pairs (V/W) from the computed reference pair energies ΔE(V/W)(ref) that are based on QCISD(T)(q-ζ) (marked with R in Figure 1). For tautomer/isomer pairs not listed here, no theoretical reference pair energy is available. Abscissa: tautomer/isomer pairs numbered according to Figure 1. Ordinate: deviation of electronic pair energies ΔE(V/W) obtained with the DFT functionals BH&HLYP, PW6B95, and BH&HLYP(G) and the QCISD(T)(d-ζ) method relative to the reference energies.

Table 2. Ground State Electronic Energy Differences ΔE0(X) electronic (kcal/mol) for 15 Tautomer/Isomer Pairs V/W from Figure 3, Computed Using Different Quantum Chemical Methods V/W

B3LYP

BH&HLYP(G)

B3LYP-LOCa

BH&HLYP

PW6B95

PW6B95-corr

QCISD(T)(d-ζ)

QCISD(T)(t-ζ)

QCISD(T)(q-ζ)

1 (A/B) 5 (A/B) 8 (A/B) 10 (A/B) 11 (A/C) 14 (A/B) 15 (A/Z) 20 (A/B) 21 (A/Y) 22 (A/B) 23 (A/Y) 29 (A̅ /B̅ ) 30 (A̅ /C̅ ) 31 (B̅ /D̅ ) 32 (C̅ /D̅ ) rmsdb

−0.56 1.24 0.60 −1.66 −9.33 −7.79 −18.56 −0.96 −1.73 1.00 −5.35 −11.24 −11.53 −1.67 −1.37 1.31

0.83 2.94 2.14 −0.34 −9.54 −6.60 −21.58 0.30 −1.62 2.04 −5.84 −11.00 −11.75 −0.78 −0.03 0.83

1.10 2.90 2.26 0.00 −9.33 −6.13 −18.70 0.70 −2.60 2.66 −6.22 −9.58 −9.87 −0.01 0.29 0.74

1.42 3.55 2.73 0.29 −9.50 −6.12 −21.00 0.85 −1.77 2.84 −5.74 −10.33 −10.91 −0.06 0.53 0.72

0.33 2.26 1.56 −0.77 −9.39 −6.99 −19.00 −0.02 −1.91 1.86 −5.63 −10.62 −10.95 −0.77 −0.43 0.67

0.78 2.73 2.03 −0.33 −9.52 −6.54 −19.05 0.35 −2.08 2.15 −5.87 −10.44 −10.83 −0.42 −0.03 0.53

1.09 3.33 1.78 −0.44 −10.12 −6.35 −20.75 1.33 0.25 2.01 −5.77 −10.93 −12.42 −0.93 0.56 0.69

1.25 3.65 2.11 −0.12 −9.43 −6.01 −19.47 1.41 −0.22 2.32 −5.45 −10.04 −11.12 −0.28 0.81 0.26

0.99 3.33 1.89 −0.34 −9.19 −6.19 −19.07 0.96 −0.48 2.13 −5.39 −10.20 −11.09 −0.44 0.44 0.00

a The correction parameters fitted to the G3 data set for the cc-pVTZ basis set were used.64 These parameters were applied to the B3LYP energies obtained with the cc-pVQZ(-g) basis set and therefore may provide less precise results. bThe rmsd values are root-mean-square deviations of the computed pair energies relative to the QCISD(T)(q-ζ) reference energies.

BH&HLYP in this application may be the lack of balance of the exchange energy terms for aromatic compounds. In the Table 1 we also present B3LYP energies that are empirically corrected using the localized orbital correction scheme (abbreviated as B3LYP-LOC).63 The parameters needed for the B3LYP-LOC approach are for instance available for the cc-pVTZ basis set using G3 as set of chemical compounds.64 This basis set is closest to cc-pVQZ(-g) used in the present study and therefore used here. Using B3LYP-LOC the rmsd between the computed and measured five tautomer energies diminishes from 1.0 to 0.8 kcal/mol. The rmsd may become even smaller, if the parameters of the B3LYP-LOC approach would be optimized using the proper basis set. The correction energies for the different types of tautomer/isomer pairs resulting from the B3LYP-LOC approach are listed in the Supporting Information (Table S2). (T) of According to Figure 2 the free energies ΔG(V/W) gas tautomer pairs (V/W) calculated with the electronic energies based on QCISD(T)(d-ζ) agree better with the experimental free energies than those calculated on the basis of the

QCISD(T)(t-ζ) electronic energies and nearly as good as the results obtained with QCISD(T)(q-ζ). This effect is due to fortuitous error compensation in the difference of the pair energies, as can be seen in Figure 3. Therefore, we preferentially compare the tautomer/isomer pair energies computed with the DFT methods with results obtained with QCISD(T)/cc-pVQZ [QCISD(T)(q-ζ)], which we use as reference energy. A closer analysis of Figure 2 reveals small but systematic deviations between the pair energies from experiment and computations using the reference method QCISD(T)(q-ζ). 3.2. Comparing DFT Pair Energies with Reference Values from QCISD(T)(q-ζ). Figure 3 provides a comparison of computed electronic energy differences (ΔE(V/W)) of 15 tautomer/isomer pairs (1, 5, 8, 10, 11, 14, 15, 20−23, 29−32; see Figure 1). Eleven pairs are oxo/hydroxy tautomers (1, 5, 8, 10, 14, 20, 22, 29−32), two are amino/imino tautomers (21, 23), and two are positional isomer pairs (11, 15). The energies of all 15 tautomer/isomer pairs are calculated with all four DFT functionals (B3LYP, BH&HLYP, BH&HLYP(G), and PW6B95) and compared with more precise quantum chemical 6889



Article

S4 of the Supporting Information. These three pairs possess less aromatic character in the fused ring when they are in states 3 (A), 7 (C), and 18 (A), compared to their corresponding states 3 (B), 7 (A), and 18 (B), respectively. For all other tautomer/isomer pairs with structures involving a fused ring (2, 4, 6, 9, 12, 13, 16, 17, 19, 24−28), the changes in NICS(1)zzb values within a pair are smaller than ±2.9 ppm. The deviations between tautomer/ isomer pair energies computed with PW6B95 and BH&HLYP(G) can be explained in part by the different amounts of exact exchange energy (28% for PW6B95 and 50% for BH&HLYP(G)) inherent in the two DFT functionals. Hence, we recommend the use of the BH&HLYP(G) rather than the PW6B95 DFT functional, if there are larger changes in aromaticity between the two states of a tautomer/isomer pair, which is in agreement with other work.23 The isomer pairs 15, 17, and 19 involve a zwitterionic state labeled by the letter Z (Figure 1). For these three isomer pairs we observe severe systematic deviations of the pair energies obtained with the PW6B95 and the BH&HLYP(G) functionals. Using PW6B95 instead of BH&HLYP(G), the pair energies of 15, 17, and 19 are larger by 2.6, 3.1, and 2.3 kcal/mol, respectively. Values of individual pair energies are given in Table S1 of the Supporting Information. For the isomer pair 15 the pair energy computed with PW6B95 agrees much better with the QCISD(T)(q-ζ) reference energy than by using BH&HLYP(G) (see Figure 3). The presence of a fused ring in the isomer pair 17 increases the discrepancy of the pair energies computed with these two DFT functionals from 2.6 to 3.1 kcal/mol as compared to the isomer pair 15 involving just the main ring (see Figure 1). Interestingly, this discrepancy decreases to 2.3 kcal/mol from isomer pair 17 to 19, even though both isomer pairs involve a fused ring. The differences of the NICS(1) zz values (ΔNICS(1)zz(A/Z) = NICS(1)zz(A) − NICS(1)zz(Z)) are small for all three isomer pairs, varying between 2 and 3 ppm only (see Table S4 in the Supporting Information). Therefore, changes in aromaticity between the two states of those three isomer pairs (15, 17, and 19) cannot be made responsible for the large deviations of pair energies obtained with the DFT functionals BH&HLYP(G) and PW6B95. On the basis of the results for the isomer pair 15 involving a zwitterionic state (Figure 3), we can conclude that PW6B95 generally outperforms BH&HLYP(G) in computing pair energies of isomer pairs involving zwitterions. The zwitterionic states are known to be stable in aqueous solution. But there is no clear experimental evidence that the three compounds (15, 17, and 19) assume a stable zwitterionic form in the gas phase. Computations of transition state energies were performed with the PWB6K functional (appropriate for kinetic studies)50 (more details are given in the Supporting Information). They yielded large activation barriers (70.0, 70.2, and 69.0 kcal/mol for the isomer pairs 15, 17, and 19, respectively; see Figure S1, Supporting Information) and electronic energies, which are larger in the zwitterionic states (Z) by 20.5, 22.0, and 18.1 kcal/mol for the pairs 15, 17, and 19, respectively. However, the gas-phase energies of these isomer pairs will be needed to compute the corresponding pair energies in aqueous solution using the thermodynamic cycle technique to connect gas and aqueous phases.65 In view of the scarce experimental data to compare with (see Figure 2), it is difficult to judge which one of the two DFT functionals (PW6B95, BH&HLYP(G)) is more appropriate to compute tautomer/isomer pair energies in the gas phase. However, if we trust in the accuracy of the pair energies obtained with the reference values calculated by QCISD(T)(q-ζ), the

computations using the QCISD(T)/cc-pVQZ method (QCISD(T)(q-ζ)) for pairs 1, 5, 8, 10, 11, 14, 15, 20−23, 29−32. In Figure 3, the tautomerization energy differences obtained with the three DFT functionals, BH&HLYP, BH&HLYP(G), and PW6B95, and the QCISD(T)(d-ζ) method, relative to the reference pair energies ΔE(V/W)(ref), are displayed for the 15 tautomer/isomer pairs, where reference energies are available. Pair energies obtained with B3LYP are not shown because of the large discrepancies from the reference energies. Numerical values of the pair energies displayed in Figure 3 are given in Table 2. According to Figure 3, the DFT functionals PW6B95 and BH&HLYP(G) yield tautomer/isomer pair energies with an rmsd of 0.67 and 0.83 kcal/mol, respectively, relative to the computed reference tautomer/isomer pair energies ΔE(V/W)(ref) that are based on QCISD(T)(q-ζ). These reference pair energies are only available for a subset of 15 tautomer/isomer pairs marked by the letter R in Figure 1. All these compounds possess a single aromatic ring only. The DFT functional PW6B95 severely underestimates (by more than 1.0 kcal/mol) the tautomer pair energies of the pairs 5, 20, and 21 by 1.1, 1.0, and 1.4 kcal/mol, respectively. BH&HLYP(G) severely underestimates (by more than 1.1 kcal/mol) the tautomer/isomer pair energies of the pairs 15 and 21 by 2.5 and 1.1 kcal/mol, respectively. If these specific tautomer/isomer pairs are not considered, the rmsds of the pair energies obtained with the PW6B95 and BH&HLYP(G) DFT functionals fall down from 0.67 and 0.83 to 0.47 and 0.45 kcal/ mol, respectively. BH&HLYP shows an rmsd of 0.72 kcal/mol, slightly smaller than the rmsd obtained with BH&HLYP(G) (Figure 3). This DFT functional severely underestimates the energies of the tautomer/isomer pairs 15 and 21, by 1.9 and 1.3 kcal/mol, respectively. Omitting these two tautomer/isomer pairs lowers the corresponding rmsd from 0.72 to 0.42 kcal/mol. The B3LYP DFT functional (results not shown in Figure 3) systematically underestimates the reference pair energies ΔE(V/W)(ref) for most of the tautomer/isomer pairs considered in this work (see Table 2). The B3LYP-LOC63,64 energies show significant improvement over the B3LYP results. It should be noted that the pair energy of the isomer pair 15 involving a zwitterionic state is well-reproduced only by the PW6B95 functional. Overall, PW6B95 performs better than the other considered DFT functionals in computing the tautomer/isomer pair energies. BH&HLYP becomes better if cases with severe deviations are not counted. However, its performance is worse than that of QCISD(T)(d-ζ) and QCISD(T)(t-ζ) in reproducing the QCISD(T)(q-ζ) reference pair energies (see also Figure 2). Hence, we can safely assume that the BH&HLYP(G) and PW6B95 DFT functionals are generally more appropriate to compute gas-phase energies of tautomeric/isomeric equilibria for the molecular systems considered in the present study. This is particularly true, since the pair energies obtained with these two functional agree better with the available five measured pair energies (see Figure 2). The tautomer/isomer pair energies computed with PW6B95 do not always agree well with the pair energies obtained with the BH&HLYP(G) functional. They deviate considerably for the tautomer/isomer pairs 3, 7, and 18 by +1.4, +1.3, and +1.1 kcal/ mol, respectively (for the individual energy values see Table S3 in the Supporting Information). These three tautomer/isomer pairs exhibit considerable changes of aromatic character of the fused rings between the pair states corresponding to NICS(1)zzb changes of 12.1, −8.5, and 11.9 ppm for the pairs 3, 7, and 18, respectively. A detailed list of NICS(1)zzb values is given in Table 6890



Article

DFT functional PW6B95 yields on average better agreement with the reference pair energies than BH&HLYP(G) (Figure 3). 3.3. Correlation between Aromaticity and Accuracy of Computed Pair Energies. For most of the 36 considered tautomer/isomer pairs the degree of aromaticity differs between the two pair states. This can be deduced by analyzing the schematic structures displayed in Figure 1. For instance, state A of the tautomer pair 1 is less aromatic than state B. However, there are also tautomer/isomer pairs where the aromaticity practically does not differ between the two states of a pair, as is the case for the pairs 11, 13, 26, and 27. This is corroborated by the computed absolute ΔNICS(1)zz values, which vanish for the tautomer/isomer pair 27 (0.3 ppm) and are small for 26, 11, and 13 (2.2, −3.8, −3.5 ppm). Figure 4 shows the correlation between ΔNICS(1)zz values and deviations of tautomer/isomer pair energies

ΔE(V/W)(PW6B95‐corr) = ΔE(V/W)(PW6B95) + 0.03ΔNICS(1)zz (kcal/mol) (5)

where eq 3 was employed and the ΔNICS(1)zz values are given in ppm units. The rmsd between the corrected tautomer/isomer pair energies based on our approach using eq 5 and the 15 computed energies based on the precise QCISD(T)(q-ζ) approach or the five measured energies are 0.53 or 0.50 kcal/mol, respectively. Not considering the tautomer pair 21 in the correlation diagram (Figure 4), the former rmsd drops down from 0.53 to 0.34 kcal/ mol.

4. CONCLUSIONS The gas-phase energies of 36 tautomer/isomer pairs of 18 different lactams have been computed with different DFT functionals and compared with high-precision quantum chemical computations based on QCISD(T)(q-ζ) as well as with experimental values if available. These energies are necessary to compute isomerization energies for these compounds in aqueous solution, which is important information to understand the interactions of bases in DNA and RNA. The accuracy of the PW6B95 functional50 is slightly higher than that of the BH&HLYP(G) functional.49 Employing quadruple-ζ basis sets, both PW6B95 and BH&HLYP(G) are more successful in computations of tautomer/isomer equilibrium energies of the studied lactams, as compared with the QCISD(T)(d-ζ) method. In contrast to BH&HLYP(G), the PW6B95 DFT functional can be safely used to compute also isomer pair energies involving zwitterions. Errors in tautomer/isomer pair energies computed with the PW6B95 DFT functional were correlated with the difference between NICS(1)zz values characterizing the change in aromaticity between tautomeric/isomeric pair states. From this correlation a correction term to the PW6B95 based pair energies was proposed. It allows computing keto−enol, amide− imide tautomerization energies of N-heterocyclic ketones and their derivatives with high precision yielding an rmsd of 0.27 kcal/mol relative to the very precise but CPU time demanding QCISD(T)(q-ζ) method. Combined with the implemented electrostatic energy approach, the quantum chemical PW6B95 DFT functional and its correction terms open the possibility to compute tautomer/ isomer equilibrium energies in aqueous solution with high precision at moderate CPU cost. The more precise quantum chemical procedures based for instance on QCISD(T) with triple or even more with quadruple-ζ basis sets are prohibitively expensive for larger molecules involving two aromatic rings. Such molecules are also considered in this study. Finally, we make a remark of caution. Although we consider a relatively large set of tautomer/isomer equilibria, the compounds are all lactams and thus of similar type. Hence, caution should be applied when transferring the recommended procedures to different types of compounds.

Figure 4. Correlation between ΔNICS(1)zz values and deviations of tautomer/isomer pair energies [ΔE(V/W)(PW6B95)] calculated with the PW6B95 DFT functional from the reference pair energies [ΔE(V/W)(ref)] based on QCISD(T)(q-ζ).

[ΔE(V/W)(PW6B95)] calculated with the PW6B95 DFT functional from the reference pair energies [ΔE(V/W)(ref)] based on QCISD(T)(q-ζ). Excluding the pair 21, where the pair energy computed with the PW6B95 DFT functional deviates the most from the reference energy, we obtain approximately a linear dependence of ΔΔE(V/W)(PW6B95) = ΔE(V/W)(PW6B95) − ΔE(V/W)(ref) as a function of the ΔNICS(1)zz values with an rmsd of 0.3 kcal/mol, if the pair 21 is ignored. Hence, we can write ΔΔE(V/W)(PW6B95) = − 0.20 − 0.03ΔNICS(1)zz (kcal/mol)

(4)

where the ΔNICS(1)zz values are given in ppm units yielding energies in units of kcal/mol. For a perfect correlation between ΔΔE(V/W) and the ΔNICS(1)zz values, the constant term should vanish. In principle, one can use ΔΔE(V/W), eq 4, to correct the tautomer/isomer pair energies computed with the PW6B95 DFT functional. However, we know from the analysis of the experimental data (see Figure 2) that the reference pair energies based on QCISD(T)(q-ζ) deviate systematically from the measured pair energies by ΔE(ref) = 0.2 kcal/mol. Hence, to get estimates of tautomer/isomer pair energies computed with the PW6B95 DFT functional, which agree well with experiments, we should correct the pair energy ΔE(V/W)(PW6B95) according to

■

ASSOCIATED CONTENT

S Supporting Information *

CAS (Chemical Abstracts Service) names of considered compounds, optimized geometries, vibrational energies and ground state electronic energies for all tautomer species, all methods used, examples of all inputs for all kinds of computations in all programs used in this work, and the full 6891



Article

(21) Géza, F. J. Phys. Chem. A 2002, 106, 1381−1390. (22) Hanus, M.; Ryjácek, F.; Kabelác, M.; Kubar, T.; Bogdan, T. V.; Trygubenko, S. A.; Hobza, P. J. Am. Chem. Soc. 2003, 125, 7678−7688. (23) Piacenza, M.; Grimme, S. J. Comput. Chem. 2004, 25, 83−99. (24) Raczynska, E. D.; Zientara, K.; Kolczynska, K.; Stepniewski, T. J. Mol. Struct. (THEOCHEM) 2009, 894, 103−111. (25) Skillman, A. G.; Geballe, M. T.; Nicolls, A. J. Comput. Aided Mol. Des. 2010, 24, 257−258. (26) Soteras, I.; Orozco, M.; Luque, F. J. J. Comput. Aided Mol. Des. 2010, 24, 281−291. (27) Nicholls, A.; Wlodek, S.; Grant, J. A. J. Comput. Aided Mol. Des. 2010, 24, 293−306. (28) Klimovich, P. V.; Mobley, D. L. J. Comput. Aided Mol. Des. 2010, 24, 307−316. (29) Ribeiro, R. F.; Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. J. Comput. Aided Mol. Des. 2010, 24, 317−333. (30) Ellingson, B. A.; Skillman, A. G.; Nicholls, A. J. Comput. Aided Mol. Des. 2010, 24, 335−342. (31) Kast, S. M.; Heil, J.; Güssregen, S.; Schmidt, K. F. J. Comput. Aided Mol. Des. 2010, 24, 343−353. (32) Kast, S. M.; Heil, J.; Güssregen, S.; Schmidt, K. F. J. Comput. Aided Mol. Des. 2010, 24, 355 (E). (33) Klamt, A.; Diedenhofen, M. J. Comput. Aided Mol. Des. 2010, 24, 357−360. (34) Meunier, A.; Truchon, J. F. J. Comput. Aided Mol. Des. 2010, 24, 361−372. (35) Purisima, E. O.; Corbeil, C. R.; Sulea, T. J. Comput. Aided Mol. Des. 2010, 24, 373−383. (36) Pople, J. A.; Head-Gordon, M.; Raghavachari, K. J. Chem. Phys. 1987, 87, 5968−5975. (37) Grimme, S. J. Chem. Phys. 2003, 118, 9095−9102. (38) Møller, C.; Plesset, M. S. Phys. Rev. 1934, 46, 618−622. (39) Cremer, D. In Encyclopedia of Computational Chemistry; Schleyer, P. v. R., Ed.; John Wiley and Sons: New York, 1998. (40) Parr, R. G.; Yang, W. Density-Functional Theory of Atoms and Molecules; Oxford University Press: Oxford, 1989. (41) Koch, W.; Holthausen, M. C. A Chemist’s Guide to Density Functional Theory; Wiley-VCH: Weinheim, 2000. (42) Dewar, M. J. S.; Zoebisch, E. G.; Healy, E. F.; Stewart, J. J. P. J. Am. Chem. Soc. 1985, 107, 3902−3909. (43) Raczynska, E. D.; Krygowski, T. M.; Duczmal, K.; Hallmann, M. Book of Abstracts, XVIII International Conference on Physical Organic Chemistry, Warsaw, 2006; p 31. (44) Vosko, S. H.; Wilk, L.; Nusair, M. Can. J. Phys. 1980, 58, 1200− 1211. (45) Lee, C.; Yang, W.; Parr, R. G. Phys. Rev. B 1988, 37, 785−789. (46) Becke, A. D. J. Chem. Phys. 1993, 98, 5648−5652. (47) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. J. Phys. Chem. 1994, 98, 11623−11627. (48) Becke, A. D. J. Chem. Phys. 1993, 98, 1372−1377. (49) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; et al. Gaussian 09, Revision A.1; Gaussian, Inc.: Wallingford, CT, 2009. (50) Zhao, Y.; Truhlar, D. G. J. Phys. Chem. A 2005, 109, 5656−5667. (51) Jaguar, version 7.7; Schrödinger, LLC, New York, NY, 2010. (52) Hehre, W. J.; Ditchfield, R.; Pople, J. A. J. Chem. Phys. 1972, 56, 2257−2261. (53) Hariharan, P. C.; Pople, J. A. Theor. Chim. Acta 1973, 28, 213− 222. (54) Francl, M. M.; Pietro, W. J.; Hehre, W. J.; Binkley, J. S.; Gordon, M. S.; DeFrees, D. J.; Pople, J. A. J. Chem. Phys. 1982, 77, 3654−3665. (55) Dunning, T. H. J. Chem. Phys. 1989, 90, 1007−1023. (56) Schleyer, P. v. R.; Maerker, C.; Dransfeld, A.; Jiao, H.; Hommes, N. J. R. v. E. J. Am. Chem. Soc. 1996, 118, 6317−6318. (57) Steiner, E.; Fowler, P. W.; Jennesskens, L. W. Angew. Chem., Int. Ed. 2001, 40, 362−366. (58) Corminboeuf, C.; Heine, T.; Seifert, G.; Schleyer, P. v. R. Phys. Chem. Chem. Phys. 2004, 6, 273−276.

citation including all authors for ref 49. This material is available free of charge via the Internet at http://pubs.acs.org/.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We are grateful to Dr. Boris Proppe from the computer center of the Freie Universität Berlin (ZEDAT) for support. We thank Drs. Arthur Galstyan and Arturo Robertazzi for useful discussions. We like to acknowledge ZEDAT, for providing generously CPU time to perform computations with Gaussian 09.

■

ABBREVIATIONS AM1, Austin model 1; BH&HLYP(G), modified BH&HLYP functional as implemented in Gaussian; d-ζ, double-ζ; t-ζ, tripleζ; q-ζ, quadruple-ζ; DFT, density functional theory; HF, Hartree−Fock theory; HOMED, harmonic oscillator model of electron delocalization; MP2, second-order Møller−Plesset perturbation theory; NICS tensor, nucleus-independent chemical shift tensor; QCISD, quadratic configuration interaction with single and double excitations; QCISD(T), quadratic configuration interaction with single and double excitations and triple excitations added as perturbation; SCS-MP2, spin-component scaled second-order Møller−Plesset perturbation theory; B3LYP-LOC, empirically corrected B3LYP using the localized orbital correction scheme

■

REFERENCES

(1) Watson, J. D.; Crick, F. H. C. Nature 1953, 171, 737−738. (2) Löwdin, P. O. Adv. Quantum Chem. 1965, 2, 213−360. (3) Pullman, B.; Pullman, A. Adv. Heterocycl. Chem. 1971, 13, 77−159. (4) Topal, M. D.; Fresco, J. R. Nature 1976, 263, 285−289. (5) Gerega, A.; Lapinski, L.; Nowak, M. J.; Furmanchuk, A.; Leszczynski, J. J. Phys. Chem. A 2007, 111, 4934−4943. (6) Geballe, M. T.; Skillman, A. G.; Nicholls, A.; Guthrie, J. P.; Taylor, P. J. J. Comput. Aided Mol. Des. 2010, 24, 259−279. (7) Fabian, W. M. F. J. Comput. Chem. 1991, 12, 17−35. (8) Ha, T.-K.; Gunthard, H. H. J. Mol. Struct. (THEOCHEM) 1992, 276, 209−249. (9) Boughton, J. W.; Pulay, P. Int. J. Quantum Chem. 1993, 47, 49−58. (10) Ha, T.-K.; Gunthard, H. H. J. Am. Chem. Soc. 1993, 115, 11939− 11950. (11) Gould, I. R.; Burton, N. A.; Hall, R. J.; Hillier, I. J. Mol. Struct. (Theochem) 1995, 331, 147−154. (12) Paglieri, L.; Corongiu, G.; Estrin, D. A. Int. J. Quantum Chem. 1995, 56, 615−625. (13) Ha, T.-K.; Keller, H. J.; Gunde, R.; Gunthard, H. H. J. Mol. Struct. 1996, 376, 375−397. (14) Gorb, L.; Leszczynski, J. J. Am. Chem. Soc. 1998, 120, 5024−5032. (15) Leszczynski, J. J. Phys. Chem. A 1998, 102, 2357−2362. (16) Dkhissi, A.; Houben, L.; Smets, J.; Adamowicz, L.; Maes, G. J. Mol. Struct. 1999, 484, 215−227. (17) Ha, T.-K.; Keller, H. J.; Gunde, R.; Gunthard, H. H. J. Phys. Chem. A 1999, 103, 6612−6623. (18) Tian, S. X.; Zhang, C. F.; Zhang, Z. J.; Chen, X. J.; Xu, K. Z. Chem. Phys. 1999, 242, 217−223. (19) Kryachko, E. S.; Nguyen, M. T.; Zeegers-Huyskens, T. J. Phys. Chem. A 2001, 105, 1288−1295. (20) Kryachko, E. S.; Nguyen, M. T.; Zeegers-Huyskens, T. J. Phys. Chem. A 2001, 105, 1934−1943. 6892



Article

(59) Fallah-Bagher-Shaidaei, H.; Wannere, C. S.; Corminboeuf, C.; Puchta, R.; Schleyer, P. v. R. Org. Lett. 2006, 8, 863−866. (60) Wolinski, K.; Hinton, J. F.; Pulay, P. J. Am. Chem. Soc. 1990, 112, 8251−8260. (61) Perdew, J. P.; Chevary, J. A.; Vosko, S. H.; Jackson, K. A.; Pederson, M. R.; Singh, D. J.; Fiolhais, C. Phys. Rev. B 1992, 46, 6671− 6687. (62) Perdew, J. P.; Chevary, J. A.; Vosko, S. H.; Jackson, K. A.; Pederson, M. R.; Singh, D. J.; Fiolhais, C. Phys. Rev. B 1993, 48, 4978 (E). (63) Friesner, R. A.; Knoll, E. H.; Cao, Y. J. Chem. Phys. 2006, 125, 124107. (64) Goldfeld, D. A.; Bochevarov, A. D.; Friesner, R. A. J. Chem. Phys. 2008, 129, 214105. (65) Busch, M. S. a.; Knapp, E. W. Chem. Phys. Chem 2004, 5, 1513− 1522.

6893


Isomer Equilibria of Different Lactams

Recommend Documents