Direct dynamics simulations using Hessian-based predictor-corrector integration algorithms

Tell Machine Learning Potentials What They Are Needed For: Simulation-Oriented Training Exemplified for Glycine

Machine learning potentials (MLPs) are widely applied as an efficient alternative way to represent potential energy surfaces (PESs) in many chemical simulations. The MLPs are often evaluated with the root-mean-square errors on the test set drawn from the same distribution as the training data. Here, we systematically investigate the relationship between such test errors and the simulation accuracy with MLPs on an example of a full-dimensional, global PES for the glycine amino acid. Our results show that the errors in the test set do not unambiguously reflect the MLP performance in different simulation tasks, such as relative conformer energies, barriers, vibrational levels, and zero-point vibrational energies. We also offer an easily accessible solution for improving the MLP quality in a simulation-oriented manner, yielding the most precise relative conformer energies and barriers. This solution also passed the stringent test by diffusion Monte Carlo simulations.

Unsupervised Machine Learning Neural Gas Algorithm for Accurate Evaluations of the Hessian Matrix in Molecular Dynamics

The Hessian matrix of the potential energy of molecular systems is employed not only in geometry optimizations or high-order molecular dynamics integrators but also in many other molecular procedures, such as instantaneous normal mode analysis, force field construction, instanton calculations, and semiclassical initial value representation molecular dynamics, to name a few. Here, we present an algorithm for the calculation of the approximated Hessian in molecular dynamics. The algorithm belongs to the family of unsupervised machine learning methods, and it is based on the neural gas idea, where neurons are molecular configurations whose Hessians are adopted for groups of molecular dynamics configurations with similar geometries. The method is tested on several molecular systems of different dimensionalities both in terms of accuracy and computational time versus calculating the Hessian matrix at each time-step, that is, without any approximation, and other Hessian approximation schemes. Finally, the method is applied to the on-the-fly, full-dimensional simulation of a small synthetic peptide (the 46 atom N-acetyl-l-phenylalaninyl-l-methionine amide) at the level of DFT-B3LYP-D/6-31G* theory, from which the semiclassical vibrational power spectrum is calculated.

Transition between [R]- and [S]-stereoisomers without bond breaking

The fifty-year old proposal of a nondissociative racemization reaction of a tetracoordinated tetrahedral center from one enantiomer to another via a planar transition state by Hoffmann and coworkers has been explored by many research groups over the past five decades. A number of stable molecules with planar tetracoordinated and higher-coordinated centers have been designed and experimentally realized; however, there has not been a single example of a molecular system that can possibly undergo such racemization. Here we show examples of molecular species that undergo inversion of stereochemistry around tetrahedral centers (Si, Al- and P+) either via a planar transition state or an intermediate state using quantum mechanical, ab initio quasi-classical dynamics calculations, and Born-Oppenheimer molecular dynamics (BOMD) simulations. This work is expected to provide potential leads for future studies on this fundamental phenomenon in chemistry.

Temporally Coarse-Grained All-Atom Molecular Dynamics Achieved via Stochastic Padé Approximants

A Padé approximant scheme for realizing the discrete-time evolution of the state of a many-atom system is introduced. This temporal coarse-graining scheme accounts for the underlying Newtonian physics and avoids the need for construction of spatially coarse-grained variables. Newtonian physics is incorporated through short molecular dynamics simulations at the beginning of each of the large coarse-grained timesteps. The balance between stochastic and coherent dynamics expressed by many-atom systems is captured via incorporation of the Ito formula into a Padé approximant for the time dependence of individual atom positions over large timesteps. Since the time for a many-atom system to express a characteristic ensemble of atomic velocity fluctuations is typically short relative to the characteristic time of large-scale atomic displacements, a computationally efficient and accurate temporal coarse-graining of the atom-resolved Newtonian dynamics is formulated, denoted all-atom Padé-Ito molecular dynamics (APIMD). Evolution of the system over a time step much longer than that required for standard molecular dynamics (MD) is achieved via incorporation of information from the short MD simulations into a Padé approximant extrapolation in time. The extrapolated atomic configuration is subjected to energy minimization and, when needed, thermal equilibration so as to avoid occasional unphysical close encounters deriving from the Padé approximant extrapolation and to represent configurations appropriate for the temperature of interest. APIMD is implemented and tested via comparison with traditional MD simulations of five phenomena: (1) pertussis toxin subunit deformation, (2) structural transition in a T = 1 capsid-like structure of HPV16 L1 protein, (3) coalescence of argon nanodroplets, and structural transitions in dialanine in (4) vacuum, and (5) water. Accuracy of APIMD is demonstrated using semimicroscopic descriptors (rmsd, radius of gyration, residue-residue contact maps, and densities) and the free energy. Significant computational acceleration relative to traditional molecular dynamics is illustrated.

Can reactions follow non-traditional second-order saddle pathways avoiding transition states?

We report here an ab initio (CASSCF/6-31+G*) trajectory simulation study on the mechanisms of the denitrogenation of 1-pyrazoline and its subsituted analogue that reveals reaction pathways via a high energy second-order saddle (SOS) region. This mechanism involves the molecule adopting a five-membered planar structure contrary to the traditional boat-like transition state. The SOS offers a trifurcation point where a pathway branches into three, different from the single pathway associated with a transitions state. We observe that the molecules following the SOS path exhibit distinctive dynamical features and form products with high translational energies and low rotational energies compared to those following the traditional pathways. In addition, the SOS pathway provides an alternative mechanism for the formation of stereo-selective products. Interestingly, although the reaction proceeds via a trimethylene diradical intermediate, the simulations show that the product cyclopropane is formed with a major single inversion of the configuration consistent with experimental observations. They also reveal mechanisms that do not follow the minimum energy paths and exhibit non-statistical dissociation dynamics.

A new (multi-reference configuration interaction) potential energy surface for H2CO and preliminary studies of roaming

We report a new global potential energy surface (PES) for H2CO, based on precise fitting of roughly 67 000 MRCI/cc-pVTZ energies. This PES describes the global minimum, the cis- and trans-HCOH isomers, and barriers relevant to isomerization, formation of the molecular (H2+CO) and radical (H+HCO) products, and the loose so-called roaming transition-state saddle point. The key features of the PES are reviewed and compared with a previous PES, denoted by PES04, based on five local fits that are 'stitched' together by switching functions (Zhang et al. 2004 J. Phys. Chem. A108, 8980-8986 (doi:10.1021/jp048339l)). Preliminary quasi-classical trajectory calculations are performed at the total energy of 36 233 cm-1 (103 kcal mol-1), relative to the H2CO global minimum, using the new PES, with a particular focus on roaming dynamics. When compared with the results from PES04, the new PES findings show similar rotational distributions, somewhat more roaming and substantially higher H2 vibrational excitation.This article is part of the themed issue 'Theoretical and computational studies of non-equilibrium and non-statistical dynamics in the gas phase, in the condensed phase and at interfaces'.

Recent advances in quantum scattering calculations on polyatomic bimolecular reactions

This review surveys quantum scattering calculations on chemical reactions of polyatomic molecules in the gas phase published in the last ten years.

A high order predictor–corrector integration algorithm for first principle chemical dynamics simulations

H. Wu submitted the paper without having obtained permissions for submission from some of the co-authors or some of the major contributors, and hereby requested to retract the publication.

Photodissociation of acryloyl chloride at 193 nm: interpretation of the product energy distributions, and new elimination pathways

The use of an automated TS search method leads to the finding of novel HCl elimination pathways.

HCN elimination from vinyl cyanide: product energy partitioning, the role of hydrogen–deuterium exchange reactions and a new pathway

A schematic diagram of HCN elimination channels from vinyl cyanide including a new CCdiss pathway.

Mechanisms and dynamics of protonation and lithiation of ferrocene

Ab initio direct dynamics studies reveal that protonation and lithiation of ferrocene follow in general an exo mechanism.

Dynamic effects dictate the mechanism and selectivity of dehydration-rearrangement reactions of protonated alcohols [Me2 (R)CCH(OH2 )Me](+) (R=Me, Et, iPr) in the gas phase

The gas-phase dehydration-rearrangement (DR) reactions of protonated alcohols [Me2 (R)CCH(OH2 )Me](+) [R=Me (ME), Et (ET), and iPr (I-PR)] were studied by using static approaches (intrinsic reaction coordinate (IRC), Rice-Ramsperger-Kassel-Marcus theory) and dynamics (quasiclassical trajectory) simulations at the B3LYP/6-31G(d) level of theory. The concerted mechanism involves simultaneous water dissociation and alkyl migration, whereas in the stepwise reaction pathway the dehydration step leads to a secondary carbocation intermediate followed by alkyl migration. Internal rotation (IR) can change the relative position of the migrating alkyl group and the leaving group (water), so distinct products may be obtained: [Me(R)CCH(Me)Me⋅⋅⋅OH2 ](+) and [Me(Me)CCH(R)Me⋅⋅⋅OH2 ](+) . The static approach predicts that these reactions are concerted, with the selectivity towards these different products determined by the proportion of the conformers of the initial protonated alcohols. These selectivities are explained by the DR processes being much faster than IR. These results are in direct contradiction with the dynamics simulations, which indicate a predominantly stepwise mechanism and selectivities that depend on the alkyl groups and dynamics effects. Indeed, despite the lifetimes of the secondary carbocations being short (<0.5 ps), IR can take place and thus provide a rich selectivity. These different selectivities, particularly for ET and I-PR, are amenable to experimental observation and provide evidence for the minor role played by potential-energy surface and the relevance of the dynamics effects (non-IRC pathways, IR) in determining the reaction mechanisms and product distribution (selectivity).

Synchronous concerted multiple-body photodissociation of oxalyl chloride explored by ab initio-based dynamics simulations

Photo-induced multiple body dissociation is of fundamental interest in chemistry and physics. A description of the mechanism associated with n-body (n ≥ 3) photodissociation has proven to be an intriguing and yet challenging issue in the field of chemical dynamics. Oxalyl chloride, (ClCO)2, is the sole molecule reported up to date that can undergo four-body dissociation following absorption of a single UV photon, with a rich history of mechanistic debate. In the present work, the combined electronic structure calculations and dynamics simulations have been performed at the advanced level, which provides convincing evidence for resolving the mechanistic debate. More importantly, synchronous and asynchronous concertedness were explored for the first time for the (ClCO)2 photodissociation, which is based on the simulated time constants for the C-C and C-Cl bond fissions. Upon photoexcitation of (ClCO)2 to the S1 state, the adiabatic C-C or C-Cl fission takes place with little possibility. The four-body dissociation to 2Cl((2)P) and 2CO((1)Σ) was determined to a dominant channel with its branch of ∼0.7, while the three-body dissociation to ClCO((2)A(')) + CO((1)Σ) + Cl((2)P) was predicted to play a minor role in the (ClCO)2 photodissociation at 193 nm. Both the four-body and three-body dissociations are non-adiabatic processes, which proceed in a synchronous concerted way as a result of the S1 → S0 internal conversion. There is a little possibility for two-body dissociation to occur in the S0 and S1 states.

Ab initio and semi-empirical Molecular Dynamics simulations of chemical reactions in isolated molecules and in clusters

Recent progress in “on-the-fly” trajectory simulations of molecular reactions, using different electronic structure methods is discussed, with analysis of the insights that such calculations can provide and of the strengths and limitations of the algorithms available.

Selectivity and mechanisms driven by reaction dynamics: the case of the gas-phase OH(-) + CH3ONO2 reaction

Well-established statistical approaches such as transition-state theory based on high-level calculated potential energy profiles are unable to account for the selectivity observed in the gas-phase OH(-) + CH(3)ONO(2) reaction. This reaction can undergo bimolecular nucleophilic displacement at either the carbon center (S(N)2@C) or the nitrogen center (S(N)2@N) as well as a proton abstraction followed by dissociation (E(CO)2) pathway. Direct dynamics simulations yield an S(N)2:E(CO)2 product ratio in close agreement with experiment and show that the lack of reactivity at the nitrogen atom is due to the highly negative electrostatic potential generated by the oxygen atoms in the ONO(2) group that scatters the incoming OH(-). In addition to these dynamical effects, the nonstatistical behavior of these reactions is attributed to the absence of equilibrated reactant complexes and to the large number of recrossings, which might be present in several ion-molecule gas-phase reactions.

Chemical dynamics simulations of X- + CH3Y → XCH3 + Y- gas-phase S(N)2 nucleophilic substitution reactions. Nonstatistical dynamics and nontraditional reaction mechanisms

Extensive classical chemical dynamics simulations of gas-phase X(-) + CH(3)Y → XCH(3) + Y(-) S(N)2 nucleophilic substitution reactions are reviewed and discussed and compared with experimental measurements and predictions of theoretical models. The primary emphasis is on reactions for which X and Y are halogen atoms. Both reactions with the traditional potential energy surface (PES), which include pre- and postreaction potential energy minima and a central barrier, and reactions with nontraditional PESs are considered. These S(N)2 reactions exhibit important nonstatistical atomic-level dynamics. The X(-) + CH(3)Y → X(-)---CH(3)Y association rate constant is less than the capture model as a result of inefficient energy transfer from X(-)+ CH(3)Y relative translation to CH(3)Y rotation and vibration. There is weak coupling between the low-frequency intermolecular modes of the X(-)---CH(3)Y complex and higher frequency CH(3)Y intramolecular modes, resulting in non-RRKM kinetics for X(-)---CH(3)Y unimolecular decomposition. Recrossings of the [X--CH(3)--Y](-) central barrier is important. As a result of the above dynamics, the relative translational energy and temperature dependencies of the S(N)2 rate constants are not accurately given by statistical theory. The nonstatistical dynamics results in nonstatistical partitioning of the available energy to XCH(3) +Y(-) reaction products. Besides the indirect, complex forming atomic-level mechanism for the S(N)2 reaction, direct mechanisms promoted by X(-) + CH(3)Y relative translational or CH(3)Y vibrational excitation are possible, e.g., the roundabout mechanism.

The need for enzymatic steering in abietic acid biosynthesis: gas-phase chemical dynamics simulations of carbocation rearrangements on a bifurcating potential energy surface

Abietic acid, a constituent of pine resin, is naturally derived from abietadiene --a process that requires four enzymes: one (abietadiene synthase) for conversion of the acyclic, achiral geranylgeranyl diphosphate to the polycyclic, chiral abietadiene (a complex process involving the copalyl diphosphate intermediate) and then three to oxidize a single methyl group of abietadiene to the corresponding carboxylic acid. In previous work (Nature Chem.2009, 1, 384), electronic structure calculations on carbocation rearrangements leading to abietadienyl cation revealed an interesting potential energy surface with a bifurcating reaction pathway (two transition-state structures connected directly with no intervening minimum), which links two products--one natural and one not yet isolated from Nature. Herein we describe direct dynamics simulations of the key step in the formation of abietadiene (in the gas phase and in the absence of the enzyme). The simulations reveal that abietadiene synthase must intervene in order to produce abietadiene selectively, in essence steering this reaction to avoid the generation of byproducts with different molecular architectures.

Electronic structure theory study of the F(-) + CH(3)I → FCH(3) + I(-) potential energy surface

MP2 and DFT electronic structure theories, with the functionals OPBE, OLYP, HCTH407, BhandH, and B97-1 for the latter, were used to investigate stationary point properties on the F(-) + CH(3)I → FCH(3) + I(-) potential energy surface (PES). The aug-cc-pVDZ and aug-cc-pVTZ basis sets for C, H, and F, with Wadt and Hay's 3s3p valence basis and an effective core potential (ECP) for iodine, were employed for both MP2 and DFT. Single-point CCSD(T) calculations were also performed to obtain the complete basis set (CBS) limit for the stationary point energies. The CCSD(T)/CBS reaction exothermicity is only 5.0 kJ/mol different than the experimental value. MP2 and DFT do not predict the same stationary points on the PES. MP2 predicts the C(3v) F(-)-CH(3)I and FCH(3)-I(-) ion-dipole complexes and traditional [F-CH(3)-I](-) central barrier as stationary points, as well as a C(s) hydrogen-bonded F(-)-HCH(2)I complex and a [F-HCH(2)-I](-) transition state connecting this latter complex to the F(-)-CH(3)I complex. A CCSD(T)/CBS relaxed potential energy curve, calculated for the MP2 structures, shows that going from the F(-)-CH(3)I complex to the [F-CH(3)-I](-) TS is a barrierless process, indicating these two structures are not stationary points. This is also suggested by the DFT calculations. The structures and frequencies for CH(3)I and CH(3)Cl given by MP2 and DFT are in overall good agreement with experiment. The calculations reported here indicate that the DFT/B97-1 functional gives the overall best agreement with the CCSD(T) energies, with a largest difference of only 7.5 kJ/mol for the FCH(3)-I(-) complex.

A direct dynamics trajectory study of F- + CH(3)OOH reactive collisions reveals a major non-IRC reaction path

A direct dynamics simulation at the B3LYP/6-311+G(d,p) level of theory was used to study the F- + CH3OOH reaction dynamics. The simulations are in excellent agreement with a previous experimental study (J. Am. Chem. Soc. 2002, 124, 3196). Two product channels, HF + CH2O + OH- and HF + CH3OO-, are observed. The former dominates and occurs via an ECO2 mechanism in which F- attacks the CH3- group, abstracting a proton. Concertedly, a carbon-oxygen double bond is formed and OH- is eliminated. Somewhat surprisingly this is not the reaction path, predicted by the intrinsic reaction coordinate (IRC), which leads to a deep potential energy minimum for the CH2(OH)2...F- complex followed by dissociation to HF + CH2(OH)O-. None of the direct dynamics trajectories followed this path, which has an energy release of -63 kcal/mol and is considerably more exothermic than the ECO2 path whose energy release is -27 kcal/mol. Other product channels not observed, and which have a lower energy than that for the ECO2 path, are F- + CO + H2 + H2O (-43 kcal/mol), F- + CH2O + H2O (-51 kcal/mol), and F- + CH2(OH)2 (-60 kcal/mol). Formation of the CH3OOH...F- complex, with randomization of its internal energy, is important, and this complex dissociates via the ECO2 mechanism. Trajectories which form HF + CH3OO- are nonstatistical events and, for the 4 ps direct dynamics simulation, are not mediated by the CH3OOH...F- complex. Dissociation of this complex to form HF + CH3OO- may occur on longer time scales.