The Marcus dimension: identifying the nuclear coordinate for electron transfer from ab initio calculations

The Marcus model forms the foundation for all modern discussion of electron transfer (ET). In this model, ET results in a change in diabatic potential energy surfaces, separated along an ET nuclear coordinate. This coordinate accounts for all nuclear motion that promotes electron transfer. It is usually assumed to be dominated by a collective asymmetric vibrational motion of the redox sites involved in the ET. However, this coordinate is rarely quantitatively specified. Instead, it remains a nebulous concept, rather than a tool for gaining true insight into the ET pathway. Herein, we describe an ab initio approach for quantifying the ET coordinate and demonstrate it for a series of dinitroradical anions. Using sampling methods at finite temperature combined with density functional theory calculations, we find that the electron transfer can be followed using the energy separation between potential energy surfaces and the extent of electron localization. The precise nuclear motion that leads to electron transfer is then obtained as a linear combination of normal modes. Once the coordinate is identified, we find that evolution along it results in a change in diabatic state and optical excitation energy, as predicted by the Marcus model. Thus, we conclude that a single dimension of the electron transfer described in Marcus–Hush theory can be described as a well-defined nuclear motion. Importantly, our approach allows the separation of the intrinsic electron transfer coordinate from other structural relaxations and environmental influences. Furthermore, the barrier separating the adiabatic minima was found to be sufficiently thin to enable heavy-atom tunneling in the ET process.


IVCT spectra
To connect the electronic structures with experimental information, the IVCT bands were reconstructed from the nuclear ensemble method as shown in Fig. 1.The spectra were computed according to: 1

ϵ(E) =
N A e 2 4 ln (10) where N A is the Avogadro constant, e is the elementary charge, c is the speed of light, m e is the mass of the electron, ϵ 0 is the permittivity of vacuum, δ is the Dirac delta function, g is the normalized Gaussian line shape with full width at half maximum 0.5 eV, f ij and E ij are oscillator strength and excitation energy for i-th electronic state and j-th geometry in the ensemble, respectively.
We compare the results for the electronic structures computed with the local hybrid density functional LH20t 2 with those obtained with the BLYP35 3 functional proposed previously that is based on a global hybrid functional (BHLYP 4 with 35 % exact exchange).The BLYP35 calculations were run with the ORCA 5 suite of programs employing the def2-TZVP 6 basis set for carbon, oxygen and nitrogen atoms and def2-SVP for hydrogens.Solvation with ACN was modelled implicitly using a polarizable continuum model.The resolution of identity approximation was employed to approximate Coulomb and exchange integrals. 7Convergence criterion for the self-consistent field procedure was set to 10 −8 E h .
Both electronic structure methods show a reasonable performance in predicting the IVCT bands.We note in particular that the relative intensities for the Class II and Class III cases are reproduced well; no scaling factors are applied in the comparison of experiment and theory.The characteristic spikes in the experimental spectra of Class III compounds are not reproduced in the computational spectra.A possible explanation might be found in the approximated treatment of vibronic transitions in the nuclear ensemble method (the Duschinsky effect and the frequency shift between ground and excited state are neglected).
Of the two density functionals, LH20t predicts the band shapes better for all systems studied and yields slightly better intensities for the Class III cases (p-DNB •− and 2,6-DNN •− ).With BLYP35, the predicted bands are too wide, which might be due to a slightly too localized electronic structure or a not ideally described vibrational structure which in turn affects the sampling of geometries.3 Marcus dimensions

Experiment
Vacuum Acetonitrile 4 Conventional approach: linear interpolation of Cartesian coordinates In this Section we show the electron transfer dimension obtained from a linear interpolation of Cartesian coordinates (LICC) as has been used in the literature so far. 10This method requires two well-defined structures in Cartesian coordinates, usually an adiabatic minimum and a totally symmetric structure representing the transition state.The reaction coordinate is obtained by linear interpolation between those structures.
interpolation minimum to symmetric molecule interpolation between minima We note that this approach cannot be applied to Class III systems where only one well-defined structure exists.LICC cannot identify the ET dimension in Class III systems and could not be then used to decouple intramolecular and solvent motion.In addition, mapping the potential curves obtained by a LICC approach for a Class III system onto the diabatic (i.e.uncoupled) states of the underlying Marcus theory is thus precluded, and hence the progression of coupling strength, nuclear motion or force constants in a series of Class III systems or for borderline Class II/Class III cases could not be evaluated.Regardless of these conceptual limitations, we tested the LICC approach on m-DNB •− in ACN.The LICC scans were performed between the adiabatic minimum and the totally symmetric molecule (Fig. 4a), and between the two adiabatic minima (Fig. 4d), since these are the only well defined structures of a class II systems.
The ET dimensions obtained by LICC (Fig. 4a,d) are not anti-symmetric motions: a twisting motion of the oxygen atoms of one nitro group dominates the overall motion, and some hydrogen and carbon atoms of the aromatic ring have non-negligible contributions.In the adiabatic minimum, where the unpaired electron is localized on one side of the molecule (in this case, the nitro group on the left-hand side), the other nitro group is further away from the aromatic bridging unit and is slightly distorted, which will decrease the potential coupling with the bridging unit.A twisting motion is expected to be an important part of the ET only in cases of very large distortions. 11For a nitro group perpendicular to the aromatic plane, the conjugation with the rest of the molecule breaks, and the potential coupling would thus tend to zero, making the ET adiabatically forbidden.
Unrelaxed scans along the two types of LICC coordinate (Fig. 4b,e) produce double-well potentials that change smoothly everywhere.For the LICC scan between the adiabatic minimum and the totally symmetric molecule, the minima do not have the same energy: the second minimum is destabilized by 741 cm −1 , which is more than half of the ET barrier.The potential obtained by scanning between the two adiabatic minima is symmetric and the barrier height agrees well with the value obtained by simulating experimentally measured EPR spectra (1000 cm −1 ). 8 The electron position changes upon overcoming the potential barrier in both cases ), with a sharper transition in the scan between the two adiabatic minima.
Clearly, the results of the LICC approach depend on the choice of the reference structures, and therefore it is unclear whether the ET coordinate obtained with a solvent model would lead to meaningful potentials when transferred into vacuum.Furthermore, a strong limitation of LICC is that it can only be applied to Class II systems, since for Class III systems no second geometric reference structure besides the adiabatic minimum exists.We will show in the following sections that by exploiting properties postulated in the Marcus model when interpreting the ab initio calculations, we arrive at a chemically intuitive dimension that drives electron transfer in Class II and leads to electron localization in Class III systems.

Multideterminant calculations
Seeking an explanation for the cusp observed in the scanned potential energy curve for m-DNB •− using DFT, we repeated the scan with state-averaged complete active space self-consistent field (SA-CASSCF) calculations 12 and a subsequent treatment with the N -electron valence state perturbation theory (SA-CASSCF/NEVPT2). 13e calculations were run with the ORCA 5 suite of programs employing the def2-TZVP 6 basis set and implicit solvation using a polarizable continuum model with acetonitrile as the solvent modeled.The resolution of identity approximation (RIJK) was employed to approximate Coulomb and exchange integrals. 7Two electronic states were considered in the state-averaged calculations: the electronic ground and first excited state at the adiabatic minimum.The active space contains nine electrons in eight orbitals, i.e. (9, 8), which are selected from the conjugated π-system comprising the two nitro groups and the benzene ring.The remaining four π bonding orbitals were significantly lower in energy and not included in the active space.The orbital shapes change during the scan, see Fig. 7.
The scan along the Marcus dimension with CASSCF and CASSCF/NEVPT2 is shown in Fig. 5.The CASSCF approach provides unreasonably high excitation energies and the minima in the ground state potential curve are significantly shifted apart from each other.These are consequences of the lack of dynamical correlation, resulting in an overstabilized localized electronic structure.The CASSCF method results in a blue-shift of the IVCT band by more than 4000 cm −1 in comparison with the experiment.Accounting for dynamical correlation with the NEVPT2 approach (Fig. 5 right) leads to a substantial improvement of the potential energy curves.The double well potential in the ground state is correctly reproduced and the excitation energies match the position of experimentally measured IVCT band.The perturbation theory correction is, however, very large which may question the validity of the perturbative treatment in this case and calls for substantially enlarging the active space.This was tested on one structure by adding three low-lying orbitals from the conjugated π-system to the active space, i.e. (15, 11), which resulted in very similar energies and energy corrections to the electronic states.
Since we are showing a proof of concept, enlarging the active space to achieve quantitative agreement is not the principal aim of this work.It can be reasonably expected that adding more orbitals besides the conjugated π-system considered here would enlarge the active space to the realms of a DMRG-SCF/NEVPT2 calculation.The differences in the description of electronic structure with DFT and a multiconfigurational method are displayed in Fig. 6.The CASSCF/NEVPT2 method provides slightly lower excitation energies and the resulting IVCT band is thus expected to be red-shifted with respect to the experiment.The adiabatic minima of the CASSCF/NEVTP2 calculations are closer to each other than for the calculations with the LH20t density functional.
Most importantly, the scan with the wavefunction method results in a smooth shape of the potential barrier instead of the presumably unrealistic cusp.In fact, the barrier exhibits a feature exactly opposite to the cusp, a dip.This behaviour was already observed for other MV systems. 14The dip can be explained as a stabilization of the wavefunction by delocalization, see Fig. 7.If the electronic structures adjacent to the dip are not as delocalized, they are expected to have higher energies.
The barrier height obtained from the CASSCF/NEVPT2 scan is, however, unreasonably small (290 cm −1 ) which would result in an ET rate that is two orders of magnitude faster than the experimental value.We expect that a faithful description of the Marcus dimension can be achieved with a strongly correlated multireference method (e.g.MRCI or DMRG-SCF/NEVPT2 with a large active space).We suggest that DFT is prone to fail only in proximity of the very top of the barrier.While CASSCF/NEVPT2 delivers a smooth potential and a reasonable shape of the barrier, the absolute values obtained with the current active space are questionable as reflected in the unreasonably large corrections from the perturbative treatment. Energy

Properties of the Marcus dimension
In practice, the Marcus dimension is described as a normalized 3N -vector of Cartesian displacements q.The mass of the fictitious particle moving along the Marcus dimension was obtained by weighting the atomic masses with the squared values of their displacement coefficients: In order to estimate the classical ET rate, a frequency along the Marcus dimension is needed.To obtain this frequency, we can make a use of the fact that the mass of the phonon in the Marcus dimension and the Hessian matrix at the adiabatic minimum are known.
As a first step, the harmonic force constant for the Marcus dimension needs to be computed.To this end, we are going to need the following 3N × 3N matrices: Ω as the diagonal matrix of harmonic frequencies, M as the diagonal matrix of atomic masses, and C as the matrix of eigenvectors of the Hessian (normal modes) collected as columns.We express the Marcus dimension as a combination of normal coordinates, a = C T M 1 2 q, which should be subsequently normalized.The force constant will be computed in a similar manner as the mass in Eq. 2, i.e., combining force constants of normal modes.The harmonic frequency can then be computed as ω = k µ .Combining all of the above results in the following equation for the frequency along the Marcus dimension: In the main text of the paper, we described two possibilities for the parameterization of the Marcus model from the ab initio scan.The procedure is described in more detail below.
We work with the two diabatic surfaces G a = f x 2 and G b = f (x − d) 2 .Here, x denotes the progress along the ET coordinate, f the force constant and d the separation of the diabatic states.Under the action of potential coupling V ab , we obtain the two adiabatic states: The properties of the adiabatic states and their spectroscopic implications are already discussed elsewhere in the literature. 15,16Here, we make use of some important properties described in the following.The excitation energy at the adiabatic minimum, also know as the reorganization energy, is λ = f d 2 .If the potential coupling 2V ab < λ, we will obtain a Class II system with a double well potential in the ground state.The position of the barrier will be d/2 and the positions of the adiabatic minima 1 2 . The excitation energy at the top of the barrier, i.e. position d/2, is 2V ab .The barrier height is .
From the ab initio scan we can extract the following quantities: the excitation energy at the adiabatic minimum E Ex. min., the excitation energy at the top of the barrier E Ex. barr., the distance between the adiabatic minimum and the top of the barrier R, and the barrier height ∆G.In the main text, we showed the two possibilities (A and B) for relating the obtained quantities to the Marcus model.
In parameterization A, we expressed the quantities E Ex.
barr. , R and ∆G using the adiabatic states G 1 and G 2 from the Marcus model, resulting in three equations: In parameterization B, we expressed the quantities E Ex. barr., E Ex. min.and R using the adiabatic states G 1 and G 2 from the Marcus model, again resulting in three equations: Either set of Eqs. 5 or 6 can then be solved for f , d and V ab which are the parameters of a symmetrical Marcus model.The solution can be obtained numerically or using any algebraic software of choice.

Heavy atom tunneling
To investigate the relevance of heavy atom tunneling, we computed the transmission coefficient using the semiclassical WKB (Wentzel-Kramers-Brillouin) approximation.The probability of tunneling through the barrier depends on the classical action accumulated between the so-called turning points: Here, T is the transmission coefficient, µ is the mass of the particle, V (x) is the potential energy, E is the total energy of the incident particle and x 1,2 are the turning points defined as solutions of equation E = V (x).
In our approach, we assume a non-stationary behaviour in the double-well potential, i.e., a particle with energy equal to the harmonic zero point energy (ZPE) is scattered by the potential barrier, see Fig. 9.The turning points x 1,2 are obtained as solutions of ZPE = V (x).Note that the region 1/2 is imaginary).As a potential V (x), the analytical potential from parameterization A, presented in the previous section, was used.The ZPE was then calculated as: , where µ is the phonon mass from Table 1.The resulting transmission coefficient are presented in Table 2.The purpose of this section is to provide an explanation why the Marcus dimension of a class II systems is difficult to obtain by fitting to the electron position.We start with the Marcus-Hush theory, the diabatic states a and b are separated on dimensionless coordinate by 1. Their energies are: The coupling of the diabatic states V ab will lead to the adiabatic wavefunctions ψ 1 and ψ 2 .The unnormalized coefficients c a and c b of the ground state (ψ 1 ) can be written as: It is clear from Eqs. 9 that the admixture of the diabatic state in the adiabatic ground state depends on the progression along the ET coordinate X but also non-trivially on the coupling element V ab .To elucidate more this dependence the squared absolute values of the normalized coefficients c a and c b are plotted in Fig. 10

Figure 3 :
Figure 3: The electron transfer dimensions resulting from the multi-component linear fit as introduced in the main text for vacuum (left) and acetonitrile (right).In vacuum, all systems behave like a Class III system, and therefore the fit used the electron position.In acetonitrile, m-DNB •− and 2,7-DNN •− belong to Class II, and therefore the fit used the excitation energy to the lowest doublet state, whereas p-DNB •− and 2,6-DNN •− belong to Class III and the fit used the electron position.

Figure 4 :
Figure 4: ET coordinates from the LICC approach and corresponding potential energy curves obtained from a scan along them.a) ET dimension obtained by LICC between the optimized geometry (adiabatic minimum) and the totally symmetrical one.b) Unrelaxed scan along the LICC dimension, with labels indicating the values of the reorganization energy, the electronic coupling, the height of the barrier, and the offset of the second minimum.c) Change of the electron position along the scanned coordinate.Panels d, e, f depict the same quantities for the ET coordinate obtained with LICC between the two adiabatic minima.

Figure 6 :
Figure 6: Comparison of state-averaged CASSCF/NEVPT2 and LH20t results for the scan along the Marcus dimension of m-DNB •− .

Figure 7 :Figure 8 :
Figure 7: Active space orbitals for the SA-CASSCF/NEVPT2 calculations.The orbitals are shown for three positions in the scan in Fig. 5.

Figure 9 : 13 10
Figure 9: Semi-classical estimation of nuclear tunneling through the potential barrier.V (x) represents the parameterized Marcus model.The energy of the incident particle is taken from harmonic ZPE at the adiabatic minimum (dashed horizontal line).Turning points x1,2 represent the boundaries of the classically forbidden region.
as a function of X for two different values of V ab .

Figure 10 :
Figure 10: Square absolute value of the diabatic state coefficient in the adiabatic ground state as a function of dimension less coordinate.Left panel depicts a class II situation where λ > 2V ab , right panel depicts a class III situation where λ ≤ 2V ab .

Table 1 :
Phonon masses and harmonic frequencies for the Marcus dimensions in ACN.