Exploring complex reaction networks using neural network-based molecular dynamics simulation

Ab initio molecular dynamics (AIMD) is an established method to reveal the reactive dynamics of complex systems. However, the computational cost of AIMD restricts the explorable length and time scales to a great extent. Here, we develop a fundamentally different approach using molecular dynamics simulations powered by a neural network potential to investigate complex reaction networks. This potential is trained via a workflow combining AIMD and interactive molecular dynamics in virtual reality (VRMD) to accelerate the sampling of a rare reactive process. The capability of the methodology is demonstrated by achieving a panoramic visualization of the complex reaction networks for decomposition of a novel high explosive (ICM-102), without any predefined reaction coordinates. The study leads to the discovery of new pathways that would be difficult to uncover employing established methods. These results highlight the power of neural network-based molecular dynamics simulations for exploration of complex reaction mechanisms under extreme conditions at the ab initio level, pushing the limit of theoretical and computational chemistry towards the realism and fidelity of experiments.


Introduction
High explosives (HEs) have contributed enormously to the prosperity of humankind since the invention of black powder in China and the era of nitroglycerin brought by Alfred Nobel [1,2]. In an HE, the chemical energy stored in the bonds of the explosive material is converted into kinetic energy of the gaseous products via chemical reactions. A quantitative as well as qualitative understanding of the complex reaction network is critical to the application of HE. Such a reaction network involves thousands of elementary reactions that take place over a wide range of time scales. In addition, the reactions usually occur in extreme conditions (i.e., > 10 GPa and > 3000 K [3]), which makes it very difficult or even impossible to obtain the detailed reaction mechanisms of HE materials using experimental or conventional simulation approaches.
In recent decades, ab initio molecular dynamics (AIMD) simulations [4] have been applied to gain atomic insights into HE materials [4][5][6][7][8][9][10]. In an AIMD simulation, an HE material is represented by an atomic model with interatomic forces determined by electronic structure calculations, and the Newtonian equations of motion are solved to obtain the dynamic trajectories. Thus, AIMD allows chemical reactions (i.e. bond breaking and forming events) to occur and accounts for electronic polarization effects.
The initiation reaction of HE involves coupled processes subject to mechanical, thermal, and chemical stimulation. With the AIMD method, the response of HE material to stimulation and the subsequent reaction mechanism can, in principle, be captured.
AIMD has been applied to investigate the shock-induced mechanical response and initiation of a detonation in numerous studies, providing fundamental insight into the microscopic mechanisms of such complex phenomena. However, the computational cost of AIMD is so high that the accessible system size and time scale in simulations are limited to several hundred atoms and dozens of picoseconds, respectively.
Recently, artificial neural networks (NNs) have been applied to construct potential energy surfaces (PESs) in a fully data-driven manner, where the PES is abstracted from a well-selected training dataset using suitable functional expressions [11]. Several formulations, such as DeepMD [12], GAP [13], sGDML [14], and SchNet [15], have been proposed to develop NN potentials and have achieved success in the modeling of water [16], small organic molecules [14], and metal materials [17]. The performance of NN potential depends on the completeness of the training dataset [18]. In other words, NN does well in finding solutions in the function space of the training dataset but might fail in configurations outside the dataset. Therefore, the quality of the dataset is critical to developing NN potentials, and it is recommended to include all the critical configurations during potential reaction processes [19]. MD sampling is the most straightforward way to construct the training dataset, and the evolution of energies and forces are recorded as training datasets [20]. Various bond dissociations and recombinations exist in a reactive system, e.g., the decomposition of HE materials. Such processes cannot be sampled appropriately in classical MD sampling due to the high energy barriers. Enhanced sampling techniques, such as metadynamics [21,22] , [19], could improve the sampling quality in the potential reaction coordinates. The concept of interactive molecular dynamics in virtual reality (VRMD) was first proposed by Glowacki et al. [23] and applied to sample the PES of a simple hydrogen abstraction reaction [24,25]. This method has been proven to accelerate the intelligent curation of high-quality datasets. In a VRMD simulation, VR forces are implemented on atoms to accelerate the reaction processes, similar to metadynamics.
In this paper, we develop a neural network-based molecular dynamics (NNMD) method for investigation into the reaction dynamics of a novel high explosive (ICM-102) [26]. The training dataset is derived from both AIMD and VRMD simulations. The configuration spaces of datasets from different sampling methods are compared to yield a good NN potential. A set of MD simulations is performed to explore the elementary reactions in the decomposition of ICM-102 molecules. A panoramic visualization of the complex reaction networks is abstracted, and new pathways are identified from atomic trajectories for the first time.

Workflow of VR-enhanced NN potential
The full dataset consists of configurations sampled from AIMD and VRMD methods. A detailed description of the training dataset is shown in Table S1. As shown in Fig. 1, the following training procedures are performed: Performing AIMD simulation in the NVT ensemble at multiple temperatures (300, 1000, 2000, 3000, 4000 K); (ii) Selecting additional configurations with an active learning strategy; (iii) Extracting the configurations of potential chemical reactions with the VRMD method. Fig. 1 Schematic illustration of the workflow for NN potential training.

DFT and AIMD Calculations
DFT calculations were performed using the CP2K package [27]. Core electrons were treated using Goedecker−Teter−Hutter (GTH) pseudopotentials and the Perdew Burke Ernzerhof generalized gradient approximation method [28,29]. The Grimme DFT-D3 method [30] was used to account for dispersion interactions. A double-zeta Gaussian basis set plus polarization (DZVP-MOLOPT) [31] was considered. AIMD calculations were performed for an ICM-102 system of 160 atoms using the Quickstep module in CP2K. In both cases, simulations were carried out at constant volume and temperature conditions with periodic boundary conditions. An integration time step of 0.5 fs was used.

NN descriptor and training
We trained NN potentials using the DeePMD-kit package [32,33]. The smooth version of the deep potential model is adopted with a cut-off radius of 6.0 Å [32]. To remove the discontinuity introduced by the cut-off, the 1/r term in the network construction is smoothly switched off by a cosine shape function from 1.0 Å to 6.0 Å.
The filter (embedding) network has three layers with (25, 50, 100) nodes/layer, and the fitting net is composed of three layers, with 240 nodes each. The network is trained with the ADAM optimizer, with an exponentially decaying learning rate from 1.0 × 10 -3 to 5.0× 10 -8 . During the optimization process, the pre-factors in the loss function change from 1 to 10 and 1000 to 1 for the energy and force terms, respectively. The final NN model used for the production run was trained for 1.0 × 10 6 steps. Additional configurations were obtained using an active learning sampling procedure as implemented in the DP-GEN package [34]. Four NN potentials were trained on the same dataset with different initializations of weights and biases. NVT-MD simulations were performed at multiple temperatures (300, 3000, and 4000 K) using the LAMMPS package [35]. By comparing the structures from MD simulations, the agreement on the force predictions made by these potentials is used to select new configurations. When the deviation of NN prediction for one configuration was in the range of [0.4, 0.8] eV/Å, the corresponding structure was labeled as a candidate for the training set. The upper limit is used to filter nonphysical configurations [19].

Reaction network
The reaction network was postprocessed using ReacNetGenerator by Zeng et al. [36] from the NNMD trajectories. All the species in the reaction network are clustered by their fingerprints using the scikit-learn package [37]. The fingerprints are defined by the molecular structure and calculated by RDKit [38]. The importance of each species is calculated according to the observed number that it occurs and is represented by the size of the dot.

Sampling the PES with VRMD
The decomposition of ICM-102 involves complex PESs, and VRMD is applied as an enhanced sampling method to potential reactions. To illustrate the effect of VR forces on the sampling processes, we first sample the reaction of a single ICM-102 molecule. As shown in Fig. 3a and 1b, the ICM-102 molecule is visualized in the helmet screen. The user interacts with atoms using handles, where the movement of handles is translated as VR forces applied on atoms. In Fig. 3c, we select a reaction that involves

Dataset exploration
The performance of NN potential depends on the quality of the training dataset, which describes the chemical space of PES [18,39]. The whole training dataset of bulk ICM-102 molecules is constructed using both AIMD and VRMD sampling methods, and the detailed configurations are listed in Table S1.  Figure 4b shows the average molecular masses along PC1. When PC1 < -2, the average molecular mass is 202 a.u., corresponding to the mass of the ICM-102 molecule. As PC1 increases, the molecular mass gradually decreases, indicating the start of the decomposition reaction. Therefore, we conclude that PC1 refers to the overall reaction coordinate, where the initial stable reactants gain energy to overcome energy barriers and reach a stable product state. Then, the distributions of subsets are compared. The first subset is, in fact, obtained from AIMD simulations at 300-4000 K, combined with an active learning sampling strategy [40]. configurations. Fig. S3 shows that the NN potentials trained on the AIMD+VRMD datasets have lower mean absolute errors (MAEs) than potentials trained on only the AIMD dataset. This demonstrates that the final NN potential can perfectly reproduce the energy and forces in AIMD trajectories (Fig. S4). The configurations sampled by VRMD can be a good supplement for AIMD samples to create a balanced training set for model training. The distribution of AIMD and VRMD datasets.

Reaction network of ICM-102 decomposition
To capture the reactive dynamics of ICM-102 molecules, 100 ps MD simulations were performed using the new NN potential (e.g., case 4 in Table S2) with 64 ICM-102 molecules at 3000 K. The species evolution during ICM-102 decomposition predicted by NN and AIMD are compared in Fig. S5. It is clear that NN agrees well with the AIMD results, where ICM-102 molecules are rapidly consumed, and H2O is the main product. With the application of NN potential, the system size could be extended to thousands or even millions of atoms [41,42]. Compared to the large oscillation in the predicted species evolution of AIMD, NNMD simulates an eight-times larger system, resulting in smooth species evolution. Such a large system allows studies on the macroscopic nature of a statistical ensemble, which is beyond the ability of the traditional AIMD method. In total, 5799 products are identified from the NNMD trajectories. Direct derivation of the reaction mechanism from thousands of products can be quite challenging. We propose an interactive reaction network to illustrate the detailed reaction mechanism. In Fig. 6, the full reaction network contains detailed pathways and intermediates. These intermediates are represented by fingerprints [38] and clustered into eight groups using the k-means clustering algorithm [43]. Groups 1 and 2 represent ICM-102 and molecules produced by releasing H, O and OH radicals, respectively. Groups 3 and 4 share similar elemental compositions as groups 1 and 2 but with broken C-N rings. The decomposition of groups 5, 6, and 7 produces intermediates with a low H/O ratio (see Figure S6). Group 8 represents small gas products, including CHNO, H2O, OH, NO, and N2. These species are frequently involved in the reaction network as the product of many reaction pathways. There are two pathways for reactants to release the final gas products: group 1 => 2 => 8; group 1 => 3 and 4 => 5 and 6 => 7 => 8. The first pathway corresponds to reactions in early initiation (<1 ps) and produces molecules/radicals such as H2O, NO, OH and H. The second pathway requires a higher energy barrier and mainly occurs during 5-50 ps.  In summary, ICM-102 decomposition reactions are explored using NN-based molecular dynamics simulations. This is the first study to achieve a panoramic view of the complex explosive reaction process with an ab initio level of accuracy. NNMD extends simulation of high explosives to the sub-nanosecond scale, which helps the development of high accuracy kinetics models and promotes the design and application of high explosives.

Conclusion
This work evaluates the datasets generated from different sampling methods: AIMD and VR enhanced sampling.  [14]. However, such a strategy is highly dependent on the wellselected training set of high-level ab initio calculations. We are working on integrating the VR-enhanced sampling algorithms to develop NN potentials with a high level of accuracy.
Gaps between the MD simulations and experiments still exist in explosive decomposition reactions. Although the landscape of ICM-102 molecules has been drawn via AIMD and VR-enhanced sampling, further work is required to assemble reaction networks into kinetics models to interpret macro-phenomena in diverse fields, such as combustion, catalysis, and atmospheric chemistry. There are many wellestablished methods, such as model reduction and sensitivity analysis, to construct kinetic models and describe the macroproperties. Importantly, the present work constructs a novel VR-enhanced sampling method to train NN potentials and enables MD simulations of a novel high explosive (i.e. ICM-102) with complex reaction networks at the level of ab initio calculations, providing atomic insights into the detailed reaction pathways that are otherwise difficult to uncover, such as ring-opening reactions and the formation of N2, CHNO, and CO2. Thus, the study opens up new possibilities to build reaction kinetics models based on high-fidelity and low-cost AI algorithms.

Table of Contents Graphic
We develop a workflow combining ab initio molecular dynamics and interactive molecular dynamics in virtual reality (VRMD) to construct a neural network potential of a novel high explosive. The VRMD helps sample the reactive process and significantly improves the model performance. Our potential enables the large-scale MD simulation with an ab initio level of accuracy, achieving a panoramic view of the complex explosive reaction process.