An energy decomposition and extrapolation scheme for evaluating electron transfer rate constants: a case study on electron self-exchange reactions of transition metal complexes

A simple approach to the analysis of electron transfer (ET) reactions based on energy decomposition and extrapolation schemes is proposed. The present energy decomposition and extrapolation-based electron localization (EDEEL) method represents the diabatic energies for the initial and final states using the adiabatic energies of the donor and acceptor species and their complex. A scheme for the efficient estimation of ET rate constants is also proposed. EDEEL is semi-quantitative by directly evaluating the seam-of-crossing region of two diabatic potentials. In a numerical test, EDEEL successfully provided ET rate constants for electron self-exchange reactions of thirteen transition metal complexes with reasonable accuracy. In addition, its energy decomposition and extrapolation schemes provide all the energy values required for activation-strain model (ASM) analysis. The ASM analysis using EDEEL provided rational interpretations of the variation of the ET rate constants as a function of the transition metal complexes. These results suggest that EDEEL is useful for efficiently evaluating ET rate constants and obtaining a rational understanding of their magnitudes.


Introduction
2][3][4][5][6] Theoretical tools are needed to elucidate and design novel reactivities and functions triggered by ET. [7][8][9][10] The ET rate constant is given by eqn (1). 9,11,12 where R is the gas constant, T is the temperature, p 0 is the standard atmosphere, q is either 1 or 0 for intermolecular or intramolecular ET, k el is the electronic transmission coefficient, n n is an effective nuclear frequency along the reaction coordinate, and DG ‡ is the Gibbs energy of activation on the adiabatic potential energy surface (PES).k el is given by eqn (2)- (4).
> > : where P LZ is the Landau-Zener transition probability, k B is the Boltzmann constant, h is the Planck constant, V 12 is the coupling constant, and l is the reorganization energy (Fig. 1).When V 12 is relatively small and 2pg ( 1, the Marcus theory equation 13 eqn (5) is obtained from Taylor series of eqn (3).
The harmonic approximation is oen used to describe the diabatic PESs of the initial and nal states, which allows one to estimate the Gibbs energy of activation DG ‡ na (subscript "na" means "non-adiabatic") by the following eqn (6) without explicitly locating the intersection between the two diabatic PESs on the multidimensional coordinate space.Note that the reaction barrier in eqn (1) is DG ‡ , not DG ‡ na (Fig. 1).
5][16][17][18][19][20][21][22][23] For example, the generalized Mulliken-Hush (GMH) theory is for ET and generates diabatic states based on the assumption that the initial and nal diabatic states are charge-localized states at different centers of the donor and acceptor. 24,25Constrained density functional theory (CDFT) is also useful for ET and forms diabatic states by imposing constraints on the partial charge (or spin) of arbitrary molecules or fragments using the Lagrange multiplier. 26,27o go beyond the framework based on the harmonic approximation, it is necessary to consider the anharmonic effect on the reaction coordinate.This requires explicit evaluations of PESs at around the non-adiabatic transition takes place.Molecular dynamics simulations can be performed to obtain accurate transition probabilities by calculating a large number of trajectories passing through such a region. 9,28-310][41][42][43][44][45] In this approach, DG ‡ na is estimated as the energy gap between the SX/ CI geometry and the initial state equilibrium (EQ) geometry.

Method
In this study, a simple diabatization scheme called energy decomposition and extrapolation-based electron localization (EDEEL) is proposed.In EDEEL, the diabatic energies for the initial and nal states V 11 and V 22 are given by eqn (7) and (8), respectively.
where R is the geometry of the donor-acceptor complex system, and R D and R A are its donor and acceptor components, respectively.The subscript for each E corresponds to the number of electrons in the corresponding system (n = m in the electron self-exchange reactions).E C n+m (R) is the adiabatic energy of the donor-acceptor complex without the moving electron.E D n (R) and E D n+1 (R) are the adiabatic energies of the donor species without and with the moving electron, respectively.E A m (R) and E A m+1 (R) are the adiabatic energies of the acceptor species without and with the moving electron, respectively.This scheme assumes that the moving electron is localized on either the donor or acceptor species and does not  interact with the corresponding counterpart.In other words, this scheme sets these assumptions as the conditions for diabatization.In the EDEEL scheme, an SX geometry between V 11 and V 22 is the critical point at which k ET is evaluated.
is the Gibbs energy gap between the SX and the initial state EQ geometries.Gibbs energy corrections at the EQ and SX are calculated by the normal mode analysis.The coupling V 12 is computed as the energy difference , where E C n+m+1 is the lower adiabatic energy between those for the two states that contain the contributions of the two diabatic states the most.Thus, at the SX between V 11 and V 22 , i.e., when V 11 = V 22 , the EDEEL scheme reproduces the diabatic-adiabatic relation of the two-state model in eqn (9).
In the electron self-exchange reactions, E C n+m+1 is the adiabatic energy of the donor-acceptor complex in its electronic ground state.The parameter l is obtained as eqn (10) based on eqn (6).
where relations DG ‡ na = G SX − G initial−EQ and DG 0 = G nal−EQ − G initial−EQ are used, G nal−EQ is Gibbs energy of the nal state EQ, and eqn (6) gives l = 4DG ‡ na for electron self-exchange reactions. 70he distance dependence of the ET reaction rate may have a distinct peak. 60This study regards adiabatic ground state at the SX geometry between V 11 and V 22 as the approximate transition state and maximizes the ET rate constant k ET within the SX region to obtain the nal k ET (=k calc ).This is done by taking the donor-acceptor distance r DA (the metal-metal distance in the transition metal examples below) as the reaction coordinate, evaluating the k ET at various r DA , and maximizing the k ET along r DA .Details of how the initial r DA (ini) is systematically determined and how the k ET is maximized are described in the ESI.† As shown in Fig. 2, the entire k calc evaluation workow consists of the following four steps.
(1) Systematic search for monomer conformations using the SC-AFIR 71,72 method implemented in the Global Reaction Route Mapping (GRRM) program. 732) Optimization of a minimum energy SX (MESX) geometry 74,75 within the SX hypersurface between V 11 and V 22 while maintaining the donor-acceptor distance to r DA (ini).
(3) k ET maximization along r DA by calculating k ET for ve r DA values near r DA (ini), where line search is done by simple quadratic curve tting using three k ET values in this study.In the calculation of k ET with eqn (1)-( 4), the Gibbs energies are extrapolated by the electronic energies obtained from single point calculations with larger basis set.

Results and discussion
The workow in Fig. 2 was applied to the electron self-exchange reactions of thirteen transition metal complexes listed in Table 1.Electronic structure calculations were performed with the Gaussian 16 program 76 and geometry optimization was performed with the GRRM program at the UuB97X-D/Def2-SV(P) 77,78 level taking into account the solvent effect of water by the conductor-like polarizable continuum model (C-PCM)., 79,80 where UuB97X-D stands for a spin unrestricted DFT calculation using the uB97X-D functional.k ET was maximized using the electronic energy calculated at the UuB97X-D/Def2-TZVP 81 level and the Gibbs energy correction at the UuB97X-D/Def2-SV(P) level; the calculation level is represented as UuB97X-D/Def2-TZVP//Def2-SV(P).In the calculation of k ET , n n in eqn (1) was approximated as k B T/h, the coefficient of transition state theory rate constant equation.The spin multiplicity was chosen to stabilize each state as much as possible and was set to the values listed in Table 1.The reaction set in Table 1 covers reactions of different timescales with rate constants in the wide range of 10 −7 to 10 9 dm 3 mol −1 s −1 , and thus would be suitable as a test set.It was assumed that both +3 and +2 charged complexes were in the electronic ground state at SX, although there were possibilities for ET via metal-to-ligand-charge transfer (MLCT) or electronic excited states of each complex. 82urther computational details are presented in Computational section.
In Table 2, the coupling values are relatively large, making k el almost unity.This is because our diabatization constraint assumes that there is no interaction between the electron to be transferred and the acceptor molecule.In other words, our constraint pushes the strong attraction between the negative charge on the moving electron and the +3 positive charge on the acceptor molecule into the coupling value, making the coupling values large.As for the non-adiabaticity, electron transfer reactions in [M(H 2 O) 6 ] 3+/2+ , which showed the largest coupling values among the systems in our calculation, are known as nonadiabatic process with smaller actual coupling values (<100 cm −1 ). 83The log 10 k calc does not include the diffusioncontrolled rate constant, which is estimated to be k diff = 3 × 10 9 dm 3 mol −1 s −1 for these reactions. 49In other words, k diff would be a source of the large error seen for fast reactions such as [Fe(bpy) 3 ] 3+/2+ and [Ru(bpy) 3 ] 3+/2+ .In our scheme, the C-PCM environments of the initial and nal states are relaxed independently, which does not satisfy the Franck-Condon principle.This would have led to the underestimation of reaction barriers.For example, the outer-sphere reorganization energy of [Ru(H 2 O) 6 ] 3+/2+ was estimated to be 0.68 eV based on dielectric continuum theory. 13This implies additional reaction barrier of  0.17 eV, reducing its log 10 k ET from 1.1 to −1.7, which is threeorder of magnitude smaller than the experimental value.The simple addition of the outers-sphere reorganization energy correction is not sufficient in our scheme.The better treatment of outer-sphere reorganization in our scheme needs further study in the future.Many other factors such as explicit solvation, dynamic and quantum motion of the atoms, and higher order (or static) electron correlation would contribute to the error.However, further improvement of the accuracy considering these factors is beyond the scope of this study.
Thanks to the energy decomposition and extrapolation scheme of EDEEL, the ASM analysis 67-69 can be performed using the energy components, i.e., E C n+m , E D n , E D n+1 , E A m , and E A m+1 , calculated during the EDEEL calculations, without any additional calculations.Note that the ASM analysis was performed using electronic energies rather than Gibbs energies for the SX geometry with the highest k ET value.In the ASM/EDEEL analysis, the strain energies in donor DE ‡ Strain(D) and acceptor DE ‡ Strain(A) and their interaction energy DE ‡ Interaction(C) are represented as follows. where m and E C(SX) n+m are the electronic energies at the SX geometry of the donor with the moving electron, that of the donor without the moving electron, that of the acceptor without the moving electron, and that of the donor-acceptor complex without the moving electron, and E D(initial−EQ) n+1 and E A(initial−EQ) m are the electronic energies at the initial state EQ geometry of the donor with the moving electron and of the acceptor without the moving electron.The donor and acceptor in the initial EQ geometry are assumed to be innitely far apart and not interacting.Fig. 4 and 5 show the results of the ASM/ EDEEL analysis at the UuB97X-D/Def2-TZVP level.The magnitudes of the strain energies of the aqua and cobalt complexes are large as shown in Fig. 4, reecting the large oxidation/redox potentials of their +2/+3 species.In these systems, to compensate for the large energy gaps between the +2 and +3 species, large structural deformations (Dl ‡ Strain(D) and Dl ‡ Strain(A) , changes in mean coordination bond length in donor and acceptor, respectively) are seen in the SXs around their central metal, as shown in Fig. 5.In pyridine and bipyridine complexes, attractive interactions are seen due to p-stacking between pyridine moieties and p-stacking between bipyridine molecules, respectively. 13,49On the other hand, aqua, ammonia, and ethylenediamine complexes show repulsive interactions, even though aqua complexes formed hydrogen bonding between water molecules.These interpretations by ASM/EDEEL would be benecial in understanding ET rate constants quantitatively.

Conclusion
In this article, we have proposed a simple diabatization scheme called EDEEL for the analysis of ET reactions.EDEEL represents the diabatic energies, V 11 and V 22 , by combining the adiabatic energies of the donor, acceptor, and their complex.Such a simple representation allows us to easily optimize the SX geometries between the two diabatic potentials.A scheme for estimating the ET rate constants at the SX geometries has also been introduced.The diabatic coupling V 12 is also estimated using the adiabatic energy of the donor-acceptor complex.In other words, the present scheme allows one to obtain all the parameters necessary to estimate ET rate constants using only the adiabatic energies.Numerical tests with electron selfexchange reactions of thirteen transition metal complexes have shown that EDEEL reproduces the trend of the experimental rate constants well and is semi-quantitative.
Two advantages of using EDEEL can be suggested.One is that EDEEL can be combined with any ab initio method and program without touching their codes.This is because all diabatic potential elements, i.e., V 11 , V 22 and V 12 , are represented by the adiabatic energies.The other is that EDEEL provides all the energy components and SX geometries needed in the ASM analysis.The latter is helpful in interpreting the ET efficiency and further designing a system with higher ET efficiency.In the present applications, the EDEEL-based ASM analysis successfully provided rational explanations for the variation of the magnitudes of the ET rate constants depending on the transition metal complexes.

Computational section
The above procedures were implemented in an in-house Python script as an interface program between the GRRM program and the Gaussian 16 program. 73,76The script takes a geometry from GRRM, performs the necessary electronic structure calculations to obtain the EDEEL PES at the geometry using Gaussian 16, and returns the energy and gradient (and Hessian if necessary) of the EDEEL PES at the geometry to GRRM.Geometry optimizations were performed by a developer version of the GRRM program at the UuB97X-D/Def2-SV(P) level.The "Stable = Opt" option was also used to identify the electronic ground state conguration of a given spin multiplicity.The "Int(Grid = 99 590)" option was used in DFT calculations.Although some of the aqua-complexes are highly acidic and may not prevail as simple [M(H 2 O) 6 ] 3+/2+ complexes in actual aqueous solution, the metal centers were assumed to be hexa-coordinated in this study as previous reports. 60,62,63,65Solvent water was modeled by the conductor-like polarizable continuum model (C-PCM). 79,80ibbs energy corrections were obtained from the standard normal mode analysis of 3N−8 dimensions (one direction orthogonal to SX and one direction along the vector between two metal atoms were removed from the full 3N−6 dimensions) at T = 298.15K and p 0 = 1 atm, with normal mode frequencies smaller than 100 cm −1 replaced by 100 cm −1 as suggested in the literature. 84The bases of Gibbs activation energies were calculated as the sum of the extrapolated Gibbs energies of the most stable donor and acceptor monomers, where the most stable structures were identied by a systematic conformation search using SC-AFIR with possible spin multiplicities.6][87] The rate constant maximization was done at the UuB97X-D/Def2-TZVP//Def2-SV(P) level.

Fig. 1
Fig. 1 Crossing of two diabatic potential energy surfaces and the relations among V 12 , l, DG 0 , DG ‡ , and DG ‡ na in the rate constant expressions.

Fig. 2
Fig.2The entire workflow for obtaining an electron transfer rate constant by EDEEL and the calculation programs used in each procedure.

Fig. 5
Fig. 5 Changes in mean coordination bond length.

Table 1
Spin multiplicity used in the calculations, experimental rate constants k expt of the thirteen electron self-exchange reactions of transition metal complexes, and their data sources Values for monomer and SX are shown outside and inside the parentheses, respectively. a

Table 2
log 10 k calc , Gibbs energies of activation DG ‡ , reorganization energy l, diabatic coupling V 12 , electronic transmission coefficient k el , and donor-acceptor distance r DA by EDEEL at the UuB97X-D/Def2-TZVP//Def2-SV(P) level Redox couple log 10 k calc DG ‡ /eV l/eV V 12 /eV k el r © 2023 The Author(s).Published by the Royal Society of Chemistry RSC Adv., 2023, 13, 32097-32103 | 32099 Paper RSC Advances (4) ASM analysis at the geometry of the largest k ET (=k calc ) obtained in step 3 (optional).