Article pubs.acs.org/accounts
First-Principles Molecular Dynamics Studies of Organometallic Complexes and Homogeneous Catalytic Processes Published as part of the Accounts of Chemical Research special issue “Computational Catalysis for Organic Synthesis”. Pietro Vidossich,* Agustí Lledós,* and Gregori Ujaque* Departament de Química, Edifici C.n., Universitat Autònoma de Barcelona, 08193 Cerdanyola del Vallès, Catalonia, Spain CONSPECTUS: Computational chemistry is a valuable aid to complement experimental studies of organometallic systems and their reactivity. It allows probing mechanistic hypotheses and investigating molecular structures, shedding light on the behavior and properties of molecular assemblies at the atomic scale. When approaching a chemical problem, the computational chemist has to decide on the theoretical approach needed to describe electron/nuclear interactions and the composition of the model used to approximate the actual system. Both factors determine the reliability of the modeling study. The community dedicated much effort to developing and improving the performance and accuracy of theoretical approaches for electronic structure calculations, on which the description of (inter)atomic interactions rely. Here, the importance of the model system used in computational studies is highlighted through examples from our recent research focused on organometallic systems and homogeneous catalytic processes. We show how the inclusion of explicit solvent allows the characterization of molecular events that would otherwise not be accessible in reduced model systems (clusters). These include the stabilization of nascent charged fragments via microscopic solvation (notably, hydrogen bonding), transfer of charge (protons) between distant fragments mediated by solvent molecules, and solvent coordination to unsaturated metal centers. Furthermore, when weak interactions are involved, we show how conformational and solvation properties of organometallic complexes are also affected by the explicit inclusion of solvent molecules. Such extended model systems may be treated under periodic boundary conditions, thus removing the cluster/continuum (or vacuum) boundary, and require a statistical mechanics simulation technique to sample the accessible configurational space. First-principles molecular dynamics, in which atomic forces are computed from electronic structure calculations (namely, density functional theory), is certainly the technique of choice to investigate chemical events in solution. This methodology is well established and thanks to advances in both algorithms and computational resources simulation times required for the modeling of chemical events are nowadays accessible, though the computational requirements use to be high. Specific applications reviewed here include mechanistic studies of the Shilov and Wacker processes, speciation in Pd chemistry, hydrogen bonding to metal centers, and the dynamics of agostic interactions.
1. INTRODUCTION
mechanical/molecular mechanical (QM/MM) approach was developed to account for these effects.2,3 In the fields of organic chemistry and homogeneous catalysis, solvent effects were recognized early to modulate the reactivity of charged or polar fragments, and their modeling required the development of specialized techniques, consisting of a mean field representation of bulk electrostatic effects.4,5 This approach was very successful in many instances. However, recent evidence suggests that microscopic effects modulate molecular properties to an extent that cannot be captured by the continuum approach. Hybrid cluster-continuum models, in which few solvent molecules are explicitly represented, are a first step to account for microscopic solvation. However, solvent molecules at the cluster boundary orient themselves to satisfy point interactions (e.g., Hbonding), hampering the comparison of cluster configurations
Molecular modeling contributes to advancing our understanding of chemical systems at the microscopic scale. Simulations proved their usefulness to rationalize experimental data (chemical reactivity or spectroscopic measurements), providing detailed atomistic information not accessible by experimental means. The challenge is to gain predictive power. This will certainly rely on the accuracy of the level of theory necessary to investigate a chemical problem (computational quantum chemistry methods in the case of chemical reactivity). But accuracy is not the only issue. The point we want to make in this Account is that the model system used in the calculations can be as important as the level of theory.1 The importance of the model was recognized early in the context of biochemistry, where it became evident that it was not possible to unravel the activity of an enzyme without accounting for the steric and electrostatic constraints imposed by the protein framework on the active site. The multiscale quantum © XXXX American Chemical Society
Received: January 31, 2016
A
DOI: 10.1021/acs.accounts.6b00054 Acc. Chem. Res. XXXX, XXX, XXX−XXX
Article
Accounts of Chemical Research differing in the number of such interactions.6 To mimic homogeneous conditions, models including solutes and solvent molecules at the experimental densities can be treated as periodically repeating units removing the explicit continuum (or vacuum) boundary. Handling such extended models requires proper statistical mechanics techniques to sample the many accessible configurations. Molecular dynamics (MD) found widespread application in the context of molecular systems in the liquid state and ab initio MD (AIMD), in which the forces acting on the nuclei are computed from electronic structure calculations, offers the required flexibility to model chemical reactions. The objective of this Account is to highlight areas in which calculations based on an extended model allow proper modeling of chemical events. In the context of this dedicated issue, focused on organometallic complexes and their use in catalytic organic transformations, these areas concern charge transfer and speciation processes. We also include a section on the investigation of the dynamic properties of organometallic complexes.
actual nuclear configuration at each step during the dynamics. In Car−Parrinello MD (CPMD), the electronic orbitals are evolved together with the ions, thus not requiring optimization of the wave function at each MD step.8 Both approaches require considerable computational power and rely on high performance computing resources. Applications reported in this Account addressed reactivity or properties on the ground state electronic surface. Extensions to perform AIMD on exited states have been developed.7 Density functional theory (DFT),9 because of accuracy and performance issues, is the electronic structure method of choice in AIMD simulations. Clearly, AIMD inherits all the limitations of DFT approximations and benefits from its developments (including the treatment of exited states, open shell systems, and dispersion interactions, the latter particularly relevant for the description of molecular assemblies).10−12 Given that force evaluation is required at each step along the simulation, AIMD applications tend to use generalized gradient approximation (GGA) functionals to estimate the exchange−correlation energy rather than the more costly hybrid functionals, which include (a fraction of) exact exchange. For large systems, the multiscale quantum mechanical/molecular mechanical (QM/ MM) approach allows reducing the computational cost compared with a full QM treatment.13 The approach is particularly suited when part of the system is not participating in the reactivity or knowledge of its electronic properties is not required. When dealing with homogeneous systems, solvent molecules may be treated at the QM level, and their diffusion should either be prevented by applying geometrical restraints14 or accounted for by adaptive QM/MM partitioning schemes.15 Temperature and pressure of the system may be controlled.7 Hence, the distribution sampled in an AIMD simulation is representative of the corresponding statistical mechanics ensembles. Because of this, high-energy states (transition states) will be observed with low probability compared with the bottom of the free energy surface, which will be sampled more extensively. Thus, chemical reactivity is most often not observed within the time span of the simulation. To overcome this limitation, enhanced sampling methods have to be used: thermodynamic integration16 and metadynamics17 are particularly suitable to study chemical reactions. In the former, a constraint is applied on a selected variable, which is meant to represent the reaction coordinate (e.g., a distance). In the latter, a time-dependent bias acts on the variable(s) used to drive the reaction. If enough sampling is performed, the potential of mean force (i.e., the free energy profile along the reactive path) may be reconstructed. In both approaches, the time evolution of the system is distorted, and thus the generated trajectories should not be interpreted in a dynamical sense. The applications reported in the following sections were performed with the highly efficient CPMD18 and CP2K19 program packages, both capable of massively parallel calculations. Other codes also include capabilities for performing AIMD simulations.
2. METHODOLOGICAL ASPECTS In the studies outlined in this Account, first-principles molecular dynamics simulations were used as a means to investigate the dynamic and structural properties of organometallic systems (Figure 1). According to this approach, a
Figure 1. Example of model system used in the studies surveyed in this Account. An organometallic complex (colored ball and sticks) is fully solvated (here by water molecules, in white and red sticks) and the primary simulation cell (lower right) is treated under periodic boundary conditions. When available, X-ray structures of the solutes are used to build the model systems and provide a valuable reference to assess the quality of the simulation setup.
numerical solution of Newton’s equations of motion (eq 1) provides the time evolution of the coordinates (Ri) of the atomic nuclei (treated as classical particles with masses mi) under the effect of the interatomic forces Fi:
3. APPLICATIONS The use of explicit solvent models for the study of organometallic complexes and their reactivity is illustrated in the following sections. In sections 3.1 and 3.2, we focus on selected steps of catalytic processes. In section 3.3, investigations of the conformational and solvation properties of organometallic systems are discussed. Most of the studies presented concern aqueous systems, and this certainly reflects
2
mi
d Ri dt 2
= Fi({R i}) = −∇i E({R i})
i = 1...N
(1)
In AIMD, forces are computed on-the-fly from firstprinciples electronic structure calculations.7 Two approaches are in use. In Born−Oppenheimer AIMD, the timeindependent electronic structure problem is solved for the B
DOI: 10.1021/acs.accounts.6b00054 Acc. Chem. Res. XXXX, XXX, XXX−XXX
Article
Accounts of Chemical Research
Figure 2. (a,b) Mechanisms of C−H activation in ref 29. (c) Mechanisms of C−H activation in ref 30. (d) Representative snapshots from an explicit water AIMD simulation showing the C−H bond breakage (first step) and transfer of the resulting proton to the medium (second step). The insets show chemical diagrams of the trans- PtCl2(H2O)(CH4) methane, the PtHCl2(H2O)(CH3) square pyramidal, and the PtCl2(H2O)(CH3) methyl complexes.
Figure 3. (a) Reaction steps of the Wacker process studied by explicit solvent AIMD simulations. (b) Model for the inner- and outer-sphere nuclephilic attack. Reproduced with permission from reference 35.
molecules in the models allows stabilization of the nascent fragments, thus improving the description of these processes. This effect was recognized for reactions in high dielectric media, and indeed aqueous oxidation reactions were among the first to be studied by means of explicit solvent AIMD. In the Shilov and Wacker processes discussed below, a proton is released to the solution. Further applications, not reviewed here, include the Fenton system,21 water oxidation,22−24 hydrogenation,25,26 and reductive dehalogenation reactions.27
the general interest in understanding how the peculiar properties of water (hydrogen-bonding and coordination abilities) influence chemical systems.20 Nevertheless, interest in exploring nonaqueous solutions is increasing, and we will show that AIMD simulations in supposedly inert solvents like toluene may indeed provide new valuable insight. 3.1. Processes Involving Charge Separation
Heterolytic breakage of chemical bonds leads to the formation of charged fragments. The inclusion of explicit solvent C
DOI: 10.1021/acs.accounts.6b00054 Acc. Chem. Res. XXXX, XXX, XXX−XXX
Article
Accounts of Chemical Research
subsequent β-hydrogen elimination took place spontaneously. Finally, formation of the aldehyde via the deprotonation of the OH group of the hydroxyethyl ligand was shown to involve the transfer of the OH proton to a water molecule from the bulk solution. These findings are in agreement with the experimental rate law.39 Interestingly, with a microkinetic model that employs the rate constants derived from simulations, the experimental rate law was reproduced (the rate law shows that the reaction is inhibited by [H+] and [Cl¯]2).40
3.1.1. The Shilov Reaction. Shilov’s system,28 based on Pt and PtIV complexes, allows for hydrocarbon oxidation in water. C−H bond activation is the rate limiting step of the catalytic cycle. Siegbahn and Crabtree considered the transPtCl2(H2O)(CH4) species with one extra water molecule and found two competitive pathways.29 In the first, the C−H breaks protonating one of the Cl ligands and maintaining the square planar arrangement at Pt (ΔE⧧ = 16.5 kcal mol−1, pathway a in Figure 2). In the other, the C−H breaks to form a Pt−hydride, while the water molecule moves from the first solvation shell to coordinate Pt, resulting in an octahedral PtIV complex (ΔE⧧ = 20.0 kcal mol−1, pathway b). Zhu and Ziegler found for the same species a barrierless process that leads to a square pyramidal trans-PtHCl2(H2O)(CH3) hydride (pathway c).30 The authors considered no solvating water molecules in their model. Recently, we reinvestigated the issue.31 Explicitly solvated cisand trans-PtCl2(H2O)(CH4) configurations were considered. CPMD simulations showed the trans-PtCl2(H2O)(CH4) species evolving to a square pyramidal complex after oxidative addition of the C−H bond to the Pt center (Figure 2d). This, however, was not stable and released the apical proton to a close by water molecule, thus forming a square planar methyl complex. The cis complex was observed to be stable on the time scale investigated. Thus, to observe the chemical transformation within a computationally accessible simulation time, we applied a biased simulation technique (namely, metadynamics) to promote the breaking of the C−H bond, which required an activation free energy of 6.4 kcal mol−1. As for the trans isomer, a Pt−H bond formed, which then spontaneously broke transferring a proton to the bulk solution. The trans-effect explains the different behavior observed for the cis and trans σcomplexes in the C−H bond breakage step.30 Remarkably, the use of explicit solvent models revealed the acidity of the hydride intermediate and showed release of a proton to the solution rather than the coordination of a water molecule to stabilize the formally PtIV center. 3.1.2. The Wacker Process. In the Wacker process, ethene is oxidized to acetaldehyde using O2 as oxidant and a Pd catalyst.32 The catalytic cycle includes a hydroxypalladation step, in which nucleophilic addition of a water molecule to the Pd-coordinated alkene takes place. Its mechanism is debated.33,34 One of the questions concerns whether the attack is intramolecular (inner-sphere mechanism), in which a ciscoordinated water attacks the carbon double bond, or intermolecular (outer-sphere mechanism), in which the attacking water comes from the bulk solution (Figure 3). Early computational studies based on cluster models could not solve the dispute. Stirling, Ujaque, and co-workers35,36 and the group of Nair37,38 independently addressed the reaction mechanism using explicit solvent CPMD simulations. The free energy profile along both outer- and inner-sphere mechanisms was estimated, showing that the former is energetically favored. Simulations revealed that in both pathways nucleophilic attack proceeds with the simultaneous transfer of a proton to the solvent. Importantly, analysis of the speciation of the reagent [PdCl4]2− showed that the formation of the trans-[PdCl2(H2O)(C2H4)] species is more stable than the cis species (the latter is required for the inner-sphere mechanism; see also section 3.2.2). The modeling study was further extended to the complete catalytic cycle, showing that the rate-determining step consists of trans-to-cis isomerization of chloride ligands whereas the II
3.2. Speciation
Catalyst speciation includes either a change in the metal coordination sphere, via ligand exchange or loss, or a change in the nature of the ligands, most notably by protonation/ deprotonation events.41,42 Identifying the catalytically active species in the reaction mixture is of great interest for catalyst optimization but, unfortunately, may constitute a challenging problem when species are present at low concentrations or are short-lived. In the simulation studies reported below, the stability of different species in solution was investigated. Further studies investigated ligand exchange reactions for cisplatin43,44 and the uranyl cation.45 3.2.1. Phosphine Dissociation from Pd0 Complexes in Toluene. The 2010 Nobel Prize in Chemistry recognized the relevance of palladium/phosphine catalytic systems in preparative organic synthesis. A highly debated issue in the field of Pd-catalyzed coupling reactions concerns the identification of the catalytically active species, because this is generated in solution from an inactive Pd0 precursor. Experimental and theoretical studies indicate that low coordinate Pd−phosphine complexes react faster in the oxidative addition step, which initiates the catalytic cycle.46 Accordingly, design of bulky phosphine ligands to favor monoligated-Pd species has been a successful strategy to increase catalyst activity.46 Thus, speciation at Pd is a key factor, which depends on the nature of the phosphine ligand but also on the presence of additives. The involvement of solvent was put forward by some authors, yet exchange equilibria involving solvent molecules were not quantitatively addressed. We investigated this issue for the Pd/ PPh3 system by means of QM/MM-MD simulations in explicit toluene.47 Simulations showed the formation of solvent adducts including [Pd(PPh3)2(Tol)], [Pd(PPh3)(Tol)], and [Pd(PPh3)(Tol)2] (Tol = toluene). Interestingly, a toluene molecule was observed to coordinate and decoordinate from Pd(PPh3)2, with the tricoordinated species lying
