Received 20 August 2013
Developing an approach for first-principles catalyst design: application to carbon-capture catalysis
Heather J. Kulik,a,b Sergio E. Wong,a Sarah E. Baker,a Carlos A. Valdez,a Joe H. Satcher Jr,a Roger D. Ainesa and Felice C. Lightstonea*
An approach to catalyst design is presented in which local potential energy surface models are first built to elucidate design principles and then used to identify larger scaffold motifs that match the target geometries. Carbon sequestration via hydration is used as the model reaction, and three- and four-coordinate sp2 or sp3 nitrogen-ligand motifs are considered for ZnII metals. The comparison of binding, activation and product release energies over a large range of interaction distances and angles suggests that four-coordinate short ZnII-Nsp3 bond distances favor a rapid turnover for CO2 hydration. This design strategy is then confirmed by computationally characterizing the reactivity of a known mimic over a range of metal-nitrogen bond lengths. A search of existing catalysts in a chemical database reveals structures that match the target geometry from model calculations, and subsequent calculations have identified these structures as potentially effective for CO2 hydration and sequestration.
The rational design of small-molecule transition-metal catalysts has an impact in fields ranging from energy science (Baldwin & Pecoraro, 1996; Balzani et al., 1998) to pharmaceuticals (Blaser et al., 2001). The sequestration of carbon dioxide (CO2) in industrial processes to reduce emission of greenhouse gases (Pachauri & Reisinger, 2007) is one example of a key challenge that can be addressed through the design of new catalysts. Industrial separation of CO2 typically relies on converting CO2 gas to carbonic acid in a liquid phase (Figueroa et al., 2008) and the subsequent removal of the liquid fraction, where the formation of carbonic acid is rate-limiting (Figueroa et al., 2008). Acceleration of carbonic acid formation through the use of catalysts would improve the effectiveness of these industrial separation processes (Bao & Trachtenberg, 2006). In particular, we focus here on the design of catalysts that could be used in industrial settings, potentially tethered close to the air-water interface to ensure that they are able to capture dissolved CO2 and convert it to carbonic acid.
The enzyme human carbonic anhydrase (CAII) rapidly converts CO2 to bicarbonate at a zinc metal center (turnover 106 s-1 at pH 9 and 298 K; Khalifah, 1971; Bao & Trachtenberg, 2006). The environment around the Zn center is tetrahedral, with three chelating histidine (His) residues and the fourth site occupied by a catalytically relevant hydroxide (Steiner et al., 1975; Lindskog, 1997; Håkansson et al., 1992; Liljas et al., 1972). Additional features of the enzyme include a hydrophobic portion responsible for CO2 binding, and a hydrophilic side that facilitates the binding and deprotonating of water. Several small-molecule CAII mimics exist (Krishnamurthy et al., 2008), which primarily aim to replicate the Zn center and His ligands of the enzyme active site with a Zn metal ligated by comparable nitrogen electron donors, such as imidazoles (Parkin, 2004) and secondary amines (Zhang & van Eldik, 1995; Zhang et al., 1993). Both carbonic anhydrase and its mimics are believed to form bicarbonate via oxidative attack of CO2 by the Zn-OH group. The newly formed bicarbonate is then displaced by a water molecule, which is deprotonated to regenerate the active hydroxide intermediate (Silverman & Lindskog, 1988).
Existing CAII mimics catalyze CO2 hydration at rates much slower than the native enzyme and, in this work, we focus on developing a new strategy for identifying CAII mimics that have the potential to exceed the turnover rates of those mimics already reported in the carbon-capture literature. These activity results suggest the exploration of alternative transition metals as an avenue to modulating carbonic anhydrase mimic behavior, but the clear advantage of Zn in particular, both natively and in mimics, is its relative inertness against oxidative stress (Powell, 2000).
The computational design of catalysts can augment experimental efforts to design and screen new catalyst mimics (Guidoni et al., 2004). While computational approaches are much cheaper than comparable experimental screens (Nørskov et al., 2009), the current limited accuracy of the electronic structure methods that can be applied to studying the system sizes typical for molecular catalysts has curtailed the application of computational approaches. Nevertheless, locality in transition metal chemistry (the strong dependence of the electronic structure of a transition metal center on the character of directly bonded atoms) has long been exploited (Orgel, 1961) to study smaller model catalysts than are typically used in industrial settings. Shortcomings in the accuracy of the chosen electronic structure method may be ameliorated somewhat by focusing on relative trends, which are more robust from error cancellation, rather than on absolute values.
We believe that methods that examine relative trends in small model systems with minimal ligands, once verified against larger realistic catalysts, should provide a wealth of information. This information may then be used for searches of existing catalysts in structural databases or for the creation of never-before-studied catalysts that have desirable characteristics. Carbon capture provides an excellent test case for computational design methods for two principle reasons: (a) there is a significant difference in the activity of previously characterized model catalysts and native enzymes, and (b) there are relatively few apparent reaction steps that will need to be studied and optimized. If optimization strategies for catalysis can be obtained from the computed electronic structure trends on small models, then computation will be better able to guide catalyst selection for experimental study in a variety of relevant chemical reactions.
Here, we present a computational approach to designing catalysts for carbon capture. We have developed extensive potential-energy-based reactivity maps for key components of the CO2 hydration reaction using both three- and fourfold-coordinated metal-nitrogen complexes (sp2 or sp3 nitrogen) with both ZnII metals. We present the results of optimization strategies for the primary aspects of the CO2 hydration reaction: hydroxyl intermediate formation, CO2 activation and conversion, and product release. We identify the metal-ligand sets that provide the best activity at equilibrium geometries, as well as the sets that are most easily optimized, as determined by characteristics all improving along the same vector of geometric change. We verify our minimal model by constraining an experimentally characterized carbonic anhydrase mimic and studying changes in its electronic structure under the constraint to check for correlation against the minimal model potential-energy surface (PES) maps. Finally, we characterize catalysts from the literature that both meet and fail to meet the geometric targets obtained from our PES and compare the resulting catalytic properties against a previously characterized carbonic anhydrase mimic.
Density functional theory calculations were carried out using the Quantum-ESPRESSO package (Giannozzi et al., 2009). An ultrasoft plane-wave pseudopotential approach was used with a cut-off of 30 Ry for the wavefunction and 300 Ry for the density, which was chosen by converging the relative energy between two intermediates to <0.01 eV. A similar test of convergence of forces and relative energetics of intermediates was used to determine that a cell size of 26 Bohrs was suitable for minimal-model PES calculations. For larger catalysts, 8 Å of vacuum was added to the total van der Waals radius of the molecule, giving a final cell size in the range of about 30-50 Bohrs. The pseudopotential for Zn included both semi-core 3d and valence 4s electrons in the valence. Plane-wave DFT is used for both small and large calculations in this study because of the favorable scaling of this approach for the large-scale calculations, thus enabling us to characterize rapidly and consistently both the PES and the larger catalysts. Additionally, with plane-wave calculations, we circumvent issues of extrapolating to the complete basis set limit and issues with basis set superposition error that are likely to plague the binding and release steps of the reaction. The Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional (Perdew et al., 1996) was used for all calculations, since previous work (Kulik & Marzari, 2010) showed that augmenting functionals with a Hubbard U term has little effect on Zn. For Co, an average U of 3.0 eV obtained from the linear response across all intermediates was applied. We note that this Hubbard U term is not a fitting parameter but a calculated measure of the extent of self-interaction error present in the system. Our previous study of DFT+U has shown that, for transition metals, DFT+U systematically reduces the energetic errors more consistently than hybrid functionals (Kulik et al., 2006), which require intermediate-specific tuning of the exact amount of exchange in order to achieve results that are consistent with experiment.
Potential-energy surfaces (PESs) were calculated via constrained relaxations in which the metal-nitrogen distances and the distance of the metal from the plane of the N atoms were fixed at different values, while all other variables were permitted to relax. In each case, stable intermediates were energy minimized under constraint of both the bond distance (M-L) and the dihedral angle1 formed by the metal with three of its ligands L, while all other degrees of freedom were relaxed. The PES calculations for M = Co or Zn and L = NH3 were carried out over a very wide range of M-L distances (1.75-2.55 Å in 0.05 Å increments) and M-3L dihedral angles (0-55° in 5° increments). Selected additional calculations for M = Zn and L = NH2 were carried out over narrower M-L distance ranges (1.85-2.35 Å in 0.05 Å increments) and M-4L dihedral-angle ranges (0-40° in 5° increments). For reference, the lowest energies or equilibrium points of most intermediates on these PESs reside around 2.1 Å for the distance and 10° for the dihedral angle. The energetics of the squeezed catalyst were determined by fixing only the metal-nitrogen distance and allowing everything else to relax. Transition states were obtained approximately through the use of the nudged elastic band (NEB) method, with climbing image and variable springs to ensure resolution of the geometry and energy near the peak of the minimum energy pathway, and forces on the seven images in the path were converged to <0.04 eV Å-1 (Henkelman et al., 2000). Neither minima nor approximate transition states were vibrationally characterized. Nevertheless, our reaction coordinate results are both quantitatively and qualitatively consistent with previous studies that have characterized the transition states (Koziol et al., 2012). Makov-Payne corrections (Makov & Payne, 1995) to the total energy were used when comparing systems of different charge [+1 overall for most reactants and intermediates, +2 for intermediates with axial water and carbonic acid or with no axial ligand, and charges are shifted net negative by 1 to allow for singlet spin in the case of ZnII-(NH2)3]; if the charge remained the same, the systems were treated without the Makov-Payne correction.
In the limit of low pH, carbonic acid is stable and this assumption likely mimics the actual mechanism for product release (Greenwood & Earnshaw, 1997). For neutral to high pH, carbonic acid species are not in high concentration at equilibrium; nevertheless, short-lived carbonic acid species that dissociate at lower energetic cost are known to form transiently, regardless of the predominance of carbonate or dissolved CO2 at equilibrium (Adamczyk et al., 2009). Additionally, carbonate salts (Greenwood & Earnshaw, 1997) may be more common at high pH as well, and the assumptions we make about the mechanism of release via exchange of neutral carbonic acid and water also hold for the exchange of carbonate salts and water. The motivation for studying neutral carbonic acid is due to the well known poor accuracy in electronic structure methods in describing the dissociation of charged species (bicarbonate is negatively charged and the catalyst is positively charged) (Bally & Sastry, 1997; Braïda et al., 1998; Zhang & Yang, 1998; Exner & von Ragué Schleyer, 2001; Gräfenstein et al., 2004; Petrie & Stranger, 2004; Tozer & De Proft, 2005; Dutoi & Head-Gordon, 2006; Mori-Sánchez et al., 2006; Cohen et al., 2008), despite the fact that carbonic acid is in low concentration in aqueous solution. Because the focus is on trends, the ease with which a water molecule may displace a carbonic acid molecule provides a sufficient relative metric of product release. We also note that some of the challenges for density functional theory (DFT) in correctly describing charge transfer may be alleviated when explicit or implicit solvent is incorporated into the simulations (Lange & Herbert, 2007). No explicit or implicit solvent has been added to the present simulations, because the trends in the reaction coordinate are not likely to be substantially affected, but additional considerations of pH-dependence and explicit solvent effects might be relevant for future work.
The candidate structures of putative carbonic anhydrase mimics were obtained by searching data available in the Cambridge Structural Database (CSD, Version 5.24; Allen, 2002) using the Conquest program (Bruno et al., 2002), with constraints on the number and character of ligands coordinating the metal center, as well as on the metal identity. Multiple searches were carried out for a variety of first-row transition metals, but the distributions of the resulting metal-nitrogen bonds were only compared across a single metal identity at a time. If candidate structures with metal centers other than Zn or Co were obtained, the metal in the crystal structure was replaced and re-optimized alongside results obtained with Zn or Co centers that were also optimized with the native metal.
We have characterized a model reaction coordinate for carbon capture via CO2 hydration in order to identify how to optimize catalysts towards facilitating this reaction. We identified potential reaction steps through what was known about carbonic anhydrases and its mimics, and by comparing reactants and products and possible intermediate species that would need to bind to the metal. For other reactions, we would also carry out a literature search to identify potential competing reaction steps. Additionally, a large focus of our optimization search is on controlling the relative binding energy of different species to the catalytic center. We can identify which species are likely to be present in order to identify which binding energies should be calculated. In this case, the characteristics of the reaction coordinate (see Fig. 1) were identified as
| || Figure 1 |
Schematic of reaction characteristics, labeled as follows: (1) EM-OH = hydroxyl stability (green), (2) Ea = activation energy (red), (3) Erxn = reaction energetics (blue), and (4) De(P) = product release (orange). Approximate structures are provided above each intermediate.
(1) stability of hydroxyl formation and binding
(2) activation energy Ea for M-OH attack on CO2
where the subscript TS denotes transition state,
(3) reaction energetics Erxn for M-OH attack on CO2
(4) product release De(P)
There are a few assumptions made in our definition of reaction characteristics. In equation (1), we assign hydroxyl formation as the relative energy of replacing a water molecule by a hydroxyl and a free hydronium ion. This approximation is justified because trends in hydroxyl formation across geometries are not particularly sensitive to the role and identity of the species that deprotonates. In equation (4), we assume that a neutral carbonic acid molecule is formed in order to aid its release from the catalyst, and that a free water molecule is present and may displace a neutral carbonic acid molecule (see §2, Computational methods).
In our strategy to optimize computationally the catalysts for CO2 hydration, we varied both metal and ligand identity and coordination number, as well as the relative distances between metals and ligands. Namely, we consider both ZnII metals either four- or threefold coordinated by ligands (L = NH3, NH2). ZnII is d10 with no net spin in any configuration. The ligands were chosen to represent sp2 (L = NH2) and sp3 (L = NH3) N atoms, and they coordinate the metal in either a threefold trigonal-planar to tetrahedral-like geometry, or a fourfold square-planar geometry. The tetra-aza scaffold we discuss in detail later is an example of a scaffold that chelates metals with an Nsp3 ligand set. The Nsp2 ligands approximately model ligand sets including imidazole rings or porphines, but we note that it is possible that some features of Nsp2 atoms in aromatic rings might be missed by our minimal model (Togni & Venanzi, 1994). The metals were chosen on the basis of the metal bound in wild-type carbonic anhydrase, as well as the activity of alternative metals when bound to carbonic anhydrase (Coleman, 1967a,b; Vallee & Galdes, 1984; Lipscomb, 1983). For each metal and ligand set, full PESs were obtained for all the intermediates and transition states that make up the reaction coordinate in Fig. 1, by varying both M-L distance (dM-L) and M-3L dihedral angle (M-3L). (Breaking symmetry via alternating short and long bonds was pursued in a preliminary fashion and was not observed to alter broad trends significantly, but further examination may be of interest in future work.) Once these PESs have been analyzed for geometric properties that optimize the reaction coordinate, structures that fit these characteristics from the literature will be fully characterized (Fig. 2). The benefits of carrying out a literature search are twofold: (a) if we confine ourselves initially to previously characterized structures, then the potential catalysts may be straightforwardly synthesized and tested; and (b) without restricting ourselves to carbon-capture catalysts, we may identify structures relevant to other reactions that may be repurposed for CO2 hydration.
| || Figure 2 |
Flow chart indicating the strategy for carbon-capture catalyst design. Starting from the top, first metals and ligands are chosen based on the enzyme, then the dependence of the reaction parameters on the metal-ligand interaction geometry is explored, thirdly reaction-oriented PESs are used as input for a database search, and finally the reaction coordinates of candidates from the search are obtained.
Using the results of the four properties outlined previously at different fixed distances and dihedral angles, we obtained the direction along which the property is optimized when moving from equilibrium (i.e. either shortening or lengthening bonds, or flattening or enhancing the dihedral angle). The results of these optimization directions for all four reaction coordinate characteristics are reported in Table 1 for ZnII-ligand sets. For each of the four reaction-oriented PESs we calculated, we obtained equilibrium values of the four reaction characteristics. Using equilibrium data, we observed that hydroxyl formation [equation (1)] is favorable when we compare the relative binding energies of a hydroxyl species and a water species. We also observe that the CO2 hydration reaction [equation (3)] is always exothermic. The other two points of consideration, the height of the CO2 hydration barrier [equation (2)] and product release [equation (4)], were instead observed to be potentially rate-limiting for all equilibrium points.
We have collected geometric trends for reactivity optimization for each reaction step with a variety of metal-ligand sets (Table 1). Steps that are not rate-limiting [e.g. steps (1) and (3) here] should be weighted less in systems where the optimization directions differ between different characteristics. Additionally, when comparing multiple metal-ligand sets, even if one set appears to be more easily optimized, the favorability of the reaction coordinates at the respective geometric equilibria of the metal-ligand sets should also be compared. Using the two key reaction characteristics, we ranked each metal-ligand set and observed that four-coordinate ZnII with L = NH3 or NH2 started from the best equilibrium characteristics.
Both ZnII-(NH3)4 and ZnII-(NH2)4 compounds are relatively optimizable. In both cases, three out of four reaction characteristics are in the same direction, with the former compound preferring short and flat bonds, while the latter prefers long and domed interactions. We focus here on using the PES and reaction coordinate optimization data obtained for ZnII-(NH3)4, both to better understand the underlying electronic structure that gives rise to these optimization trends and to verify that our minimal model predicts trends that are meaningful for larger catalysts. We reserve comparison of ZnII-(NH2)4 compounds for future study because we are better able to validate the results of the ZnII-(NH3)4 PES optimization against a previously computationally characterized catalyst mimic (Wong et al., 2011). By comparing against previous results on 1,4,7,10-tetraazacyclodedecane (tetra-aza), we are able to verify our minimal ligand approach. Using the tetra-aza scaffold with a ZnII metal center, we carried out structural relaxations for the larger scaffold structures under constrained Zn-N bond distances, and examined how the two principal reaction characteristics, activation energy and product release, change with varying Zn-N bond distance (the Zn-N dihedral angle is not constrained here). Without constraints, the energy-minimized tetra-aza structure has an equilibrium Zn-N bond distance of 2.22 Å. We find that compressing the Zn-N bond in the scaffold from its equilibrium value leads to a decrease in activation energy for CO2 hydration, in excellent agreement with the computational results on model ligands (Fig. 3). This strong agreement verifies our minimal-ligand PES approach. Expansion of the Zn-N bonds to values larger than in the unconstrained relaxed tetra-aza scaffold instead increases the activation energy. The large energy range (5 kcal mol-1; 1 kcal mol-1 = 4.184 kJ mol-1) and distance range (0.5 Å), as well as the fact that trends are being compared, ensures that this test is a robust one for our approach, despite any residual errors in the absolute values of the energetics. We find that the dihedral angles, which we did not constrain, change from 35 to 12° as the scaffold is squeezed from a Zn-N bond distance of 2.3 to 1.8 Å. The geometry of the scaffold corresponds to the flattening of the Zn-N dihedral angle as the lowest energy pathway to shortened Zn-N bonds.
| || Figure 3 |
Comparison of activation energy (kcal mol-1) for squeezed tetra-aza scaffolds (black circles) with the model ligand potential energy surface over the range 1.8-2.3 Å for the Zn-N bond distances d. Blue contour lines represent different dihedral angles, with light-blue shading covering the whole range of PES contours.
As an additional test, we considered whether our minimal-ligand PESs could reproduce trends in product release from the squeezed tetra-aza scaffold. Our ZnII-(NH3)4 results suggested that a shorter Zn-N bond would facilitate product release (see Table 1). Product release has recently been proposed as the rate-limiting step in these small-molecule CAII mimics (Koziol et al., 2012), especially when using a definition of product release that involves the release of bicarbonate (instead of the carbonic acid presented here). Comparing product-release characteristics for squeezed tetra-aza structures revealed a trend similar to that from the minimal-ligand structures (Fig. 4). While the absolute agreement is a little weaker than for the activation energies, the trend is consistent across a range of Zn-N distances. We note that, indeed, at very short Zn-N distances, product release becomes increasingly unfavorable, but for intermediately short distances, product release is still weakly favorable in both our minimal model and in the squeezed scaffold. This justifies the focus on short bonds later in our database search of catalyst scaffolds. Since the optimization characteristics observed in the minimal-ligand set were upheld when studying the tetra-aza scaffold under constraint, we next considered what distance-dependent electronic structure properties were giving rise to these optimization characteristics.
| || Figure 4 |
Comparison of product release (kcal mol-1) for squeezed tetra-aza scaffolds (black circles) with the model ligand potential energy surface over the range 1.8-2.3 Å for the Zn-N bond distances d. Blue contour lines represent different dihedral angles, with light-blue shading covering the whole range of PES contours.
In order to identify the electronic-structure underpinnings of these trends, we have analyzed the population and molecular orbital character in our minimal-ligand ZnII-(NH3)4 set. When the Zn-N bond is squeezed, the largest changes are a decline in the Zn 4s population and a commensurate increase in the N 2s population (Fig. 5). Additionally, the 2s and 2p populations for the axial O atom in the reactant state increase for shorter Zn-N bond distances, while the Zn-N 3d population decreases. The overall population trends suggest that charge is transferred primarily from the Zn metal center to N, likely through increased bonding interactions between Zn 4s and N 2s/2p as the species are squeezed together. We also consider the change in the eigenvalues and eigenstates as the model system is squeezed (Fig. 5). Not surprisingly, the energy window for all states with Zn 3d character increases as the molecule is squeezed. The most strongly Zn-N distance-dependent levels (Fig. 5) have large Zn 3d contributions. Bonding orbitals that have strong interactions between N 2p and Zn 3d are pushed further down in energy upon squeezing, enhancing the energetic overlap with deeper N 2p and N 2s levels. Of the levels that are pushed higher in energy upon squeezing, they all appear to have strong antibonding character between in-plane Zn 3d and N 2p orbitals. Overall, interactions in the plane between Zn 3d/4s and N 2p/2s are strengthened through both charge transfer and better energetic and geometric overlap when the Zn-N bond is squeezed. This strengthening in the plane weakens the axial interactions, facilitating turnover through localizing more charge on the hydroxyl for interaction with CO2 in the transition state of the hydration reaction. Notably, the increase in 2s/2p populations on the axial O atom that is observed upon squeezing provides quantitative support for the chemical picture that CO2 insertion through nucleophilic attack by an increasingly nucleophilic hydroxyl is becoming more favorable.
| || Figure 5 |
Electronic structure dependence on Zn-N bond distance. (a) Relative atomic population change across a range of Zn-N distances in a four-coordinate Nsp3 model compound (Zn s in black circles has been divided by five to show on the same scale as the 3d electrons in red squares). Each curve is set to one at the point where the relative occupancy is the lowest across the range of distances. (b) Variation in eigenvalues across a range of Zn-N distances in a model compound. Rapidly varying levels are highlighted in orange, red, green, and blue, and labeled with their respective atomic contributions.
Following verification of the minimal-ligand PESs, we searched the literature to identify both the spread of Zn-N bonds in experimentally characterized catalysts and the outliers on the tails of this Zn-N distribution. The results of the reaction-oriented PES analysis were used as the basis for identifying structures from the CSD that are potentially suitable candidate catalysts, as well as those that oppose the optimization characteristics. A CSD search may be straightforwardly carried out for a variety of metal-ligand combinations using the coordination number of the metal and the ligand character. In our own searches, we used both native-metal structures (i.e. scaffolds with ZnII centers) and structures that contained other 3d MII metals. We identified outliers as structures with metal-ligand distances and dihedral angles that are two standard deviations below (Fig. 6a) or above (Fig. 6b) the mean. Whether or not the structures were obtained with ZnII experimentally, the metal center was replaced with zinc, and a full structural optimization was carried out in all cases. We verified that these structures still had short or long bonds after optimization and that having bond lengths far from the mean was not simply a result of error from the experimental crystallography. After this point, we then characterized the reaction coordinates fully with no constraints on the CSD-derived catalysts.
| || Figure 6 |
Four-coordinate Nsp3 structures obtained from the database search, with Zn-N distances at least (a) two standard deviations below the mean and (b) two standard deviations above the mean Zn-N separation in the distribution of search results. The structures are labeled with their CSD refcode.
Structures that fulfil the optimization vector for ZnII-(NH3)4 were identified in the distribution of structures obtained from the CSD search, while the opposing tail of the distribution was used for control structures. From each of these tails, we selected two structures based on the correspondence between experimentally observed Zn-N bond lengths being within 1% of the calculated bond length (Table 2).
As a reference, we use the previously characterized tetra-aza scaffold (Satcher et al., 2011) that we also considered in our constrained Zn-N bond studies. We note again that we focus our CSD search on Nsp3 compounds largely because of the available reference data, but Nsp2 compounds will be considered in future work. We compare all four reaction characteristics, including the hydroxyl stability and reaction exothermicity steps. We observe that the short Zn-N bond structures [CSD refcodes QOVZEH (Lin et al., 2008) and FEQYOP (Choi et al., 1998)] are comparable with or improve upon the tetra-aza barrier height characteristics, though product release is slightly less favorable. For the long-bonded structures [OBOJOE (Handley et al., 2001) and MECVEV (Kong et al., 2000)], the activation energies are noticeably higher. In one long-bond Zn-N structure (OBOJOE), the product release becomes favorable, while in the other, product release is unfavorable (MECVEV). This result suggests that steric and electrostatic factors not considered in our PES search or fully modeled by our minimal ligands may come into play in determining these reaction characteristics.
The difference in the results for the two long-bond Zn-N CSD structures becomes apparent in observing the differences between the structures (Fig. 6b). In one structure (OBOJOE), the methyl groups on the coordinating N atoms strain the ring, contributing to long Zn-N bonds, while in the other (MECVEV), bulk is added by the scaffold and functionalizing CN groups. In the MECVEV case, the addition of steric forces from axial CN groups might stabilize carbonic acid binding, disfavoring product release. In fact, we note that we calculate the product release in terms of the relative binding energy of a water molecule and a carbonic acid molecule, and analysis of the structure of MECVEV suggests steric clashing would be more important with the larger carbonic acid molecule than with the relatively smaller water molecule. In the other long-bond Zn-N case (OBOJOE), the methyl groups responsible for the longer Zn-N bonds possibly interfere sterically with an axially bound carbonic acid.
In both short-bonded CSD structures (QOVZEH and FEQYOP), the short bonds are achieved through a planar orientation of the ligands with accompanying scaffold features that favor short bonds. For instance, in both QOVZEH and FEQYOP, two five-membered rings are formed with the Zn center, two nitrogen ligands and two Csp3 atoms from the scaffold. In the case of FEQYOP, the other two rings formed are seven-membered, but this geometry appears to minimize strain on the short Zn-N bonds formed in satisfying the geometric constraints of the other two five-membered rings. Correlated with the observance of short Zn-N bonds, the QOVZEH and FEQYOP structures have longer axial Zn-OH interactions, lower activation energies and weakly less favorable product-release energetics. However, observation of the effect of the bulky axial groups in one long-bond Zn-N case (MECVEV) suggests that product release might be further facilitated through some combination of both short Zn-N bonds and axial steric effects, where possible. Future directions will include a focus on trends in relative binding energies when adding bulky axial groups to scaffold structures, either by identifying structures in the literature with bulky axial ligands or through identifying scaffolds that may be straightforwardly functionalized.
We have introduced and tested a computational approach for the design and optimization of catalysts. We have carried out an optimization scheme for carbon-capture catalysts based on the four key characteristics of the CO2 hydration reaction in carbonic anhydrase and synthetic mimics: (1) reactant hydroxyl formation, (2) CO2 hydration activation energy, (3) CO2 hydration reaction step energy, and (4) product release. Using both ZnII centers in combination with either four- or threefold coordination by Nsp3 or Nsp2 ligands, we observed that steps (2) and (4) are key for optimization. We identified the ZnII-(NH3)4 and ZnII-(NH2)4 metal/ligand sets as both starting from favorable reactivity at their equilibrium geometries and having reactivities that can be optimized through variation of structural properties. By constraining a previously characterized tetra-aza scaffold over a range of Zn-N bond lengths, we verified our initial observation that a short bond/shallow dihedral angle optimization direction improves the reaction characteristics around the CO2 insertion barrier. This approach of squeezing scaffolds to reduce activation energies, while simply a test of our minimal-ligand PES approach, could have significant application in embedding catalysts in compressive media (Bamford, 1962) or for the use of applied pressure to modify reaction kinetics (le Noble, 1967).
Finally, we searched the literature for outlying structures that were both consistent with and opposed to our geometric target. We found that the reaction properties, specifically the activation and product-release energies, for the short-bond Zn-N structures were consistent with or improved upon the tetra-aza scaffold reaction energetics. In the long-bond Zn-N structures, the activation energies were significantly higher, but product release appeared to depend also on the presence of axial scaffold components. Now that our approach has been validated for Nsp3 atoms, future study is also merited in searching for Nsp2 catalysts.
Overall, we have presented an approach for identifying potential catalysts. We employed existing literature on carbonic anhydrase and its mimics to identify potentially rate-limiting steps in our reaction and optimized putative reaction coordinates through constraining the coordination of metal centers. By seeking out geometric properties that optimize energy-based reaction characteristics, we may now use this information to build scaffolds combinatorially, combine existing literature scaffolds or design materials that constrain scaffolds away from equilibrium geometries to enhance their reactivity. This approach should easily be extended to have significant promise for catalyst design in other reactions, in which the number of steps may be significantly greater or more complex. In all cases, the relative binding energies of the species known to be present in the reaction provide an excellent starting point for characterizing relative reaction trends with geometric properties. Additionally, enumerating possible ways in which molecules may interconvert through transition states is computationally affordable due to our use of minimal models in the PES scan. Our results thus far have shown that reaction characteristics are relatively monotonic with varying geometry, justifying the use in the future of a sparser PES to characterize a greater number of metal and ligand combinations more rapidly. A key focus for our approach will be examining how we can employ geometric searches to identify how to control the relative energetics of different closely spaced spin surfaces, which are often present in few-site mid-row transition-metal catalysts.
This work was partly performed under the auspices of the US Department of Energy by the Lawrence Livermore National Laboratory under Contract No. DE-AC52-07NA27344. This project was funded by Laboratory Directed Research and Development grant No. 10-ERD-035. HJK holds a Career Award at the Scientific Interface from the Burroughs-Wellcome Fund. The use of computer resources at the Lawrence Livermore National Laboratory is gratefully acknowledged. This document has been reviewed and released with IM No. LLNL-JRNL-579437.
Adamczyk, K., Prémont-Schwarz, M., Pines, D., Pines, E. & Nibbering, E. T. J. (2009). Science, 326, 1690-1694.
Allen, F. H. (2002). Acta Cryst. B58, 380-388.
Baldwin, M. J. & Pecoraro, V. L. (1996). J. Am. Chem. Soc. 118, 11325-11326.
Bally, T. & Sastry, G. N. (1997). J. Phys. Chem. A, 101, 7923-7925.
Balzani, V., Campagna, S., Denti, G., Juris, A., Serroni, S. & Venturi, M. (1998). Acc. Chem. Res. 31, 26-34.
Bamford, C. R. (1962). Phys. Chem. Glasses, 3, 189-201.
Bao, L. & Trachtenberg, M. C. (2006). J. Membr. Sci. 280, 330-334.
Blaser, H. U., Spindler, F. & Studer, M. (2001). Appl. Catal. A, 221, 119-143.
Braïda, B., Hiberty, P. C. & Savin, A. (1998). J. Phys. Chem. A, 102, 7872-7877.
Bruno, I. J., Cole, J. C., Edgington, P. R., Kessler, M., Macrae, C. F., McCabe, P., Pearson, J. & Taylor, R. (2002). Acta Cryst. B58, 389-397.
Choi, K. Y., Suh, I. H. & Park, J. R. (1998). Main Group Met. Chem. 21, 783-788.
Cohen, A. J., Mori-Sánchez, P. & Yang, W. (2008). Science, 321, 792-794.
Coleman, J. E. (1967a). J. Biol. Chem. 242, 5212-5219.
Coleman, J. E. (1967b). Nature, 214, 193-194.
Dutoi, A. D. & Head-Gordon, M. (2006). Chem. Phys. Lett. 422, 230-233.
Exner, K. & von Ragué Schleyer, P. (2001). J. Phys. Chem. A, 105, 3407-3416.
Figueroa, J. D., Fout, T., Plasynski, S., McIlvried, H. & Srivastava, R. D. (2008). Int. J. Greenhouse Gas Control, 2, 9-20.
Giannozzi, P. et al. (2009). J. Phys. Condens. Matter, 21, 395502.
Gräfenstein, J., Kraka, E. & Cremer, D. (2004). J. Chem. Phys. 120, 524-539.
Greenwood, N. N. & Earnshaw, A. (1997). Chemistry of the Elements. 2nd revised ed, pp. 261-327. London: Butterworth-Heinemann Ltd.
Guidoni, L., Spiegel, K., Zumstein, M. & Röthlisberger, U. (2004). Angew. Chem. Int. Ed. 43, 3286-3289.
Håkansson, K., Carlsson, M., Svensson, L. A. & Liljas, A. (1992). J. Mol. Biol. 227, 1192-1204.
Handley, D. A., Hitchcock, P. B. & Leigh, G. J. (2001). Inorg. Chim. Acta, 314, 1-13.
Henkelman, G., Uberuaga, B. P. & Jónsson, H. (2000). J. Chem. Phys. 113, 9901-9904.
Khalifah, R. G. (1971). J. Biol. Chem. 246, 2561-2573.
Kong, D.-Y., Meng, L.-H., Ding, J., Xie, Y.-Y. & Huang, X.-Y. (2000). Polyhedron, 19, 217-223.
Koziol, L., Valdez, C. A., Baker, S. E., Lau, E. Y., Floyd, W. C., Wong, S. E., Satcher, J. H., Lightstone, F. C. & Aines, R. D. (2012). Inorg. Chem. 51, 6803-6812.
Krishnamurthy, V. M., Kaufman, G. K., Urbach, A. R., Gitlin, I., Gudiksen, K. L., Weibel, D. B. & Whitesides, G. M. (2008). Chem. Rev. 108, 946-1051.
Kulik, H. J., Cococcioni, M., Scherlis, D. A. & Marzari, N. (2006). Phys. Rev. Lett. 97, 103001.
Kulik, H. J. & Marzari, N. (2010). J. Chem. Phys. 133, 114103.
Lange, A. & Herbert, J. M. (2007). J. Chem. Theory Comput. 3, 1680-1690.
Liljas, A., Kannan, K. K., Bergstén, P. C., Waara, I., Fridborg, K., Strandberg, B., Carlbom, U., Järup, L., Lövgren, S. & Petef, M. (1972). Nature New Biol. 235, 131-137.
Lin, B.-Z., He, L.-W., Xu, B.-H., Li, X.-L., Li, Z. & Liu, P.-D. (2008). Cryst. Growth Des. 9, 273-281.
Lindskog, S. (1997). Pharmacol. Ther. 74, 1-20.
Lipscomb, W. N. (1983). Annu. Rev. Biochem. 52, 17-34.
Makov, G. & Payne, M. C. (1995). Phys. Rev. B, 51, 4014-4022.
Mori-Sánchez, P., Cohen, A. J. & Yang, W. T. (2006). J. Chem. Phys. 125, 201102.
Nørskov, J. K., Bligaard, T., Rossmeisl, J. & Christensen, C. H. (2009). Nature Chem. 1, 37-46.
Orgel, L. E. (1961). An Introduction to Transition Metal Chemistry: Ligand-Field Theory, pp. 1-184. New York: Wiley.
Pachauri, R. K. & Reisinger, A. (2007). Editors. Climate Change 2007 - Synthesis Report. Contribution of Working Groups I, II and III to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change. Geneva: IPCC.
Parkin, G. (2004). Chem. Rev. 104, 699-767.
Perdew, J. P., Burke, K. & Ernzerhof, M. (1996). Phys. Rev. Lett. 77, 3865-3868.
Petrie, S. & Stranger, R. (2004). Inorg. Chem. 43, 2597-2610.
Powell, S. R. (2000). J. Nutr. 130, 1447S-1454S.
Satcher, J. H. Jr, Baker, S. E., Kulik, H. J., Valdez, C. A., Krueger, R. L., Lightstone, F. C. & Aines, R. D. (2011). Energy Procedia, 4, 2090-2095.
Silverman, D. N. & Lindskog, S. (1988). Acc. Chem. Res. 21, 30-36.
Steiner, H., Jonsson, B. H. & Lindskog, S. (1975). Eur. J. Biochem. 59, 253-259.
Noble, W. J. le (1967). Progress in Physical Organic Chemistry, Vol. 5, edited by A. Streitwieser & R. W. Taft, Kinetics of Reactions in Solutions Under Pressure, pp. 207-330. New York: John Wiley and Sons Inc.
Togni, A. & Venanzi, L. M. (1994). Angew. Chem. Int. Ed. 33, 497-526.
Tozer, D. J. & De Proft, F. (2005). J. Phys. Chem. A, 109, 8923-8929.
Vallee, B. L. & Galdes, A. (1984). Adv. Enzymol. Relat. Areas Mol. Biol. 56, 283-430.
Wong, S. E., Lau, E. Y., Kulik, H. J., Satcher, J. H., Valdez, C., Worsely, M., Lightstone, F. C. & Aines, R. (2011). Energy Procedia, 4, 817-823.
Zhang, X. P. & van Eldik, R. (1995). Inorg. Chem. 34, 5606-5614.
Zhang, X. P., van Eldik, R., Koike, T. & Kimura, E. (1993). Inorg. Chem. 32, 5749-5755.
Zhang, Y. K. & Yang, W. T. (1998). J. Chem. Phys. 109, 2604-2608.