^{1}, Marialore Sulpizi

^{1}and Michiel Sprik

^{1,a)}

### Abstract

The density functional theory based molecular dynamics (DFTMD) method for the computation of redoxfree energies presented in previous publications and the more recent modification for computation of acidity constants are reviewed. The method uses a half reaction scheme based on reversible insertion/removal of electrons and protons. The proton insertion is assisted by restraining potentials acting as chaperones. The procedure for relating the calculated deprotonation free energies to Brønsted acidities and the oxidationfree energies to electrode potentials with respect to the normal hydrogen electrode is discussed in some detail. The method is validated in an application to the reduction of aqueous 1,4-benzoquinone. The conversion of hydroquinone to quinone can take place via a number of alternative pathways consisting of combinations of acid dissociations, oxidations, or dehydrogenations. The free energy changes of all elementary steps (ten in total) are computed. The accuracy of the calculations is assessed by comparing the energies of different pathways for the same reaction (Hess’s law) and by comparison to experiment. This two-sided test enables us to separate the errors related with the restrictions on length and time scales accessible to DFTMD from the errors introduced by the DFT approximation. It is found that the DFT approximation is the main source of error for oxidationfree energies.

J.C. thanks Ali Shah and Chris Adriaanse for kindly offering assistance at the beginning of the project. Ruth Lynden-Bell is acknowledged for a number of helpful discussions. M.S. and J.C. are grateful for EPSRC for financial support. The calculations for this work have been performed using an allocation of computer time on HECToR, the U.K.’s high-end computing resource funded by the Research Councils, as part of a grant to the UKCP consortium.

I. INTRODUCTION

II. METHODS

A. Free energies for half reactions from vertical energy gaps

B. Restraining potentials and chemical species definition

C. Finite system size and periodic boundary effects

D. and redox potentials from half reaction energies

E. Corrections for restraints and zero point motion

F. Computational setup

III. RESULTS AND DISCUSSION

A. Vertical energy gaps and free energy changes

B. Solvent reorganization and hydration structure

C. Comparison to experiment

D. Summary of error analysis and evaluation of DFT

IV. CONCLUSION AND OUTLOOK

## Figures

Decomposition of the reduction of quinone to hydroquinone in reduction (horizontal arrows) and protonation half reactions (vertical arrows). The quinone, semiquinone, and hydroquinone are indicated by Q, HQ, and , respectively. Diagonal arrows correspond to (de)hydrogenation (proton coupled electron transfer). Note that the usual dot indicating that a species is open shell (a radical) has been suppressed.

Decomposition of the reduction of quinone to hydroquinone in reduction (horizontal arrows) and protonation half reactions (vertical arrows). The quinone, semiquinone, and hydroquinone are indicated by Q, HQ, and , respectively. Diagonal arrows correspond to (de)hydrogenation (proton coupled electron transfer). Note that the usual dot indicating that a species is open shell (a radical) has been suppressed.

Accumulative averages of vertical energy gaps as a function of MD time for the four oxidation reactions of Fig. 1: (a) ; (b) ; (c) ; and (d) . The three curves for each reaction are the result for three different PESs as defined by the coupling parameter [see Eq. (9)]. In the order of decreasing value of the energy gap the coupling parameter values are (reactant), (mixed state), and (product). Each curve is followed by the overall time average of the vertical energy gap in eV (see also Table III).

Accumulative averages of vertical energy gaps as a function of MD time for the four oxidation reactions of Fig. 1: (a) ; (b) ; (c) ; and (d) . The three curves for each reaction are the result for three different PESs as defined by the coupling parameter [see Eq. (9)]. In the order of decreasing value of the energy gap the coupling parameter values are (reactant), (mixed state), and (product). Each curve is followed by the overall time average of the vertical energy gap in eV (see also Table III).

Accumulative averages of vertical energy gaps as a function of MD time for the deprotonation and dehydrogenation reactions of Fig. 1. Labeling is continued from Fig. 2 (see also Table III): (e) ; (f) ; (g) ; (h) ; (i) ; (j) . For reactions (e), (f), (i), and (j) only the results for the two end states and the half way mixed state are given similar to the reactions (a)–(d) in Fig. 2. For reactions (g) and (h) a more closely spaced set of coupling parameter values has been used, namely, , for reaction (h) an additional state at . The gap energies are strictly monotonously decreasing in [see plot of charging curves for reactions (g) and (h) in Fig. 5]. Energies at the end of the accumulative average plots are the final values used in the calculation of the free energies.

Accumulative averages of vertical energy gaps as a function of MD time for the deprotonation and dehydrogenation reactions of Fig. 1. Labeling is continued from Fig. 2 (see also Table III): (e) ; (f) ; (g) ; (h) ; (i) ; (j) . For reactions (e), (f), (i), and (j) only the results for the two end states and the half way mixed state are given similar to the reactions (a)–(d) in Fig. 2. For reactions (g) and (h) a more closely spaced set of coupling parameter values has been used, namely, , for reaction (h) an additional state at . The gap energies are strictly monotonously decreasing in [see plot of charging curves for reactions (g) and (h) in Fig. 5]. Energies at the end of the accumulative average plots are the final values used in the calculation of the free energies.

Free energy diagram for quinone reduction. The numbers alongside the arrows give the calculated reaction free energies (units in eV) in the reverse direction, i.e., hydroquinone oxidation to quinone. The free energy changes from experimental (Ref. 10) NHE reduction potentials and are given in parentheses. The numbers in bold are the free energy changes of four full reactions. The three red (blue) arrows indicate three alternative pathways for the , respectively reactions. The energies marking the arrows are the corresponding DFTMD reaction free energies which according to Hess’s law should be the same for each reaction.

Free energy diagram for quinone reduction. The numbers alongside the arrows give the calculated reaction free energies (units in eV) in the reverse direction, i.e., hydroquinone oxidation to quinone. The free energy changes from experimental (Ref. 10) NHE reduction potentials and are given in parentheses. The numbers in bold are the free energy changes of four full reactions. The three red (blue) arrows indicate three alternative pathways for the , respectively reactions. The energies marking the arrows are the corresponding DFTMD reaction free energies which according to Hess’s law should be the same for each reaction.

Correlation between vertical energy gap and solvation structure for half reactions [panels (a) and (d)], [(b) and (e)], and [(c) and (f)]. Black squares in (a)–(c) give the vertical energy gaps as determined by averaging over MD trajectories at the indicated values of the coupling parameter . Red circles are the corresponding coordination numbers obtained by integrating over the first solvation shell of the rdf between quinone O atoms and solvent H atoms (denoted by ). The rdfs at the end points (black solid curves), (blue dashed-dotted curves), and one intermediate point (red dashed curves) are given in panels (d)–(f) on the right. The blue dashed line in (a) is the linear fit of the vertical energy gap as a function of coupling parameter with the regression coefficient .

Correlation between vertical energy gap and solvation structure for half reactions [panels (a) and (d)], [(b) and (e)], and [(c) and (f)]. Black squares in (a)–(c) give the vertical energy gaps as determined by averaging over MD trajectories at the indicated values of the coupling parameter . Red circles are the corresponding coordination numbers obtained by integrating over the first solvation shell of the rdf between quinone O atoms and solvent H atoms (denoted by ). The rdfs at the end points (black solid curves), (blue dashed-dotted curves), and one intermediate point (red dashed curves) are given in panels (d)–(f) on the right. The blue dashed line in (a) is the linear fit of the vertical energy gap as a function of coupling parameter with the regression coefficient .

## Tables

Force constants and structural parameters for the restraining potentials defined in Eq. (19). , , and are the force constants for bonds, bends, and torsions, respectively, in atomic unit. , , and are the corresponding equilibrium values in atomic unit for bond lengths and radians for angles and dihedrals. For weak acids . For strong acids is applied to the OH bond only (for both O–H bonds are restrained).

Force constants and structural parameters for the restraining potentials defined in Eq. (19). , , and are the force constants for bonds, bends, and torsions, respectively, in atomic unit. , , and are the corresponding equilibrium values in atomic unit for bond lengths and radians for angles and dihedrals. For weak acids . For strong acids is applied to the OH bond only (for both O–H bonds are restrained).

Vibrational frequencies , vibrational temperatures , moments of inertia , and rotational temperatures used in the correction for restraining potentials and zero point motion (see Sec. II E). The molecule is the hydronium ion with one proton replaced by a dummy of the same mass but without charge. Similarly is a quinone with the acid proton replaced by a dummy proton. The three vibrational frequencies given for are the normal modes of the dummy atom with all of the quinone anion atoms fixed.

Vibrational frequencies , vibrational temperatures , moments of inertia , and rotational temperatures used in the correction for restraining potentials and zero point motion (see Sec. II E). The molecule is the hydronium ion with one proton replaced by a dummy of the same mass but without charge. Similarly is a quinone with the acid proton replaced by a dummy proton. The three vibrational frequencies given for are the normal modes of the dummy atom with all of the quinone anion atoms fixed.

DFTMD free energy changes , variances of energy gaps , and reorganization energies for the half reactions of Fig. 1. Subscripts LR and TI indicate the numerical method used to calculate free energy changes. LR stands for linear response (end point) approximation [Eq. (13)] and TI for thermodynamic integration using additional intermediate coupling parameter states (the number of integration points can be deduced from Figs. 2 and 3). Subscripts 0 and 1 denote the value of the coupling parameter corresponding to begin and end states [see Eq. (9)]. The reorganization energy has been calculated according to Eq. (15) allowing for nonlinear solvent response. All simulations have been performed in a cubic periodic cell box with edge except the results of redox reactions in parentheses, which were calculated in a smaller box . Energies are in eV and variances are in .

DFTMD free energy changes , variances of energy gaps , and reorganization energies for the half reactions of Fig. 1. Subscripts LR and TI indicate the numerical method used to calculate free energy changes. LR stands for linear response (end point) approximation [Eq. (13)] and TI for thermodynamic integration using additional intermediate coupling parameter states (the number of integration points can be deduced from Figs. 2 and 3). Subscripts 0 and 1 denote the value of the coupling parameter corresponding to begin and end states [see Eq. (9)]. The reorganization energy has been calculated according to Eq. (15) allowing for nonlinear solvent response. All simulations have been performed in a cubic periodic cell box with edge except the results of redox reactions in parentheses, which were calculated in a smaller box . Energies are in eV and variances are in .

Full reactions for the determination of the oxidation free energy versus NHE [(a)-(d)], Brønsted acidity [(e)–(h)], and dehydrogenation free energies [(i) and (j)] of the half reactions of Table III. The data marked repeat the free energies of the column of Table III. are the experimental free energies as collected in Ref. 10 (same data are given again in Fig. 4). is the free energy cost of deprotonating a hydronium ion estimated by aligning simulation and experimental free energies. The average obtained in this way from the acid dissociation reactions [(e)–(h)] have been used to compute the oxidation free energies and acidities denoted by (see Sec. III C for more details). Figures in parentheses include the corrections for the restraining potentials and zero point motion (see Sec. II E). All energies are in eV. The mean deviation from experiment of the restraint corrected oxidation free energies is −0.52 eV.

Full reactions for the determination of the oxidation free energy versus NHE [(a)-(d)], Brønsted acidity [(e)–(h)], and dehydrogenation free energies [(i) and (j)] of the half reactions of Table III. The data marked repeat the free energies of the column of Table III. are the experimental free energies as collected in Ref. 10 (same data are given again in Fig. 4). is the free energy cost of deprotonating a hydronium ion estimated by aligning simulation and experimental free energies. The average obtained in this way from the acid dissociation reactions [(e)–(h)] have been used to compute the oxidation free energies and acidities denoted by (see Sec. III C for more details). Figures in parentheses include the corrections for the restraining potentials and zero point motion (see Sec. II E). All energies are in eV. The mean deviation from experiment of the restraint corrected oxidation free energies is −0.52 eV.

Quinone energetics in vacuum as computed in this work using BLYP compared with the results taken from the literature. The DFT functional used in Refs. 43, 48, and 49 is B3LYP and B3PW91 in Ref. 97. The unit of energy is eV.

Quinone energetics in vacuum as computed in this work using BLYP compared with the results taken from the literature. The DFT functional used in Refs. 43, 48, and 49 is B3LYP and B3PW91 in Ref. 97. The unit of energy is eV.

Article metrics loading...

Full text loading...

Commenting has been disabled for this content