Toward Accurate QM/MM Reaction Barriers with Large QM Regions Using Domain Based Pair Natural Orbital Coupled Cluster Theory.

The hydroxylation reaction catalyzed by p-hydroxybenzoate hydroxylase and the Baeyer-Villiger reaction catalyzed by cyclohexanone monooxygenase are investigated by means of quantum mechanical/molecular mechanical (QM/MM) calculations at different levels of QM theory. The geometries of the stationary points along the reaction profile are obtained from QM/MM geometry optimizations, in which the QM region is treated by density functional theory (DFT). Relative energies are determined from single-point QM/MM calculations using the domain-based local pair natural orbital coupled cluster DLPNO-CCSD(T) method as QM component. The results are compared with single-point DFT/MM energies obtained using popular density functionals and with available experimental and computational data. It is found that the choice of the QM method strongly affects the computed energy profiles for these reactions. Different density functionals provide qualitatively different energy barriers (variations of the order of 10 kcal/mol in both reactions), thus limiting the confidence in DFT/MM computational predictions of energy profiles. On the other hand, the use of the DLPNO-CCSD(T) method in conjunction with large QM regions and basis sets makes it possible to achieve high accuracy. A critical discussion of all the technical aspects of the calculations is given with the aim of aiding computational chemists in the application of the DLPNO-CCSD(T) methodology in QM/MM calculations.

Multilevel Approaches within the Local Pair Natural Orbital Framework.

The linear-scaling local coupled cluster method DLPNO-CCSD(T) allows calculations on systems containing hundreds of atoms to be performed while reproducing canonical CCSD(T) energies typically with chemical accuracy (<1 kcal/mol). The accuracy of the method is determined by two main truncation thresholds that control the number of electron pairs included in the CCSD iterations and the size of the pair natural orbital virtual space for each electron pair, respectively. While the results of DLPNO-CCSD(T) calculations converge smoothly toward their canonical counterparts as the thresholds are tightened, the improved accuracy is accompanied by a fairly steep increase of the computational cost. Many applications study events that are confined to a relatively small region of the system of interest. Hence, it is viable to develop methods that allow the user to treat different parts of a large system at various levels of accuracy. In this work we present an extension to the native DLPNO method that fragments the system into preselected molecular parts and uses different thresholds or even different levels of theory for the interaction between individual fragments. Thereby chemical intuition can be used to focus computational resources on a more accurate evaluation of the properties at the center of interest, while permitting a less demanding description of the surrounding moieties. The strategy was implemented within the DLPNO-CCSD(T) framework. We tested the scheme for a series of realistic quantum chemical applications such as the calculation of the dimerization energies, potential energy surfaces, enantiomeric excess in organometallic catalysis, and the binding energy of the anticancer drug ellipticine to DNA. This work demonstrates the power of the approach and offers guidance to its setup.

Decomposition of Intermolecular Interaction Energies within the Local Pair Natural Orbital Coupled Cluster Framework.

The local coupled cluster method DLPNO-CCSD(T) allows calculations on systems containing hundreds of atoms to be performed while typically reproducing canonical CCSD(T) energies with chemical accuracy. In this work, we present a scheme for decomposing the DLPNO-CCSD(T) interaction energy between two molecules into physical meaningful contributions, providing a quantification of the most important components of the chemical interaction. The method, called Local Energy Decomposition (LED), is straightforward and requires negligible additional computing time. Both the Hartree-Fock and the correlation energy are decomposed into contributions from localized or pairs of localized occupied orbitals. Assigning these localized orbitals to fragments allows one to differentiate between intra- and intermolecular contributions to the interaction energy. Accordingly, the interaction energy can be decomposed into electronic promotion, electrostatic, exchange, dynamic charge polarization, and dispersion contributions. The LED scheme is applied to a number of test cases ranging from weakly, dispersively bound complexes to systems with strong ionic interactions. The dependence of the results on the one-particle basis set and various technical aspects, such as the localization scheme, are carefully studied in order to ensure that the results do not suffer from technical artifacts. A numerical comparison between the DLPNO-CCSD(T)/LED and the popular symmetry adapted perturbation theory (DFT-SAPT) is made, and the limitations of the proposed scheme are discussed.

Exploring the Accuracy Limits of Local Pair Natural Orbital Coupled-Cluster Theory.

The domain based local pair natural orbital coupled cluster method with single-, double-, and perturbative triple excitations (DLPNO­CCSD(T)) is an efficient quantum chemical method that allows for coupled cluster calculations on molecules with hundreds of atoms. Because coupled-cluster theory is the method of choice if high-accuracy is needed, DLPNO­CCSD(T) is very promising for large-scale chemical application. However, the various approximations that have to be introduced in order to reach near linear scaling also introduce limited deviations from the canonical results. In the present work, we investigate how far the accuracy of the DLPNO­CCSD(T) method can be pushed for chemical applications. We also address the question at which additional computational cost improvements, relative to the previously established default scheme, come. To answer these questions, a series of benchmark sets covering a broad range of quantum chemical applications including reaction energies, hydrogen bonds, and other noncovalent interactions, conformer energies, and a prototype organometallic problem were selected. An accuracy of 1 kcal/mol or better can readily be obtained for all data sets using the default truncation scheme, which corresponds to the stated goal of the original implementation. Tightening of the three thresholds that control DLPNO leads to mean absolute errors and standard deviations from the canonical results of less than 0.25 kcal/mol (<1 kJ/mol). The price one has then to pay is an increased computational time by a factor close to 3. The applicability of the method is shown to be independent of the nature of the reaction. On the basis of the careful analysis of the results, three different sets of truncation thresholds (termed "LoosePNO", "NormalPNO", and "TightPNO") have been chosen for "black box" use of DLPNO­CCSD(T). This will allow users of the method to optimally balance performance and accuracy.

Chemical applications carried out by local pair natural orbital based coupled-cluster methods.

The scope of this review is to provide a brief overview of the chemical applications carried out by local pair natural orbital coupled-electron pair and coupled-cluster methods. Benchmark tests reveal that these methods reproduce, with excellent accuracy, their canonical counterparts. At the same time, the speed up achieved by exploiting the locality of the electron correlation permits us to tackle chemical systems that, due to their size, would normally only be addressable with density functional theory. This review covers a broad variety of the chemical applications e.g. simulation of transition metal catalyzed reactions, estimation of weak interactions, and calculation of lattice properties in molecular crystals. This demonstrates that modern implementations of wavefunction-based correlated methods are playing an increasingly important role in applied computational chemistry.

Mechanism of Olefin Asymmetric Hydrogenation Catalyzed by Iridium Phosphino-Oxazoline: A Pair Natural Orbital Coupled Cluster Study.

Since the development of chiral phosphino-oxazoline iridium catalysts, which hydrogenate unfunctionalized alkenes enantioselectively, the asymmetric hydrogenation of prochiral olefins has become important in the production of chiral compounds. For the last 10 years, details of the mechanism, including formal oxidation state assignment of the metal center and the nature of intermediates and transition states have been debated. Various contributions have been given from a theoretical point of view, but due to the size of the structures, these have been forced to rely on density functional theory (DFT) methods. In our investigation of the catalytic cycle, we employ both DFT and a correlated ab initio method, namely, the newly implemented domain-based local pair natural orbital coupled-cluster theory with single and double excitations and the inclusion of perturbative triples correction (DLPNO-CCSD(T)). Our results show that the most likely active paths involve the formation of an intermediate Ir(V) species. Furthermore, we have been able to predict the absolute configuration of the major products, and where comparison to experiment is possible, the results of our calculations agree with the enantiomeric excess obtained from hydrogenating five prochiral substrates. This work also shows that it is now possible to study catalytic reactions with untruncated models (having up to 88 atoms) at the CCSD(T) level of theory.

Metal-dependent activity of Fe and Ni acireductone dioxygenases: how two electrons reroute the catalytic pathway.

Two virtually identical acireductone dioxygenases, ARD and ARD', catalyze completely different oxidation reactions of the same substrate, 1,2-dihydroxy-3-keto-5-(methylthio)pentene, depending exclusively on the nature of the bound metal. Fe(2+)-dependent ARD' produces the α-keto acid precursor of methionine and formate and allows for the recycling of methionine in cells. Ni(2+)-dependent ARD instead produces methylthiopropionate, CO, and formate, and exits the methionine salvage cycle. This mechanistic difference has not been understood to date but has been speculated to be due to the difference in coordination of the substrate to Fe(2+)versus Ni(2+): forming a five-membered ring versus a six-membered ring, respectively, thus exposing different carbon atoms for the attack by O2. Here, using mixed quantum-classical molecular dynamics simulations followed by the density functional theory mechanistic investigation, we show that, contrary to the old hypothesis, both metals preferentially bind the substrate as a six-membered ring, exposing the exact same sites to the attack by O2. It is the electronic properties of the metals that are then responsible for the system following different reaction paths, to yield the respective products. We fully explain the puzzling metal-induced difference in functionality between ARD and ARD' and, in particular, propose a new mechanism for ARD'. All results are in agreement with available isotopic substitution and other experimental data.

The Role of the Flexible L43-S54 Protein Loop in the CcrA Metallo-ß-lactamase in Binding Structurally Dissimilar ß-Lactam Antibiotics.

The CcrA di-Zn ß-lactamase is a bacterial enzyme capable of efficiently hydrolyzing and thus disabling a diverse set of ß-lactam antibiotics. Understanding the factors that contribute to the efficiency of CcrA is essential for the design of new CcrA-resistant antibiotics and enzyme inhibitors. The efficacy of CcrA has been speculated to be partially attributable to the flexible protein loop located above the active site (L43-S54), which would mold around structurally different substrates, for snag binding. Confirmation of this hypothesis about the role of the loop has been a challenge, from both an experimental and a theoretical point of view. We employed our newly developed method that combines extensive sampling of the protein structure via discrete molecular dynamics (DMD) and quantum mechanical (QM) treatment of the active site, QM/DMD, to investigate the structural role of the L43-S54 loop in binding three different ß-lactam antibiotics: imipenem, ampicillin, and cephalorodine. QM/DMD sampling was followed by high level ab initio calculations for the assessment of the energy contributions to loop-substrate interactions. We show that upon binding of all three antibiotic molecules, the loop comes in direct contact with the substrates and adopts distinctly different conformations depending on the bound substrate. The loop contributes to the binding affinity of CcrA to antibiotics. The primary component of the loop-substrate interaction is hydrophobic, and nonspecific, except for cephalorodine that is capable of π-stacking with W49 via one of the two competing modes.

How metal substitution affects the enzymatic activity of catechol-o-methyltransferase.

Catechol-O-methyltransferase (COMT) degrades catecholamines, such as dopamine and epinephrine, by methylating them in the presence of a divalent metal cation (usually Mg(II)), and S-adenosyl-L-methionine. The enzymatic activity of COMT is known to be vitally dependent on the nature of the bound metal: replacement of Mg(II) with Ca(II) leads to a complete deactivation of COMT; Fe(II) is slightly less than potent Mg(II), and Fe(III) is again an inhibitor. Considering the fairly modest role that the metal plays in the catalyzed reaction, this dependence is puzzling, and to date remains an enigma. Using a quantum mechanical / molecular mechanical dynamics method for extensive sampling of protein structure, and first principle quantum mechanical calculations for the subsequent mechanistic study, we explicate the effect of metal substitution on the rate determining step in the catalytic cycle of COMT, the methyl transfer. In full accord with experimental data, Mg(II) bound to COMT is the most potent of the studied cations and it is closely followed by Fe(II), whereas Fe(III) is unable to promote catalysis. In the case of Ca(II), a repacking of the protein binding site is observed, leading to a significant increase in the activation barrier and higher energy of reaction. Importantly, the origin of the effect of metal substitution is different for different metals: for Fe(III) it is the electronic effect, whereas in the case of Ca(II) it is instead the effect of suboptimal protein structure.

Hybrid dynamics simulation engine for metalloproteins.

Quality computational description of metalloproteins is a great challenge due to the vast span of time- and lengthscales characteristic of their existence. We present an efficient new method that allows for robust characterization of metalloproteins. It combines quantum mechanical (QM) description of the metal-containing active site, and extensive dynamics of the protein captured by discrete molecular dynamics (DMD) (QM/DMD). DMD samples the entire protein, including the backbone, and most of the active site, except for the immediate coordination region of the metal. QM operates on the part of the protein of electronic and chemical significance, which may include tens to hundreds of atoms. The breathing quantum-classical boundary provides a continuous mutual feedback between the two machineries. We test QM/DMD using the Fe-containing electron transporter protein, rubredoxin, and its three mutants as a model. QM/DMD can provide a reliable balanced description of metalloproteins' structure, dynamics, and electronic structure in a reasonable amount of time. As an illustration of QM/DMD capabilities, we then predict the structure of the Ca(2+) form of the enzyme catechol O-methyl transferase, which, unlike the native Mg(2+) form, is catalytically inactive. The Mg(2+) site is ochtahedral, but the Ca(2+) is 7-coordinate and features the misalignment of the reacting parts of the system. The change is facilitated by the backbone adjustment. QM/DMD is ideal and fast for providing this level of structural insight.

B13+ : a photodriven molecular Wankel engine.

Revved-up rotary: A molecular Wankel motor, the dual-ring structure B(13)(+), is driven by circularly-polarized infrared electromagnetic radiation. Calculations show that this illumination leads to a guided unidirectional rotation of the outer ring, which is achieved with rotational frequency of the order of 300 GHz.

Vibrational coupled cluster response theory: a general implementation.

The calculation of vibrational contributions to molecular properties using vibrational coupled cluster (VCC) response theory is discussed. General expressions are given for expectation values, linear response functions, and transition moments. It is shown how these expressions can be evaluated for arbitrary levels of excitation in the wave function parameterization as well as for arbitrary coupling levels in the potential and property surfaces. The convergence of the method is assessed by benchmark calculations on formaldehyde. Furthermore, excitation energies and infrared intensities are calculated for the fundamental vibrations of furan using VCC limited to up to two-mode and up to three-mode excitations, VCC[2] and VCC[3], as well as VCC with full two-mode and approximate three-mode couplings, VCC[2pt3].

New Formulation and Implementation of Vibrational Self-Consistent Field Theory.

A new implementation of the vibrational self-consistent field (VSCF) method is presented on the basis of a second quantization formulation. A so-called active terms algorithm is shown to be a significant improvement over a standard implementation reducing the computational effort by one order in the number of degrees of freedom. Various types of screening provide even further reductions in computational scaling and absolute CPU time. VSCF calculations on large polyaromatic hydrocarbon model systems are presented. Further, it is demonstrated that in cases where distant modes are not directly coupled in the Hamiltonian, down to linear scaling of the required CPU time with respect to the number of vibrational modes can be obtained. This is illustrated with calculations on simple model systems with up to 1 million degrees of freedom.

Using Electronic Energy Derivative Information in Automated Potential Energy Surface Construction for Vibrational Calculations.

The availability of an accurate representation of the potential energy surface (PES) is an essential prerequisite in an anharmonic vibrational calculation. At the same time, the high dimensionality of the fully coupled PES and the adverse scaling properties with respect to the molecular size make the construction of an accurate PES a computationally demanding task. In the past few years, our group tested and developed a series of tools and techniques aimed at defining computationally efficient, black-box protocols for the construction of PESs for use in vibrational calculations. This includes the definition of an adaptive density-guided approach (ADGA) for the construction of PESs from an automatically generated set of evaluation points. Another separate aspect has been the exploration of the use of derivative information through modified Shepard (MS) interpolation/extrapolation procedures. With this article, we present an assembled machinery where these methods are embedded in an efficient way to provide both a general machinery as well as concrete computational protocols. In this framework we introduce and discuss the accuracy and computational efficiency of two methods, called ADGA[2gx3M] and ADGA[2hx3M], where the ADGA recipe is used (with MS interpolation) to automatically define modest sized grids for up to two-mode couplings, while MS extrapolation based on, respectively, gradients only and gradients and Hessians from the ADGA determined points provides access to sufficiently accurate three-mode couplings. The performance of the resulting potentials is investigated in vibrational coupled cluster (VCC) calculations. Three molecular systems serve as benchmarks: a trisubstituted methane (CHFClBr), methanimine (CH2NH), and oxazole (C3H3NO). Furthermore, methanimine and oxazole are addressed in accurate calculations aiming to reproduce experimental results.

Potential energy surfaces for vibrational structure calculations from a multiresolution adaptive density-guided approach: implementation and test calculations.

A multiresolution procedure to construct potential energy surfaces (PESs) for use in vibrational structure calculations is developed in the framework of the adaptive density-guided approach. The implementation of the method allows the construction of hybrid PESs with different mode-coupling terms calculated with a variety of combinations of electronic structure methods and basis sets. Furthermore, the procedure allows the construction of hybrid PESs that incorporate a variety of contributions and corrections to the electronic energy, such as infinite basis set extrapolation and core correlation effects. A full account of the procedure is given together with a rather large set of benchmark calculations on a set of 20 small molecules, from diatomics to tetratomics.

Activity of rhodium-catalyzed hydroformylation: added insight and predictions from theory.

We have performed a density functional theory investigation of hydroformylation of ethylene for monosubstituted rhodium-carbonyl catalysts, HRh(CO)3L, where the modifying ligand, L, is a phosphite (L = P(OMe)3, P(OPh)3, or P(OCH2CF3)3), a phosphine (L = PMe3, PEt3, PiPr3, or PPh3), or a N-heterocyclic carbene (NHC) based on the tetrahydropyrimidine, imidazol, or tetrazol ring, respectively. The study follows the Heck and Breslow mechanism. Excellent correspondence between our calculations and existing experimental information is found, and the present results constitute the first example of a realistic quantum chemical description of the catalytic cycle of hydroformylation using ligand-modified rhodium carbonyl catalysts. This description explains the mechanistic and kinetic basis of the contemporary understanding of this class of reaction and offers unprecedented insight into the electronic and steric factors governing catalytic activity. The insight has been turned into structure-activity relationships and used as guidelines when also subjecting to calculation phosphite and NHC complexes that have yet to be reported experimentally. The latter calculations illustrate that it is possible to increase the electron-withdrawing capacity of both phosphite and NHC ligands compared to contemporary ligands through directed substitution. Rhodium complexes of such very electron-withdrawing ligands are predicted to be more active than contemporary catalysts for hydroformylation.

Structure and stability of networked metallofullerenes of the transition metals.

A DFT investigation of substitutionally doped fullerenes MC59 of second- and third-row transition metals shows that their stability increases toward the right-hand side of the d-block. Whereas the structural deviation from that of C60 depends on the size of the metal atom, stability is governed by electronic properties of the transition metal atom. A range of MC59 compounds of group 6-8 metals are predicted to have sufficient stability for experimental observation.