Dynamics with Explicit Solvation Reveals Formation of the Prereactive

Aug 7, 2018 - (40) We would like to note that this diffusion is much faster compared to the diffusion of the pyrene substituted [RuVO(bda)L2]+ catalys...
0 downloads 0 Views 689KB Size
Subscriber access provided by UNIVERSITY OF TOLEDO LIBRARIES

Article

Dynamics with explicit solvation reveals formation of the prereactive dimer as sole determining factor for the efficiency of Ru(bda)L2 catalysts Shaoqi Zhan, Rongfeng Zou, and Marten S. G. Ahlquist ACS Catal., Just Accepted Manuscript • DOI: 10.1021/acscatal.8b02519 • Publication Date (Web): 07 Aug 2018 Downloaded from http://pubs.acs.org on August 7, 2018

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 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 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.

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 23 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 Catalysis

Dynamics with explicit solvation reveals formation of the prereactive dimer as sole determining factor for the efficiency of Ru(bda)L2 catalysts Shaoqi Zhan, Rongfeng Zou, Mårten S. G. Ahlquist* Department of Theoretical Chemistry & Biology, School of Engineering Sciences in Chemistry Biotechnology and Health, KTH Royal Institute of Technology, 10691 Stockholm, Sweden.

ABSTRACT: This report describes all key steps in the O-O bond formation from two separated [RuV=O(bda)L2]+ cations to form the dinuclear [(bda)L2RuIV-O-O-RuIV(bda)L2]2+ in explicit solvent. The three steps involve the diffusion of the catalysts in the water phase, formation of the prereactive dimer, and the bond formation between the two catalysts. Based on the calculated parameters, we compute the rate constant of two catalysts with different L-ligands, isoquinoline and picoline, and the computed values are in excellent agreement with the experimental ones. The interaction of the axial ligands is key to the improved rates of the larger ligand, mainly by facilitating the formation of the prereactive dimer from the solvated monomer. By comparing the binding free energy of hydrophilic RuIV-OH and hydrophobic RuV=O, the hydrophobic driving force of RuV=O in this system has been estimated to 1 kcal mol-1.

ACS Paragon Plus Environment

1

ACS Catalysis 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

Page 2 of 23

KEYWORDS: water oxidation, binding free energy, diffusion rate, O-O bond formation, rate constant, molecular dynamics, empirical valence bond

INTRODUCTION Water oxidation (2H2O → O2 + 4e- + 4H+) is one of the key reactions in artificial photosynthesis with requirements for both four electron transfer reactions from two water molecules and O-O bond formation.1,2 Among the huge variety of catalytic systems for water oxidation,3-9 the Ru(bda)L2 (bda = 2,2’-bipyridine-6,6’-dicarboxylate, L = typically nitrogen containing heterocycle e.g. pyridine) type catalysts, using CeIV as driving force, stand out in activity, with values surpassing the activity as photosystem II.10-13 Unlike most mononuclear catalysts where the O-O bond is formed via the water nucleophilic attack (WNA) mechanism,14,15 the Ru(bda)L2 catalysts have been proposed to react via an interaction of two metal centers (I2M) mechanism.16,17 This mechanism involves a bimolecular coupling of two RuV=O radicals,11,18 where the key RuV=O intermediate was spectroscopically recently confirmed (Figure 1).19 The bimolecular coupling has been shown experimentally to work even at extremely low concentrations, down to 100 nM, which indicates extremely low activation energies for the O-O bond formation. In our first study, which included explicit solvation, we showed that the oxo fragment in the RuV=O species is hydrophobic, which explained the tendency of the complex for forming of the prereactive dimer with the oxos pointing at each other.20 From the prereactive dimer we found that water molecules do not interfere with the reaction, and the activation energy was entrirely due to small rearrangements of the prereactive dimer to the transition state, while the two oxos have a slightly attractive interaction.21 The kinetic measurements of the complexes at low concentration showed that the rate limiting step is

ACS Paragon Plus Environment

2

Page 3 of 23 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 Catalysis

associated with the coupling step and the second order rate constants are 7.83 × 105 M-1s-1 of complex 1 and 2.39 × 108 M-1s-1 of complex 2 (Figure 2).11,18 These results show that the choice of axial ligand has a huge effect on the reactivity, where bicyclic aromatic ligands such as isoquinoline (isoq) showed extreme reactivity, which was attributed to their π-stacking ability.22 In our initial study we found no correlation between the activation free energy of the O-O bond forming step alone, with the measures activities. This observation in combination with the fact that second order kinetics is observed with respect to the Ru concentration suggests that the formation of the prereactive dimer needs careful consideration for a complete description of the reaction.

Figure 1. Proposed the full reaction of the O-O bond formation, including diffusion, collision and binding, and O-O bond coupling in the water phase. As illustrated in Figure 1, the O-O bond formation between two [RuV=O(bda)L2]+ complexes is a combination of three components 1) random diffusion 2) occasional encountering between two catalysts forming the prereactive dimer and finally 3) the actual O-O coupling step. Until this year studies have only focused on the O-O coupling step,11,15,16,21,23 and no study give even qualitative agreement with the experimental rate constants from the free solvated monomers. The

ACS Paragon Plus Environment

3

ACS Catalysis 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

Page 4 of 23

diffusion and formation of prereactive dimer from [RuV=O(bda)L2]+ monomer is still not well understood. In a recent report we studied the reaction between two catalysts on a carbon nanotube-water interface.24 The approach to model all the three components of the reaction was to 1) study the diffusion of a monomer at the surfaces using classical molecular dynamics (MD) and from that extract the diffusion coefficient 2) use unbiased MD to let two catalysts form a prereactive complex followed by potential of mean force (PMF) simulations to extract the interaction free energy and 3) study the reaction at the interface using empirical valence bond (EVB) simulation.25 We found that the reaction could be studied with high accuracy by using the calculated diffusion coefficient for estimation of the collision frequency of the catalyst, and the activation free energy with the binding free energy subtracted as the total activation energy. The collision frequency was assumed to be equal to the Arrhenius prefactor (A) (at a given concentration), and the total activation energy equal to the Arrhenius activation energy Ea. This allowed us to describe a reaction of very high complexity with very high accuracy. Herein, we extend our studies to describe the full reaction of the rate-determining step, O-O bond formation from two separated [RuV=O(bda)L2]+ cations to form the dinuclear [(bda)L2RuIV-O-ORuIV(bda)L2]2+ in explicit solvent. We have focused on the two complexes 1 and 2 since accurate data is available for the rate constants of the O-O bond formation, where other steps of the catalysis do not interfere significantly with the measurements. Herein, we therefore compare directly to the second order rate constants of the reaction of two RuV=O complexes, rather than attempt to reproduce turnover frequencies (TOF), since TOFs often are measured at conditions where other steps could limit the rate.

ACS Paragon Plus Environment

4

Page 5 of 23 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 Catalysis

Figure 2. Complexes RuV(bda)(pic)2 1 and RuV(bda)(isoq)2 2 included in this study. COMPUTATIONAL DETAILS All DFT calculations for the estimation of Gibbs free energies were carried out with the Jaguar 8.3 program package by Schrödinger LLC.26 Molecular geometries were optimized using Becke’s three-parameter hybrid functional and the LYP correlation functional (B3LYP)27 with D3 correction of Grimme et al.28,29 with the LACVP** basis set.30 To identify the transition states for O-O bond formation, we searched the potential energy surface by scanning the terminal O-O bond distance [RuV=O•••O=RuV] of the antiferromagnetic open shell singlet. The thermochemical corrections for estimations of the Gibbs free energy barrier from free separated monomers (RuV=O) were calculated by at the B3LYP-D3/LACVP** level for both the monomer and transition state structures. The Gibbs free energy were defined as the following equation G = E(B3LYP-D3/LACVP**) + ZPE + H298 - TS298. The geometries of the complexes reactant monomer and product dimer have been optimized with the quantum mechanics (QM) method described in the previous section. Force field parameters of the complexes are from previous work20 to augment with standard OPLS-AA (allatom optimized molecular potential for liquid simulation)31 bonded and van der Waals

ACS Paragon Plus Environment

5

ACS Catalysis 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

Page 6 of 23

parameters for the picoline groups. Partial charge parameters were calculated in the same way using ESP charges as previously described.20 In order to run molecular dynamics (MD) simulations, the parameters were then transformed into the Gromacs topology format, where the Fourier coefficients of the dihedral potential term were transformed into the Ryckaerd-Bellemans type.32 MD simulations were performed with the Gromacs 5.1.2 MD software package.33 Firstly, we performed a 100 ns MD simulation in vacuum, to test the stability. Then we performed a 100 ns MD run in a 42 × 43 × 54 Å3 periodic box filled with TIP3P34 water molecules and a chloride ion to neutralize the charge. In MD simulations, the resulting systems were subject to 100000 steps of steepest descent minimization. The periodic boundary condition was applied in the simulation. The cut-off radius for the Lennard Jones and electrostatic interactions were set to be 10 Å. For accurate evaluation of the long range Coulomb interactions, Particle Mesh Ewald (PME)35 summation method is used for electrostatic interactions beyond the cut-off. The system was heated to 300 K in 100 ps by using a v-rescale36 thermostat for the canonical ensemble (NVT) simulations. During this process, the Linear Constraint Solver (LINCS) algorithm37 was used to constrain all the bond lengths. The isothermal isobaric ensemble (NPT) was used in the subsequent simulations, with the pressure set to 1 bar in 100 ps, controlled using a v-rescale thermostat and Parrinello-Rahman barostat38. Thereafter, the systems were simulated for 100 ns and 1 us at different cases separately. Three repeated MD simulations with different initial velocities were also performed. Throughout the simulations a time step of 2.0 fs was used. The trajectory of each simulation was recorded every 100000 time steps (200.0 ps), with the last 50 ns of the trajectory chosen for mean square displacement (MSD) calculation to obtain the diffusion coefficient of the catalyst.

ACS Paragon Plus Environment

6

Page 7 of 23 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 Catalysis

RESULTS AND DISCUSSIONS Formation of the prereactive dimer from two [RuV=O(bda)L2]+ monomers. We recently concluded that the formation of the preractive dimer from the free solvated monomers must be key to understanding the difference in reactivity between complexes with different axial ligands, e.g. 1 and 2.21 To understand this step in more detail we performed PMF simulations using umbrella sampling approach to estimate the binding free energy of encountering reaction between two mono-catalysts. Our initial structure was the DFT optimized prereactive dimers (Figure 3a, 3c), where the O-O distance is 3.5 Å and the Ru-Ru distance is 5.2 Å. The PMF simulations for the Ru-Ru distance were preceded by equilibrations for 100 ps under an NVT ensemble. The two Ru complexes were pulled apart for total 10 Å from the initial structure over 200 ps, using a spring constant of 10 kJ mol-1 Å-2 and a pull rate of 0.05 Å ps-1. From these trajectories, snapshots were taken to generate the starting configurations for the umbrella sampling windows. We then sampled the Ru-Ru distances using roughly a 0.5 Å spacing. Such spacing allowed for increasing detail at smaller Ru-Ru distance, and resulted in 21 windows. In each window, 10 ns of MD was performed for a total simulation time of 210 ns utilized for umbrella sampling. Analysis of results was performed with the weighted histogram analysis method (WHAM)39. These PMF simulations resulted in smooth dissociation curves (Figure S9) and the mean values from 5 repeated simulations were 3.2 kcal mol-1 for complex 1 and 6.3 kcal mol-1 for complex 2 (Table 1). The much larger binding energy for formation of the prereactive dimer for complex 2 should favorably impacts its tendency to form the dimeric structure, and thereby, possibly enhance the rate. When we compare the intermediate structures of two complexes along the PMF energy profile, a striking difference was found in the interaction between the molecules. Complex 2 appear to interact between the two isoq groups even at longer

ACS Paragon Plus Environment

7

ACS Catalysis 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

Page 8 of 23

distances. At the longer distances the isoq-isoq interaction appears to be more T-shaped than πstacking, as it is in the prereactive dimer (Figure 3b). In complex 1, however, instead of forming interactions between the axial picoline (pic) ligand, one pic form noncovalent interaction with the bda ligands of another [RuVO(bda)L2]+ species, mainly CH···O electrostatic interactions (Figure 3c). The interaction completely disappears when the distance is increased to around 9.5 Å.

a

c

b

d

Figure 3. Structure of PMF simulation. a, Prereactive dimers of complex 2 optimized with the DFT. b, Snapshots of complex 2 from PMF simulation, isoq-isoq groups tend to form T-shaped arrangement with a distance around 3.9 Å. c, Prereactive dimer structure of complex 1 optimized from DFT. d, Snapshots of complex 1 from PMF simulation. All PMF simulations are performed with TIP3P water, which have been removed in the figures for clarity.

ACS Paragon Plus Environment

8

Page 9 of 23 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 Catalysis

The binding free energy includes both the different intermolecular forces between the two catalysts, as well as the catalyst-solvent and solvent-solvent interactions. We found evidence for a hydrophobic nature of the oxo in RuV=O in our previous report, while the RuIV-OH was found to interact more strongly with water and seemed to be more hydrophilic. The oxo of complex 1 is also hydrophobic from hydrogen bonds analysis (Figure S6) and radial distribution functions (RDF) analysis (Figure 4). We therefore performed PMF simulations for the RuIV-OH complex to quantify the additional driving force the change from hydrophilic hydroxide to hydrophobic oxo induces.20 The difference in binding free energy of RuIV-OH is merely 1 kcal mol-1 smaller than that of RuV=O complex in 5 repeated simulations. This shows that it is a quite fine balance between hydrophilic and hydrophobic nature. From the intermediate structures of RuV=O and RuIV-OH complexes, we find different behavior at the initial configurations. In the RuIV-OH complex, one of the RuIV-OH units turns around from pointing to the other RuIV-OH, to pointing towards the water phase. This interaction with water indicates the prereactive dimer is not a strongly favored configuration (Figure 5a, 5b). In RuV=O complex on the other hand, prereactive dimer structure stays stable longer and the two RuV=O keep pointing towards each other until the distance is increased to around 8 Å (Figure 3c, 5c).

ACS Paragon Plus Environment

9

ACS Catalysis 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

Page 10 of 23

Figure 4. RDF analysis for O7 of complex 1 in the water phase with water. A first obvious solvation shell at 3.1-3.6 Å which is approximately the sum of the vdW radii, a result which indicates that the oxo is hydrophobic.

a

b

c

Figure 5. Structure from the PMF simulations. a, Prereactive dimers of RuIV-OH complex optimized with the DFT, the O-O distance is 3.5 Å. b, Snapshots of complex RuIV-OH from PMF simulation. c, Snapshots of complex 1 from PMF simulation. All PMF simulations are performed with TIP3P water, which have been removed in the figures for clarity. Diffusion of [RuV=O(bda)L2]+ in Solvent. The collisions rate of the reactants monomer is one factor that could influence rate constant and the full reaction of the O-O bond formation step. To study the diffusion we ran MD simulations in a periodic box filled with TIP3P water molecules and a chloride ion to balance the charge. From the MD trajectories of the complexes we calculated the diffusion coefficient in solvent by the mean square displacement calculation that is a part of the Gromacs 5.1.2 suite.32 To gain stable diffusion rates, we used long MD simulation time of 1 µs. The calculated diffusions are 5.1 × 10-6 cm2 s-1 of complex 1 and 4.5 × 10-6 cm2 s-1 of complex 2, which are in the same range as experimental value by Meyer and co-workers (1.6 × 10-6 cm2 s-1 in water solution).40 We would like to note that this diffusion is much faster compared to the diffusion of the pyrene substituted [RuVO(bda)L2]+ catalysts at the carbon

ACS Paragon Plus Environment

10

Page 11 of 23 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 Catalysis

nanotube (CNT) water interface.24 The complexes diffuse more freely in solvent than on the CNT surface, where the large contact area by pyrene-CNT π-stacking interaction could decrease the diffusion rate. This faster diffusion of complexes in solution is one contributing factor to the higher rate constant in water than catalyst at the CNT-water interface. The values of the two catalysts 1 and 2 are very similar, indicating that the diffusion rate is not causing the higher rate constant of complex 2. O-O Bond Formation in Solvent. Complex 2 was part of our previous report20, and the approach for 1 was identical to use EVB to simulate O-O bond coupling reaction both in the vacuum and aqueous phases. The initial structure was the same one from DFT optimized prereactive dimer. Before EVB simulation, equilibrations with distance restraint of 5 kcal/mol between the oxyl radicals, and the harmonic restraint of 0.2 kcal/mol on each atom were preceded for 1 ns. The EVB simulations then were performed in 51 EVB mapping frames of 10 ps in length each with the same restraints used in the equilibrations, using a 1 fs time step. In these frames, the two Ru complexes were gradually approached to form the O-O bond with a distance of 1.36 Å. The 10 simulations were started from equilibrations starting with ten different random seeds, which give rise to ten independent free energy profiles. All EVB calculations were performed using the standard EVB free energy perturbation procedure as implemented into the Q software package41 with the force fields for the reactant and products which were based on the OPLSAA force field.20 The EVB simulation mapping procedure requires the activation free energy and reaction energy as input data in order to fit the EVB parameters. This can be fitted either to high-level quantum chemical calculations, or experimental data. Due to the lack of relevant experimental data for these systems, we used reference values from DFT to parametrize EVB simulations in vacuum. We then used these parameters to simulate the reaction with EVB

ACS Paragon Plus Environment

11

ACS Catalysis 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

Page 12 of 23

in a 20 Å radius droplet of TIP3P water molecules centered on the centroid of the reacting atoms (Figure 6 shows the reaction coordinate of complex 1). Reactions in vacuum and in water proceeded smoothly from reactants to products with complex 1. In gas phase, complex 2 has a higher activation energy than complex 1 (by 3.5 kcal mol-1), which is likely due to higher favorable interaction energy in the prereactive dimer that need to be distorted to reach the transition state.21 In water phase, on the other hand, the activation free energy of complex 1 with pic ligands increased slightly compared to that in vacuum phase, while it decrease slightly for 2 with the isoq ligands. The different solvent effect on the complexes could be explained by the different π-stacking interaction. We found in our previous report that the π-stacking was largely solvent induced, meaning that the solvent can favor structures where the π-stacking tendency increases from reactant to product.20 This is exactly what we observe for 2, but much less so for 1. In solvent we calculate activation free energies of complex 1 and 2 for the O-O radical coupling step to be very close, 7.9 kcal mol-1 and 7.1 kcal mol-1, respectively (Table 1). This indicates that the energy barrier of the intrinsic O-O coupling step does not, in itself, contribute to the much higher rate constant of complex 2.

ACS Paragon Plus Environment

12

Page 13 of 23 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 Catalysis

Figure 6. The EVB energy profiles of complex 1 in the vacuum phase and solvent phase, energy profiles of complex 2 was in previous report20. The profile describes the relative energy of two complexes species as a function of the reaction coordinate. The O-O distance at transition state is about 1.89 Å in the gas phase and 1.92 Å in the water phase. Table 1. Activation free energy for forward (∆∆G) and reverse (∆∆G*) reaction of complexes. All values are in kcal mol-1.

∆∆G[a]

∆∆G*[a]

∆∆G[b]

∆∆G*[b]

gas

gas

water

water

1

4.4

21.4

7.9 (6.9)

20.2 (20.4)

3.2

2

7.9

24.8

7.1 (6.4)

28.9 (24.6)

6.3

Complex

[a]

All values calculated with DFT in the gas phase.

[b]

∆Gdiss

All values are calculated with EVB an

explicit TIP3P water environment. DFT results with the PBF solvation model21 are provided in parenthesis for comparison. Estimation of Rate Constant Based on Simulated Data. In the rate-limiting step the reaction rate constant depends upon the collision rate of the reactants monomer, a proper orientation of the reactants that could form prereactive dimer easily and enough energy to overcome the activation barrier. When the catalysts diffuse in the solvent, there is a chance of collision that leads to chemical reaction. The hydrophobic oxo of [RuVO(bda)L2]+ adds a driving force for the collision of the catalysts to form the prereactive dimer with the proper geometry, where two RuV

ACS Paragon Plus Environment

13

ACS Catalysis 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

Page 14 of 23

oxos pointing towards each other, hence, we ignored the steric factor in this work. With the analysis of diffusion rate, O-O bond formation step from free solvated monomer, we could use same methodology presented in our recent report24 to estimate rate constant, except that for the collision rate, we use three dimensional collision theory instead of two dimensional as in that report. We then make the assumption that the collision rate divided by the concentration is equal to the prefactor in the Arrhenius equation (eq 1). To estimate the activation energy Ea in the Arrhenius equation, weuse the values from the PMF and EVB simulations. Using the ∆GOO‡ from the EVB simulation as the forward activation energy, and the ∆Gdiss from the dissociation of the two catalysts from the PMF simulation, the combined value is then Ea = ∆GOO‡ - ∆Gdiss. We calculate this to be 4.7 kcal mol-1 of complex 1 and 0.8 kcal mol-1 of complex 2. The collision rate Φ for 3-dimensional diffusion have been derived by Hardt42 (2), and the prefactor A is assumed to be Φ/C2:

kobs = Ae− Ea /( RT )

(1)

Φ = 8π NaC 2 D

(2)

where D is the diffusion coefficient. C is the concentration catalysts in the medium, N is Avogadro’s number and a is the encounter radius of the reaction partners, which was assumed to be 10 Å with these catalysts. The calculated collision frequencies of 3-dimensional diffusion and the activation energies Ea, were then inserted in equation (1) yields a second-order rate constants are 1.38 × 106 M-1s-1 of complex 1 and 8.82 × 108 M-1s-1 of complex 2 which agrees very well with experimental results, 7.83 × 105 M-1s-1 and 2.39 × 108 M-1s-1 for complex 1 and 2, respectively. From studying the diffusion coefficient in the pre-exponential, the activation energy (the value of energy barrier minus binding free energy) in the exponential, we could analyze

ACS Paragon Plus Environment

14

Page 15 of 23 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 Catalysis

effect of different parts on the rate constant. Both the diffusion coefficient and the O-O bond activation energy from the respective prereactive dimers are close to identical, and therefore, these two steps cannot explain the difference in the constants of complex 1 and 2. Instead the only step of the process where a significant difference is seen is in the binding free energy of two [RuVO(bda)L2]+ monomers when the prereactive dimer is formed, indicating that the prereactive dimer formation is the main contributor for the difference in reactivity of the different complexes CONCLUSIONS In summary, we have studied the key aspects of the O-O bond formation between two Ru(V)oxo catalysts in a realistic environment, using dynamic simulations. In this work, we moved beyond the simple bond forming step and studied the formation of prereactive dimer from two free [RuVO(bda)L2]+ catalysts. By calculating diffusion of individual catalysts and extracting interaction energies between two complexes using PMF simulations, and EVB for the bond formation, we extracted all the necessary data for a full description of the reaction, including collision frequency, tendency to form the prereactive dimer, and activation energy of O-O bond formation. Combined these three parameters made it possible to calculate the rate constant of the catalysts, in excellent agreement with experimental values, to less than one order of magnitude. The binding energy of the isoq catalyst 2 is 3.1 kcal mol-1 larger than the pic catalyst 1, which is the only significant differing parameter of higher rate constant of the catalyst 2. From the prereactive dimer the activation free energy is close to identical, and therefore, the increased πstacking interaction of 2 is mainly contributing to the binding free energy rather than stabilizing the transition state of the radical coupling step. By comparing the binding free energy of hydrophilic RuIV-OH and hydrophobic RuV=O, we estimate the hydrophobic driving force of RuV=O in this system, to 1 kcal mol-1. While this value is small, this driving force, together with

ACS Paragon Plus Environment

15

ACS Catalysis 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

Page 16 of 23

the hydrophobicity of the rest of the molecule assists the formation of the prereactive dimer – the key step of the reaction. This report illustrates how complex processes where solvent behavior affects the reaction can be simulated by a combination of quantum and classical methodologies. ASSOCIATED CONTENT Supporting Information. The Supporting information is available free of charge. Details of the computational study, force field parameterization, molecular dynamics simulations and empirical valence bond simulations, complexes structures in vacuum and water phase, stability analysis, H-bond analysis, radial distribution functions analysis and potential of mean force analysis and EVB free energy profiles (PDF). AUTHOR INFORMATION Corresponding Author * [email protected] Funding Sources This work was supported by Vetenskapsrådet and the China Scholar-ship Council (CSC). Notes The authors declare no competing financial interest. ACKNOWLEDGMENT This work was supported by Vetenskapsrådet and the China Scholarship Council (CSC). All calculations were performed on resources provided by the Swedish National Infrastructure for

ACS Paragon Plus Environment

16

Page 17 of 23 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 Catalysis

Computing (SNIC) at PDC Centre for High Performance Computing (PDC-HPC), Uppsala Multidisciplinary Center for Advanced Computational Science (UPPMAX) under the project number SNIC2017-1-339 and the National Supercomputing Center in Linköping, Sweden. REFERENCES (1) Francàs L.; Matheu R.; Pastor E.; Reynal A.; Berardi S.; Sala X.; Llobet A.; Durrant J. R. Kinetic Analysis of an Efficient Molecular Light-Driven Water Oxidation System. ACS Catal. 2017, 7, 5142-5150. (2) Chen Z.; Concepcion J. J.; Hu X.; Yang W.; Hoertz P. G.; Meyer T. J. Concerted O Atom– Proton Transfer in the O–O Bond Forming Step in Water Oxidation. Proc. Natl. Acad. Sci. USA 2015, 107, 7225-7229. (3) Gersten, S. W.; Samuels, G. J.; Meyer, T. J. Catalytic Oxidation of Water by an Oxo-Bridged Ruthenium Dimer. J. Am. Chem. Soc. 1982, 104, 4029-4230. (4) Kärkäs, M. D.; Verho, O.; Johnston, E. V.; Åkermark, B. Artificial Photosynthesis: Molecular Systems for Catalytic Water Oxidation. Chem. Rev. 2014, 114, 11863-12001. (5) Garrido-Barros P.; Gimbert-Suriñach C.; Matheu R.; Sala X.; Llobet A. How to Make an Efficient and Robust Molecular Catalyst for Water Oxidation. Chem. Soc. Rev. 2017, 46, 60886098. (6) Zong R.; Thummel R. P. A New Family of Ru Complexes for Water Oxidation. J. Am. Chem. Soc. 2005, 127, 12802-12803.

ACS Paragon Plus Environment

17

ACS Catalysis 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

Page 18 of 23

(7) Ellis, W. C.; McDaniel, N. D.; Bernhard, S.; Collins, T. J. Fast Water Oxidation Using Iron. J. Am. Chem. Soc. 2010, 132, 10990-10991. (8) Barnett, S. M.; Goldberg, K. I.; Mayer, J. M. A Soluble Copper-Bipyridine Water-Oxidation Electrocatalyst. Nature Chem. 2012, 4, 498-502. (9) Sartorel, A.; Miró, P.; Salvadori, E.; Romain, S.; Carraro, M.; Scorrano, G.; Valentin, M. D.; Llobet, A.; Bo, C.; Bonchio, M. Water Oxidation at a Tetraruthenate Core Stabilized by Polyoxometalate Ligands: Experimental and Computational Evidence To Trace the Competent Intermediates. J. Am. Chem. Soc. 2009, 131, 16051-16053. (10) Duan, L.; Wang, L.; Inge, A. K.; Fischer, A.; Zou, X.; Sun, L. Insights into Ru-Based Molecular Water Oxidation Catalysts: Electronic and Noncovalent-Interaction Effects on Their Catalytic Activities. Inorg. Chem. 2013, 52, 7844-7852. (11) Duan, L.; Bozoglian, F.; Mandal, S.; Stewart, B.; Privalov, T.; Llobet, A.; Sun, L. A Molecular Ruthenium Catalyst with Water-Oxidation Activity Comparable to That of Photosystem II. Nature Chem. 2012, 4, 418-423. (12) Duan, L.; Wang, L.; Li, F.; Li, F.; Sun, L. Highly Efficient Bioinspired Molecular Ru Water Oxidation Catalysts with Negatively Charged Backbone Ligands. Acc. Chem. Res. 2015, 48, 2084-2096. ( 1 3) Schulze M.; Kunz V.; Frischmann P, D.; Würthner F. A Supramolecular Ruthenium Macrocycle with High Catalytic Activity for Water Oxidation That Mechanistically Mimics Photosystem II. Nature Chem. 2016, 8, 576-583.

ACS Paragon Plus Environment

18

Page 19 of 23 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 Catalysis

(14) Meyer T. J.; Sheridan M. V.; Sherman B. D. Mechanisms of Molecular Water Oxidation in Solution and on Oxide Surfaces. Chem. Soc. Rev. 2017, 46, 6148-6169. (15) Shaffer D. W.; Xie Y.; Concepcion J. J. O–O Bond Formation in Ruthenium-Catalyzed Water Oxidation: Single-Site Nucleophilic Attack vs. O–O Radical Coupling. Chem. Soc. Rev. 2017, 46, 6170-6193. (16) Nyhlén, J.; Duan, L.; Åkermark, B.; Sun, L.; Privalov, T. Evolution of O2 in a SevenCoordinate RuIV Dimer Complex with a [HOHOH]− Bridge: A Computational Study. Angew. Chem. Int. Ed. 2010, 49, 1773-1777. (17) Hetterscheid, D. G. H.; Reek, J. N. H. Mononuclear Water Oxidation Catalysts. Angew. Chem. Int. Ed. 2012, 51, 9740-9747. (18) Duan, L.; Fischer, A.; Xu, Y.; Sun, L. Isolated Seven-Coordinate Ru(IV) Dimer Complex with [HOHOH]− Bridging Ligand as an Intermediate for Catalytic Water Oxidation. J. Am. Chem. Soc. 2009, 131, 10397-10399. (19) Lebedev, D.; Pineda-Galvan, Y.; Tokimaru, Y.; Fedorov, A.; Kaeffer, N.; Copéret, C.; Pushkar, Y. The Key RuV=O Intermediate of Site-Isolated Mononuclear Water Oxidation Catalyst Detected by in Situ X-ray Absorption Spectroscopy. J. Am. Chem. Soc. 2018, 140, 451458. (20) Zhan, S.; Mårtensson, D.; Purg, M.; Kamerlin, S. C. L.; Ahlquist, M. S. G. Capturing the Role of Explicit Solvent in the Dimerization of RuV(bda) Water Oxidation Catalysts. Angew. Chem. Int. Ed. 2017, 56, 6962-6965.

ACS Paragon Plus Environment

19

ACS Catalysis 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

Page 20 of 23

(21) Fan, T.; Zhan, S.; Ahlquist, M. S. G. Why Is There a Barrier in the Coupling of Two Radicals in the Water Oxidation Reaction? ACS Catal. 2016, 6, 8308-8312. (22) Wang, L.; Duan, L.; Wang, Y.; Ahlquist, M. S. G.; Sun, L. Highly Efficient and Robust Molecular Water Oxidation Catalysts Based on Ruthenium Complexes. Chem. Commun. 2014, 50, 12947-12950. (23) Richmond, C. J.; Matheu, R.; Poater, A.; Falivene, L.; Benet-Buchholz, J.; Sala, X.; Cavallo, L.; Llobet, A. Supramolecular Water Oxidation with Ru–bda-Based Catalysts. Chem. Eur. J. 2014, 20, 17282-17286. (24) Zhan, S.; Ahlquist, M. S. G. Dynamics and Reactions of Molecular Ru Catalysts at Carbon Nanotube–Water Interfaces. J. Am. Chem. Soc. 2018, 140, 7498-7503. (25) Warshel, A. Computer Modelling of Chemical Reactions in Enzymes and Solutions, WileyInterscience, New York, 1997. (26) Bochevarov, A. D.; Harder, E.; Hughes, T. F.; Greenwood, J. R.; Braden, D. A.; Philipp, D. M.; Rinaldo, D.; Halls, M. D.; Zhang, J.; Friesner, R. A. Jaguar: A High-Performance Quantum Chemistry Software Program with Strengths in Life and Materials Sciences. Int. J. Quantum Chem. 2013, 113, 2110-2121. (27) Zhao, Y.; Trular, D. G. The M06 Suite of Density Functionals for Main Group Thermochemistry, Thermochemical Kinetics, Noncovalent Interactions, Excited States, and Transition Elements: Two New Functionals and Systematic Testing of Four M06-class Functionals and 12 Other Functionals. Theor. Chem. Acc. 2006, 120, 215-241.

ACS Paragon Plus Environment

20

Page 21 of 23 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 Catalysis

(28) 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 HPu. J. Chem. Phys. 2010, 132, 154104-154109. (29) Goerigk, L.; Grimme, S. A Thorough Benchmark of Density Functional Methods for General Main Group Thermochemistry, Kinetics, and Noncovalent Interactions. Phys. Chem. Chem. Phys. 2011, 13, 6670-6688. (30) Hay, P. J.; Wadt, W. R. Ab Initio Effective Core Potentials for Molecular Calculations. Potentials for K to Au Including the Outermost Core Orbitals. J. Chem. Phys. 1985, 82, 299310. (31) Jorgensen, W.; Maxwell, D.; Tirado-Rives, J. Development and Testing of the OPLS AllAtom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 118, 11225-11236. (32) Berendsen, H. J. C.; Vanderspoel D.; Vandrunen R. GROMACS: A Message-passing Parallel Molecular Dynamics Implementation. Comput. Phys. Commun. 1995, 91, 43-56. (33) Abraham, M. J.; Murtola, T.; Schulz, R.; Pall, S.; Smith, J. C.; Hess, B.; Lindahl, E. GROMACS: High Performance Molecular Simulations Through Multi-level Parallelism from Laptops to Supercomputers. SoftwareX 2015, 1-2, 19-25. (34) 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, 926935.

ACS Paragon Plus Environment

21

ACS Catalysis 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

Page 22 of 23

(35) Essmann, U.; Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577-8593. (36) Bussi, G.; Donadio, D.; Parrinello, M. Canonical Sampling Through Velocity Rescaling. J. Chem. Phys. 2007, 126, 014101. (37) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS: A Linear Constraint Solver for Molecular Simulations. J. Comput. Chem. 1997, 18, 1463-1472. (38) Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52, 7182-7190. (39) Kumar, S.; Rosenberg, J. M.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A. THE Weighted Histogram Analysis Method for Free-Energy Calculations on Biomolecules. I. The Method. J. Comput. Chem. 1992, 13, 1011-1021. (40) Song, N.; Concepcion, J. J.; Binstead, R. A.; Rudd, J. A.; Vannucci, A. K.; Dares, C. J.; Coggins, M. K.; Meyer, T. J. Base-Enhanced Catalytic Water Oxidation by a Carboxylate– Bipyridine Ru(II) Complex. Proc. Natl. Acad. Sci. USA 2015, 112, 4935-4940. (41) Marelius, J.; Kolmodin, K.; Feierberg, I.; Åqvist, J. Q: a Molecular Dynamics Program for Free Energy Calculations and Empirical Valence Bond Simulations in Biomolecular Systems. J. Mol. Graph. Model 1998, 16, 213-225. (42) Hardt, S. L. Rates of Diffusion Controlled Reactions in One, Two and Three Dimensions. Biophys. Chem. 1979, 10, 239-243.

ACS Paragon Plus Environment

22

Page 23 of 23 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 Catalysis

TOC

ACS Paragon Plus Environment

23