Redox Levels through Constant Fermi-Level ab Initio Molecular

Article pubs.acs.org/JCTC

Redox Levels through Constant Fermi-Level ab Initio Molecular Dynamics Assil Bouzid* and Alfredo Pasquarello Chaire de Simulation à l’Echelle Atomique (CSEA), Ecole Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland ABSTRACT: We develop a method to determine redox levels of half reactions through the use of ab initio molecular dynamics evolving at constant Fermi energy. This scheme models the effect of an electrode by controlling the charge transfer between the single-particle energy levels of the system and an electron reservoir set at a given potential during the dynamics. Like the thermodynamic integration method, our scheme does not require a priori knowledge of the products of the reaction, which can simply be obtained by driving the reaction through the variation of the Fermi level. The simulations are performed subject to periodic boundary conditions in the absence of any counterelectrode. We extract the redox level from the evolution of the Kohn−Sham energies upon charging (or discharging) the system. Using Janak’s theorem and assuming a quadratic evolution of the Kohn−Sham energy of the highest occupied state upon charging, we demonstrate that our scheme for the determination of redox levels is equivalent to the thermodynamic integration method. The approach is illustrated for the redox couples Fe2+/Fe3+, HO•2 /HO−2 , and MnO−4 /MnO−2 4 in aqueous solution and yields redox potentials differing by less than 0.1 eV from respective ones achieved with the thermodynamic integration method. was described by a continuum model.19,20 These methods are widely used and are available in various quantum-chemistry codes. The second approach consists in treating explicitly the solvent molecules around the solute through density functional based molecular dynamics (DFT-MD). Albeit computationally expensive, this approach provides an accurate description of the solute structure and accounts for the structural evolution of the solvent during the reaction. In the following, we mainly focus on the DFT-MD level of theory. Within the DFT-MD approach, the theoretical determination of redox levels and their related electrochemical properties have gone through various steps of improvement. In an early paper, Tavernelli et al. modeled redox reactions by allowing the exchange of electrons between the system and an electron reservoir, which mimics the presence of an electrode in the solution.21 In their formulation, the potential energy surfaces (PES) of the oxidized and the reduced species are superimposed in a grand-canonical ensemble. Both PESs correspond to an integer number of electrons, and their weight in the superposition is governed by their vertical ionization energies. While this technique does not provide a clear connection to the redox level, it nevertheless achieves a successful chemical description of the charge transition in several aqueous redox reactions.22−26 Next, the thermodynamic integration (TI) method provides a clear simulation procedure to access the redox level. In this method, the Hamiltonians of the reactants and the products are continuously connected through a series of fictitious Hamil-

1. INTRODUCTION Redox potentials are relevant in a variety of research areas, which include fuel cells, solar-energy conversion, energy storage technology, and biological systems. The control and design of redox-based catalysts and the development of redox-based devices requires a clear understanding of the underlying oxidoreduction mechanism as well as the related electrochemical properties. Therefore, methodologies allowing for the determination and prediction of such mechanisms and properties are of invaluable importance. In the last decades, several computational methods have been developed to study the thermodynamical properties of redox half reactions. In particular, most of the methods developed to determine redox potentials are based on the thermodynamic Born−Haber cycle linking the gas and liquid phases. Following this cycle, the redox potentials can be obtained by evaluating free energy differences associated with the reaction. The accurate determination of redox potentials in solution depends crucially on the description of solvation effects. Hence, several methods have been developed to model the solvent, either in an implicit or explicit fashion. In practice, two approaches are possible. The first approach consists in replacing the solvent by a homogeneous dielectric continuum, as for example in the polarizable continuum model (PCM).1 In such approaches, it is assumed that the dynamical evolution of the solvent does not affect the reaction in a relevant way. Despite its success in reproducing the redox potentials of many systems,2−15 this technique proved unsatisfactory when applied to transition metal complexes.16−18 To overcome this limitation, the first two hydration shells around the complex were treated explicitly while the remaining part of the solvent © XXXX American Chemical Society

Received: December 21, 2016 Published: March 2, 2017 A

DOI: 10.1021/acs.jctc.6b01232 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation tonians /η depending on a Kirkwood coupling parameter η ranging between 0 and 1.27 The path between the reactants and the products is sampled through molecular dynamics simulations governed by the Hamiltonian /η of the fictious intermediate systems. The free-energy difference between the reactants and products is obtained as a thermodynamic integral of the vertical energy gap over η.28−30 In this context, the development of the standard hydrogen electrode (SHE) by Sprik and co-workers has been instrumental in facilitating the comparison between calculated and measured redox potentials.31−36 Because of their high computational cost, DFT-MD simulations have mostly been performed through the use of density functionals based on the generalized gradient approximation (GGA) for the treatment of the exchange and correlation energy. Very recently, a higher accuracy has been achieved by combining these methods with hybrid functionals33,36,37 or with the random phase approximation.38 Recently, Ambrosio et al.39 made available a formulation for defining redox levels in aqueous solutions in analogy to the definition of charge transition levels of defects in solids. In this work, the authors introduced an electron chemical potential when considering redox half-reactions and calculated the associated free energy differences through the TI method. The redox level is then obtained as the electron chemical potential at which the Gibbs free energies of formation of the oxidized and reduced systems are equal. Subsequently, this formulation has also been applied to determine structural and electrochemical properties of aqueous species at the hybrid functional level of theory.37 In experiments, redox levels are measured through cyclic voltammetry. This technique requires the presence of an explicit electrode in the solution. Upon application of a voltage, the induced current is measured. The redox potential corresponds to the voltage at which a rapid increase of the electric current is observed. This increase results from the release of free electrons into the solution upon occurrence of oxidation.40 In spite of the success of computational methods in computing the electrochemical properties of redox pairs, the simulation of a redox reaction at the surface of an explicit electrode and the consequent structural modifications are still out of reach. As discussed above, the redox level is obtained through a thermodynamic cycle in which the various transformations are studied separately. In this case, the reaction is simulated in the absence of an electrode/electrolyte interface and by controlling the number of electrons rather than the voltage. A major step toward modeling the presence of an explicit electrode has been achieved by Lozovoi et al.41 These authors performed simulations at constant bias potential and considered an electric field at the metal/water interface. The slab was allowed to exchange electrons with a reference electrode under the constraint of maintaining a constant preset chemical potential. This supercell approach made possible the simulation of charged slabs as in electrochemical set-ups. Charge neutrality was preserved through the introduction of counterions in the supercell. Albeit useful, Lozovoi’s method is limited by the strong fluctuations of the electron density which make the electronic structure calculations difficult.41,42 More sophisticated techniques were introduced by Otani et al. 43 and Jinnouchi et al. 44 to overcome the charge compensation problem. However, they could study the dynamics at the solvent/electrode interface only at fixed excess charge.45,46 In a later paper, Bonnet et al. developed a scheme

to perform molecular dynamics simulations at constant electrode potential42 in which the reactions at the metal/ electrolyte interface could be simulated in a nonrepeated supercell. Similarly to the way the Nosé−Hoover thermostat controls the temperature of a system,47−50 this technique treats the electronic charge as a fictitious dynamical variable and allows for charge exchanges with an external potentiostat at a preset potential. The equations of motion of the electronic charge and of the atomic positions are derived from the potential energy of the extended system, which includes the atoms and the potentiostat. The sequential treatment of the electronic structure and of the charge dynamics ensures a robust and efficient algorithm. The authors adopt the effective screening medium (ESM) technique43 to compensate for the excess charge in the supercell. While this compensation technique allows one to gain insight into the reaction mechanism at the metal/water interface, it also affects the energetics, and the redox potential cannot trivially be accessed. In this work, we develop a scheme for modeling redox halfreactions through the use of first-principles molecular dynamics at constant Fermi level. Following ref 42, we connect the system to an external reservoir of electrons acting like a potentiostat at a preset chemical potential. The present constant Fermi energy scheme allows us to control the charge transfer between the electrode potential and the single-particle energy levels of the system. The charge dynamics are governed by a set of fictitious Newton-like equations of motion. In our simulations, the investigated systems evolve under periodic boundary conditions in the absence of any counterelectrode altering the electronic properties of the system. In this work, we first show that we can drive redox half-reactions by varying the Fermi level, thereby revealing reaction mechanisms without any prior knowledge of the reaction products, just like with thermodynamic integration methods. Next, on the basis of Janak’s theorem,51 we derive a scheme for determining redox potentials from the evolution of single-particle energy levels upon charging (or discharging) the system. The developed scheme is here illustrated through the study of redox halfreactions of various aqueous species within a semilocal density functional scheme augmented to include van der Waals interactions. For the considered systems, the calculated redox levels are found to agree closely with those obtained via the thermodynamic integration method. The agreement with experiment is consistent with the accuracy of the adopted level of theory.

2. CONSTANT FERMI-LEVEL SIMULATION 2.1. Formulation. In this section, we first review the dynamical equations for achieving constant Fermi-level simulations, as introduced by Bonnet et al.42 We consider a system of particles described by a set of atomic positions ri and a total electronic charge Ne. To control the Fermi level, this system is connected to an external potentiostat at a fixed Fermi level ϵF̅ acting like an electron reservoir or a fictitious external electrode. The electronic charge Ne is considered as a dynamical variable with inertia Me, and the extended system is driven by the grand canonical potential:42 Ω = Etot(ri, Ne) − NeϵF̅ , where Etot(ri, Ne) is the potential energy in the absence of the electron reservoir and NeϵF̅ the corresponding energy of Ne electrons in the reservoir. According to this definition, the forces acting on the atoms and the electronic charges are given by42 B


Article

Journal of Chemical Theory and Computation Fi = −

∂Etot(ri , Ne) and Fe = −(ϵF − ϵF̅ ) ∂ri

occupations. In our simulations, this choice of Gaussian smearing ensures that only the occupation f of the singleparticle state g falling within the band gap carries a fractional value and does not affect the occupation of the other states. Because our simulated systems are generally charged, a mechanism to ensure charge neutrality should be adopted when using periodic boundary conditions. Several compensation mechanisms have been proposed for constant potential simulations.21,43,58,59 In our setup, we did not encounter any instability in resorting to a uniform neutralizing background, as the charge dynamics and the electronic structure calculation are carried out sequentially in the constant Fermi-level scheme. At each MD step, the electronic charge is first computed via eq 2. Then, the electronic structure is minimized at this fixed charge. Next, the atomic positions and the electronic charge are updated via the forces Fi and Fe defined in section 2.1. The varying electronic charge within the simulation cell involves electrostatic finite size effects that need to be accounted for. In this work, we adopted the finite-size corrections introduced by Freysoldt, Neugebauer, and Van de Walle (FNV).60,61 In the present work, we used the dielectric constant of liquid water at 300 K (78.3, ref 62). The resulting corrections are negligible in our system for the singly charged state (0.03 eV) but amount to 0.10 and 0.23 eV for +2 and +3 charge states, respectively. The use of the dielectric constant at 350 K (62.4, ref 62) would change these corrections by less than 0.03 eV. In this work, the energy levels were referred to the standard hydrogen electrode, for which we took the level calculated in ref 39 at the same level of theory. With this choice of alignment, the conduction band minimum and the valence band maximum of liquid water occur at −2.5 and 1.7 eV, respectively. We validate our method in the case of the reduction of the Fe3+ aqua ion:

(1)

where ϵF is the instantaneous Fermi energy. Fe drives electrons into the system when ϵF < ϵF̅ and drives electrons into the reservoir when ϵF > ϵF̅ so that the instantaneous Fermi energy ϵF averages out to the preset ϵF̅ . The dynamical equations for the charge evolution then read: Nė =

Pe and Pė = Fe = −(ϵF − ϵF̅ ) Me

(2)

where Pe is a fictitious momentum associated with the dynamical variable Ne. In this technique, the electronic charge in the system is controlled in an analogous way as the temperature and the pressure in constant-temperature47−50 and constant-pressure52 molecular dynamics, respectively. When additionally a thermostat setting the atomic degrees of freedom ri to a temperature T is used, it has been demonstrated that the grand-canonical distribution is sampled (cf. Supporting Information in ref 42). To control the fluctuations of Ne and consequently those of ϵF, it is also possible to couple the charge dynamics to a separate thermostat set at the temperature Te.42 This does not create any significant thermal flow in the system due to the weak coupling between the atoms and the charges.42 For a good choice of Me and Te, the instantaneous Fermi level of the system is governed by the preset ϵF̅ . In the case of redox reactions, performing MD at constant Fermi level corresponds to defining the net charge in the system. Hence, by controlling ϵF̅ , one can drive the reaction between the reactant in a charge state q and the products in a charge state q′ and get the structural transformations induced by the charge transfer. 2.2. Application to the Aqueous Fe2+/Fe3+ Redox Couple. In the present work, the electronic structure was described within density functional theory (DFT) through the rVV10 functional,53,54 which is based on a generalized gradient approximation for the exchange-correlation energy and augmented to include nonlocal van der Waals interactions. Within this DFT scheme, the short-range interactions are controlled via an empirical parameter b. We adopted the value of b = 9.3, as this has been shown to reproduce the experimental mass density in the case of liquid water.55 The core−valence interactions were described by normconserving pseudopotentials according to the prescription of Troullier and Martins.56 The wave functions of the valence electrons were expanded on a plane-wave basis set defined by a kinetic energy cutoff of 80 Ry. The Brillouin zone was sampled at the Γ point. We performed Born−Oppenheimer molecular dynamics using a time step of Δt = 0.48 fs to integrate the equations of motion. We used a velocity rescaling method to set the temperature at T = 350 K. We used this temperature to ensure liquid-like behavior.55 The Fe aqua ion was simulated with 31 water molecules in a cubic supercell with a side of 9.85 Å. We used the suite of codes provided in the Quantum-ESPRESSO package.57 In the fictitious charge dynamics, we set the fictitious electronic mass to Me = 1000 au The temperature of the electronic charge was controlled by a velocity rescaling method and was set to Te = 50 K. We used Gaussian-smeared occupations with widths of at most 0.15 eV to prevent numerical instabilities associated with the use of fractional

Fe3 + + e− → Fe 2 +

(3) 3+

After 3 ps of equilibration of Fe at 350 K, the last configuration has been taken as the starting point of several runs at different Fermi energies ϵF̅ , as shown in Figure 1. The choice of the target Fermi energies is guided by the energy level

Figure 1. Time evolution of the instantaneous Fermi level ϵF (magenta lines) and total electronic charge q (green lines) in constant Fermilevel molecular dynamics of the aqueous Fe2+/Fe3+ redox couple. The panels correspond to four different preset values of the Fermi energy ϵ̅F: −0.3, +0.2, +1.0, and +1.3 eV. The average values after equilibration are given for both ϵF and q̅ (dashed lines). The energies are referred to the SHE. C


Article

Journal of Chemical Theory and Computation of the lowest unoccupied single-particle state of Fe3+. The occupation number of this state goes from 0 to 1 as the reduction proceeds. In practice, the instantaneous Fermi energy goes within ∼1 ps to the target value ϵF̅ by taking on a fractional number of electrons from the potentiostat. After this period, the time averages of the instantaneous Fermi energy are taken over at least 10 ps and are found to agree with the target value within at most 0.02 eV. However, the charge dynamics are slower and require longer equilibration times (at most ∼6 ps). The charge fluctuations depend directly on the choice of the fictitious mass parameter Me and the fictitious electronic temperature Te. Me controls the frequency of the fluctuations while Te their amplitude. In our setup, we chose a sufficiently large Me to allow for structural rearrangements as the change in the electronic charge is accommodated. Statistical averages are taken over the equilibrated parts of the trajectories and correspond to periods of at least 10 ps. Figure 2a shows the charge transition as a function of the preset

the use of a quadratic model for capturing this dependence. This aspect is further discussed in section 3.2. It is interesting to study the structural properties as the system evolves from Fe3+ to Fe2+. In all simulations, the Fe aqua ion is octohedrally coordinated by six water molecules in the first coordination shell, in agreement with experimental observations.63 In Figure 3, we give the partial pair correlation

Figure 3. Partial pair correlation functions gFe−O(r) and gFe−H(r) for Fermi energies preset at −0.3, +0.2, +1.0, and +1.3 eV with respect to SHE.

functions Fe−O and Fe−H as calculated for our simulations at fixed ϵF̅ . As the total electronic charge of the system decreases from +3 to +2, the main peak in the gFe−O(r) shifts steadily from 1.98 to 2.09 Å due to electrostatic effects. These values are in excellent agreement with the estimates of 1.97 and 2.10 Å derived from extended X-ray absorption fine structure (EXAFS) experiments.63 Our findings are also consistent with results of previous simulations.64−67 Similarly, the first peak in the gFe−H(r) shifts toward shorter distances by 0.04 Å. Hence, we conclude that the present simulations at constant Fermi energy yield a very good structural description of the reduced and oxidized species.

3. REDOX POTENTIAL FROM CONSTANT FERMI-LEVEL SIMULATIONS 3.1. Thermodynamic Integration Method. Before going into the determination of redox potentials from constant Fermi-level simulations, we briefly review the thermodynamic integration method. In this work, we compute redox potentials following the formulation introduced by Ambrosio et al.39 Let us consider the following redox reaction:

Figure 2. Constant Fermi-level molecular dynamics of the aqueous Fe2+/Fe3+ redox couple: (a) evolution of the electronic charge (left axis) and of the occupation number (right axis) as a function of the preset Fermi energy ϵF̅ ; (b) evolution of the single-particle energy level ϵg as a function of ϵF̅ ; (c) ϵg as a function of occupation number η. The energies are referred to the SHE.

X − → X • + e−

(4)

The Gibbs free energy of formation of a solute X in the charge state q is given by

value of the Fermi energy ϵF̅ . Through the appearance of fractional charges, the redox reaction is driven toward the formation of Fe2+. We remark that the charge of the system is fully determined by the occupation η of the single-particle energy level in the gap, which is the only state undergoing an occupation change in our simulation. During this evolution, the single-particle energy level ϵg practically corresponds to ϵF̅ (cf. Figure 2b), the small differences resulting from the width of the employed Gaussian smearing. From Figures 2a,b, it is possible to extract the dependence of ϵg on η (Figure 2c). The dependence of ϵg on η shows a concave pattern and suggests

Gfq[X ] = Gq[X ] − G[bulk] − +

q Ecorr

∑ niμi + q(ϵv + μ) (5)

q

where G [X] and G[bulk] are the Gibbs free energies of the solute X in the charge state q and of the pristine bulk system, respectively, μi is the chemical potential of species i, μ the electron chemical potential, ϵv the valence band maximum of the bulk system, and Eqcorr an electrostatic finite-size correction term. The free-energy difference pertaining to the reaction in eq 4 is given by D


Article

Journal of Chemical Theory and Computation

that during the reaction a single state g is involved and that its occupation goes from 0 to 1. The electronic charge then corresponds to the occupation η of state g. In this description, the TI method gives the redox level through the calculation of the thermodynamic integral ΔA, as given in eq 8. Using Janak’s theorem in eq 11, the thermodynamic integral can be expressed as

ΔG(q/q′) = Gfq ′[X •] − Gfq[X −] = Gq ′[X •] − Gq[X −] + (ϵv + μ)(q′ − q) q′ q + Ecorr − Ecorr

(6)

The redox level then corresponds to the electron chemical potential for which ΔG(q/q′) = 0. By applying this condition to eq 6, we derive the following expression for the redox level defined with respect to ϵv: q

μ(q/q′) =

q′

−

•

G [X ] − G [X ] + q′ − q −

q Ecorr

q′ Ecorr

− q′ − q

− ϵv

ΔA =

ΔA = A(X −) − A(X •) =

∫0

q

(7)

•

where η is the Kirkwood coupling parameter and ⟨ΔE⟩η is the total energy difference between charge states q and q′ averaged over configurations achieved at fixed η. Substituting eq 8 into eq 7, we obtain μ(q/q′) =

ΔA + q′ − q

− q′ − q

− ϵv

ϵ(f , η) df dη

(12)

(13)

∫0 ∫0

1

⎛1 1⎞ ϵ(f , η) df dη ≅ ϵ⎜ , ⎟ = ϵS ⎝2 2⎠

(9)

Figure 4. Schematic representation of the function ϵ(f, η) in the twodimensional (f, η) space. Assuming linear dependence for ϵ( f, η) results in considering the function in the sole

∫0

( 12 , 12 ) point (black dot).

For capturing the quadratic dependence of ϵ( f, η), the dots shown in red are sampled in the present work.

(10)

As seen in section 2.2, ϵ( f, η) is not necessarily linear. A more accurate approximation then consists in assuming a quadratic functional form in both η and f:

1

ϵg (f ) df

(14)

A single molecular dynamics simulation at half filling of the state g would thus be sufficient at this level of approximation (cf. Figure 4).

In its integral form, Janak’s theorem gives the change in energy of the system upon addition of one electron at fixed geometry, EN+1 − EN =

ϵg (f ) df ⟩η dη

1

1

ΔA =

In practice, various simulations are performed with values of the coupling parameter η ranging between 0 and 1. For each system, the total energy difference between charge states q and q′ (vertical energy gap) is obtained as a time average. Consequently, ΔA is obtained as an integral of the vertical energy gap over η. When referred to the SHE, the redox levels calculated with eq 9 correspond to redox potentials and can directly be compared with experiment. 3.2. Redox Potential from Single-Particule Eigenvalues. According to Janak’s theorem,51 the derivative of the DFT energy functional E[N] with respect to the occupation f of the highest occupied Kohn−Sham (KS) state g gives the KS eigenvalue ϵg corresponding to that state: ∂E[N ] = ϵg ∂f

1

In the latter expression, the average is taken over structures corresponding to a state g with occupation η, while f indicates the occupation of state g at fixed geometry. In analogy to Slater’s approximation,68,69 the simplest approximation to the thermodynamic integral consists in assuming a linear dependence of ϵ( f, η) in η and f. This gives

⟨ΔE⟩η dη

q′ Ecorr

∫0 ∫0

∫0 ⟨∫0

ϵ(f , η) = ⟨ϵg (f )⟩η

(8)

q Ecorr

1

⟨ΔE⟩η dη =

where we have used the following definition:

1

= Gq ′[X −] − Gq[X •]

1

1

=

The free-energy difference G [X ] − G ′[X ] can be calculated from the thermodynamic integral (ΔA) of the vertical energy gap as follows: q

∫0

(11)

The constant-Fermi-level technique consists in driving the electronic charge from state q to state q′ through the realization of several intermediate systems at fixed Fermi level. We remark that semilocal energy functionals generally give a poor physical description of fractionally occupied states. Nevertheless, the use of these states acquires significance through Janak’s theorem (eq 11), as they establish a continuous connection between physical states carrying an integer number of electrons. Although it is the Fermi level that is controlled during the simulations, the parameters of the simulation can be set to allow only for minimal fluctuations of the electronic charge of the system, as seen in section 2.2. Then, taking the electronic charge of the system as Kirkwood coupling parameter allows one to draw a clear similarity with the TI method. We assume

⎛ ⎛ ⎛ 1⎞ 1⎞ 1 ⎞2 ϵ(f , η) = ϵS + α⎜η − ⎟ + β ⎜f − ⎟ + α′⎜η − ⎟ ⎝ ⎝ ⎝ 2⎠ 2⎠ 2⎠ 2 ⎛ ⎛ 1⎞ 1 ⎞⎛ 1⎞ + β′⎜f − ⎟ + γ ⎜η − ⎟⎜f − ⎟ ⎝ ⎝ (15) 2⎠ 2 ⎠⎝ 2⎠

where ϵS, α, β, α′, β′, and γ are constants. Substituting eq 15 into eq 12 and carrying out the integrals leads to the following expression for the thermodynamic integral: 1

ΔA =

∫0 ∫0

1

ϵ(f , η) df dη = ϵS +

α′ + β ′ 12

(16)

To determine ϵS, α′, and β′, we propose calculating ϵ( f, η) in the following cases. First, it is natural to determine ϵ( f, η) for f E


Article


− Fe2+, O2H−/O2H•, and MnO2− 4 /MnO4 . We carry out a detailed comparison between the redox potentials as obtained through the TI method and as extracted from the constant Fermi level simulations. The simulation details described in section 2.2 for the Fe3+/Fe2+ pair were also adopted for the other redox pairs. However, to account for the larger solute sizes, we used 31 and 28 water molecules in the simulation cell for O2H−/O2H• and − MnO2− 4 /MnO4 , respectively. For each redox pair, we performed constant Fermi-level simulations at fixed ϵF̅ , corresponding to various values η of occupation of the state g. These simulations could be used directly for determining the redox potentials through the TI method. For the configurations obtained during the molecular dynamics at a given η, we calculated the average total energy difference ⟨ΔE⟩η between occupations f = 1 and f = 0 (cf. Figure 5). Integration of ⟨ΔE⟩η yielded the thermodynamic

= η as this follows directly from the constant Fermi level simulations: ⎛ ⎛ 1⎞ 1 ⎞2 ϵ(η , η) = ϵS + ( α + β) ⎜η − ⎟ + ( α ′ + β′  + γ ) ⎜η − ⎟ ⎝ ⎝ 2⎠ 2⎠ A

⎛ ⎛ 1⎞ 1⎞ = ϵS + A⎜η − ⎟ + B⎜η − ⎟ ⎝ ⎝ 2⎠ 2⎠

B 2

(17)

where the parameters ϵS, A, and B are easily obtained from a quadratic fit of ϵ(η, η). However, this fit is not sufficient to determine ΔA in eq 16. Therefore, we consider two additional points in the two-dimensional ( f, η) space: ϵ(0, 1) = ϵS −

1 1 (α − β ) + (α′ + β ′ − γ ) 2 4

(18)

ϵ(1, 0) = ϵS +

1 1 (α − β ) + (α′ + β ′ − γ ) 2 4

(19)

Combining eqs 17, 18, and 19 with eq 16, we can express the thermodynamical integral as 1

ΔA = =

∫0 ∫0

1

ϵ(f , η) df dη

ϵ(0, 1) + ϵ(1, 0) 5 B ϵS + + 6 24 12

(20)

If ϵ(f, η) is given by a quadratic function, the present expression for the thermodynamic integral leads to the exact value of the redox level defined in eq 9. Figure 4 schematically illustrates the function ϵ(f, η) in the two-dimensional ( f, η) space and the points sampled in the present scheme. We remark that the presently adopted second-order form (eq 15) goes beyond Marcus’ linear regime.70,71 Indeed, on the basis of Janak’s theorem, the integration of ϵ(f, η) over f produces the vertical energy gaps ⟨ΔE⟩η required in the TI method (eq 8): ⟨ΔE⟩η =

∫0

1

Δϵ(f , η) df

⎛ ⎛ β′ 1⎞ 1 ⎞2 = ϵS + α⎜η − ⎟ + α′⎜η − ⎟ + ⎝ ⎝ 2⎠ 2⎠ 12

(21)

which by construction accounts for one higher order more than Marcus’ linear regime. However, we note that having recourse to a quadratic form for ϵ( f, η) is not inherent to the constant Fermi level approach. Indeed, when strong nonlinearities occur, higher-order forms should be fitted to the data. Consequently, this entails a more detailed sampling of the function ϵ( f, η), in analogy to the more detailed sampling required for the vertical energy gaps in the TI method. It is interesting to compare the amount of calculations required in the present scheme compared to the TI method. The required number of intermediate states should be set with the motivation of properly describing the evolution of the Kohn−Sham energy level upon Fermi energy variation, just like the evolution of vertical energy gaps upon varying the Kirkwood parameter in TI. The TI method additionally requires the calculation of vertical energy gaps for each considered state, while such kind of calculations are performed only for the initial and final states in the present scheme. 3.3. Determination of Redox Potentials. In this section, we focus on the redox potentials of three redox pairs: Fe3+/

Figure 5. Energetic diagram as a function of the occupation η of the − − state g for the redox couples Fe2+/Fe3+, MnO2− 4 /MnO4 , and O2H / O2H•. The diagram represents the KS eigenvalue of state g, ϵ(η, η), as obtained during the constant Fermi-level simulation (black lines with dots), the KS eigenvalues ϵ(0, η) (blue lines with triangles) and ϵ(1, η) (red lines with triangles), and the total energy difference ⟨ΔE⟩η (magenta lines with stars).

integral ΔA, from which the redox potential was derived through eq 9. This result is taken as reference value in the following. For comparison, we also estimated the thermodynamic integral through Marcus’ linear picture,70,71 in which the vertical energy differences at the end points (at η = 0 and η = 1) are connected linearly. The obtained values are listed in Table 1. Following the scheme outlined in section 3.2, we also determined the redox potentials via the averaged Kohn−Sham F


Article


ensured through a uniform background. The technique allows one to directly control the charge transfer from an electron reservoir at a preset potential to the single-particle energy levels of the system, thereby modeling the effect of an electrode. In particular, the present scheme removes any simulation bias due to compensating ions or to a counterelectrode. Furthermore, the redox reaction can be driven through the control of the Fermi level without any prior knowledge of the reaction products. The charge evolution of the system is governed by fictitious dynamical equations of motion, which are solved separately with respect to the electronic structure calculations. For aqua ions, we demonstrated molecular dynamics with stable evolutions of the Fermi level and of the electronic charge, without any divergence or instability. We showed that the structural properties of the hydrated ions agree well with experiment. Thus, constant Fermi level simulations constitute a valuable tool for describing structural mechanisms associated with charge transfer. In the second part, we focused on the determination of redox potentials from constant Fermi level simulations. We monitored the evolution of a single particle energy level as a function of its occupation. Through Janak’s theorem, a relation with the redox potential could be established. Introducing a quadratic function to model the dependence of the singleparticle energy level on occupation, we found agreement within less than 0.1 eV with redox potentials calculated through the thermodynamic integration method. We also compared the calculated redox potentials with their experimental counterparts, finding agreement within errors typical of the level of DFT adopted in this work. To yield higher accuracy,37,38 the present scheme can be carried over without impediment to higher levels of theory, such as hybrid functionals. Constant Fermi level simulations are not only limited to the quantitative study of redox potentials in aqueous solution but could also be used to determine redox potentials at electrode/ water interfaces. Through the control of the Fermi level, the oxidation mechanism as well as the associated structural changes could realistically be modeled. Such processes will shed light in an unprecedented fashion on the electrochemistry occurring at electrode surfaces.

Table 1. Redox Potentials vs SHE (in eV) As Obtained through Constant Fermi Level Simulations in the Linear and Quadratic Approximations, Compared to Those Achieved with the Thermodynamic Integration (TI) Method, Both in the Linear End-Point Approximation and through a Full Sampling (Reference Value)a constant Fermi energy

a

TI

expt

system

linear

quadratic

linear

full

O2H•/O2H− Fe2+/Fe3+ MnO−4 /MnO−2 4

0.58 1.10 −0.06

0.47 0.85 −0.31

0.22 0.65 −0.30

0.42 0.77 −0.28

0.75 0.77 0.56

Experimental values are shown for comparison.

eigenvalues ϵ( f, η) of state g. The constant Fermi-level simulation provided us directly with ϵ(η, η), as shown in 1 1 Figure 5. For the linear approximation, only ϵ 2 , 2 is

(

)

required. The obtained values are found to differ from the reference TI results up to 0.33 eV (cf. Table 1). In the more accurate quadratic approximation, we additionally calculated ϵ(0, 1) and ϵ(1, 0). These values are reported in Figure 5, where we additionally show the full dependence of ϵ(0, η) and ϵ(1, η) on η for the sole purpose of illustration. Through eq 20, the quadratic approximation yields an improved description of the thermodynamic integrals. The corresponding redox potentials are given in Table 1. Comparison with the reference TI values shows agreement within 0.08 eV, validating the proposed scheme for the calculation of redox potentials. The present redox couples also illustrate that the quadratic form in the present approach goes beyond Marcus’ linear end-point approximation. Indeed, the redox potentials evaluated in the latter scheme generally show larger deviations with respect to the reference result, reaching up to 0.20 eV in the case of the redox couple O2H•/ O2H−. All redox couples give redox levels falling inside the reduced band gap achieved in the present DFT calculations. The calculated redox potentials can thus be compared with their experimental counterparts at 300 K (Table 1). The use of a temperature of 350 K in our simulations is not expected to invalidate this comparison.39 Considering the reference results obtained with the TI method, we found an agreement with experiment within 0.3 eV for O2H•/O2H− and Fe2+/Fe3+ and a larger deviation of 0.84 eV for MnO−4 /MnO2− 4 . This level of agreement can be considered typical for DFT schemes at the present level of theory, as errors up to 0.9 eV have been found in redox potentials calculated with GGA functionals in previous studies.36,38 In particular, the large deviation observed for the permanganate redox potential could in part also arise from the choice of the volume that we allocated to the MnO−4 ion in the absence of any estimate of its partial molar volume. However, we remark that these deficiencies do not affect the main achievement in our work, which consists in validating constant Fermi-level simulations as a practical tool for the calculation of redox potentials.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: assil.bouzid@epfl.ch. ORCID

Assil Bouzid: 0000-0002-9363-7240 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This project has received funding from the European Union’s Seventh Framework Programme for research, technological development, and demonstration under grant agreement no. 291771 (call 2014). This work has been performed in the context of the National Center of Competence in Research (NCCR) “Materials’ Revolution: Computational Design and Discovery of Novel Materials (MARVEL)” of the Swiss National Science Foundation. We used computational resources of CSCS and CSEA-EPFL.

4. CONCLUSIONS In this work, we developed a computational scheme to compute redox potentials through constant Fermi-level simulations. In the first part of our study, we adapted the constant Fermi-level simulation technique to be used with three-dimensional periodic boundary conditions and with charge neutrality G


Article


■

(21) Tavernelli, I.; Vuilleumier, R.; Sprik, M. Ab initio molecular dynamics for molecules with variable numbers of electrons. Phys. Rev. Lett. 2002, 88, 213002. (22) Blumberger, J.; Bernasconi, L.; Tavernelli, I.; Vuilleumier, R.; Sprik, M. Electronic structure and solvation of copper and silver ions: a theoretical picture of a model aqueous redox reaction. J. Am. Chem. Soc. 2004, 126, 3928−3938. (23) Blumberger, J.; Sprik, M. Free energy of oxidation of metal aqua ions by an enforced change of coordination. J. Phys. Chem. B 2004, 108, 6529−6535. (24) Tateyama, Y.; Blumberger, J.; Sprik, M.; Tavernelli, I. Densityfunctional molecular-dynamics study of the redox reactions of two anionic, aqueous transition-metal complexes. J. Chem. Phys. 2005, 122, 234505. (25) Blumberger, J.; Sprik, M. Ab initio molecular dynamics simulation of the aqueous Ru2+/Ru3+ redox reaction: The Marcus perspective. J. Phys. Chem. B 2005, 109, 6793−6804. (26) VandeVondele, J.; Ayala, R.; Sulpizi, M.; Sprik, M. Redox free energies and one-electron energy levels in density functional theory based ab initio molecular dynamics. J. Elect. Chem. 2007, 607, 113− 120. (27) Kirkwood, J. G. Statistical mechanics of fluid mixtures. J. Chem. Phys. 1935, 3, 300−313. (28) Cheng, J.; Sulpizi, M.; Sprik, M. Redox potentials and pKa for benzoquinone from density functional theory based molecular dynamics. J. Chem. Phys. 2009, 131, 154504. (29) Costanzo, F.; Sulpizi, M.; Della Valle, R. G.; Sprik, M. The oxidation of tyrosine and tryptophan studied by a molecular dynamics normal hydrogen electrode. J. Chem. Phys. 2011, 134, 244508. (30) Blumberger, J.; Tavernelli, I.; Klein, M. L.; Sprik, M. Diabatic free energy curves and coordination fluctuations for the aqueous Ag+/ Ag2+ redox couple: A biased Born-Oppenheimer molecular dynamics investigation. J. Chem. Phys. 2006, 124, 064507. (31) Sulpizi, M.; Sprik, M. Acidity constants from vertical energy gaps: density functional theory based molecular dynamics implementation. Phys. Chem. Chem. Phys. 2008, 10, 5238−5249. (32) Cheng, J.; Liu, X.; Kattirtzi, J. A.; VandeVondele, J.; Sprik, M. Aligning electronic and protonic energy levels of proton-coupled electron transfer in water oxidation on aqueous TiO2. Angew. Chem., Int. Ed. 2014, 53, 12046−12050. (33) Adriaanse, C.; Cheng, J.; Chau, V.; Sulpizi, M.; VandeVondele, J.; Sprik, M. Aqueous redox chemistry and the electronic band structure of liquid water. J. Phys. Chem. Lett. 2012, 3, 3411−3415. (34) Cheng, J.; Sprik, M. Alignment of electronic energy levels at electrochemical interfaces. Phys. Chem. Chem. Phys. 2012, 14, 11245− 11267. (35) Sulpizi, M.; Sprik, M. Acidity constants from DFT-based molecular dynamics simulations. J. Phys.: Condens. Matter 2010, 22, 284116. (36) Cheng, J.; Liu, X.; VandeVondele, J.; Sulpizi, M.; Sprik, M. Redox potentials and acidity constants from density functional theory based molecular dynamics. Acc. Chem. Res. 2014, 47, 3522−3529. (37) Ambrosio, F.; Miceli, G.; Pasquarello, A. Structural, dynamical, and electronic properties of liquid water: A hybrid functional study. J. Phys. Chem. B 2016, 120, 7456−7470. (38) Cheng, J.; VandeVondele, J. Calculation of electrochemical energy levels in water using the random phase approximation and a double hybrid functional. Phys. Rev. Lett. 2016, 116, 086402. (39) Ambrosio, F.; Miceli, G.; Pasquarello, A. Redox levels in aqueous solution: Effect of van der Waals interactions and hybrid functionals. J. Chem. Phys. 2015, 143, 244508. (40) Agostiano, A. Cyclic voltammetry: A powerful tool to follow redox processes in biologically relevant systems. J. Photochem. Photobiol., B 1993, 17, 294−297. (41) Lozovoi, A. Y.; Alavi, A.; Kohanoff, J.; Lynden-Bell, R. M. Ab initio simulation of charged slabs at constant chemical potential. J. Chem. Phys. 2001, 115, 1661−1669.

REFERENCES

(1) Mennucci, B.; Tomasi, J.; Cammi, R.; Cheeseman, J.; Frisch, M.; Devlin, F.; Gabriel, S.; Stephens, P. Polarizable continuum model (PCM) calculations of solvent effects on optical rotations of chiral molecules. J. Phys. Chem. A 2002, 106, 6102−6113. (2) Baik, M.-H.; Schauer, C. K.; Ziegler, T. Density functional theory study of redox pairs: 2. Influence of solvation and ion-pair formation on the redox behavior of cyclooctatetraene and nitrobenzene. J. Am. Chem. Soc. 2002, 124, 11167−11181. (3) Dutton, A. S.; Fukuto, J. M.; Houk, K. N. Theoretical reduction potentials for nitrogen oxides from CBS-QB3 energetics and (C) PCM solvation calculations. Inorg. Chem. 2005, 44, 4024−4028. (4) Yu, A.; Liu, Y.; Li, Z.; Cheng, J.-P. Computation of pKa values of substituted aniline radical cations in dimethylsulfoxide solution. J. Phys. Chem. A 2007, 111, 9978−9987. (5) Shimodaira, Y.; Miura, T.; Kudo, A.; Kobayashi, H. DFT method estimation of standard redox potential of metal ions and metal complexes. J. Chem. Theory Comput. 2007, 3, 789−795. (6) Eckert, F.; Leito, I.; Kaljurand, I.; Kütt, A.; Klamt, A.; Diedenhofen, M. Prediction of acidity in acetonitrile solution with COSMO-RS. J. Comput. Chem. 2009, 30, 799−810. (7) Roy, L. E.; Jakubikova, E.; Guthrie, M. G.; Batista, E. R. Calculation of one-electron redox potentials revisited. Is it possible to calculate accurate potentials with density functional methods? J. Phys. Chem. A 2009, 113, 6745−6750. (8) Gómez-Bombarelli, R.; González-Pérez, M.; Pérez-Prior, M. T.; Calle, E.; Casado, J. Computational study of the acid dissociation of esters and lactones. A case study of diketene. J. Org. Chem. 2009, 74, 4943−4948. (9) Zeng, Y.; Qian, H.; Chen, X.; Li, Z.; Yu, S.; Xiao, X. Thermodynamic estimate of pKa values of the carboxylic acids in aqueous solution with the density functional theory. Chin. J. Chem. 2010, 28, 727−733. (10) Eckert, F.; Diedenhofen, M.; Klamt, A. Towards a first principles prediction of pKa: COSMO-RS and the cluster-continuum approach. Mol. Phys. 2010, 108, 229−241. (11) Jono, R.; Sumita, M.; Tateyama, Y.; Yamashita, K. Redox reaction mechanisms with non-triiodide mediators in dye-sensitized solar cells by redox potential calculations. J. Phys. Chem. Lett. 2012, 3, 3581−3584. (12) Guerard, J. J.; Arey, J. S. Critical evaluation of implicit solvent models for predicting aqueous oxidation potentials of neutral organic compounds. J. Chem. Theory Comput. 2013, 9, 5046−5058. (13) Yang, C.; Xue, X.-S.; Jin, J.-L.; Li, X.; Cheng, J.-P. Theoretical study on the acidities of chiral phosphoric acids in dimethyl sulfoxide: hints for organocatalysis. J. Org. Chem. 2013, 78, 7076−7085. (14) Sviatenko, L. K.; Gorb, L.; Hill, F. C.; Leszczynski, J. Theoretical study of ionization and one-electron oxidation potentials of Nheterocyclic compounds. J. Comput. Chem. 2013, 34, 1094−1100. (15) Yan, L.; Lu, Y.; Li, X. A density functional theory protocol for the calculation of redox potentials of copper complexes. Phys. Chem. Chem. Phys. 2016, 18, 5529−5536. (16) Baik, M.-H.; Friesner, R. A. Computing redox potentials in solution: Density functional theory as a tool for rational design of redox agents. J. Phys. Chem. A 2002, 106, 7407−7412. (17) Srnec, M.; Chalupský, J.; Fojta, M.; Zendlová, L.; Havran, L.; Hocek, M.; Kývala, M.; Rulisek, L. Effect of spin-orbit coupling on reduction potentials of octahedral ruthenium (II/III) and osmium (II/ III) complexes. J. Am. Chem. Soc. 2008, 130, 10947−10954. + 2+ /Cuaq redox reaction exhibits strong (18) Blumberger, J. Cuaq nonlinear solvent response due to change in coordination number. J. Am. Chem. Soc. 2008, 130, 16065−16068. (19) Keith, T. A.; Frisch, M. J. Inclusion of explicit solvent molecules in a self-consistent-reaction field model of solvation. Modeling the hydrogen bond 1994, 569, 22−35. (20) Uudsemaa, M.; Tamm, T. Density-functional theory calculations of aqueous redox potentials of fourth-period transition metals. J. Phys. Chem. A 2003, 107, 9997−10003. H


Article

Journal of Chemical Theory and Computation (42) Bonnet, N.; Morishita, T.; Sugino, O.; Otani, M. First-principles molecular dynamics at a constant electrode potential. Phys. Rev. Lett. 2012, 109, 266101. (43) Otani, M.; Sugino, O. First-principles calculations of charged surfaces and interfaces: A plane-wave nonrepeated slab approach. Phys. Rev. B: Condens. Matter Mater. Phys. 2006, 73, 115407. (44) Jinnouchi, R.; Anderson, A. B. Electronic structure calculations of liquid-solid interfaces: Combination of density functional theory and modified Poisson-Boltzmann theory. Phys. Rev. B: Condens. Matter Mater. Phys. 2008, 77, 245417. (45) Otani, M.; Hamada, I.; Sugino, O.; Morikawa, Y.; Okamoto, Y.; Ikeshoji, T. Structure of the water/platinum interface - A first principles simulation under bias potential. Phys. Chem. Chem. Phys. 2008, 10, 3609−3612. (46) Otani, M.; Hamada, I.; Sugino, O.; Morikawa, Y.; Okamoto, Y.; Ikeshoji, T. Electrode dynamics from first principles. J. Phys. Soc. Jpn. 2008, 77, 024802. (47) Nosé, S. A molecular dynamics method for simulations in the canonical ensemble. Mol. Phys. 1984, 52, 255−268. (48) Nosé, S. A unified formulation of the constant temperature molecular dynamics methods. J. Chem. Phys. 1984, 81, 511−519. (49) Hoover, W. G. Canonical dynamics: equilibrium phase-space distributions. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 31, 1695−1697. (50) Blöchl, P. E.; Parrinello, M. Adiabaticity in first-principles molecular dynamics. Phys. Rev. B: Condens. Matter Mater. Phys. 1992, 45, 9413−9416. (51) Janak, J. F. Proof that ∂E/∂ni=ϵ in density-functional theory. Phys. Rev. B: Condens. Matter Mater. Phys. 1978, 18, 7165. (52) Andersen, H. C. Molecular dynamics simulations at constant pressure and/or temperature. J. Chem. Phys. 1980, 72, 2384−2393. (53) Vydrov, O. A.; Van Voorhis, T. Nonlocal van der Waals density functional: The simpler the better. J. Chem. Phys. 2010, 133, 244103. (54) Sabatini, R.; Gorni, T.; de Gironcoli, S. Nonlocal van der Waals density functional made simple and efficient. Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 87, 041108. (55) Miceli, G.; de Gironcoli, S.; Pasquarello, A. Isobaric firstprinciples molecular dynamics of liquid water with nonlocal van der Waals interactions. J. Chem. Phys. 2015, 142, 034501. (56) Troullier, N.; Martins, J. L. Efficient pseudopotentials for planewave calculations. Phys. Rev. B: Condens. Matter Mater. Phys. 1991, 43, 1993−2006. (57) Giannozzi, P.; Baroni, S.; Bonini, N.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Chiarotti, G. L.; Cococcioni, M.; Dabo, I.; Dal Corso, A.; de Gironcoli, S.; Fabris, S.; Fratesi, G.; Gebauer, R.; Gerstmann, U.; Gougoussis, C.; Kokalj, A.; Lazzeri, M.; Martin-Samos, L.; Marzari, N.; Mauri, F.; Mazzarello, R.; Paolini, S.; Pasquarello, A.; Paulatto, L.; Sbraccia, C.; Scandolo, S.; Sclauzero, G.; Seitsonen, A. P.; Smogunov, A.; Umari, P.; Wentzcovitch, R. M.; et al. QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials. J. Phys.: Condens. Matter 2009, 21, 395502. (58) Fu, C. L.; Ho, K. M. External-charge-induced surface reconstruction on Ag(110). Phys. Rev. Lett. 1989, 63, 1617. (59) Bohnen, K. P.; Kolb, D. M. Charge- versus adsorbate-induced lifting of the Au(100)-(hex) reconstruction in an electrochemical environment. Surf. Sci. 1998, 407, L629−L632. (60) Freysoldt, C.; Neugebauer, J.; Van de Walle, C. G. Fully ab initio finite-size corrections for charged-defect supercell calculations. Phys. Rev. Lett. 2009, 102, 016402. (61) Komsa, H.-P.; Rantala, T. T.; Pasquarello, A. Finite-size supercell correction schemes for charged defect calculations. Phys. Rev. B: Condens. Matter Mater. Phys. 2012, 86, 045112. (62) Malmberg, C. G.; Maryott, A. A. Dielectric constant of water from 0° to 100°C. J. Res. Nat. Bureau Stand. 1956, 56, 1−8. (63) Brunschwig, B. S.; Creutz, C.; Macartney, D. H.; Sham, T. K.; Sutin, N. The role of inner-sphere configuration changes in electronexchange reactions of metal complexes. Faraday Discuss. Chem. Soc. 1982, 74, 113−127. (64) Guàrdia, E.; Padró, J. Molecular dynamics simulation of ferrous and ferric ions in water. Chem. Phys. 1990, 144, 353−362.

(65) Drechsel-Grau, C.; Sprik, M. The temperature dependence of the symmetry factor for a model Fe3+(aq)/Fe2+(aq) redox half reaction. Mol. Phys. 2015, 113, 2463−2475. (66) Ando, K. Solvent nuclear quantum effects in electron transfer reactions. III. Metal ions in water. Solute size and ligand effects. J. Chem. Phys. 2001, 114, 9470−9477. (67) Amira, S.; Spångberg, D.; Probst, M.; Hermansson, K. Molecular dynamics simulation of Fe2+(aq) and Fe3+(aq). J. Phys. Chem. B 2004, 108, 496−502. (68) Slater, J. C.; Johnson, K. H. Self-consistent-field Xα cluster method for polyatomic molecules and solids. Phys. Rev. B 1972, 5, 844. (69) Slater, J. C. Statistical exchange-correlation in the self-consistent field. Adv. Quantum Chem. 1972, 6, 1−92. (70) Marcus, R. A. On the theory of electron-transfer reactions. VI. Unified treatment for homogeneous and electrode reactions. J. Chem. Phys. 1965, 43, 679−701. (71) Marcus, R. A. On the theory of oxidation-reduction reactions involving electron transfer. I. J. Chem. Phys. 1956, 24, 966−978.

I


Redox Levels through Constant Fermi-Level ab Initio Molecular

Recommend Documents