Benchmarking Semiempirical Methods to Compute Electrochemical

Subscriber access provided by UNIV OF DURHAM

A: New Tools and Methods in Experiment and Theory

Benchmarking Semiempirical Methods to Compute Electrochemical Formal Potentials Rebecca L. M. Gieseking, Mark A. Ratner, and George C. Schatz J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.8b05143 • Publication Date (Web): 19 Jul 2018 Downloaded from http://pubs.acs.org on July 21, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Benchmarking Semiempirical Methods to Compute Electrochemical Formal Potentials

Rebecca L. M. Gieseking, Mark A. Ratner, George C. Schatz*

Department of Chemistry, Northwestern University 2145 Sheridan Road, Evanston, Illinois 60208, United States

* Corresponding author: [email protected]; +1(847)-491-5657

1 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract

Computational methods to predict and tune electrochemical redox potentials are important for the development of energy technologies. Here, we benchmark several semiempirical models to compute reduction potentials of organic molecules, comparing approaches based on (1) energy differences between the S0 and one-electron-reduced D0 states of the isolated molecules and (2) an orbital energy shift approach based on tuning the charge-transfer triplet energy of the Ag20molecule complex; the second model enables explicit modeling of electrode-molecule interactions. For molecules in solution, the two models yield nearly identical results. Both PM7 and PM6 predict formal potentials with only a slight loss of accuracy compared to standard density functional theory models, and the results are robust across several choices of geometries and implicit solvent models. PM6 and PM7 show dramatically improved accuracy over older semiempirical Hamiltonians (MNDO, AM1, PM3, and INDO/S). However, our recently developed INDO parameters model the electronic properties of our Ag20 model electrode much more accurately than does PM7. These results demonstrate the need for further development of semiempirical models to accurately model molecules on surfaces.


Page 2 of 31

Page 3 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1. Introduction

Electrochemical electron transfer processes are central to the functioning of many energy technologies, such as batteries,1–4 fuel cells,4,5 and electrocatalysis.6–8 Of particular interest is the ability to predict redox potentials for electrochemical reactions at interfaces and tune those potentials via structural modifications to either the molecule or the electrode surface. To this end, the development of computational techniques that can predict site-specific reduction potentials by modeling electron transfer events in the explicit presence of an electrode surface is important. However, the development of theoretical approaches to describe these processes is challenging due to the complexity of the experimental environment, including the near-continuum of electrode states, the applied potential, and the electrical double layer formed in an electrolyte solution. Significant simplifications must be made to develop a computationally tractable model.

There has been significant effort to develop computational approaches to model the formal potentials E° of molecules in solution.9,10 The most commonly used approaches have been based °, on the Born-Haber cycle, where the free energy of the redox reaction in solution ∆ is

computed based on the gas-phase free energy difference between the oxidized and reduced species and the solvation free energies of each species.9–14 However, recent work has shown that the vibrational contributions to the free energies are in many cases constant enough between the various species that they can be neglected.15 As a result, a DFT model based solely on computing the energies of the oxidized and reduced species in solution has been shown to predict E° with mean residual errors on the order of 70 mV for a broad range of organic π-conjugated molecules.16,17 Due to systematic errors in DFT, it is challenging to compute absolute potentials



but possible to have strong predictive power based on correlation between computational and experimental results.18 Although these approaches have been effective for predicting redox potentials of molecules in solution, these approaches are generally not straightforward to translate to molecules interacting with the electrode surface.

There has also been much work in recent years to develop computational models for electrochemical reactions with an electrode surface explicitly present; however, to date this work has focused largely on reactions involving small molecules such as H2,19–21 OH-,22–25 and CO2.26– 29

Although varying the electrode potential in principle requires dynamically modulating the

electrode charge,30–32 in practice most models involve either applying potential-specific thermodynamic corrections to computations at a single electrode charge33–35 or converting computations using several fixed charges to a series of effective electrode potentials.36–38 In these computations, the electrode is typically modeled as either a finite cluster or as an infinite periodic slab.23 Periodic slab models require maintaining an overall neutral unit cell, and several approaches to balance the electrode charge have been developed;36–38 these computations have been largely limited by computational expense to typical DFT functionals that are well-known to provide spuriously low charge-transfer energies,39,40 which limits their reliability for these applications. In contrast, cluster models have the advantages of lower computational cost and the ability to model a system with overall non-zero charge, but electrode finite size effects may limit accuracy.23,41–43 These models are also typically restricted to unit charges, which limits voltage changes to large steps; for an Ag20 cluster electrode, addition or removal of a single electron shifts the potential by nearly 3 V,44 much larger than the differences on the order of tens to hundreds of mV that are often relevant when comparing molecules or surface sites.


Page 4 of 31

Page 5 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


To address these limitations, we have recently introduced an approach to compute the formal potentials of molecules at electrode surfaces based on the orbital energy shift approximation (OESA).44 This model is based on using a semiempirical Hamiltonian such as INDO, which has been shown to yield more accurate charge-transfer energies than typical DFT approaches.45,46 The electrode potential is tuned by applying a fixed shift to the parameters corresponding to the metal atomic orbital energies.24,47 We have used a tetrahedral Ag20 cluster as a model electrode; the faces of the this cluster are equivalent to the Ag(111) crystal face commonly used as an electrode experimentally, and the edge and vertex sites enable investigation of site-specific effects. We have focused on systems where both the Ag20 cluster model for the electrode and the oxidized form of the molecule of interest are closed-shell, so that the most stable [Ag20molecule] charge-localized state is a closed-shell singlet and the most stable [Ag20+-molecule-] charge-transfer state is a triplet. Thus, the site-specific formal potential can be determined by applying a shift to the Ag orbital energies such that the charge-localized singlet and chargetransfer triplet are equal in energy. We have shown that this approach predicts significant sitespecific variation in molecular formal potentials,44 highlighting the importance of explicitly modeling the electrode surface. However, to date the OESA model has not been benchmarked to experimental data.

Here, we benchmark several semiempirical Hamiltonians for computing reduction potentials of molecules in acetonitrile (ACN), considering both our OESA approach and an energy difference approach (EDA) based on typical DFT methods for molecules in solution.16,17 We show that for molecules in solution, the EDA and OESA approaches yield nearly identical formal potentials.



The PM7 Hamiltonian predicts formal potentials with only a slight loss of accuracy compared to DFT, and the results are robust across several choices of geometries and implicit solvent models. The semiempirical models reduce the computational cost of these computations by 1-3 orders of magnitude relative to DFT with minimal loss of accuracy and avoid the problems of selfinteraction error that plague charge-transfer energies in typical DFT functionals. Our results also highlight the dramatic improvement of semiempirical methods of the 2000s and later (PM6 and PM7) over those of the 1980s and earlier (MNDO, AM1, PM3, and INDO/S). However, our recently developed INDO parameters48,49 model the electronic properties of the Ag20 model electrode much more accurately than does PM7. These results demonstrate the need for further development and parameterization of semiempirical models to improve the predictive ability of these methods for the properties of molecules on surfaces.

2. Methodology

Benchmarking of computational methods was performed using a set of 74 molecules for which experimental formal potentials for reduction in ACN have been experimentally measured relative to the saturated calomel electrode (SCE), as compiled by Lynch et al.16,17 This set includes 13 cyanoaromatic molecules, 13 quinones, 9 N-methyl heteroaromatic cations, 9 flexible πconjugated molecules, 12 polycyclic aromatic hydrocarbons (PAHs), and 18 heterocyclic amines, as shown in the SI. As in previous DFT studies,16,17 we focus on reduction of these molecules from the S0 states to their one-electron-reduced D0 states.


Page 6 of 31

Page 7 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The geometries of the isolated tetrahedral Ag20 cluster and of the isolated molecules in both their S0 and one-electron-reduced D0 forms were optimized using density functional theory (DFT) with the BP86 functional50,51 and a double-zeta (DZ) Slater basis set; this approach has been shown to perform well for the geometries of Ag clusters.52,53 A frozen-core approximation was used for Ag atoms for all electrons in the 4p and lower orbitals (DZ.4p). All optimizations were performed using the COSMO model of implicit solvation for ACN (ε = 37.5).54,55 The solvent radius was set to 1.0 Å, and the screening energy was scaled by a factor of = − 1⁄ + 0.5 relative to a perfect conductor for consistency with the original COSMO algorithm.55 Scalar relativistic effects were incorporated using the zero-order regular approximation (ZORA).56 All DFT calculations were performed using the Amsterdam Density Functional (ADF) 2016 program.57,58

Semiempirical calculations to compute the formal potentials were computed using several approaches. All semiempirical calculations were performed using a home-built code incorporating parts of the INDO/CI code from Jeffrey Reimers59 into MOPAC2016.60,61 To benchmark the energy difference approach (EDA), we optimized the S0 and one-electronreduced D0 states of each molecule using a series of semiempirical Hamiltonians, including the PM7,62 PM6,63 PM3,64,65 AM1,66 and MNDO67 Hamiltonians. The open-shell systems were treated using a UHF scheme. The COSMO solvent model was used, with van der Waals radii for all atoms,68,69 the solvent radius set to 1.3 Å, and the screening energy scaled by a factor of = − 1⁄ + 0.5 relative to a perfect conductor54 as is the default in MOPAC2016.



For the PM7 Hamiltonian, further benchmarking of the EDA was performed. Single-point energies of the S0 states were computed using the S0 DFT-optimized geometries, and single-point energies of the D0 states were computed using both the S0 and D0 DFT-optimized geometries. The PM7 S0 and D0 state geometries were also optimized using an alternate set of COSMO parameters using Bondi radii70 scaled by a factor of 1.2 and the solvent radius set to 1.0 Å.

The formal potentials were also computed using our recently developed Orbital Energy Shift Approximation (OESA)44 using the PM7 and INDO Hamiltonians. The INDO calculations used our recently developed parameters for Ag48,49,71 and Zerner’s INDO/S parameters72 for all other atoms. Using the S0 DFT-optimized geometries, each molecule was oriented roughly parallel to one of the triangular faces of the Ag20 cluster, with its closest atom 8 Å from the Ag20 face as shown in Figure 1. Our previous results show that this distance is sufficient for the distance dependence of the computed formal potential to be negligible in implicit solvent.44 The S0 state of the Ag20-molecule complex was first computed, and then the triplet states were computed as excited states using a single configuration interaction (SCI) approach. Within the voltage range of interest, the first triplet state is a charge-transfer state from Ag20 to the molecule, which is effectively the one-electron-reduced D0 state of the molecule weakly interacting with the D0 state of [Ag20]+. As we are interested in the triplet states with a fully relaxed configuration of the solvent,44 COSMO corrections to the triplet states73,74 were computed using a dielectric constant equal to that used for the ground states (ε = 37.5). Within the INDO model, all possible oneelectron excitations of triplet multiplicity were generated, and the lowest 1000 configurations were included in the CI matrix; this matrix was then diagonalized to obtain the lowest 100 triplet states. Within the PM7 model, the CI matrix was constructed using all one-electron excitations


Page 8 of 31

Page 9 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


within the first 12 HOMOs and first 8 LUMOs. As only the first triplet state is of interest here, the computed properties are fairly insensitive to the size of the CI matrix.

Figure 1. Depiction of the geometry of the Ag20-molecule complex selected for the OESA approach to model formal potentials.

To compute the electrochemical properties, the PM7 and INDO parameters corresponding to the Ag s, p, and d orbital energies were shifted by a value corresponding to the applied potential. To determine the OESA computed formal potential, we find the shift in Ag orbital energies required to obtain equal energies for the first charge-localized singlet and charge-transfer triplet states. The potential was varied using linear interpolation based on the singlet–triplet energy gap until the singlet and triplet state energies converged to within 0.001 eV.

3. Results and Discussion

3.1. PM7 Hamiltonian



Page 10 of 31

To evaluate the accuracy of semiempirical methods for computing the formal potentials E° for reduction of organic molecules, we consider a set of 74 molecules in acetonitrile (ACN) that has been previously used in DFT benchmarking studies (see SI for molecular structures).16,17 This reference set is suitable for this study because of its focus on data in solution, which is much more sparse in typical benchmarking data sets.75 We focus initially on the PM7 Hamiltonian. We start by using an energy difference approach (EDA), examining the correlation of the energy difference between the S0 and one-electron-reduced D0 states with the experimental E°. Similar to previous DFT results,17 the computed D0-S0 energy difference at the PM7 level with COSMO implicit solvation is strongly correlated with the experimental E°, with the correlation coefficient R2 > 0.95 when considering geometries for the S0 and D0 states optimized at the PM7 level as shown in Figure 2a. The formal potential E° of an unknown molecule can thus be computed as: , = −1.327 #, − 7.040

1

where #, is the energy difference between the S0 and D0 states in eV. The mean absolute deviation (MAD) of the residual errors is < 100 mV, representing a loss of accuracy of ~ 30 mV from the previously reported DFT results at the B3LYP/6-311G(d,p) level (residuals < 70 mV) as shown in

Table 1.17 We note that our PM7 geometry optimizations are approximately three orders of magnitude faster than our DFT optimizations using a double-ζ basis set, and the time savings are more significant when comparing to the triple-ζ basis set data in Ref. 17. This substantial reduction in computational cost enables rapid screening of far more molecules with little loss of accuracy.


Page 11 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 1. Correlations of computed D0-S0 energy differences with experimental reduction potentials. The slope and intercept reported here define the predicted E° value based on the computed energy difference. Method

Geometries

B3LYP/6311G(d,p)17 PM7 PM7

B3LYP/6311G(d,p) PM7 PM7 BP6/DZ (neutral) BP6/DZ (neutral/anion)

PM7 PM7

Solvent radii

R2

MAD residual (V)

Slope (V/eV)

Intercept (V)

UAKS

0.9889

0.063

-0.838

-3.931

VDW Bondi

0.9581 0.9522

0.093 0.100

-1.327 -1.328

-7.040 -7.003

VDW

0.9540

0.097

-1.242

-6.676

VDW

0.9520

0.117

-1.322

-6.839



Figure 2. Correlation between D0-S0 energy differences computed at the PM7 level with COSMO solvation and experimental formal potentials. Panels a-d compare several choices of COSMO parameters and molecular geometries, as detailed in each panel.

Overall, this correlation is quite robust across all six chemical classes of molecules in this reference set. The largest residual errors of -500 and -370 mV occur for 1,4-cyano-2,3,5,6fluorobenzene (3 as shown in the SI) and 1,4-cyano-2,3,5,6-chlorobenzene (4), respectively (E° of -1.10 and -0.95 V, respectively), which may be related to limitations of the PM7 parameters for halides. The largest positive residual is 320 mV for benzophenone (38; E° = -1.86 V), which is the only molecule in this set with a large out-of-plane twist between two identical conjugated moieties. All other molecules have residuals < 300 mV, suggesting fairly consistent performance of PM7 across a range of chemical functionalities in conjugated organic molecules. The small increase in residuals relative to the previous DFT results suggests that the PM7 parameters are slightly less transferable between classes of molecules, so care should be taken to ensure that the predictions are benchmarked to a test set involving reasonably similar molecules.

This strong correlation is robust across small methodological variations in the choice of COSMO parameters and molecular geometries. The results are nearly identical when constructing the solvent-accessible surface using either the Bondi radii or the van der Waals atomic radii (Figure 2a,b; see also

Table 1); we focus on results using the van der Waals radii for the remainder of this paper since these radii lead to a slightly larger R2 value. We also compare results using the PM7-optimized geometries of the S0 and D0 states to single-point PM7 computations using DFT-optimized geometries at the BP86/DZ level; this functional was chosen for its good performance for metal 12 ACS Paragon Plus Environment

Page 12 of 31

Page 13 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


clusters, which will become relevant later when we consider the OESA approach. Two choices of DFT-optimized geometries were chosen: (1) DFT geometries of the S0 state were used to compute both PM7 S0 and D0 energies (Figure 2c), and (2) DFT geometries of the S0 and D0 states, respectively, were used for the PM7 S0 and D0 energies (Figure 2d). In both cases, the R2 values remain > 0.95; the MAD remains at 100 mV when only DFT S0 geometries are used but increases slightly to 117 mV when the DFT S0 and D0 geometries are used. These results show that PM7 has adequate accuracy to predict the formal potentials of molecules and that minor methodological changes do not significantly change this predictive ability.

Since the correlation of the B3LYP computed formal potentials with experimental formal potentials can be slightly improved to a MAD of < 60 mV by a simpler computational approach considering only the LUMO energies of the S0 species,76 we perform an analogous analysis for the PM7 results (Figure 3). Although the LUMO energies of the S0 species at the PM7 level are reasonably well correlated with the experimental formal potentials, the MAD increases to 155 mV. Thus, in contrast with the DFT results, the S0-D0 energy difference is a more accurate approach to predict formal potentials at the PM7 level than is the LUMO energy.



Figure 3. Correlation of PM7 computed LUMO energies with experimental formal potentials.

We now turn to our OESA approach,44 which allows us to model formal potentials in the presence of an electrode, modeled as an Ag20 cluster. Since we are interested here in evaluating the accuracy of this approach for molecules in solution, we place the molecules with their nearest atom 8 Å from the electrode (Figure 1), using DFT-optimized S0 geometries of the isolated Ag20 cluster and the molecule. Our previous results show that this distance is large enough to fully solvate the molecules such that the computed E° is nearly distance-independent.44 Although here we focus on benchmarking this approach for well-solvated molecules, we have shown that this excited-state-based scheme enables computation of site-specific formal potentials for adsorbed molecules.44 This scheme also explicitly includes the Coulomb interaction between the electron transferred to the molecule and the hole remaining on the electrode, which is neglected in the previous models with no explicit electrode. Even though we explicitly include an Ag20 cluster, these computations are still 1-2 orders of magnitude faster than our DFT computations for the isolated molecules; these calculations can also be extended to include explicit interactions of the molecule with the electrode surface, which would be much more computationally expensive using a DFT-based approach.

The correlation between the OESA orbital shifts and the experimental E° values is nearly as good as for the EDA model described previously, with R2 of 0.9388 and a MAD residual of 120 mV as shown in Figure 4a, showing that both models have comparable ability to predict E°. This similar predictive ability is due to the extremely strong internal consistency between the OESA and EDA results. As shown in Figure 4b, the computed OESA orbital energy shifts are nearly 14 ACS Paragon Plus Environment

Page 14 of 31

Page 15 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


perfectly correlated with the EDA energy differences, with a slope of 1.0048 and R2 of 0.985. The formal potential E° of an unknown molecule can thus be computed as: ,%& = −1.237 ',%& + 0.586

2

where ',%& is the orbital energy shift required to obtain equal energies of the singlet and triplet states. Despite the apparent differences between the two computational approaches, this strong correlation is unsurprising given the details of the models. The OESA formal potential is formally based on computing a charge-transfer energy as the T1 excited state of the S0 ground state. However, in our model, the Ag20 and molecular moieties are sufficiently far apart that the interactions are negligible. Thus, the S0 energy can be approximated as a sum of the energies of the S0 states of the isolated Ag20 and molecule with a small interaction term; the T1 energy can similarly be approximated as a sum of the energies of the D0 states of the isolated Ag20+ and oneelectron-reduced molecule with a small interaction term. The residual differences of < 60 mV are related to the interactions between the Ag20 and molecular moieties.



Figure 4. Correlation of PM7 computed OESA orbital energy shifts with experimental formal potentials, using (a) Ag20 and (c) Ag2 model electrodes. (b) Correlation of PM7 computed OESA orbital energy shifts (Ag20 electrode) with the EDA D0-S0 energy differences (using DFT S0 geometries).


Page 16 of 31

Page 17 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Although here we focus on Ag20, these results can be generalized to other electrode structures. As shown in the SI, the orbital shift computed for each molecule using the smaller Ag2 as an electrode is larger than the value computed using Ag20 by 1.703 eV as shown in Figure 4c; the equation to predict E° is shown in Table 2. This difference corresponds to the difference between the ionization potentials of the clusters. Thus, predicting the formal potential of a molecule near an arbitrary electrode can be done by performing a single benchmarking computation to determine the energetic offset and does not require re-benchmarking the full series of molecules with a new electrode.

Table 2. Correlations of computed orbital energy shifts with experimental reduction potentials. The slope and intercept reported here define the predicted E° value based on the computed energy difference. Method

Electrode

R2

PM7 PM7 INDO (all molecules) INDO (planar aromatic only)

Ag20 Ag2

0.9388 0.9405

MAD residual (V) 0.125 0.124

Ag20

0.8100

Ag20

0.8201

Slope (V/eV)

Intercept (V)

-1.237 -1.237

0.586 2.694

0.307

-1.076

1.312

0.133

-1.629

3.091

Interestingly, the slope of the predicted E° value with computed orbital shift is -1.237, significantly different from the slope of -1 assumed in our previous work developing the OESA model44 and in similar earlier models.24,47,77 This slope is nearly identical to that seen for PM7 in the EDA model but significantly larger than the slopes of -0.8 to -0.9 seen in previous DFT models.16,17 Since the current benchmarking is limited to molecules in ACN, this suggests the need for benchmarking to high-quality experimental data in other solvents to learn more about general trends and enable prediction of formal potentials in a broader range of environments. 17 ACS Paragon Plus Environment


In addition to showing that the PM7 Hamiltonian is able to predict E° values with little loss of accuracy relative to DFT-based predictions, these results also show that the OESA model has predictive ability equivalent to that of the EDA model for well-solvated molecules. The OESA model also provides the advantage of explicitly including an electrode. This will enable us to directly model molecule-electrode interactions and evaluate their effect on the molecular formal potential,44 an important capability for understanding single-molecule electrochemical experiments.

3.2. Comparison to INDO Hamiltonian

We now turn to a comparison of PM7 with the INDO Hamiltonian, which we have previously used for electrochemical computations.44,48 We have recently re-parameterized INDO to accurately model the excited-state properties of Ag clusters, showing that the INDO electronic structure of the Ag20 cluster is comparable to TDDFT results.48,49 The INDO/S parameters used for all other atoms72 were developed specifically for excited-state properties and are known to have limited accuracy for ground-state energies. Because of technical limitations of our code for open-shell systems, our discussion here focuses solely on the OESA model.

The INDO Hamiltonian shows some correlation between the computed orbital shift and the experimental E°, but with significantly poorer correlation than for PM7 (Figure 5); R2 decreases to only 0.8100 and the MAD of the residuals increases to 307 mV as shown in Table 2. This poorer correlation is unsurprising given the known limitations of INDO/S for ground-state


Page 18 of 31

Page 19 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


energies. As described in the previous section, even though the orbital shift is computed as an excited-state property, the computation effectively probes the relative energies of the S0 and D0 states of the molecules.

Figure 5. Correlation of INDO computed OESA orbital energy shifts with experimental formal potentials, showing results for (a) all molecules and (b) planar aromatic molecules.

In contrast with the PM7 results, the INDO model shows significant deviations among the predicted formal potentials for molecules with different chemical functionalities. In particular, most of the quinones have significant negative deviations from the overall trend line; in addition,



the flexible-π molecules show significant scatter and large deviations from the overall trend line. Interestingly, these two sets of molecules are the only molecules in this reference set that are not planar aromatic molecules. To limit the effect of variations in chemical functionality, we show the correlation for only planar aromatic molecules in Figure 5b. With this more limited data set, the MAD residual drops from 307 mV to 132 mV, showing a substantial improvement in correlation when comparing molecules that are chemically more similar. Even when comparing chemically similar molecules, INDO has some predictive power but is somewhat more limited than PM7 or DFT.

Since we are interested in modeling the effects of molecule-electrode interactions on formal potentials, it is important to accurately model the electronic structures of both the molecule and the electrode. Thus, we compare the electronic structure of the Ag20 cluster at the INDO and PM7 levels. We note that the semiempirical methods that will be discussed in the next section either do not have parameters for Ag, or (in the case of PM6) do not give SCF convergence for Ag20-molecule complexes in implicit solvent.

As has been previously shown,49,52 at the BP86 and INDO levels, the molecular orbitals of the Ag20 cluster comprise a fully occupied d band and a partially occupied sp band, commonly referred to as the valence and conduction bands; some mixing between these bands also occurs. In contrast, at the PM7 level, there are discrete d, s, and p bands with little to no mixing among bands. This lack of mixing dramatically increases the excited-state energies and the energy of the first absorption peak. Both BP86 and INDO/SCI yield first excited states for the Ag20 complex around 1.8 eV; the low-energy excited states are weakly absorbing, and strong absorption peaks


Page 20 of 31

Page 21 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


are seen at 3.58 eV and 3.45 eV, respectively (Figure 6). The BP86 absorption peaks have been shown to red-shift with increasing cluster size and extrapolate to values consistent with classical electrodynamics and experimental results.52 In contrast, at the PM7/SCI level, the first excited state of the Ag20 complex is at 5.24 eV, and the first strong absorption peak is around 8 eV. The failure to describe the electronic structure of Ag20 is unsurprising given the limitations of the parameterization of PM7. Ag is specifically mentioned as an example where the benchmark data for PM7 was sparse enough that the parameters are not expected to be reliable, and the p orbital energies in metals were artificially raised in PM7 to improve performance for solids.62

Figure 6. Excited-state energies of Ag20 computed at several levels of theory.

These results point to the need for further development of semiempirical methods to understand electrochemistry on surfaces at a much lower computational cost than DFT and without the concerns of self-interaction error. The PM7 Hamiltonian is sufficiently accurate for ground-state properties of organic molecules but is limited in its ability to describe the electrode surfaces, despite developments that improve its performance for solids over PM6.62 In contrast, INDO is 21 ACS Paragon Plus Environment


well-known to perform reasonably well for excited-state properties of organic molecules,78 and we have demonstrated its ability to model the electronic and excited state properties of metals;49 however, the limitations of the INDO/S parameters for ground-state energies limit its accuracy for electrochemical properties.

3.3. Comparison to other semiempirical models

Since semiempirical methods have a long history, we now turn to a comparison of the PM7 results to those using the older semiempirical Hamiltonians PM6 (introduced 2007),63 PM3 (1989),64 AM1 (1985),66 and MNDO (1977),79 focusing solely on the EDA model. The energy differences computed at the PM6 level are nearly identical to those at the PM7 level (Figure 7). This is unsurprising, given the relatively minor modifications in the parameterization procedure and reference data for organic elements between these two models.62,63 The MAD residuals are ~100 mV for both methods (Table 3), and there are only a handful of molecules where the PM6 and PM7 energy differences differ by more than 150 mV.


Page 22 of 31

Page 23 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 7. Correlation between D0-S0 energy differences computed at the (a) PM6, PM7, (b) PM3, MNDO, and AM1 levels and experimental formal potentials.

Table 3. Correlations of computed D0-S0 energy differences with experimental reduction potentials. The slope and intercept reported here define the predicted E° value based on the computed energy difference. Method

R2

PM7 PM6 PM3 MNDO AM1

0.9581 0.9587 0.6566 0.7068 0.7657

MAD residual Slope (V/eV) (V) 0.093 -1.3282 0.089 -1.2741 0.220 -1.7549 0.279 -1.5637 0.198 -1.4811


Intercept (V) -7.0026 -6.8228 -8.3892 -7.7769 -7.2960


In contrast, the older semiempirical models MNDO, AM1, and PM3 clearly show much poorer ability to predict formal potentials. Similar to the INDO results, these methods show large deviations in the trend lines between different chemical groups, in particular significant differences between the quinones and planar aromatic molecules. Even if the quinones are excluded, all three of these methods have residual errors greater than 1 eV for some molecules. This highlights the significant advances in semiempirical methods from the methods of the 1970s and 1980s with relatively poor predictive power to the methods of the 2000s and later (PM6 and PM7) which are competitive with DFT in their accuracy at substantially lower computational cost.

4. Conclusions

Computational modeling is a key tool to understand electrochemical reactions in complex environments; however, it is important to benchmark the accuracy of these computational methods to enable direct comparison to experimental results and prediction of electrochemical properties of novel systems. We have benchmarked the accuracy of the EDA model, based on energy differences between two charge states of the isolated molecule, and the OESA model, where the bias is tuned by changing energy levels of a small model electrode, for semiempirical modeling of electrochemical formal potentials.

The PM7 results are fairly accurate and consistent across the chemical functionalities in our data set; the MAD of residual errors is

Benchmarking Semiempirical Methods to Compute Electrochemical

Recommend Documents