Neural Network Potential Energy Surface for the low temperature Ring Polymer Molecular Dynamics of the H2CO + OH reaction
NNeural Network Potential Energy Surface for the low temperature Ring Polymer Molecular Dynamics of theH CO + OH reaction
Pablo del Mazo-Sevillano, a) Alfredo Aguado, and Octavio Roncero b)1) Unidad Asociada UAM-CSIC, Departamento de Química Física Aplicada, Facultad de Ciencias M-14,Universidad Autónoma de Madrid, 28049, Madrid, Spain Instituto de Física Fundamental (IFF-CSIC), C.S.I.C., Serrano 123, 28006 Madrid,Spain
A new potential energy surface (PES) and dynamical study are presented of the reactive process between H CO + OHtowards the formation of HCO + H O and HCOOH + H. In this work a source of spurious long range interactions insymmetry adapted neural network (NN) schemes is identified, what prevents their direct application for low temperaturedynamical studies. For this reason, a partition of the PES into a diabatic matrix plus a NN many body term hasbeen used, fitted with a novel artificial neural networks scheme that prevents spurious asymptotic interactions. Quasi-classical trajectory and ring polymer molecular dynamics (RPMD) studies have been carried on this PES to evaluate therate constant temperature dependence for the different reactive processes, showing a good agreement with the availableexperimental data. Of special interest is the analysis of the previously identified trapping mechanism in the RPMDstudy, which can be attributed to spurious resonances associated to excitations of the normal modes of the ring polymer.Accepted in J. of Chem. Phys. (2021)
I. INTRODUCTION
The raise of the reaction rate constant at low tempera-tures measured in different laboratories for several organicmolecules (OM) with OH and other reactants opened a pos-sible solution for the enigma of how organic molecules aregenerated in cold dense molecular clouds, below 10 K, in theinterstellar medium. These reactions usually present reactionbarriers, which make impossible the reaction at low tempera-tures in gas phase. For that reason, it was assumed that thesemolecules were formed on cosmic ices, in cold molecularclouds at about 10 K , and then released to gas phase as theyevolve to hotter structures , such as hot cores and corinos,where they were observed . However, recently severalobservations of OMs were made in UV-shielded cold coresat 10 K , what introduces the need of changing the model.Several hypotheses were introduced, such as the possibility ofchemidesorption , the incidence of cosmic rays and ofsecondary UV photons, that may induce the photodesorption.However, OMs usually present dipoles and strong binding en-ergies, and recent experiments performed by illuminating cos-mic ice precursors showed that only their photofragments aredesorbed into gas phase . This may suggest that the ori-gin of OMs in gas phase should be a combined mechanism, inwhich they are formed on ices, then their photofragments re-leased to gas phase, where they suffer a reprocessing to finallyend in the formation of the OM.The reprocessing could be accelerated in gas phase by thehuge increase of the reaction rate observed in the laboratorybelow 100 K . This increase is explained by the formationof a complex between the reactants, from where they tunnel a) Electronic mail: [email protected] b) Electronic mail: [email protected] to form the products . However, several transition state the-ory (TST) studies of the CH OH + OH reaction did notget a so steep raise of the rate constant at low temperature,but it was much slower than the experimental results. They allattributed the difference to pressure effects, which at micro-scopic level may be attributed to the formation of complexes,between the reactants and of them with the buffer gas used inthe expansions. These complexes could experience secondarycollisions adding the extra energy needed to overcome the bar-rier. However, due to the low densities of molecular cloudsthe zero-pressure rate constant is needed in their models. Itis therefore mandatory to check if TST methods can describequantitatively the tunneling in the deep regime.In this line, we have studied the dynamics of the reactionsof OH with H CO and CH OH . For this purpose, fulldimensional potential energy surfaces (PESs) were fitted toaccurate ab initio calculations . For H CO + OH reactionthe QCT rate constant yield a semiquantitative agreement withexperimental results , while for CH OH + OH the simulatedrate constant was far too low . To check if quantum effectscould explain the difference, ring polymer molecular dynam-ics (RPMD) calculations were performed on the two PESs ,finding an excellent agreement in the case of methanol, buta similar semiquantitative agreement for the case of H CO.These good agreement could explain the experimental resultsand provide a good estimate of the zero-pressure rate constantat low temperature, but several problems still persist. First, be-low 100 K RPMD leads to a high trapping probability, form-ing complexes with extremely long lifetimes (larger than sev-eral hundreds of nanoseconds), whose trajectories cannot befinished. Therefore, the reaction rate has to be evaluated bymultiplying the trapping rate constant by a ratio to productsthat is inferred from TST method or QCT calculations . Thiscould also be due to the PES, and the second problem is toimprove the accuracy of the fit. This is the goal of this paper,by applying Artificial Neural Networks (ANN) techniques to1 a r X i v : . [ phy s i c s . c h e m - ph ] F e b t a new PES to check if the high trapping rate persist, in oneside, and to improve the comparison with experimental data.In this work we shall explore different alternatives formultidimensional fitting with ANN techniques, including thesymmetry either with permutationally invariant polynomials,PIP-NN , or with fundamental invariants, FI-NN . Inaddition, the channel towards HCOOH+H products, absent inthe previous PES , will be included. We shall focus the at-tention in the problem of including long-range interactions,which are crucial to simulate the dynamics at low tempera-ture. We shall then use the new PES to perform QCT andRPMD calculations to compare with the experimental resultsand to analyze the formation of long-lived complexes.This paper is organized as follows. Firstly, a descriptionof ab initio calculations is provided in section 2. A generaldescription of artificial neural networks is then provided insection 3, focusing on the description of a new many-bodyscheme to accurately include the long-range interactions cru-cial in the dynamics at low temperature in section 4. With this,an overall description of the PES will be provided in section5, showing the most important features. In section 6, we willfocus on the dynamical results on this new PES, employingboth QCT and RPMD methodologies, comparing with previ-ous theoretical and experimental results. Finally some conclu-sions will be extracted together with several crucial questionsto be addressed in the future. II.
AB INITIO
CALCULATIONS
The reactions H H + H ( Π ) → HC + H (1a) → HCH + H (1b)are studied using the explicitly correlated coupled clustermethod, RCCSD(T)-F12a/cc-pVDZ , as implemented in theMOLPRO package . CCSD(T) calculations are considered abenchmark, and it has recently been shown that CCSD(T)-F12is able to reach results close to CCSD(T)/CBS, even with rel-atively small basis sets . The effect over reactivity froman excited electronic state has been despise from the results ofMRCI/cc-pVDZ calculations along the RCCSD(T)-F12 min-imum energy path, since it becomes highly repulsive as thesystems approach the transition states of either HCO + H Oor HCOOH + H reaction channels. The reader is referred tothe Supplementary Information for a further description.The calculated stationary points at this level of theory arerather similar to those previously reported . The energyof the transition state to form HCO + H O is very low (of ≈ O, called TS1 (see figure 3), with an energy of 27 . − . CO plane. Then, the H CO bends and, at the sametime, the OH gets closer to the carbon atom. Finally, one ofthe H CO hydrogens leaves, forming HCOOH. The formationof c − HCOOH or t − HCOOH depends on whether the OH ro-tates or not from the stationary point RC4.For the present fit, around 180000 ab initio points wereadded to the previously calculated for the PES of Ref. . Thiswas needed to increase the accuracy of reaction 1a, already ex-plored in Ref. , and to describe, for the first time, the chan-nel towards HCOOH + H, in reaction 1b. The new pointswere calculated iteratively including points for the minimumenergy paths (MEPs) and normal modes along them evalu-ated on intermediate fits. Also QCT and RPMD trajecto-ries were ran on intermediate versions of the fits to popu-late physically accessible configurations. For each new point,the Euclidean distance to all the other points is evaluated as d ab = (cid:113) ∑ i = ( d ai − d bi ) , where d ai is the i th interatomic dis-tance of geometry a . III. ARTIFICIAL NEURAL NETWORKS
The use of Permutationally Invariant Polynomials (PIP)with Rydberg-type functions has been widely used to describethe PES of polyatomic systems with no many atoms .In these methods, few non-linear exponential parameters areused, while the linear coefficients of each term are fitted with aleast square method. In order to obtain a PES that is invariantwith respect to translation, rotation and inversion, the fit of thepotential is carried out as a polynomial in terms of the inter-nuclear distances . As the number of internuclear distancesincreases more than the number of independent coordinate itmeans that, even considering low orders of the fitted polyno-mial, the number of coefficients necessary to carry out the fitincreases enormously with a small increase of the number ofatoms in the system. This implies a limitation in the appli-cation of these methods to systems of more than 6-8 atoms.Instead, in this work we make use of neural networks (NN) tofit a new PES.Feed forward neural networks (FFNN) are a type of func-tions able to transform a vectorial input into a vectorial output,in contrast to other architectures devoted to work on matrix or graph representations. This kind of neural networks hasbeen of great interest in the field of PES fitting, where a molec-ular representation of the system is given to the neural net-work, returning its potential energy. This algorithm works ina supervised way, since the neural network is trained againsta set of target energies, e.g. ab initio , usually minimizing theroot mean squared error (RMSE) with the predictions.A FFNN is characterized by the number of layers and thenumber of neurons at each layer, increasing the flexibility of2 .. ... ... Eb b b I I I I m H H n H H H o InputLayer 1st HiddenLayer 2nd HiddenLayer OutputLayer
FIG. 1. Computational graph of a FFNN with two hidden layers. the function as more layers and neurons are included. Thecomputational graph for such a FFNN with two hidden lay-ers is presented in figure 1. A simple rule is given to moveforward in the layer structure, from the input (molecular ge-ometry) to output (energy): H lj = σ (cid:32) ∑ i w lji H l − i + b lj (cid:33) , (2)where H lj is the j th value of the l th hidden layer, w li j and b lj are the trainable weights and bias in that layer and σ is thetransfer function which is, in general, not linear.Typical choices of σ are hyperbolic tangent or sigmoidfunctions between the hidden layers in the context of PES fit-ting , and a linear function between the last hidden layer andoutput layer to map the whole range of real numbers. In thiswork we used a logistic function, σ ( x ) = / ( + exp ( − x )) ,between the hidden layers. The input representation for eachANN is defined as the negative exponential of each inter-atomic distances. The corresponding PIP of FI polynomialsare evaluated on them whenever the permutational symmetryis considered. Finally, all input values are standardize to zeromean and unit standard deviation. The PIP and FI functionsused are presented in the Supplementary Information.The NN structure will be expressed in the usual way, I-H - · · · -H n -O, with I and O the number of neurons in the inputand output layers, and n hidden layers, with H i neurons in the i th hidden layer. PIP-NN and FI-NN schemes, including thepermutational symmetry of identical atoms, have many prop-erties in common so artificial neural networks (ANN) will beused to refer indistinctly to either PIP-NN or FI-NN in thispaper.For each permutational symmetry group, the FI are cal-culated with the King’s algorithm implemented in thecomputer algebra system Singular . For many systemsthese FI have already been evaluated and FORTRAN rou-tines for their evaluation are provided in https://github.com/pablomazo/FI , expanded from the original repository https://github.com/kjshao/FI .For training ANNs we have developed NeuralPES, aPython code based on PyTorch library, meant to assist at each step of the training process, beginning with data prepro-cessing, hyperparameter tuning and training itself. IV. LONG RANGE INTERACTION WITH NEURALNETWORKS
The long range behaviour of a PIP PES was analyzed inRef. , concluding that certain terms of the PIP expansion leadto spurious interactions in the asymptotic region due to poly-nomials involving distances between unconnected fragments.This was solved removing some of the polynomials that con-nect the fragments asymptotically in what is called “PurifiedBasis” . An equivalent solution has been recently proposedby J. Li et al. to be used in the fit of many body (MB) termswith PIP-NN by removing those PIPs relating unconnecteddistances.In this paper we demonstrate that, on top of these terms ofunconnected distances, there exist another source of spuriousinteractions due to the non linearity of the transfer function inthe ANN. For doing so, a simplified A BC system is consid-ered, focusing on the asymptote A · · · BC, such that there isno interaction between them. A PIP PES is defined as V = ∑ i c i PIP ( i ) (3)where c i are the coefficients to fit and PIP(i) is the i th permu-tational invariant polynomial. The generator of all the PIPs( G ) for this particular permutational symmetry is: G = h l A A h l BC (cid:104) h l A B · h l A C · h l A B · h l A C + h l A B · h l A C · h l A B · h l A C (cid:105) (4)where h i j = f ( d i j ) , with d i j the distance between atoms i and j , and l i j is and exponent. It is common to define h i j = exp ( − α d i j ) , such that h i j tends to zero as the dis-tance tends to infinity. In this situation only the PIPs where l = l = l = l = G ( ∞ ) = h l A A h l BC (5)Setting max ( l i ) = PIP = { h A A , h BC , h A A · h BC } (6)Therefore, the PES in the asymptote is evaluated as V ( ∞ ) = c · h A A + c · h BC + c · h A A · h BC (7)Notice that, through the third term in the expansion a force be-tween A and BC fragments arises, even though they were setat an infinite distance. This PIP set could be purified by justremoving the third function, hence disconnecting both frag-ments. At this point, the asymptotic potential energy surfacewould depend linearly on both fragment distances, so no spu-rious interactions are introduced up to now.The next step is to show that even when the purified PIP setis used, PIP-NN introduces spurious interactions. In a PIP-NN3ith only one hidden layer, the energy is evaluated as V = ∑ j w E j H j + b E (8) H j = σ (cid:32) ∑ i w ji PIP ( i ) + b j (cid:33) , (9)where w and b are the weight matrices and bias vectors, H is the hidden layer and σ is the transfer function, which ingeneral is not linear. The values of the neurons of the hiddenlayer evaluated on the purified PIP set become H j = σ (cid:0) w j · h A A + w j · h BC + b j (cid:1) . (10)Since σ is not linear, a connection between both fragmentsis introduced and, with that, an spurious force between themarises. ∂ V / ∂ h A A , which ultimately relates with the force A and A experience is ∂ V ∂ h A A = ∑ j w E j w j σ (cid:48) (cid:0) w j · h A A + w j · h BC + b j (cid:1) , (11)being clear that there is always a dependence with the BCfragment.This same derivation can be develop with a set of FI, arriv-ing to the same conclusion. Moreover, this problem is mag-nified when more than one hidden layer is introduced. Thisspurious interactions lead to an energy transfer between thefragments, that will have a critical effect in dynamical studies,specially at low collision energies, where trajectories couldeven be reflected. We may then conclude at this point, that themain source of spurious interactions in a PIP-NN or FI-NN isthe non linearity of the transfer function.Although it would yield to a very low fitting error, the abovearguments show that it is not possible to fit the PES with oneANN since a residual interaction persists between the reac-tants at long distances. This causes small changes in the inter-nal energies of each fragment, introducing errors that stronglyaffect the dynamics at low energies, making inconvenient theuse of PIP or FI-NN PES that expand the whole configurationspace.One possible solution is to separate the PES in two regions,short and long range, connected by a switching function .The long range term would be a well behaved separable func-tion and the short range an ANN PES. The main difficulty withthis functional form is to tune the switching function to pro-duce a smooth change between both regions, specially as thedimensionality of the system increases. This can also havea huge effect when studying low energy or temperature dy-namics of the system, due to the possible presence of spuriousmatching problems.In this work we propose to use a partition of the potentialenergy as the one used in the previous PESs , werea diabatic matrix, V diab , is used to produce a first order PEScorrected with a six-body term, as expressed by V = E diab + V C , (12)where E diab is the lowest eigenvalue of V diab . This matrix isof dimension N × N , being N the number of rearrangement channels, each of them describing reactans or products frag-ments, and the interaction between them, including long rangeinteractions. This matrix is analogous to those previouslyused , and is described in detail in the SI. In the presentimplementation of V diab , the accuracy of each fragment of thePES has been increased by replacing them by ANN fittings.The main advantage of this form is that the diabatic matrixcaptures the basic features, namely the long range interactionsin each rearrangement channel and a first approximation toshort range interactions, only using PESs for each of the frag-ments and their interactions on each rearrangement channel.This idea is based on triatomics-in-molecules (TRIM) anddiatomic-in-molecules (DIM) , that allows an accurate de-scription of the fragments of physical relevance as well as theinteractions among them, both short and long range. V. POTENTIAL ENERGY SURFACE: NN SIX-BODYTERM
For the six-body term, V C , we propose to use the partition V C = S · V ANN , (13)where S is the switching function, designed to be zero as thesystem tends to any of the asymptotic regions, removing thespurious interaction from the V ANN , the PIP-NN or FI-NNfunction, trained to the difference between the ab initio en-ergies and E diab , i . e . it should tend to zero at long distances.In this way, the switching function ensures once more that theMB term will go to zero in the desired regions, and the frag-ments in the asymptotic region will not be affected by spuriousinteractions. On top of that, notice that it is not necessary topurify the set of invariants anymore, since is the S functionwhich removes the spurious connections.In order to stress the importance of removing the long rangespurious interaction in ANN PES, the energy of the H COand OH fragments has been followed along two QCT trajec-tories, being both fragments well separated so that there isno physical interaction between them. In the first trajectory, V C = V ANN and in the second one V C = S · V ANN . The resultsare shown in figure 2. It is clear that when the S function isnot included to build the six-body term, both fragments expe-rience an spurious interaction that makes their energy changealong the trajectory. This transfer is of about 10 meV, com-parable to the relative kinetic energy at low collisional energy.This interaction is completely removed by the switching func-tion as it becomes evident by the constant energy of the frag-ments in dashed lines.The asymptotic channels to consider in this system areH CO + OH, HCO + H O and HCOOH + H. The switchingfunction ( S ) is a product of hyperbolic tangents depending onthe distances d CO , d CH and d CH , so permutational symmetryof H CO is preserved. V ANN is fitted with NeuralPES program through a FI-NNwith a structure 21-80-80-1, meaning that 21 fundamental in-variants are required for taking into account the permutationalsymmetry of the two formaldehyde hydrogen atoms. The FIsare evaluated over the negative exponential of the interatomic4 t /ps E / c m H CO energyOH energy
FIG. 2. Energy of H CO and OH fragments in two QCT trajecto-ries with the fragments set at a large distance so there is no physicalinteraction between them. With dashed lines a trajectory in which V C = S · V ANN and with solid lines one where V C = V ANN . distances. The reader is referred to the Supplementary Infor-mation for the hyperparameter description used here. TABLE I. RMSE evaluated over each energy range for the diabaticterm and the full PES. E / eV < Points E diab / meV V / meV . .
38 39 . . .
41 51 . . .
64 65 . . .
92 86 . . .
00 92 . . .
60 105 . The RMSE of E diab and the fitted PES are presented in tableI. The ANN V C term reduces the fitting error by almost afactor of 10 for the lower energy ranges.In figure 3 the ab initio and the PES minimum energy pathsfor reactions to HCO + H O and HCOOH + H are shown.In both cases we find a very good agreement between thepaths, with the only exception of a shoulder that appears inthe HCOOH + H products channel. This is not expected tohave any effect on the dynamics since it is well below the re-actants asymptote. We find an improvement with respect tothe previous PES in the description of the region of stationarypoint RC1, which was geometrically not well located. Also,another local minima on the reactants rearrangement channel,not related to the minimum energy path, have been improved,what is of importance for a reaction with such a roaming be-haviour.The short range region of the PES is compared against abinitio results in figure 4. In general, the PES is very well be-haved, showing a perfect agreement for both attractive andrepulsive regions.A very good agreement is also found in the normal modefrequencies on the different transition states as can be checkedin the table presented in the Supplementary Information. Highfrequency modes have been considerably improved, in partdue to a better description of the fragments PESs. It is alsofound that, in general, low frequency modes present greater differences, since they relate with flatter regions of the PESthat can be more challenging to fit.
VI. DYNAMICAL RESULTSA. QCT results
QCT calculations were carried with an extension of miQCTcode for N atom systems. The reaction cross section σ ( E ) is calculated at fixed collision energy with the reactants intheir ground vibrational and rotational states as: σ ( E ) = π b max P r ( E ) (14)where b max ( E ) is the maximum impact parameter and P r ( E ) the reaction probability at a given collision energy. The initialconditions of the internal degrees of freedom of the fragmentshave been calculated with the adiabatic switching method .The remaining initial conditions are set by random numbersaccording to the usual Monte Carlo method . The initialmaximum impact parameter is set according to the capturemodel for a dipole-dipole interaction shown in Fig. 5along with the converged b max .The reactive cross section for the formation of HCO + H Ois presented in figure 5, calculated from more than 5 × trajectories at each collisional energy within b max ( E ) . Thepresent QCT results have the same qualitative behaviour asthose obtained with a previous PES , and reactive cross sec-tion experiences a huge increase as the collisional energy de-creases below 100 meV, even at energies lower than the tran-sition state energy barrier, at 27 meV. This increase in the re-active cross section is explained by an enormous increase ofthe maximum impact parameter, while the reaction probabil-ity remains constant or slightly increases.The reaction at high relative collisional energy is character-ized by a direct mechanism in which OH and H CO collideand either react or fly apart. At low collisional energies thereaction is dominated by the long range, dipole-dipole, inter-action which rotates and orients the reactants as to maximizetheir interaction, capturing the trajectory for high values of theimpact parameter. This leads to a rotational excitation of thesystem as the reactants become closer that traps the system forlong times, since it is improbable to stop the rotation turning itinto translation. Direct and trapped QCT trajectories are sim-ilar to the RPMD ones, shown in figure 9, with the differencethat trapped trajectories live much less, as will be discussedbelow.The reaction probability does not got to zero for collisionalenergies below the potential barrier, since the zero-point en-ergy (ZPE) is enough to overcome it. Still, there must ex-ist couplings between the orthogonal degrees of freedom thatpromote the energy transfer among them. This coupling be-comes favoured by the roaming-like mechanism and the hugecomplex lifetimes experienced at low collisional energies.The QCT reaction rate constants have been calculated atconstant temperature, with the reactants in their ground vibra-tional states and considering their rotational distribution. Only5
S1RC1 PC1
RC1 RC4 TS3 TS5t-HCOOH + Hc-HCOOH + HTS2
FIG. 3. Minimum energy paths for the formation of HCO + H O (left) and HCOOH + H (right), both cis and trans rotamers. The geometriesof the stationary points are represented along the path.FIG. 4. Energy of the system when OH approaches H CO in theplane. On the top (below) panels, the H (O) of OH always facesthe H CO molecule. Black and blue isolines represent energies of0 meV and -50 meV, respectively. Distance is expressed in Å. TheH CO and OH fragments are only meant to serve as a reference, theyare not at scale. the two states that correlate with the ground spin-orbit state,OH( Π / ), react, so the electronic partition function is q e ( T ) = + exp ( − . / T ) , (15)and the reaction rate constant is evaluated as k ( T ) = q e ( T ) (cid:115) k B T π µ π b max ( T ) P r ( T ) , (16)where b max ( T ) and P r ( T ) are the maximum impact parameterand reaction probability at constant temperature, and µ is thereduced mass of the OH, H CO system.The calculated reaction rate constants are shown in figure 6along with the experimental results in the temperature rangefrom 22 K to 1200 K, together with rates calculated earlierwith another PES .The reaction rate constant behaviour with temperature isbasically the same for the two PESs above 100 K, but be-low this temperature there is a quantitative difference between E /meV10 ( E ) / a E /meV0255075100125 b m a x / a Capture model 10 E /meV012345 P r FIG. 5. QCT results at fixed collisional energy for the formation ofHCO + H O. On the top left panel, the maximum impact parameter,on the top right panel the reaction probability and in the lower panel,the reactive cross section. both of them, although both reproduce the increase experi-mentally shown. It is not surprising that the reaction rate con-stants calculated with the two PESs is not well reproduced athigh temperatures, since only the vibrational ground state ofthe fragments is being considered. Still, there is a factor ofabout 1.5 between our results and the experimental ones fromWang et al. at 300 K. This difference increases as the tem-perature lowers, with respect to experimental results and alsobetween both PES. The discrepancies between the two PESscould be attributed to variations in the transition state region.Due to the relatively small height of the TS and the high an-harmonicity of this region, tiny differences, below the fittingerror, can lead to differences in dynamical behaviour.The reaction to form HCOOH + H is secondary with respectto the formation of HCO + H O. In particular, R. A. Yetter6
200 400 600 800 1000 1200 1400 T /K10 k ( T ) / c m s H CO + OH HCO + H O PES v1This workOcaña et al . fitWang et al . fitWang 15% uncertainty FIG. 6. QCT reaction rate constant at constant temperature for HCO+ H O formation. With black points, the results obtained in thiswork. The error bars are two times the standard error of the MonteCarlo integration . With red points, the results from our previousPES , With solid lines, the recommended fit from Ocaña et al. inred and from Wang et al. in blue. et al. measured the kinetic rate constants for the formationof HCO + H O to be ( . ± . ) × − cm molecule − s − , in good agreement with the present results, and 0 . + . − . × − cm molecule − s − for HCOOH + H at 300 K.No formation of HCOOH + H was observed in the RPMDtrajectories described below, because of the low reaction prob-ability of this channel and the impossibility to enrich theRPMD statistics. At fixed collisional energy, the only reactiv-ity was obtained in QCT calculations above 300 meV, which isconsistent with a barrier height of 276.0 meV, including ZPE.This reaction follows a direct mechanism, where the OH ap-proaches the H CO in a conformation close to the one in TS2.Indeed, the maximum impact parameter for this process isabout 3 times smaller than the one for HCO + H O at col-lisional energies of 300 meV, decreasing the gap at 1 eV to1.5. As the collisional energy increases so does the reactionprobability, since the system has more energy in the reactioncoordinate. These results are summarized in figure 7.At fixed temperature, reactivity was obtained above 900 K.It is important to remember that only reactivity from reactantsin their ground vibrational states is being taken into account.Our results are in good agreement with the ones calculated byG. de Souza et al. with the canonical variational transitionstate method (CVTST) and are shown in Fig. 8. B. RPMD results
Quantum effects, such as tunneling and zero-point energy,are important at the low temperatures of interest here. Ex-act quantum calculations are not feasible in this system, andone of the most successful alternative methods is Ring Poly-mer Molecular Dynamics (RPMD) , which is used in thiswork. RPMD is a semiclassical formalism based on Path Inte-gral Molecular Dynamics (PIMD), that includes quantum ef-
300 400 500 600 700 800 900 1000 E /meV10 ( E ) / a
400 600 800 1000 E /meV2.02.53.03.54.0 b m a x / a
400 600 800 1000 E /meV01234 P r FIG. 7. QCT results at fixed collisional energy for the formation ofHCOOH + H. On the top left panel, the maximum impact parameter,on the top right panel the reaction probability and in the lower panel,the reactive cross section.
800 900 1000 1100 1200 1300 T /K10 k ( T ) / c m s H CO + OH HCOOH + H
This workG. de Souza et al . FIG. 8. QCT reaction rate constant at constant temperature forHCOOH + H formation. The error bars are two times the standard er-ror of the Monte Carlo integration . In red line the results obtainedby G. de Souza et al. . fects such as ZPE and tunneling . RPMD has been suc-cessfully applied to many reactions, with barriers or deepwells . Here we use the dRPMD program , a direct ver-sion of the RPMD method born as an alternative to the RPM-Drate code to deal with reactions with no barriers. dRPMDis parallelized to allow long propagations at low temperatures,requiring many replicas or beads, N b . The direct RPMD con-sists in two steps, thermalization and real-time dynamics.The thermalization is a constrained PIMD simulation, usingthe Andersen thermostat and freezing the distance betweenthe two reactants at 120 a . This propagation is performedduring 10 -10 steps, depending on the temperature, to war-rant convergence of the initial ZPE of each fragment. At theend of the thermalization the relative position and velocities of7
10 0 10 20 30-6 0 6 12
T = 50 K R z ( B oh r) R x (Bohr)R vector 5 7 920 30 40 50 60 70 R ( B oh r) time (ps) |R|06012018090 93 96 ang l e ( deg r ee ) time (ps) θ R-OH θ CO-OH φ R -10 0 10 20 30-15 -10 -5 0 5 10 15 E coll = 10 meV R z ( B oh r) R x (Bohr)R vector 0 10 20 30 40 50 60 0 10 20 30 40 50 60 70 R ( B oh r) time (ps)|R|060120180 50 52 54 56 58 60 ang l e ( deg r ee ) time (ps) θ R-OH θ CO-OH φ R FIG. 9. Typical trapped RPMD trajectory (left 3 panels) obtained at 50 K. Right panels show a typical QCT trajectory for high (red) initialimpact parameter, at fixed collisional energy of 0.1 meV. In each case, the left panel shows the components R x and R z of the vector, R ,connecting the centers of mass of the reactants along the trajectory. The right bottom panels show the modulus of R as a function of time.The right top panels show the evolution of the angles: between R and OH internuclear vector, θ R − OH , between the CO and OH internuclearvectors, θ CO − OH and the azimuthal angle associated to R , φ . the two reactants are reoriented, imposing a maximum initialimpact parameter, similar to QCT calculations.In the second step, the constrain and the thermostat are re-moved, leading to the classical evolution of a system con-sisting on N atoms × N b particles, which are propagated us-ing a second order simplectic operator. In previous RPMDcalculations , the simplectic propagation is separated in twosteps: first the free polymer and second the real potentialterms. The free polymer propagation was done using a Fouriertransform (FT) , and, being nearly exact, allows a large timestep. In the previous RPMD calculations it was also foundthat, for temperatures below 100 K, many trajectories weretrapped in the complex well in the reactants channel. Thosetrapped trajectories lived for more than 500 ns, and were toolong to be finished. So long-lived complexes, as comparedto those obtained in QCT calculations (of several ps), couldbe explained by the reduction of the available reactants phasespace produced in RPMD, in which ZPE is taken into ac-count. However, it is also known that RPMD may show spuri-ous resonances , which can be attributed to polymer nor-mal mode excitations probably accessed when the frequenciesof these modes are close to the frequencies associated to the“physical” system. These spurious resonances seem to be en-hanced when the free ring polymer is propagated as a separateentity and several alternatives have been proposed . The choice made in this work is not to propagate the freering polymer. Instead, we separate the ring polymer Hamilto-nian in physical kinetic energy and total potential, containingthe PES among the N atoms atoms and the harmonic oscillatorsamong the beads of the same atoms. This implies to reducethe time step to approximately ∆ t / √ N b , with ∆ t the time stepused in the FT propagation. This reduction of the time step in-creases the computational time, but when using parallel com-puters the present factorization reduces the communication toonly neighbor processors, thus allowing a higher speed up ofthe parallelization.The parameters used for the calculation of the RPMD rateconstants for different temperatures are listed in table II. Atthe end of each RPMD trajectory, the product channel is de-termined by analyzing the centroid, similarly to what is donein QCT calculations. What is different in RPMD calculations,is that there is an increasing probability of trapped trajectoriesas temperature decreases below 300 K. Those trapped trajec-tories are formed following an orbit similar to the QCT orbitof Fig. 9 as those shown in previous RPMD calculations inthis system : the long range interaction deviates the initialtrajectories for rather long impact parameters, increasing theend-over-end angular momentum and the rotation of each ofthe two reactants. These trajectories get then trapped orbitingcontinuously keeping the relative orientation of the two reac-8ants fixed, close to the geometry of reactant well, RC1. Whatis different from QCT calculations, is that these trapped orbitslive very long times, longer than t max = 10 ns in this case, whenthey are stopped. Thus we give separate probabilities for di-rect (reactive RPMD trajectories) and trapped trajectories intable II. The rates for direct reaction and trapping processesare then evaluated through equation (16), as listed in table IIThe trapping mechanism becomes significant at tempera-tures below 200 K, where trajectories finished after 10 nswithout dissociation or reaction. This results perfectly mimicthe ones from our previous work as shown in figure 11. Ascompared to our previous PES, the trapping mechanism be-gins at temperatures slightly higher, around 300 K. With re-spect to the direct mechanism rate, the difference found inthe QCT study seems to magnify, what on top of what hasalready been said in the previous section, it is worth notingthat some of the errors may emerge from a poorer statistic,due to the huge computational cost of RPMD trajectories atthese low temperatures. A typical trapped RPMD trajectory isshown in the right panels of Fig. 9, compared with a direct andcomplex-forming QCT trajectories shown in the right panelsof Fig. 9. The first difference in both trajectories is that RPMDtrajectory keeps trapped for much longer times than QCT tra-jectory, that actually manages to break the reactive complex inthe propagation time. Once the two reactants collide, φ R os-cillates quasi periodically in both cases, but the angle θ R − OH vector keeps almost constant. This behaviour is endorsed bythe θ CO − OH angle, and the only appreciable difference is thatQCT amplitudes have larger variations.With respect to our previous work , not only the partitionof the RP Hamiltonian has been varied, but also the time step, ∆ t , has been reduced by a factor of 10 and the number of beadsused in the low temperature calculations has been consider-ably increased, in order to check whether or not long-livedcomplex lifetimes are artificially due to poor convergence.The computational cost of this PES has been considerably re-duced, what has enabled to perform this kind of study. How-ever, the formation of extremely long trapped trajectories per-sist here, as can be seen in figure 11. Therefore, we discardthe possibility of this results being a lack of convergence ofthe calculations. Similar behaviour was also found for otherreactions at low temperature, OH + CH OH , D + H + andH + H + . We therefore conclude that the trapping is a rathergeneral outcome of RPMD calculations at low temperature.An increase of complex lifetimes when including quantumeffects is expected and can be explained according to the sta-tistical theory RRKM by the descent of accesible disso-ciative states in the RPMD calculations with respect to theQCT ones, because ZPE is included. However, the impossi-bility to end the trapped trajectories, dominant at low tempera-tures, seems to indicate that RPMD method fails to reproducethe real lifetime of the reaction complexes. This is attributedto the appearance of spurious resonances as those reportedpreviously , that appear when the frequency of the com-plexes are of the order of those of the ring polymer modes.The normal modes of the free ring polymers are given by ω k = N b k B T / ¯ h sin ( | k | π / N b ) , with k B being the Boltzmannconstant and − N b / ≤ k ≤ N b /
2. The lowest frequency mode, F e r quen cy ( c m - ) Temperature (K)RC1 NMsfree RP w i FIG. 10. RC1 complex normal mode frequencies (blue) and first five ω i free ring-polymer frequencies versus temperature. k =
1, for large N b can be approximated by ω ≈ π T k B / ¯ h ,and when this quantity is of the order of the lower frequencyof the reaction complex, the energy flow between the “physi-cal” and “free ring-polymer” modes may be enhanced, givingrise to spurious resonances, in which the energy is stored inthe free ring-polymer normal modes. The lower ω i are plottedas a function of temperature and compared with the physi-cal frequencies of the H CO-complex in Fig. 10. Clearly at ≈
100 K, the first RP normal mode is of the order of thelower RC1 normal modes. As temperature decreases, thereare Ring-polymer modes in near resonance with several phys-ical modes of the RC1 complex, favoring the energy transfer.Since the density of the RP modes is large, the energy is storedthere for very long times, giving rise to artificial trapping orspurious resonances. T /K10 k ( T ) / c m s PES v1, k dir
PES v1, k trap
This work, k dir
This work, k trap
FIG. 11. Reaction and trapping rates for the present PES, comparedwith the ones from out previous fit.
Several solutions have been proposed to solve the problemof spurious resonances in RPMD calculations of spectra .9 ABLE II. Parameters of the direct RPMD and rate constants. T /K N b N total t max /ns b dirmax / a b trapmax / a P dir P trap k dir / cm s − k trap / cm s − . .
44 0 .
00 0 .
025 0 .
00 1 . × − . . .
82 0 .
00 0 .
048 0 .
00 1 . × − . . .
93 0 .
00 0 .
019 0 .
00 7 . × − . . .
76 0 .
00 0 .
026 0 .
00 6 . × − . . .
93 0 .
00 0 .
016 0 .
00 3 . × − . . .
93 7 .
86 0 .
003 2 . × − . × − . × −
200 384 20000 1 . .
78 13 .
82 0 .
005 2 . × − . × − . × −
100 768 1000 1 . .
95 16 .
60 0 .
003 1 . × − . × − . × −
50 1536 400 1 . .
39 18 .
96 0 .
015 5 . × − . × − . × −
20 3072 143 1 . .
53 21 .
43 0 .
048 6 . × − . × − . × − Among them, the inclusion of a thermostat to the internalmodes of the ring polymer during the dynamics is consid-ered here that may affect the reaction dynamics, because theextra energy of the beads is expected to flow to the physi-cal bonds, producing their fragmentation. A detailed analy-sis needs to be done before applying it to reactive collisionsat low temperatures. Instead, in order to get the total reac-tion rate constant we have to look for an alternative methodto consider the fragmentation ratio of the trapped trajectories:back to reactants (redissociation) and tunneling through reac-tion barriers to products (tunnel). Under this assumption, thetotal reaction rate becomes: k ( T ) = k dir ( T ) + k CF ( T ) (17) k CF ( T ) = k trap ( T ) k tunnel ( T ) k tunnel ( T ) + k rediss ( T ) , (18)where k dir is the rate constant for the direct mechanism and k CF is the product of the trapping rate constant and the ratioof tunneling trajectories towards products. This approach wassuggested before , and the ratios were obtained either fromTST or QCT methods. Here we adopted the QCT reactionprobability as k tunnel ( T ) k tunnel ( T ) + k rediss ( T ) ≈ P QCTr ( T ) , (19)that is, ≈
1% for collision energies below 1 meV. This clas-sical estimate is justified by the small height of the TS, thatcan even be reduced due to anharmonic effects. The resultingRPMD reaction rate constants are shown in figure 12 com-pared with previous results. The behaviour of both results aresimilar and always below the experimental values.Black points and blue triangles represent the RPMD andQCT results obtained in this work, respectively. Red pointsare the RPMD results from our previous PES . Solid linesshow the recommended fit from Ocaña et al. in red and fromWang et al. in blue.The difference between simulated and experimental rateconstants for T >
400 K is attributed to the appearance of othermechanisms, such as non-adiabatic contributions of excitedelectronic state, and/or a lack of accuracy to describe otherchannels such as H + HCOOH product channel. In this workwe focus on the lower temperature behaviour. Here the prob-lem is to determine the tunneling/redissociation fraction of T /K10 k ( T ) / c m s PES v1This workThis work, QCTOcaña et al . fitWang et al . fitWang 15% uncertainty FIG. 12. RPMD reactive rate constant for HCO + H O formation.Black points represent the RPMD results obtained in this work. Bluetriangles represent the QCT results obtained in this work. Red pointsare the results from our previous PES . Solid lines show the rec-ommended fit from Ocaña et al. in red and from Wang et al. inblue. the trapped RPMD trajectories. Such study is fundamentalfor astrochemistry, since it is necessary to use the pure zero-pressure rate constant for the formation of OM in cold envi-ronments.Also, it is important to determine the life-time of the reac-tion complexes to establish the role of complexes (or pressure)in the measured rate constants in CRESU experiments. Workin this direction is now on the way, following the preliminaryarguments already published . VII. CONCLUSIONS
In this work a new full dimensional potential energy surfacehas been developed for the title reaction based on artificialneural networks. The long-range behaviour of the ANN hasbeen analyzed in detail. A source of spurious interactions hasbeen identified due to the non linearity of the transfer func-tion in the ANN, that mixes the internal degrees of freedomof each of the polyatomic fragments. These spurious interac-tions change the energy of the fragments at very long distances10ntroducing artifacts in the dynamics at low energies, thus be-coming inappropriate to study reactivity at low temperatures.To solve this problem, a new factorization of the potentialis proposed consisting in two terms, a diabatic matrix and afull dimensional term, V MB , expanded using neural networks.In the diabatic matrix, the diagonal elements describe eachrearrangement channel, including the potential of each inde-pendent fragment plus their interaction among them, with thelong range interaction properly set. Thus, V MB fits the dif-ference and tends to zero at intermediate distance. This termis further multiplied by a switching function to fully removethe spurious long range interaction introduced by the artificialneural network function. The present potential is more accu-rate than the previous one and also incorporates the channeltowards HCOOH + H product channel.This PES has been used to calculate the reaction rate con-stants using QCT and RPMD methods. It is found that theHCOOH + H products present a near negligible contributionat the energy range considered. The HCO + H O reaction rateconstant presents a non-Arrhenius V-shape as a function oftemperature, dominated by a direct mechanism at high tem-peratures and a complex-forming at low energies, as it wasalso found in the previous PES and in qualitative agreementwith the available experimental data .In spite of the changes done in the calculations, the RPMDresults show an increasing trapping probability as tempera-ture decreases, as reported before . This is attributed to thepresence of spurious resonances occurring when the free ring-polymer normal mode frequencies enter in near resonancewith the low intermolecular normal modes frequencies. It iscrucial to determine the life-time of these complexes and thefragmentation ratio in order to properly determine the zero-pressure reaction rate constant, needed in astrophysical mod-els, and to determine the role of complexes in the measure-ments made in laval expansions. Some work in these direc-tions are now-a-days under way. VIII. SUPPLEMENTARY MATERIAL
See supplementary material for the multi-reference calcu-lations of ground and first electronic states along the MEP,the diabatic matrix, V diab construction, the hyper-parametersused in the NN fits and the normal mode frequencies in thetransition states.A Fortran 90 implementation of this potential energy sur-face is provided in the Supplementary Information. IX. ACKNOWLEDGEMENTS
The research leading to these results has received fund-ing from MICIU (Spain) under grant FIS2017-83473-C2. Wealso acknowledge computing time at Finisterre (CESGA) andMarenostrum (BSC) under RES computational grants ACCT-2019-3-0004 and AECT-2020-1-0003.
X. DATA AVAILABILITY
The data that support the findings of this study are availablefrom the corresponding author upon reasonable request.
REFERENCES R. J. Shannon, M. A. Blitz, A. Goddard, and D. E. Heard,“Accelerated chemistry in the reaction between the hydroxilradical and methanol at interstellar temperatures facilitatedtunnelling,” Nature Chem. , 745 (2013). J. C. Gómez-Martín, R. L. Caravan, M. A. Blitz, D. E.Heard, and J. M. C. Plane, “Low temperature kinetics ofthe CH OH + OH reaction,” J. Phys. Chem. A. , 2693(2014). R. L. Caravan, R. J. Shannon, T. Lewis, M. A. Blitz, andD. E. Heard, “Measurements of rate coefficients for reac-tions of oh with ethanol and propan-2-ol at very low tem-peratures,” J. Phys. Chem. A , 7130 (2015). M. Antiñolo, M. Agúndez, E. Jiménez, B. Ballesteros,A. Canosa, G. E. Dib, J. Albadalejo, and J. Cernicharo,“Reactivity of OH and CH OH between 22 and 64 K: mod-eling the gas phase production of CH O in Barnard 1b,” As-troPhys. J. , 25 (2016). A. J. Ocaña, E. Jiménez, B. Ballesteros, A. Canosa,M. Antiñolo, J. Albadalejo, M. Agúndez, J. Cernicharo,A. Zanchet, P. del Mazo, O. Roncero, and A. Aguado, “Isthe gas phase OH+H CO reaction a source of hco in inter-stellar cold dark clouds? a kinetic, dynamics and modellingstudy,” AstroPhys. J. , 28 (2017). D. E. Heard, “Rapid acceleration of hydrogen atom ab-straction reactions of oh at very low temperatures throughweakly bound complexes and tunneling,” Accounts ofChemical Research , 2620–2627 (2018). A. J. Ocaña, S. Blázquez, A. Potapov, B. Ballesteros,A. Canosa, M. Antiñolo, L. Vereecken, J. Albaladejo, andE. Jiménez, “Gas-phase reactivity of CH OH toward OH atinterstellar temperatures (11.7-177.5 K): Experimental andtheoretical study,” PCCP , 6942 (2019). R. L. Hudson and M. H. Moore, “Laboratory Studies of theFormation of Methanol and Other Organic Molecules byWater+Carbon Monoxide Radiolysis: Relevance to Comets,Icy Satellites, and Interstellar Ices,” Icarus , 451 (1999). N. Watanabe, O. Mouri, A. Nagaoka, T. Chigai, A. Kouchi,and V. Pirronello, “Laboratory simulation of competitionbetween hydrogenation and photolysis in the chemical evo-lution of H O-CO ice mixtures,” AstroPhys. J. , 1001(2007). R. T. Garrod and E. Herbst, “Formation of methyl for-mate and other organic species in the warm-up phase of hotmolecular cores,” A&A , 927–936 (2006). L. E. Snyder, D. Buhl, B. Zuckerman, and P. Palmer, “Mi-crowave detection of interstellar Formaldehyde,” Phys. Rev.Lett. , 679 (1969). J. A. Ball, C. A. Gottlieb, A. E. Lilley, and H. E. Radford,“Detection of Methyl Alcohol in Sagittarius,” AstroPhys. J. , L203 (1970).11 C. A. Gottlieb, P. Palmer, L. J. Rickard, and B. Zuck-erman, “Studies of Interstellar Formamide,” Astrophys. J. , 699–710 (1973). P. D. Godfrey, R. D. Brown, B. J. Robinson, andM. W. Sinclair, “Discovery of Interstellar Methanimine(Formaldimine),” AstroPhys lett. , 119 (1973). A. J. Remijan, J. M. Hollis, L. E. Snyder, P. R. Jewell, andF. J. Lovas, “Methyltriacetylene (CH C H) toward TMC-1:The Largest Detected Symmetric Top,” Astrophys. lett. ,L37–L40 (2006). A. Bacmann, V. Taquet, A. Faure, C. Kahane, and C. Cecca-relli, “Detection of complex organic molecules in a prestel-lar core: a new challenge for astrochemical models,” Astron.Astrophys. , L12 (2012). J. Cernicharo, N. Marcelino, E. Roueff, M. Gerin,A. Jiménez-Escobar, and G. M. Muñoz Caro, “Discoveryof the Methoxy Radical, CH O, toward B1: Dust Grain andGas-phase Chemistry in Cold Dark Clouds,” AstroPhys. J.lett. , L43 (2012). C. Vastel, C. Ceccarelli, B. Lefloch, and R. Bachiller, “TheOrigin of Complex Organic Molecules in Prestellar Cores,”AstroPhys. J. lett. , L2 (2014). S. Cazaux, M. Minissale, F. Dulieu, and S. Hocuk, “Dustas interstellar catalyst. ii. how chemical desorption impactsthe gas,” Astron. Astrophys. , A55 (2016). Y. Oba, T. Tomaru, T. Lamberts, A. Kouchi, and N. Watan-abe, “An infrared measurement of chemical desorption frominterstellar ice analogues,” Nature Astronomy , 228–232(2018). M. Mainitz, C. Anders, and H. M. Urbassek, “Irradiation ofastrophysical ice grains by cosmic-ray ions: a REAX simu-lation study,” A&A , A35 (2016). G. A. Cruz-Diaz, R. Martín-Doménech, G. M. Muñoz-Caro, and Y.-J. Chen, “The negligible photodesorption ofmethanol ice and the active photon-induced desorption ofits irradiation products,” Astron. Astrophys. , 68 (2016). M. Bertin, C. Romanzin, M. Doronin, L. Philippe, P. Jeseck,N. Ligterink, H. Linnartz, X. Michaut, and J.-H. Fillion,“UV photodesorption of methanol in pure and CO-rich ices:desorption rates of the intact molecule and of the photofrag-ments,” AstroPhys. J. Lett. , L12 (2016). W. Siebrand, Z. Smedarchina, E. Martínez-Núñez, andA. Fernández-Ramos, “Methanol dimer formation drasti-cally enhances hydrogen abstraction from methanol by OHat low temperature,” Phys. Chem. Chem. Phys. , 22712(2016). L. G. Gao, J. Zheng, A. Fernández-Ramos, D. G. Truh-lar, and X. Xu, “Kinetics of the methanol reaction with ohat interstellar, atmospheric, and combustion temperatures,”Journal of the American Chemical Society. , 2906–2918(2018). T. L. Nguyen, B. Ruscic, and J. F. Stanton, “A master equa-tion simulation for the oh + ch oh reaction,” The Journal ofChemical Physics , 084105 (2019). A. Zanchet, P. del Mazo, A. Aguado, O. Roncero,E. Jiménez, A. Canosa, M. Agúndez, and J. Cernicharo,“Full dimensional potential energy surface and low tem-perature dynamics of the H CO + OH → HCO + H O re- action,” PCCP , 5415 (2018). O. Roncero, A. Zanchet, and A. Aguado, “Low temperaturereaction dynamics for CH OH + OH collisions on a new fulldimensional potential energy surface,” Phys. Chem. Chem.Phys. , 25951 (2018). P. del Mazo-Sevillano, A. Aguado, E. Jiménez, Y. V.Suleimanov, and O. Roncero, “Quantum roaming in thecomplex forming mechanism of the reactions of OH withformaldehyde and methanol at low temperature and zeropressure: a ring polymer molecular dynamics approach,” J.Phys. Chem. Lett. , 1900 (2019). F. Naumkin, P. del Mazo-Sevillano, A. Aguado,Y. Suleimanov, and O. Roncero, “Zero- and high-pressure mechanisms in the complex forming reactions ofOH with methanol and formaldehyde at low temperatures,”ACS Earth and Space Chem. , 1158 (2019). B. Jiang and H. Guo, “Permutation invariant polynomialneural network approach to fitting potential energy sur-faces,” The Journal of Chemical Physics , 054112(2013). J. Li, B. Jiang, and H. Guo, “Permutation invariant poly-nomial neural network approach to fitting potential energysurfaces. II. Four-atom systems,” The Journal of ChemicalPhysics , 204103 (2013). J. Li and H. Guo, “Communication: An accurate full 15dimensional permutationally invariant potential energy sur-face for the OH + CH → H O + CH reaction,” The Journalof Chemical Physics , 221103 (2015). B. Jiang, J. Li, and H. Guo, “Potential energy surfacesfrom high fidelity fitting of ab initio points: the permutationinvariant polynomial - neural network approach,” Interna-tional Reviews in Physical Chemistry , 479–506 (2016). B. Jiang, J. Li, and H. Guo, “High-Fidelity Potential En-ergy Surfaces for Gas-Phase and Gas–Surface ScatteringProcesses from Machine Learning,” The Journal of Physi-cal Chemistry Letters , 5120–5131 (2020). K. Shao, J. Chen, Z. Zhao, and D. H. Zhang, “Commu-nication: Fitting potential energy surfaces with fundamentalinvariant neural network,” Journal of Chemical Physics ,071101 (2016). X. Lu, K. Shao, B. Fu, X. Wang, and D. H. Zhang, “Anaccurate full-dimensional potential energy surface and qua-siclassical trajectory dynamics of the H + H O two-channelreaction,” Physical Chemistry Chemical Physics , 23095–23105 (2018). G. Knizia, T. B. Adler, and H.-J. Werner, “SimplifiedCCSD(T)-F12 methods: Theory and benchmarks,” TheJournal of Chemical Physics , 054104 (2009). H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby,M. Schütz, et al. , “Molpro, version 2015.1, a package ofab initio programs,” (2015). T. B. Adler, G. Knizia, and H. J. Werner, “A simple and ef-ficient CCSD(T)-F12 approximation,” Journal of ChemicalPhysics (2007), 10.1063/1.2817618. Y. Ajili, K. Hammami, N. E. Jaidane, M. Lanza, Y. N.Kalugina, F. Lique, and M. Hochlaf, “On the accuracy ofexplicitly correlated methods to generate potential energysurfaces for scattering calculations and clustering: applica-12ion to the HCl–He complex,” Phys. Chem. Chem. Phys. ,10062–10070 (2013). J. S. Francisco, “An examination of substituent effects onthe reaction of OH radicals with HXCO (where X=H, F,and Cl),” The Journal of Chemical Physics , 7597–7602(1992). J. R. Alvarez-Idaboy, N. Mora-Diez, R. J. Boyd, andA. Vivier-Bunge, “On the importance of prereactive com-plexes in molecule-radical reactions: Hydrogen abstractionfrom aldehydes by OH,” Journal of the American ChemicalSociety , 2018–2024 (2001). B. D’Anna, V. Bakken, J. Are Beukes, C. J. Nielsen,K. Brudnik, and J. T. Jodkowski, “Experimental and theo-retical studies of gas phase NO and OH radical reactionswith formaldehyde, acetaldehyde and their isotopomers,”Phys. Chem. Chem. Phys. , 1790–1805 (2003). M. Akbar Ali and J. R. Barker, “Comparison of Three Iso-electronic Multiple-Well Reaction Systems: OH + CH O,OH + CH CH , and OH + CH NH,” The Journal of Physi-cal Chemistry A , 7578–7592 (2015). Y. Zhao, B. Wang, H. Li, and L. Wang, “Theoretical stud-ies on the reactions of formaldehyde with OH and OH − ,”Journal of Molecular Structure: THEOCHEM , 155–161 (2007). G. de Souza Machado, E. M. Martins, L. Baptista, and G. F.Bauerfeldt, “Prediction of Rate Coefficients for the H CO +OH → HCO + H O Reaction at Combustion, Atmosphericand Interstellar Medium Conditions,” The journal of physi-cal chemistry. A , 2309–2317 (2020). A. Aguado and M. Paniagua, “A new functional form to ob-tain analytical potentials of triatomic molecules,” J. Chem.Phys. , 1265 (1992). B. J. Braams and J. M. Bowman, “Permutationally invariantpotential energy surfaces in high dimensionality,” Int. Rev.Phys. Chem. , 577 (2009). A. Aguado, C. Suarez, and M. Paniagua, “Accurate globalfit of the h potential energy surface,” J. Chem. Phys. ,4004 (1994). F. Murtagh, “Multilayer perceptrons for classification andregression,” Neurocomputing , 183 – 197 (1991). D. C. Ciresan, U. Meier, J. Masci, L. M. Gambardella,and J. Schmidhuber, “Flexible, high performance convolu-tional neural networks for image classification,” in
Twenty-second international joint conference on artificial intelli-gence (2011). D. K. Duvenaud, D. Maclaurin, J. Iparraguirre, R. Bom-barell, T. Hirzel, A. Aspuru-Guzik, and R. P. Adams, “Con-volutional networks on graphs for learning molecular finger-prints,” Advances in neural information processing systems , 2224–2232 (2015). S. A. King, “Minimal generating sets of non-modular in-variant rings of finite groups,” Journal of Symbolic Compu-tation , 101–109 (2013), 0703035. W. Decker, G.-M. Greuel, G. Pfister, and H. Schönemann,“S
INGULAR (2020). A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury,G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga,A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison,A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai,and S. Chintala, “Pytorch: An imperative style, high-performance deep learning library,” in
Advances in NeuralInformation Processing Systems 32 , edited by H. Wallach,H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, andR. Garnett (Curran Associates, Inc., 2019) pp. 8024–8035. Y. Paukku, K. R. Yang, Z. Varga, and D. G. Truhlar, “Globalab initio ground-state potential energy surface of N ,” TheJournal of Chemical Physics , 044309 (2013). Q. Yu and J. M. Bowman, “Ab Initio Potential for H O + → H + + H O: A Step to a Many-Body Representation of theHydrated Proton?” Journal of Chemical Theory and Com-putation , 5284–5292 (2016). J. Li, Z. Varga, D. G. Truhlar, and H. Guo, “Many-BodyPermutationally Invariant Polynomial Neural Network Po-tential Energy Surface for N ,” Journal of Chemical Theoryand Computation , 4822–4832 (2020). A. Li and H. Guo, “A full-dimensional global poten-tial energy surface of H O + (ã A) for the OH + ( ˜X Σ − ) +H ( ˜X Σ + g ) → H( S) + H O + ( ˜X B ) reaction,” Journal ofPhysical Chemistry A , 11168–11176 (2014). A. Aguado, P. Barragan, R. Prosmiti, G. Delgado-Barrio,P. Villarreal, and O. Roncero, “A new accurate full dimen-sional potential energy surface of H + based on triatomics-in-molecule analytical functional form,” J. Chem. Phys. , 024306 (2010). C. Sanz-Sanz, O. Roncero, M. Paniagua, and A. Aguado,“Full dimensional potential energy surface for the groundstate of H + system based on triatomic-in-molecules formal-ism,” J. Chem. Phys. , 184302 (2013). F. O. Ellison, “A Method of Diatomics in Molecules. I. Gen-eral Theory and Application to H O,” Journal of the Amer-ican Chemical Society , 3540–3544 (1963). F. O. Ellison, N. T. Huff, and J. C. Patel, “A Method of Di-atomics in Molecules. II. H and H + ,” Journal of the Amer-ican Chemical Society , 3544–3547 (1963). A. Dorta-Urra, A. Zanchet, O. Roncero, and A. Aguado,“A comparative study of the Au + H , Au + + H , and Au − + H systems: Potential energy surfaces and dynamics ofreactive collisions,” The Journal of Chemical Physics ,154301 (2015). A. Zanchet, O. Roncero, and N. Bulut, “Quantumand quasi-classical calculations for the S + + H ( v , j ) → SH + ( v (cid:48) , j (cid:48) ) + H reactive collisions,” Phys. Chem. Chem.Phys. , 11391–11400 (2016). P. Ehrenfest, “ XLVIII. Adiabatic invariants and the theoryof quanta ,” The London, Edinburgh, and Dublin Philosoph-ical Magazine and Journal of Science , 500–513 (1917). T. Nagy and G. Lendvay, “Adiabatic Switching Extended toPrepare Semiclassically Quantized Rotational-VibrationalInitial States for Quasiclassical Trajectory Calculations,”Journal of Physical Chemistry Letters , 4621–4626 (2017). M. Karplus, R. N. Porter, and R. D. Sharma, “Exchange Re-actions with Activation Energy. I. Simple Barrier Potentialfor (H, H2),” The Journal of Chemical Physics , 3259–13287 (1965). R. D. Levine and R. B. Bernstein,
Molecular Reaction Dy-namics and Chemical Reactivity (Oxford University Press,1987). S. Wang, D. F. Davidson, and R. K. Hanson, “High temper-ature measurements for the rate constants of C1-C4 aldehy-des with OH in a shock tube,” Proceedings of the Combus-tion Institute , 473–480 (2015). R. A. Yetter, H. Rabitz, F. L. Dryer, R. G. Maki, andR. B. Klemm, “Evaluation of the rate constant for the re-action OH+H CO: Application of modeling and sensitivityanalysis techniques for determination of the product branch-ing ratio,” The Journal of Chemical Physics , 4088–4097(1989). I. R. Craig and D. E. Manolopoulos, “Quantum statisticsand classical mechanics: Real time correlation functionsfrom ring polymer molecular dynamics,” J. Chem. Phys. , 3368 (2004). I. R. Craig and D. E. Manolopoulos, “Chemical reactionrates from ring polymer molecular dynamics,” J. Chem.Phys. , 084106 (2005). I. R. Craig and D. E. Manolopoulos, “A refined ring poly-mer molecular dynamics theory of chemical reaction rates,”J. Chem. Phys. , 034102 (2005). Y. V. Suleimanov, R. Collepardo-Guevara, and D. E.Manolopoulos, “Bimolecular reaction rates from ring poly-mer molecular dynamics: application to H + CH → H +CH ,” J. Chem. Phys. , 044131 (2011). Y. V. Suleimanov, F. J. Aoiz, and H. Guo, “Chemical reac-tion rate coefficients from ring polymer molecular dynam-ics: theory and practical applications,” J. Phys. Chem. A , 8488 (2016). R. Pérez de Tudela, F. J. Aoiz, Y. V. Suleimanov, and D. E.Manolopoulos, “Chemical reaction rates from ring polymermolecular dynamics: zero point energy conservation in Mu+ H → MuH + H,” The Journal of Physical Chemistry Let-ters , 493 (2012). R. Pérez de Tudela, Y. V. Suleimanov, J. O. Richardson,V. Sáez Rábanos, W. H. Green, and F. J. Aoiz, “Stress testfor quantum dynamics approximations: deep tunneling in the muonium exchange reaction D + HMu → DMu + H,”The Journal of Physical Chemistry Letters , 4219 (2014). Y. V. Suleimanov, A. Aguado, S. Gómez-Carrasco, andO. Roncero, “Ring polymer molecular dynamics approachto study the transition between statistical and direct mech-anisms in the H +H + → H + + H reaction,” J. Phys. Chem.Lett. , 2133 (2018). Y. Suleimanov, J. Allen, and W. Green, “RPMDrate: Bi-molecular chemical reaction rates from ring polymer molec-ular dynamics,” Computer Physics Communications ,833 – 840 (2013). H. C. Andersen, “Molecular dynamics simulations at con-stant pressure and/or temperature,” The Journal of ChemicalPhysics , 2384–2393 (1980). A. Witt, S. D. Ivanov, M. Shiga, H. Forbert, and D. Marx,“On the applicability of centroid and ring polymer path in-tegral molecular dynamics for vibrational spectroscopy,” J.Chem. Phys. , 194510 (2009). M. Rossi, M. Ceriotti, and D. E. Manolopoulos, “How toremove the spurious resonances from ring polymer molecu-lar dynamics,” J. Chem. Phys. , 234116 (2014). S. Jang, A. V. Sinitskiy, and G. A. Voth, “Can the ring poly-mer molecular dynamics method be interpreted as real timequantum dynamics,” J. Chem. Phys. , 154103 (2014). R. Korol, N. Bou-Rabee, and T. F. M. III, “Cayley modi-fications for strongly stable path-integral and ring-polymermolecular dynamics,” J. Chem. Phys. , 124103 (2019). N. Bulut, A. Aguado, C. Sanz-Sanz, and O. Roncero,“Quantum effects on the D + H + → H D + deuteration re-action and isotopic variants,” J. Phys. Chem. A , 8766(2019). R. A. Marcus, “Lifetimes of active molecules. i,” The Jour-nal of Chemical Physics , 352–354 (1952). R. A. Marcus, “Lifetimes of active molecules. ii,” The Jour-nal of Chemical Physics , 355–359 (1952). J. O. Richardson and S. C. Althorpe, “Ring-polymer molec-ular dynamics rate-theory in the deep-tunneling regime:Connection with semiclassical instanton theory,” The Jour-nal of Chemical Physics131