Exhaustive state-to-state cross sections for reactive molecular collisions from importance sampling simulation and a neural network representation

There are numerous situations in physical chemistry which involve a potentially large number of states and transitions between them. Examples include complete line lists for polyatomic molecules in hot environments [e.g., high-resolution transmission molecular absorption (HITRAN)¹ database] or state-to-state (STS) cross sections in reactive and nonreactive molecular collisions. Exhaustively probing and enumerating all relevant combinations or creating high-dimensional analytical representations is usually impossible. On the other hand, it has been shown for spectroscopic applications that omission of certain crucial states makes prediction or modeling of the spectroscopic band difficult or even impossible.² Another example is hypersonic flow around a space vehicle reentering the atmosphere. Temperatures can easily reach 20 000 K, for which reliable experimental data are sparse and the energy distributions are out of equilibrium. In such an environment, the space vehicle is exposed to a collisionally dense environment which generates an immensely diverse population of rovibrational states between which collisions take place.³ An accurate, fast, and reliable method is required to include this information in more coarse grained models to study the associated dynamics. The question thus is how to best probe and represent a function with >10⁷ values for discrete input data, without explicitly computing each entry which may require thousands to millions of samples to statistically converge each of the entries.

The present work attempts to develop such a model for state-to-state cross sections σ_{$v$⁠,j→$v$′j′}(E_t) between initial (⁠$v$⁠, j) and final (⁠$v$′, j′) rovibrational states at given relative translational energy E_t. For this, the N(⁴S) + NO(²Π)(⁠$v$⁠, j) → O(³P) + N₂(⁠$X1Σg+$⁠)(⁠$v$′, j′) reaction is considered because (a) it is relevant in the hypersonic flight regime characterized by temperatures T ≤ 20 000 K at which a multitude of the available rovibrational states are occupied and accessible and (b) an accurate, fully dimensional, and reactive potential energy surface (PES) is available.⁴ Specifically, N₂-formation rates from simulations on the ³A′ and ³A″ PESs are in favorable agreement with experiments for temperatures T ≥ 2000 K and N₂-formation below 5000 K is dominated by processes on the ³A″ surface.⁴ State-to-state cross sections are typically required when modeling the reactive flow around a re-entry object with macroscopic dimensions using techniques such as direct simulation Monte Carlo (DSMC)⁵ as the flows are usually not in thermal equilibrium. However, it should be noted that the present ansatz will be applicable to other relevant simulations and hypersonic flow is merely chosen because the occupation number of the available states is high and the number of states is large as well.

Scattering cross sections can be determined from quasiclassical trajectory (QCT) simulations.^6–8 Here, Hamilton’s equations of motion are solved using a fourth order Runge-Kutta numerical integration in reactant Jacobi coordinates. The analytical representation of the ³A′ potential energy surface (PES) is a reproducing kernel Hilbert space (RKHS)^9,10 based on MRCI+Q/aug-cc-pVTZ calculations. The two possible reactive channels are (I) nitrogen exchange (NO + N′ → N′O + N) and (II) N₂ molecule formation (NO + N → N₂ + O). Besides these two channels, the elastic or inelastic collisions can also lead to kinetically or internally excited reactants.

For the ³A′ surface of N₂O, a maximum of 47 and 57 vibrational states are available for NO and N₂, respectively, while the maximum rotational quantum numbers for NO and N₂ are 241 and 273. Altogether, there are 6329 rovibrational states for the N + NO channel and 8733 states for the O + N₂ channel. For one defined initial state (⁠$v$⁠, j), there are >10⁴ possible final states (⁠$v$′, j′). To converge each of the σ_{$v$⁠,j→$v$′j′}(E_t) for given E_t with a statistical significance of ∼10%, an estimated ≥10¹³ trajectories would be required (10⁷ trajectories to converge one cross section, see below; ∼10⁴ initial states; ∼10⁴ final states per initial state). Hence, using QCT simulations to directly sample all possible rovibrational initial states is computationally impractical, even for this simple 3-body system which also ignores the complexity of electronic states.¹¹ Thus, alternative approaches need to be explored.

To better define the QCT sampling problem, state-to-state cross sections for σ_{$v$=6,j=30→$v$′,j′}(E_t = 2.5 eV) were considered. A total of 2 × 10⁷ trajectories were run initially, which is considered as the reference. Out of the 3784 energetically accessible states, 3420 final states are found as products, i.e., 90.4%. Compared to this, running fewer trajectories (8 × 10⁴, 1.6 × 10⁵, 1.6 × 10⁶, and 9.6 × 10⁶) finds 37%, 54%, 83%, and 90% of the final states, respectively, see Fig. S1, which converges at approximately the cube root of the number of trajectories. In addition, running too few trajectories leads to highly oscillatory cross sections due to the limited statistics of the final state.

Some computational gain can be obtained from exploring the fact that cross sections often vary smoothly with $v$ and j. This allows us to locally average computed cross sections according to

where n_$v$ and n_j are the bin widths for vibration and rotation over which the state to state cross sections are averaged. Figures 1(a) and 1(b) report the raw data from 1.6 × 10⁶ trajectories and locally averaged cross sections from the same data set with n_$v$ = 2 and n_j = 3, respectively. Comparison with the unaveraged result from N_tot = 2 × 10⁷ trajectories in Fig. 1(d) illustrates the benefit of local averaging; see also Fig. S4, for example (⁠$v$ = 10, j = 60 and E_t = 1.8 eV). Convergence of the state-to-state cross sections with respect to N_tot and the effect of local averaging are shown in Figs. S2 and S3. Averaging over neighboring states reduces the noise and also decreases the number of trajectories required approximately by the square root of the number of states averaged over, i.e., $(2nv+1)(2nj+1)=35$⁠. Besides that, local averaging also provides values for σ_{$v$⁠,j→$v$′,j′} > 0 for unsampled transitions and reduces local oscillations while correctly describing the broad features (see Fig. 1 and Fig. S4). However, sharp resonances, for which the width is comparable to the size of the averaging window, are washed out in this approach, and depending on the application, this point would need to be considered separately.

Concerning the sampling strategy for choosing initial conditions, it is worthwhile to note that large impact parameters b mostly lead to nonreactive collisions. A straightforward approach to sampling b is to draw it from $b2/bmax2$ with 0 ≤ b ≤ b_max (see Fig. S5) and the number of trajectories sampled in the interval b + db increases with increasing b, which leads to a larger fraction of nonreactive trajectories when b increases. For such situations, importance sampling (IS)¹² can provide a more advantageous protocol as those values of b for which reactive trajectories are more likely to occur are chosen with higher probability, which causes the cross section of all of the exit channels to converge at the same rate.

In order to determine the necessary weighting function $w$(b) for a particular trajectory, the following strategy is used. For given (⁠$v$⁠, j, E_t), first a few thousand trajectories are run by uniformly sampling 0 ≤ b ≤ b_max (see Fig. S5). The number of reactive trajectories as a function of b is fitted to

From this, the distribution of the cross sections, g(b) = 2πbn(b), is determined. Next, 1.6 × 10⁵ initial conditions are sampled from g(b) (see Fig. S5), the trajectories are explicitly run, and the weight of each trajectory is calculated according to $w$(b) = f(b)/g(b), where $f(b)=2b/bmax2$⁠.

The performance of IS is illustrated in Fig. 1(c) and Fig. S4(C) which report the locally averaged cross sections from 1.6 × 10⁵ samples for σ_{$v$=6,j=30→$v$′,j′}(E_t = 2.5 eV) and σ_{$v$=10,j=60→$v$′,j′}(E_t = 1.8 eV), respectively. For these two examples, 42 193 and 38 665 reactive trajectories are found using IS from N_tot = 1.6 × 10⁵ compared with 44 737 and 44 979 reactive trajectories for N_tot = 1.6 × 10⁶ from conventional sampling, respectively. This leads to an efficiency increase by one order of magnitude when IS is used. The effect of IS on the convergence of state-to-state cross section can also be seen in Figs. S2 and S3.

Overall, IS and local averaging lead to an estimated reduction of the required number of QCT trajectories by a factor of ∼60 which will be explored next to cover the entire state space for state-to-state cross sections for different selected initial states. Those will then be used for training a neural network (NN) to construct a model to compute state-to-state cross sections.

To compute the necessary reference data to train the NN, 10 initial $v$-states (⁠$v$ = 0, 3, 6, 9, 12, 15, 19, 23, 28, and 34) and 12 initial j−states (j = 0, 25, 50, 75 100, 125, 145, 160, 175, 190, 200, and 210) were sampled. The relative translational energies (E_t) for the N + NO collision were 0.05 0.1, 0.25, 0.5, 0.8, 1.2, 1.6, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, and 5.5 eV. Thus, for a total of 1232 initial conditions in the (⁠$v$⁠, j, E_t) space, for each of them, 1.6 × 10⁵ QCT trajectories were run with IS of b. The cross sections from these ∼10⁸ trajectories constitute the training set for learning the NN to predict σ_{$v$⁠,j→$v$′j′}(E_t) for any initial (⁠$v$⁠, j, E_t).

The network architecture used here is inspired by ResNet.¹³ The NN transforms its input through four identical residual blocks,¹³ after which a linear transformation followed by a scaled sigmoid function is used to obtain the final output; see Fig. S6. Before they can be used as input to the residual blocks, the raw input features $f∈RNin$ are first transformed by a linear transformation x = Wf + b to a vector $x∈RF$ with the same dimensionality F as the hidden features in the residual blocks in order to simplify the formulation of skip connections.¹³ Here, $W∈RNin×F$ (weight matrix) and $b∈RF$ (bias vector) are parameters to be optimized. The residual blocks consist of two dense layers with the same number of nodes F (see Fig. S6) and transform their input x^l according to

where the superscript l denotes parameters or feature representations of layer l and $W∈RF×F$ and $b∈RF$ are parameters. Two different activation functions, one for rectified linear units (ReLU)¹⁴ and a self-normalizing¹⁵ inverse hyperbolic sine (snasinh),¹⁶ are used in the residual blocks. Use of residual blocks improves the flow of gradients between layers^13,17 and helps alleviate the “vanishing gradients problem.”¹⁸ After the last residual block, the final output is obtained from

where C is a scaling constant, $sig(x)=(1+e−x)−1$ denotes the sigmoid function, and the superscripts o and l denote parameters $Wo∈RF×1$ and $bo∈R$ corresponding to the output layer and the hidden features x^l obtained after the last residual block, respectively.

The initial states of the reactants are described by N_in = 12 input features f, namely: (i) internal energy, (ii) vibrational energy, (iii) vibrational quantum number, (iv) rotational energy, (v) rotational quantum number, (vi) angular momentum of the diatom, (vii) relative translational energy, (viii) relative velocity, [(ix) and (x)] turning periods of the diatom, (xi) rotational barrier height, and (xii) vibrational time period of the diatom. For state-to-state cross sections, the same 12 features for the final states of the products are also included as input (i.e., N_in = 24). The loss functions (L_f) were defined as

where y′ and y are the reference (QCT) and predicted values (NN), respectively. This loss function penalizes all relative errors in the prediction of cross sections approximately equally irrespective of the absolute magnitude of the reference value. This is important because for small values of the cross section, a small absolute error corresponds to a large relative error. All parameters of the NN are initialized according to the Glorot initialization scheme¹⁸ and optimized using Adam optimization¹⁹ with an exponentially decaying learning rate. From all state-to-state cross sections (∼8 × 10⁶), different numbers (N_train) of training samples were randomly chosen for training and N_train/10 data were used for validation. From the remaining data, 2 × 10⁵ values constitute the test data set. The convergence of the model with training set size is reported in Fig. S7. Training the NN for state-to-state cross sections (referred to as NN-STS) was done with N_in = 24, F = 24, and C = 0.4.

Since trajectories were calculated for different initial reactant states (⁠$v$⁠, j, E_t), the total cross sections (σ_tot(⁠$v$⁠, j, E_t)) for the N + NO → N₂ + O reaction are also available. On these data, another NN (with N_in = 12, F = 8, and C = 85.0) was trained using 1122 training data and 110 validation data to get a model for σ_tot(⁠$v$⁠, j, E_t) (referred to as NN-Tot). Both networks were trained using TensorFlow.²⁰ The models with the lowest validation losses [Eq. (5)] are selected and used to predict different observables in the following and compared with values from explicit QCT simulations.

The entire reference state-to-state data set (∼10⁷) has an average cross section of 0.00246 $a02$⁠. To test the quality of model NN-STS, 2 × 10⁵ state-to-state cross sections were randomly chosen from the test data. The root-mean-square error (RMSE) is 0.000 328 $a02$ for those data with a maximum deviation of 0.008 $a02$ (for σ_QCT = 0.211 $a02$⁠, σ_NN−STS = 0.219 $a02$⁠), and the correlation coefficient is R² = 0.993; see Fig. S8. A similar analysis was carried out for model NN-Tot which has an average cross section of 16.18 $a02$ for all 1232 data points. The RMSE is found to be 0.1552 $a02$⁠, and the correlation between NN-Tot and QCT is R² = 0.9997, see Fig. S9, indicative of the high quality of the fit.

To quantify the accuracy of the NN, additional QCT calculations were performed for independent initial conditions at fixed E_t. Total QCT cross sections are then compared with the NN predictions; see Figs. 2 and 3 and Table S1. Again, the NN results describe the explicit QCT simulations which validate the use of such a model to predict microscopic information for such a reaction.

According to the method described in Ref. 11 (see also the supplementary material), initial state specific rates for the reaction were also calculated at temperatures between 2000 and 20 000 K for a few selected reactant states using QCT simulations and compared with the rates obtained from the NN models in Fig. 4 and Fig. S10. The NN models successfully capture the trends as well as the magnitudes of the rates from QCT. Although maximum relative errors of ∼17% are found for $v$ = 5, j = 85, and T = 2000 K from model NN-STS and ∼13% for $v$ = 5, j = 25, and T = 2000 K from NN-Tot, in most cases, the relative errors are <5%. Good agreement between the QCT and NN rates can be seen in the correlation diagrams shown in Fig. S11.

As another test, total thermal rates k(T) were calculated from QCT simulations and compared with the results obtained from the NN models; see Fig. 5 and the supplementary material. The NN rates are calculated from integrating the NN cross sections over the (⁠$v$⁠, j)-state and translational energy phase space using Monte Carlo integration. For NN-Tot, the agreement with QCT is particularly striking, whereas it is still good for NN-STS except for temperatures 1000 K and 20 000 K. Here, it should be mentioned that only a range of 0.05–5.5 eV in translational energy is covered by the QCT calculations to generate the state-to-state data. Sampling a broader range of energies will most likely further improve rates determined from the NN model in the low and high temperature regions.

As a final validation of the robustness of NNs, the distribution of the final vibrational and rotational states and the rovibrational energies of N₂ after N + NO collisions at different temperatures are calculated from QCT and compared with those from NN-STS; see Fig. 6. The NN correctly captures the shape of all distributions but lacks the oscillatory features, in particular, for the rotational distribution. Thus, the NN provides a physically robust model based on validated, microscopic data from which information about nonequilibrium systems can be obtained, obviating the construction of models based on simple, empirical expressions.²¹ 

In this work, an NN-based model for state-to-state cross sections has been constructed. For this purpose, a total of ∼8 × 10⁶ state-to-state cross sections have been explicitly determined from QCT simulations for selectively chosen 1232 initial states. Local averaging over (⁠$v$′, j′) reduces the noise of the data set, and IS to sample impact parameters for the trajectories further accelerates the convergence of the cross sections. Typical training times for the NN-STS models are a few days on a 64 bit 2.40 GHz Intel E5-2620 v3 central processing unit (CPU) workstation using 8 processors, and the evaluation time of the NN for 10⁶ state-to-state cross sections is 24.4 s on the same computer using only a single processor. This makes the technique suitable for direct use in DSMC simulations where a large number of collision cross sections are required to model hypersonic air flow. The average error from the NN compared with the reference QCT data is ∼5%. This compares with errors ranging from 25% to 60% for vibrational relaxation rates and state-specific dissociation rates from a maximum entropy model for O₂ + O.²² Unfortunately, the error for the state-to-state cross sections is not reported and cannot be compared here.

In summary, the state-to-state cross sections for a reactive collision relevant to the hypersonic flight regime has been successfully modeled using a neural network based on explicitly calculated data from QCT simulations on an accurate, fully dimensional reactive PES. Such an approach is general and demonstrates that for situations in which large amounts of data constitute the relevant state space, subsampling and subsequent machine learning can provide a viable, accurate, and computationally tractable alternative to exhaustive computations.

See the supplementary material for the methodologies used to compute QCT and NN rate coefficients, Figs. S1–S11 and Table S1.

Support from the Swiss National Science Foundation through Grant No. 200021-117810, the NCCR MUST (to M.M.), and the University of Basel is also acknowledged. Part of this work was supported by the United State Department of the Air Force which is gratefully acknowledged (to M.M.).

The codes for training the NN-model and computing the cross sections from the NN-model are freely available from https://github.com/MeuwlyGroup/NNcross together with a sample training data set.