Ionization Energies and Aqueous Redox Potentials of Organic

Article pubs.acs.org/JCTC

Ionization Energies and Aqueous Redox Potentials of Organic Molecules: Comparison of DFT, Correlated ab Initio Theory and Pair Natural Orbital Approaches Miho Isegawa, Frank Neese,* and Dimitrios A. Pantazis* Max Planck Institute for Chemical Energy Conversion, Stiftrasse 34-38, 45470 Mülheim and der Ruhr, Germany S Supporting Information *

ABSTRACT: The calculation of redox potentials involves large energetic terms arising from gas phase ionization energies, thermodynamic contributions, and solvation energies of the reduced and oxidized species. In this work we study the performance of a wide range of wave function and density functional theory methods for the prediction of ionization energies and aqueous one-electron oxidation potentials of a set of 19 organic molecules. Emphasis is placed on evaluating methods that employ the computationally efficient local pair natural orbital (LPNO) approach, as well as several implementations of coupled cluster theory and explicitly correlated F12 methods. The electronic energies are combined with implicit solvation models for the solvation energies. With the exception of MP2 and its variants, which suffer from enormous errors arising at least partially from the poor Hartree−Fock reference, ionization energies can be systematically predicted with average errors below 0.1 eV for most of the correlated wave function based methods studies here, provided basis set extrapolation is performed. LPNO methods are the most efficient way to achieve this type of accuracy. DFT methods show in general larger errors and suffer from inconsistent behavior. The only exception is the M06-2X functional which is found to be competitive with the best LPNO-based approaches for ionization energies. Importantly, the limiting factor for the calculation of accurate redox potentials is the solvation energy. The errors in the predicted solvation energies by all continuum solvation models tested in this work dominate the final computed reduction potential, resulting in average errors typically in excess of 0.3 V and hence obscuring the gains that arise from choosing a more accurate electronic structure method.

■

and properties studied. From the field of wave function based methods, coupled cluster with single and double excitations with perturbative triples, CCSD(T), has long been considered among the most reliable practical electronic structure methods,2,3 at least for main-group chemistry,4 though the associated computational cost is very high with a formal O(N7) scaling with system size. Moreover, coupled cluster methods are very slowly convergent with respect to the one-electron basis set, necessitating the use of basis set extrapolation schemes5−12 or of explicitly correlated approaches13,14 such as R12 and F12.15−17 Drastic reductions in the computational cost of these correlated methods can however be achieved with techniques such as the recently developed local pair natural orbital (LPNO) approach.18−20 These methods exploit the locality of electron correlation in a compact way which limits the correlation problem to orbital pairs that contribute to the correlation energy above a given threshold. LPNO methods converge toward the canonical coupled cluster results at a much reduced cost and yield about 99.9% of the correlation energy.18 This percentage is often sufficient to ensure that in chemical applications the results of LPNO calculations reach the high

INTRODUCTION The prediction of accurate ionization energies and redox potentials is crucial for many fields of molecular science that deal with reactions involving exchange of electrons. Owing to the importance of these quantities in the study of redox processes, many theoretical studies have addressed the question of how to best obtain reliable results, focusing on diverse systems and applications.1 In the context of a thermodynamic cycle, the accurate description of the redox process depends on a quantitatively correct description of the changes in electronic structure that accompany the redox event, as well as of the interaction of the molecule (solute) with the solvent. Given the high requirements for achieving results of predictive quality, ionization energies, electron affinities, and reduction potentials are thus good target properties to evaluate both electronic structure methods and solvation models. In terms of electronic structure methods, the most widely used and cost-effective approach is any of the many flavors of density functional theory (DFT). This approach is most often adopted because of the favorable scaling, low computational cost, and often sufficient accuracy. However, the latter is by no means a given; a “careful choice” of functional is often crucial due to the inconsistent behavior and occasionally even erratic performance of existing functionals for different types of system © 2016 American Chemical Society

Received: March 7, 2016 Published: April 11, 2016 2272

DOI: 10.1021/acs.jctc.6b00252 J. Chem. Theory Comput. 2016, 12, 2272−2284

Article

Journal of Chemical Theory and Computation

Figure 1. Organic molecules included in the present test set.

redox potentials in aqueous solution. We have two parallel aims in this study. The first one is to evaluate and benchmark a number of modern wave function based methods that have not been systematically evaluated so far for such properties, especially methods that employ the LPNO approach, in comparison with a collection of traditional ab initio approaches and a range of DFT functionals. The second is to identify and quantify the errors in the calculation of ionization energies and redox potentials for these one-electron oxidations, thus separately determining the accuracy expected from various electronic structure methods and from common implicit solvation models.

accuracy expected from coupled cluster methods. In most cases, the results are significantly more accurate than DFT approaches,21,22 thus extending the applicability of coupled cluster methods to molecules of size traditionally associated with DFT.23 Considering the contribution of the solvent, required for calculating redox potentials, implicit solvation models24,25 have been of great utility in incorporating solvation effects without significant increase in computational cost compared with the gas phase calculation. As a consequence, they are widely used in combination with various electronic structure methods as well as with multiscale modeling approaches. Limitations of these models arise from certain assumptions regarding the solvation free energy contributions26 and from their inability to properly account for specific interactions such as hydrogen bonding or charge transfer between solvent and solute. On the other hand, explicit solvation is not always a practical alternative because the calculation of solvation free energies requires extensive sampling and may result in calculations that are simply impractical when coupled to expensive electronic structure methods. The prediction of ionization energies, electron affinities, and redox potentials is such a historically relevant and extensively pursued topic in computational chemistry that it is impossible to provide any concise overview of the literature here. Many reports on such calculations have appeared over the years for a broad range of chemically diverse systems, using various approaches that encompass implicit solvation with several quantum mechanical (QM) methods, hybrid approaches with an explicit molecular mechanics treatment of the solvent (QM/ MM), and classical or DFT molecular dynamics. Valuable summaries of the many approaches and pointers to numerous applications are provided in excellent comprehensive reviews1,27,28 and in recent papers that focus on specific classes of compounds and on various methodological aspects.29−46 In the present study, we have collected a set of organic molecules to test a variety of modern wave function based approaches as well as DFT methods for the calculation of ionization energies and subsequently, with the use of two popular implicit solvation models, for obtaining one-electron

■

METHODOLOGY AND COMPUTATIONAL DETAILS Definition of the Test Set. The molecules that comprise the present test set are shown in Figure 1. The test set includes both aromatic and aliphatic molecules. The selection was designed such that a sufficiently broad range of different classes of compounds and functional groups is included but was also guided by the availability of experimentally determined redox potentials and solvation energies. Relevant experimental data with respect to ionization energies (IE), redox potentials (E0), and Henry’s constants (HC), which are used to derive solvation energies, are listed in Table 1.47−63 We note that the experimental uncertainties for the redox potentials are reported in the range 0.01−0.05 eV. For the errors in the experimentally derived Henry’s constants there are large differences depending on the methods employed; therefore, the choices were made according to the recommendations in ref 64: Henry’s constant is first converted to the dimensionless constant sometimes called “water-air partitioning coefficient” PW/A = RTHC

(1)

where R is the ideal gas constant, and T is temperature. Then the solvation energy is calculated as ° = −2.303RT log PW/A ΔG W/A

(2)

Electronic Structure Methods. All computations were performed using the ORCA program system.65 Geometries of 2273


Article


of the correlation energy using the LPNO-CCSD/cc-pV[Q/ 5]Z level as the best method to obtain reference electronic energies. The combination of canonical CCSD(T) with LPNOCCSD (or LPNO-CEPA/1) has been previously evaluated and was found to be superior to extrapolation methods that use MP2 for the basis set incompleteness correction.75 The tested density functionals involve different formalisms and different percentages of Hartree−Fock exchange. They include GGA (PBE76), meta-GGA (TPSS66 and M06-L77), hybrid-GGA (B3LYP78 and PBE079), hybrid-meta-GGA (TPSSh,80 M06, and M06-2X81,82), and double-hybrid GGA functionals (B2PLYP83 and mPW2LYP84). DFT integration grids were increased to “Grid 6” in ORCA notation. Multiple basis sets were used in this study to obtain single point energies at various levels of theory. For all density functionals the def2-TZVPP basis was used to obtain singlepoint energies. With the exception of the explicitly correlated methods, all wave function based methods employed extrapolation to the complete basis set limit using the family of correlation-consistent basis sets cc-pVnZ, with n = D, T, Q, 5, and 6. The cc-pVnZ-F12 basis sets85 were used for explicitly correlated F12 calculations for both second-order perturbation and for coupled cluster methods. It is noted that s and p functions in cc-pVDZ-F12 are the same as those in aug-ccpVTZ with the exception that the most diffuse p function is decontracted. Therefore, the F12 version of this basis set is considerably larger than the conventional valence-double-ζ basis. Since the complementary auxiliary basis set (CABS) approach86 is employed in the formalism of the F12, the ccpVnZ-F12-CABS basis sets, which have been produced so as to minimize the RI error,87 was utilized. Further details about basis sets will be provided as necessary in the discussion of each individual method. Calculation of Redox Potentials. For the molecules under study the oxidation potential is calculated by Nernst equation as

Table 1. Experimental Values for Gas Phase Ionization Energies IE (eV), Redox Potentials in Aqueous Solution E0 (eV), and Henry’s Constants HC (mol/m3pa) aniline p-methoxyaniline p-methylaniline p-chloroaniline N-methylaniline dimethylamine diethylamine pyrrolidine piperidine anisole indole phenol p-methoxyphenol p-methylphenol p-chlorophenol p-cyanophenol dimethylsulfide dimethyldisulfide thioanisole

IE

E0

HC

7.72 7.12 7.46 7.74 8.43 8.24 8.01 8.00 8.20 8.20 7.76 8.49 7.71 8.17 8.44 9.01 8.69 8.18 7.92

1.02 0.79 0.92 0.675 0.95 1.27 1.36 1.26 1.34 1.62 1.08 1.5 1.23 1.38 0.653 1.71 1.39 1.68 1.45

5.2 150 13 10 0.87 0.3 0.39 4.2 2.8 0.029 7.1 28 7.7 10 1400 14000 0.0056 0.0058 0.04

closed-shell (reduced) and singly oxidized species were optimized using the meta-GGA functional TPSS66 with the def2-TZVP(-f) basis set.67 The geometries in aqueous solution are determined self-consistently using the SMD model.68 Vibrational frequency calculations for the optimized geometries are carried out to obtain zero point vibrational energy (ZPVE) and thermochemical contributions. Both wave function based electronic structure methods and a broad range of density functionals are evaluated. Among the wave function based methods are CCSD and CCSD(T), explicitly correlated F12-CCSD(T), as well as a range of methods from the MP2 family, specifically conventional MP2, spin-component-scaled MP2 (SCS-MP2), and orbital optimized MP2 (OO-MP2). For all MP2 calculations the resolution of the identity approximation (RI-MP2) was employed.69−73 Explicitly correlated F12-MP2 calculations were also performed. A range of methods coupled to the LPNO approach were tested, specifically LPNO-CCSD, LPNO-pCCSD, and LPNO-CEPA (coupled electron pair approximation) in the variant known as LPNO-CEPA/1.19,74 For the LPNO calculations, the truncation parameters involved in the calculation (TCutPNO, TCutPairs, and TCutMKN), which define cutoffs for the occupation numbers in the pair natural orbitals, for the estimated pair correlation energies, and for the fitting domain selection, are chosen according to the default “NormalPNO” settings in ORCA (3.33 × 10−7, 1 × 10−4, and 1 × 10−3, respectively). The LPNO calculations were performed with the open-shell code developed by Hansen et al.20 The closed-shell species were also treated with the openshell code since in this specific method, the open-shell LPNOCCSD method only precisely reproduces the restricted LPNOCCSD result in the limit of zero-thresholds. A set of calculations on smaller systems was carried out to determine the optimal reference method that most closely approximates the complete-basis-set (CBS) extrapolated CCSD(T) results but is still applicable to the compounds under study. As will be described in the following, this led to the choice of CCSD(T)/cc-pVTZ combined with extrapolation

° = ESHE ° (Ox) − E R/O

° ΔG R/O (3)

nF

where ΔG°R/O is the Gibbs free energy difference between the reductant and the oxidant in aqueous solution, n is the number of electrons which attribute to the reaction, and F is Faraday constant. In the present study, only one electron transfer is considered, and n is equal to 1 for all the systems. If the employed unit of energy is eV, F is equal to the unit charge e. The superscript circle indicates the standard state of 1 atom ideal gas for gases and a standard state of 1 molal ideal solution for solute. The difference of the Gibbs free energy is calculated by using a thermodynamic (Born−Haber) cycle as ° = IE(0K) + ΔGtrans,vib,rot, ° ° + ESHE ° ΔG R/O ‐TS + ΔΔGsolv (4)

where IE(0K) is the adiabatic ionization energy at 0 K with the consideration of zero point vibrational energy (ZPVE) IE(0K) = ΔUelec + ZPVE

(5)

with ΔUelec being the difference of total electronic energy between reductant and oxidant. The second term in eq 4 involves the thermal contribution of translational, vibrational, and rotational motion at 298 K and the stabilization by entropic contribution. The fourth term is given as a difference of the solvation energy of oxidant and reductant 2274


Article

Journal of Chemical Theory and Computation Table 2. Differences in Predicted Adiabatic Ionization Energies (eV) Obtained with Various Methods from Reference CCSD(T)/cc-pV[5/6]Z Extrapolated Results for a Test Set of Small Molecules CCSD(T)/cc-pV[5/6]Z acetonitrile ammonia formaldehyde hydrogen cyanide hydroxylamine isocyanic acid methanimine methanol MAE

12.345 10.224 10.986 13.645 9.076 11.676 10.056 10.994

° = ΔGsolv ° (Ox) − ΔGsolv ° (Red) ΔΔGsolv

F12-CCSD(T)/ cc-pVTZ-F12

CCSD(T)/cc-pVTZ + LPNO-CCSD/[T/Q]

CCSD(T)/cc-pVTZ + LPNO-CCSD/[Q/5]

−0.040 −0.032 −0.046 −0.038 −0.049 −0.054 −0.033 −0.047 0.042

−0.014 −0.023 −0.013 0.000 −0.026 −0.015 −0.020 −0.017 0.016

−0.014 −0.020 −0.014 −0.006 −0.021 −0.015 −0.015 −0.013 0.015

following the second approach can be obtained if a CCSD(T) calculation with a triple-ζ basis seta level of theory usually feasible on small to medium size moleculesis combined with extrapolation of the correlation energy at the LPNO-CCSD level using basis sets of higher cardinality.75 In this case for example the correlation energy at the basis set limit is approximately represented by

(6)

where ΔGsolv ° represents the Gibbs free energy associated with moving of a single ion from the gas phase into the solution. The last term in eq 4 represents the free energy change for the reference reaction in the standard hydrogen electrode (SHE). In this study, the value −4.28 V is consistently used in all calculations.88,89 Solvation Free Energies. Two implicit solvation models, the conductor-like screening model (COSMO)90 and the SMD68 model that contains a number of empirical parameters that were fitted to experimental data, are tested to find feasible combinations with employed electronic structure methods. The standard definitions of the models were used in each case, with the corresponding dielectric constants for bulk water (80.4 for COSMO and 78.355 for SMD) and a refractive index of 1.33. The reference value for the difference of solvation energies between reduced and oxidized form is calculated from the experimental ionization potential and the oxidation potential based on eq 4. Once the difference in solvation energy and the experimentally derived value of solvation energy for the neutral form is known, the solvation energy for the other species (cation) can be calculated from eq 6.

E corr (CCSD(T)/CBS) ≈ E corr (CCSD(T)/cc‐pVTZ) + E corr (LPNO‐CCSD/cc‐pV[T/Q]Z) − E corr (LPNO‐CCSD/ccpVTZ)

(7)

where Ecorr(CCSD(T)/cc-pV[T/Q]Z) denotes the extrapolated energy for correlation part of LPNO-CCSD with ccpVTZ and cc-pVQZ basis sets. Here we compare F12-CCSD(T)/cc-pVTZ-F12 and CCSD(T)/cc-pVTZ with LPNO-CCSD/cc-pV[Q/5]Z extrapolation, since the latter combination of basis sets is still feasible owing to the computational efficiency of the LPNO approach. To perform this comparison we constructed a separate test set of small molecules, for which CCSD(T) calculations with up to sextuple-zeta basis sets are feasible. These small molecules are listed in Table 2, along with the reference CCSD(T)/cc-pV[5/ 6]Z values and a comparison of other methods. In terms of absolute deviations, all methods slightly underestimate the ionization energies compared to the reference. Overall, the comparison in terms of average errors is in favor of the composite extrapolation method. It must be noted that the F12-CCSD(T) calculations do recover more correlation energy in absolute terms compared to LPNOCCSD/cc-pV[Q/5]Z extrapolated results, but they do not appear to recover a consistently similar amount of correlation energy for both the reduced and the oxidized species. Specifically, for the small molecules of the above test set, the absolute F12-CCSD(T)/cc-pVTZ-F12 correlation energies approach, on average, the CCSD(T)/cc-pV5Z energies for the reduced forms and the CCSD(T)/cc-pV6Z energies for the oxidized forms. Although the absolute deviations from the reference values are small for most intents and purposes, the above behavior results in somewhat inferior performance for ionization energies compared to the method that employs correlation energy extrapolation with LPNO-CCSD. In the latter case the ionization energies are already converged with cc-pV[T/Q]Z extrapolation, since only marginal improvement is achieved with the more demanding cc-pV[Q/5]Z extrapolation. Based on the above results, reference values for the organic molecules in the following will be obtained with

■

RESULTS AND DISCUSSION Geometry Optimizations. As described above, each species is independently optimized for the reduced and oxidized forms and those with implicit solvation. Table S1 lists root-mean-square deviations (RMSD) between reduced (neutral) and oxidized forms in gas phase and in solution compared with those of the reduced form in gas phase for bond lengths, bond angles, and dihedral angles. The structures of reduced and oxidized forms are presented superimposed in Figure S1. The individual results and average deviations suggest that in general the change of oxidation state has a greater effect on geometries (primarily in dihedral angles) than the modeling of the environment, that is, gas phase or implicit solvation. Definition of Reference Electronic Structure Method. The electronic structure method that will provide the reference values for the present study is required to reproduce CCSD(T)/CBS ionization energies as close as possible. Among the methods that can be used for this purpose and that are applicable to the size of the molecules included in the present test set are the explicitly correlated F12-CCSD(T) with a medium-size F12-optimized basis set and approaches that combine a CCSD(T) calculation performed with a mediumsize basis set with correlation energy extrapolation using a computationally cheaper method such as MP2 or LPNOCCSD. Based on a recent evaluation study, optimum results 2275


Article


Table 3. Performance of Wave Function Based Methods for the Adiabatic Ionization Energy (eV) against the Reference Values of the Present Test Set, Ordered According to Increasing Mean Unsigned Error (MUE) method

main basis set

MUE

MSE

UEmax

UEmax species

CCSD(T) LPNO-pCCSD OO-RI-MP2 CCSD LPNO-CCSD LPNO-CEPA F12-CCSD(T) SCS-MP2 F12-MP2 MP2 HF

cc-pV[D/T]Z cc-pV[T/Q]Z cc-pV[T/Q]Z cc-pV[D/T]Z cc-pV[T/Q]Z cc-pV[T/Q]Z cc-pVDZ-F12 cc-pV[T/Q]Z cc-pVTZ-F12 cc-pV[T/Q]Z cc-pVTZ

0.02 0.05 0.07 0.08 0.09 0.12 0.13 0.54 0.63 0.65 1.42

0.02 −0.02 0.06 −0.08 −0.06 −0.08 −0.13 0.54 0.63 0.65 −1.42

0.05 0.48 0.13 0.12 0.51 0.55 0.17 1.19 1.32 1.39 1.63

thioanisole chlorophenol pyrrolidine aniline chlorophenol chlorophenol chloroaniline cyanophenol cyanophenol aniline N-methylaniline

CCSD(T)/cc-pVTZ + LPNO-CCSD/cc-pV[Q/5]Z extrapolation. Ionization Energies. Statistical metrics (mean unsigned error MUE, mean signed error MSE, and maximum unsigned error UEmax) for the performance of the various wave function and DFT methods compared to the reference electronic ionization energies obtained by CCSD(T)/cc-pVTZ + LPNOCCSD/cc-pV[Q/5]Z extrapolated energies are provided in Tables 3 and 4. Note that here only the electronic component

Table 5. Demonstration of Basis Set Dependence of Selected Wave Function Based Methods for the Mean Unsigned Error MUE (eV) in the Prediction of Adiabatic Ionization Energies of the Molecules under Study method CCSD(T) CCSD LPNO-CCSD LPNOpCCSD LPNO-CEPA MP2 SCS-MP2 OO-MP2

Table 4. Performance of Density Functional Theory Methods with the def2-TZVPP Basis Set for the Adiabatic Ionization Energy (eV) against the Reference Electronic Energy Values of the Present Test Set, Ordered According to Increasing MUE functional

MUE

MSE

UEmax

UEmax species

M06-2X mPW2PLYP B2PLYP M06 PBE0 B3LYP PBE TPSSh M06-L TPSS

0.11 0.24 0.25 0.25 0.32 0.38 0.40 0.41 0.42 0.45

−0.01 −0.15 −0.16 −0.18 −0.24 −0.30 −0.32 −0.34 −0.35 −0.38

0.94 0.88 0.88 0.68 0.78 0.68 0.73 0.69 0.68 0.66

phenol phenol phenol phenol phenol phenol phenol phenol phenol phenol

of the adiabatic ionization energy is considered and is defined as the energy required to generate the most stable oxidized species. Table 5 shows the basis set dependence of selected methods, while Figure 2 provides an overall graphical representation of the performance of the methods. A more extended version of Table 3, in which all calculations for extrapolated and nonextrapolated energies for the specific basis set are listed, is given as Table S2, and the ionization energies of all test molecules for the methods listed in Table 3 as well as for all DFT methods are summarized in Tables S3 and S4 of the Supporting Information. For all wave function based methods where this was possible, the results reported above were obtained with [T/Q] basis set extrapolation. We note however that the conclusions regarding each method, both in terms of individual errors and in terms of average performance, are identical when the considerably cheaper [D/T] extrapolation is used instead. On the topic of basis set convergence (Table 5), it should be pointed out that results obtained with the cc-pVTZ basis set are always inferior

ccpVDZ

ccpVTZ

ccpVQZ

cc-pV[D/T]Z

cc-pV[T/Q]Z

0.41 0.45 0.43 0.41

0.14 0.22 0.21 0.17

0.10 0.06

0.02 0.08 0.09 0.05

0.09 0.05

0.47 0.29 0.23 0.39

0.23 0.45 0.35 0.11

0.10 0.57 0.46 0.04

0.11 0.62 0.52 0.06

0.12 0.65 0.54 0.07

Figure 2. Bar graph of MUEs in electronic adiabatic ionization energies (eV) for wave function based methods (blue bars) and DFT methods (orange bars) tested in the present study. The values are obtained with the basis sets and basis set extrapolations indicated in Tables 3 and 4 (def2-TZVPP for all density functionals), except where explicitly indicated otherwise by “[D]” or “[T]” for cc-pVDZ and ccpVTZ.

to extrapolated results, even using the simplest possible [D/T] extrapolation. LPNO-based methods benefit to a significant extent from basis set extrapolation because the extrapolation 2276


Article


over MP2 or SCS-MP2. Given that each OO-MP2 iteration has a comparable cost to that of an RI-MP2 relaxed density calculation, the OO-MP2 approach is much more expensive than RI-MP2. In fact, the CPU times required to obtain ccpV[T/Q]Z extrapolated energies with OO-RI-MP2 range from 5 to 10 times more than those of corresponding LPNO-CCSD calculations on the same hardware. Therefore, with the exception of specific applications, it is difficult to see how this approach, despite its robustness for the present test set, can be competitive in general applications with computationally more convenient and better scaling LPNO approaches, in particular once linear scaling open-shell LPNO-CCSD(T) variants become available. As can be seen from the results presented above, all LPNO methods perform extremely well. This is particularly true for pCCSD, the parametrized variant of CCSD proposed by Huntington and Nooijen.93 The good performance of pCCSD has been previously demonstrated in other applications.94,95 LPNO-CCSD itself has somewhat larger MUE of 0.09 eV, practically converging to the canonical CCSD value, while LPNO-CEPA reaches a MUE of 0.12 eV. Focusing on more specific aspects, an interesting comparison is that between the bare CCSD(T)/cc-pVTZ values and those obtained with the explicitly correlated F12-CCSD(T)/ccpVDZ-F12 (note that the double-ζ nomenclature is used here rather liberally and the F12 basis set is actually similar in size to cc-pVTZ). As expected, the F12 method recovers more correlation energy in all species involved in the test set. However, these absolute gains do not translate to significant improvements in the ionization energies. Selected F12CCSD(T) calculations with the larger cc-pVTZ-F12 basis set, a substantially more expensive task that borders the limits of applicability for the larger molecules of the present set, demonstrate a nontrivial difference between the DZ-F12 and the TZ-F12 results in the final ionization energies: the larger basis set reduces the F12-CCSD(T) errors in ionization energies by up to a factor of 3 (see the SI for details). It therefore becomes clear that the F12-CCSD(T) calculations with the cc-pVDZ-F12 basis set are not yet sufficiently converged to deliver the same level of accuracy as that of the [D/T] extrapolated CCSD(T) calculations, or of the much cheaper extrapolated LPNO methods, which can be thus considered as the optimum combination of accuracy and efficiency for ionization energies in the present test set. The best performance among the density functionals is shown by M06-2X, which represents the only DFT method that achieves a MUE close to 0.1 eV. The M06-2X method is a hybrid meta-GGA functional with 54% of Hartree−Fock exchange and 32 empirically optimized parameters.81 This functional was shown to perform well for ground state properties and also for valence and vertical and Rydberg electronic transition energies96,97 for benchmark sets of molecules consisting of main group elements. For ionization energies of the present test set the MUE and MSE of M06-2X are respectively 0.11 eV and −0.01 eV, thus this level of performance is directly comparable with the best-performing wave function based methods. The two double hybrid functionals, B2PLYP and mPW2PLYP, are in second place with respect to their mean deviations from the reference and perform almost identically with respect to each other. The M06 functional, a version of the M06 family with “normal” percentage of HF exchange (27%), essentially matches the performance of the double hybrids. All

tends to reduce the PNO specific errors. This can be traced back to the fact that the LPNO errors are highly systematic and approach the correlation energy always from below (in absolute value).75 In line with previous observations,22 it must be emphasized that the use of a double-ζ basis set as demonstrated here for example with the behavior of CCSD(T)/cc-pVDZ, essentially guarantees results that are inferior in terms of mean unsigned error and error spread to most DFT methods. As an overall impression, there is a very widespread of results among the methods tested, which range all the way from catastrophic failures to consistent and quantitatively accurate methods that yield mean errors of less than 0.1 eV on average. MP2 fares so badly for the ionization energies considered here, that the results can only be described as unusable. Various problems of MP2 with the description of the electronic structure, energetics, and other properties of radical species are documented in the literature.91 As expected from the fact that the ionization process studied here involves a change from closed-shell to open-shell species, the inability of Hartree−Fock to account correctly for the differential correlation energy in the two cases leads to gross systematic underestimation of the ionization energies by more than 1.4 eV on average. Beyond rendering HF results of no practical use whatsoever, this has important implications for the performance and reliability of post-Hartree−Fock correlation methods that cannot sufficiently compensate for such a gross shortcoming of the reference wave function. The reason for the failure of MP2 is therefore the extremely poor Hartree−Fock reference. This is also why the spincomponent-scaled variant (SCS-MP2), which typically corrects shortcomings of canonical MP2 with respect to the correlation energy, offers somewhat improved performance over MP2 but still falls far short of eliminating the enormous errors observed in the present set since it cannot address the underlying enormous HF errors. In contrast to conventional MP2, in the SCS-MP2 method the opposite-spin and the same-spin correlation energy are scaled differently. Specifically, Grimme proposed the coefficient 1.2 for the same-spin and 0.3 for the opposite-spin correlation, demonstrating superior performance of SCS-MP2 for ground state energies compared with conventional MP2.92 The basis set limit is reached quite fast, both by the explicitly correlated F12 approach and by conventional extrapolation, but the errors of MP2 and SCSMP2 become even larger at this limit. Therefore, the performance of these methods is dominated by the extremely poor HF reference, rendering second-order perturbation methods practically unusable for calculations of this type. In stark contrast, the orbital-optimized MP2 method gives a staggering improvement in accuracy over conventional MP2. In fact, the orbital-optimized approach surpasses in accuracy even the CCSD(T)/cc-pVTZ result (but not the extrapolated CCSD(T) values). The basis set dependence is also quite small, with the limit being reached already with a [D/T] extrapolation. The above discussion helps to better understand the superior performance of OO-MP2 over canonical MP2, since the fundamental idea here is not only to minimize the Hylleraas functional with respect to the MP2 amplitudes but also to additionally minimize the total energy with respect to orbital changes. More details on the singles amplitudes as derived from the variational OO-MP2 procedure and an overview of the HF, MP2, and OO-MP2 errors are provided in the Supporting Information (Table S5). The major shortcoming of OO-MP2 is, however, its significantly increased cost 2277


Article


Table 6. Errors in Ionization Energies (eV) Calculated by Wave Function Methods against Experimental Values Listed in Table 1 method

main basis set

MUE

MSE

UEmax

UEmax species

CCSD(T) LPNO-CCSD LPNO-pCCSD LPNO-CEPA F12-CCSD(T) OO-RI-MP2 F12-MP2 MP2 HF

cc-pV[D/T]Z cc-pV[T/Q]Z cc-pV[T/Q]Z cc-pV[T/Q]Z cc-pVDZ cc-pV[T/Q]Z cc-pVTZ-F12 cc-pV[T/Q]Z cc-pVTZ

0.07 0.08 0.08 0.11 0.11 0.11 0.65 0.68 1.40

0.04 −0.04 0.00 −0.06 −0.11 0.08 0.65 0.67 −1.40

0.13 0.47 0.44 0.51 0.26 0.20 1.38 1.45 1.63

piperidine chlorophenol chlorophenol chlorophenol piperidine DMA cyanophenol aniline methylaniline

local density functionals tested in this study give larger errors, although it should be pointed out that even the worstperforming DFT functionals fair substantially better than MP2. The lowest MUE among GGA functionals is given by PBE. Focusing on individual molecules, the largest deviation is observed for phenol for all functionals, whereas no single molecule was identified as particularly problematic for the wave function methods. The S2 expectation values for the oxidized compounds (radical cations) in gas phase and solution are respectively listed in Tables S6 and S7 for the ten tested density functionals. Relatively large spin contamination in the gas phase, up to ∼0.90 against the ideal value S(S+1) = 0.75, is observed for the two double-hybrid density functionals, B2PLYP and mPWP2LYP. Overall, the spin contamination seems approximately proportional to the percentage of Hartree−Fock exchange; PBE0, M06, M06-2X, B2LPYP, and mPW2PLYP, with 25%, 27%, 54%, 53%, and 55% of Hartree−Fock exchange, respectively, show a clear trend toward larger spin contamination with increasing the Hartree−Fock exchange. Therefore, the double hybrids have similar problems with the “reference” determinant as MP2, and hence the nonlocal, perturbative, MP2-like correlation component of the double hybrids suffers somewhat from the spin contamination in the Kohn−Sham determinant (we note in passing that the S2 expectation values for the double hybrids were calculated only from the Kohn− Sham determinant and hence do not include the two-electron contribution from the perturbative correction). So far we have compared the performance of all methods against the highly accurate electronic structure reference, CCSD(T)/cc-pVTZ with LPNO-CCSD/cc-pV[Q/5]Z extrapolation. The comparison to available experimental ionization energies is also important, although it should be kept in mind that this evaluation is necessarily less rigorous because of the experimental uncertainties, the different sources of experimental data, and the need to correct the electronic energies with an estimate of the zero-point vibrational energy (ZPVE), for experimental IEs that are referenced to 0 K. For the systems of the current study, the ZPVE can only be calculated with realistic effort using DFT or MP2 methods. Based on the poor performance of MP2, DFT is hence the method of choice. For the purposes of the present comparison with available experimental values, we have consistently applied to all electronic energies the ZPVE correction calculated at the TPSS/def2-TZVP(-f) level used for geometry optimizations. Table S8 gives the contributions of ZPVE, demonstrating that the contribution is small but certainly not negligible, with absolute values of up to 0.06 eV. Comparisons with experimental values are provided in Tables 6 and 7.

Table 7. Errors in Ionization Energies (eV) Calculated by Ten Density Functionals against Experimental Values functional

MUE

MSE

UEmax

UEmax species

M06-2X B2PLYP mPW2PLYP M06 PBE0 B3LYP PBE TPSSh M06-L TPSS

0.11 0.23 0.23 0.24 0.31 0.36 0.38 0.40 0.41 0.43

0.01 −0.14 −0.13 −0.16 −0.22 −0.29 −0.30 −0.32 −0.33 −0.36

1.00 0.94 0.94 0.74 0.84 0.74 0.79 0.75 0.74 0.72


Although the statistical errors change slightly from those calculated against the reference electronic values, the order of error for the tested methods does not change except for methods which already showed a low MSE from both computational and experimental references. Overall, the performance of individual methods is dominated by the description of the electronic structure, and our conclusions regarding their relative performance with respect to the extrapolated CCSD(T)/cc-pVTZ+LPNO-CCSD/cc-pV[Q/ 5]Z reference do not change when we compare with the available experimental values. As an aside, we mention a methodological aspect relating to what is usually referred to in DFT as the “IP-theorem” (equivalent to Koopmans’ theorem in Hartree−Fock theory), according to which the vertical ionization energy would correspond to the negative of the energy of the highest occupied molecular orbital (HOMO) in an “exact” Kohn− Sham DFT calculation. The relationship between the HOMO energy and oxidation potential has been examined previously (for example, ref 98). We compare three types of IE: (i) adiabatic IEs used in the present work, (ii) vertical IEs (Table S9), and (iii) the negative of the HOMO energy (Table S10). A summary of these values for the DFT methods of the present study is plotted against the reference values in Figure S2. When considering vertical IEs, all functionals except M06-2X give lower MUEs with respect to the reference compared with the adiabatic IE values, whereas the MUE of M06-2X is 0.32 eV for the vertical compared to 0.11 eV for the adiabatic IEs. Values obtained from the HOMO energies severely underestimate the IEs. Even though the trends are reasonably reproduced, the absolute values can be so far off that the approximation is not a practical alternative. This is to be expected from the work of van Leeuwen and Baerends,99 who blame the erroneous longrange behavior of the currently employed functionals for 2278


Article

Journal of Chemical Theory and Computation Table 8. Difference of Solvation Energy (eV) of Reduced and Oxidized Species SMD PBE TPSS M06-L B3LYP PBE0 TPSSh M06 M06-2X B2PLYP mPW2PLYP

COSMO

MUE

MSE

UEmax

UEmax species

MUE

MSE

UEmax

UEmax species

0.22 0.22 0.22 0.22 0.23 0.22 0.22 0.22 0.22 0.22

0.15 0.16 0.14 0.12 0.15 0.15 0.13 0.14 0.12 0.12

1.13 1.13 1.11 1.09 1.12 1.12 1.1 1.11 1.1 1.09

chlorophenol chlorophenol chlorophenol chlorophenol chlorophenol chlorophenol chlorophenol chlorophenol chlorophenol chlorophenol

0.40 0.40 0.39 0.38 0.39 0.39 0.38 0.39 0.38 0.38

0.39 0.38 0.37 0.36 0.37 0.38 0.36 0.37 0.37 0.36

1.65 1.68 1.67 1.63 1.66 1.67 1.56 1.63 1.68 1.68


Table 9. Comparison of COSMO and SMD for Errors in the Solvation Energy (eV) of the Neutral (Reduced) Species Obtained by Different Functionals SMD PBE TPSS M06-L B3LYP PBE0 TPSSh M06 M06-2X B2PLYP mPW2PLYP

COSMO

MUE

MSE

UEmax

UEmax species

MUE

MSE

UEmax

UEmax species

0.07 0.07 0.09 0.07 0.05 0.07 0.07 0.06 0.07 0.07

0.06 0.07 0.09 0.06 0.04 0.06 0.07 0.04 0.06 0.06

0.18 0.18 0.20 0.17 0.16 0.17 0.18 0.16 0.17 0.17

chlorophenol chlorophenol methoxyaniline chlorophenol chlorophenol chlorophenol chlorophenol N-methylaniline N-methylaniline N-methylaniline

0.06 0.06 0.06 0.06 0.06 0.06 0.05 0.06 0.06 0.06

0.00 0.01 0.03 0.01 −0.01 0.00 0.01 −0.01 0.00 0.00

0.11 0.11 0.13 0.11 0.10 0.11 0.11 0.11 0.10 0.10

pyrrolidine pyrrolidine chlorophenol pyrrolidine piperidine pyrrolidine pyrrolidine methoxyphenol chlorophenol chlorophenol

Table 10. Comparison of COSMO and SMD for Errors in the Solvation Energy (eV) of the Oxidized Species Obtained by Different Functionals SMD PBE TPSS M06-L B3LYP PBE0 TPSSh M06 M06-2X B2PLYP mPW2PLYP

COSMO

MUE

MSE

UEmax

UEmax species

MUE

MSE

UEmax

UEmax species

0.27 0.27 0.28 0.26 0.26 0.27 0.27 0.26 0.26 0.26

0.21 0.22 0.23 0.19 0.19 0.21 0.2 0.18 0.18 0.18

1.30 1.29 1.30 1.25 1.27 1.28 1.27 1.25 1.26 1.25

chlorophenol chlorophenol chlorophenol chlorophenol chlorophenol chlorophenol chlorophenol chlorophenol chlorophenol chlorophenol

0.41 0.41 0.42 0.39 0.39 0.4 0.4 0.39 0.39 0.39

0.39 0.39 0.39 0.37 0.36 0.38 0.37 0.36 0.37 0.36

1.67 1.70 1.72 1.65 1.67 1.69 1.59 1.63 1.70 1.69


standard statistical thermodynamics on the basis of harmonic vibrational frequencies, are shown in Table S8. Both thermal and entropic contributions are small for the ionization energies but comparable to the ZPE corrections. The absolute entropic contribution of ca. 0.03 eV is on average larger than the thermal contribution (ca. 0.01 eV) in the present systems. These contributions may become more important when one considers larger systems where the cancellation may not be as complete. The computed ionization energies at 298 K are given in Table S11. As can been seen from the magnitude of thermal and entropic contributions, there are in practice no large changes from those at 0 K. Overall, it is concluded that the major determinant for the quality of ionization energies is the electronic structure method employed. Fast implementations of wave function based

HOMO energies that are far too high (for the exact functional Janak’s theorem suggests that the HOMO energy coincides with the vertical IE100). Since the Hartree−Fock exchange potential has the correct asymptotic behavior, it is not surprising that the functionals with the highest percentage of HF exchange fare best in this comparison. Similar conclusions were reached in recent studies of ionization energies and electron affinities in a series of acceptor molecules, where it was shown that long-range corrected hybrid functionals in which the range-separation parameter was nonempirically tuned to obey the IP-theorem perform very reliably for vertical IEs and EAs.101,102 Coming back to the actual values used in the present work, the ionization energies at 298 K are obtained by further adding the thermal and entropic contributions to the 0 K ionization energies. Each of these contributions, as determined by 2279


Article


Table 11. Errors in Redox Potentials (V), Computed by Wave Function Based Methods in Combination with TPSS Thermal Contributions and TPSS-SMD Solvation Energies method

main basis set

MUE

MSE

UEmax

UEmax species

LPNO-CCSD LPNO-pCCSD CCSD(T) F12-CCSD(T) LPNO-CEPA OO-RI-MP2 SCS-MP2 F12-MP2 MP2 HF

cc-pV[T/Q]Z cc-pV[T/Q]Z cc-pV[D/T]Z cc-pVDZ-F12 cc-pV[T/Q]Z cc-pV[T/Q]Z cc-pV[T/Q]Z cc-pVTZ-F12 cc-pV[T/Q]Z cc-pVTZ

0.28 0.28 0.28 0.30 0.30 0.30 0.69 0.78 0.80 1.29

0.07 0.12 0.16 0.01 0.06 0.19 0.67 0.76 0.79 −1.29

1.02 1.05 1.09 0.91 0.95 1.12 2.25 2.38 2.38 1.73

cyanophenol cyanophenol cyanophenol cyanophenol cyanophenol cyanophenol cyanophenol cyanophenol cyanophenol N-methylaniline

Table 12. Errors in Redox Potentials (V), Obtained for the Density Functionals Included in the Present Study, with the Two Different Implicit Solvation Models SMD M06-2X M06 mPW2PLYP B2PLYP PBE0 PBE B3LYP TPSSh M06-L TPSS

COSMO

MUE

MSE

UEmax

UEmax species

MUE

MSE

UEmax

UEmax species

0.28 0.30 0.33 0.33 0.36 0.39 0.40 0.41 0.41 0.43

0.12 −0.06 −0.04 −0.04 −0.10 −0.18 −0.19 −0.20 −0.22 −0.23

1.39 1.10 1.30 1.31 1.23 1.16 1.10 1.13 1.09 1.10


0.34 0.27 0.32 0.32 0.33 0.34 0.34 0.35 0.35 0.36

0.29 0.11 0.14 0.14 0.07 0.00 −0.01 −0.03 −0.05 −0.06

1.54 1.26 1.47 1.47 1.38 1.32 1.27 1.29 1.24 1.26


correlated methods can ensure affordable and “black box” access to highly consistent and accurate predictions. Solvation Energies. Table 8 lists the calculated difference of solvation energies between the neutral and the cationic species ΔΔGsolv for the density functionals tested in the present study in combination with the two implicit solvation models, SMD and COSMO (we do not report such values for wave function based methods because in this case no fully relaxed correlation contribution to the solvation energy is available). This difference of solvation energy is one of the dominant properties that control the accuracy in the calculation of redox potentials. As described in the Methodology section, the reference ΔΔGsolv values are derived using the experimental redox potentials. For both implicit solvation methods examined here, SMD and COSMO, the magnitude of error does not depend significantly on the choice of density functional. The calculated ΔΔGsolv with COSMO is almost always overestimated compared to the reference values, except for diethylamine and dimethyldisulfide. On the other hand, SMD shows smaller errors for ΔΔGsolv in general and does not appear to have a bias, as reflected by the low MSE values compared to COSMO. The details of errors of individual molecules in COSMO and SMD are summarized in Table S12 and Table S13. The difference of solvation energy, ΔΔGsolv, is decomposed into the contribution from the neutral species and the oxidized cationic species, and these contributions are respectively given in Tables 9 and 10. The details of solvation energies of reduced forms are summarized in Tables S14 (SMD) and S15 (COSMO) and oxidized forms in Tables S16 (SMD) and S17 (COSMO). The COSMO and SMD solvation models show comparable performance for the neutral species of the

present test set, but SMD shows consistently better performance compared to COSMO for all density functionals for the cationic species. It is precisely this better performance for the oxidized form that results in overall better performance of the SMD model for the differential solvation energies ΔΔGsolv reported in Table 8. The maximum error UEmax is quite dispersed among the molecules of the test set for the neutral species, but it is also consistently low and there is no obvious dramatic failure. In stark contrast, the oxidized (cationic) radical species display much larger errors, and UEmax is consistently observed for chlorophenol in the case of SMD and phenol in the case of COSMO, regardless of the specific functional employed. As can be seen from Tables S16 and S17, significant deviations are observed in substituted phenols for the oxidized species. The largest error is in phenol with MUE over 1.5 eV on average over ten density functionals; this magnitude of error corresponds to roughly twice the error in ionization energies. The large error in the solvation energy of the oxidized form is presumably related to the treatment of the environment, and improvements can be achieved by explicitly including hydrogen bonded water molecules, as discussed in a later section. Redox Potentials. The final step is the combination of all individual energetic components (electronic energies, thermodynamic corrections, and solvation contributions) to obtain the values for the redox potentials. Table 11 shows the statistical errors of redox potentials compared to experiment for wave function based methods. Here the ionization potential is estimated by the respective wave function method, the thermal corrections are adopted from the TPSS values (as for the comparison with experimental IEs, since this is used for geometry optimizations and frequencies calculations), and 2280


Article


approach when the explicitly solvated cluster is used in combination with an implicit solvent model) has indeed been reported to improve the accuracy of predicted redox potentials.108,109 Such a treatment is also applicable for a QM/MM approach to solvation,110 which can be combined with a force field that accounts for electronic polarization. We note that an approach to describe charge transfer111 has been suggested in which the charge transfer is realized in the framework of the molecular fragmentation method.112 Here we test this cluster−continuum approach in combination with SMD to see how much the description is improved for the DFT methods used in the present work in the case of phenolwater clusters in which a hydrogen bonding network is established in the vicinity of the hydroxyl group.113,114 For a higher level and more complete treatment of this problem that couples MD sampling with the equation-of-motion ionization potential coupled cluster and effective fragment potential methods, we refer the reader to the excellent study of Krylov and co-workers.35 Figure 3 shows the optimized structures for phenol−water clusters at the TPSS/def2-TZVP(-f) level for the neutral form;

similarly the solvation energies are those computed by the TPSS functional in conjunction with SMD solvation model. These are expected to be the best practical choices for the additional nonelectronic energy contributions given that the corresponding thermodynamic corrections and solvation energies cannot presently be obtained directly with the correlated electronic structure methods listed in Table 11. The increase in the MUE for the redox potential over that of ionization energies becomes conspicuous in highly accurate methods such as CCSD(T)/cc-pV[D/T]Z; however, for methods that suffer from very large errors in ionization energies (for example, HF), lower errors may be observed due to error cancellation. The range of the MUE is similar between two solvation models for ten density functionals; 0.28−0.43 V for SMD and 0.27−0.36 V for COSMO (see Table 12). Some methods show a lower MUE than the ionization energies, meaning that a cancellation of errors occurs in the process of computing the redox potentials. Notice that despite the very low MUE of 0.08 eV for the ionization energies obtained with M06-2X, an excellent result that compared very well even with much more expensive wave function based methods, a mean absolute error close to 0.3 V is still observed for the redox potential, indicating that the dominant errors originate from the solvation energy. It is of course expected that the modeling of aqueous solution in the present case can present particular difficulties and expose weaknesses of the implicit solvation modelswhich is precisely why this choice was made. In terms of comparison to experimental values, we would like to emphasize that we do not suggest the present values exhaust the capabilities of computational modeling, because case-by-case consideration of tautomerism, acid−base equilibria, and hydrogen bonding (as shown using the example of phenol in the following section) might selectively improve the results. In any case, the best performance is seen here for M06 combined with COSMO and M06-2X combined with SMD. The largest MUEs are observed for TPSS and M06-L, for both SMD and COSMO. This does not necessarily imply a requirement for Hartree−Fock exchange, since, for example, TPSSh and B3LYP do not perform better on average than the nonhybrid functionals. Finally, we note that the MSEs are smaller than MUEs for most cases, indicating that the deviations are not systematic. Cluster Approach for Phenol. As pointed out in the preceding sections, the errors for phenol and some substituted phenols for both the ionization energy and the oxidation potential from reference values were particularly pronounced compared to other compounds. Phenol and its derivatives have been the subject of numerous theoretical studies,35,103−105 and the challenges in obtaining accurate ionization energies and oxidation potentials are well documented. For example, Klein and Lukeš studied 37 substituted phenols and found that DFT considerably underestimates the ionization energy with an average deviation of 2.1 eV.105 Large errors in redox potentials have been observed for a set of substituted phenols, where the MUE for 18 species had been calculated as 2.44 V at the BPW91/cc-pVDZ level,103 an error significantly larger than that for substituted anilines (0.09 V).106 It is well-known that hydrogen bonding interactions are not sufficiently well described by implicit solvation models.107 Therefore, it is expected that a simple way to improve the accuracy of predictions would be to include one or more explicit solvent molecules that establish specific interactions with the solute. This cluster approach (or “cluster−continuum”

Figure 3. Neutral (reduced) forms of optimized phenol−water clusters.

starting geometries were the ones previously reported as stable structures in the gas phase.113,114 For single water molecule containing clusters, the reduced and oxidized forms adopt similar geometries; however, there are obvious differences between the two species for clusters containing two and three water molecules. Table 13 provides calculated ionization energies and redox potentials for phenol-water clusters. The IE is reduced by 1.54−2.1 eV for C6H5OH-H2O, 1.67−2.29 eV for C6H5OH(H2O)2, and 1.64−2.23 eV for C6H5OH-(H2O)3 compared to the naked phenol, as can be expected from the stabilization of the highest molecular orbital. The MUE in the oxidation potential is reduced from 1.19 to 0.15 V for C6H5OH-H2O, to 0.34 V for C6H5OH-(H2O)2, and to 0.25 V for C6H5OH(H2O)3. Thus, the improvements offered by the cluster approach are dramatic. The larger error for C6H5OH-(H2O)2 and C6H5OH-(H2O)3 compared with the single water cluster may be indicative of complications that arise when a larger number of solvent molecules are explicitly included in the model and could be related to the imbalance of explicit and implicit modeling of hydration, as well as the absence of configurational sampling for the explicitly included water molecules. The ZPVE, thermodynamic, and entropic contribution for the redox potentials are respectively 0.01, − 0.02, 0.05 eV for C6H5OH-H2O, − 0.05, 0.01, − 0.04 eV for C6H5OH-(H2O)2, 2281


Article


Table 13. Ionization Energies (eV) and Redox Potentials (V) for the Bare Phenol and the Three Phenol−Water Clusters E0

IE PBE TPSS M06-L B3LYP PBE0 TPSSh M06 M06-2X B2PLYP mPW2PLYP MUE

no water

H2O

(H2O)2

(H2O)3

no water

H2O

(H2O)2

(H2O)3

9.27 9.20 9.22 9.22 9.32 9.23 9.22 9.48 9.42 9.42

7.58 7.51 7.55 7.59 7.65 7.55 7.68 7.86 7.32 7.35

7.39 7.31 7.42 7.40 7.46 7.34 7.55 7.71 7.13 7.17

7.44 7.36 7.46 7.45 7.51 7.40 7.58 7.74 7.19 7.23

2.66 2.60 2.59 2.60 2.73 2.63 2.60 2.89 2.81 2.80 1.19

1.42 1.34 1.37 1.40 1.50 1.39 1.50 1.69 1.13 1.16 0.15

1.17 1.08 1.19 1.16 1.25 1.13 1.32 1.48 0.91 0.95 0.34

1.29 1.20 1.28 1.27 1.36 1.24 1.40 1.57 1.00 1.05 0.25

and −0.10, − 0.01, 0.02 eV for C6H5OH-(H2O)3. Contributions of C6H5OH-(H2O)2 and C6H5OH-(H2O)3 are substantially larger than those of naked phenol (0.01, 0.00, and 0.02 eV), especially the effects of ZPVE and entropic contribution are substantial for the specific structures shown in Figure S3. Table S21 lists the solvation energies for the neutral and cationic species of the three types of phenol-water clusters. The solvation energies of the neutral species are increased by more than 0.3 eV for all three types of clusters. On the other hand, the solvation of the radical cation does not show significant differences between the bare and the cluster radical, which implies that some cancellations occur between the increase of the solvation energy by water molecules and the reduction of the energy by the delocalization of the charge distribution in the cluster models.

the solvation of the oxidized form (cationic species). Overall, the solvation contribution emerges as the limiting factor in achieving a level of accuracy for the redox potentials similar to what can be achieved for ionization energies. In practice, therefore, the errors in solvation energy remove much of the advantage in employing more consistently accurate wave function based methods since they overshadow the improved electronic energies that these methods yield compared to DFT. Further improvement in this direction can be achieved with a fundamental improvement in the performance of the implicit solvation models or through explicit treatment of solvation. However, the latter approach comes with its own set of challenges, but its advantages are clearly seen already with a very simple treatment of minimally hydrated phenol. A main target for the near future is to wield the power of the new DLPNO methods, which promise to deliver the same level of accuracy as the best performing electronic structure methods presented in this work, but for a cost comparable to that of DFT,23 so that the applicability of accurate correlation methods can eventually be extended to the treatment of explicitly solvated systems.

■

CONCLUSIONS A test set of organic molecules was used to evaluate wave function and density functional methods for ionization energies and, in combination with implicit solvation models, for aqueous oxidation potentials. The results confirm that accuracy in electronic energies better than 0.1 eV on average can only be obtained with wave function based methods, although the M062X functional comes very close to this level of performance and is the only DFT method that can be considered directly competitive with more expensive wave function based correlated methods. Among the latter, MP2 and its variants (with the exception of orbital-optimized MP2) fail badly owing to the poor Hartree−Fock reference and are unusable for the study of such processes. LPNO-CCSD and LPNO-CEPA are among the top performers compared with the canonical CCSD(T) and the extrapolated reference values for electronic ionization energies. Their excellent performance and relatively low cost make them an obvious choice even for routine use when accurate energetics and tight error control are desired. Further improvements are expected from a consistent open shell implementation of the domain-based LPNO (DLPNO) methods115 that is currently under development in our laboratory. The calculation of redox potentials using the Born−Haber thermodynamic cycle involves as principal contributions the gas phase ionization energy and the solvation energy. Two continuum solvation models were tested for the latter, SMD and COSMO, in combination with DFT methods. Both COSMO and SMD perform similarly for the solvation energy of the reduced forms of the molecules (neutral species), but SMD provides distinctively significantly better performance for

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.6b00252.

■

Figures S1−S3, Tables S1−S21, Cartesian coordinates (PDF)

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. *E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work is supported by the Cluster of Excellence RESOLV (EXC 1069) funded by the Deutsche Forschungsgemeinschaft. Support by the Max Planck Society is gratefully acknowledged. The authors thank Robert Izsak and Dimitrios Liakos for discussions on explicitly correlated methods and extrapolation techniques and Marius Retegan and Vera Krewald for insightful comments. 2282


Article


■

(40) Blinco, J. P.; Hodgson, J. L.; Morrow, B. J.; Walker, J. R.; Will, G. D.; Coote, M. L.; Bottle, S. E. J. Org. Chem. 2008, 73, 6763−6771. (41) Roy, L. E.; Batista, E. R.; Hay, P. J. Inorg. Chem. 2008, 47, 9228−9237. (42) Roy, L. E.; Jakubikova, E.; Guthrie, M. G.; Batista, E. R. J. Phys. Chem. A 2009, 113, 6745−6750. (43) Rulíšek, L. J. Phys. Chem. C 2013, 117, 16871−16877. (44) Castro, L.; Bühl, M. J. Chem. Theory Comput. 2014, 10, 243− 251. (45) Cheng, J.; Liu, X.; VandeVondele, J.; Sulpizi, M.; Sprik, M. Acc. Chem. Res. 2014, 47, 3522−3529. (46) Krewald, V.; Neese, F.; Pantazis, D. A. Phys. Chem. Chem. Phys. 2016, 18, 10739−10750. (47) CRC Handbook of Chemistry and Physics, 90th ed.; Lide, D. R., Ed.; CRC Press: Boca Raton, FL, 2009; p 2804. (48) NIST Chemistry WebBook, NIST Standard Reference Database Number 69; Linstrom, P. J., Mallard, W. G., Eds.; National Institute of Standards and Technology: Gaithersburg, MD, 2005. (49) Lin, J.; Lin, J. L.; Tzeng, W. B. Chem. Phys. Lett. 2003, 370, 44− 51. (50) Lin, J. L.; Tzeng, W. B. J. Chem. Phys. 2001, 115, 743−751. (51) Lin, J. L.; Lin, K. C.; Tzeng, W. B. J. Phys. Chem. A 2002, 106, 6462−6468. (52) Wu, R. H.; Lin, J. L.; Lin, J.; Yang, S. C.; Tzeng, W. B. J. Chem. Phys. 2003, 118, 4929−4937. (53) Johnstone, R. A. W.; Mellon, F. A. J. Chem. Soc., Faraday Trans. 2 1973, 69, 36−42. (54) Jonsson, M.; Wayner, D. D. M.; Lusztyk, J. J. Phys. Chem. 1996, 100, 17539−17543. (55) Boogaarts, M. G. H.; Hinnen, P. C.; Meijer, G. Chem. Phys. Lett. 1994, 223, 537−540. (56) Boogaarts, M. G. H.; von Helden, G.; Meijer, G. J. Chem. Phys. 1996, 105, 8556−8568. (57) Li, C.; Pradhan, M.; Tzeng, W. B. Chem. Phys. Lett. 2005, 411, 506−510. (58) Li, C.; Su, H.; Tzeng, W. B. Chem. Phys. Lett. 2005, 410, 99− 103. (59) Shirhatti, P. R.; Wategaonkar, S. Indian J. Phys. 2012, 86, 159− 164. (60) Vondrák, T.; Sato, S.-i.; Kimura, K. J. Phys. Chem. A 1997, 101, 2384−2389. (61) Choi, S.; Choi, K.-W.; Kim, S. K.; Chung, S.; Lee, S. J. Phys. Chem. A 2006, 110, 13183−13187. (62) Li, W. K.; Chiu, S. W.; Ma, Z. X.; Liao, C. L.; Ng, C. Y. J. Chem. Phys. 1993, 99, 8440−8444. (63) Vondrák, T.; Sato, S.-i.; Špirko, V.; Kimura, K. J. Phys. Chem. A 1997, 101, 8631−8638. (64) Sander, R. Atmos. Chem. Phys. Discuss. 2014, 14, 29615−30521. (65) Neese, F. WIREs Comput. Mol. Sci. 2012, 2, 73−78. (66) Tao, J.; Perdew, J. P.; Staroverov, V. N.; Scuseria, G. E. Phys. Rev. Lett. 2003, 91, 146401. (67) Weigend, F.; Ahlrichs, R. Phys. Chem. Chem. Phys. 2005, 7, 3297−3305. (68) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem. B 2009, 113, 6378−6396. (69) Feyereisen, M.; Fitzgerald, G.; Komornicki, A. Chem. Phys. Lett. 1993, 208, 359−363. (70) Eichkorn, K.; Weigend, F.; Treutler, O.; Ahlrichs, R. Theor. Chem. Acc. 1997, 97, 119−124. (71) Weigend, F.; Häser, M. Theor. Chem. Acc. 1997, 97, 331−340. (72) Weigend, F.; Häser, M.; Patzelt, H.; Ahlrichs, R. Chem. Phys. Lett. 1998, 294, 143−152. (73) Weigend, F.; Köhn, A.; Hättig, C. J. Chem. Phys. 2002, 116, 3175−3183. (74) Neese, F.; Hansen, A.; Wennmohs, F.; Grimme, S. Acc. Chem. Res. 2009, 42, 641−648. (75) Liakos, D. G.; Neese, F. J. Phys. Chem. A 2012, 116, 4801−4816. (76) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865−3868.

REFERENCES

(1) Marenich, A. V.; Ho, J.; Coote, M. L.; Cramer, C. J.; Truhlar, D. G. Phys. Chem. Chem. Phys. 2014, 16, 15068−15106. (2) Raghavachari, K.; Trucks, G. W.; Pople, J. A.; Head-Gordon, M. Chem. Phys. Lett. 1989, 157, 479−483. (3) Bartlett, R. J.; Musiał, M. Rev. Mod. Phys. 2007, 79, 291−352. (4) Xu, X.; Zhang, W.; Tang, M.; Truhlar, D. G. J. Chem. Theory Comput. 2015, 11, 2036−2052. (5) Feller, D.; Peterson, K. A.; Hill, J. G. J. Chem. Phys. 2011, 135, 044102. (6) Feller, D. J. Chem. Phys. 2013, 138, 074103. (7) Feller, D. J. Chem. Phys. 1992, 96, 6104−6114. (8) Nyden, M. R.; Petersson, G. A. J. Chem. Phys. 1981, 75, 1843− 1862. (9) Petersson, G. A.; Bennett, A.; Tensfeldt, T. G.; Al-Laham, M. A.; Shirley, W. A.; Mantzaris, J. J. Chem. Phys. 1988, 89, 2193−2218. (10) Helgaker, T.; Klopper, W.; Koch, H.; Noga, J. J. Chem. Phys. 1997, 106, 9639−9646. (11) Halkier, A.; Helgaker, T.; Jørgensen, P.; Klopper, W.; Koch, H.; Olsen, J.; Wilson, A. K. Chem. Phys. Lett. 1998, 286, 243−252. (12) Okoshi, M.; Atsumi, T.; Nakai, H. J. Comput. Chem. 2015, 36, 1075−1082. (13) Hättig, C.; Klopper, W.; Köhn, A.; Tew, D. P. Chem. Rev. 2012, 112, 4−74. (14) Kong, L.; Bischoff, F. A.; Valeev, E. F. Chem. Rev. 2012, 112, 75−107. (15) Kutzelnigg, W. Theor. Chim. Acta 1985, 68, 445−469. (16) May, A. J.; Manby, F. R. J. Chem. Phys. 2004, 121, 4479−4485. (17) Ten-no, S. Chem. Phys. Lett. 2004, 398, 56−61. (18) Neese, F.; Hansen, A.; Liakos, D. G. J. Chem. Phys. 2009, 131, 064103. (19) Neese, F.; Wennmohs, F.; Hansen, A. J. Chem. Phys. 2009, 130, 114108. (20) Hansen, A.; Liakos, D. G.; Neese, F. J. Chem. Phys. 2011, 135, 214102. (21) Liakos, D. G.; Hansen, A.; Neese, F. J. Chem. Theory Comput. 2011, 7, 76−87. (22) Sameera, W. M. C.; Pantazis, D. A. J. Chem. Theory Comput. 2012, 8, 2630−2645. (23) Liakos, D. G.; Neese, F. J. Chem. Theory Comput. 2015, 11, 4054−4063. (24) Cramer, C. J.; Truhlar, D. G. Chem. Rev. 1999, 99, 2161−2200. (25) Tomasi, J.; Mennucci, B.; Cammi, R. Chem. Rev. 2005, 105, 2999−3094. (26) Harris, R. C.; Pettitt, B. M. J. Chem. Theory Comput. 2015, 11, 4593−4600. (27) Orozco, M.; Luque, F. J. Chem. Rev. 2000, 100, 4187−4226. (28) Arumugam, K.; Becker, U. Minerals 2014, 4, 345−387. (29) Guerard, J. J.; Arey, J. S. J. Chem. Theory Comput. 2013, 9, 5046−5058. (30) Fu, Y.; Liu, L.; Yu, H.-Z.; Wang, Y.-M.; Guo, Q.-X. J. Am. Chem. Soc. 2005, 127, 7227−7234. (31) Ho, J. Phys. Chem. Chem. Phys. 2015, 17, 2859−2868. (32) Psciuk, B. T.; Lord, R. L.; Munk, B. H.; Schlegel, H. B. J. Chem. Theory Comput. 2012, 8, 5107−5123. (33) Baik, M.-H.; Friesner, R. A. J. Phys. Chem. A 2002, 106, 7407− 7412. (34) Thapa, B.; Schlegel, H. B. J. Phys. Chem. A 2015, 119, 5134− 5144. (35) Ghosh, D.; Roy, A.; Seidel, R.; Winter, B.; Bradforth, S.; Krylov, A. I. J. Phys. Chem. B 2012, 116, 7269−7280. (36) Sviatenko, L.; Isayev, O.; Gorb, L.; Hill, F.; Leszczynski, J. J. Comput. Chem. 2011, 32, 2195−2203. (37) Sviatenko, L. K.; Gorb, L.; Hill, F. C.; Leszczynski, J. J. Comput. Chem. 2013, 34, 1094−1100. (38) Hodgson, J. L.; Namazian, M.; Bottle, S. E.; Coote, M. L. J. Phys. Chem. A 2007, 111, 13595−13605. (39) Namazian, M.; Coote, M. L. J. Phys. Chem. A 2007, 111, 7227− 7232. 2283


Article

Journal of Chemical Theory and Computation (77) Zhao, Y.; Truhlar, D. G. J. Chem. Phys. 2006, 125, 194101. (78) Becke, A. D. J. Chem. Phys. 1993, 98, 5648−5652. (79) Adamo, C.; Barone, V. J. Chem. Phys. 1999, 110, 6158−6170. (80) Staroverov, V. N.; Scuseria, G. E.; Tao, J.; Perdew, J. P. J. Chem. Phys. 2003, 119, 12129−12137. (81) Zhao, Y.; Truhlar, D. G. Theor. Chem. Acc. 2008, 120, 215−241. (82) Zhao, Y.; Truhlar, D. G. Acc. Chem. Res. 2008, 41, 157−167. (83) Grimme, S. J. Chem. Phys. 2006, 124, 034108. (84) Schwabe, T.; Grimme, S. Phys. Chem. Chem. Phys. 2006, 8, 4398−4401. (85) Peterson, K. A.; Adler, T. B.; Werner, H.-J. J. Chem. Phys. 2008, 128, 084102. (86) Valeev, E. F. Chem. Phys. Lett. 2004, 395, 190−195. (87) Yousaf, K. E.; Peterson, K. A. J. Chem. Phys. 2008, 129, 184108. (88) Truhlar, D. G.; Cramer, C. J.; Lewis, A.; Bumpus, J. A. J. Chem. Educ. 2004, 81, 596. (89) Isse, A. A.; Gennaro, A. J. Phys. Chem. B 2010, 114, 7894−7899. (90) Klamt, A.; Schüürman, D. J. Chem. Soc., Perkin Trans. 2 1993, 799−805. (91) Byrd, E. F. C.; Sherrill, C. D.; Head-Gordon, M. J. Phys. Chem. A 2001, 105, 9736−9747. (92) Grimme, S. J. Chem. Phys. 2003, 118, 9095−9102. (93) Huntington, L. M. J.; Nooijen, M. J. Chem. Phys. 2010, 133, 184109. (94) Huntington, L. M. J.; Hansen, A.; Neese, F.; Nooijen, M. J. Chem. Phys. 2012, 136, 064101. (95) Schwabe, T. J. Comput. Chem. 2012, 33, 2067−2072. (96) Isegawa, M.; Peverati, R.; Truhlar, D. G. J. Chem. Phys. 2012, 137, 244104. (97) Isegawa, M.; Truhlar, D. G. J. Chem. Phys. 2013, 138, 134111. (98) Mendez-Hernandez, D. D.; Tarakeshwar, P.; Gust, D.; Moore, T. A.; Moore, A. L.; Mujica, V. J. Mol. Model. 2013, 19, 2845−2848. (99) van Leeuwen, R.; Baerends, E. J. Phys. Rev. A: At., Mol., Opt. Phys. 1994, 49, 2421−2431. (100) Parr, R. G.; Yang, W. Density-Functional Theory of Atoms and Molecules; Oxford University Press: Oxford, 1989; p 352. (101) Richard, R. M.; Marshall, M. S.; Dolgounitcheva, O.; Ortiz, J. V.; Brédas, J.-L.; Marom, N.; Sherrill, C. D. J. Chem. Theory Comput. 2016, 12, 595−604. (102) Gallandi, L.; Marom, N.; Rinke, P.; Körzdörfer, T. J. Chem. Theory Comput. 2016, 12, 605−614. (103) Winget, P.; Cramer, C. J.; Truhlar, D. G. Theor. Chem. Acc. 2004, 112, 217−227. (104) Le, H. T.; Flammang, R.; Gerbaux, P.; Bouchoux, G.; Nguyen, M. T. J. Phys. Chem. A 2001, 105, 11582−11592. (105) Klein, E.; Lukeš, V. Chem. Phys. 2006, 330, 515−525. (106) Winget, P.; Weber, E. J.; Cramer, C. J.; Truhlar, D. G. Phys. Chem. Chem. Phys. 2000, 2, 1231−1239. (107) Marten, B.; Kim, K.; Cortis, C.; Friesner, R. A.; Murphy, R. B.; Ringnalda, M. N.; Sitkoff, D.; Honig, B. J. Phys. Chem. 1996, 100, 11775−11788. (108) Li, J.; Fisher, C. L.; Chen, J. L.; Bashford, D.; Noodleman, L. Inorg. Chem. 1996, 35, 4694−4702. (109) Uudsemaa, M.; Tamm, T. J. Phys. Chem. A 2003, 107, 9997− 10003. (110) Gao, J. Acc. Chem. Res. 1996, 29, 298−305. (111) Isegawa, M.; Gao, J.; Truhlar, D. G. J. Chem. Phys. 2011, 135, 084107. (112) Gao, J. J. Phys. Chem. B 1997, 101, 657−663. (113) Watanabe, H.; Iwata, S. J. Chem. Phys. 1996, 105, 420−431. (114) Parthasarathi, R.; Subramanian, V.; Sathyamurthy, N. J. Phys. Chem. A 2005, 109, 843−850. (115) Riplinger, C.; Neese, F. J. Chem. Phys. 2013, 138, 034106.

2284


Ionization Energies and Aqueous Redox Potentials of Organic

Recommend Documents