Redox Potentials from Ab Initio Molecular Dynamics and Explicit

Subscriber access provided by Columbia University Libraries

Article

Redox potentials from ab initio molecular dynamics and explicit entropy calculations: application to transition metals in aqueous solution Miguel A. Caro, Olga Lopez-Acevedo, and Tomi Laurila J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b00314 • Publication Date (Web): 17 Jul 2017 Downloaded from http://pubs.acs.org on July 24, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 13

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Redox potentials from ab initio molecular dynamics and explicit entropy calculations: application to transition metals in aqueous solution Miguel A. Caro,1, 2, ∗ Olga Lopez-Acevedo,2, 3 and Tomi Laurila1 1

Department of Electrical Engineering and Automation, Aalto University, Espoo, Finland 2 COMP Centre of Excellence in Computational Nanoscience, Department of Applied Physics, Aalto University, Espoo, Finland 3 Departamento de Ciencias B´ asicas, Universidad de Medell´ın, Colombia (Dated: June 14, 2017)

We present a complete methodology to consistently estimate redox potentials strictly from first principles, without any experimental input. The methodology is based on (i) ab initio molecular dynamics (MD) simulations, (ii) all-atom explicit solvation, (iii) the two-phase thermodynamic (2PT) model and (iv) the use of electrostatic potentials as references for the absolute electrochemical scale. We apply the approach presented to compute reduction potentials of the following redox +2/+3 couples: Cr+2/+3 , V+2/+3 , Ru(NH3 )6 , Sn+2/+4 , Cu+1/+2 , FcMeOH0/+1 and Fe+2/+3 (in aque0/+1 ous solution) and Fc (in acetonitrile). We argue that fully quantum-mechanical simulations are required to correctly model the intricate dynamical effects of the charged complexes on the surrounding solvent molecules within the solvation shell. Using the proposed methodology allows for a computationally efficient and statistically stable approach to compute free energy differences, yielding excellent agreement between our computed redox potentials and the experimental references. The root-mean-square deviation with respect to experiment for the aqueous test set and the two exchange-correlation density functionals used, PBE and PBE with van der Waals corrections, are 0.659 V and 0.457 V, respectively. At this level of theory, the functional employed and its ability to correctly describe each particular molecular complex seems to be the factor limiting the accuracy of the calculations. I.

INTRODUCTION

Free-energy differences drive chemical reactions and therefore determine the characteristics of a vast array of physical and chemical processes. From electrochemistry to catalysis, through phase stability and equilibrium properties, the change in a system’s free energy is the key quantity to be evaluated. In electrochemistry, the position of redox potentials and reaction rates involved in the electrochemical detection of important biomolecules directly depend on the free energy difference between the oxidized and reduced forms of the biomolecule, both of which could be either adsorbed on the electrode’s surface or in solution, as a function of electrode potential. Unfortunately, atomistic simulations, either based on quantum mechanical approaches or classical force fields, can only easily access potential (or internal) energy values, whereas free-energy calculations are computationally significantly more demanding. In addition, simulating the effect of tuning the electrode’s potential remains a challenging issue, in particular within the context of densityfunctional theory (DFT). Computational electrochemistry is a relatively underdeveloped (and fairly recent) branch of computational physical chemistry.1–7 This is mostly because it combines two classical challenges of general atomistic simulations. On the one hand, computing free energies and free energy differences relies on long molecular dynamics (MD) simulations that are required to sample a statisticallysignificant portion of the momentum-position parameter space.8 On the other hand, explicit interaction between electrolyte and electrode is often limited to charge-

neutral systems where only one value of all the possible electrode potentials is surveyed.1–3 Other possibilities to model the effect of varying electrochemical potentials involve removal of electrons7 and indirect references, such as a “computational hydrogen electrode”6 and electrostatic potential differences.5 Directly tuning electrode potentials while retaining the full atomistic detail remains both challenging and computationally expensive.9 The possibility to apply methodologies based on nonequilibrium Green’s functions to treat electrode/solution interfaces in a similar way to semiconductor/metal interfaces10 also emerges as an attractive prospect. However, the community needs to work around the difficulties involved in affordable computation of free energies within simulation frameworks which are intrinsically computationally demanding. This is an effort worth pursuing, since being able to accurately simulate electrochemical processes from first principles would open the door to new applications in fields of great technological importance. Computational electrochemistry with predictive power will enable computer-aided design for the next generation of materials with targeted electrochemical and electrocatalytical applications: enhanced electrochemical detection, more efficient fuel cells, chemical energy conversion, photoelectrochemical solar cells, and so on. In this paper, to accurately compute redox potentials strictly from first principles only, we tackle the two aforementioned challenges. In a first step, we focus on addressing directly the issue of computationally-affordable free energy calculations. We apply a modified two-phase thermodynamic (2PT) approach11 to explicitly calculate entropies (and thus free energies) from ab initio MD simula-

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Page 2 of 13

Page 3 of 13

Journal of Chemical Theory and Computation 3

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

tablished, the connection between the absolute potential Vabs and the equivalent potential under PBC is Veq (Vabs ) = Vabs + Vref + ∆Vref/vac .

(3)

Within the simulation, the energy of an electron in an energy state corresponding to electrode potential Vabs is −eVeq (Vabs ) rather than the usual −eVabs . Only when Vvac = 0 in the simulation do both values coincide. The average electrostatic (Hartree) potential of the solvent becomes an obvious candidate to serve as reference, as schematically shown in Fig. 1. Then, the correction term from Eq. (1) becomes: 1 sol sol/vac δV = − (EHartree + ∆EHartree ). e sol/vac

(4)

∆EHartree is a solvent-dependent quantity, and in general will also show a small dependence on the DFT functional used. Lucking et al. proposed a strategy to estimate this value from MD simulations.5 Their PBE value for water/vac water is ∆EHartree = 2.97 eV. We have estimated the value for acetonitrile using the same approach, yielding acetonitrile/vac ∆EHartree = 3.15 eV. Having established the corrections that need to be applied to the computed energies and potentials, the standard absolute potential for the reaction can then be calculated as V0 such that GR −GO (Vabs = V0 ) = 0. In other words, the standard potential is the potential value for which O and R lie at the same energy, i.e., there is an equal probability of finding the molecule in either oxidized or reduced form. In this paper we work under the assumption that PR VR ≈ PO VO (c.f. Eq. (1)), which means that we can neglect pressure-volume effects and approximate Gibbs free energy differences as Helmholtz free energy differences, ∆G ≈ ∆A. This brings us to the pressing matter of how to estimate free energies, or rather, the entropic term in Eq. (1). Several strategies exist to tackle the issue of free energy calculations from MD simulations.8 Because of computational costs and difficulties in the implementation, many of these strategies have not been consistently applied to ab initio molecular dynamics (AIMD) and are usually circumscribed to the realm of classical MD. A popular strategy which has been used in combination with AIMD is thermodynamic integration (TI). Within the framework of TI, the initial (i) and final (f) states of the system (e.g., two different oxidation states of an analyte molecule) are coupled via the parameterization of a combined Hamiltonian, Hcomb = λHf + (1 − λ)Hi . The change in free energy can then be computed by evaluating the configurational average h∂Hcomb /∂λiλ via molecular dynamics, which means that the method can be more efficient than other approaches based on stochastic techniques, such as Monte Carlo algorithms. For ab initio molecular dynamics, constructing a combined Hamiltonian is far from trivial and computationally-demanding simulations limit the amount of configurational sampling that can be performed. Wang and Van Voorhis recently

performed multiscale quantum-mechanical/molecularmechanical (QM/MM) simulations to estimate standard potentials of popular redox couples using TI.4 Because of the use of DFT, they could only simulate the initial and final states (λ = 0, 1) and had to approximate the configurational average integral in λ by a finite difference (this is known as the “linear approximation”). The main contribution from Ref. 4 is the realization that, to obtain accurate redox potentials from computational simulations, it is vital to include the explicit interaction of the solvent molecules with the solute. We will confirm that result in the present study. Work on this topic using thermodynamic integration has also been done by the Sprik group, who rely on a “computational hydrogen electrode”, which allows direct comparison of redox potentials with experiment.6 However, the problem remains regarding the computational cost of sampling configuration space within thermodynamic integration and its applicability to quantum-mechanical systems where the coupling scheme cannot be properly applied. A recent approach particularly geared towards cheap calculation of thermodynamic properties from MD is the so-called two-phase thermodynamic (2PT) model.11 The main idea behind 2PT is that the 3N degrees of freedom of a liquid system can be expressed as a combination of 3f N “gas-like” degrees of freedom and 3(1 − f )N “solidlike” degrees of freedom. 2PT relies on a vibrationalmode analysis based on the “density of states” (DoS),19 calculated as the velocity spectrum of the atoms (the Fourier transform of the velocity auto-correlation function). This analysis allows to compute the contributions to thermodynamic properties arising from solid- and gaslike degrees of freedom separately, as a functional of the DoS. The model was extended from monoatomic liquids in 200311 to polyatomic liquids in 201020 to liquid mixtures in 2012.21,22 We recently assessed the performance of the method with a variety of liquid mixtures and force fields.23 The great advantage of 2PT is to obtain converged thermodynamic properties from short MD simulations (a rule of thumb being 20 ps of dynamics), which immediately points at the utilization of 2PT to compute free energies in computationally-demanding situations, including very large molecular systems, high sampling (for statistical purposes) and ab initio MD. Surprisingly though, except for some exceptions,24 most of the applications of 2PT to date have been limited to classical MD, and the combination of 2PT and AIMD has not been widely employed. We have very recently created our own implementation of 2PT, called DoSPT,29 where several improvements have been made, including enhanced convergence of results23 and improved suitability for AIMD. A crucial difference between AIMD and classical MD is the possibility of bond breaking and bond formation occurring during the course of the dynamics. This can become an issue for instance when a charged redox couple induces dissociation of water molecules into hydroxyl and hydronium groups, which we observe in several of our simulations.




Page 4 of 13

Page 5 of 13




Page 6 of 13 6

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

TABLE I. Solvation shells for different redox complexes and probability to find a particular configuration during the MD. Only probabilities larger than 1% are shown (values are rounded up to nearest 1%). These numbers are for the 20 ps of MD used to compute the redox potentials: the solvation shells change dynamically and the solvation structure at equilibrium, especially for Sn+2 , may differ. See text. PBE

PBE+TS

Most likely

2nd most likely

3rd most likely

Most likely

2nd most likely

3rd most likely

Cr+2

(H2 O)5 [100%]

n/a

n/a

(H2 O)4 [51%]

(H2 O)5 [49%]

n/a

Cr+3

(H2 O)5 OH−1 [83%]

(H2 O)6 [17%]

n/a

(H2 O)6 [61%]

(H2 O)5 OH−1 [39%]

n/a

V+2

(H2 O)5 [100%]

n/a

n/a

(H2 O)6 [100%]

n/a

n/a

V+3

(H2 O)3 (OH)−2 [93%] 2

(H2 O)4 OH−1 [7%]

n/a

(H2 O)3 (OH)−2 [97%] 2

(H2 O)2 (OH)−3 [2%] 3

(H2 O)4 OH−1 [1%]

Sn

+2

(H2 O)3 OH

Sn

+4

(H2 O)2 (OH)−3 3

−1

[48%]

(H2 O)2 OH

[81%]

(H2 O)3 (OH)−2 2

−1

[37%] [11%]

(H2 O)3 [12%] H2 O(OH)−4 4

[8%]

(H2 O)2 OH

−1

(H2 O)2 (OH)−3 3

[84%] [61%]

−1

(H2 O)3 OH

[11%]

H2 O(OH)−4 4

H2 O(OH)−2 [5%] 2

[39%]

n/a

Cu+1

(H2 O)3 [79%]

(H2 O)2 [21%]

n/a

(H2 O)2 [58%]

(H2 O)3 [38%]

H2 OOH−1 [5%]

Cu+2

(H2 O)4 [81%]

(H2 O)3 OH−1 [19%]

n/a

(H2 O)4 [78%]

(H2 O)3 OH−1 [21%]

n/a

+2

(H2 O)5 [100%]

n/a

n/a

(H2 O)6 [87%]

(H2 O)5 [12%]

n/a

Fe+3

(H2 O)6 [93%]

(H2 O)5 OH−1 [7%]

n/a

(H2 O)6 [87%]

(H2 O)5 OH−1 [13%]

n/a

Fe

Note that experimental conditions differ somewhat from those imposed in our simulation. Experimentally, standard redox potentials are measured at 298.15 K and 1 mol/L concentrations. We carry out our simulations at 300 K rather than 298.15 K. This small difference has negligible impact on the results. Regarding concentration, from our 12 ˚ A and 14 ˚ A boxes we have that 1 solute molecule corresponds to approximate concentrations of 0.96 and 0.61 mol/L, respectively. At these levels, given the outer-sphere nature of the probed couples, the experimental potentials should be comparable to the results of our simulations with no concentration correction. To gain further insight into the factors that impact the simulation results, we have conducted a detailed analysis of the solvation shell structure for all the systems studied. The most likely configurations are summarized in Table I. Radial distribution functions and more detailed information on solvation shell structures are given in Ref. 43. These results indicate a sizable functional dependence on the predicted solvation shell, with either 5or 6-fold coordination for the +2/+3 metal cation couples. This includes partial hydroxylation of the water molecules surrounding some of the +3 cations, particularly V+3 . Cu+1/+2 shows lower coordination, between 2- and 4-fold, which increases with oxidation number, and partial hydroxylation in all cases. Sn+2/+4 shows the most complex shell structure, with 3- to 5-fold coordination and partial hydroxylation, especially in the +4 state, where up to 4 out of 5 surrounding oxygens are part of OH groups rather than H2 O groups. Solvation shell structures across functionals seem to make a notable difference only for V+2 , where PBE+TS predicts a (H2 O)6 shell (rather than (H2 O)5 with PBE) and further stabilization of the complex; this is likely the reason for the improved performance of PBE+TS for redox potential prediction of the V+2/+3 couple (Fig. 2). It is crucial to note that the solvation shells change dynamically, especially for Sn. The numbers in Table I are those calculated for the 20 ps of MD which were used to compute the redox potentials and may not fully rep-

resent the most stable solvation shell configuration. For instance, the Sn+2 shell, for which we have carried out extended simulations (see the discussion on statistical convergence at the end of this section) settles eventually after 60 ps of MD at a predominantly (H2 O)4 shell configuration, with some probability for (H2 O)3 shell and (H2 O)3 OH− shell. The latter solvation shells lie of the order of 0.5 eV lower than those shown in Table I, in terms of potential energy; in terms of free energy, they only lie 0.15 eV below. This means that the solvation shell dynamics can be entropy driven, although it is possible that the free energy barrier becomes too large to overcome once the Sn+2 complex has settled in the lower potential energy configuration. On the other hand, the Sn+4 complex, for which we also have extended MD, remains quite similar to the configuration initially observed, with the (H2 O)2 (OH)−3 and H2 O(OH)−4 shells 3 4 dominating (the relative probabilities fluctuate throughout the MD). These results further highlight the importance of including explicit interaction between molecule and solvent, to a level of detail which can only be achieved from either fully quantum-mechanical simulations or perhaps QM/MM schemes with a sizable number of solvent molecules included in the QM part. Another question that merits discussion in the context of solvation shell structure is the role of pH in the simulations, in particular given supercell sizes where one single proton already corresponds to pH = 0. In fact, most of the listed redox couples are not stable at neutral pH, as can be seen from their experimental Pourbaix diagrams where, e.g., Sn+2/+4 or V+2/+3 are only stable at very low pHs.44,45 This means that any experimental results utilizing these couples are also from this low pH region. If one looks in more detail for example at the Pourbaix diagram of Sn, the equilibrium between Sn+4 and Sn(OH)4 occurs at very low pH (approximately 0): below pH = 0 Sn+4 is experimentally more stable than Sn(OH)4 . Experimentally, the protons in solution will compete with the Sn+4 cation to bind the OH− groups, and the removal of hydroxyls takes place gradually, i.e.,


Page 7 of 13




Page 8 of 13 8

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

compare different scales.48 Since all the values we compute within our methodological framework are naturally given in the absolute scale, we propose that the present approach can provide a definitive link between external and internal references via placing the thermodynamic redox potential of the Fc0/+1 couple in acetonitrile on the absolute scale. The values we obtain are 5.63 V and 6.24 V with PBE and PBE+TS, respectively. To obtain a highly accurate reference in this fashion, the calculations could be repeated with a larger simulation box and/or a hybrid functional (see related discussion in Sec. V). To assess the statistical stability of our approach, we carried out extended simulations for FcMeOH, Sn+2 and Sn+4 . We chose these systems because they were simultaneously the least computationally expensive to simulate and also representative examples. A total of 180 ps of dynamics was sliced into nine 20 ps pieces, that were subsequently analyzed in the same fashion as the previous results. The statistical analysis of the results for potential energy, entropy (separately for solvent and molecule) and free energy is shown in Fig. 5. Our statistical test assumes that the simulations are fully converged after 180 ps, and computes the statistical error incurred by using only 20 ps of dynamics. For Sn+2 the maximum errors incurred are 327 meV, 319 meV and 138 meV for internal energy, solvent entropy and total free energy, respectively. The corresponding typical errors by using only 20 ps of MD (i.e., the standard deviation of the results) are 209 meV, 181 meV and 73 meV, as indicated in Fig. 5. We attribute the larger statistical error shown in Fig. 5 for FcMeOH compared to Sn+2 to the relative size of the molecule and the solvent box that contains it. The results for Sn+4 are similar to Sn+2 , with typical deviations of 156 meV, 125 meV and 135 meV for potential energy, solvent entropy and free energy, respectively. The maximum deviations are 272 meV, 220 meV and 287 meV, respectively. The error in free energy is dominated by the last 20 ps slice from t = 160 ps to t = 180 ps. This means that in the worst case scenario, by using 20 ps MD slices, the maximum statistical error associated to the Sn+2 /Sn+4 redox potential calculation would be 155 mV + 287 mV√= 425 mV, whereas the typical statistical error would be 732 + 1352 mV = 153 mV. Therefore, within this treatment, we can claim that the typical error in the calculated free energies is about 100 meV to 150 meV for solvated cations and 200 meV for solvated molecules, given the supercell size and time sampling that we are employing. Statistically, these uncertainties propagate to the calculated redox potentials with a multiplicative √ factor of 2. These statistical uncertainties can be reduced by increasing either sampling, system size, or both. Therefore, although our 20 ps results are not fully converged, we can provide reliable estimates of what kind of errors should be expected by employing the methodology presented. In addition, at this level of theory (PBE or PBE+TS) the error in the energies arising from the choice of functional will be, at the very least, of the same order of magnitude as the statistical error and, possibly,

even larger. Given the high computational cost of these simulations, one could argue that converging the free energies well below the accuracy of the functional would be a waste of computational resources.

In all three cases studied in Fig. 5, the calculated free energy values are significantly more stable than the potential energy values, which highlights an opposite correlation between entropic deviations and potential energy deviations from the configurational average (i.e., when the entropic term is overestimated the potential energy term is underestimated, and vice versa). For instance, one can observe an apparent drift in potential energy for Sn+2 , where the difference between the potential energy between the first and last 20 ps slices is approximately 550 meV. This drift could be due to either insufficient equilibration or a slow potential energy oscillation. Insufficient equilibration is always bound to be a possible issue for AIMD. However, the drift in free energy, which is our main quantity of interest, is only 150 meV. Since the evolution of the system is driven by finding the minimal free energy configuration (rather than potential energy), this drift is not significant. This analysis of free energy statistical stability offers strong support for the suitability of the present methodology for free energy calculations from AIMD simulations.

To finalize our discussion, it is worth adding a note on computational performance of the method. The 2PT computations are post-processing only, and of negligible cost compared to the MD simulations. The CPU cost associated with the MD is typically about 17 000 CPU hours per calculated redox potential (usually one spin-polarized and one spin-paired simulation per couple). This performance was achieved on CSC’s Taito supercluster with Intel Sandy Bridge 2.6 GHz processors (Intel E5-2670, 64bits). Higher quality simulations at the same level of theory (in terms of removing statistical uncertainty) can be carried out by increasing the sampling time, the system size, or a combination of both, at the associated increase in CPU time. The cost of the simulation grows linearly with time sampling for fixed size and (generally for DFT, although it depends on the code) as the cube of the system size for fixed duration of the dynamics. The use of a different DFT functional within the GGA family should not significantly impact the CPU time requirements. Going up one level of theory, for instance running MD with a hybrid functional, is still quite computationally expensive, albeit doable. For instance, based on testing we carried out, we estimate to be able to do HSE0649 AIMD with settings similar to those presented here at a cost of about 250 000 CPU hours per redox couple. This is certainly an expendable amount of CPU resources by today’s standards for any particular system whose technological or scientific importance warrants the investment on high quality simulations.


Page 9 of 13


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

V.

ACKNOWLEDGMENTS

CONCLUSIONS

From the results presented here we extract the following conclusions. a) Including explicit solvation effects is crucial in order to accurately estimate free energy differences between different oxidation states of a redox couple. The complexity of solvation shells that we have observed in some of these simulations, in particular for the aqua ions, point to the conclusion that full quantummechanical simulations are necessary, at least to model the first few solvation layers, because some of the metal cations will induce O-H bonds to break during the dynamics (especially Sn and V, and to a lesser extent Cu and Cr). Bond breaking within the solvation shell cannot be accurately (or at all) captured by modeling with classical force fields. b) Inclusion of van der Waals corrections can significantly improve the quality of the predicted energy differences in some cases, although not systematically since the intrinsic inability of the base functional used to correctly describe strongly correlated electrons and electronic self-interaction is also present in the corresponding van der Waals-corrected functional. Note that at the small level of disagreement observed, the intrinsic error in the energies predicted by the functional will dominate the overall error. c) The statistical error associated to employing the modified 2PT formalism is quite small: 190 meV, 73 meV and 135 meV for the calculated free energies of the tested FcMeOH, Sn+2 and Sn+4 , respectively. This makes the methodology presented here highly attractive to compute free energies and free-energy differences in ab initio MD simulations, not only in the realm of computational electrochemistry but also in general. With the method presented, excellent statistical stability can be achieved for free energy calculations, with uncertainties below 200 meV. This uncertainty can be further reduced by a combination of increased system size and extended MD. The accumulated error arising from the different contributions: entropy, internal energy, finite-size corrections and electrostatic potential alignment, amount to RMSDs with respect to experiment of 0.659 V and 0.457 V for PBE and PBE+TS, respectively, for the test set studied. The main source of error at this point is the DFT functional used, and its ability to correctly describe challenging systems, for instance Cr and V in this paper. This is evidenced by the fact that previous QM/MM simulations with a hybrid functional4 achieved a better description of Cr and V aqua ions despite only including one solvation layer within the QM region. The use of more sophisticated DFT exchangecorrelation functionals could lead to improved overall accuracy of the approach. All in all, we conclude that high quality estimates of unknown redox potentials can be conducted using the methodology presented in this paper at affordable CPU cost. This opens the door to first-principles electrochemistry with predictive power, with innumerable applications in industry and research.

Computational resources for this project were provided by CSC–IT Center for Science through the Taito supercluster and by Aalto University’s Science-IT Project through the Triton cluster. Funding from the Academy of Finland through grants no. 284621, 285015 and 285526 is gratefully acknowledged. Appendix: Details of the 2PT calculation

The two-phase thermodynamics (2PT) model of Lin and collaborators11 builds up on the idea introduced by Berens et al.19 of computing thermodynamic properties of a molecular system from its density of states. For more details, the reader is referred to previous literature on this topic;11,19,20,22,23 here we simply outline the procedure to complement the discussion of our manuscript and to make the explanation of the theoretical framework more self contained. The expression of the partition function of a (quantum) harmonic oscillator is known analytically: q(ν) =

exp(−βhν/2) , 1 − exp(−βhν)

(A.1)

where β = 1/(kB T ) is the inverse of the product of Boltzmann’s constant times the temperature, and h is Planck’s constant. For a system made up of an ensemble of 3N such harmonic oscillators, each at frequency νi , the total partition function Q and its logarithm ln Q can be expressed as Q=

3N Y

q(νi ),

ln Q =

i

3N X

ln[q(νi )].

(A.2)

i

In the limit where the frequencies are continuously distributed, one can extend the summation to an integral, provided that an appropriate weighting function, giving the distribution of “oscillator density” as a frequencydependent quantity, can be found. This density is the density of states S(ν), which allows to express Eq. (A.2) as Z ∞ dν S(ν) ln [q(ν)] , (A.3) ln Q = 0

and must of course integrate to the total number of degrees of freedom: Z ∞ dν S(ν) = 3N. (A.4) 0

As shown by Berens et al.,19 this density of states can be obtained as a sum over the individual density of states of each atom (or particle) in the system: 2 X S(ν) = mj skj (ν), (A.5) kB T


j,k


Page 10 of 13 10

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

where mj are the atomic masses and skj (ν) is given by the squared modulus of the Fourier transform of the atomic velocities vjk : skj (ν)

2 Z 1 τ k −2πiνt = lim dt vj (t) e . τ →∞ τ 0

(A.6)

From the thermodynamic relation between entropy and partition function S = kB T

∂ ln Q + kB ln Q, ∂T

(A.7)

and the fact the we know the analytical form of q(ν), we can differentiate under the integral sign and obtain a weighting function WS (ν) that allows to compute entropy as a functional of the DoS: Z ∞ dν S(ν) WS (ν), (A.8) S= 0

where WS (ν) = kB

βhν −βhν − ln (1 − e ) . eβhν − 1

(A.9)

Because we are constrained here to the harmonic approximation, this theoretical framework works remarkably well for solids but fails for liquids. For liquids, lowfrequency components of the DoS, related to diffusion, are important and contribute heavily to the entropy. Lin et al.11 realized that one could conceptually partition the DoS into a solid-like part S s (ν) and a gas-like part S g (ν). That is, one ends up with a separate description of solidlike and gas-like degrees of freedom, and the entropy is thus calculated as the sum of two integrals: Z ∞ Z ∞ dν S g (ν) WSg (ν), dν S s (ν) WSs (ν) + S= 0

∗ 1

2

3

4

0

(A.10)

[email protected] Rossmeisl, J.; Nørskov, J. K.; Taylor, C. D.; Janik, M. J.; Neurock, M. Calculated phase diagrams for the electrochemical oxidation and reduction of water over Pt (111). J. Phys. Chem. B 2006, 110, 21833. Sk´ ulason, E.; Karlberg, G. S.; Rossmeisl, J.; Bligaard, T.; Greeley, J.; J´ onsson, H.; Nørskov, J. K. Density functional theory calculations for the hydrogen evolution reaction in an electrochemical double layer on the Pt(111) electrode. Phys. Chem. Chem. Phys. 2007, 9, 3241. Rossmeisl, J.; Sk´ ulason, E.; Bj¨ orketun, M. E.; Tripkovic, V.; Nørskov, J. K. Modeling the electrified solid– liquid interface. Chem. Phys. Lett. 2008, 466, 68. Wang, L.-P.; Van Voorhis, T. A Polarizable QM/MM Explicit Solvent Model for Computational Electrochemistry in Water. J. Chem. Theory Comput. 2012, 8, 610.

where the gas DoS is calculated from the DoS of a hardsphere gas with the same diffusive properties as the liquid being described, via the “fluidicity” parameter f : that is, S g (ν) depends parametrically on f . This parameter determines how many degrees of freedom are treated as solid and how many as liquid, with the obvious requirement that the total number of degrees of freedom must remain unchanged: 3N s = 3(1 − f )N , 3N g = 3f N . In practice, f is computed from the ratio of real diffusion constant (calculated from the MD) and ideal diffusion constant (obtained from the hard-sphere formalism), f = D/DHS . The solid weighting function in Eq. (A.10) is taken from the harmonic treatment outlined earlier. The gas weighting function and the form of the gas DoS are retrieved from the hard-sphere theory of gases. All the details are given in Ref. 11. In addition, the descriptions of translational, rotational and vibrational degrees of freedom in molecular systems are also done separately: translational degrees of freedom contribute much more significantly to the entropy and treating all the types of degrees of freedom on the same footing would result in a gross oversimplification of the problem. The different weighting functions, including the type of degree of freedom, have been obtained by Lin et al.20 Finally, one needs to incorporate the entropy of mixing contributions that arise from having more than one substance in the system. The analysis of entropy of mixing in the context of 2PT has been done by Lai et al.,21 Pascal and Goddard22 and, in somewhat more detail, in our recent paper.23 However, because in this manuscript we are in the dilute limit (only one solute molecule) and only interested in entropy differences, mixing entropy contributions do not significantly affect our results. The 2PT method is implemented in the DoSPT code, together with several improvements, such as accelerated convergence23 and the topology reconstruction feature already discussed in this manuscript, Eq. (5), which was required in the present case because of the dynamical effects of protonation/deprotonation induced by our charged metal complexes.

5

6

7

8

9

Lucking, M.; Sun, Y.-Y.; West, D.; Zhang, S. Absolute redox potential of liquid water: a first-principles theory. Chem. Sci. 2014, 5, 1216. 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. Holmberg, N.; Laasonen, K. Ab Initio Electrochemistry: Exploring the Hydrogen Evolution Reaction on Carbon Nanotubes. J. Phys. Chem. C 2015, 119, 16166. Christ, C. D.; Mark, A. E.; van Gunsteren, W. F. Basic ingredients of free energy calculations: a review. J. Comput. Chem. 2010, 31, 1569. Aarva, A.; Laurila, T.; Caro, M. A. Doping as a means to probe the potential dependence of dopamine asdorption on carbon-based surfaces: A first-principles study. J. Chem.


Page 11 of 13


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

Phys. 2017, accepted for publication. Stradi, D.; Martinez, U.; Blom, A.; Brandbyge, M.; Stokbro, K. General atomistic approach for modeling metalsemiconductor interfaces using density functional theory and nonequilibrium Green’s function. Phys. Rev. B 2016, 93, 155302. Lin, S.-T.; Blanco, M.; Goddard III, W. A. The two-phase model for calculating thermodynamic properties of liquids from molecular dynamics: Validation for the phase diagram of Lennard-Jones fluids. J. Chem. Phys. 2003, 119, 11792. Bard, A. J.; Faulkner, L. R. Electrochemical methods: Fundamentals and applications, 2nd ed.; Wiley, New York, 2001. Makov, G.; Payne, M. C. Periodic boundary conditions in ab initio calculations. Phys. Rev. B 1995, 51, 4014. Rusnak, A. J.; Pinnick, E. R.; Calderon, C. E.; Wang, F. Static dielectric constants and molecular dipole distributions of liquid water and ice-Ih investigated by the PAWPBE exchange-correlation functional. J. Chem. Phys. 2012, 137, 034510. Caro, M. A.; M¨ aa ¨tt¨ a, J.; Lopez-Acevedo, O.; Laurila, T. Energy band alignment and electronic states of amorphous carbon surfaces in vacuo and in aqueous environment. J. Appl. Phys. 2015, 117, 034502. Wu, Y.; Chan, M. K. Y.; Ceder, G. Prediction of semiconductor band edge positions in aqueous environments from first principles. Phys. Rev. B 2011, 83, 235301. Van de Walle, C. G.; Martin, R. M. Theoretical calculations of heterojunction discontinuities in the Si/Ge system. Phys. Rev. B 1986, 34, 5621. Zhang, S. B.; Tom´ anek, D.; Louie, S. G.; Cohen, M. L.; Hybertsen, M. S. Quasiparticle calculation of valence band offset of AlAs-GaAs(001). Solid State Comms. 1988, 66, 585. Berens, P. H.; Mackay, D. H. J.; White, G. M.; Wilson, K. R. Thermodynamics and quantum corrections from molecular dynamics for liquid water. J. Chem. Phys. 1983, 79, 2375. Lin, S.-T.; Maiti, P. K.; Goddard III, W. A. Two-Phase Thermodynamic Model for Efficient and Accurate Absolute Entropy of Water from Molecular Dynamics Simulations. J. Phys. Chem. B 2010, 114, 8191. Lai, P.-K.; Hsieh, C.-M.; Lin, S.-T. Rapid determination of entropy and free energy of mixtures from molecular dynamics simulations with the two-phase thermodynamic model. Phys. Chem. Chem. Phys. 2012, 14, 15206. Pascal, T. A.; Goddard III, W. A. Hydrophobic segregation, phase transitions and the anomalous thermodynamics of water/methanol mixtures. J. Phys. Chem. B 2012, 116, 13905. Caro, M. A.; Laurila, T.; Lopez-Acevedo, O. Accurate schemes for calculation of thermodynamic properties of liquid mixtures from molecular dynamics simulations. J. Chem. Phys. 2016, 145, 244504. Zhang, C.; Spanu, L.; Galli, G. Entropy of Liquid Water from Ab Initio Molecular Dynamics. J. Phys. Chem. B 2011, 115, 14190. Cannes, C.; Kanoufi, F.; Bard, A. J. Cyclic voltammetry and scanning electrochemical microscopy of ferrocenemethanol at monolayer and bilayer-modified gold electrodes. J. Electroanal. Chem. 2003, 547, 83. Tsierkezos, N. G. Cyclic Voltammetric Studies of Ferrocene in Nonaqueous Solvents in the Temperature Range from

27

28

29 30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

248.15 to 298.15 K. J. Sol. Chem. 2007, 36, 289. Wang, Y.; Rogers, E. I.; Compton, R. G. The measurement of the diffusion coefficients of ferrocene and ferrocenium and their temperature dependence in acetonitrile using double potential step microdisk electrode chronoamperometry. J. Electroanal. Chem. 2010, 648, 15. Tsierkezos, N. G.; Ritter, U. Electrochemical impedance spectroscopy and cyclic voltammetry of ferrocene in acetonitrile/acetone system. J. Appl. Electrochem. 2010, 40, 409. http://dospt.org. Kresse, G.; Furthm¨ uller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B 1996, 54, 11169. Bl¨ ochl, P. E. Projector augmented-wave method. Phys. Rev. B 1994, 50, 17953. Kresse, G.; Joubert, D. From ultrasoft pseudopotentials to the projector augmented-wave method. Phys. Rev. B 1999, 59, 1758. Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865. Tkatchenko, A.; Scheffler, M. Accurate Molecular Van Der Waals Interactions from Ground-State Electron Density and Free-Atom Reference Data. Phys. Rev. Lett. 2009, 102, 073005. Bjelkmar, P.; Larsson, P.; Cuendet, M. A.; Hess, B.; Lindahl, E. Implementation of the CHARMM Force Field in GROMACS: Analysis of Protein Stability Effects from Correction Maps, Virtual Interaction Sites, and Water Models. J. Chem. Theory Comput. 2010, 6, 459. Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79, 926. Van Der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. GROMACS: Fast, flexible, and free. J. Comput. Chem. 2005, 26, 1701. Nos´e, S. A unified formulation of the constant temperature molecular dynamics methods. J. Chem. Phys. 1984, 81, 511. Hoover, W. G. Canonical dynamics: equilibrium phasespace distributions. Phys. Rev. A 1985, 31, 1695. Gillan, M. J.; Alf`e, D.; Michaelides, A. Perspective: How good is DFT for water? J. Chem. Phys. 2016, 144, 130901. Lee, C.; Yang, W.; Parr, R. G. Development of the ColleSalvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B 1988, 37, 785. Dudarev, S. L.; Botton, G. A.; Savrasov, S. Y.; Humphreys, C. J.; Sutton, A. P. Electron-energy-loss spectra and the structural stability of nickel oxide: An LSDA+U study. Phys. Rev. B 1998, 57, 1505. Caro, M. A. Solvation shells and radial distribution functions of transition metal complexes in aqueous solution: results from ab initio molecular dynamics. Zenodo 2017, Post, K.; Robins, R. G. Thermodynamic diagrams for the vanadium-water system at 298.15 K. Electrochim. Acta 1976, 21, 401. House, C. I.; Kelsall, G. H. Potential-pH diagrams for the Sn/H2 O–Cl system. Electrochim. acta 1984, 29, 1459. Pavlishchuk, V. V.; Addison, A. W. Conversion constants for redox potentials measured versus different reference electrodes in acetonitrile solutions at 25 C. Inorg. Chim. Acta 2000, 298, 97.



Page 12 of 13 12

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

47

48

Trasatti, S. The absolute electrode potential: an explanatory note (Recommendations 1986). Pure Appl. Chem. 1986, 58, 955. Gagne, R. R.; Koval, C. A.; Lisensky, G. C. Ferrocene as an

49

internal standard for electrochemical measurements. Inorg. Chem. 1980, 19, 2854. Heyd, J.; Scuseria, G. E.; Ernzerhof, M. Hybrid functionals based on a screened Coulomb potential. J. Chem. Phys. 2003, 118, 8207.


Page 13 of 13



Redox Potentials from Ab Initio Molecular Dynamics and Explicit

Recommend Documents