Normal Mode Sampling Research Articles

Transition state search and geometry relaxation throughout chemical compound space with quantum machine learning.

We use energies and forces predicted within response operator based quantum machine learning (OQML) to perform geometry optimization and transition state search calculations with legacy optimizers but without the need for subsequent re-optimization with quantum chemistry methods. For randomly sampled initial coordinates of small organic query molecules, we report systematic improvement of equilibrium and transition state geometry output as training set sizes increase. Out-of-sample SN2 reactant complexes and transition state geometries have been predicted using the LBFGS and the QST2 algorithms with an root-mean-square deviation (RMSD) of 0.16 and 0.4 Å-after training on up to 200 reactant complex relaxations and transition state search trajectories from the QMrxn20 dataset, respectively. For geometry optimizations, we have also considered relaxation paths up to 5'595 constitutional isomers with sum formula C7H10O2 from the QM9-database. Using the resulting OQML models with an LBFGS optimizer reproduces the minimum geometry with an RMSD of 0.14 Å, only using ∼6000 training points obtained from normal mode sampling along the optimization paths of the training compounds without the need for active learning. For converged equilibrium and transition state geometries, subsequent vibrational normal mode frequency analysis indicates deviation from MP2 reference results by on average 14 and 26 cm-1, respectively. While the numerical cost for OQML predictions is negligible in comparison to density functional theory or MP2, the number of steps until convergence is typically larger in either case. The success rate for reaching convergence, however, improves systematically with training set size, underscoring OQML's potential for universal applicability.

Computational Spectroscopy of the Cr-Cr Bond in Coordination Complexes.

We report the accurate computational vibrational analysis of the Cr-Cr bond in dichromium complexes using second-order multireference complete active space methods (CASPT2), allowing direct comparison with experimental spectroscopic data both to facilitate interpreting the low-energy region of the spectra and to provide insights into the nature of the bonds themselves. Recent technological development by the authors has realized such computation for the first time. Accurate simulation of the vibrational structure of these compounds has been hampered by their notorious multiconfigurational electronic structure that yields bond distances that do not correlate with bond order. Some measured Cr-Cr vibrational stretching modes, ν(Cr2), have suggested weaker bonding, even for so-called ultrashort Cr-Cr bonds, while others are in line with the bond distance. Here, we optimize geometries and compute ν(Cr2) with CASPT2 for three well-characterized complexes, Cr2(O2CCH3)4(H2O)2, Cr2(mhp)4, and Cr2(dmp)4. We obtain CASPT2 harmonic ν(Cr2) modes in good agreement with experiment at 282 cm-1 for Cr2(mhp)4 and 353 cm-1 for Cr2(dmp)4, compute 50Cr and 54Cr isotope shifts, and demonstrate that the use of the so-called IPEA shift leads to improved Cr-Cr distances. Additionally, normal mode sampling was used to estimate anharmonicity along ν(Cr2), leading to an anharmonic mode of 272 cm-1 for Cr2(mhp)4 and 333 cm-1 for Cr2(dmp)4.

Beyond the normal mode picture: the importance of the reactant’s intramolecular mode coupling in quasiclassical trajectory simulations

In quasiclassical trajectory simulations of bimolecular reactions of polyatomic molecules, the sets of initial coordinates and momenta of the colliding molecules generated by the widely used normal mode sampling (NMS) are generally nonstationary and evolve during the initial free flight of the reactants. For several isotopologues of the CH4 + H abstraction reaction, the calculated cross sections were found to oscillate as a function of the initial reactant separation. The arising ambiguity can be removed by averaging over an entire period of such oscillations (1PA). Adiabatic switching (AS) produces stationary ensembles of initial states, allowing one to generate unambiguous reaction probabilities and cross sections regardless of the choice of coordinate system used in the zeroth-order harmonic Hamiltonian. In this work the accuracy of the 1PA method was tested against AS-based reactivity parameters. Overall, 1PA provides a reasonably good correction to the non-stationarity of the initial states generated by NMS, but in extreme cases the cross sections are overestimated by up to 70%, thus the accuracy of 1PA varies currently unpredictably and is often sensitive to the choice of coordinates used for NMS. Therefore, it seems advisable to use adiabatic switching to prepare initial states in quasiclassical trajectory simulations.

Theoretical kinetic analysis and molecular dynamic simulations for exploring the dissociation pathways of CF3OH

The unimolecular dissociation of trifluoromethanol (CF3OH) had been studied by static and dynamic calculations. The electronic structure calculations and molecular dynamic (MD) simulations were performed at the second order Moller-Plesset perturbation theory (MP2) theory. The kinetic behavior was investigated in the temperature range from 298 to 1400 K and pressure range between 10−9 to 104 kPa using Troe’s formalism. In addition to this, Born-Oppenheimer molecular dynamic (BOMD) model was employed for exploring the dynamic insights of the reaction. Reaction trajectories were estimated by microcanonical normal mode sampling, wherein selected vibrational modes were pumped with additional energy of 95, 155, and 208 kcal mol−1. Our calculations suggest that, the unimolecular decomposition of CF3OH results in the generation of molecular HF with concomitant generation of CF2O via intramolecular reaction. The MD simulations suggest the critical role of vibrational modes for such intramolecular reaction.

Metadynamics for training neural network model chemistries: A competitive assessment.

Neural network model chemistries (NNMCs) promise to facilitate the accurate exploration of chemical space and simulation of large reactive systems. One important path to improving these models is to add layers of physical detail, especially long-range forces. At short range, however, these models are data driven and data limited. Little is systematically known about how data should be sampled, and "test data" chosen randomly from some sampling techniques can provide poor information about generality. If the sampling method is narrow, "test error" can appear encouragingly tiny while the model fails catastrophically elsewhere. In this manuscript, we competitively evaluate two common sampling methods: molecular dynamics (MD), normal-mode sampling, and one uncommon alternative, Metadynamics (MetaMD), for preparing training geometries. We show that MD is an inefficient sampling method in the sense that additional samples do not improve generality. We also show that MetaMD is easily implemented in any NNMC software package with cost that scales linearly with the number of atoms in a sample molecule. MetaMD is a black-box way to ensure samples always reach out to new regions of chemical space, while remaining relevant to chemistry near kbT. It is a cheap tool to address the issue of generalization.

Adiabatic Switching Extended To Prepare Semiclassically Quantized Rotational-Vibrational Initial States for Quasiclassical Trajectory Calculations.

An approximation-free adiabatic switching method to generate semiclassically quantized ensembles of rovibrational states of polyatomic molecules for use as initial conditions in quasiclassical trajectory calculations is presented. Vibrational states are prepared, starting from an ensemble of classical states corresponding to the desired quantum state of the normal-mode Hamiltonian by slowly switching on the anharmonicity in internal coordinates, thereby avoiding rotational contamination. To generate rovibrational states, an extension is proposed: The vibrationally quantized molecules are slowly spun up to the desired quantized angular momentum. The ensembles obtained with adiabatic switching for CH4 are insensitive to the choice of internal coordinates and stationary; furthermore, their mean energies agree remarkably well with the quantum mechanical values: The zero-point energy and 15 vibrational levels of the first three polyads are within 20 cm-1, the rotational levels are between J = 1 and 50 within 1%, and the standard deviation is always <1%. Adiabatic switching produces classical state ensembles with significantly better properties than normal-mode sampling, making them more appropriate in quasiclassical trajectory calculations.

ANI-1: an extensible neural network potential with DFT accuracy at force field computational cost.

Deep learning is revolutionizing many areas of science and technology, especially image, text, and speech recognition. In this paper, we demonstrate how a deep neural network (NN) trained on quantum mechanical (QM) DFT calculations can learn an accurate and transferable potential for organic molecules. We introduce ANAKIN-ME (Accurate NeurAl networK engINe for Molecular Energies) or ANI for short. ANI is a new method designed with the intent of developing transferable neural network potentials that utilize a highly-modified version of the Behler and Parrinello symmetry functions to build single-atom atomic environment vectors (AEV) as a molecular representation. AEVs provide the ability to train neural networks to data that spans both configurational and conformational space, a feat not previously accomplished on this scale. We utilized ANI to build a potential called ANI-1, which was trained on a subset of the GDB databases with up to 8 heavy atoms in order to predict total energies for organic molecules containing four atom types: H, C, N, and O. To obtain an accelerated but physically relevant sampling of molecular potential surfaces, we also proposed a Normal Mode Sampling (NMS) method for generating molecular conformations. Through a series of case studies, we show that ANI-1 is chemically accurate compared to reference DFT calculations on much larger molecular systems (up to 54 atoms) than those included in the training data set.

Oscillatory reaction cross sections caused by normal mode sampling in quasiclassical trajectory calculations.

The quasiclassical trajectory (QCT) method is an efficient and important tool for studying the dynamics of bimolecular reactions. In this method, the motion of the atoms is simulated classically, and the only quantum effect considered is that the initial vibrational states of reactant molecules are semiclassically quantized. A sensible expectation is that the initial ensemble of classical molecular states generated this way should be stationary, similarly to the quantum state it is supposed to represent. The most widely used method for sampling the vibrational phase space of polyatomic molecules is based on the normal mode approximation. In the present work, it is demonstrated that normal mode sampling provides a nonstationary ensemble even for a simple molecule like methane, because real potential energy surfaces are anharmonic in the reactant domain. The consequences were investigated for reaction CH4 + H → CH3 + H2 and its various isotopologs and were found to be dramatic. Reaction probabilities and cross sections obtained from QCT calculations oscillate periodically as a function of the initial distance of the colliding partners and the excitation functions are erratic. The reason is that in the nonstationary ensemble of initial states, the mean bond length of the breaking C-H bond oscillates in time with the frequency of the symmetric stretch mode. We propose a simple method, one-period averaging, in which reactivity parameters are calculated by averaging over an entire period of the mean C-H bond length oscillation, which removes the observed artifacts and provides the physically most reasonable reaction probabilities and cross sections when the initial conditions for QCT calculations are generated by normal mode sampling.

QM/MM Protocol for Direct Molecular Dynamics of Chemical Reactions in Solution: The Water-Accelerated Diels-Alder Reaction.

We describe a solvent-perturbed transition state (SPTS) sampling scheme for simulating chemical reaction dynamics in condensed phase. The method, adapted from Truhlar and Gao's ensemble-averaged variational transition state theory, includes the effect of instantaneous solvent configuration on the potential energy surface of the reacting system (RS) and allows initial conditions for the RS to be sampled quasiclassically by TS normal mode sampling. We use a QM/MM model with direct dynamics, in which QM forces of the RS are computed at each trajectory point. The SPTS scheme is applied to the acceleration of the Diels-Alder reaction of cyclopentadiene (CP) + methyl vinyl ketone (MVK) in water. We explored the effect of the number of SPTS and of solvent box size on the distribution of bond lengths in the TS. Statistical sampling of the sampling was achieved when distribution of forming bond lengths converged. We describe the region enclosing the partial bond lengths as the transition zone. Transition zones in the gas phase, SMD implicit solvent, QM/MM, and QM/MM+QM (3 water molecules treated by QM) vary according to the ability of the medium to stabilize zwitterionic structures. Mean time gaps between formation of C-C bonds vary from 11 fs for gas phase to 25 fs for QM/MM+QM. Mean H-bond lengths to O(carbonyl) in QM/MM+QM are 0.14 Å smaller at the TS than in MVK reactant, and the mean O(carbonyl)-H(water)-O(water) angle of H-bonds at the TS is 10° larger than in MVK reactant.

Ab initio classical trajectory calculations of 1,3-cyclobutanedione radical cation dissociation

The dissociation of 1,3-cyclobutanedione radical cation was studied by ab initio direct classical trajectory calculations at the BH&HLYP/6-31G(d) level of theory. A microcanonical ensemble using normal mode sampling was constructed by distributing 10 kcal/mol of excess energy above the transition state for the tautomerization of the keto-enol to the diketo form. A total of 210 trajectories were run starting from this transition state, yielding chemically activated 1,3-cyclobutanedione radical cation. The majority of the trajectories resulted in CH2CO+· + CH2CO, with the activated CC bond breaking nearly twice as often as the spectator CC bond. The non-statistical behavior is observed because the rate of energy redistribution within the molecule is comparable to or slower than the dissociation rates. In addition to the expected products, dissociation to CH2COCH 2 +· + CO and formation of a proton-transferred product HCCO· + CH3CO+ were also seen in some of the trajectories.

Influence of Structural Dynamics on Polarization Energies in Anthracene Single Crystals

Normal mode sampling and molecular dynamics simulations are coupled to a valence-bond/Hartree−Fock approach to evaluate the impact of the lattice and molecular vibrations on site energies in anthracene single crystals. The calculations are conducted in the temperature range 0−400 K and show substantial contributions from high-frequency modes, which calls for a quantum-mechanical model even at room temperature. External reorganization energies are also obtained from these modeling studies and found to be much smaller than their internal counterparts. Implications for charge transport in organic single crystals are discussed.

Ab Initio Classical Trajectory Study of the Dissociation of Neutral and Positively Charged Methanimine (CH2NHn+ n = 0−2)

The structures and energetics of the reactants, intermediates, transition states, and products for the dissociation of methanimine neutral, monocation, dication, and trication were calculated at the CBS-APNO level of theory. The dissociations of the neutral, monocation, and dication were studied by ab initio direct classical trajectory calculations at the B3LYP/6-311G(d,p) level of theory. A microcanonical ensemble using quasiclassical normal mode sampling was constructed by distributing 200, 150, and 120 kcal/mol of excess energy above the local minima of the neutral, singly, and doubly charged species, respectively. Many of the trajectories dissociate directly to produce H+, H atom, or H2. However, for a fraction of the cases, substantial migration of the hydrogen occurs within the molecule before dissociation. The preferred dissociation product for the neutral and the monocation is hydrogen atom. Elimination of H(2) was seen in 20% of the trajectories for the neutral and in 5% of the trajectories for the monocation. Dissociations of the dication and trication produced H+ rather than H atom. HCNH+ was formed in 85-90% of the dissociating trajectories for the monocation and dication.

Large Nonstatistical Branching Ratio in the Dissociation of Pentane-2,4-dione Radical Cation: An Ab Initio Direct Classical Trajectory Study

The dissociation of pentane-2,4-dione radical cation has been studied by ab initio direct classical trajectory calculations at the MP2/6-31G(d) level of theory. A bond additivity correction has been used to improve the MP2 potential energy surface (BAC-MP2). A microcanonical ensemble was constructed using quasiclassical normal-mode sampling by distributing 10 kcal/mol of excess energy above ZPE for the transition state for the tautomerization of the enol with a terminal double bond, 4-hydroxypent-4-en-2-one radical cation, to the diketo form. A total of 244 trajectories were run starting from this transition state, yielding pentane-2,4-dione radical cation and depositing energy in the terminal CC bond. As a result, the branching ratio for dissociation of the terminal CC bond versus the interior CC bonds is significantly larger than expected from RRKM theory. The branching ratio for the dissociation of the two interior CC bonds is ∼20:1, with the one closest to the activated methyl breaking more often. Since the two interior bonds are equivalent and should dissociate with equal probability, this branching ratio represents a very large deviation from statistical behavior. A simple kinetic scheme has been constructed to model the dissociation rates. The nonstatistical behavior is seen because the rate of energy flow within the molecule is comparable to or less than the rates of dissociation for the activated system. In addition to the expected dissociation products, some of the trajectories also lead to the formation of an ester-like product, prop-1-en-2-yl acetate radical cation.

RRKM and direct MP2/6-31G(d,p) quasiclassical trajectory study of the H 2 elimination in the photodissociation of vinyl chloride at 193 nm

RRKM calculations and MP2/6-31G(d,p) quasiclassical trajectories (QCT) were performed to investigate the H 2 elimination in the photodissociation of vinyl chloride at 193 nm. The kinetic study involved nine elementary steps and the QCT investigation was made for the relevant reaction channels. The trajectories were initiated at the transition states using a microcanonical, quasiclassical normal-mode sampling. Comparison between the calculated product energy distributions and the available experimental data is discussed. Overall, the present results (RRKM and QCT) suggest that the direct three-center elimination, which may occur via two similar transition states, is the dominant mechanism for H 2 production.

Rovibrational distributions of HF in the photodissociation of vinyl fluoride at 193 nm: a direct MP2 quasiclassical trajectory study.

Quasiclassical trajectory calculations were performed to calculate rovibrational distributions of the nascent HF fragment in the photodissociation of vinyl fluoride at 193 nm. The trajectories were initiated at the transition states of the four-center (4C) and three-center (3C) HF elimination channels, using a microcanonical, quasiclassical normal-mode sampling. In general, the calculated distributions are in reasonably good agreement with experiment. In particular, the trajectory distributions show bimodal character, although not as pronounced as that observed experimentally. The calculations predict that the 3C and 4C distributions are rather similar to each other, which suggests that the low-J and high-J components of the rotational distributions cannot be specifically assigned to each of these channels.

Classical trajectory studies of the photodissociation reaction of sym-triazine

Classical trajectory simulations were performed for the photodissociation of sym-triazine on an interpolated potential energy surface, for energies equivalent to 193, 248, 266, 285, and 295 nm light. Initial conditions were chosen using microcanonical normal mode sampling. The calculated average translational and rotational energies of the HCN products were in good agreement with some experimental results. The shapes of the product translational energy distributions are also in good agreement with experiment for low energy light, but a little narrower than the experimental distributions at high energies.

Dissociation of acetone radical cation (CH3COCH3+˙ → CH3CO++CH3˙): An ab initio direct classical trajectory study

The dissociation of acetone radical cation (CH3COCH3+˙ → CH3CO+ + CH3˙) has been studied by ab initio direct classical trajectory calculations at the MP2/6-31G(d) level of theory. A total of 247 trajectories were initiated at the transition state of the keto–enol tautomerization. A microcanonical ensemble using quasiclassical normal mode sampling was constructed by distributing 10 kcal mol−1 of excess energy above the barrier. The dissociation is found to favor the loss of the newly formed methyl group in agreement with experiments. The branching ratio of methyl loss was calculated to be 1.53 ± 0.20 which is fortuitously in very good agreement with the experimental ratio of 1.55 for 8–12 kcal mol−1 excess energy. Nearly 50% of the available energy is retained by the acetyl fragment as vibrational energy. The methyl fragment has very little vibrational energy but receives ca. 25% of the available energy in translation. The translational energy distribution of the methyl radicals is bimodal with the newly formed methyl having higher average translational energy than that the methyl derived from the spectator group.

A Direct Classical Trajectory Study of HCl Elimination from the 193 nm Photodissociation of Vinyl Chloride

Direct AM1 with specific reaction parameters (AM1-SRP) trajectory calculations were performed to investigate the photodissociation of vinyl chloride at 193 nm. Marked differences were found between the trajectories run from the transition states to the products for the four-center (4C) and three-center (3C) dissociation channels. The average internal energy of HCl is found to be about 2.5 times larger for the 4C channel. The calculated translational and HCl vibrational distributions, especially those obtained with a quasi-classical normal-mode sampling method, are in good agreement with the experimental data. The results are consistent with a concerted but nonsynchronous 3C mechanism, producing a lifetime for vinylidene to acetylene isomerization in line with the experimental observations. The theoretical HCl rotational distributions are in reasonable agreement with experiment, although they do not support the experimental interpretation of Lin et al. [J. Chem. Phys. 2001, 114, 160], who attributed low-J and high-J components to the 4C and 3C channels, respectively.

Initial conditions and paths in semiclassical tunneling

The choice of initial conditions and tunneling paths for semiclassical tunneling are investigated for a two-dimensional double-well system which is often used for modeling hydrogen tunneling in unimolecular isomerization processes. Tunneling effects are incorporated in a standard trajectory simulation. The exact quantum mechanical calculations are performed for comparison. We find that tunneling splittings are in close agreement for initial conditions obtained by using adiabatic switching and normal-mode sampling methods. These results suggest that the approximate normal-mode sampling can be used to compute accurate level splittings. We also find that simple straight-line tunneling paths work well from the sudden to intermediate regimes. However, the method is less accurate for the intermediate cases.

Sampling of semiclassically quantized polyatomic molecule vibrations by an adiabatic switching method: Application to quasiclassical trajectory calculations

We apply the adiabatic switching (AS) method to determine the polyatomic classical motions that correspond to selected vibrational quantum states on multidimensional, anharmonic potential energy surfaces, and use these semiclassically quantized motions as initial conditions for quasiclassical trajectory (QCT) calculations of state-to-state reaction dynamics. Specifically, we calculate the classical motion corresponding to the quantum mechanical zero-point vibration of deuterated methane, CD4, and run QCT calculations on the H+CD4→HD(v′,j′)+CD3 reaction. The distribution of CD4 vibrational zero-point energy (ZPE) associated with the AS-sampled motions is compared with that from normal-mode-sampled motions. The spread of total zero-point energy in the AS calculations is much narrower than with normal-mode sampling, and the ZPE’s are appropriately shifted to lower energy due to anharmonic effects. Reverse adiabatic switching is used as an indirect check of the quantum numbers of the adiabatically sampled motion, but numerical limitations made this test inconclusive. The AS method thus appears to be superior to normal-mode sampling, but this superiority cannot be demonstrated conclusively for the fully anharmonic CD4 potential. However, the AS method is shown to perform very well for transformation from one CD4 harmonic potential to another and for transformation from an harmonic to an anharmonic, but decoupled potential in which CD4 is described by Morse oscillators. Evidence is presented that suggests the AS calculations are limited by numerical inaccuracies or intrinsic features of the potential energy surface, both of which are unavoidable. H+CD4→HD(v′,j′)+CD3 QCT calculations of state-to-state dynamics using CD4 with no ZPE, the ZPE from AS sampling, and the ZPE from normal-mode sampling are reported and compared.