Exploring Solvation Effects in Ligand-Exchange Reactions via Static

Article pubs.acs.org/JCTC

Exploring Solvation Effects in Ligand-Exchange Reactions via Static and Dynamic Methods Florian H. Hodel,† Peter Deglmann,‡ and Sandra Luber*,† †

Department of Chemistry, University of Zurich, Winterthurerstrasse 190, CH-8057 Zurich, Switzerland BASF SE, Carl-Bosch-Straße 38, 67056 Ludwigshafen am Rhein, Germany

‡

S Supporting Information *

ABSTRACT: We investigate ligand-exchange reactions of a biomimetic Co(II)-based heterocubane complex in aqueous solution by means of various approaches for consideration of solvent effects. Static calculations based on geometry optimizations carried out in vacuum, with solvent continuum models, or with several explicit solvent molecules have been carried out as well as density functional theory (DFT)-based molecular dynamics simulations. In addition, reaction pathways and barriers have been elucidated via nudged elastic band calculations and metadynamics. The results show that static approaches with approximate consideration of the solvent environment lead to reaction energies, which may change drastically depending on the method employed. A more sophisticated approach is DFTmolecular dynamics at ambient conditions with full solvation, i.e. enough solvent molecules to retain bulk water properties far from the solute, which, however, comes with a much higher computational cost. The investigated example of the exchange of an acetate ligand by a hydroxide demonstrates that entropic contributions can be vital and consideration of electronic energies alone may be a rather rough approximation.

■

(1) (see Figure 1) and [CoII3Er(hmp)4(OAc)5H2O].2 The goal of this study had been to determine the most probable structures of the ground states of the catalytic cycles paving the way for further investigations regarding the energetics of catalysis itself.3,4 Calculating differences between average electronic energies from sampling with Born−Oppenheimer molecular dynamics (MD) with Kohn−Sham density functional theory (DFT)5,6 as the electronic structure method, we had looked at the exchange of an acetate ligand with hydroxide (see Figure 1 for their naming): one of the monodentate acetate ligands of 1 has been replaced by hydroxide (resulting in 1-OH). In addition, we had detached some of the bridging ligands from only one metal center, breaking the colored bonds labeled “a” and “b” in Figure 1, and had placed a water molecule at the now vacant coordination site (resulting in 1μ1(a) and 1-μ1(b)). Among other, more subtle effects, we had found far higher ligand mobility for the heterocubane [CoII3Er(hmp)4(OAc)5H2O] than for 1.2 One issue that had become apparent during this study had been the strong dependence of the ligand exchange reaction energies on solvation. In this present work, we therefore set out to further elucidate on that matter. To this end, we explored different levels of solvation ranging from hundreds of explicit DFT-water

INTRODUCTION Proper description of environmental effects is a key point for reliable modeling of chemical processes. In fact, many chemical reactions take place in the condensed phase. This is also true for catalytic reactions. Homogeneous catalysis, for instance, typically occurs in the liquid phase, whereas heterogeneous catalysis often involves a liquid−solid interface. Photocatalytic water splitting is a prominent example and of paramount importance for environmentally friendly production of molecular hydrogen needed for sustainable energy storage and conversion. With the water oxidation half-reaction constituting the bottleneck of the process,1 the development of efficient and cheap water oxidation catalysts (WOCs) is highly desirable. It not only takes place in the realm of classical chemistry but also is driven forward with ab initio computational studies. In such calculations, static approaches are mostly employed where solvent effects are standardly included via computationally cheap, approximate solvent continuum models and/or consideration of several explicit solvent molecules. Investigating a WOC mechanism in solution requires knowledge of the structures occurring during the catalytic cycle, which are often experimentally unknown. Calculations can be of great help but need a sound description of solvent effects to provide meaningful results. In a previous study, we had investigated the ligand mobility of two Co(II)-based biomimetic WOCs, [CoII4(hmp)4(μOAc)2(μ2-OAc)2(H2O)2] (hmp = 2-(hydroxymethyl)pyridine) © XXXX American Chemical Society

Received: March 1, 2017

A

DOI: 10.1021/acs.jctc.7b00214 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

dielectric continuum depends on the parametrization that was used to ensure that solvation free energies reproduce experimental values. These energies therefore contain certain not directly quantifiable thermal corrections and do not exactly represent the electronic energy of the solute in solution.22 Especially for the calculation of reaction energies, people have thus attempted to model their systems with at least those water molecules they deemed crucial for the reaction mechanism explicitly represented.23 Computational cost is unfortunately not the only drawback of explicit solvent models. One also needs to be very careful when calculating energy differences or reaction energies to have both systems in the “same” local energy minimum with respect to the solvent coordinates, otherwise a big contribution to the energy difference will stem simply from the different solvation shell structures. In general, due to the limited sizes of the systems, large fluctuations in electronic energy differences can occur, even for similar solute configurations or reactions with small structural changes. Implicit solvation models can be thought of as providing an average or “configurationally sampled” solvent effect.24 The configuration of the solute is however still a static one and is not necessarily representative of the statistical ensemble at the desired temperature.25 In this lies also another important difference between explicit and implicit solvation models, which we broached above but want to stress here again: Dielectric continuum solvent models are parametrized to give accurate free energies of solvation (compared to experimental solvation energies) at a certain temperature,22 whereas the explicit solvent approach delivers electronic/free energies within the limits of DFT. DFT-based MD calculations, albeit at a far higher computational cost and relatively rarely used, have some distinct advantages over “static” ones.26,27 Given the size of the box is large enough and the duration of equilibration and production runs long enough, and assuming that the generated trajectory is ergodic, sampling at a finite temperature, containing an average over many different configurations, should yield a more realistic picture than calculating electronic energies of frozen systems. Of course, in practice, considering the obtained time averages of certain properties equal to the corresponding ensemble averages is only an approximation. Additionally, through the use of periodic boundary conditions, the artificial phase boundaries vanish, and at a certain distance from the solutes, bulk properties of water are replicated in large enough boxes. Furthermore, the need to “intelligently guess” the structures of important conformations introduces some degree of bias in static calculations, which is absent in MD. Methods. All calculations were carried out using either the CP2K28 or the Turbomole29 program package (version 6.5). With the former, the electronic energies (and forces) were calculated using spin unrestricted DFT with the mixed Gaussian and plane waves method as implemented in the QUICKSTEP program.30 We chose the BP86 exchangecorrelation (XC) density functional,31,32 which has been reported to deliver reliable structures of transition metal complexes in static as well as dynamic calculations (including solvation).33−39 Additionally, we employed the D3 dispersion correction by Grimme et al.,40 which has been shown to provide an improved description of transition metal complexes,41 as well as to counteract the overstructuring of water by the GGA XC functional (see the SI, section Bulk Properties of Water).42−45 The energy cutoff for the auxiliary plane wave

Figure 1. Molecular structure of 1 with one of the two monodentate acetate ligands visible in the lower left corner. The two bidentate, bridging acetate ligands are located on the “top” and “bottom” of the cubane cage, and the bonds broken in order to obtain 1-μ1(a) and 1μ1(b) are labeled a and b, respectively.

molecules in a periodic DFT-based MD framework to static approaches in vacuo, with explicit solvation shells and/or solvent continuum models as well as quantum mechanics/ molecular mechanics (QM/MM) calculations. In addition, barriers and free energies are investigated using the nudged elastic band method (NEB)7 and metadynamics,8 respectively. Since computational investigations of ligand-exchange reactions, especially for systems containing metals, are rare in literature,9−11 this report may serve as an important guideline for future studies.

■

THEORY Solvation Models. Energy differences and other properties of solvated molecules are often obtained from geometry optimizations in one of two ways:12 by investigating either molecules in an implicit solvent environment or systems with not only the solutes but also a number of solvent molecules explicitly represented (perhaps surrounded again by a dielectric).13,14 What that number might be depends on the specific system and the effects to be captured and can range from a few molecules that might be necessary to stabilize a certain ligand (microsolvation)15,16 to a full first17 and even second solvation shell18 describing short-range effects such as hydrogen-bonding in more detail. When applying an implicit solvent model, the solute is placed inside a cavity in a polarizable medium, and electrostatic interactions are calculated. Van der Waals interactions can be incorporated via parametrization.19 The treatment of the interactions between the cavity and the dielectric continuum is the main distinction between the available methods.20 This continuum representation of the solvent is computationally cheap and is therefore more commonly found in literature than the explicit solvent model. Based mainly on the only small structural differences between conformations obtained with implicit and explicit solvation, it has been suggested that modeling of actual water molecules is not essential for the calculation of certain properties.21 It has however been shown that reaction energies can vary greatly depending on the continuum model employed as well as the shape and size of the cavity.20 The electronic energy of the solute within the B


Article


pitfalls of ab initio MD can be found in the SI, section Bulk Properties of Water. The starting configurations for the geometry optimizations of 1 including a shell of 68 explicit water molecules were extracted from equilibrated structures by deleting water molecules at a distance greater than 9 Å from the heterocubane’s center of mass. The initial coordinates of all other systems were derived from those two to ensure the structural similarity of the solvation shells. In contrast to the DFT-MD simulations, these systems did not contain solvated OH− or Ac− ions, the energies of which were instead calculated separately from geometry optimizations containing 25 water molecules. Electronic energy differences of ligand exchanges were obtained as e.g. ΔE(1−1OH) = E(1) + E(OH−) − E(1-OH) − E(Ac−). The initial geometries for the structural optimizations with COSMO or in vacuo were obtained by deleting all water molecules of the respective explicit solvation shell systems. Furthermore, for structures optimized within the solvation description by the COSMO model, we obtained free energy differences using a protocol involving the solvation model COSMO-RS,58 which is, compared to COSMO, a more sophisticated description based on statistical thermodynamics and better takes into account hydrogen-bonding by explicitly optimized terms.59 Although each species in solution is described in this model only by its so-called σ-profile, which represents a histogram of surface charge (screening) density, i.e. the possibility of specific local interactions is not taken into account, the accuracy of COSMO-RS Gibbs free energies of solvation of smaller molecules in solution was found to be quite high (0.32 kcal/mol, i.e. 1.34 kJ/mol in a rather old parametrization60). That COSMO-RS solvation treatment can yield also reasonable differences in Gibbs free energies for aqueous systems involving species of strongly differing ionicity has been demonstrated for e.g. aqueous acid−base equilibria,61 so that also for the equilibria studied in this work a reliable treatment of solvation effects is expected. The primary outcome of COSMO-RS are chemical potentials or vapor pressures, which directly translate into the aforementioned Gibbs free energies of solvation. For the computation of enthalpies and entropies, the Gibbs−Helmholtz equation can be employed,62 which involves the computation of derivatives of the Gibbs free energy with respect to temperature. Although it has been shown that for systems involving ionic species with different coordination power, enthalpies and entropies are of lower accuracy than Gibbs free energies;61 at least an idea of the dimension of entropic effects can be obtained in this way. The COSMO-RS calculations of Gibbs free energies of solvation (and enthalpies of solvation by forming the temperature derivative according to the Gibbs−Helmholtz equation62) were performed as recommended based on singlepoint calculations at the BP86/def-TZVP and BP86+COSMO/ def-TZVP level. For certain solvation shell systems, we included additionally a dielectric continuum using a method implemented in CP2K. It has to be noted, however, that it differs significantly from COSMO in several aspects: For example, whereas the cavity in COSMO is described by the solvent accessible surface which is constructed from atom-centered spheres with radii proportional to the respective van der Waals radii and also depends on the assumed hard-sphere radius of the solvent molecules, the selfconsistent reaction field (SCRF) model in CP2K simply assumes a spherical cavity with a user-defined radius. We used a

expansion of the charge density was set to 400 Ry, and for the evaluation of the XC functional and its derivative (which uses the same grid), a nearest-neighbor smoothing approach (NN50 smoothing) was employed,30 a setup that has already successfully been applied to systems of 1 and other heterocubanes,2−4 as well as in an ab initio MD study of liquid water with a GGA XC functional.45 The core electrons were replaced in all calculations with the norm conserving Goedecker-Teter-Hutter (GTH) pseudopotentials.46 For cobalt, we used DZVP-MOLOPT-SR-GTH, and for all other atom types we used DZVP-MOLOPT-GTH basis sets.47 The Co(II) centers were assumed to be in the quartet spin state, coupling ferromagnetically.48,49 The calculations with Turbomole either used the conductorlike screening model COSMO50 or were carried out in vacuo. We employed the BP86 XC functional and Grimme’s D3 dispersion correction. def2-TZVP basis sets51,52 were used as well as the resolution-of-the-identity density-fitting technique53 with corresponding auxiliary basis sets.54 DFT-MD sampling, the results of which we already reported in a previous study,2 was performed with CP2K. We refer the reader to the above-mentioned work for a detailed description of the computational setup. In short, we first ran 100 ps of classical MD in a periodic cubic box using the CHARMM force field55 with the TIP3P water model56 in the NPT ensemble with positional constraints on all solute atoms. Next, DFT-MD trajectories were generated with a time step of 0.5 fs in the NVT ensemble without any constraints. The temperature was kept around 300 K using the canonical sampling through velocity rescaling thermostat.57 After discarding the first 5 to 8 ps as equilibration, we averaged the electronic energy calculated at every time step over another 2 to 4 ps. With this setup, the 95% confidence intervals for the averages of the electronic energies could be kept smaller than ±8 kcal/mol. Due to limited computational resources, longer production runs had not been possible. Unless otherwise noted in the Results and Discussion section, each system is solvated with 328 water molecules and contains either a heterocubane molecule with one monodentate acetate ligand coordinated to it and OH− at some distance in solution or a heterocubane molecule with the aforementioned acetate in solution and the hydroxide ligand being attached to the heterocubane, instead. For the systems with partly detached ligands, additionally, one oxygen atom of a chelating acetate ligand was detached from cobalt by pivoting the whole acetate molecule around the other oxygen atom, which thus remained at the same position. Next, the now vacant coordination site on Co was occupied by the water molecule closest to it. The electronic energy differences were then calculated as the differences between average electronic energies over MD runs of two systems. Finally, we generated DFT-MD trajectories with OH−, Ac−, 1-OH, and 1, each solvated in its own periodic box containing 191 water molecules in the case of the two anions and again 328 for the two cubanes, thus enabling us to obtain the electronic reaction energy by simply subtracting the sums of the average electronic energies of a ligand and the corresponding complex. The equilibration and sampling of these systems was carried out analogously to the other MD simulations described above. We investigated some (electronic) structure properties to ensure sufficient solvation and the retention of bulk properties at large distances from the solutes. The results of these investigations and a short overview of some related possible C


Article


in Table 1). In agreement with our findings from the aforementioned simulations where we placed each cubane and corresponding ion in the same box, the method yielded a positive electronic energy difference, deeming 1 + OH− thermodynamically more favorable than 1-OH + Ac−. This pattern changes drastically when solvation is reduced, and instead of sampling, static calculations are performed. As above, OH−, Ac−, 1-OH, and 1 were each solvated separately; this time, however, only with 68 water molecules for the cubane complexes and 25 for the ligand molecules. We then carried out geometry optimizations and calculated the electronic energy differences between the local minimum energy structures, which resulted in an energetically almost neutral picture (ΔE = −0.7 kcal/mol, see Table 1). Repeating the calculations without any explicitly represented solvent, i.e. using COSMO and in vacuo, further lowered the energy of 1−OH + Ac−, with respect to the one of 1 + OH−, rendering the complex with hydroxide ligand now clearly thermodynamically more favorable. The only significant difference in the conformation of the solutes between explicit and implicit solvation is the orientation of the hydroxide ligand of 1-OH. There is a clear trend visible going from full solvation to explicit solvation shell, to implicit solvation, and finally to in vacuo calculations. The reaction f ree energy obtained with COSMO-RS complements this picture by a value of −6.9 kcal/ mol. However, the entropic contribution was found to account for altogether almost 20 kcal/mol, as the enthalpy difference exhibits a positive value of 12.6 kcal/mol, which agrees better with the results obtained from DFT-MD. It is worth mentioning that the enthalpy difference is approximately given by the electronic energy difference within an NVT ensemble for the DFT-MD runs (neglecting contributions from pressure fluctuations). Enthalpy differences are, however, in general more difficult to converge than free energies.66,67 Setting aside the perhaps insufficient treatment of solvation effects, the COSMO-RS result nevertheless suggests that entropic corrections may have a decisive contribution to reaction free energies, which are not covered by considering only electronic energies as given in Table 1. For the partial dissociation of chelating acetate ligands and the subsequent attachment of water, we found that for 1 there was at least for three out of four systems again a similar trend visible when going from DFT-MD simulations to explicit solvation shell to COSMO (see Table 2). The four systems we investigated were created by breaking bonds a or b (see Figures 1 and 2) of 1 or 1-OH and attaching a water molecule to the

radius of 11 Å around the center of mass of the system under study. For selected single point energy calculations using explicit solvent, as well as for all systems with COSMO or in vacuo, we had compared the electronic energies to the ones obtained employing the B3LYP XC functional63,64 instead of BP86 with all other settings unchanged.2 For the sake of consistency, we focus on the BP86 results in this manuscript. Furthermore, except for the metadynamics studies and the free energies obtained with COSMO-RS, we only investigated electronic energy differences based on the assumption that the electronic energy is the by far dominating contribution to the free energy,2 which is often made in static calculations. All reaction path studies were performed using settings similar to those chosen for the geometry optimizations. A short summary of the theoretical background of the metadynamics can be found in the SI, section Theory. The NEB calculations consisted of 16 replicas along the pathways optimized by a climbing image NEB every fifth optimization step, the others utilizing the improved tangent method.65 For the metadynamics calculations, we used only coordination numbers as collective variables (CVs). Gaussian potentials were added every 40 MD steps. Their heights and widths for every metadynamics simulation performed in this study can be found in the SI, section NEB and Metadynamics Parameters.

■

RESULTS AND DISCUSSION 1. Reaction Energies. For the dissociation of one of the monodentate acetate ligands from 1, i.e. [CoII 4 (hmp)4 (μ − OAc)2 (μ2 − OAc)2 (H 2O)2 ] + OH− → [CoII 4 (hmp)4 (μ − OAc)(μ2 − OAc)2 (H 2O)2 OH] + Ac−

we compared electronic energy differences obtained with different solvation models (see Table 1). Table 1. Electronic Energy Differences for the Exchange of a Monodentate Acetate of 1 with OH− sampling

ΔE [kcal/mol]

geometry optimization

ΔE [kcal/mol]

DFT-MD 1 DFT-MD 2 fragments

34.3 47.8 20.5

solvation shell COSMO vacuum

−0.7 −5.8 −25.2

We carried out two DFT-MD runs with the acetate ligand in solution. The initial part of the first trajectory was generated by heating a geometry-optimized system; the initial structure of the second one was obtained from an earlier DFT-MD run. In addition, another DFT-MD run was performed with the acetate bound to 1. The average electronic energy difference between the latter and the two DFT-MD runs with acetate in solution are labeled DFT-MD 1 and DFT-MD 2 in Table 1. The difference between the DFT-MD electronic energies is most likely rooted in the fact that the two simulations did not end up sampling exactly the same potential energy basin. This issue also manifests itself in the abstraction of H+ from water by the coordinated hydroxide ligand leading to a [CoII4(hmp)4(μOAc)(μ2-OAc)2(H2O)3] species followed by hydrogen bonding the newly created hydroxide in solution for the remaining duration of DFT-MD 1. This contrasts with the second DFTMD run, where OH− remained inertly attached. Next, we generated DFT-MD trajectories with OH−, Ac−, 1OH, and 1, each solvated in its own periodic box (“fragments”

Table 2. Electronic Energy (for COSMO-RS: Free Energy) Differences in kcal/mol for the Partial Removal of Bidentate Ligands from 1 1-μ1(a) vs 1 DFT-MD sampling static calculations: solv. shell 1 solv. shell 2 solv. shell 3 COSMO COSMO-RS free energy (in brackets: enthalpy) D

1-μ1(b) vs 1

1-OHμ1(a) vs 1-OH

1-OHμ1(b) vs 1-OH

30.1

17.5

−2.0

1.1

10.9 11.6 5.4 −6.9 7.3 (2.3)

12.4 16.8 −2.3 −8.8 10.2 (7.8)

−4.5 −3.2 8.5 −14.7 4.4 (5.9)

19.8 −1.1 10.5 −16.4 4.3 (8.0)


Article


underlying principles of this continuum model differ from those of COSMO. Furthermore, we performed geometry optimizations using a QM/MM approach with 68 water molecules around the cubane represented quantum mechanically and 260 additional water molecules calculated with the same setup we used for the classical equilibration. Finally, we employed COSMO in addition to a first solvation shell of 45 water molecules. Whereas the results for systems with explicit solvent molecules (except for the one with the simple SCRF solvent model) (and from the COSMO-RS free energy) suggest that this reaction is thermodynamically not favorable, COSMO predicts the opposite. In order to investigate the changes in electronic energy difference with the number of explicitly represented solvent molecules, several water molecules (for exact numbers, see Table 4) were randomly selected and removed starting from

Figure 2. Partial dissociation of bridging ligands from 1-OH. Either bond a or bond b (see Figure 1) is broken resulting in 1-OH-μ1(a) (green) or 1-OH-μ1(b) (orange). For the DFT-MD runs, one monodentate acetate ligand of 1 was displaced by hydroxide into solution (background, top).

Table 4. Electronic Energy Differences in kcal/mol between 1-μ1(a) and 1 with Decreasing Number of Explicit Water Molecules Included in the System solvent fixed

now vacant coordination site. The structures of the solutes were again similar for explicit and implicit solvation. Only the orientation of the hydroxide and the partially detached bridging acetate ligands differs slightly. Compared to the COSMO results, the COSMO-RS values, which, however, are free energy differences instead of electronic energy differences, show an improved agreement with the DFT-MD results. In all investigated cases shown in Table 2, neglecting the entropic contributions (between 1.5 and 5 kcal/mol in the COSMO-RS results) to the free energy differences leads to less agreement with the average electronic energy differences obtained from DFT-MD but does not result in a change of sign of the corresponding values. In order to be able to assess the dependency of the electronic energy differences on the choice of initial atom positions, we repeated the first solvation shell calculations for another two randomly chosen different starting configurations also extracted from DFT-MD equilibration runs (labeled “solv. shell 2” and “solv. shell 3” in Table 2). These systems differ not only in the positions of the water molecules but also in the conformation and orientation of the hydroxide and bidentate acetate ligands (whether partially detached or not). The standard deviations of the electronic energy differences over the three starting configurations of the four systems are between 3.4 and 10.4 kcal/mol demonstrating that the energy differences do in fact strongly depend on the chosen initial configurations. In Table 2, the values calculated with COSMO were significantly lower than all others. For the example “1-μ1(a) vs 1”, we therefore added the simple, rather approximate SCRF continuum model implemented in CP2K to the first solvation shell (see Table 3). As mentioned in the Methods section, the

68 67 60 55 45 35 25 15

no constraints

ΔE

ΔΔE per water molecule

ΔE


10.9 7.5 −1.0 0.2 −23.0 −19.3 −7.7

0.0 −0.4 0.1 0.6 −3.3 −1.5 −0.3

10.9 10.4 −1.7 −5.4 9.9 −4.4 −5.0 9.7

−0.5 −1.3 −0.9 1.0 1.9 1.4 1.7

the optimized structure of 1-μ1(a) including 68 H2O. We repeated this with the corresponding water molecules surrounding 1. The remaining ones were fixed in space, the geometries of the two cubane-structures were optimized (columns 2 and 3 of Table 4), and finally, the whole systems were optimized with no constraints (columns 4 and 5 of Table 4). This procedure was repeated until the solvation shells consisted of only 15 water molecules. Also reported in Table 4 are the changes in electronic energy difference between 1-μ1(a) and 1 per deleted water molecule (ΔΔE) when releasing the constraint (column 5) or deleting solvent molecules (column 3). Close to the site of the partial detaching of the bidentate acetate ligand, the water shells of the two systems differ the most. If a water molecule in this region and its counterpart in the other system are deleted, the resulting change in electronic energy difference will be bigger than the one for the deletion of water molecules having almost exactly the same positions in both systems. The large change in electronic energy difference when releasing the constraint on 60 water molecules (−1.3 kcal mol−1/H2O) is rooted mainly in the positions of 3 water molecules for which the difference in positional change between going from 1 + 60 H2O (fixed) to 1 + 60 H2O (without constraints) and going from 1-μ1(a) + 60 H2O (fixed) to 1-μ1(a) + 60 H2O (without constraints) is between 2 and 3 Å. All 3 water molecules are within a very short distance of the bidentate acetate ligand. Hence, removing certain water molecules from the systems containing 67 water molecules had a similar effect on the heterocubane structures of both 1

Table 3. Electronic Energy Differences in kcal/mol between 1-μ1(a) and 1 Calculated with Different Solvation Models from Geometry Optimization ΔE 68 H2O 68 H2O + SCRF 68 H2O (QM) + 260 H2O (MM) 45 H2O + COSMO COSMO, no water

no. of H2O

10.9 −13.1 14.9 8.7 −6.9 E


Article


accompanying the deletion of another 10 water molecules far from the bidentate ligand under investigation, which is also reflected in the structures of the systems: While the molecular geometry of 1 changes only slightly (hmp-ligand), the one of 1μ1(a) is more strongly affected by the reduction of the number of solvent molecules. The position and orientation of one monodentate acetate ligand, the hmp-ligands, as well as the partially detached bidentate acetate ligand are significantly changed. The positions of the solvent molecules are (except for those in the vicinity of the partially detached bidentate acetate) similar for 1-μ1(a) and 1 when surrounded by 45 water molecules. Removing 10 of them destabilizes 1-μ1(a) more than 1, and when we relaxed the geometries of the systems without constraints, the water shells of both of them underwent significant rearrangements and diverged structurally from each other going along with the above-mentioned changes of the ligand positions. With less “dense” water shells (i.e., fewer hydrogen bonds per water molecule) making water molecules more “mobile”, it appears that the removal of certain ones of them leads to changes in the hydrogen bond network that accumulate and propagate to the solvent molecules close to the bidentate acetate ligand enabling them to transfer their positional differences between 1-μ1(a) and 1 to the whole solvation shell. All this strongly points toward the fact that the solvation shells we tested are not sufficient to accurately calculate energy differences even of a system involving relatively small structural changes. An improved description of the reaction energies may be obtained if an average over a large number of optimized shapshots is taken. Due to the high computational effort, we did however not investigate this approach further. 2. Reaction Pathways and Barriers for a Selected Example. To gain further insight and calculate not only electronic energy differences but also barriers and free energies, we investigated the ligand exchange of a monodentate acetate ligand of 1 employing NEB and metadynamics. We set up 4 different initial guesses for minimum energy paths: In the first, a water molecule is held in place at a distance of 10 Å from the cubane’s center of mass, while one of the monodentate acetate ligands is moved 11 Å away from the cubane, followed by an attack of the water on the empty coordination site of cobalt (labeled “H2O”, in Figure 3). Next, not only did we set up the same path with hydroxide instead of water (labeled “OH− (a)” in Figure 3), but we also generated an initial guess for a third pathway containing the attack of OH− on Co subsequently displacing the acetate ligand (labeled

and 1-μ1(a) but resulted in different changes of the minimum electronic energy conformations of the water shells. The same is true for most systems, except the ones with 35 or 25 water molecules, where another effect can be seen: Deleting water molecules of the systems including 45 and 35 of them and optimizing the resulting structures with the solvent fixed in space leads to large changes in electronic energy difference of −3.3 kcal mol−1/H2O and −1.5 kcal mol−1/H2O, respectively. For 1 + 35 H2O (fixed solvent coordinates), the geometrical change of the cubane-structure after deleting the 10 water molecules is located mainly on the bidentate acetate and consists of a rotation of the methyl group and slight displacement of the whole acetate ligand, whereas no such changes can be observed for 1-μ1(a) + 35 H2O (fixed solvent coordinates). Thus, removal of a certain water molecule close to the bidentate acetate ligand leads to a change in electronic energy difference already for the heterocubane structure alone. The same is true for the systems with 25 water molecules with fixed coordinates, where this time, the difference in configurational change is located mostly on the hmp-ligand coordinated to the cobalt atom from which the bidentate acetate ligand is detached in 1-μ1(a). To summarize, it is evident that the presence or absence of water molecules close to the site of the largest structural difference has a significant impact possibly not only on the difference in configuration of the water shells but also on one of the cubane structures. To investigate, whether it might be possible to find a more obvious trend for a reasonable number and position of explicit solvent molecules, we repeated the calculations going from 55 to 35 solvent molecules, but this time always deleting the 10 water molecules with the biggest positional difference between the 1 and the 1-μ1(a) system, which are, as mentioned above, the solvent molecules closest to the bidentate acetate ligand. Next, again going from 55 to 35 solvent molecules, we randomly selected water molecules for removal but only among those far away from the site of the biggest structural difference (see Table 5). Table 5. Electronic Energy Differences in kcal/mol between 1-μ1(a) and 1 with Water Molecules Deleted Close to and Far from the Bidentate Acetate Liganda solvent fixed close/ far

close far

no constraints

no. of H2O

ΔE


ΔE


55 45 35 45 35

7.1 −0.9 −8.0 12.9

1.2 −0.4 −0.3 1.8

−5.4 3.3 −12.4 −5.5 9.4

−0.9 −0.4 −1.2 0.3 −0.4

In the first column, it is reported whether water molecules close or far from the (partially detached) bidentate acetate were deleted. The second column gives the number of remaining water molecules. a

Deleting solvent molecules close to the bidentate acetate influences the structures of 1-μ1(a) and 1 unequally and significantly changes the electronic energy difference. Nevertheless, this change is not necessarily larger than the one observed for the random selection of water molecules (compare Table 4). Removing 10 of 55 water molecules far from the partially detached acetate entails only minor structural rearrangements and no change in electronic energy difference. This picture is changed when looking at the changes

Figure 3. NEB energy profiles. F


Article

Journal of Chemical Theory and Computation “OH− (b)” in Figure 3). Finally, we displaced Ac− without filling the resulting vacancy on Co (labeled “Ac−“). None of these systems contained any solvent molecules. The electronic energy difference between the first and last replica of the ligand exchange of acetate with water is similar to the electronic energy difference obtained with the DFT-MD simulation DFT-MD 1 (See Table 1 and Figure 3.) where in fact, as already mentioned in the beginning of the Results and Discussion section, for the most part of the trajectory, water and not OH− was attached. We also calculated the energy difference for the exchange of a monodentate acetate with water solvated with COSMO and found 1 to be only 10.1 kcal/mol more stable than 1-H2O. The energy barrier of almost 60 kcal/ mol would make it virtually impossible for the ligand exchange to occur at room temperature. However, it should be kept in mind that solvation is most certainly a very important issue we did not account for in our NEB calculations, especially for the exchange of charged species, as well as shortcomings of the exchange-correlation density functional used. For our systems at hand, one can expect an approach thoroughly considering solvation and its dynamics as well as entropic effects to yield exact results. As usual, such methods come at a significant computational cost. The drawback of the neglect of solvent effects becomes more evident when looking at the two pathways involving hydroxide. Their electronic energy difference between starting and end frame is similar to the electronic energy difference obtained with COSMO (see Table 1). The huge discrepancy between them and the electronic energy difference between the first and last replica of the ligand exchange of acetate with water presumably stems from the Coulomb interaction between the positively charged {1-H2O}+ and the acetate anion in the last frames of the “H2O-pathway” and the absence of charge screening by solvent molecules or a solvent continuum. The rather good agreement of DFT-MD sampling and in vacuo electronic energy difference between 1-H2O and 1 is presumably due to error cancelation. The reaction barrier for OH− (a) is similar in height to the one for replacing acetate with water, since in both cases the increase in electronic energy with respect to the first frame corresponds mainly to the dissociation of Ac −, which is further corroborated by investigating the minimum energy path of removing acetate without replacing it with any other ligand (see Figure 3). The beginning of the NEB pathway OH− (b) goes along with a decrease in energy. Since frame 1 was obtained from a geometry optimization and should therefore represent a local electronic energy minimum (the “OH− (a)” pathway, which uses the same first frame continuously increases in electronic energy to a maximum), there must be some small barriers invisible on the scale of the NEB calculation. This is however not the only issue with this pathway: after the hydroxide has moved in and attacked the Co, the acetate is slightly displaced, now only hydrogen bonding with one oxygen atom to the nearby water ligand and with the other one to the newly attached OH-ligand. Before it moves further away, hydrogenbonding with its oxygen atom previously close to the OHligand to the water ligand and the other oxygen sticking out into vacuum, the system goes through a high-energy transition state structure with a stretched C−O bond and a distorted O− C−O angle. The inability of the NEB procedure to arrive at a reasonable transition state structure points toward the issue of solvation: A few water molecules around the ligand exchange site could have facilitated the removal of the acetate and kept it

from binding so tightly to the water ligand or the newly attached hydroxide. This system was further investigated by a NEB consisting of 32 images between the third replica before and the third one after the saddle point of the above-mentioned NEB (see Figure 3). The observations stayed qualitatively the same, with the height of the barrier going from 111.5 kcal/mol for the “coarse approach” to 159.9 kcal/mol for the one exploring the saddle point region in more detail. In order to investigate the free energy surface corresponding to the ligand exchange of acetate with water, we prepared a metadynamics simulation with 8 water molecules placed in close proximity of the monodentate acetate ligand. We defined two CVs: The coordination of the cobalt center to the oxygen of the acetate initially bound to it (CV1) and the coordination of oxygen of any water molecule to that Co (CV2). To ensure that the 8 water molecules did not move beyond a certain distance from the cubane, we applied a quartic wall potential to CV2 in the direction of a smaller coordination number at position 0.1. The exact setup parameters can be found in the SI, Table S1. Starting in the local free energy minimum a (see Figure 4), CV1 decreased from 0.9 to 0.1, describing the loss of acetate (region b), an almost barrier-free process. Next, and hence following a two-step mechanism, CV2 increased from 0.2 to 1

Figure 4. Metadynamics free energy surface and corresponding selected minimum structures of 1 solvated with eight water molecles. G


Article

Journal of Chemical Theory and Computation with again almost no free energy barrier, and a water molecule occupied the vacant coordination site (c). The obtained free energy difference between regions a and c (−24 kcal/mol) differs remarkably from the electronic energy difference obtained from NEB or DFT-MD simulations (see Table 1 and Figure 3). This can surely be attributed not only to entropy and sampling effects but also to the level of solvation. For example, in all NEB calculations, after breaking the Co−OAc bond, the acetate first hydrogen-bonded with one oxygen to the water ligand close by, next bonded with the other, and only then was able to drift off into vacuum. During the metadynamics trajectory, on the other hand, it broke its bond to cobalt and hydrogen-bonded with one oxygen to the water ligand as above; with the other one, however, it coordinated to a solvent water molecule thus facilitating its departure and the charge separation going along with it (Figure 4). During the 8 ps necessary to adequately sample the relevant parts of phase space, a fourth local free energy minimum was discovered (d) corresponding to the coordination of a second water and the partial disintegration of the cubane cage. This region of phase space is only 1 kcal/mol higher in free energy than 1 but is accessed via a barrier of 28 kcal/mol. The huge discrepancy between free energy difference and electronic energy difference obtained from DFT-MD simulations sampling only the initial and final state warrants closer scrutiny. To this end, we carried out another metadynamics study, this time, however, in a periodic box containing (in addition to the solute) not only 8 but also 328 water molecules plus a hydroxide anion. This is the same system as the ones used in DFT-MD1 and DFT-MD2 (see Table 1). The simulation was run for 49 ps with the setup parameters given in Table S1 in the SI and further described in the following paragraphs. In the final setup of the full solvation metadynamics study, we defined as CV1 the coordination of both oxygen atoms of the monodentate acetate ligand to the Co center and chose a “step-function-like” coordination function (see SI, Table S1). Reasoning that the attachment of water (or the hydroxide ion) should take place without the need to enhance sampling in this dimension, we chose as CV2 the coordination of the hydrogen atoms of any solvent water molecule to the oxygen atom of the acetate initially forming the bond with Co thereby “describing” its solvation and through that perhaps facilitating the acetate’s departure. We decided not to include the coordination to the other acetate oxygen to avoid excessive fluctuations. For the remainder of this paper, we will label this atom “O2’, and we will label the one which is part of CV2 “O1”. The local free energy minimum corresponding to the attached acetate is labeled a in Figure 5. CV1 takes a value of approximately 1.7, and CV2 is small since O1 is not easily accessible to the solvent. A rearrangement of the water shell takes place with a free energy barrier of 6 kcal/mol, and CV2 is slightly decreased; whereas the Co−O coordination is virtually unchanged. The acetate changes however its orientation, and O2 no longer interacts with the water ligand. From this point, which has approximately the same energy as the initial configuration, the system drops into basin b. In this region, O1 interacts with cobalt as well as the water ligand, and the Co−O bond is strongly elongated. O1 is even less accessible to the solvent, and CV2 still takes small values. O2, on the other hand, is (partially) solvated. This free energy basin is 10 kcal/ mol lower than region a.

Figure 5. Metadynamics free energy surface of the exchange of a monodentate acetate with water in 1. The system includes 1,328 solvent water molecules and one hydroxide ion. (a) acetate attached to the Co center via O1; (b) acetate interacting with water ligand and the Co center via O1; (c) acetate hydrogen-bonding to water ligand via O1; (d) acetate hydrogen-bonding to water ligand via O2; (e) acetate in solution; (f) high energy conformation of acetate attached to Co via O2 (for details, see text).

From here, crossing a barrier of 4 kcal/mol, the free energy well c is accessed, where O1 is only interacting with the water ligand close by and is already far (approximately 5 Å) from the cobalt center. This thermodynamically very stable configuration (24 kcal/mol lower than the initial structure) is somewhat reminiscent of basin c of the above-mentioned metadynamics study with only 8 solvent molecules. It has to be mentioned, though, that at this point no water molecule has occupied the empty coordination site on cobalt, yet. Next, the system crosses a barrier of 26 kcal/mol into region d, where the acetate is “flipped” and now interacts with the water ligand via O2, with O1 pointing into the solution. Since O1 is now (partially) surrounded by solvent molecules, CV2 is increased. Finally, over a 48 kcal/mol barrier, the acetate goes into solution, or, more precisely, the interactions between the water ligand and O1 or O2 are broken; and the acetate faces the cubane with its methyl group. Both oxygen atoms have a distance from the cobalt center of approximately 7 Å. This part of the free energy surface (e) is a very shallow minimum, which is, however, not due to insufficient sampling or premature abortion of the metadynamics simulation. When we stopped the calculations after 49 ps, the behavior of CV1 had already started to become diffusive, and the system was at a point in phase space corresponding again to an attached acetate (i.e., a large value of CV1). Unfortunately, during the entire simulation, water (or hydroxide) did not coordinate to cobalt. Contrary to our initial assumption, there is probably a barrier associated with that process, which therefore could not take place without enhancing sampling in the corresponding dimension. The bulky methyl group pointing toward the cobalt center might make it sterically difficult for water to insert itself between the acetate and the heterocubane, and the negative charge of the acetate keeps the hydroxide ion from approaching. Occupying the empty coordination site might certainly change the free energy of the final state corresponding to the detached acetate. The H


Article

Journal of Chemical Theory and Computation thermodynamic energy difference between the initial state a and that final state e, 34 kcal/mol, is similar to the electronic energy difference found with nonenhanced DFT-MD runs and reported in Table 1. However, this free energy difference is expected to be lowered further by coordination of hydroxide to the cobalt center. Finally, it should be noted that basin d would in principle (not following a minimum free energy path) be directly accessible from a via a barrier of 12 kcal/mol. This would correspond to O2 maintaining its interaction with the water ligand and O1 detaching from Co and being solvated. Region f corresponds to the acetate being coordinated to cobalt via O2. In these configurations, O1 does not interact with the water ligand but is surrounded by the solvent environment. To summarize, the results and high barrier from the metadynamics are in line with previous results,2 which suggested that the monodentate acetate ligand of 1, even though it might detach from the cobalt center and interact with the cubane via a water ligand, is not likely to leave the catalyst.

desirable to arrive at a very accurate description of such complex heterocubanes.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.7b00214. Analysis of the bulk properties of water, summary of the theoretical background of the metadynamics method, and metadynamics parameters used in this work (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Sandra Luber: 0000-0002-6203-9379 Funding

We thank the National Center of Competence in Research − Materials Revolution: Computational Design and Discovery of Novel Materials (NCCR-MARVEL) for financial support and the Swiss National Supercomputing Center (project ID: s502) for computing resources. The work has been supported by the University Research Priority Program “Light to Chemical Energy Conversion” (LightChEC). Funding by Forschungskredit of the University of Zurich, grant no. FK-15-101, is gratefully acknowledged.

■

CONCLUSIONS The calculations of ligand exchange reaction energies presented in this work clearly demonstrate the general importance of adequate solvent representation, especially for reactions involving significant structural changes. Energetics obtained from geometry optimizations with a popular solvent continuum model (COSMO) or explicit solvation shells appeared to be unable to sufficiently capture all effects determining the reaction energies, which change drastically depending on the method employed. Even for large solvation shells and reactions involving a comparably small rearrangement of the molecular geometry (e.g., partial dissociation of a bridging acetate ligand), the results show a strong dependence on the number of explicit solvent molecules and the initial guesses for the geometry optimizations. Hence, DFT-MD simulations at ambient conditions with full solvation, i.e. enough solvent molecules to retain bulk water properties far from the solute, provide in general a better description. Due to the high computational cost of this approach and the large size of the systems under investigation, it can however be difficult to obtain sufficiently converged averages. Beyond that, the COSMO-RS calculations demonstrate that entropy may have a decisive influence on the thermodynamics of ligand-exchange reactions, which are not captured by solely considering electronic energies. For calculations of barriers and reaction paths, the situation is similar. Neither static approaches, such as NEB calculations in vacuo, nor metadynamics with only a few water molecules included can provide a complete picture of the investigated systems. The results from a sufficiently solvated system used for metadynamics showed the vital influence of solvation on free energies and corresponding entropic contributions. DFT-MD calculations are too expensive for routine calculations leaving static calculations with solvent continuum models as the most common method to approximately include solvent effects. A sound, yet computationally inexpensive way to include solvent effects into the modeling of complex cubane catalysts is thus preferable. Future work may additionally include the assessment of the self-consistent Direct-COSMORS68 method, which has, among others, been successfully used for Diels−Alder reactions 69 and organometallic compounds.70,71 Moreover, approaches beyond DFT for an improved description of the electronic structure are highly

Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Peter Deglmann acknowledges hospitality and many fruitful discussions upon a sabbatical with the group of Jürg Hutter at the University of Zurich.

■

REFERENCES

(1) Young, K. J.; Martini, L. A.; Milot, R. L.; Snoeberger, R. C.; Batista, V. S.; Schmuttenmaer, C. A.; Crabtree, R. H.; Brudvig, G. W. Light-driven water oxidation for solar fuels. Coord. Chem. Rev. 2012, 256, 2503−2520. (2) Evangelisti, F.; Moré, R.; Hodel, F.; Luber, S.; Patzke, G. R. 3d− 4f {CoII3Ln(OR)4} Cubanes as Bio-Inspired Water Oxidation Catalysts. J. Am. Chem. Soc. 2015, 137, 11076−11084. (3) Hodel, F. H.; Luber, S. What Influences the Water Oxidation Activity of a Bioinspired Molecular CoII4O4 Cubane? An In-Depth Exploration of Catalytic Pathways. ACS Catal. 2016, 6, 1505−1517. (4) Hodel, F. H.; Luber, S. Redox-Inert Cations Enhancing Water Oxidation Activity: The Crucial Role of Flexibility. ACS Catal. 2016, 6, 6750−6761. (5) Hohenberg, P.; Kohn, W. Inhomogeneous Electron Gas. Phys. Rev. 1964, 136, B864−B871. (6) Kohn, W.; Sham, L. J. Self-Consistent Equations Including Exchange and Correlation Effects. Phys. Rev. 1965, 140, A1133− A1138. (7) Jónsson, H.; Mills, G.; Jacobsen, K. W. Nudged Elastic Band Method for Finding Minimum Energy Paths of Transitions. In Classical and Quantum Dynamics in Condensed Phase Simulations; Berne, B. J., Cinotti, G., Coker, D. F., Eds.; World Scientific: NJ, 1998; pp 385−404, DOI: 10.1142/9789812839664_0016. (8) Laio, A.; Parrinello, M. Escaping free-energy minima. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 12562−12566. (9) Jacobsen, H.; Cavallo, L. On the Accuracy of DFT Methods in Reproducing Ligand Substitution Energies for Transition Metal Complexes in Solution: The Role of Dispersive Interactions. ChemPhysChem 2012, 13, 562−569. I


Article

Journal of Chemical Theory and Computation (10) Jensen, J. H. Predicting accurate absolute binding energies in aqueous solution: thermodynamic considerations for electronic structure methods. Phys. Chem. Chem. Phys. 2015, 17, 12441−12451. (11) Matosziuk, L. M.; Holbrook, R. J.; Manus, L. M.; Heffern, M. C.; Ratner, M. A.; Meade, T. J. Rational Design of [Co(acacen)L2]+ Inhibitors of Protein Function. Dalton Trans. 2013, 42, 4002−4012. (12) Eilmes, A. Solvatochromic probe in molecular solvents: implicit versus explicit solvent model. Theor. Chem. Acc. 2014, 133, 1538. (13) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Sorting Out the Relative Contributions of Electrostatic Polarization, Dispersion, and Hydrogen Bonding to Solvatochromic Shifts on Vertical Electronic Excitation Energies. J. Chem. Theory Comput. 2010, 6, 2829−2844. (14) Caricato, M.; Lipparini, F.; Scalmani, G.; Cappelli, C.; Barone, V. Vertical Electronic Excitations in Solution with the EOM-CCSD Method Combined with a Polarizable Explicit/Implicit Solvent Model. J. Chem. Theory Comput. 2013, 9, 3035−3042. (15) Bachrach, S. M. Microsolvation of Glycine: A DFT Study. J. Phys. Chem. A 2008, 112, 3722−3730. (16) Huang, Z.; Dai, Y.; Wang, H.; Yu, L. Microsolvation of aminoethanol: a study using DFT combined with QTAIM. J. Mol. Model. 2011, 17, 2781−2796. (17) Amaro-Estrada, J. I.; Maron, L.; Ramirez-Solis, A. Aqueous solvation of HgClOH. Stepwise DFT solvation and Born-Oppenheimer molecular dynamics studies of the HgClOH-(H2O)24 complex. Phys. Chem. Chem. Phys. 2014, 16, 8455−8464. (18) De, S.; Ali, Sk. M.; Ali, A.; Gaikar, V. G. Micro-solvation of the Zn2+ ion − a case study. Phys. Chem. Chem. Phys. 2009, 11, 8285− 8294. (19) Sundararaman, R.; Gunceler, D.; Arias, T. A. Weighted-density functionals for cavity formation and dispersion energies in continuum solvation models. J. Chem. Phys. 2014, 141, 134105. (20) Ozkanlar, A.; Clark, E. A. Sensitivity of the properties of ruthenium “blue dimer” to method, basis set, and continuum model. J. Chem. Phys. 2012, 136, 204104. (21) Bandaru, S.; English, N. J.; MacElroy, J. M. D. Density functional theory calculations of catalytic mechanistic pathways for the formation of O2 involving triazolylidene iridium complexes. New J. Chem. 2014, 38, 4060−4070. (22) Ho, J.; Klamt, A.; Coote, M. L. Comment on the Correct Use of Continuum Solvent Models. J. Phys. Chem. A 2010, 114, 13442− 13444. (23) Bianco, R.; Hay, P. J.; Hynes, J. T. Theoretical Study of O−O Single Bond Formation in the Oxidation of Water by the Ruthenium Blue Dimer. J. Phys. Chem. A 2011, 115, 8003−8016. (24) Mennucci, B. Continuum Solvation Models: What Else Can We Learn from Them? J. Phys. Chem. Lett. 2010, 1, 1666−1674. (25) Kolar, M.; Fanfrlik, J.; Lepsik, M.; Forti, F.; Luque, F. J.; Hobza, P. Assessing the accuracy and performance of implicit solvent models for drug molecules: conformational ensemble approaches. J. Phys. Chem. B 2013, 117, 5950−5962. (26) Brüssel, M.; Zahn, S.; Hey-Hawkins, E.; Kirchner, B. Theoretical Investigation of Solvent Effects and Complex Systems: Toward the calculations of bioinorganic systems from ab initio molecular dynamics simulations and static quantum chemistry. Adv. Inorg. Chem. 2010, 62, 111−141. (27) Marx, D.; Hutter, J. Ab Initio Molecular Dynamics: Basic Theory and Advanced Methods, 1st ed.; Cambridge University Press: Cambridge, NY, 2009. (28) CP2K Developers Group. http://www.cp2k.org (accessed June 7, 2017). (29) Ahlrichs, R.; Bär, M.; Häser, M.; Horn, H.; Kölmel, C. Electronic structure calculations on workstation computers: The program system Turbomole. Chem. Phys. Lett. 1989, 162, 165−169. (30) VandeVondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J. Quickstep: Fast and accurate density functional calculations using a mixed Gaussian and plane waves approach. Comput. Phys. Commun. 2005, 167, 103−128.

(31) Becke, A. D. Density-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. A: At. Mol. Opt. Phys. 1988, 38, 3098−3100. (32) Perdew, J. P. Density-functional approximation for the correlation energy of the inhomogeneous electron gas. Phys. Rev. B: Condens. Matter Mater. Phys. 1986, 33, 8822−8824. (33) Bill, E.; Bothe, E.; Chaudhuri, P.; Chlopek, K.; Herebian, D.; Kokatam, S.; Ray, K.; Weyhermüller, T.; Neese, F.; Wieghardt, K. Molecular and Electronic Structure of Four- and Five-Coordinate Cobalt Complexes Containing Two o-Phenylenediamine- or Two oAminophenol-Type Ligands at Various Oxidation Levels: An Experimental, Density Functional, and Correlated ab initio Study. Chem. - Eur. J. 2005, 11, 204−224. (34) Ames, W.; Pantazis, D. A.; Krewald, V.; Cox, N.; Messinger, J.; Lubitz, W.; Neese, F. Theoretical Evaluation of Structural Models of the S2 State in the Oxygen Evolving Complex of Photosystem II: Protonation States and Magnetic Interactions. J. Am. Chem. Soc. 2011, 133, 19743−19757. (35) Orio, M.; Pantazis, D. A.; Neese, F. Density functional theory. Photosynth. Res. 2009, 102, 443−453. (36) Neese, F. J. A critical evaluation of DFT, including timedependent DFT, applied to bioinorganic chemistry. JBIC, J. Biol. Inorg. Chem. 2006, 11, 702−707. (37) Neese, F. Prediction of molecular properties and molecular spectroscopy with density functional theory: From fundamental theory to exchange-coupling. Coord. Chem. Rev. 2009, 253, 526−563. (38) Neese, F.; Ames, W.; Christian, G.; Kampa, M.; Liakos, D. G.; Pantazis, D. A.; Roemelt, M.; Surawatanawong, P.; Shengfa, Y. E. Dealing with complexity in open-shell transition metal chemistry from a theoretical perspective: reaction pathways, bonding, spectroscopy, and magnetic properties. Adv. Inorg. Chem. 2010, 62, 301−349. (39) Bühl, M. Molecular Dynamics of a Vanadate−Dipeptide Complex in Aqueous Solution. Inorg. Chem. 2005, 44, 6277−6283. (40) 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. (41) Siegbahn, P. E. M.; Blomberg, M. R. A.; Chen, S.-C. Significant van der Waals Effects in Transition Metal Complexes. J. Chem. Theory Comput. 2010, 6, 2040−2044. (42) Wang, J.; Roman-Perez, G.; Soler, J. M.; Artacho, E.; FernandezSerra, M.-V. Density, structure, and dynamics of water: the effect of van der Waals interactions. J. Chem. Phys. 2011, 134, 024516. (43) Schmidt, J.; VandeVondele, J.; Kuo, I.-F. W.; Sebastiani, D.; Siepmann, J. I.; Hutter, J.; Mundy, C. J. Isobaric−Isothermal Molecular Dynamics Simulations Utilizing Density Functional Theory: An Assessment of the Structure and Density of Water at Near-Ambient Conditions. J. Phys. Chem. B 2009, 113, 11959−11964. (44) 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. (45) Jonchiere, R.; Seitsonen, A. P.; Ferlat, G.; Saitta, A. M.; Vuilleumier, R. Van der Waals effects in ab initio water at ambient and supercritical conditions. J. Chem. Phys. 2011, 135, 154503. (46) Goedecker, S.; Teter, M.; Hutter, J. Separable dual-space Gaussian pseudopotentials. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 1703−1710. (47) VandeVondele, J.; Hutter, J. Gaussian basis sets for accurate calculations on molecular systems in gas and condensed phases. J. Chem. Phys. 2007, 127, 114105. (48) Evangelisti, F.; Güttinger, R.; Moré, R.; Luber, S.; Patzke, G. R. Closer to Photosystem II: A Co4O4 Cubane Catalyst with Flexible Ligand Architecture. J. Am. Chem. Soc. 2013, 135, 18734−18737. (49) Wang, P.; Shannigrahi, S.; Yakovlev, N. L.; Hor, T. S. General One-Step Self-Assembly of Isostructural Intermetallic CoII3LnIII Cubane Aggregates. Inorg. Chem. 2012, 51, 12059−12061. (50) Klamt, A.; Schüürmann, G. COSMO: a new approach to dielectric screening in solvents with explicit expressions for the J


Article

Journal of Chemical Theory and Computation screening energy and its gradient. J. Chem. Soc., Perkin Trans. 2 1993, 2, 799−805. (51) Weigend, F.; Häser, M.; Patzelt, H.; Ahlrichs, R. RI-MP2: optimized auxiliary basis sets and demonstration of efficiency. Chem. Phys. Lett. 1998, 294, 143−152. (52) Weigend, F.; Ahlrichs, R. Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy. Phys. Chem. Chem. Phys. 2005, 7, 3297−3305. (53) Bauernschmitt, R.; Häser, M.; Treutler, O.; Ahlrichs, R. Calculation of excitation energies within time-dependent density functional theory using auxiliary basis set expansions. Chem. Phys. Lett. 1997, 264, 573−701. (54) Weigend, F. Accurate Coulomb-fitting basis sets for H to Rn. Phys. Chem. Chem. Phys. 2006, 8, 1057−1065. (55) Brooks, B. R.; Brooks, C. L.; Mackerell, A. D.; Nilsson, L.; Petrella, R. J.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch, S.; Caflisch, A.; Caves, L.; Cui, Q.; Dinner, A. R.; Feig, M.; Fischer, S.; Gao, J.; Hodoscek, M.; Im, W.; Kuczera, K.; Lazaridis, T.; Ma, J.; Ovchinnikov, V.; Paci, E.; Pastor, R. W.; Post, C. B.; Pu, J. Z.; Schaefer, M.; Tidor, B.; Venable, R. M.; Woodcock, H. L.; Wu, X.; Yang, W.; York, D. M.; Karplus, M. CHARMM: the biomolecular simulation program. J. Comput. Chem. 2009, 30, 1545−1614. (56) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79, 926−935. (57) Bussi, G.; Donadio, D.; Parrinello, M. Canonical sampling through velocity rescaling. J. Chem. Phys. 2007, 126, 014101. (58) Eckert, F.; Klamt, A. Fast solvent screening via quantum chemistry: COSMO-RS approach. AIChE J. 2002, 48, 369−385. (59) Klamt, A. The COSMO and COSMO-RS solvation models. Wiley Interdiscipl. Rev. 2011, 1, 699−709. (60) Klamt, A.; Jonas, V.; Bürger, T.; Lohrenz, J. C. W. Refinement and Parametrization of COSMO-RS. J. Phys. Chem. A 1998, 102, 5074−5085. (61) Deglmann, P.; Schenk, S. Thermodynamics of chemical reactions with COSMO-RS: The extreme case of charge separation or recombination. J. Comput. Chem. 2012, 33, 1304. (62) Deglmann, P.; Müller, I.; Becker, F.; Schäfer, A.; Hungenberg, K.-D.; Weiß, H. Prediction of Propagation Rate Coefficients in Free Radical Solution Polymerization Based on Accurate Quantum Chemical Methods: Vinylic and Related Monomers, Including Acrylates and Acrylic Acid. Macromol. React. Eng. 2009, 3, 496−515. (63) Becke, A. D. Density-functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648−5652. (64) 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: Condens. Matter Mater. Phys. 1988, 37, 785−789. (65) Henkelman, G.; Jónsson, H. Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points. J. Chem. Phys. 2000, 113, 9978−9985. (66) Pearlman, D. A.; Rao, B. G. Free Energy Calculations: Methods and Applications. In Encyclopedia of Computational Chemistry; Schleyer, P. v. R., Ed.; John Wiley & Sons: New York, 1998. (67) Lu, N.; Kofke, D. A.; Woolf, T. B. Staging Is More Important than Perturbation Method for Computation of Enthalpy and Entropy Changes in Complex Systems. J. Phys. Chem. B 2003, 107, 5598−5611. (68) Sinnecker, S.; Rajendran, A.; Klamt, A.; Diedenhofen, M.; Neese, F. Calculation of solvent shifts on electronic g-tensors with the conductor-like screening model (COSMO) and its self-consistent generalization to real solvents (direct COSMO-RS). J. Phys. Chem. A 2006, 110, 2235−2245. (69) Theilacker, K.; Buhrke, D.; Kaupp, M. Validation of the DirectCOSMO-RS Solvent Model for Diels−Alder Reactions in Aqueous Solution. J. Chem. Theory Comput. 2015, 11, 111−121. (70) Parthey, M.; Kaupp, M. Quantum-chemical insights into mixedvalence systems: within and beyond the Robin−Day scheme. Chem. Soc. Rev. 2014, 43, 5067−5088.

(71) Renz, M.; Kaupp, M. Predicting the Localized/Delocalized Character of Mixed-Valence Diquinone Radical Anions. Toward the Right Answer for the Right Reason. J. Phys. Chem. A 2012, 116, 10629−10637.

K


Exploring Solvation Effects in Ligand-Exchange Reactions via Static

Recommend Documents