A general statistical mechanical approach for modeling redox thermodynamics: the reaction and reorganization free energies

By using rigorous statistical mechanical derivations we obtain a general theoretical model providing the thermodynamics of redox processes, with a focus on the reaction and reorganization free energies and on the relationship between these key thermodynamic quantities. In particular, we define two distinct reorganization free energies, l P and l R , for the reactants (R) A products (P) reaction and for the inverse process, respectively. We first derive in principle exact relationships, then gradually introduce different levels of approximation to obtain more and more simplified, though less general, working equations. The results of the calculation of thermodynamic properties for two model systems are then used to compare general and more approximated expressions and critically assess their applicability to the description of redox processes. Finally, we obtain specific relationships that can be used as a diagnostic tool to test the actual reliability of the assumption of Gaussian fluctuations, a priori accepted within Marcus theory, for any redox system under investigation. For both benchmark molecules studied in the present paper, the Gaussian approximation turns out to be inappropriate to describe the redox thermodynamics.


Introduction
Electron transfer (ET) processes play a fundamental role in energy transduction pathways of living cells. 1,24][5] Thorough investigation of the molecular determinants to ET would lead not only to elucidation of several biological processes but also to severe improvements in the design of the aforementioned devices. 6][9][10][11][12][13][14][15][16][17][18][19][20] The usual framework for the description of ET reactions is provided by the semiclassical Marcus theory, 21 furnishing information on the kinetics and thermodynamics of redox processes, which is actually based on the severe assumption of Gaussian fluctuations of the energy change due to the electron transition (transition energy). 22,23In this paper we describe a new, general approach for modeling the thermodynamics of redox processes without assuming such Gaussian fluctuations.In particular, our attention is focused on two key quantities essential in ET reactions and, more in general, in redox processes: the reaction free energy and the reorganization free energy, where the latter is typically defined as the free energy change necessary to relax the product state right after the vertical electron transition.

General definitions and derivations
We consider a system composed by a single redox center (i.e. the solute, defined by a single acceptor or donor molecule or by the donor-acceptor complex) embedded in a large amount of solvent molecules in order to describe the statistical mechanical behavior of real systems at infinite solute dilution.We define the chemical state of the system via the solute reactant (R) or product (P) redox state: the chemical states characterized at each solute and environment configuration by the solute electronic state corresponding to the proper R or P donor/acceptor charge and properties, with the environment in its electronic ground state.Therefore, the R A P chemical reaction may correspond to an electron transfer between the donor and the acceptor as well as a pure oxidation or reduction of a single chemical group.We define the system energy when the redox center is in the lth vibrational state of its reactant or product electronic condition with the environment in the jth vibrational state of its electronic ground state as R,l, j $ R (j,p j ) + E v,R,l + E v,e, j and P,l, j $ P (j,p j ) + E v,P,l + E v,e, j , respectively, where j, p j are the semi-classical coordinates and conjugated momenta defining the system phase space, R and P include the system electronic energy and classical kinetic energy with the redox center in the R or P electronic condition, E v,R,l , E v,P,l are the lth vibrational state energies of the redox center in the R and P condition, respectively, and E v,e, j is the jth vibrational state energy of the environment electronic ground state which we assume as independent of the solute R and P electronic conditions.Hence, the corresponding energy change due to the R A P transition at each phase space position (transition energy 24 ) is DH l, j ~H P,l, j {H R,l, j %U P (j,p j ){U R (j,p j )zE v,P,l {E v,R,l DU (j,p j )~U P (j,p j ){U R (j,p j ) (2) where we have assumed that the vibrational energies are independent of the phase space position (such approximation, which is not really appropriate for every possible phase space region, becomes very accurate when we deal with statistical mechanical calculations 25 as in the present work, see the appendix).Note that D can be in general accurately approximated by the solute R A P electronic energy change at each phase space position, 24 as it follows assuming that the environment electronic energy just like the classical kinetic energy is independent of the solute R and P electronic conditions, thus neglecting such weak higher order perturbation effects. 24Within such usual approximations the canonical partition functions of the reactant Q R and product Q P redox states are where dC= Hdjdp j /h D is the quantum state differential with h the Planck's constant, D the number of the semi-classical j coordinates and H a constant providing the quantum corrections needed to express the partition function via a phase space integral, e.g. the correction for the permutations of identical particles.Q v,R , Q v,P are the solute vibrational partition functions for the R and P electronic condition, respectively, Q v,e is the environment electronic ground state vibrational partition function and l/b = kT with k the Boltzmann constant.It must be remarked that the phase space integrals in eqn (4) and ( 5) are in principle meant over all the phase space regions where the R and P electronic conditions correspond to specific solute electronic Hamiltonian eigenstates.Therefore, when considering as solute the donor-acceptor complex, such integrals cannot involve the phase space tiny regions where the solute electronic eigenstates, corresponding to the R or P condition elsewhere, are characterized by a mixed R and P electronic distribution, thus defining the transition regions where the electron transfer occurs. 24However, such transition regions where the diabatic states (i.e. the electronic eigenstates obtained fixing the donor/ acceptor charge to either the R or P condition) do not coincide anymore with the adiabatic states (i.e. the true electronic eigenstates as obtained without constraining the donor/ acceptor charge) are typically negligible when calculating phase space integrals, as in most of the phase space positions adiabatic and diabatic states are virtually identical (for a detailed discussion on this matter for electron transfer reactions see ref. 24).Hence the phase space integrals in eqn ( 4) and ( 5) can be considered over virtually the whole available phase space disregarding, when dealing with the donoracceptor complex, the presence of transition regions (note that for a single acceptor or donor solute redox center by definition these integrals are exactly defined over the whole phase space).
For the reactant ensemble we can calculate the Landau free energy A R as a function of D via and analogously for the product state the corresponding Landau free energy A P is which leads to and hence It must be noted that eqn (8) necessarily implies with n ¢ 2.Moreover, expansion of the Landau free energies around their minima provides where D 0 P and D 0 R are the transition energy values corresponding to the A P and A R minima, respectively.From the last equations it follows that if A R were exactly defined by a purely quadratic function, i.e. h 2 A R /hD 2 = c for whatever transition energy value and hence then indicating that in such a particular case the two Landau free energy curves would only differ for the minimum position.
From the definition of the partition function and the Landau free energy, we can express the probability density as a function of D for the reactant r R (D ) and product r P (D ) state respectively providing Therefore, from where DA = 2kTlnQ P /Q R stands for the R A P Helmholtz free energy difference (i.e.DA = A P 2 A R ) which exactly corresponds to the solute R A P standard chemical potential change.Note that the R and P Helmholtz free energies are defined involving all the vibrational states of the R and P electronic conditions and not just a single R or P vibronic state.The exact eqn ( 8) and ( 17)- (20) provide the fundamental relations between the Landau free energies and the probability densities as well as between the reactant and product states (note that from eqn ( 17) and ( 18) quadratic Landau free energies provide Gaussian probability distributions).Moreover, in eqn (20) e 2b(D 2DA) Q v,P /Q v,R is the operator transforming the reactant probability density into the product one.
At this stage we need to introduce the concept of the Landau relaxation free energy involved in our general definition of the reorganization free energy for the R A P reaction, considering that both the reactant and product chemical states are characterized by the equilibrium distributions given by the corresponding Landau free energy curves.For each possible initial reactant and final product position (i.e.reactant and product local states) along the A R and A P surfaces we may define the corresponding Landau free energy change as composed by a first vertical Landau free energy variation followed by the relaxation along the A P curve to reach the final product local state (see Fig. 1).Therefore, we should properly define the product reorganization free energy l P from the non-equilibrium product distribution provided by the vertical transitions to the final equilibrium product distribution, as the average Landau relaxation free energy provided by (see eqn ( 8)) where the R and P angle brackets subscripts stand for averaging in the R and P ensemble, respectively.Similarly, for the P A R reaction we can express the reactant reorganization free energy l R as and hence we also obtain We can proceed further by expressing the Helmholtz free energy change via the Landau free energies (see eqn ( 6) and ( 7)) g~Ð e {b(A P {SA P T P ) dDU which can be used to substitute SA P T P 2 SA R T R in eqn ( 21) and ( 22) to give where the reaction free energy DA can be obtained by the general expression of the Helmholtz free energy change Subtracting eqn (27) from eqn (26) we obtain or rearranging to express DA which explicitly relate the reorganization free energies to the reaction free energy.
All the statistical mechanical relations derived so far can be considered, in principle, exact relations which are then completely general and always accurate to treat any redox system.
In the next subsections, we will introduce progressive levels of simplifications based on the use of an increasing number of approximations.

First simplification: vibrational partition function invariance and local quadratic approximation
Eqn (28) can be simplified by assuming Q v,R $ Q v,P which then provides Note that within the approximation of invariant vibrational partition function, eqn (31) coincides with the expression of the free energy change we gave in a previous paper 24 for the electron transfer process defined by the transition from a R vibronic state to a P vibronic state, when considering the R and P vibronic states as involving the ground vibrational states.The approximation of assuming the vibrational partition functions of the R and P chemical states as virtually identical is rather general and accurate as it follows from the fact that the vibrational energy variations due to the redox state transition are typically negligible with respect to the corresponding electronic energy changes.In this subsection, together with the approximation just described, we will also use the local quadratic behaviour of the Landau free energies close to their minima.In fact, assuming that SA P T P , SA R T R and g can be almost exactly evaluated by using relatively narrow D ranges around the Landau free energy minima and considering that A R and A P close to their minima are almost indistinguishable from quadratic curves (i.e. the corresponding probability densities close to their maxima are almost indistinguishable from Gaussian distributions), we obtain within such narrow ranges A P %A P (DU 0 P )z and hence we can write, assuming accurate the quadratic approximation to evaluate the second D moment,  g% By using Q v,P $ Q v,R and eqn (40) and (41) into eqn ( 24), ( 26), ( 27) and (30), we obtain A P (DU 0 P ){A R (DU 0 R )%DAzkT ln where in eqn ( 42)-( 44) DA is given by eqn (31).
It must be now remarked that eqn (40)-( 45) can be still considered as general accurate approximations as they only require, beyond Q v,P $ Q v,R , that A R and A P can be properly described by quadratic curves within relatively narrow D ranges around the Landau free energy minima.

Second simplification: indistinguishable reorganization free energies
If we further assume l P $ l R , i.e. we assume that the mean Landau relaxation free energies are very close, that is we can write from eqn (45) Eqn ( 46) and (47) stemming from l P $ l R , although might be reasonably accurate in many redox systems, must be considered as less general approximations than the ones described in the previous subsection, applicable only to a subset of redox reactions.It should be also noted that all the approximated relations derived so far were obtained without assuming a global quadratic shape of the Landau free energies (i.e.Gaussian shape of the probability densities) which hence in our general approach, differently from Marcus theory, is not at all a required assumption.In the next, final, theory subsection we will extend the quadratic approximation to the whole transition energy interval relevant in the redox process, hence assuming a global quadratic behaviour of both the Landau free energies.

Third simplification: global quadratic approximation
In case of quadratic Landau free energies, that is the Landau free energies are properly described by quadratic curves in the whole relevant D range, as required in Marcus theory, we would have (see eqn (15) and ( 16)) and from the definition of the moments generating function of Gaussian distributions 25 (54) providing when using also eqn ( 31), ( 42)-( 44) (55) Therefore, from leading to (when used into eqn (55) and ( 56)) providing the usual relations utilized in Marcus theory.However, such last equations are not at all necessarily due to the quadratic behavior of the Landau free energy as they correspond to eqn (46) and (47) (for the latter when considering the special case s R = s P ) which were obtained by using much less demanding approximations than assuming virtually exact Gaussian distributions for D in the R and P chemical states.In fact, in order to address the accuracy of the quadratic approximation we should consider eqn (55) and the relation, following from eqn (48)-( 52), which, when fulfilled, are specifically indicating that a redox system is properly described by the global quadratic approximation, i.e. the R (P) probability density is properly modeled by a Gaussian distribution in the whole relevant D range.In Fig. 2 we provide a simple scheme summarizing the main simplifications-approximations described in the Theory section.

Results and discussion
We will now focus on the application of the previously derived theoretical relationships to the calculation of the reaction and reorganization free energies of two redox systems we recently studied: namely, the relatively small organic molecule 2FSN 24 and a protein, more precisely yeast cytochrome c (cytc hereafter). 26For both systems, the transition energy was calculated by using a mixed quantum mechanics/molecular dynamics theoretical approach, based on the Perturbed Matrix Method 27,28 (PMM).As for others QM/MM methods, a portion of the system is treated explicitly at the electronic level (the quantum centre, QC), while the rest of the system is described at a classical atomistic level providing the electrostatic perturbation on the electronic states of the QC.The first step of this method is the calculation of an orthonormal set of unpeturbed electronic Hamiltonian H ˜0 eigenfunctions of the QC.Then, for each configuration generated by the explicit solvent MD simulation, a perturbed electronic Hamiltonian matrix H ˜is generated as: where q T is the total charge of the QC, V and E are the perturbing electric potential and field, respectively, exerted by the environment on the QC (typically obtained by the environment atomic charge distribution and evaluated in the QC center of mass), I ˜is the identity matrix, DV is a short range potential energy term which approximates all the terms higher than the dipolar one, m ^is the dipole operator and W 0 m , W 0 m the m and n QC unperturbed electronic eigenstates.The diagonalization of H ˜at each MD frame yields a set of eigenvectors and eigenvalues which represent the perturbed electronic eigenstates and energies of the QC.The detailed description of the molecular dynamics (MD) simulations, as well as the ab initio quantum chemical calculations and the mixed quantum mechanics/molecular dynamics theoretical approach to evaluate the perturbed quantum energy and properties 27,28 for these two systems, can be found elsewhere. 24,26or both systems the estimates of the reaction free energy DA using the most general relation given by eqn (31), based on the approximation Q v,P $ Q v,R , are listed in Table 1 where they are compared with the corresponding experimental values.A very good agreement, within the estimated noise of a few kJ mol 21 , can be observed confirming that the only assumption made in that equation, i.e. invariant vibrational partition function, is rather general and accurate.As shown in the Theory section, assuming indistinguishable reorganization free energies allows for the calculation of DA using the further simplified expression given by eqn (47).The redox free energies for both 2FSN and cytc as obtained via such a simplified relation, are also reported in Table 1 showing that the results of the two different levels of approximation, i.e. eqn (31) and (47), are so close to be indiscernible within the noise (a few kJ mol 21 ).This finding implies that the expression provided in eqn ( 47) is likely to be often reliable and accurate, yielding a computationally efficient shortcut for the calculation of DA.Note that for the systems considered in this paper the reorganization free energies for the R A P (l P ) and P A R (l R ) transitions, reported in Table 2, are very close, with the relative errors between l R , l P and their average value l of only 0.3% and 4.2% for 2FSN and cytc, respectively.
The assessment of the trustworthiness of a priori assuming Gaussian fluctuations of the transition energy (i.e.global quadratic behavior of the Landau free energy) for any redox system is another key issue of the present work.In the Theory section we have derived specific relations which are diagnostic Fig. 2 Simplified scheme of the approximation levels derived in the Theory section.
Table 1 Calculated and experimental redox free energies DA of solvated 2FSN 24 and yeast cytochrome c. 26 All data are expressed in kJ mol 21 .Calculated DA have been determined at two levels of approximation: by using the most general relation given by eqn (31) or by means of the approximation given by eqn (47) DA a DA b DA exp c of the Gaussian behavior: among them, eqn (61) provides a proper tool for addressing the reliability of the Gaussian approximation, which can be accepted when SD T R 2 SD T P = 2l $ s 2 /(kT); if s 2 /(kT) differs significantly from SD T R 2 SD T P , then the system under investigation cannot be properly described by the global quadratic approximation.As it is shown in Table 3, for both redox systems we considered in this paper such an approximation is largely inaccurate.It is important to stress that the Gaussian approximation must be rejected for both systems although s R and s P are fairly close (see Table 2), explicitly showing that the accuracy of the approximation given by eqn (47) is indeed not necessarily associated to the presence of Gaussian fluctuations of the transition energy.

Appendix
The most general expression of the canonical partition function for a system as defined by a single solute molecule embedded into a large amount of solvent molecules, i.e. a model system representative of solvated solute molecules at infinite dilution, is given for the R and P chemical states by Q R ~X n ð e {b½eR(j)zE R,n (j)zK (j,p j ) dC Q P ~X n ð e {b½eP(j)zE P,n (j)zK (j,p j ) dC where the sum runs over the vibrational states, e R (j),e P (j) and E R,n (j),E P,n (j) are the system electronic and vibrational energies corresponding to the solute R and P electronic conditions and (j, p j ) is the system classical kinetic energy.By introducing a phase space independent reference vibrational energy for each nth vibrational state of the R (E ref,R,n ) and P (E ref,P,n ) electronic condition, we obtain In liquid state systems (the systems we consider in this paper) infrared (IR) spectra of vibrational modes provide widths of at most a few hundreds of Joules per mole with almost temperature independent maxima along an isochore, thus indicating that compared to thermal energy only small deviations of the vibrational energies as a function of the phase space position are allowed and that the average vibrational energies can be considered as virtually constant within the whole temperature range of interest.Therefore, in typical liquid state systems, once properly defined the reference vibrational energies, at virtually any liquid state temperature along the isochore considered to define the canonical partition functions, we can reasonably assume ~X l e {bE v,P,l the R and P solute vibrational partition functions and Q v,e = S j e 2bE v,e, j the environment vibrational partition function.

Fig. 1
Fig.1Example of a R A P transition path, composed by a vertical transition (red) from the reactant to the product Landau free energy curve, followed by a relaxation along the product curve (green).

s 2 P~SDU T P z b 2 s 2
53) kT ln Se bDU T P ~SDU T P z b 2 Q R ~X n e {bE ref,R,n ð e {b(eRzK ) e {b(E R,n {E ref,R,n ) dC Q P ~X n e {bE ref,P,nð e {b(ePzK ) e {b(E P,n {E ref,P,n ) dC EP,n2Eref,P,n) $ 1 2 b(E P,n 2 E ref,P,n ) {bE ref,P,n ð e {b(ePzK ) ½1{b(E P,n {E ref,P,n )dC %Q v,P Q v,e n e {bE ref,R,n ð e {b(eRzK ) ½1{b(E R,n {E ref,R,n )dC %Q v,R Q v,e ð e {bU R dC Q P % X n e n e {bE ref,R,n %Q v,R Q v,e X n e {bE ref,P,n %Q v,P Q v,e with Q v,R ~X l e {bE v,R,l Q v,P