Next Article in Journal
First Principles Study of the Vibrational and Thermal Properties of Sn-Based Type II Clathrates, CsxSn136 (0 ≤ x ≤ 24) and Rb24Ga24Sn112
Previous Article in Journal
Photocracking Silica: Tuning the Plasmonic Photothermal Degradation of Mesoporous Silica Encapsulating Gold Nanoparticles for Cargo Release
Previous Article in Special Issue
Modeling the Catalyst Activation Step in a Metal–Ligand Radical Mechanism Based Water Oxidation System
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Determination of pKa Values via ab initio Molecular Dynamics and its Application to Transition Metal-Based Water Oxidation Catalysts

Department of Chemistry, University of Zurich, Winterthurerstrasse 190, CH-8057 Zurich, Switzerland
*
Author to whom correspondence should be addressed.
Submission received: 11 May 2019 / Revised: 31 May 2019 / Accepted: 6 June 2019 / Published: 12 June 2019
(This article belongs to the Special Issue Recent Advances in Water Oxidation Catalysis)

Abstract

:
The p K a values are important for the in-depth elucidation of catalytic processes, the computational determination of which has been challenging. The first simulation protocols employing ab initio molecular dynamics simulations to calculate p K a values appeared almost two decades ago. Since then several slightly different methods have been proposed. We compare the performance of various evaluation methods in order to determine the most reliable protocol when it comes to simulate p K a values of transition metal-based complexes, such as the here investigated Ru-based water oxidation catalysts. The latter are of high interest for sustainable solar-light driven water splitting, and understanding of the underlying reaction mechanism is crucial for their further development.

1. Introduction

Nowadays, transition metal-based catalysts are employed for a broad range of applications including pharmaceuticals, industrial scale synthesis and the development of renewable energy sources. Among the latter is the class of water-splitting catalysts that is subdivided into water oxidation catalysts which oxidize water and evolve molecular oxygen, and water reduction catalysts, which reduce protons in order to release molecular hydrogen. Analogously to nature’s photosynthesis process, the goal of those catalysts is either to oxidize water to molecular oxygen thereby generating reduction equivalents or to use the latter to form energy rich chemical bonds, for example in the form of molecular hydrogen.
An in-depth understanding of the underlying reaction mechanisms for those processes is crucial in order to further improve or rationally design novel catalysts. A prerequisite for mechanistic studies is the knowledge of the chemical speciation of the transition metal complex under catalytic conditions. When working under aqueous conditions, special attention has to be given to functional groups which might undergo protonation/deprotonation reactions. Unfortunately it is often difficult to determine physical properties such as acidity constants ( p K a ) of catalytic intermediates due to their elusive nature. In this study we apply density functional theory (DFT)-based molecular dynamics, so-called ab initio molecular dynamics (AIMD), to reliably determine p K a values of transition metal complexes used as catalysts for water oxidation. The here discussed methodologies have been successfully applied to a variety of compounds up to now, among them are organic molecules [1,2,3], amino acids and peptides [4,5] as well as aqua complexes of transition metals [6].
While all the studies mentioned above are based on the same protocol, namely the Bluemoon methodology, there are significant differences when it comes to the post-processing. In the following we will shortly introduce the general protocol, as well as the different flavors of post-processing. Then we compare those approaches for a set of benchmarking molecules, as well as for the system of interest—a ruthenium-based water oxidation catalyst [Ru(II)Py5OMe(OH2)]2+ (where Py5OMe = 6,6 -(methoxy(pyridin-2-yl)methylene)di-2,2 -bipyridine) [7]. The said system, and in particular the thermodynamics and kinetics of the water oxidation process, were already investigated in-depth by Luber and co-workers employing DFT simulations [7,8]. Those studies gave important insights with respect to the water oxidation mechanism and further helped to come up with some design guidelines on how to further improve those catalysts. However, the limitations imposed by the implicit solvation led to the desire for a more sophisticated description of the solvation shell. In this context, also this study serves as a mini-review and benchmark study to validate the methodology for p K a determination to be applied to the same or similar systems in the future.

2. Methodology

The calculation of p K a values is a common task in computational chemistry—it is therefore not surprising that there are many protocols available for various levels of theory [9]. The most common approach used for small molecules is based on the calculation of the difference in free energy between the protonated and the deprotonated species using a thermodynamic cycle scheme, whereby the interaction between the solvent and solute was approximated by an electrostatic potential by means of an implicit solvation model [10,11]. Alternatively, the electrostatic contribution of solvation can be obtained in a two-step procedure. In the first step, the point charge distribution was generated by the restricted electrostatic potential (RESP) procedure [12]. Then the Poisson equation is solved in order to obtain the electrostatic energies of solvation [13,14]. In combination with sampling of the conformational space, the second approach is in particular useful for large systems such as proteins which possess multiple protonation sides [15,16]. Even though many protocols are able to reliably reproduce experimental p K a values, their performance is still strongly system dependent. For example, large, flexible or highly charged species species stand in conflict with the underlying approximations of some of those protocols. There are numerous protocols and correction schemes depending on the system of interest which can be used to account for such shortcomings. A full review of them is beyond the scope of the current work, and we refer the interested reader to a number of selected articles [17,18,19,20,21,22,23,24,25].
An obvious approach to improve the previously presented protocol was to describe the solvent at an atomistic level. However, if both the solute and the solvent are treated explicitly using AIMD not only the computational cost rises drastically but also the simulation protocols become more elaborated [26,27,28,29,30,31,32].
Among those protocols was the so-called Bluemoon ensemble, where the free energy difference Δ F between the protonated and deprotonated state is calculated by a thermodynamic integration scheme. In this case the integrand is the average force ( f ξ ) acting on a system to impose the constraint ( ξ ), while the discrete values of the constraint ( ξ ) define the range of the integration:
Δ F = ξ 0 ξ 1 f ξ d ξ .
The latter is also referred to as potential of mean force (PMF). The average force f ξ is derived from the Lagrange multiplier λ according to
f ξ = Z 1 / 2 [ λ k B T G ] ξ Z 1 / 2 ξ ,
where k B is the Boltzmann constant, T the temperature, and Z and G are correction factors associated with the transformation from generalized to Cartesian coordinates. In the case of distance constraint (d(A−H)) Equation (2) simplifies to λ ξ = f ξ .
A detailed derivation of the Bluemoon methodology and in particular Equation (2) can be found in the corresponding original literature by Sprik and Ciccotti [33,34,35].
The application of the Bluemoon methodology to determine p K a values is in principle straight forward, however certain aspects deserve attention.

2.1. Choice of Constraint

From a chemical point of view the two states of interest, i.e., the protonated and deprotonated species, are well defined. Either the proton is bound to the acidic functional group (A) or it is (infinitely) far away from the latter stabilized by an extensive hydrogen bonding network. Simulating the deprotonation state in a chemical sense is restricted due to the limitations with regard to the size of the simulation cell. There are sophisticated proton insertion schemes which avoid this problem [3,26,27,28,29]. However their high demand in terms of computational resources renders this approach unsuitable for complex systems.
A crucial choice within the Bluemoon methodology is the nature of the constraint that describes the two states. The most simple one is the A−H distance (d(A−H)). While it is often applicable, it suffers from some intrinsic problems. First, this constraint does not prevent the re-protonation of the acidic group. Secondly, upon deprotonation, proton hopping might take place from the proton accepting molecule to other solvent molecules according to the Grotthuss mechanism. This is problematic since for large distances, in principle the intermolecular distance between the acidic group and a water molecule somewhere in solution is constraint. The latter is different from simulating a hydronium ion at infinite separation from the acidic group. Besides that, the important question arises at which distance is a covalent O−H bond broken, respectively, which distances have to be sampled in order to reach the deprotonated state. In particular proton hopping makes this decision ambiguous. Nevertheless there are several systems known for which the simple distance constrained turned out to work reasonably well [2,5,36].
Those issues might be circumvented by not only constraining a single A−H bond but the coordination number (CN) of all the protons to the acidic group [37]. The coordination number is commonly represented as the sum over sigmoid functions such as a Fermi–Dirac distribution. The latter has been successfully applied to several systems [1,3,31,32,38].
While a CN prevents the reprotonation of the acidic group it still does not prevent proton hopping. De Meyer et al. resolved this problem by constraining the difference in the CN of all the protons (index i, total number of protons: N) to the acidic group (index j) and a selected solvent molecule (index k) in its proximity:
C N ( r i j , r i k ) = i N 1 ( r i j r 0 ) n i N 1 ( r i j r 0 ) m i N 1 ( r i k r 0 ) n i N 1 ( r i k r 0 ) m ,
where r i j is the length of the vector r i j , describing the distance of proton i to the acidic group j, analogously r i k is the distance of proton i to the solvent molecule k. The inflection point is defined by r 0 , and exponential factors n and m define the overall shape of the switch function. This constraint in principle guaranties a smooth transition from the protonated acid to a solvent molecule without further proton hopping [3]. However, the said constraint introduces an additional empirical parameter, namely the choice of the proton accepting solvent molecule. Nonetheless, similar p K a values were obtained when the accepting molecule was either part of the first or second solvation shell [3].

2.2. Estimation of p K a Values from the Free Energy Differences

Independent of the constraint used for the simulation, the free energy differences obtained from the Bluemoon ensemble might be interpreted as equilibrium constants from which p K a values can be determined. There are several proposed ways, all of which have been shown to be able to reproduce experimental results reasonably well. However, to the best of our knowledge, those approaches have not be directly compared for the very same system. In the following we will present the different techniques, highlight their requirements in terms of simulations protocols and the necessary empiric parameters for the evaluation. For the detailed derivation of those methods, we refer the reader to the indicated references.

2.2.1. Absolute p K a

The most straight forward approach is to use the relation between the equilibrium constant p K a and the difference in free energy ( Δ F ):
p K a = β Δ F ln ( 10 ) ,
with β = 1 / k B T .
In order to reproduce experimental p K a values one has to assure that the difference in free energy sufficiently describes the two states. This boils down to the nature of the constraint as well as to the range of sampling the corresponding phase space. In any case a key requirement is the convergence of Δ F ( ξ ) to a constant value towards the limits of the sampled range ( ξ m i n to ξ m a x ), i.e., lim ξ ξ m a x Δ F ( ξ ) = constant. Further, in order to obtain meaningful results from the thermodynamic integration the free energy at the bound state F ( ξ 0 ) has to be set to unity. This method has been successfully applied to several systems using either the d(A−H)−d(O−H), the CN(A−H) or CN(A−H)−CN(O−H) (see Equation (3)) constraints [2,3,39].
Whether the calculated free energy difference Δ F can directly be related to the equilibrium constant is disputed in literature [4,6,31]. Authors who disagree with the previously discussed method commonly refer to Chandler’s derivations of the equilibrium constant based on a classical statistical mechanical description. The basic principle described there is the relation of the free energy difference Δ F and the radial distribution function (RDF) according to the reversible work theorem [40]. The RDF itself might be interpreted as the probability to find a proton within a certain radius of the acidic group. The probability distribution is related to the (inverse) acidity constant according to [31]:
K a 1 = c 0 0 e x p [ β Δ F ( r ) ] 4 π r 2 d r ,
where c 0 is the standard concentration. Note the difference in free energy Δ F ( r ) is a function of radius of a sphere around the acidic group and not of the constraint ξ . The latter are only equal in case of the distance constraint d(A−H). In principle Δ F ( r ) has to be known for infinite separation, in practice however only a finite separation R m a x L / 2 (L is the length of the cubic cell simulation) is accessible. Since Δ F ( r ) asymptotically approaches a constant value, one often defines R c R m a x , where R c is the radius which distinguishes A−H from A + H+, i.e., the distance at which the covalent bond is broken [31]. The limitations with respect to the simulation cell result in an uncertainty in the p K a value which Davies et al. quantified as Δ F ( R m a x ) / ( 2.3 k B T ) [31].
Based on the RDF there are two common approaches to derive p K a values, both of which require an additional set of simulations in order to reduce the potential errors describe above.

2.2.2. Relative p K a

The approach presented by Ivanov et al.—in the later referred to as “relative p K a ”—takes advantage of error-cancellation upon reporting the p K a value relative to a p K a value of a reference system (REF):
K a H A K a R E F = 0 R c e x p [ β Δ F R E F ( r ) ] r 2 d r 0 R c e x p [ β Δ F H A ( r ) ] r 2 d r .
Compared to Equation (5) the integration is from 0 to R c ( R c R m a x ) in order to account for the fact that lim r R m a x p K a ( r ) asymptotically approaches p K a ( R m a x ) , i.e., p K a ( r ) quickly becomes a constant value [4]. Further, the lower bound of the integral can be approximated by the value slightly smaller than the average A−H bond length, since the contributions of large values of Δ F ( r ) in the exponent are negligible [4]. The reference system (see Equation (6)) has to be calculated within the same computational framework, which in fact doubles the cost of the approach compared to the previously described method. However, the reference system might be used for the determination of p K a values of multiple species. Further in order to equate Equation (6), Δ F ( R m a x ) is set to unity for both the acid and the reference system [4,5].
This method has been successfully applied to several isomers of histidine as well as a histidine-tryptophan dimer employing a simple d(A−H) constraint [4,5].

2.2.3. Probabilistic p K a

The second approach was introduced by Davies et al., later referred to as “probabilistic p K a ” [31]. The main idea is to define the acid dissociation constant K a by the probabilities to find protons within a certain volume of the acidic group. Then the normalized probability to find a proton within a cutoff radius R c is
P ( R c ) = 0 R c e x p [ β Δ F ( r ) ] r 2 d r 0 R m a x e x p [ β Δ F ( r ) ] r 2 d r ,
where R c is the cutoff radius defining the protonated and the deprotonated state and R m a x the limit for infinite separation.
In the limit for weak acids the K a value then becomes
K a ( R c ) = ( 1 P ( R c ) ) 2 P ( R c ) N c 0 V ,
where c 0 is the standard concentration, V the volume of the simulation cell and N the number of acidic sides—here N = 1. The latter is converted to mol by dividing it by the Avogadro constant.
Here again the choice of the cutoff radius R c is crucial. However, unlike the relative p K a scheme, there is no obvious asymptotic behavior within the range of the deprotonation. This led Davies et al. to calculate the p K w of liquid water using the following relation:
K w ( R c ) = ( 1 P ( R c ) ) 2 N w c 0 V ,
where N w is the number of water molecules in the simulation cell. R c is the radius at which p K w ( r ) = 14 is true [31]. The relation between the Equation (8) and (9) is that the activity of the undissociated reactant is set to unity [31].
As for the previous methods, the current one has been successfully applied to several systems employing either a simple d(A−H) constraint or more commonly a CN constraint [1,6,31,32,38].
By introducing the three approaches named ”absolute”, ”relative” and ”probabilistic” p K a we have laid the foundation for the following study.

3. Computational Settings

AIMD simulations were performed employing the CP2K program package [41]. All atoms were described by the DZVP-MOLOPT-SR-GTH basis sets [42] as well as the corresponding GTH-BLYP pseudo potentials [43]. In order to enlarge the time step the mass of all hydrogen atoms was set to 2 a.m.u. in accordance with literature [3]. The influence of the latter is discussed in Section 4. The BLYP [44,45] exchange-correlation functional together with Grimme’s D3 dispersion correction [46], and a cutoff of 800 Ry for the auxiliary plane wave basis set were used.
In order to simulate liquid water, we used a cubic simulation cell with a side-length of 15.6404 Å containing 128 water molecules. This box size corresponds to liquid water at 1 bar and 300 K using the TIP5P force field [47]. The same simulation cell was used for solvated molecules, by deleting several water molecules we assured that the pressure remained approximately the same as the one of clean liquid water. For some model systems a larger simulation cell with side-length 19.7340 Å containing 256 water molecules was employed.
The simulations were performed in the NVT ensemble with a time-step of 0.5 fs. The temperature was kept constant at 320 K by a Nosé-Hoover chain thermostat [48,49]. The slightly elevated temperature is required in order to avoid the glassy behavior of BLYP water [50].
The general settings mentioned above closely resemble the protocols that have been employed previously by other groups in order to determine p K a values [3,4,5].
Each model system was equilibrated in the protonated state for 5–10 ps without any constraint. Starting from those simulations, constrained AIMD runs were performed each for an additional 15–20 ps (30,000–40,000 steps (see Supplementary Materials Tables S11 and S12)). The first 5 ps (10,000 steps) of each individual run were neglected in order to give the system time to equilibrate e.g., adopt to the imposed constraint. The convergence of each model system with respect to simulation time is given in the supporting information (see Supplementary Materials Tables S3–S5).

3.1. Model Systems

Simulations were conducted for the systems given in Table 1. p K a values were not only calculated for transition metal complexes, but also for two small organic molecules with p K a values in the same range as the molecules of interest. Thereby phenol serves as an internal standard which allows for direct comparison with p K a values obtained by the same methodology [3].
In our previous study we have investigated the mechanism of the water oxidation reaction catalyzed by [Ru(II)Py5Me(H2O)]2+ and [Ru(II)Py5OMe(H2O)]2+ listed in Table 1, from a kinetic and thermodynamic point of view employing state of the art DFT simulations [7,8]. We found that both ligand frameworks Py5OMe and Py5Me were virtually identical in terms of their thermodynamics and kinetics of the water-oxidation reaction, which is a not too surprising a result as the replacement of a methyl-group by a methoxy-group at an sp 3 carbon is not expected to remarkably alter electronics or sterics at the metal center. However, experimentally, the catalytic activity between the two ligands was found to be rather different, which was attributed to a rapid halide substitution at the catalyst bearing a Py5Me ligand. The latter leads to a deactivation of the catalyst [7]. As mechanistic studies for the more active catalyst Py5OMe are still underway, we decided to choose it as a model system even though no experimental p K a values are currently available. Based on our previous study we would not expect the thermodynamics of the two ligands and related p K a values to differ significantly [7,53]. Therefore comparing the calculated p K a value of the Py5OMe and Py5Me systems serves as a further validation of the method.

3.2. Error Analysis

The standard deviation of the average forces used to calculated the free energy difference according to Equation (1) is calculated by block averaging methods [54]. An upper limit for the standard deviation ( σ ) of the p K a is obtained by calculating the free energy Δ F of both, the average force ( λ ) and the average force plus its standard deviation i.e., ( λ + σ λ ) [36].

4. Results and Discussion

4.1. Convergence of the AIMD Simulations

For our comparison of the post-processing methods we ran constrained AIMDs for the systems described in Section 3.1. As a constraint we chose the distance of the acidic proton from the acidic group (d(A−H)), which we scanned in steps of 0.1 Å over a range from 0.9 Å to 1.6 Å. The constraint was primarily chosen for its simplicity, which avoids additional parameters, i.e., the explicit definition of the proton accepting molecule. Further, Equation (1) holds only for the distance constraint, otherwise correction terms are necessary (see Equation (2)) [33]. For most of the systems, proton hopping is observed for d(A−H) distances larger than 1.5 Å, which made the scanning of distances larger than 1.6 Å obsolete. In particular since large values of d(A−H) do not necessarily describe the distance between the acidic group and the initially formed hydronium, since the latter potentially loses its one of its proton to other solvent molecules.
An exemplary PMF profile of PhOH is shown in Figure 1. Error bars on the average forces were obtained by block averaging methods. The absolute value of those standard deviations is in the range of 0.5 to 1.3 kcal/mol, which corresponds to 0.3 to 0.9 p K a units. The PMF profiles of all the other model systems can be found in the Supplementary Materials Figures S2–S9.
The convergence of the force acting on the constraint can be seen in Figure 2. Each of the traces represents a single point in the top part of Figure 1. Towards larger values of the constraint the fluctuations start to oscillate around zero, highlighting the fact that the A−H bond has been broken.
The autocorrelation function of the force acting on the constraint is used to determine the optimal block size for the average which is in the range of 0.5 to 1 ps (see Figure 3). The force appears to be heavily correlated for the constraints ≤ 1.3 Å which corresponds to the region where the A−H bond is broken (see Figure 2). Ivanov et al. further analyzed those autocorrelation functions in order to elucidate the bond breaking process [4]. Since the thermostat modulates the autocorrelation function, simulations in the micro-canonical ensemble (NVE) would be required which are beyond the scope of the current work.

4.2. Reference System

In order to calculate the p K a value according to the “relative” or “probabilistic” method, the simulation of water is required. We determined a cutoff radius R c of 1.24 Å according to Equation (9) (see Figure 4). The same value was obtained for the two simulation cells containing either 128, or 256 water molecule (see Supplementary Materials Figure S17), which is in good agreement with 1.22 Å, respectively, 1.28 Å reported in literature [6,31,32]. As already reported by others, strong bases such as OH tend to be reprotonated by the solvent. This is also the case in our simulations. Starting from a d(O−H) of 1.5 Å, we were able to observe the reprotonation of the OH moiety (see Figure 5). As discussed earlier, choosing a different collective variable could circumvent this issue.
In the following we will present the p K a values calculated based on the protocol described beforehand.

4.3. Overview of Calculated p K a Values

In Table 2 the p K a values obtained with the three simulation protocols are given. In general we find that the “absolute” protocol (see Equation (5)) underestimates the p K a values compared to the experiment by about 1–2 p K a units. p K a values obtained by the relative protocol on the other hand are an overestimation by about two p K a units. The best agreement is achieved with the probabilistic protocol where in particular p K a values in the range of 10–11 are accurately reproduced. Hereby [Ru(III)Py5OMe(H2O)]3+ is somewhat of an outliner with the largest difference compared to the experiment by about two p K a units. As the only difference between [Ru(II)Py5OMe(H2O)]2+ and [Ru(III)Py5OMe(H2O)]3+ is the charge of the system, we reasoned that the simulation cell might be too small for such highly charged species. The latter was verified by a set of simulations in a bigger simulation cell (see Table 3).
The larger simulation cell slightly improves the agreement between the calculated p K a values using the relative protocol and the experiment. However, there is no systematic improvement of the p K a values with the system size, independent of the overall charge. Due to the significantly higher computational cost, we could not obtain the same level of convergence as compared to the smaller system. This can been seen from the shorter trajectories (see Supplementary Materials Tables S7 and S8).

4.4. Deuterated Solvent

The use of deuterated solute and solvent i.e., setting the mass of hydrogen atoms to 2 a.m.u may have a profound impact on calculated p K a values. In principle all values presented in Table 2 and Table 3 are p K a D values, i.e., p K a values in D2O instead of p K a H values, i.e., p K a values in H2O. Based on experiments, a correlation between p K a D and p K a H was reported already decades ago [55]. At first, the correlation was suspected to be linear only for p K a values > 7 and constant for more acidic species [55]. Later, Delgado et al. experimentally determined a linear relation between p K a D and p K a H over the whole range of p K a values [56]:
p K a D = 1.044 p K a H + 0.32 .
In a more recent study a detailed derivation of the linear relation between p K a D and p K a H has been presented by Krȩżel and Bal [57]:
p K a H = 0.929 p K a H + 0.41 ,
where p K a H is the p K a value determined in a D2O solution by a pH-meter which was calibrated by H2O. Conversion to p K a D is achieved by adding the empirically determined constant of 0.4 to p K a H [57]:
p K a D = p K a H + 0.4 .
Combining Equations (11) and (12) and solving for p K a D results in:
p K a D = 1.076 p K a H 0.041 ,
which is of the same mathematical form as Equation (10).
Before employing either Equation (10) or (13) to convert p K a D to p K a H values, some additional alterations to the post-processing procedure have to be made, since some of them include an explicit reference to H2O respectively its p K w value.
For the probabilistic method, the cut-off radius R c has to be redetermined since the p K w D of D2O is 14.951 (25 C) [58]. The obtained values for R c are 1.26 Å, respectively 1.25 Å for the two simulation cells (see Supplementary Materials Figures S18 and S19). The p K a D values obtained with the cut-off R c determined for D2O were converted to p K a H values (see Table 4 and Table 5).
The probabilistic p K a H values obtained by applying the conversion schemes discussed above are very similar to p K a values determined by referencing our calculations to H2O instead of D2O. The difference when converting the p K a D values either according to Equation (10) or Equation (13) is negligible.
For p K a D values, respectively p K a H values derived from the latter, calculated either by the absolute or relative method see Supplementary Materials Tables S7–S10. For the relative protocol, the overestimation p K a H values is reduced by about 1–2 units. However, no systematic improvement is achieved as in particular the p K a H values for the acidic compounds were still significantly overestimated. Further, the results for the large simulation cell deteriorate. The agreement between the experimental and calculated p K a H values using the absolute protocol decreased by 0.1 to 0.6 units upon converting p K a D to p K a H values.
Taking into account the multitude of empirical factors required to convert the p K a D values led to the conclusion that the referencing p K a values to H2O is acceptable. This is in particular true for the probabilistic method. The validity of this conclusion could in principle be checked by repeating all the simulations with H2O instead of D2O. With a large enough test set, it would also be possible to adjust equation of the linear relation to the employed methodology. However this is beyond the scope of the current work.
In the following we are going to highlight the dependence of the three methods on the cut-off radius ( R c ).

4.5. Absolute and Probabilistic p K a —Dependence on R c

In Figure 6 both the absolute as well as the probabilistic p K a values are shown as a function of the constraint (see Supplementary Materials Figures S11–S17 for plots for the other model systems). The choice of R c is obviously crucial (see Figure 6). In order to illustrate the influence of R c , p K a values were also calculated with R c values reported in literature (see Supplementary Materials Tables S1 and S2). Some systems are remarkably independent of R c while for others they spread over a range of almost 3 p K a units. This suggests that it is crucial to determine R c with exactly the same settings as the systems of interest, rather than to rely on an previously published value.

4.6. Relative p K a

As for the probabilistic protocol, the dependency of the relative p K a on R c was investigated. As expected, the p K a values calculated by Equation (6) show an asymptotic behavior towards larger values of R c [4] (see Figure 7). The latter suggests that large values of R c might be neglected in the evaluation. Indeed the results can be tuned i.e., lowered by 1 p K a unit, if only constraints between 0.9 and 1.4 Å are considered (see Supplementary Materials Table S6). Nonetheless the results for the acidic compounds are still overestimated with 5.8 and 6.2 for HCOOH, respectively [Ru(III)Py5OMe(H2O)]3+. Similar inconsistencies when applying the relative p K a protocol to lumiflavins have also been reported by Kiliç et al. [2].

5. Summary and Conclusions

In summary, we have investigated the Bluemoon methodology and its various post-processing methods in order to determine p K a values previously only applied to small organic molecules and transition metal aqua complexes. Our simulation cells with 128 respectively 256 water molecules were considerably larger than the ones previously presented in literature. This highlights the robustness and applicability of said protocol to larger chemically interesting systems, in particular in the context of water oxidation catalysts such as the here presented Ru(Py5Me) and Ru(Py5OMe) catalysts, for which we found, as expected, both qualitative and quantitative similar p K a values independent of the applied post-processing method.
When comparing the three post-processing methods suggested in literature (1) absolute, (2) relative and (3) probabilistic p K a (all referenced to H2O), we find that method 1 and 3 are able to quantitatively reproduce experimental values with an accuracy of about 1 p K a unit. In case of method 1 this is rather surprising as there is no guarantee that the free energy levels off within the scanned range of the constraint. Method 2 appears to at least qualitatively reproduce experimental p K a values, however there seems to be no necessity to use said approach, as it requires exactly the same set of calculations as method 3 which preformed the best within our test set. This conclusion also holds if one assumes that p K a D instead of p K a H values were calculated. For the sake of consistency we have only applied a very simple constraint i.e., the A−H distance. While this choice was fine for our test-set, there might be cases where more complex constraints are necessary, in particular if the conjugated base is very strong.
The overall accuracy of p K a values calculated by the Bluemoon methodology using a simple distance constraint lies between 1–3 p K a units. However if the computational settings are tailored in order to reproduce relevant experimental reference values, the accuracy might be increased. Overall it appears that both methods 2 and 3 slightly overestimate small p K a values while large p K a values are reproduced with high fidelity. This issue should be addressed in an extended benchmark study before one attempts to accurately predict low absolute p K a values.
Further we want to highlight the importance of a sufficiently large simulation cell when investigating highly charged systems such as water oxidation catalysts with different oxidation states. While the influence on the calculated p K a values was found to be minor, only a large enough simulation cell can guarantee that there are no spurious interactions between the solute and its mirror image in the neighboring simulation cells.

Supplementary Materials

The following are available online at https://0-www-mdpi-com.brum.beds.ac.uk/2304-6740/7/6/73/s1, Manuscript.pdf (Figures S1 (Lewis structure of [Ru(II)Py5R(H2O)]2+), Figures S2–S9 (Potentials of mean force as a function of the constraint - for all systems), Figures S10–S19 (Absolute and probabilistic p K a values as a function of the constraint - for all systems), Tables S1–S5 (Convergence of p K a with respect to simulation time), Tables S7–S10 ( p K a D values converted to p K a H values for the absolute and relative protocol), Tables S11 and S12 (Summary of simulation times per system/constraint), the following CP2K files: .inp, .ener, .xyz (only start structure), .LagrangeMultLog (forces) are available for each calculation.

Author Contributions

M.S. and S.L.; methodology, M.S.; software, M.S.; validation, M.S. and S.L.; formal analysis, M.S.; investigation, M.S.; resources, S.L.; writing—original draft preparation, M.S.; writing—review and editing, S.L.; visualization, M.S.; supervision, S.L.

Funding

The work has been supported by the University of Zurich Research Priority Program “Solar Light to Chemical Energy Conversion” (LightChEC) and the Swiss National Science Foundation (grant no. PP00P2_170667). We thank the Swiss National Supercomputing Center and the Platform for Advanced Computing in Europe for computing resources (project ID: s745 and s788).

Acknowledgments

Mauro Schilling would like to state his gratitude towards Richard A. Cunha for the fruitful discussions.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
MDPIMultidisciplinary Digital Publishing Institute
DOAJDirectory of open access journals
TLAThree letter acronym
LDLinear dichroism

References

  1. Chen, Y.L.; Doltsinis, N.L.; Hider, R.C.; Barlow, D.J. Prediction of Absolute Hydroxyl pKa Values for 3-Hydroxypyridin-4-ones. J. Phys. Chem. Lett. 2012, 3, 2980–2985. [Google Scholar] [CrossRef]
  2. Kiliç, M.; Ensing, B. Acidity constants of lumiflavin from first principles molecular dynamics simulations. Phys. Chem. Chem. Phys. 2014, 16, 18993–19000. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  3. De Meyer, T.; Ensing, B.; Rogge, S.M.J.; De Clerck, K.; Meijer, E.J.; Van Speybroeck, V. Acidity Constant (pKa) Calculation of Large Solvated Dye Molecules: Evaluation of Two Advanced Molecular Dynamics Methods. ChemPlusChem 2016, 17, 3447–3459. [Google Scholar] [CrossRef]
  4. Ivanov, I.; Chen, B.; Raugei, S.; Klein, M.L. Relative pKa Values from First-Principles Molecular Dynamics: The Case of Histidine Deprotonation. J. Phys. Chem. B 2006, 110, 6365–6371. [Google Scholar] [CrossRef] [PubMed]
  5. Bankura, A.; Klein, M.L.; Carnevale, V. Proton affinity of the histidine-tryptophan cluster motif from the influenza A virus from ab initio molecular dynamics. Chem. Phys. 2013, 422, 156–164. [Google Scholar] [CrossRef] [Green Version]
  6. Bernasconi, L.; Baerends, E.J.; Sprik, M. Long-Range Solvent Effects on the Orbital Interaction Mechanism of Water Acidity Enhancement in Metal Ion Solutions: A Comparative Study of the Electronic Structure of Aqueous Mg and Zn Dications. J. Phys. Chem. B 2006, 110, 11444–11453. [Google Scholar] [CrossRef] [PubMed]
  7. Gil-Sepulcre, M.; Böhler, M.; Schilling, M.; Bozoglian, F.; Bachmann, C.; Scherrer, D.; Fox, T.; Spingler, B.; Gimbert-Suriñach, C.; Alberto, R.; et al. Ruthenium Water Oxidation Catalysts based on Pentapyridyl Ligands. ChemSusChem 2017, 10, 4517–4525. [Google Scholar] [CrossRef] [Green Version]
  8. Schilling, M.; Böhler, M.; Luber, S. Towards the rational design of the Py5-ligand framework for ruthenium-based water oxidation catalysts. Dalton Trans. 2018, 47, 10480–10490. [Google Scholar] [CrossRef] [Green Version]
  9. Shields, G.S.; Seybold, P.G. CRC Handbook of Chemistry and Physics; CRC Press: Boca Raton, FL, USA; Taylor & Francis Group: Abingdon, UK, 2014. [Google Scholar]
  10. Tomasi, J. Thirty years of continuum solvation chemistry: A review, and prospects for the near future. Theor. Chem. Acc. 2004, 112, 184–203. [Google Scholar] [CrossRef]
  11. Tomasi, J.; Cancès, E.; Pomelli, C.S.; Caricato, M.; Scalmani, G.; Frisch, M.J.; Cammi, R.; Basilevsky, M.V.; Chuev, G.N.; Mennucci, B. Modern Theories of Continuum Models; John Wiley & Sons, Ltd.: Hoboken, NJ, USA, 2007. [Google Scholar]
  12. Bayly, C.I.; Cieplak, P.; Cornell, W.; Kollman, P.A. A well-behaved electrostatic potential based method using charge restraints for deriving atomic charges: The RESP model. J. Chem. Phys. 1993, 97, 10269–10280. [Google Scholar] [CrossRef]
  13. Am Busch, M.S.; Knapp, E.W. Accurate pKa Determination for a Heterogeneous Group of Organic Molecules. ChemPlusChem 2004, 5, 1513–1522. [Google Scholar] [CrossRef]
  14. Galstyan, G.; Knapp, E.W. Computing pKa Values of Hexa-Aqua Transition Metal Complexes. J. Comput. Chem. 2015, 36, 69–78. [Google Scholar] [CrossRef] [PubMed]
  15. Simonson, T.; Carlsson, J.; Case, D.A. Proton Binding to Proteins: pKa Calculations with Explicit and Implicit Solvent Models. J. Am. Chem. Soc. 2004, 126, 4167–4180. [Google Scholar] [CrossRef] [PubMed]
  16. Meyer, T.; Kieseritzky, G.; Knapp, E.W. Electrostatic pKa computations in proteins: Role of internal cavities. Proteins 2011, 79, 3320–3332. [Google Scholar] [CrossRef] [PubMed]
  17. Cramer, C.J.; Truhlar, D.G. Implicit Solvation Models: Equilibria, Structure, Spectra, and Dynamics. Chem. Rev. 1999, 99, 2161–2200. [Google Scholar] [CrossRef] [PubMed]
  18. Ho, J.; Coote, M.L. A universal approach for continuum solvent pKa calculations: Are we there yet? Theor. Chem. Acc. 2009, 125, 3. [Google Scholar] [CrossRef]
  19. Zhan, C.G.; Dixon, D.A. Absolute Hydration Free Energy of the Proton from First-Principles Electronic Structure Calculations. J. Phys. Chem. A 2001, 105, 11534–11540. [Google Scholar] [CrossRef]
  20. Klamt, A.; Eckert, F.; Diedenhofen, M.; Beck, M.E. First Principles Calculations of Aqueous pKa Values for Organic and Inorganic Acids Using COSMO-RS Reveal an Inconsistency in the Slope of the pKa Scale. J. Phys. Chem. A 2003, 107, 9380–9386. [Google Scholar] [CrossRef]
  21. Nielsen, J.E.; Gunner, M.R.; García-Moreno, E.B. The pKa Cooperative: A collaborative effort to advance structure-based calculations of pKa values and electrostatic effects in proteins. Proteins 2011, 79, 3249–3259. [Google Scholar] [CrossRef]
  22. Alexov, E.; Mehler, E.L.; Baker, N.M.; Baptista, A.; Huang, Y.; Milletti, F.; Erik Nielsen, J.; Farrell, D.; Carstensen, T.; Olsson, M.H.M.; et al. Progress in the prediction of pKa values in proteins. Proteins 2011, 79, 3260–3275. [Google Scholar] [CrossRef]
  23. Marenich, A.V.; Ho, J.; Coote, M.L.; Cramer, C.J.; Truhlar, D.G. Computational electrochemistry: Prediction of liquid-phase reduction potentials. Phys. Chem. Chem. Phys. 2014, 16, 15068–15106. [Google Scholar] [CrossRef] [PubMed]
  24. Ho, J. Are thermodynamic cycles necessary for continuum solvent calculation of pKas and reduction potentials? Phys. Chem. Chem. Phys. 2015, 17, 2859–2868. [Google Scholar] [CrossRef] [PubMed]
  25. Gunner, M.; Baker, N. Chapter One—Continuum Electrostatics Approaches to Calculating pKas and Ems in Proteins. In Computational Approaches for Studying Enzyme Mechanism Part B; Voth, G.A., Ed.; Academic Press: Cambridge, MA, USA, 2016; Volume 578, pp. 1–20. [Google Scholar] [CrossRef]
  26. Sulpizi, M.; Sprik, M. Acidity constants from vertical energy gaps: Density functional theory based molecular dynamics implementation. Phys. Chem. Chem. Phys. 2008, 10, 5238–5249. [Google Scholar] [CrossRef] [PubMed]
  27. Cheng, J.; Sulpizi, M.; Sprik, M. Redox potentials and pKa for benzoquinone from density functional theory based molecular dynamics. J. Chem. Phys. 2009, 131, 154504. [Google Scholar] [CrossRef] [PubMed]
  28. Sulpizi, M.; Sprik, M. Acidity constants from DFT-based molecular dynamics simulations. J. Phys. Condens. Matter. 2010, 22, 284116. [Google Scholar] [CrossRef] [PubMed]
  29. 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–3529. [Google Scholar] [CrossRef] [Green Version]
  30. Sprik, M. Computation of the pK of liquid water using coordination constraints. Chem. Phys. 2000, 258, 139–150. [Google Scholar] [CrossRef]
  31. Davies, J.E.; Doltsinis, N.L.; Kirby, A.J.; Roussev, C.D.; Sprik, M. Estimating pKa Values for Pentaoxyphosphoranes. J. Am. Chem. Soc. 2002, 124, 6594–6599. [Google Scholar] [CrossRef]
  32. Doltsinis, N.L.; Sprik, M. Theoretical pKa estimates for solvated P(OH)5 from coordination constrained Car-Parrinello molecular dynamics. Phys. Chem. Chem. Phys. 2003, 5, 2612–2618. [Google Scholar] [CrossRef]
  33. Sprik, M.; Ciccotti, G. Free energy from constrained molecular dynamics. J. Chem. Phys. 1998, 109, 7737–7744. [Google Scholar] [CrossRef]
  34. Ciccotti, G.; Ferrario, M. Blue Moon Approach to Rare Events. Mol. Sim. 2004, 30, 787–793. [Google Scholar] [CrossRef]
  35. Ciccotti, G.; Kapral, R.; Vanden-Eijnden, E. Blue Moon Sampling, Vectorial Reaction Coordinates, and Unbiased Constrained Dynamics. ChemPlusChem 2005, 6, 1809–1814. [Google Scholar] [CrossRef]
  36. Brüssel, M.; Di Dio, P.J.; Muñiz, K.; Kirchner, B. Comparison of Free Energy Surfaces Calculations from Ab Initio Molecular Dynamic Simulations at the Example of Two Transition Metal Catalyzed Reactions. Int. J. Mol. Sci. 2011, 12, 1389–1409. [Google Scholar] [CrossRef] [PubMed]
  37. Sprik, M. Coordination numbers as reaction coordinates in constrained molecular dynamics. Faraday Discuss. 1998, 110, 437–445. [Google Scholar] [CrossRef]
  38. Liu, X.; Lu, X.; Wang, R.; Zhou, H. In Silico Calculation of Acidity Constants of Carbonic Acid Conformers. J. Phys. Chem. A 2010, 114, 12914–12917. [Google Scholar] [CrossRef] [PubMed]
  39. Sinha, V.; Govindarajan, N.; de Bruin, B.; Meijer, E.J. How Solvent Affects C–H Activation and Hydrogen Production Pathways in Homogeneous Ru-Catalyzed Methanol Dehydrogenation Reactions. ACS Catal. 2018, 8, 6908–6913. [Google Scholar] [CrossRef]
  40. Chandler, D. Introduction to Modern Statistical Mechanics; Oxford University Press, Inc.: Oxford, UK, 1987. [Google Scholar]
  41. CP2K Developers Group. CP2K Program Package. Available online: https://www.cp2k.org/ (accessed on 11 June 2019).
  42. VandeVondele, J.; Hutter, J. Gaussian basis sets for accurate calculations on molecular systems in gas and condensed phases. J. Chem. Phys. 2007, 127, 114105. [Google Scholar] [CrossRef] [Green Version]
  43. Goedecker, S.; Teter, M.; Hutter, J. Separable dual-space Gaussian pseudopotentials. Phys. Rev. B 1996, 54, 1703–1710. [Google Scholar] [CrossRef] [Green Version]
  44. Becke, A.D. Density-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. A 1988, 38, 3098–3100. [Google Scholar] [CrossRef]
  45. Lee, C.; Yang, W.; Parr, R.G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B 1988, 37, 785–789. [Google Scholar] [CrossRef] [Green Version]
  46. Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H–Pu. J. Chem. Phys. 2010, 132, 154104. [Google Scholar] [CrossRef] [PubMed]
  47. Mahoney, M.W.; Jorgensen, W.L. A five-site model for liquid water and the reproduction of the density anomaly by rigid, nonpolarizable potential functions. J. Chem. Phys. 2000, 112, 8910–8922. [Google Scholar] [CrossRef]
  48. Nosé, S. A unified formulation of the constant temperature molecular dynamics methods. J. Chem. Phys. 1984, 81, 511–519. [Google Scholar] [CrossRef] [Green Version]
  49. Nosé, S. A molecular dynamics method for simulations in the canonical ensemble. Mol. Phys. 1984, 52, 255–268. [Google Scholar] [CrossRef]
  50. VandeVondele, J.; Mohamed, F.; Krack, M.; Hutter, J.; Sprik, M.; Parrinello, M. The influence of temperature and density functional models in ab initio molecular dynamics simulation of liquid water. J. Chem. Phys. 2005, 122, 014515. [Google Scholar] [CrossRef] [PubMed]
  51. Gil-Sepulcre, M.; Axelson, J.C.; Aguiló, J.; Solà-Hernández, L.; Francàs, L.; Poater, A.; Blancafort, L.; Benet-Buchholz, J.; Guirado, G.; Escriche, L.; et al. Synthesis and Isomeric Analysis of Ru(II) Complexes Bearing Pentadentate Scaffolds. Inorg. Chem. 2016, 55, 11216–11229. [Google Scholar] [CrossRef] [PubMed]
  52. Braude, A.E.; Nachod, F.C. Determination of Organic Structures by Physical Methods; Academic Press., Inc.: Cambridge, MA, USA, 1955; Volume 1. [Google Scholar]
  53. Schilling, M.; Luber, S. Computational Modeling of Cobalt-Based Water Oxidation: Current Status and Future Challenges. Front Chem. 2018, 6, 100. [Google Scholar] [CrossRef] [Green Version]
  54. Flyvbjerg, H.; Petersen, H.G. Error estimates on averages of correlated data. J. Chem. Phys. 1989, 91, 461–466. [Google Scholar] [CrossRef]
  55. Robinson, R.A.; Paabo, M.; Bates, R.G. Deuterium Isotope Effect on the Dissociation of Weak Acids in Water and Deuterium Oxide. J. Res. Natl. Bur. Stand. 1969, 73A, 299–308. [Google Scholar] [CrossRef]
  56. Delgado, R.; Silva, J.D.; Amorim, M.; Cabral, M.; Chaves, S.; Costa, J. Dissociation constants of Brønsted acids in D2O and H2O: Studies on polyaza and polyoxa-polyaza macrocycles and a general correlation. Anal. Chim. Acta 1991, 245, 271–282. [Google Scholar] [CrossRef]
  57. Krȩżel, A.; Bal, W. A formula for correlating pKa values determined in D2O and H2O. J. Inorg. Biochem. 2004, 98, 161–166. [Google Scholar] [CrossRef] [PubMed]
  58. Lide, D.R. CRC Handbook of Chemistry and Physics; Number 8–85; CRC Press: Boca Raton, FL, USA, 2005. [Google Scholar]
Figure 1. (top): Average force ( λ ) acting on the constraint for the PhOH model system, constraining the A−H bond. (bottom): Potential of mean force, i.e., free energy obtained by integrating the average forces.
Figure 1. (top): Average force ( λ ) acting on the constraint for the PhOH model system, constraining the A−H bond. (bottom): Potential of mean force, i.e., free energy obtained by integrating the average forces.
Inorganics 07 00073 g001
Figure 2. Force ( λ ) acting on the constraint—here for PhOH at a constrained d(A−H). The vertical line at 5.0 ps marks the equilibration time in relation to the production run.
Figure 2. Force ( λ ) acting on the constraint—here for PhOH at a constrained d(A−H). The vertical line at 5.0 ps marks the equilibration time in relation to the production run.
Inorganics 07 00073 g002
Figure 3. Autocorrelation function of the force acting on the constraint. From top to bottom the constraint increases by 0.1 Å, starting from 0.9 Å.
Figure 3. Autocorrelation function of the force acting on the constraint. From top to bottom the constraint increases by 0.1 Å, starting from 0.9 Å.
Inorganics 07 00073 g003
Figure 4. Determination of R c from simulation of water, by fitting to the experimental value. The simulations were carried out in a cubic box with a side length of 15.6406 Å.
Figure 4. Determination of R c from simulation of water, by fitting to the experimental value. The simulations were carried out in a cubic box with a side length of 15.6406 Å.
Inorganics 07 00073 g004
Figure 5. Autodissociation of H2O (top): Average force acting on the constraint. (bottom): Free energy obtained by integrating the average forces. Note the drop in free energy at 1.5 Å is caused by the reprotonation of OH.
Figure 5. Autodissociation of H2O (top): Average force acting on the constraint. (bottom): Free energy obtained by integrating the average forces. Note the drop in free energy at 1.5 Å is caused by the reprotonation of OH.
Inorganics 07 00073 g005
Figure 6. p K a values calculated for PhOH using the absolute and probabilistic method. The latter is strongly dependent on the choice of R c , even a change by 0.01 Å might result in change of up to 0.5 p K a unit.
Figure 6. p K a values calculated for PhOH using the absolute and probabilistic method. The latter is strongly dependent on the choice of R c , even a change by 0.01 Å might result in change of up to 0.5 p K a unit.
Inorganics 07 00073 g006
Figure 7. Relative p K a as a function of d(A−H) according to Equation (6), the value asymptotically approaches a constant value.
Figure 7. Relative p K a as a function of d(A−H) according to Equation (6), the value asymptotically approaches a constant value.
Inorganics 07 00073 g007
Table 1. Model systems used in this study. N w stands for the number of water molecules in the simulation cell. The calculated p K a value is only given if the simulations found in the literature were obtained employing the Bluemoon methodology. The Ru complex bears a pentapyridine (Py5) ligand that is composed of two bipyridyl fragments linked to a fifth pyridyl via an sp 3 carbon. The fourth fragment connected to the latter is either a methyl (Me) or methoxy (OMe) group resulting in Py5Me or Py5OMe (see [7,8,51] for more details). Experimental p K a values denoted with “*” were only available for the Py5Me ligand framework (see Supplementary Materials Figure S1 for a graphical representation of the catalysts).
Table 1. Model systems used in this study. N w stands for the number of water molecules in the simulation cell. The calculated p K a value is only given if the simulations found in the literature were obtained employing the Bluemoon methodology. The Ru complex bears a pentapyridine (Py5) ligand that is composed of two bipyridyl fragments linked to a fifth pyridyl via an sp 3 carbon. The fourth fragment connected to the latter is either a methyl (Me) or methoxy (OMe) group resulting in Py5Me or Py5OMe (see [7,8,51] for more details). Experimental p K a values denoted with “*” were only available for the Py5Me ligand framework (see Supplementary Materials Figure S1 for a graphical representation of the catalysts).
Molecule N w Side-Length [Å] p K a (exp.) p K a (calc.)
H2O12815.614.0
H2O25619.714.0
HCOOH12615.63.8 [52]
PhOH12315.610.0 [52]9.7 [3]
[Ru(II)Py5Me(H2O)]2+11215.6∼11 [7]
[Ru(II)Py5OMe(H2O)]2+11215.6∼11 [7] *
[Ru(II)Py5OMe(H2O)]2+23419.7∼11 [7] *
[Ru(II)Py5OMe(H2O)]3+11215.6∼3 [7] *
[Ru(II)Py5OMe(H2O)]3+23419.7∼3 [7] *
Table 2. All results presented here are calculated at 320 K, for a cubic box with side length of 15.6406 Å and a cutoff R C = 1.24 Å, over a trajectory of about 20 ps (where the first 5 ps were not included in the evaluation). The standard deviation is calculated using the block average method with a block size of 1 ps; (a) absolute, (b) relative, (c) probabilistic protocols.
Table 2. All results presented here are calculated at 320 K, for a cubic box with side length of 15.6406 Å and a cutoff R C = 1.24 Å, over a trajectory of about 20 ps (where the first 5 ps were not included in the evaluation). The standard deviation is calculated using the block average method with a block size of 1 ps; (a) absolute, (b) relative, (c) probabilistic protocols.
Molecule p K a (exp.) p K a ( a ) p K a ( b ) p K a ( c )
H2O14.0
HCOOH3.8 [52] 2.7 ± 0.5 6.7 ± 0.5 4.2 ± 0.6
PhOH10.0 [52] 8.7 ± 0.3 12.5 ± 0.5 10.7 ± 0.4
[Ru(II)Py5Me(H2O)]2+∼11 [51] 9.8 ± 0.4 13.7 ± 0.3 11.2 ± 0.4
[Ru(II)Py5OMe(H2O)]2+∼11 [51] 9.3 ± 0.4 13.3 ± 0.6 11.1 ± 0.4
[Ru(III)Py5OMe(H2O)]3+∼2.5 [51] 3.1 ± 0.4 7.1 ± 0.4 4.5 ± 0.5
Table 3. All results presented here are calculated at 320 K, for a cubic box with a side length of 19.7340 Å and a cutoff R C = 1.24 Å over a trajectory of 10–15 ps (where the first 5 ps were not included in the evaluation). The standard deviation is calculated using the block average method with a block size of 1 ps; (a) absolute, (b) relative, (c) probability protocols.
Table 3. All results presented here are calculated at 320 K, for a cubic box with a side length of 19.7340 Å and a cutoff R C = 1.24 Å over a trajectory of 10–15 ps (where the first 5 ps were not included in the evaluation). The standard deviation is calculated using the block average method with a block size of 1 ps; (a) absolute, (b) relative, (c) probability protocols.
Molecule p K a (exp.) p K a ( a ) p K a ( b ) p K a ( c )
H2O14.0
[Ru(II)Py5OMe(H2O)]2+∼11 [51] 10.1 ± 0.5 11.5 ± 0.5 12.7 ± 0.7
[Ru(III)Py5OMe(H2O)]3+∼2.5 [51] 3.1 ± 0.3 4.6 ± 0.3 4.3 ± 0.4
Table 4. All results presented here are calculated at 320 K, for a cubic box with a side length of 15.6406 Å, over a trajectory of about 20 ps (where the first 5 ps were not included in the evaluation). The p K a values were calculated using the probabilistic method for a R c value of 1.26 Å determined for D2O. The p K a H (a) was obtained referencing the calculations to H2O i.e., a R c of 1.24 Å. The p K a D values were converted to p K a H values according to Equation (10) (b), respectively Equation (13) (c).
Table 4. All results presented here are calculated at 320 K, for a cubic box with a side length of 15.6406 Å, over a trajectory of about 20 ps (where the first 5 ps were not included in the evaluation). The p K a values were calculated using the probabilistic method for a R c value of 1.26 Å determined for D2O. The p K a H (a) was obtained referencing the calculations to H2O i.e., a R c of 1.24 Å. The p K a D values were converted to p K a H values according to Equation (10) (b), respectively Equation (13) (c).
Molecule p K a (exp.) p K a H (a) p K a D p K a H (b) p K a H (c)
H2O14.0
HCOOH3.8 [52]4.24.33.84.0
PhOH10.0 [52]10.711.210.410.5
[Ru(II)Py5Me(H2O)]2+∼11 [51]11.212.011.211.2
[Ru(II)Py5OMe(H2O)]2+∼11 [51]11.111.811.011.0
[Ru(III)Py5OMe(H2O)]3+∼2.5 [51]4.54.64.14.3
Table 5. All results presented here are calculated at 320 K, for a cubic box with a side length of 19.7340 Å, over a trajectory of about 20 ps (where the first 5 ps were not included in the evaluation). The p K a values were calculated using the probabilistic method for a R c value of 1.25 Å determined for D2O. The p K a H (a) was obtained referencing the calculations to H2O i.e., a R c of 1.24 Å. The p K a D values were converted to p K a H values according to Equation 10 (b), respectively Equation (13) (c).
Table 5. All results presented here are calculated at 320 K, for a cubic box with a side length of 19.7340 Å, over a trajectory of about 20 ps (where the first 5 ps were not included in the evaluation). The p K a values were calculated using the probabilistic method for a R c value of 1.25 Å determined for D2O. The p K a H (a) was obtained referencing the calculations to H2O i.e., a R c of 1.24 Å. The p K a D values were converted to p K a H values according to Equation 10 (b), respectively Equation (13) (c).
Molecule p K a (exp.) p K a H (a) p K a D p K a H (b) p K a H (c)
H2O14.0----
[Ru(II)Py5OMe(H2O)]2+∼11 [51]12.713.112.212.2
[Ru(III)Py5OMe(H2O)]3+∼2.5 [51]4.34.43.94.1

Share and Cite

MDPI and ACS Style

Schilling, M.; Luber, S. Determination of pKa Values via ab initio Molecular Dynamics and its Application to Transition Metal-Based Water Oxidation Catalysts. Inorganics 2019, 7, 73. https://0-doi-org.brum.beds.ac.uk/10.3390/inorganics7060073

AMA Style

Schilling M, Luber S. Determination of pKa Values via ab initio Molecular Dynamics and its Application to Transition Metal-Based Water Oxidation Catalysts. Inorganics. 2019; 7(6):73. https://0-doi-org.brum.beds.ac.uk/10.3390/inorganics7060073

Chicago/Turabian Style

Schilling, Mauro, and Sandra Luber. 2019. "Determination of pKa Values via ab initio Molecular Dynamics and its Application to Transition Metal-Based Water Oxidation Catalysts" Inorganics 7, no. 6: 73. https://0-doi-org.brum.beds.ac.uk/10.3390/inorganics7060073

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop