{'fraction_non_alphanumeric': 0.023187052598817304, 'fraction_numerical': 0.007236227824463119, 'mean_word_length': 4.907837278786486, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Instanton rate constant calculations using interpolated potential energy surfaces in non-redundant, rotationally and translationally invariant coordinates

A trivial flaw in the utilization of artificial neural networks in interpolating chemical potential energy surfaces (PES) whose descriptors are Cartesian coordinates is their dependence on simple translations and rotations of the molecule under consideration. A different set of descriptors can be chosen to circumvent this problem, internuclear distances, inverse internuclear distances or z-matrix coordinates are three such descriptors. The objective is to use an interpolated PES in instanton rate constant calculations, hence information on the energy, gradient and Hessian is required at coordinates in the vicinity of the tunneling path. Instanton theory relies on smoothly fitted Hessians, therefore we use energy, gradients and Hessians in the training procedure. A major challenge is presented in the proper back-transformation of the output gradients and Hessians from internal coordinates to Cartesian coordinates. We perform comparisons between our method, a previous approach and on-the-fly rate constant calcuations on the hydrogen abstraction from methanol and on the hydrogen addition to isocyanic acid.


Introduction
In spite of neural networks gaining ever more attention for their ability to interpolate potential energy surfaces (PES) at levels of accuracy approaching the highest levels of theory,  simple improvements can be made to enhance their robustness and range of applications, especially for the purpose of instanton rate calculations, an area where little work has been conducted to date, though with some exceptions. 28,29 External degrees of freedom, i.e. the rotation or translation of a molecule do nothing to alter the internal energy of the molecule, neither does the exchange of the same atoms within a molecule. Neural networks may be trained using a variety of different descriptors, for example Cartesian coordinates , internuclear distances, z-matrix coordinates or projections on to normal mode displacements. 28 Though this approach produces highly accurate results, it is more cumbersome to account for geometrical symmetries. Ideally, neural networks would be constructed or adapted to account for these symmetries, thus allowing for more user-friendly interfacing with quantum chemistry programs. In this contribution we briefly introduce instanton theory in section 2.1 and thereafter focus on how to account for the geometrical symmetries when using neural networks in predicting reaction rates using instanton theory. Like-atom exchange symmetry may, in specific cases, be treated explicitly, for example using symmetry adapted internals, 30,31 but is in general better treated using atomistic neural networks. 32,33 Because we will use the interpolated PES to determine rate constants in instanton theory, information on the gradient and the Hessian of the surface is also required. To train the neural network, we use data obtained from ab initio methods. The coordinates, gradients and Hessians of the training data are then converted into the desired internal coordinate system.
A test set is also produced using the same approach. Our program offers three different input descriptors: internuclear distances, inverse internuclear distances and z-matrix coordinates.
Usage of internuclear distance based coordinates as descriptors is well established in neural networks for small molecules, 6,12,19,[34][35][36] while z-matrix coordinates, or some combination of internuclear distances, angles and dihedrals, also have some pedigree, particularly in atomic chain neural networks. 2,37,38 After the training is complete, the back-end program, which will use the trained network, reconverts the internal gradients and Hessians back into their Cartesian counterparts.
A number of significant differences exist between neural networks trained using the descriptors mentioned in comparison to those trained on Cartesian coordinates. The training can either be done by minimising the error in Cartesian derivatives or in derivatives with respect to internal coordinates. The latter is more natural, because the input data are converted into internal coordinates. However, the final interpolated PES is used in Cartesian coordinates. Therefore, the Cartesian derivatives must be as accurate as possible in the end. The transformation of gradients and Hessians from Cartesian coordinates to internals introduces a certain amount of numerical noise to the training data, or conversely, to the back-converted gradients and Hessians from their interpolated internal counterparts. This noise comes as a consequence of using pseudoinverses in the transformation process which is necessary because in general the transformation matrix is highly singular, but also because usually there are more Cartesian components than non-redundant internal coordinate components. The noise can be suppressed to acceptably low levels, the details of which we outline in section 2.3. The next major difference is the choice of which internal descriptors to use. Z-matrix coordinates require the user to determine a build list, usage of internuclear distance based descriptors usually means disregarding a large number of potential internuclear distances in order to ensure the set chosen is non-redundant. A non-redundant set of descriptors is desirable though not necessary. The minimum number of descriptors is given by min ( where N is the number of atoms in the molecule and N 0 are the number of rotational and translational degrees of freedom. Using a larger number of descriptors than this minimum will not result in a poorer interpolation but will increase the time required to train the neural network.
The last major difference is the minimisation scheme of the residual with respect to the weights and biases. The weights and biases are exported to an external program that requires the gradients and Hessians be given in Cartesian coordinates. It is important that the weights and biases are optimised in such a manner as to ensure that the back-transformation results in an as small a precision loss as possible. We describe a method to achieve this by making a few modifications to the traditional residual minimisation scheme.
The PESs to be interpolated are taken from the following reactions: hydrogen abstraction from methanol and the hydrogen addition to isocyanic acid. These reactions are of interest to the astrochemical community [39][40][41][42] and are good candidates for interpolation by neural networks due to the relatively low number of atoms involved.
The paper is organised as follows. We describe in detail the architecture of the neural network used and the procedure of transforming the ab initio training data to the desired internal coordinate. Next we show how we choose a non-redundant set of the internal coordinates and then explain the residual minimisation procedure. Finally, we export the interpolated surface to DL-FIND 43 and perform the required rate constant calculations using instanton theory and compare with on-the-fly calculations.

Instanton theory
Before expounding on the particulars of the theory of neural networks, it is necessary to explain why we need an interpolator capable of predicting energies, gradients and Hessians.
Our rate constants are calculated using instanton theory, for the finer details thereof we refer the reader to the literature [44][45][46][47][48][49][50][51][52][53] and simply reproduce the equation for the rate constant where there are N degrees of freedom, N 0 translational and rotational degrees of freedom, P discretisation points of the Feynman path (images), β is the inverse temperature (β = where q a k is the mass-weighted coordinate component a of image k. It is worth noting that transforming to mass-weighted coordinates is a simple step in Cartesian coordinates, hence the reason why DL-FIND uses this coordinate system and why modifying it to accept Hessians and gradients in any number of different internal coordinate systems is a far more strenuous task than modifying the output of the neural network to conform with this preexisting structure. The last term in equation (2) contains the second derivative of the energy (Hessian) along the instanton path, hence the reason why we seek to interpolate the PES up to the second order spatial derivative. The gradient of the PES does not appear in either of the above equations, but it too is a necessary part of the procedure to locate the instanton. The most probable instanton path is defined such that δS E [y] = 0, in other words, it is a path of stationary Euclidean action, localising this path requires the gradient of the PES. 54,55

Network architecture
The architecture used in the neural network presented here combines two multi-layer perceptrons, one with a single hidden layer (1HL) and one with two hidden layers (2HL). The outputs from each network's final layer connect to a single node representing the target energy for the input geometry.
When speaking of the descriptors, we refer to the internal coordinates fed to the 2HL network. A different, and twofold larger, set of descriptors are sent to the 1HL network. If we let x be an I dimensional vector representing the descriptors of the 2HL network, then x −1 i and x 2 i are the i th and i + I th descriptors of the 1HL network respectively. When referencing the 1HL network in formulas, roman upright script will be used.
The respective formulas for the j, k, and p th nodes of the hidden layers y (1) j , y (2) k , y (3) p and the energy are given by The symbols w and b are the usual weights and biases connecting the layers and f are the transfer functions. A useful notation for this type of architecture is The superscripts on the weights, biases and transfer functions represent their location in the neural network, whereas the superscripts on the descriptors x are powers. As can be seen, the 2HL network has J and K nodes in the first and second hidden layer respectively, while the 1HL network has 2I descriptor nodes, and P hidden nodes. The transfer functions f (1) and f (2) are chosen as f (1) is chosen as The gradients and Hessians are determined from equation (6) via the chain rule with respect to the descriptors. With this in mind, the cost function, or residual, used here accounts for the inaccuracies in the energies, gradients, Hessians and Hessian eigenvalues and is given by where lower-case greek symbols represent the interpolated values obtained by the neural network and upper-case greek are the reference values in the training set. It is assumed that the reference quantities in the training and test sets have been provided in Cartesian coordinates. The sum of all the elements squared is used as the norm for matrices, equivalent to vectors. A set of weights has also been included {A E , A G , A H , A L } should there be a need to emphasize accuracy in one or more parameters against any other. Furthermore, energies, gradients and Hessians are usually generated for one common molecular configuration, hence The reason for the unusual choice of network architecture and usage of differing powers of the descriptors is to explicitly provide to the network some of the geometrical dependencies that influence the properties we wish to interpolate, rather than relying solely on the weights and biases to intuit these dependencies. For instance, the energy of the system exhibits some dependence on the inverse of the internuclear distance, this property can be more easily an internal residualR (a tilde will be used henceforth to signify any quantities in the internal coordinate system). One would then minimize this residual with respect to the weights and biases and, once a convergence criteria has been reached, export the optimised weights and biases to the external program where the internal quantities Γ ,H,Λ are back-transformed into their Cartesian counterparts.
The alternative is of course to leave the reference quantities in their Cartesian representation, back-transform the interpolated quantities γ,η,λ → {γ, η, λ} and obtain a residual in Cartesian coordinates. It is in fact this approach which is adopted here due to the fact that the training of the neural network, based on a Cartesian residual, is able to actively adapt the weights and biases to counteract uncertainties introduced by the coordinate transform, we give an in depth explanation of this property in the following section .

Coordinate transformation and training of the neural network
We use the following algorithm for coordinate, gradient, and Hessian transformation: 56 with B, Wilson's B matrix, is in general non-square and potentially highly singular, therefore P = I. In the reactions tested in section 3, P ∼ I and therefore is set equal to I, we include it here and in the following formulas for completeness. Owing to the nature of B, one utilises pseudoinverses, indicated by † , to derive the reference gradients and Hessians in internal coordinates. The Hessians are transformed as where x a and x b refer to the Cartesian coordinates of atoms a and b, respectively. Hence subscript i maps to a unique pair of atoms a and b.
The nature of the pseudoinverse means that some accuracy is lost in performing the transformation, applying the back-transformation may produce results which are significantly different compared with the original gradient or Hessian. This phenomenon is governed by the condition number κ of B, a quantity proportional to the product B · B † .
A high condition number would mean, for instance, that in spite of obtaining a quite accurate approximationγ ∼Γ after training, the magnitude of uncertainty δγ in the internal representation, upon back-transformation, is magnified κ-fold in the uncertainty δγ ∼ κδγ of the Cartesian representation. In this work, we make use of Tikhonov regularisation 57 in creating pseudoinverses.
In the case of an internal residual, this addition of uncertainty to the interpolated Cartesian quantities occurs after the training phase, there is no way for the neural network to actively compensate for this during the training phase. If however a Cartesian residual is used, each step of the training phase, i.e. the minimum search on the residual hypersurface, is actively guiding the weights and biases to correct for the error incurred during transformation.
A more insightful form for the residual, where the interpolated, internal quantities are transformed back to Cartesian coordinates, has the following form and X is the matrix of eigenvectors of η. From equation (14) it is relatively straightforward to construct ∂R/∂ω i and further proceed with the optimisation of the residual. The reason for the replacement of λ in equation (14) is for the sake of the numerical stability of the derivative of R with respect to the weights and biases. If we stack all the weights and biases in a vector ω, take the derivative of η with respect to ω i and diagonalise, it is clear that both η and ∂η/∂ω i must share the same set of eigenvectors since the derivative with respect to the i th weight/bias ω i has no dependence on any spatial coordinate. This allows one to use ∂η/∂ω i in both the derivative of the third and fourth sums of equation (14) ∂R where • refers to the Hadamard product. Equation (15) is then used to minimise the residual by means of the L-BFGS algorithm. 58 To obtain ∂λ l /∂ω i another possibility is to directly diagonalise ∂η l /∂ω i for every i. However this approach is encumbered by larger numerical instability, and increased computational expenditure than taking ∂λ l /∂ω i = diag X T ∂η l ∂ω i X .

Finding non-redundant descriptors
When working with internuclear distances as descriptors, for the purpose of finding a non- The approach we have adopted is firstly to examine the total number C of non-redundant combinations of internuclear distances. If this number is less than a certain threshold T (10 7 is a reasonable figure), all combinations will be tested. The formula for C is If C > T we randomly choose a combination of 3N − N 0 internuclear distances from the available N (N − 1)/2 a total of T times, importantly however, only those combinations which join each atom to at least three other atoms are accepted. Once the set of internuclear distance combinations is established, the next step is to examine the sum of all internuclear distances for this particular combination and average this sum over all the configurations in the test and training set. Denoting this quantityS, for a particular combination k of internuclear distances, it isS k 1≤k≤min(C,T ) where it is assumed that N > 4 and r is an internuclear distance and ν = N train + N test .
The combination producing the smallestS k is that which is then used as our descriptors because the change in distances, which are short to start with, can be assumed to have the largest influence on the total energy of a molecule. Other approaches exist for finding the descriptors, a quasi non-redundant set of descriptors can be found by choosing more than the minimum 3N − N 0 internuclear distances. Another approach might be to build Wilson's B matrix for each internuclear distance combination k and molecular configuration j, yielding a set of singular values, the variance of which, averaged over j, is as small as possible. To that end, it is also possible to select the combination k that minimises the average variance inγ or η. In section 3 we present results based on the smallestS k approach as, for the cases tested, this gave fast convergence of the residual and produced B matrices which are reasonably non-singular. Moreover, in the cases tested, C is maximal for the hydrogen abstraction from methanol where it is only 54264, hence we test all possible combinations to find the one producing the smallestS k .

Applications
We have applied the neural network described in section 2 to the hydrogen addition to isocyanic acid H + HCNO → H 2 CNO (R1) and to the hydrogen abstraction from methanol The reaction rate constant for R1 is shown in Fig. 1. The black curve, representing rate constants calculated on the interpolated neural network PES is never more than one half of an order of magnitude away from the rate constants calculated on the fly (red curve).
Since the rate constant depends strongly on many details of the potential energy surface along the instanton path, a higher interpolation accuracy is difficult to achieve. We expect the intrinsic errors of the underlying quantum chemical methods as well as the semiclassical approximation inherent to instanton theory to lead to errors in a similar range. 63 The main source of deviation in the rate constants from the neural network for R1 comes from the fact that S contains more than the expected 6 modes with a non-negligible projection on the translation and rotation modes and hence more than the expected 6 zero (or nearly zero) eigenvalues. This makes the identification of which eigenvalues to exclude in equation (1)  The rate constant of a single instanton optimised on-the-fly at 65 K is shown in Fig. 2.
The network architecture for R2 is 15-30-30→1←60-30 and the weight parameters were The accuracy for R2, compared with the red curve, for which elongations along normal modes have been used as descriptors, 28 is very favourable. For all temperatures the deviation is less than the half an order of magnitude. We note that the red curve is the result of an average of multiple neural network calculations, while our result required only a single network to be trained. More training runs with different starting weights lead to similar rate constants.
In both reactions there seems to be very little noise in the rate constant curve, indicating that the Hessians are also quite noise-free. This is, of course, the desired outcome and shows the usefulness of equation (15) in the training procedure, and the tendency of Tikhonov regularisation to minimise uncertainties generated as part of the pseudo-inversion operation.

Discussion
Perhaps the most important aspect of any machine-learning-based approach to modelling (not simply in quantum chemistry) is to provide high quality training data. It is not difficult to imagine how error-prone it could be when attempting to create a training database for a PES if one is required to first orient and displace the molecule of interest in precisely the form prescribed by the neural network. The main objective of this contribution is to prove that a neural network interpolator, capable of accurately predicting energies, gradients and Hessians, can be created that is independent of trivial rotations and translations of the molecule being modelled. As the results in section 3 have shown, this objective is indeed achievable.
Another goal for PESs interpolated by neural networks is for the set up to be as free of user-input as possible. For the network shown here, this is only partially the case, a number of parameters must be chosen by the user such as network architecture, layer sizes, residual weights {A E , A G , A H , A L }, transfer functions, pseudoinversion scheme and internal coordinate system. The system sizes used here are small enough that one may make a number of reasonable guesses for these parameters, initiate a training run for each guess, and simply export the weights and biases for that run which gives the smallest residual. This does not, however, transfer easily to reactions involving ten or more atoms as the time required for each training step increases sharply with system size.
For larger systems, the first computational hurdle is reached in finding a non-redundant set of internuclear distances (cf. equation (16)). As the system size increases, a smaller fraction of the number of potential non-redundant combinations can be tested, increasing the chance that more favourable combinations are not utilized, hence increasing the number of steps required in the training phase before an acceptably small residual is reached, but for systems of less than or approximately 10 atoms it is not too problematic to simply test all possible non-redundant internuclear distances.
The second computational restriction is the increase in the size of the layers as the system size increases. A good rule of thumb for the extent of the hidden layers is to have twice or three times as many hidden nodes as input nodes, if we apply this to a system    It is clear then that naively applying deep neural networks to large systems is impractical at this stage. Importantly, however, the architecture of our system may make this issue

The residual and usage of eigenvalues therein
In order to give the reader some sense of the learning rate of the neural network, figure 3 shows the reduction of the test set residual in equation (14) with each optimization step. It is noteworthy that the residual for both reactions stabilizes, given that no regularization was used as one typically expects the test-set residual to begin increasing at some point while the training set residual continues to decrease. As can be seen from equation (1)  elements require greater emphasis on accuracy. This extra optimisation option may not always be necessary (as was the case with R2), but as we found with R1, it was necessary to set A L = 2, this had a self-reinforcing, positive influence on the overall residual and on the accuracy of the resulting Hessians and eigenvalues, a result not otherwise possible.

Other coordinate systems
It is worth highlighting the fact that while the procedure for obtaining non-redundant coordinates outlined in section 2.4 is appropriate and manageable for the chemical systems treated here, the combinatorial scaling of equation (16), should one wish to test all possible combinations, is prohibitive. We therefore acknowledge that the procedure outlined in 29,64 obviates this problem, at the expense of complicating somewhat the expressions for the coordinate, gradient and Hessian transformations in equation (10) and the subsequent equations for the residual.
Finally, we note that the results presented in section 3 were produced using the inverse internuclear distance coordinate system. The other coordinate systems mentioned, internuclear distances and z-matrix coordinates, were also tested, yet neither coordinate system was able to yield a residual as small as that obtained using inverse internuclear distances. As a guide, Fig. 4 shows the performance of internuclear distances as descriptors applied to R2, the rate constants calculated are overestimated by slightly more than one order of magnitude for higher temperatures. Moreover, the instanton search only reaches convergence at higher temperatures, below 210 K, no valid instanton could be found.

Conclusions
A neural network PES interpolation scheme has been presented capable of calculating energies, gradients and Hessians and which is intrinsically independent of the rotation and translation of the molecule being modelled. Our design can be readily applied to small molecules and is able to produce results which are comparable in accuracy to other established methods. The architecture of the neural network used here allows for a broad scope of refinements in the future with the aim of increasing the range of applicability of deep neural For this reaction, the rate constants are overestimed and instantons can only be located at higher temperatures.
networks to reactions involving more than approximately ten atoms.

Acknowledgements
This work was financially supported by the European Union's Horizon 2020 research and in-