An Automated and Systematic Transition Structure Explorer in Large Flexible Molecular Systems Based on Combined Global Reaction Route Mapping and Microiteration Methods

Nanoscale chemical reaction exploration with a quantum magnifying glass

Nanoscopic systems exhibit diverse molecular substructures by which they facilitate specific functions. Theoretical models of them, which aim at describing, understanding, and predicting these capabilities, are difficult to build. Viable quantum-classical hybrid models come with specific challenges regarding atomistic structure construction and quantum region selection. Moreover, if their dynamics are mapped onto a state-to-state mechanism such as a chemical reaction network, its exhaustive exploration will be impossible due to the combinatorial explosion of the reaction space. Here, we introduce a "quantum magnifying glass" that allows one to interactively manipulate nanoscale structures at the quantum level. The quantum magnifying glass seamlessly combines autonomous model parametrization, ultra-fast quantum mechanical calculations, and automated reaction exploration. It represents an approach to investigate complex reaction sequences in a physically consistent manner with unprecedented effortlessness in real time. We demonstrate these features for reactions in bio-macromolecules and metal-organic frameworks, diverse systems that highlight general applicability.

Variations on the Bergman Cyclization Theme: Electrocyclizations of Ionic Penta-, Hepta-, and Octadiynes

The Bergman cyclization of (Z)-hexa-3-ene-1,5-diyne to form the aromatic diradical p-benzyne has garnered attention as a potential antitumor agent due to its relatively low cyclization barrier and the stability of the resulting diradical. Here, we present a theoretical investigation of several ionic extensions of the fundamental Bergman cyclization: electrocyclizations of the penta-1,4-diyne anion, hepta-1,6-diyne cation, and octa-1,7-diyne dication, leveraging the spin-flip formulation of the equation-of-motion coupled cluster theory with single and double substitutions (EOM-SF-CCSD). Though the penta-1,4-diyne anion exhibits a large cyclization barrier of +66 kcal mol-1, cyclization of both the hepta-1,6-diyne cation and octa-1,7-diyne dication along a previously unreported triplet pathway requires relatively low energy. We also identified the presence of significant aromaticity in the triplet diradical products of these two cationic cyclizations.

Calculation of the permeability coefficients of small molecules through lipid bilayers by free-energy reaction network analysis following the explicit treatment of the internal conformation of the solute

Biomembrane permeation represents a major barrier to pharmacokinetics. During preclinical drug discovery, the coefficients of the permeation of molecules through lipid bilayers account for a valuable property of such molecules. Therefore, the control of the permeation of molecules through lipid bilayers is an essential factor in drug design, and the estimation of the permeation phenomena is a crucial study in pharmacy. Thus, there are many published studies on the theoretical estimations of permeation coefficients. Here, we propose a molecular dynamics (MD) simulation method for estimating the permeation of small molecules through lipid bilayers based on the free-energy reaction network (FERN) analysis. In this method, the collective variables (CVs) of the free energy calculations explicitly include the conformational changes in the rotational bonds of the solute molecules. The advantages of the present method over the other method are that it is possible to estimate reaction pathways and their reaction rates, i.e., permeation coefficients or passage times, in multidimensional space spanned by CVs though conventional methods such as the umbrella sampling method and target MDs often dealt with a few degrees of freedom. To demonstrate the efficacy of our method, we calculate the coefficients of the permeation of three small aromatic peptides, namely N-acetylphenylalanineamide (Ac-Phe-NH2 or NAFA), N-acetyltyrosineamide (Ac-Tyr-NH2 or NAYA), and N-acetyltryptophanamide (Ac-Trp-NH2 or NATA), through a 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC) lipid bilayer. In these cases we adopted one CV for the permeation direction and four CVs for the internal rotational coordinates. The results reveal that the permeation coefficients of NAFA, NAYA, and NATA are 1.7 × 10^-2, 0.51 × 10^-4, and 5.7 × 10^-4 cm s^-1, respectively. Compared with the experimental data, our simulation results followed the same trend, i.e., NAFA > NATA > NAYA. By analyzing the structures of metastable points of the solute molecules, our simulation result reveals that the aforementioned trend is caused by the differences in stability among their rotamers. Furthermore, we evaluate the statistical fluctuation of the rotamers, and the time scale of flipping the side chain reveals that the structures rigidify as the ligand moves deeper into the membrane.

Automated input structure generation for single-ended reaction path optimizations

The exploration of a reaction network requires highly automated workflows to avoid error-prone and time-consuming manual steps. In this respect, a major bottleneck is the search for transition-state (TS) structures, which frequently fails and, therefore, makes (manual) revision necessary. In this work, we present a technique for obtaining suitable input structures for automated TS searches based on single-ended reaction path optimization algorithms, which makes subsequent TS searches via this method significantly more robust. First, possible input structures are generated based on the spatial alignment of the reactants. The appropriate orientation of reacting groups is achieved via stepwise rotations along selected torsional degrees of freedom. Second, a ranking of the obtained structures is performed according to selected geometric criteria. The main goals are to properly align the reactive atoms, to avoid hindrance within the reaction channel and to resolve steric clashes between the reactants. The developed procedure has been carefully tested on a variety of examples and provides suitable input structures for TS searches within seconds. The method is in daily use in an industrial setting.

Multistructural microiteration combined with QM/MM-ONIOM electrostatic embedding

Multistructural microiteration (MSM) is a method to take account of contributions of multiple surrounding structures in a geometrical optimization or reaction path calculation using the quantum mechanics/molecular mechanics (QM/MM) ONIOM method. In this study, we combined MSM with the electrostatic embedding (EE) scheme of the QM/MM-ONIOM method by extending its original formulation for mechanical embedding (ME). MSM-EE takes account of the polarization in the QM region induced by point charges assigned to atoms in the multiple surrounding structures, where the point charges are scaled by the weight factor of each surrounding structure determined through MSM. The performance of MSM-EE was compared with that of the other methods, i.e., ONIOM-ME, ONIOM-EE, and MSM-ME, by applying them to three chemical processes: (1) chorismate-to-prephenate transformation in aqueous solution, (2) the same transformation as (1) in an enzyme, and (3) hydroxylation in p-hydroxybenzoate hydroxylase. These numerical tests of MSM-EE yielded barriers and reaction energies close to experimental values with computational costs comparable to those of the other three methods.

DFT and AFIR study on the copper(i)-catalyzed mechanism of 5-enamine-trisubstituted-1,2,3-triazole synthesis via C-N cross-coupling and the origin of ring-opening of 2H-azirines

Understanding the synthesis mechanism of substituted 1,2,3-triazoles is an important and state-of-the-art research area of contemporary copper(i)-catalyzed terminal alkyne and organic azide click reaction (CuAAC), which has invoked increasing close collaborations between experiment and theory including copper catalyzed interrupted click reaction. In this study, the mechanism of Cu(i)-catalyzed 5-enamine-functionalized fully substituted 1,2,3-triazole synthesis was rationalized via density functional theory (DFT) and multicomponent artificial force-induced reaction (MC-AFIR) methods. The reasonable reaction route consists of (a) di-copper catalyzed ring-opening of 2H-azirines, (b) alkyne hydrogen atom transfer, (c) [3 + 2] ring cycloaddition, and (d) C-N bond formation through reductive elimination. The MC-AFIR method was used for the systematic determination of transition states for the C/N-Cu bond formation, C-N bond coupling and crossing points between singlet and triplet states. Our survey on the prereactant complexes suggested that the dicopper-catalyzed 2H-azirine ring-opening and alkyne hydrogen activation are both thermodynamically feasible via a singlet/triplet crossing point. This explains why Et3N is critical for alkyne hydrogen transfer (HT) before the [3 + 2] cycloaddition reaction, and the C-N cross-coupling product instead of the click product (byproduct). Our DFT results indicate that the transmetalation process is the rate determination step along the triplet state potential energy surface. This study provides important mechanistic insights for the interrupted CuAAC reaction to form 5-enamine-fully-substituted-1,2,3-triazoles. Further insight prediction interprets that solvent and extra strong ligand coordination play a certain role in competitive reactions.

Trialkylphosphines Having a Bulky Phosphacyclopentane Backbone: Structural and Redox Properties Depending on the Exocyclic Alkyl Groups and EPR Observation of a Persistent Trialkylphosphine Radical Cation

Bulky phosphines and their redox properties have received increased attention in the view of useful auxiliary ligands for transition metal catalysts and Lewis-base components of frustrated Lewis pairs for chemical transformations. Herein we report the synthesis, structure, and properties of a series of trialkylphosphines 2R (R = methyl, ethyl, isopropyl, tert-butyl, 1-adamantyl) that possess the bulky 2,2,5,5-tetrakis(trimethylsilyl)-1-phosphacyclopentane as a structural backbone. Among these phosphines, 2Ad, which contains an adamantyl moiety, has a very large buried volume (%V_bur) for a trialkylphosphine (62.0) and shows a quasi-reversible oxidative wave at a lower oxidation potential (-0.12 V in CH₂Cl₂, vs ferrocene/ferrocenium couple) by cyclic voltammetry. The reaction of 2Ad with AgPF₆ afforded a cationic silver aquo complex [Ag(2Ad)(H₂O)]⁺[PF₆]^-, whereas the reaction with NOSbF₆ gave a persistent phosphine radical cation [2Ad]^•+. Based on the EPR spectra and DFT studies, the spin and positive charge of [2Ad]^•+ are localized on the phosphorus atom.

Analytical Method Using a Scaled Hypersphere Search for High-Dimensional Metadynamics Simulations

Metadynamics (MTD) is one of the most effective methods for calculating the free energy surface and finding rare events. Nevertheless, numerous studies using MTD have been carried out using 3D or lower dimensional collective variables (CVs), as higher dimensional CVs require costly computational resources and the obtained results are too complex to understand the important events. The latter issue can be conveniently solved by utilizing the free energy reaction network (FERN), which is a graph structure consisting of edges of minimum free energy paths (MFEPs), nodes of equation (EQ) points, and transition state (TS) points. In the present article, a new method for exploring FERNs on high-dimensional CVs using MTD and the scaled hypersphere search (SHS) method is described. A test calculation based on the MTD-SHS simulation of met-enkephalin in explicit water with 7 CVs was conducted. As a result, 889 EQ points and 1805 TS points were found. The MTD-SHS approach can find MFEPs exhaustively; therefore, the FERNs can be estimated without any a priori knowledge of the EQ and TS points.

Analysis of degradation products of nitrogen mustards via hydrophilic interaction liquid chromatography-tandem mass spectrometry

A hydrophilic interaction liquid chromatography-tandem mass spectrometry (HILIC-MS/MS) method was developed for qualitative and quantitative analysis of ethanolamines (EAs), which are nitrogen mustard degradation products. With this method, the retention times of the highly hydrophilic EAs on the HILIC column were sufficient (retention times: methyl diethanolamine, 12.2 min; ethyl diethanolamine, 11.2 min; and triethanolamine, 9.5 min) and the EAs were analyzed more efficiently than with reported HILIC-MS/MS methods. The detection limits of methyl diethanolamine and ethyl diethanolamine in serum and urine using this approach were 15-20 ng/mL. The suitability of the method for real samples was evaluated via recovery tests involving urine and serum, and the method was validated. The MS/MS fragmentation of EAs was discussed based on density functional theory calculations.

Dearomative cycloadditions of a silylene with pyrazine and quinoxaline

The dearomative cycloaddition of a dialkylsilylene with pyrazine affords an eight-membered heterocyclic compound with a trans-cycloalkene moiety.

Implementation and performance of the artificial force induced reaction method in the GRRM17 program

This article reports implementation and performance of the artificial force induced reaction (AFIR) method in the upcoming 2017 version of GRRM program (GRRM17). The AFIR method, which is one of automated reaction path search methods, induces geometrical deformations in a system by pushing or pulling fragments defined in the system by an artificial force. In GRRM17, three different algorithms, that is, multicomponent algorithm (MC-AFIR), single-component algorithm (SC-AFIR), and double-sphere algorithm (DS-AFIR), are available, where the MC-AFIR was the only algorithm which has been available in the previous 2014 version. The MC-AFIR does automated sampling of reaction pathways between two or more reactant molecules. The SC-AFIR performs automated generation of global or semiglobal reaction path network. The DS-AFIR finds a single path between given two structures. Exploration of minimum energy structures within the hypersurface in which two different electronic states degenerate, and an interface with the quantum mechanics/molecular mechanics method, are also described. A code termed SAFIRE will also be available, as a visualization software for complicated reaction path networks. © 2017 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.

Automated Transition State Theory Calculations for High-Throughput Kinetics

A scarcity of known chemical kinetic parameters leads to the use of many reaction rate estimates, which are not always sufficiently accurate, in the construction of detailed kinetic models. To reduce the reliance on these estimates and improve the accuracy of predictive kinetic models, we have developed a high-throughput, fully automated, reaction rate calculation method, AutoTST. The algorithm integrates automated saddle-point geometry search methods and a canonical transition state theory kinetics calculator. The automatically calculated reaction rates compare favorably to existing estimated rates. Comparison against high level theoretical calculations show the new automated method performs better than rate estimates when the estimate is made by a poor analogy. The method will improve by accounting for internal rotor contributions and by improving methods to determine molecular symmetry.

Multistructural microiteration technique for geometry optimization and reaction path calculation in large systems

We propose a multistructural microiteration (MSM) method for geometry optimization and reaction path calculation in large systems. MSM is a simple extension of the geometrical microiteration technique. In conventional microiteration, the structure of the non-reaction-center (surrounding) part is optimized by fixing atoms in the reaction-center part before displacements of the reaction-center atoms. In this method, the surrounding part is described as the weighted sum of multiple surrounding structures that are independently optimized. Then, geometric displacements of the reaction-center atoms are performed in the mean field generated by the weighted sum of the surrounding parts. MSM was combined with the QM/MM-ONIOM method and applied to chemical reactions in aqueous solution or enzyme. In all three cases, MSM gave lower reaction energy profiles than the QM/MM-ONIOM-microiteration method over the entire reaction paths with comparable computational costs. © 2017 Wiley Periodicals, Inc.

Reliable and efficient reaction path and transition state finding for surface reactions with the growing string method

The computational challenge of fast and reliable transition state and reaction path optimization requires new methodological strategies to maintain low cost, high accuracy, and systematic searching capabilities. The growing string method using internal coordinates has proven to be highly effective for the study of molecular, gas phase reactions, but difficulties in choosing a suitable coordinate system for periodic systems has prevented its use for surface chemistry. New developments are therefore needed, and presented herein, to handle surface reactions which include atoms with large coordination numbers that cannot be treated using standard internal coordinates. The double-ended and single-ended growing string methods are implemented using a hybrid coordinate system, then benchmarked for a test set of 43 elementary reactions occurring on surfaces. These results show that the growing string method is at least 45% faster than the widely used climbing image-nudged elastic band method, which also fails to converge in several of the test cases. Additionally, the surface growing string method has a unique single-ended search method which can move outward from an initial structure to find the intermediates, transition states, and reaction paths simultaneously. This powerful explorative feature of single ended-growing string method is demonstrated to uncover, for the first time, the mechanism for atomic layer deposition of TiN on Cu(111) surface. This reaction is found to proceed through multiple hydrogen-transfer and ligand-exchange events, while formation of H-bonds stabilizes intermediates of the reaction. Purging gaseous products out of the reaction environment is the driving force for these reactions. © 2017 Wiley Periodicals, Inc.

Study of Potential Energy Surfaces towards Global Reaction Route Mapping

The potential energy surface (PES) is just a theoretical construct based on the Born-Oppenheimer approximation, but it underlies various phenomena, including molecular vibrations, collisional ionizations, and chemical reactions. This account describes how a new idea for global reaction route mapping (GRRM), which had seemed to be impossible for chemical systems with more than three atoms, was born and has been developed during the course of the study of the PES. GRRM has pioneered new fields of chemistry. Furthermore, techniques for GRRM are still developing, and GRRM is further extending its application to various areas of chemistry and chemical physics.

Methods for Efficiently and Accurately Computing Quantum Mechanical Free Energies for Enzyme Catalysis

Enzyme activity is inherently linked to free energies of transition states, ligand binding, protonation/deprotonation, etc.; these free energies, and thus enzyme function, can be affected by residue mutations, allosterically induced conformational changes, and much more. Therefore, being able to predict free energies associated with enzymatic processes is critical to understanding and predicting their function. Free energy simulation (FES) has historically been a computational challenge as it requires both the accurate description of inter- and intramolecular interactions and adequate sampling of all relevant conformational degrees of freedom. The hybrid quantum mechanical molecular mechanical (QM/MM) framework is the current tool of choice when accurate computations of macromolecular systems are essential. Unfortunately, robust and efficient approaches that employ the high levels of computational theory needed to accurately describe many reactive processes (ie, ab initio, DFT), while also including explicit solvation effects and accounting for extensive conformational sampling are essentially nonexistent. In this chapter, we will give a brief overview of two recently developed methods that mitigate several major challenges associated with QM/MM FES: the QM non-Boltzmann Bennett's acceptance ratio method and the QM nonequilibrium work method. We will also describe usage of these methods to calculate free energies associated with (1) relative properties and (2) along reaction paths, using simple test cases with relevance to enzymes examples.

Automated Transition State Searches without Evaluating the Hessian

Accurate and speedy determination of transition structures (TSs) is essential for computational studies on reaction pathways, particularly when the process involves expensive electronic structure calculations. Many search algorithms require a good initial guess of the TS geometry, as well as a Hessian input that possesses a structure consistent with the desired saddle point. Among the double-ended interpolation methods for generation of the guess for the TS, the freezing string method (FSM) is proven to be far less expensive compared to its predecessor, the growing string method (GSM). In this paper, it is demonstrated that the efficiency of this technique can be improved further by replacing the conjugate gradient optimization step (FSM-CG) with a quasi-Newton line search coupled with a BFGS Hessian update (FSM-BFGS). A second crucial factor that affects the speed with which convergence to the TS is achieved is the quality and cost of the Hessian of the energy for the guessed TS. For electronic structure calculations, the cost of calculating an exact Hessian increases more rapidly with system size than the energy and gradient. Therefore, to sidestep calculation of the exact Hessian, an approximate Hessian is constructed, using the tangent direction and local curvature at the TS guess. It is demonstrated that the partitioned-rational function optimization algorithm for locating TSs with this approximate Hessian input performs at least as well as with an exact Hessian input in most test cases. The two techniques, FSM and approximate Hessian construction, therefore can significantly reduce costs associated with finding TSs.

Exploring Multiple Potential Energy Surfaces: Photochemistry of Small Carbonyl Compounds

In theoretical studies of chemical reactions involving multiple potential energy surfaces (PESs) such as photochemical reactions, seams of intersection among the PESs often complicate the analysis. In this paper, we review our recipe for exploring multiple PESs by using an automated reaction path search method which has previously been applied to single PESs. Although any such methods for single PESs can be employed in the recipe, the global reaction route mapping (GRRM) method was employed in this study. By combining GRRM with the proposed recipe, all critical regions, that is, transition states, conical intersections, intersection seams, and local minima, associated with multiple PESs, can be explored automatically. As illustrative examples, applications to photochemistry of formaldehyde and acetone are described. In these examples as well as in recent applications to other systems, the present approach led to discovery of many unexpected nonadiabatic pathways, by which some complicated experimental data have been explained very clearly.

On the perturbation of the H-bonding interaction in ethylene glycol clusters upon hydration

Ab initio and density functional methods have been employed to study the structure, stability, and spectral properties of various ethylene glycol (EG(m)) and ethylene glycol-water (EG(m)W(n)) (m = 1-3, n = 1-4) clusters. The effective fragment potential (EFP) approach was used to explore various possible EG(m)W(n) clusters. Calculated interaction energies of EG(m)W(n) clusters confirm that the hydrogen-bonding interaction between EG molecules is perturbed by the presence of water molecules and vice versa. Further, energy decomposition analysis shows that both electrostatic and polarization interactions predominantly contribute to the stability of these clusters. It was found from the same analysis that ethylene glycol-water interaction is predominant over the ethylene glycol-ethylene glycol and water-water interactions. Overall, the results clearly illustrate that the presence of water disrupts the ethylene glycol-ethylene glycol hydrogen bonds.

Finding Minimum Structures on the Seam of Crossing in Reactions of Type A + B → X: Exploration of Nonadiabatic Ignition Pathways of Unsaturated Hydrocarbons

A new theoretical approach is proposed for finding automatically minimum structures on the seam of crossing (MSX) in reactions of type A + B → X, where the artificial-force-induced reaction (AFIR) method is combined with the seam model function (SMF) approach. Its application to reactions between triplet dioxygen and unsaturated hydrocarbons provided many MSX structures. In addition to known ignition pathways, we discovered a pathway through a new type of MSX in the reaction of dioxygen with aromatic hydrocarbons; for benzene, this new pathway requires a lower energy than those of three known ignition pathways and is likely to be the most important. This demonstrates that the AFIR-SMF approach has the ability to discover unknown/unexpected MSX structures without prejudice for presumed pathways or mechanisms.