Article pubs.acs.org/JPCA
Phenol−Quinone Tautomerism in (Arylazo)naphthols and the Analogous Schiff Bases: Benchmark Calculations S. Tahir Ali,†,§ Liudmil Antonov,‡ and Walter M. F. Fabian*,§ †
Department of Chemistry, Federal Urdu University of Arts, Science & Technology, Karachi, Sindh, Pakistan Institute of Organic Chemistry with Centre of Phytochemistry, Bulgarian Academy of Sciences, Sofia 1113, Bulgaria § Institut für Chemie, Karl-Franzens-Universität Graz, Heinrichstrasse 28, A-8010 Graz, Austria ‡
S Supporting Information *
ABSTRACT: Tautomerization energies of a series of isomeric [(4-R-phenyl)azo]naphthols and the analogous Schiff bases (R = N(CH3)2, OCH3, H, CN, NO2) are calculated by LPNO-CEPA/1-CBS using the def2-TZVPP and def2-QZVPP basis sets for extrapolation. The performance of various density functionals (B3LYP, M06-2X, PW6B95, B2PLYP, mPW2PLYP, PWPB95) as well as MP2 and SCS-MP2 is evaluated against these results. M06-2X and SCS-MP2 yield results close to the LPNO-CEPA/1CBS values. Solvent effects (CCl4, CHCl3, CH3CN, and CH3OH) are treated by a variety of bulk solvation models (SM8, IEFPCM, COSMO, PBF, and SMD) as well as explicit solvation (Monte Carlo free energy perturbation using the OPLSAA force field).
■
INTRODUCTION Reliable predictions of tautomerization energies still present a substantial problem for computational chemistry as, e.g., exemplified by the SAMPL2 challenge.1 Although certain types of tautomeric equilibria could be described quite successfully, for others, significant disagreement with experiment had been observed.2−7 These discrepancies were attributed either to an inadequate treatment of electron correlation or, alternatively, to the solvent model used to describe tautomerization in aqueous solution. Besides the importance of tautomerism in medicinal chemistry, drug design, and chemical information systems,7−11 tautomeric equilibria, e.g., those of (arylazo)naphthols or the analogous Schiff bases, can also be exploited for molecular switches and sensing.12−20 For a rational design of such advanced materials a reliable prediction of the position of the tautomeric equilibrium, including the effect of external stimuli, e.g., change of solvent polarity or complexation of guests, would be of great help.17,21−23 Previously, we had calculated the tautomerization energies of a series of substituted isomeric [(4-R-phenyl)azo]naphthols 1−3 and analogous Schiff bases 4 and 5 (R = N(CH3)2, OCH3, H, CN, and NO2) by B3LYP.22 Though trends of experimental tautomerization constants could be reasonably well described, in a quantitative sense the results were quite poor. Meanwhile, it has become clear that B3LYP generally overestimates the stability of keto forms in DNA-base tautomers including the 2-hydroxypyridine−pyrid-2(1H)-one lactam−lactim tautomerism whereas, in contrast, the stability of lactim forms is overestimated by MP2 calculations.24−28 Hence, the analogous results found by us for the phenol−quinone © 2014 American Chemical Society
(azo−hydrazo, hydroxyiminomethyl−enaminone, enol−keto) tautomerism A−H in (arylazo)naphthols and the corresponding Schiff bases (Scheme 1) is in line with these general observations. In both types of compounds an aromatic phenol A is converted to a quinoid structure H (Scheme 1). Consequently, in the following the notation phenol−quinone tautomerism will be used throughout to denote these equilibria. Recently, a quite systematic study with respect to density functional and basis set of the azo (A)−hydrazo (H) equilibrium 1cA−1cH of compound 1c has been published.29,30 The present study describes an extension to the whole series of substituted derivatives 1a−e−5a−e (Scheme 2). The assessment of the reliability of calculated tautomerization energies is frequently hampered by the lack of reliable experimental data. Moreover, experiments usually quote tautomerization constants KT (ΔG) in some solvent. Differences between calculated and experimental values thus might not only be caused by the electronic structure method but also by the solvent model used or the rigid-rotor harmonic-oscillator approximation for calculating entropic contributions. Moreover, frequently different experimental methods (UV/vis vs NMR) result in slightly up to substantially differing KT values. For validation of approximate computational procedures it has become common practice to use benchmark calculations.31−33 Hence, the primary goal of the present study is to provide benchmark quality tautomerization energies for the phenol− Received: November 22, 2013 Revised: January 8, 2014 Published: January 13, 2014 778
dx.doi.org/10.1021/jp411502u | J. Phys. Chem. A 2014, 118, 778−789
The Journal of Physical Chemistry A
Article
Scheme 1. Azo−Hydrazo and Iminoenol−Enaminone Tautomerism in (Arylazo)naphthols and the Analogous Schiff Bases
Scheme 2. Structures of the Investigated Compounds
■
COMPUTATIONAL DETAILS All geometries were optimized in the gas phase by B3LYP35,36 using the 6-31+G(d,p) basis set37−39 and characterized as true minima by frequency calculations. These geometries were then used for single point calculations with a variety of wave function (MP2,40 SCS-MP2,41 LPNO-CEPA/142−45) and DFT (M062X, 33,46,47 PW6B95, 48 B2PLYP, 49 mPW2PLYP, 50 and PWPB9551). In SCS-MP2, the parallel- and antiparallel-spin pair components to the correlation energy are scaled separately, Ecorr,SCS‑MP2 = pSES(2) + pTET(2). Originally optimized for a benchmark set of 51 reaction energies, in general this scaling resulted in a significant improvement of thermodynamic data for a large set of chemical reactions compared with the standard MP2 method, approaching QCISD or even QCISD(T)
quinone equilibrium of compounds 1a−e−5a−e. In addition, with these benchmark ΔET values, the performance of various bulk solvation models as well as considering explicit solvent by Monte Carlo free energy perturbation will be considered through comparison with experimental ΔGT values. Third, the influence of using different geometries [B3LYP/6-31+G(d,p), M06-2X/6-31+G(d,p), MP2/cc-pVDZ, IEFPCM-B3LYP/631+G(d,p), CAM-B3LYP/6-31+G(d,p)] and dispersion correction will be discussed. Finally, when experimental and calculated Gibbs free energies of tautomerization are compared, the problem of low-frequency modes in calculating entropies will be addressed by the recently proposed quasi-RRHO approach.34 779
dx.doi.org/10.1021/jp411502u | J. Phys. Chem. A 2014, 118, 778−789
The Journal of Physical Chemistry A
Article
quality.41,52 Most importantly, a significant improvement was also found for the energetics of several tautomeric equilibria.24 The scaling parameters applied were those described previously, pS = 6/5 and pT = 1/3.41 A detailed description of the various coupled electron pair methods has been given by Neese et al.43,45 Basically, one starts with a wave function of the CISD type. Variation with respect to the CI coefficients yields the CISD equations. Although size consistency is “simulated” by replacing the correlation energy Ecorr by an excitation dependent shift Δijab, the variational property of the equations as well as unitary invariance with respect to rotations of the occupied orbitals are kept by using localized occupied MOs.43,45 The shifts Δia and Δijab applied in the various electron pair approximations (CPF, ACPF, CEPA/n, etc.) have been summarized by Wennmohs and Neese.45 Basis sets used were def2-TZVPP53,54 for PW6B95 and the double hybrid functionals; 6-311++G(2d,2p)38,55,56 and aug-cc-pVnZ (n = D, T)57,58 for MP2; 6-31+G(d,p) and aug-cc-pVTZ for M06-2X. Complete basis set extrapolation using def2-TZVPP and def2QZVPP has been done with LPNO-CEPA/1 using the modified42,59,60 scheme of Helgaker et al.61 for Ecorr and an exponential extrapolation for ESCF.62 Def2-TZVPP/def2QZVPP basis sets have been shown to give excellent extrapolations to the SCF and correlation basis set limits.60 Zero point energy and thermal corrections to Gibbs free energies were obtained from B3LYP/6-31+G(d,p) calculations using unscaled frequencies. To address the problem of lowfrequency modes in calculating entropies, the method proposed by Grimme34 where the contribution of these modes (ν < 100 cm−1) is replaced by a corresponding rotational entropy (quasiRRHO approach using the published parameters α = 4 and Bav = 10−44 kgm−1), was also applied. Solvent effects (CCl4, CHCl3, CH3CN, and CH3OH) on tautomeric equilibria were estimated by bulk solvent models: IEFPCM63 (B3LYP/6-31+G(d,p)), COSMO64 (B3LYP/6-31+G(d,p)), PBF65 (M06-2X/aug-ccpVTZ(-f)), SM866 (B3LYP/6-31+G(d,p) and M06-2X/631+G(d,p)), and SMD67 (M06-2X/6-31G(d)) as well as explicit solvation by Monte Carlo free energy perturbation (MC-FEP) calculations in NpT (isobaric-isothermic) ensembles (T = 298 K, p = 1 atm) using the OPLS-AA force field.68−71 According to the proposal of Nagy et al.,72 for these latter simulations geometries and atomic charges (CHELPG73) resulting from IEFPCM-B3LYP/6-31+G(d,p) optimizations in the respective solvent (CCl4, CHCl3, CH3CN, or CH3OH) were used as in our previous publication.21 The solution models consisted of one solute and 257 (CH3OH)−263 (CCl4) solvent molecules. The reaction coordinate parameter λ was varied between 0 (phenol) and 1 (quinone tautomer) in steps of 0.025 with linear interpolation of geometric and potential parameters between the corresponding starting and end values. Both solvent−solvent as well as solute−solvent cutoffs were set to 12 Å. In the equilibration and averaging phases, 3500 and 5000 K configurations, respectively, were considered. For DFT calculations the effect of dispersion was estimated by Grimme’s DFTD3 correction.74−76 To estimate the influence of geometry, optimizations were also done by CAM-B3LYP/631+G(d,p),77 M06-2X/6-31+G(d,p), and MP2/cc-pVDZ. In analogy to a previous study,21 for selected tautomeric pairs, IEFPCM-QCISD(T)/6-31G(d) calculations78 in the respective solvent using geometries optimized in the respective solvent by IEFPCM-B3LYP/6-31+G(d,p) were also performed and combined with solvation energies from MC-FEP simulations. Programs used were ORCA,79,80 Gaussian 03,81 GAMESS and
GAMESSPLUS,82−84 NWChem,85 Jaguar,86 BOSS,87,88 and MOLDEN89 for visualization.
■
RESULTS AND DISCUSSION I. Tautomerization Energies. Tautomerization energies obtained by various computational procedures will be compared with those resulting from CEPA/1-CBS//B3LYP/ 6-31+G(d,p) calculations. The choice of this procedure to provide benchmark data is based on the following considerations: (i) The CEPA/1 method has been found to yield energies somewhere intermediate between those obtained by CCSD and CCSD(T) calculations.42−45 Moreover, the def2TZVPP/def2-QZVPP basis sets, although optimized for HF and DFT, show excellent convergence toward the basis set limit.60 (ii) Except in special cases, e.g., [2,2]paracyclophanes,33 the B3LYP density functional yields not only very good geometries90 but also ZPE corrections.91 A more detailed discussion concerning the effect of the geometry will be presented in section III. Some pertinent structural parameters are provided in Table S1 of the Supporting Information. Comparison with experimental X-ray structures92 gives mean unsigned errors (MUE) for the quoted distances of 0.013 [B3LYP/6-31+G(d,p)], 0.012 [CAM-B3LYP/6-31+G(d,p)], 0.013 [M06-2X/6-31+G(d,p)], and 0.016 (MP2(cc-pVDZ); for bond angles and the two dihedral angles describing the twist of the azo (1−3) or iminomethyl (4, 5) and the phenyl group, MUE = 1.5 (B3LYP), 1.7 (CAM-B3LYP), 1.5 (M06-2X), and 1.8 (MP2). Specifically, for fixed model compounds of (4phenylazo)-1-naphthols [4-(methoxynaphthalen-1-yl)phenyl]diazene and [(4-methylphenyl)hydrazono]-4H-naphthalen-1one], better agreement between calculated and experimental geometries was found for B3LYP than for MP2.29 Moreover, although in conjugated oligomers among a number of various density functionals considered BHHLYP, M06-2X, and CAMB3LYP gave the best geometries or, more precisely, bond length alternations for long chain systems compared with benchmark CCSD(T) results, for smaller systems B3LYP or B2PLYP was found to be preferred.93 Obviously, for the molecules treated here, B3LYP is a good choice. (iii) Probably one of the most thoroughly investigated tautomeric equilibria is the lactam−lactim tautomerism of 2-hydroxypyridine−pyrid2(1H)-one,24−28 with a benchmark QCISD(T)/TZV(2df,2dp) calculated value24 of ΔE = 1.0 kcal mol−1. In close agreement, CEPA/1 calculations43,44 using the def2-TZVP basis set resulted in ΔE = 0.9 kcal mol−1. Our LPNO-CEPA/1 CBS (def2-TZVPP/def2-QZVPP) estimate gives ΔE = 1.1 kcal mol−1. Hence, we consider this procedure, LPNO-CEPA-1 (def2-TZVPP/def2-QZVPP) CBS extrapolation using B3LYP/ 6-31+G(d,p) geometries, to provide tautomerization energies ΔET = E(quinone) − E(phenol) of benchmark quality (Table 1). Though electron-withdrawing substituents (R = CN, NO2) preferentially stabilize the hydrazo form of azonaphthols 1−3, eventually even leading to a greater stability of this latter tautomer (2d,e and 3d,e), in the case of the analogous azomethines 4 and 5 no such substituent effect is observable and irrespective of R a substantial (>2 kcal mol−1) preference for the phenol form is evident. The least preference for the aromatic phenol tautomer is obtained for 2-[(4-R-phenyl)azo]1-naphthols 3. Tautomerization energies obtained by a variety of density functionals (B3LYP, 35,36 M06-2X, 33,46,47 PW6B95, 48 B2PLYP,49 mPW2PLYP,50 and PWPB9551 as well as wave 780
dx.doi.org/10.1021/jp411502u | J. Phys. Chem. A 2014, 118, 778−789
The Journal of Physical Chemistry A
Article
substantially better agreement with experiment and/or benchmark calculations.24 With respect to the CEPA/1 def2-CBS tautomerization energies the largest errors are obtained with B3LYP. As expected from previous experience,22 invariably the stability of the quinone tautomers is overestimated. An analogous result, albeit to a lesser extent, is found for PW6B95 as well as the double hybrid functionals. Among the density functionals tested, M06-2X by far leads to the smallest errors. Hence, if a DFT method is used for the treatment of tautomeric equilibria, M06-2X appears to be a good choice as also exemplified by previous results on substituted isomeric hydroxy naphthaldehyde Schiff bases 4 and 5.13 On the contrary, also in agreement with previous experience on these type of compounds,21 MP2 calculations invariably overestimate the stability of the aromatic phenol tautomers by almost 2 kcal mol−1. Spin-component scaling leads to a significant improvement with errors barely exceeding 0.5 kcal mol−1. Hence, if MP2 calculations are used for calculating tautomerization energies, scaling of the spin components is strongly recommended. Fortunately, CEPA/1 calculations using the def2-TZVP basis set give excellent results and hence can be used to avoid the rather costly QZVPP calculations. In view of the very good cost/performance ratio, LPNO-CEPA/1-def2-TZVP calculations appear to be the best choice for treating this type of tautomeric equilibria. II. Solvent Effects. Solvents have a profound influence on tautomeric equilibria with polar solvents usually stabilizing the quinone (hydrazo, enaminone) tautomers H in azonaphthols or the analogous Schiff bases.13,22,23,94−100 Table 3 reports the differences in Gibbs free energies of solvation ΔGsolv = Gsolv(quinone) − Gsolv(phenol) obtained by SMD-M06-2X/631G(d) calculations as used by Truhlar et al. in the SAMPL2 tautomerism challenge.5 Results of other bulk solvation models as well as considering explicit solvent by Monte Carlo free
Table 1. Calculated (CEPA-1 CBS) Tautomerization Energies ΔE for Compounds 1a−e−5a−ea
1 2 3 4 5
N(CH3)2
OCH3
H
CN
a
b
c
d
NO2 e
3.5 2.5 1.5 3.3 2.6
3.2 1.7 1.4 3.4 2.7
1.5 0.7 0.0 3.1 2.5
0.8 −1.2 −1.1 3.1 2.4
0.7 −0.7 −1.3 3.3 2.5
a
For extrapolation the def2-TZVPP/def2-QZVPP basis sets and B3LYP/6-31+G(d,p) geometries were used. ΔE = E(quinone) − E(phenol) in kcal mol−1. A negative sign indicates greater stability of the quinone tautomer H.
function based methods (MP2,40 SCS-MP2;41 LPNO-CEPA/1 using the def2-TZVP, def2-TZVPP, and def2-QZVPP basis sets) are provided in Table S2 of the Supporting Information. Also there are the errors ΔET(CEPA/1-CBS) − ΔET(x) with respect to the CEPA/1 CBS results. Mean signed (MSE) and unsigned (MUE) as well as root-mean-square (RMSE) errors are given in Table 2. According to the above definition of the Table 2. Mean Signed (MSE), Mean Unsigned (MUE), and Root-Mean-Square (RMSE) Errors (kcal mol−1) for Various Computational Procedures with Respect to CEPA-1 CBS ΔE Values B3LYP/6-31+G(d,p) M06-2X/6-31+G(d,p) M06-2X/aug-cc-pVTZ PW6B95/def2-TZVPP B2PLYP/def2-TZVPP mPW2PLYP/def2-TZVPP PWPB95/def2-TZVPP MP2/6-311++G(2d,2p) MP2/aug-cc-pVDZ MP2/aug-cc-pVTZ SCS-MP2/6-311++G(2d,2p) SCS-MP2/aug-cc-pVDZ SCS-MP2/aug-cc-pVTZ LPNO-CEPA-1/def2-TZVP LPNO-CEPA-1/def2-TZVPP LPNO-CEPA-1/def2-QZVPP
MSEa
MUE
RMSE
2.0 −0.5 −0.1 1.6 1.4 1.3 1.5 −1.8 −1.7 −1.8 0.3 0.5 0.4 0.0 −0.2 −0.1
2.0 0.6 0.4 1.6 1.4 1.3 1.5 1.8 1.7 1.8 0.3 0.5 0.5 0.2 0.2 0.1
2.1 0.8 0.5 1.7 1.5 1.4 1.6 1.9 1.8 1.9 0.4 0.6 0.5 0.2 0.2 0.1
Table 3. Calculated [SMD-M06-2X/6-31G(d)//B3LYP/631+G(d,p)] Differences in Solvation Energies ΔGsolv (kcal mol−1)a
Error = ΔET(CEPA/1-CBS) − ΔET(x). Hence, a negative sign indicates overestimation of the stability of the aromatic hydroxy tautomer A compared with CEPA/1-CBS results. a
errors, a negative sign indicates a greater stability of the phenol form calculated by the respective method as compared with the CEPA/1-CBS result. We have chosen the M06-2X and PW6B95 density functionals because the former has been found to provide quite reliable tautomerization energies for a variety of systems;7,25 the latter turned out to be one of the more robust functionals, including tautomer pair energies.28,31 In a previous publication concerning the 1cA−1cH tautomeric equilibrium, B2PLYP was found to give good results,30 and hence this double hybrid functional and variants were included in the present investigation. Although MP2 was found to significantly overestimate the stability of aromatic tautomers,21,24,27 the spin-component scaled version led to
R
solvent
1
2
3
4
5
a
N(CH3)2
b
OCH3
c
H
d
CN
e
NO2
CH3OH CH3CN CHCl3 CCl4 CH3OH CH3CN CHCl3 CCl4 CH3OH CH3CN CHCl3 CCl4 CH3OH CH3CN CHCl3 CCl4 CH3OH CH3CN CHCl3 CCl4
−0.8 0.0 −0.2 −0.4 −1.2 −0.5 −0.6 −0.6 −1.1 −0.4 −0.5 −0.6 −1.9 −1.3 −1.2 −1.0 −2.1 −1.3 −1.3 −1.0
−0.4 0.0 −0.1 −0.1 −1.0 −0.7 −0.6 −0.4 −0.7 −0.5 −0.4 −0.3 −1.1 −1.0 −0.8 −0.6 −1.2 −1.0 −0.7 −0.6
−0.4 −0.2 −0.2 −0.2 −0.5 −0.4 −0.3 −0.3 −0.5 −0.4 −0.3 −0.3 −0.7 −0.7 −0.5 −0.4 −0.8 −0.8 −0.6 −0.5
−2.3 −1.5 −1.3 −0.8 −2.1 −1.4 −1.2 −0.8 −2.0 −1.4 −1.1 −0.7 −1.5 −1.1 −0.9 −0.6 −1.4 −1.1 −0.9 −0.5
−1.6 −1.2 −1.0 −0.6 −1.9 −1.5 −1.2 −0.8 −1.4 −1.1 −0.9 −0.6 −1.8 −1.5 −1.2 −0.8 −1.4 −1.1 −0.9 −0.6
ΔGsolv = Gsolv(quinone) − Gsolv(phenol). Hence, a negative sign indicates greater stabilization of the quinone tautomers H.
a
781
dx.doi.org/10.1021/jp411502u | J. Phys. Chem. A 2014, 118, 778−789
The Journal of Physical Chemistry A
Article
is very good. Compared with M06-2X geometries, B3LYP bond lengths usually are longer but quite similar (mean signed error MSE = −0.004 Å; mean unsigned error MUE = 0.007 Å); deviations are in the range −0.019 to +0.009 Å. In contrast, MP2/cc-pVDZ bond lengths are usually longer than the corresponding B3LYP values (mean signed error MSE = 0.006 Å; mean unsigned error MUE = 0.006 Å), and deviations are in the range −0.007 to +0.022 Å. Deviations of M06-2X from B3LYP bond angles are in the range −2.3° to +1.5° (MSE = 0.0, MUE = 0.3); those for MP2 are in the range −3.1° to +1.2° (MSE = −0.3, MUE = 0.8). Larger differences are found for the (substituted) phenyl torsional angle N(CH)−N−Cipso− Cortho.29,30 In the case of azo derivatives, B3LYP/6-31+G(d,p) calculations result in completely planar structures irrespective of the substituent or the tautomeric form (Table S1, Supporting Information). Similar results are obtained with M06-2X/631+G(d,p) for both tautomers of 2a−2e and 3a−3e as well as 1a and 1b, i.e., donor-substituted [(4-R-phenyl)azo]-1-naphthols, whereas for 1c−1e, i.e., unsubstituted and acceptorsubstituted derivatives, only the hydrazo tautomer adopts a planar aryl conformation. MP2/cc-pVDZ calculations result in planar structures for both tautomers of 3a−3e as well as the hydrazo tautomers of 1a−1e and 2a−2e. For azo tautomers 1aA−1eA slightly twisted structures are obtained irrespective of the substituent as well as for the unsubstituted and acceptor derivatives 2cA−2eA; 2aA and 2bA are predicted to be planar (Table S1, Supporting Information). Considerably larger aryl twists τ(CH−N−Cipso−Cortho) are calculated for the corresponding Schiff bases 4a−4e and 5a−5e, especially the phenol tautomers. Transformation to the quinone tautomers significantly reduces these dihedral angles, especially in the case of R = acceptor (4dH, 4eH, 5dH, and 5eH), eventually leading to planar B3LYP and M06-2X structures for these derivatives. With MP2/cc-pVDZ these derivatives still are characterized by an aryl twist but the trend of increasing planarization in acceptor-substituted derivatives is also clearly evident (Table S1, Supporting Information). Inclusion of long-range correlation as in CAM-B3LYP has little influence on calculated geometries compared with those obtained with B3LYP. Mean unsigned errors are 0.006 Å for bond lengths and 0.2° for bond angles. Concerning the aryl twists, B3LYP and CAM-B3LYP optimizations lead to completely analogous results. In the case of Schiff bases 4 and 5 with the latter density functional somewhat larger deviations from planarity are calculated (Table S1, Supporting Information). Solvation not only will influence relative energies but also the structures of individual tautomers.101 Consequently, all molecules were also optimized in the respective solvents CH3OH, CH3CN, CHCl3, and CCl4 using a continuum solvation model [IEFPCM-B3LYP/631+G(d,p)]. As examples, selected geometrical parameters obtained for CH3OH as the solvent are collected in Table S1, Supporting Information (for 5 the corresponding data for CH3CN are also given). The in-solution geometries only marginally differ from gas phase structures, usually
Phenol−Quinone Tautomerism in (Arylazo)naphthols and the Analogous Schiff Bases: Benchmark Calculations S. Tahir Ali,†,§ Liudmil Antonov,‡ and Walter M. F. Fabian*,§ †
Department of Chemistry, Federal Urdu University of Arts, Science & Technology, Karachi, Sindh, Pakistan Institute of Organic Chemistry with Centre of Phytochemistry, Bulgarian Academy of Sciences, Sofia 1113, Bulgaria § Institut für Chemie, Karl-Franzens-Universität Graz, Heinrichstrasse 28, A-8010 Graz, Austria ‡
S Supporting Information *
ABSTRACT: Tautomerization energies of a series of isomeric [(4-R-phenyl)azo]naphthols and the analogous Schiff bases (R = N(CH3)2, OCH3, H, CN, NO2) are calculated by LPNO-CEPA/1-CBS using the def2-TZVPP and def2-QZVPP basis sets for extrapolation. The performance of various density functionals (B3LYP, M06-2X, PW6B95, B2PLYP, mPW2PLYP, PWPB95) as well as MP2 and SCS-MP2 is evaluated against these results. M06-2X and SCS-MP2 yield results close to the LPNO-CEPA/1CBS values. Solvent effects (CCl4, CHCl3, CH3CN, and CH3OH) are treated by a variety of bulk solvation models (SM8, IEFPCM, COSMO, PBF, and SMD) as well as explicit solvation (Monte Carlo free energy perturbation using the OPLSAA force field).
■
INTRODUCTION Reliable predictions of tautomerization energies still present a substantial problem for computational chemistry as, e.g., exemplified by the SAMPL2 challenge.1 Although certain types of tautomeric equilibria could be described quite successfully, for others, significant disagreement with experiment had been observed.2−7 These discrepancies were attributed either to an inadequate treatment of electron correlation or, alternatively, to the solvent model used to describe tautomerization in aqueous solution. Besides the importance of tautomerism in medicinal chemistry, drug design, and chemical information systems,7−11 tautomeric equilibria, e.g., those of (arylazo)naphthols or the analogous Schiff bases, can also be exploited for molecular switches and sensing.12−20 For a rational design of such advanced materials a reliable prediction of the position of the tautomeric equilibrium, including the effect of external stimuli, e.g., change of solvent polarity or complexation of guests, would be of great help.17,21−23 Previously, we had calculated the tautomerization energies of a series of substituted isomeric [(4-R-phenyl)azo]naphthols 1−3 and analogous Schiff bases 4 and 5 (R = N(CH3)2, OCH3, H, CN, and NO2) by B3LYP.22 Though trends of experimental tautomerization constants could be reasonably well described, in a quantitative sense the results were quite poor. Meanwhile, it has become clear that B3LYP generally overestimates the stability of keto forms in DNA-base tautomers including the 2-hydroxypyridine−pyrid-2(1H)-one lactam−lactim tautomerism whereas, in contrast, the stability of lactim forms is overestimated by MP2 calculations.24−28 Hence, the analogous results found by us for the phenol−quinone © 2014 American Chemical Society
(azo−hydrazo, hydroxyiminomethyl−enaminone, enol−keto) tautomerism A−H in (arylazo)naphthols and the corresponding Schiff bases (Scheme 1) is in line with these general observations. In both types of compounds an aromatic phenol A is converted to a quinoid structure H (Scheme 1). Consequently, in the following the notation phenol−quinone tautomerism will be used throughout to denote these equilibria. Recently, a quite systematic study with respect to density functional and basis set of the azo (A)−hydrazo (H) equilibrium 1cA−1cH of compound 1c has been published.29,30 The present study describes an extension to the whole series of substituted derivatives 1a−e−5a−e (Scheme 2). The assessment of the reliability of calculated tautomerization energies is frequently hampered by the lack of reliable experimental data. Moreover, experiments usually quote tautomerization constants KT (ΔG) in some solvent. Differences between calculated and experimental values thus might not only be caused by the electronic structure method but also by the solvent model used or the rigid-rotor harmonic-oscillator approximation for calculating entropic contributions. Moreover, frequently different experimental methods (UV/vis vs NMR) result in slightly up to substantially differing KT values. For validation of approximate computational procedures it has become common practice to use benchmark calculations.31−33 Hence, the primary goal of the present study is to provide benchmark quality tautomerization energies for the phenol− Received: November 22, 2013 Revised: January 8, 2014 Published: January 13, 2014 778
dx.doi.org/10.1021/jp411502u | J. Phys. Chem. A 2014, 118, 778−789
The Journal of Physical Chemistry A
Article
Scheme 1. Azo−Hydrazo and Iminoenol−Enaminone Tautomerism in (Arylazo)naphthols and the Analogous Schiff Bases
Scheme 2. Structures of the Investigated Compounds
■
COMPUTATIONAL DETAILS All geometries were optimized in the gas phase by B3LYP35,36 using the 6-31+G(d,p) basis set37−39 and characterized as true minima by frequency calculations. These geometries were then used for single point calculations with a variety of wave function (MP2,40 SCS-MP2,41 LPNO-CEPA/142−45) and DFT (M062X, 33,46,47 PW6B95, 48 B2PLYP, 49 mPW2PLYP, 50 and PWPB9551). In SCS-MP2, the parallel- and antiparallel-spin pair components to the correlation energy are scaled separately, Ecorr,SCS‑MP2 = pSES(2) + pTET(2). Originally optimized for a benchmark set of 51 reaction energies, in general this scaling resulted in a significant improvement of thermodynamic data for a large set of chemical reactions compared with the standard MP2 method, approaching QCISD or even QCISD(T)
quinone equilibrium of compounds 1a−e−5a−e. In addition, with these benchmark ΔET values, the performance of various bulk solvation models as well as considering explicit solvent by Monte Carlo free energy perturbation will be considered through comparison with experimental ΔGT values. Third, the influence of using different geometries [B3LYP/6-31+G(d,p), M06-2X/6-31+G(d,p), MP2/cc-pVDZ, IEFPCM-B3LYP/631+G(d,p), CAM-B3LYP/6-31+G(d,p)] and dispersion correction will be discussed. Finally, when experimental and calculated Gibbs free energies of tautomerization are compared, the problem of low-frequency modes in calculating entropies will be addressed by the recently proposed quasi-RRHO approach.34 779
dx.doi.org/10.1021/jp411502u | J. Phys. Chem. A 2014, 118, 778−789
The Journal of Physical Chemistry A
Article
quality.41,52 Most importantly, a significant improvement was also found for the energetics of several tautomeric equilibria.24 The scaling parameters applied were those described previously, pS = 6/5 and pT = 1/3.41 A detailed description of the various coupled electron pair methods has been given by Neese et al.43,45 Basically, one starts with a wave function of the CISD type. Variation with respect to the CI coefficients yields the CISD equations. Although size consistency is “simulated” by replacing the correlation energy Ecorr by an excitation dependent shift Δijab, the variational property of the equations as well as unitary invariance with respect to rotations of the occupied orbitals are kept by using localized occupied MOs.43,45 The shifts Δia and Δijab applied in the various electron pair approximations (CPF, ACPF, CEPA/n, etc.) have been summarized by Wennmohs and Neese.45 Basis sets used were def2-TZVPP53,54 for PW6B95 and the double hybrid functionals; 6-311++G(2d,2p)38,55,56 and aug-cc-pVnZ (n = D, T)57,58 for MP2; 6-31+G(d,p) and aug-cc-pVTZ for M06-2X. Complete basis set extrapolation using def2-TZVPP and def2QZVPP has been done with LPNO-CEPA/1 using the modified42,59,60 scheme of Helgaker et al.61 for Ecorr and an exponential extrapolation for ESCF.62 Def2-TZVPP/def2QZVPP basis sets have been shown to give excellent extrapolations to the SCF and correlation basis set limits.60 Zero point energy and thermal corrections to Gibbs free energies were obtained from B3LYP/6-31+G(d,p) calculations using unscaled frequencies. To address the problem of lowfrequency modes in calculating entropies, the method proposed by Grimme34 where the contribution of these modes (ν < 100 cm−1) is replaced by a corresponding rotational entropy (quasiRRHO approach using the published parameters α = 4 and Bav = 10−44 kgm−1), was also applied. Solvent effects (CCl4, CHCl3, CH3CN, and CH3OH) on tautomeric equilibria were estimated by bulk solvent models: IEFPCM63 (B3LYP/6-31+G(d,p)), COSMO64 (B3LYP/6-31+G(d,p)), PBF65 (M06-2X/aug-ccpVTZ(-f)), SM866 (B3LYP/6-31+G(d,p) and M06-2X/631+G(d,p)), and SMD67 (M06-2X/6-31G(d)) as well as explicit solvation by Monte Carlo free energy perturbation (MC-FEP) calculations in NpT (isobaric-isothermic) ensembles (T = 298 K, p = 1 atm) using the OPLS-AA force field.68−71 According to the proposal of Nagy et al.,72 for these latter simulations geometries and atomic charges (CHELPG73) resulting from IEFPCM-B3LYP/6-31+G(d,p) optimizations in the respective solvent (CCl4, CHCl3, CH3CN, or CH3OH) were used as in our previous publication.21 The solution models consisted of one solute and 257 (CH3OH)−263 (CCl4) solvent molecules. The reaction coordinate parameter λ was varied between 0 (phenol) and 1 (quinone tautomer) in steps of 0.025 with linear interpolation of geometric and potential parameters between the corresponding starting and end values. Both solvent−solvent as well as solute−solvent cutoffs were set to 12 Å. In the equilibration and averaging phases, 3500 and 5000 K configurations, respectively, were considered. For DFT calculations the effect of dispersion was estimated by Grimme’s DFTD3 correction.74−76 To estimate the influence of geometry, optimizations were also done by CAM-B3LYP/631+G(d,p),77 M06-2X/6-31+G(d,p), and MP2/cc-pVDZ. In analogy to a previous study,21 for selected tautomeric pairs, IEFPCM-QCISD(T)/6-31G(d) calculations78 in the respective solvent using geometries optimized in the respective solvent by IEFPCM-B3LYP/6-31+G(d,p) were also performed and combined with solvation energies from MC-FEP simulations. Programs used were ORCA,79,80 Gaussian 03,81 GAMESS and
GAMESSPLUS,82−84 NWChem,85 Jaguar,86 BOSS,87,88 and MOLDEN89 for visualization.
■
RESULTS AND DISCUSSION I. Tautomerization Energies. Tautomerization energies obtained by various computational procedures will be compared with those resulting from CEPA/1-CBS//B3LYP/ 6-31+G(d,p) calculations. The choice of this procedure to provide benchmark data is based on the following considerations: (i) The CEPA/1 method has been found to yield energies somewhere intermediate between those obtained by CCSD and CCSD(T) calculations.42−45 Moreover, the def2TZVPP/def2-QZVPP basis sets, although optimized for HF and DFT, show excellent convergence toward the basis set limit.60 (ii) Except in special cases, e.g., [2,2]paracyclophanes,33 the B3LYP density functional yields not only very good geometries90 but also ZPE corrections.91 A more detailed discussion concerning the effect of the geometry will be presented in section III. Some pertinent structural parameters are provided in Table S1 of the Supporting Information. Comparison with experimental X-ray structures92 gives mean unsigned errors (MUE) for the quoted distances of 0.013 [B3LYP/6-31+G(d,p)], 0.012 [CAM-B3LYP/6-31+G(d,p)], 0.013 [M06-2X/6-31+G(d,p)], and 0.016 (MP2(cc-pVDZ); for bond angles and the two dihedral angles describing the twist of the azo (1−3) or iminomethyl (4, 5) and the phenyl group, MUE = 1.5 (B3LYP), 1.7 (CAM-B3LYP), 1.5 (M06-2X), and 1.8 (MP2). Specifically, for fixed model compounds of (4phenylazo)-1-naphthols [4-(methoxynaphthalen-1-yl)phenyl]diazene and [(4-methylphenyl)hydrazono]-4H-naphthalen-1one], better agreement between calculated and experimental geometries was found for B3LYP than for MP2.29 Moreover, although in conjugated oligomers among a number of various density functionals considered BHHLYP, M06-2X, and CAMB3LYP gave the best geometries or, more precisely, bond length alternations for long chain systems compared with benchmark CCSD(T) results, for smaller systems B3LYP or B2PLYP was found to be preferred.93 Obviously, for the molecules treated here, B3LYP is a good choice. (iii) Probably one of the most thoroughly investigated tautomeric equilibria is the lactam−lactim tautomerism of 2-hydroxypyridine−pyrid2(1H)-one,24−28 with a benchmark QCISD(T)/TZV(2df,2dp) calculated value24 of ΔE = 1.0 kcal mol−1. In close agreement, CEPA/1 calculations43,44 using the def2-TZVP basis set resulted in ΔE = 0.9 kcal mol−1. Our LPNO-CEPA/1 CBS (def2-TZVPP/def2-QZVPP) estimate gives ΔE = 1.1 kcal mol−1. Hence, we consider this procedure, LPNO-CEPA-1 (def2-TZVPP/def2-QZVPP) CBS extrapolation using B3LYP/ 6-31+G(d,p) geometries, to provide tautomerization energies ΔET = E(quinone) − E(phenol) of benchmark quality (Table 1). Though electron-withdrawing substituents (R = CN, NO2) preferentially stabilize the hydrazo form of azonaphthols 1−3, eventually even leading to a greater stability of this latter tautomer (2d,e and 3d,e), in the case of the analogous azomethines 4 and 5 no such substituent effect is observable and irrespective of R a substantial (>2 kcal mol−1) preference for the phenol form is evident. The least preference for the aromatic phenol tautomer is obtained for 2-[(4-R-phenyl)azo]1-naphthols 3. Tautomerization energies obtained by a variety of density functionals (B3LYP, 35,36 M06-2X, 33,46,47 PW6B95, 48 B2PLYP,49 mPW2PLYP,50 and PWPB9551 as well as wave 780
dx.doi.org/10.1021/jp411502u | J. Phys. Chem. A 2014, 118, 778−789
The Journal of Physical Chemistry A
Article
substantially better agreement with experiment and/or benchmark calculations.24 With respect to the CEPA/1 def2-CBS tautomerization energies the largest errors are obtained with B3LYP. As expected from previous experience,22 invariably the stability of the quinone tautomers is overestimated. An analogous result, albeit to a lesser extent, is found for PW6B95 as well as the double hybrid functionals. Among the density functionals tested, M06-2X by far leads to the smallest errors. Hence, if a DFT method is used for the treatment of tautomeric equilibria, M06-2X appears to be a good choice as also exemplified by previous results on substituted isomeric hydroxy naphthaldehyde Schiff bases 4 and 5.13 On the contrary, also in agreement with previous experience on these type of compounds,21 MP2 calculations invariably overestimate the stability of the aromatic phenol tautomers by almost 2 kcal mol−1. Spin-component scaling leads to a significant improvement with errors barely exceeding 0.5 kcal mol−1. Hence, if MP2 calculations are used for calculating tautomerization energies, scaling of the spin components is strongly recommended. Fortunately, CEPA/1 calculations using the def2-TZVP basis set give excellent results and hence can be used to avoid the rather costly QZVPP calculations. In view of the very good cost/performance ratio, LPNO-CEPA/1-def2-TZVP calculations appear to be the best choice for treating this type of tautomeric equilibria. II. Solvent Effects. Solvents have a profound influence on tautomeric equilibria with polar solvents usually stabilizing the quinone (hydrazo, enaminone) tautomers H in azonaphthols or the analogous Schiff bases.13,22,23,94−100 Table 3 reports the differences in Gibbs free energies of solvation ΔGsolv = Gsolv(quinone) − Gsolv(phenol) obtained by SMD-M06-2X/631G(d) calculations as used by Truhlar et al. in the SAMPL2 tautomerism challenge.5 Results of other bulk solvation models as well as considering explicit solvent by Monte Carlo free
Table 1. Calculated (CEPA-1 CBS) Tautomerization Energies ΔE for Compounds 1a−e−5a−ea
1 2 3 4 5
N(CH3)2
OCH3
H
CN
a
b
c
d
NO2 e
3.5 2.5 1.5 3.3 2.6
3.2 1.7 1.4 3.4 2.7
1.5 0.7 0.0 3.1 2.5
0.8 −1.2 −1.1 3.1 2.4
0.7 −0.7 −1.3 3.3 2.5
a
For extrapolation the def2-TZVPP/def2-QZVPP basis sets and B3LYP/6-31+G(d,p) geometries were used. ΔE = E(quinone) − E(phenol) in kcal mol−1. A negative sign indicates greater stability of the quinone tautomer H.
function based methods (MP2,40 SCS-MP2;41 LPNO-CEPA/1 using the def2-TZVP, def2-TZVPP, and def2-QZVPP basis sets) are provided in Table S2 of the Supporting Information. Also there are the errors ΔET(CEPA/1-CBS) − ΔET(x) with respect to the CEPA/1 CBS results. Mean signed (MSE) and unsigned (MUE) as well as root-mean-square (RMSE) errors are given in Table 2. According to the above definition of the Table 2. Mean Signed (MSE), Mean Unsigned (MUE), and Root-Mean-Square (RMSE) Errors (kcal mol−1) for Various Computational Procedures with Respect to CEPA-1 CBS ΔE Values B3LYP/6-31+G(d,p) M06-2X/6-31+G(d,p) M06-2X/aug-cc-pVTZ PW6B95/def2-TZVPP B2PLYP/def2-TZVPP mPW2PLYP/def2-TZVPP PWPB95/def2-TZVPP MP2/6-311++G(2d,2p) MP2/aug-cc-pVDZ MP2/aug-cc-pVTZ SCS-MP2/6-311++G(2d,2p) SCS-MP2/aug-cc-pVDZ SCS-MP2/aug-cc-pVTZ LPNO-CEPA-1/def2-TZVP LPNO-CEPA-1/def2-TZVPP LPNO-CEPA-1/def2-QZVPP
MSEa
MUE
RMSE
2.0 −0.5 −0.1 1.6 1.4 1.3 1.5 −1.8 −1.7 −1.8 0.3 0.5 0.4 0.0 −0.2 −0.1
2.0 0.6 0.4 1.6 1.4 1.3 1.5 1.8 1.7 1.8 0.3 0.5 0.5 0.2 0.2 0.1
2.1 0.8 0.5 1.7 1.5 1.4 1.6 1.9 1.8 1.9 0.4 0.6 0.5 0.2 0.2 0.1
Table 3. Calculated [SMD-M06-2X/6-31G(d)//B3LYP/631+G(d,p)] Differences in Solvation Energies ΔGsolv (kcal mol−1)a
Error = ΔET(CEPA/1-CBS) − ΔET(x). Hence, a negative sign indicates overestimation of the stability of the aromatic hydroxy tautomer A compared with CEPA/1-CBS results. a
errors, a negative sign indicates a greater stability of the phenol form calculated by the respective method as compared with the CEPA/1-CBS result. We have chosen the M06-2X and PW6B95 density functionals because the former has been found to provide quite reliable tautomerization energies for a variety of systems;7,25 the latter turned out to be one of the more robust functionals, including tautomer pair energies.28,31 In a previous publication concerning the 1cA−1cH tautomeric equilibrium, B2PLYP was found to give good results,30 and hence this double hybrid functional and variants were included in the present investigation. Although MP2 was found to significantly overestimate the stability of aromatic tautomers,21,24,27 the spin-component scaled version led to
R
solvent
1
2
3
4
5
a
N(CH3)2
b
OCH3
c
H
d
CN
e
NO2
CH3OH CH3CN CHCl3 CCl4 CH3OH CH3CN CHCl3 CCl4 CH3OH CH3CN CHCl3 CCl4 CH3OH CH3CN CHCl3 CCl4 CH3OH CH3CN CHCl3 CCl4
−0.8 0.0 −0.2 −0.4 −1.2 −0.5 −0.6 −0.6 −1.1 −0.4 −0.5 −0.6 −1.9 −1.3 −1.2 −1.0 −2.1 −1.3 −1.3 −1.0
−0.4 0.0 −0.1 −0.1 −1.0 −0.7 −0.6 −0.4 −0.7 −0.5 −0.4 −0.3 −1.1 −1.0 −0.8 −0.6 −1.2 −1.0 −0.7 −0.6
−0.4 −0.2 −0.2 −0.2 −0.5 −0.4 −0.3 −0.3 −0.5 −0.4 −0.3 −0.3 −0.7 −0.7 −0.5 −0.4 −0.8 −0.8 −0.6 −0.5
−2.3 −1.5 −1.3 −0.8 −2.1 −1.4 −1.2 −0.8 −2.0 −1.4 −1.1 −0.7 −1.5 −1.1 −0.9 −0.6 −1.4 −1.1 −0.9 −0.5
−1.6 −1.2 −1.0 −0.6 −1.9 −1.5 −1.2 −0.8 −1.4 −1.1 −0.9 −0.6 −1.8 −1.5 −1.2 −0.8 −1.4 −1.1 −0.9 −0.6
ΔGsolv = Gsolv(quinone) − Gsolv(phenol). Hence, a negative sign indicates greater stabilization of the quinone tautomers H.
a
781
dx.doi.org/10.1021/jp411502u | J. Phys. Chem. A 2014, 118, 778−789
The Journal of Physical Chemistry A
Article
is very good. Compared with M06-2X geometries, B3LYP bond lengths usually are longer but quite similar (mean signed error MSE = −0.004 Å; mean unsigned error MUE = 0.007 Å); deviations are in the range −0.019 to +0.009 Å. In contrast, MP2/cc-pVDZ bond lengths are usually longer than the corresponding B3LYP values (mean signed error MSE = 0.006 Å; mean unsigned error MUE = 0.006 Å), and deviations are in the range −0.007 to +0.022 Å. Deviations of M06-2X from B3LYP bond angles are in the range −2.3° to +1.5° (MSE = 0.0, MUE = 0.3); those for MP2 are in the range −3.1° to +1.2° (MSE = −0.3, MUE = 0.8). Larger differences are found for the (substituted) phenyl torsional angle N(CH)−N−Cipso− Cortho.29,30 In the case of azo derivatives, B3LYP/6-31+G(d,p) calculations result in completely planar structures irrespective of the substituent or the tautomeric form (Table S1, Supporting Information). Similar results are obtained with M06-2X/631+G(d,p) for both tautomers of 2a−2e and 3a−3e as well as 1a and 1b, i.e., donor-substituted [(4-R-phenyl)azo]-1-naphthols, whereas for 1c−1e, i.e., unsubstituted and acceptorsubstituted derivatives, only the hydrazo tautomer adopts a planar aryl conformation. MP2/cc-pVDZ calculations result in planar structures for both tautomers of 3a−3e as well as the hydrazo tautomers of 1a−1e and 2a−2e. For azo tautomers 1aA−1eA slightly twisted structures are obtained irrespective of the substituent as well as for the unsubstituted and acceptor derivatives 2cA−2eA; 2aA and 2bA are predicted to be planar (Table S1, Supporting Information). Considerably larger aryl twists τ(CH−N−Cipso−Cortho) are calculated for the corresponding Schiff bases 4a−4e and 5a−5e, especially the phenol tautomers. Transformation to the quinone tautomers significantly reduces these dihedral angles, especially in the case of R = acceptor (4dH, 4eH, 5dH, and 5eH), eventually leading to planar B3LYP and M06-2X structures for these derivatives. With MP2/cc-pVDZ these derivatives still are characterized by an aryl twist but the trend of increasing planarization in acceptor-substituted derivatives is also clearly evident (Table S1, Supporting Information). Inclusion of long-range correlation as in CAM-B3LYP has little influence on calculated geometries compared with those obtained with B3LYP. Mean unsigned errors are 0.006 Å for bond lengths and 0.2° for bond angles. Concerning the aryl twists, B3LYP and CAM-B3LYP optimizations lead to completely analogous results. In the case of Schiff bases 4 and 5 with the latter density functional somewhat larger deviations from planarity are calculated (Table S1, Supporting Information). Solvation not only will influence relative energies but also the structures of individual tautomers.101 Consequently, all molecules were also optimized in the respective solvents CH3OH, CH3CN, CHCl3, and CCl4 using a continuum solvation model [IEFPCM-B3LYP/631+G(d,p)]. As examples, selected geometrical parameters obtained for CH3OH as the solvent are collected in Table S1, Supporting Information (for 5 the corresponding data for CH3CN are also given). The in-solution geometries only marginally differ from gas phase structures, usually
