An automatic approach to explore multi-reaction mechanism for medium-sized bimolecular reactions via collision dynamics simulations and transition state searches

We develop a broadly-applicable computational method for the automatic exploration of the bimolecular multi-reaction mechanism. The current methodology mainly involves the high-energy Born-Oppenheimer molecular dynamics (BOMD) simulation and the successive reaction pathway construction. Several computational tricks are introduced, which include the selection of the reactive regions based on the electronic-structure calculations and the employment of the virtual collision-dynamics simulations with monitoring atomic distance before BOMD. These prescreening steps largely reduce the number of trajectories in the BOMD simulations and significantly save computational cost. The hidden Markov model combined with modified atomic connectivity matrix is taken for the detection of reaction events in each BOMD trajectory. Starting from several geometries close to reaction events, the further intermediate optimization and transition-state searches are conducted. The proposed method allows us to build the complicated multi-reaction mechanism of medium-sized bimolecular systems automatically. Here we examine the feasibility and efficiency of the current method by its performance in searching the mechanisms of two prototype reactions in environmental science, which are the penicillin G anion + H2O and the penicillin G anion + OH radical reactions. The result indicates that the proposed theoretical method is a powerful protocol for the automatic searching of the bimolecular reaction mechanisms for medium-sized compounds.


Introduction
The computational studies on the reaction mechanism provide the crucial understanding of chemical reactions, and the representative progresses of this topic may be found in some review articles. [1][2][3][4] At the same time, the studies of reaction mechanism also play essential roles in various chemical-relevant research fields, such as the understanding of chemical conversions of pollutant compounds in environmental sciences [5][6][7][8] and the investigation of the biological catalysis processes in biochemistry. 9,10 Thus, the development of the novel user-friendly automatic computational approach to find the reaction mechanism is an essential topic in theoretical chemistry.
One of the center tasks in the computational studies of reaction mechanism is to find transition states (TSs). The successful TS search is highly dependent on the chemical knowledge and intuition of individual researchers, and this brings difficulty for non-experts in practices. It is also possible to employ some automatic TS searching approaches, such as Nudged Elastic Band (NEB), 11 growing-string method (GSM), 12,13 and various string methods. [14][15][16] However, all above approaches heavily rely on the proper initial information of reactant and product of each elementary reaction, i.e., local minimum-energy structures in reaction pathways, while the generation of these structures may not be an easy task for unknown reactions. Therefore, it is great interesting to develop effective theoretical approaches for the automatic exploration of complicated reaction mechanism.
The development of black-box theoretical approaches for automatic reaction mechanism discovery is a great challenging task in theoretical chemistry, and several approaches were proposed in the last decades, as summarized by available reviews. 2,3,17,18 Among them, some works tried to generate the reaction mechanism by using the knowledge of the structure-reactivity information based on databases or chemical heuristic rules. [19][20][21][22][23][24][25] Differently, several groups tried to generate the possible involved species firstly, and then used them as the starting points to perform TS search in the automatic reaction mechanism construction. The possible products may be generated by the graph representation of molecular connectivity, and these species are connected by various string approaches in the automatic reaction pathway generation. 2,[26][27][28] Alternatively, the reactive intermediates may be generated by heuristic-guided approaches and then the constrained optimization may give the initial guesses of the reaction pathways. [29][30][31] Starting from reactants, it is also possible to employ the reaction-pathway guided TS search in mechanism construction by adding the bias potentials in the trajectory propagation [32][33][34][35][36][37] or distorting molecules along some critical coordinates (such as the normal-mode coordinates with strongest anharmonicity). 4,[38][39][40] In theoretical approaches of the automatic discovery of reaction mechanism, one important idea is to use the Born-Oppenheimer molecular dynamics (BOMD). Our proposed approach also heavily relies on such idea. Given reactants, the on-the-fly BOMD simulation should in principle give all possible chemical reactions, once enough energy is added into systems. [41][42][43][44][45][46][47][48] This implies that the reaction mechanism may be directly extracted from the BOMD simulations. According to such idea, Martínez-Núñez and co-workers developed a method called "the transition state search using chemical dynamics simulation (TSSCDS)". 41,42,45,46 The high-energy BOMD is run at the semi-empirical level, which provides the basic information of chemical reactions.
When the bond cleavage and formation in BOMD were detected by monitoring the atomic connectivity, the geometries near the reaction events were chosen to perform the TS optimization. The TSSCDS method was effective in the construction of unimolecular reaction network. This approach was extended to study catalysis reactions, in which initial reactant complexes were generated by putting two randomly-orientated reactive species together. 44 Recently, Martínez-Núñez, Shalashilin, Glowacki and their co-workers tried to develop a theoretical approach for the automatic discovery of reaction network by using the boxed molecular dynamics [49][50][51] in energy space, 47 which is a powerful approach to sample rare events in molecular dynamics. Later, they When several reactants are involved, the whole reaction network in principle becomes very complicated. To address this challenge, Martínez and co-workers developed an ab initio nanoreactor simulation tool. 53,54 The long-time BOMD including many compounds was run under high temperature and pressure, which allows the generation of various chemical reactions. To manage the huge computational cost, the electronic structure calculations were speeded up using graphics processor units. The hidden Markov model 55 (HMM) analysis of the atomic connectivity matrix was employed to extract reaction events in the BOMD trajectories. The geometries in adjacent to the reaction steps were taken to conduct the intermediate and TS optimization. They also developed the trajectory smoothing method to provide better initial guesses for the reaction pathway construction. 54 The nanoreactor method was further extended by Lei and their coworkers, who used DFTB 56,57 (density functional tight-binding) in BOMD. 58 Alternatively, the chemical dynamics involving many compounds may be efficiently simulated by using ReaxFF 59-61 (a widely-used reactive force field) and even neural-network potentials, and such simulations may produce the very complicated reaction network, as shown by Zhu, Zhang and co-workers. 62,63 Inspired by the above works of the reaction-mechanism discovery relying on BOMD followed by the TS optimization or reaction pathway construction, we made our initial effort to develop a theoretical method for the automatic exploration of reaction mechanism. In current stage, we focus on bimolecular reactions of mediumsized compounds. Roughly speaking, our approach includes high initial collisionenergy BOMD simulations and successive reaction pathway constructions. Several useful technical tricks were applied in our implementation, which largely reduce computational cost and significantly improve optimization convergence. Initially, the reactive regions are selected for each reactant according to the reactivity properties derived from electronic-structure calculations. As we assume that chemical reactions should only happen between chosen reactive regions, the collision dynamics between each pair of reactive centers are simulated by BOMD. To induce chemical reactions, high initial collision energies between two compounds are added in BOMD. In this implementation, we introduced the prescreening step before BOMD, by using the virtual collision dynamics with atomic distance monitoring. This largely reduces the number of non-reactive trajectories in BOMD simulation and saves computational cost substantially in the treatment of reactions involving medium-sized compounds. After BOMD, we locate reaction events in trajectories using the HMM analysis of modified atomic connectivity matrix. Starting from several geometries close to reaction events, we try to build the elementary reaction pathway. In the TS searching step, different TS optimization manners are considered. We either directly take a few geometries nearby reaction events to conduct the TS optimization, or use the optimized intermediates (local minimum-energy geometries along reaction pathways) as the starting points in the TS searching. In the latter case, the climbing image NEB (CI-NEB) 64 approach is employed for the automatic TS location. Then, the vibrational analyses of TS geometries and intrinsic reaction coordinate (IRC) 65,66 calculations help us to identify proper TS structures and construct corresponding reaction pathways. Finally, the multireaction mechanism is constructed.
We employ the current theoretical approach to study two complicated reactions, i.e., a penicillin G anion + a H2O molecule and a penicillin G anion + a OH radical, which play important roles in environmental chemistry. The presence of antibiotics in wastewater is one of the crucial environmental problems and many researchers have paid considerable attention to their chemical conversions and removals. [67][68][69][70][71][72][73][74] Penicillin G is one of the antibiotics widely used for human beings and animals in the treatment of various infections caused by bacteria. However, because of heavy use, it has been frequently detected in natural aqueous environment due to wastewater pollution, and thus its chemical reaction also received research interests. 75,76 Previous study mentioned that most penicillin G compounds in water exist as anion forms once pH value is higher than its pKa (pKa = 2.75). 76 This indicates that the study of the chemical reactions of penicillin G anions in aqueous environment is significant. Thus, we study two typical reactions, namely a penicillin G anion + a H2O molecule and a penicillin G anion + a OH radical. The importance of the first reaction is indubitable. The second reaction is also essential, as OH radical plays a critical role in photocatalysis decomposition of penicillin G in aqueous environment. As the additional tests of the proposed theoretical approach, we also investigated two simple reaction systems, NH3 + CH3CHO reactions and NH3 + HCHO reactions, which were studied by the work of Zimmerman. 26,77 The current theoretical approach allows us to explore the multi-reaction mechanism of bimolecular systems. Our long-time goal is to provide a highly-automatic black-box approach to predict the unknown reaction products for given reaction species.
As our first effort toward this goal, we show that the complicated bimolecular reaction mechanism may be identified by the current theoretical approach. In the current stage, even if we still cannot get the whole reaction network in a fully automatic manner, the proposed approach is still rather efficient. At the same time, this should be useful for researchers in chemical-relevant fields, improving the impact of theoretical chemistry on other research fields.

Theoretical Methods
A theoretical approach is developed for the automatic exploration of the bimolecular reaction mechanism, which combines various theoretical approaches, including on-the-fly reactive molecular dynamics, HMM, TS optimization and IRC.
Our purpose is to study bimolecular reactions and herein two reactant species are labelled as MOLE_A and MOLE_B in below illustration. Based on natural bond orbital (NBO) charges 78 and condensed Fukui function (CFF) 79 values, the reactive regions of reactant compounds are located. After that, the collision dynamics between two species is simulated using the on-the-fly BOMD. In this step, the dynamics is set to enforce the collision between the selected reactive regions of MOLE_A and MOLE_B. To find as many reactions as possible, different initial collision energies are taken in the BOMD simulation. For a balance of computational cost and accuracy, the BOMD simulation is run at a low-cost electronic-structure theory. To further reduce the computational cost of the dynamics simulations, we apply the preliminary virtual dynamics purely based on interatomic distances before the BOMD simulation. This allows us to get rid of the significant amount of non-reactive trajectories, and thus largely reduces the number of BOMD trajectories. After the BOMD simulation, the reaction events are identified by the HMM analysis of the atomic connectivity matrix. Starting from several geometries nearby the identified reaction time, intermediates are optimized first, which in fact correspond to the local minimum-energy structures along the reaction pathways. If two intermediates are connected by a barrier, they become the reactant and product of an elementary reaction. In this sense, there is no strict distinguishment between intermediates and reactants/products. Next, the TS search is performed using the standard TS optimization or the CI-NEB calculation. After the location of possible TS, the vibrational frequency analysis is conducted for validation. The IRC-driven pathways are built to connect reactants, TS and products.
All theoretical frameworks are discussed below and the computational details are given in the next section. The additional implementation details are also given in Supporting Information (SI).

Reactive region selection
When the medium-sized compounds are taken into account, the brute-force search of possible bimolecular reactions through collisions by BOMD may be extremely time consuming, because many trajectories are non-reactive. Thus, some tricks must be taken to avoid the non-reactive trajectories in the BOMD step, and at least we should make sure that the BOMD simulations are mainly run for the trajectories experiencing Certainly, it is also possible to apply the chemical intuition to choose reactive regions.

Reactive dynamics simulations
After the selection of reactive regions for MOLE_A and MOLE_B, we run the dynamics that gives the collision between them.

Initial sampling
In this step, the proper sampling of initial inter-molecular velocity is necessary to make sure that effective collisions appear in the dynamics simulation. Different initial collision energies are sampled to find relevant reactions as many as possible. More technical details on sampling inter-molecular motions and adding collision energies are given in SI.
For initial molecular orientation, we first fix the coordinates of MOLE_A. The collision dynamics should be aligned along with the direction connecting two individual centers of different reactive regions. All possible orientations of MOLE_B are sampled and the angular momentum is not considered. We wish to emphasize that molecular rotation generally takes place at a much longer time scale, for instance at picosecond time domain, while the current collision dynamics happens in the much shorter time scale. Thus, the sampling of angular momentum is not necessary. At the same time, we know that the orientation of molecule is extremely important for the reaction dynamics, thus we perform the sampling of different molecular orientations. In this way, the orientation-dependent reactions can be discovered. More details on how to perform the sampling of different molecular orientations are given in SI. Because the collision dynamics is used only for the generation of possible reactions, the current sampling approach should be acceptable.
The normal modes are obtained for both MOLE_A and MOLE_B individually.
The Wigner sampling 82 of the lowest vibrational level is performed to create initial coordinates and velocities.

Virtual dynamics
In principle, it is possible to run the BOMD from all sampled initial conditions. However, this gives many non-reactive trajectories that require a huge amount of computational cost. Thus, additional efforts are taken to avoid them in the BOMD stimulation and here we introduce the virtual dynamics for prescreening purpose.
The virtual dynamics simulation is realized as follows. Starting from an initial sampling condition with only collision energy included, we run the virtual dynamics simulations without potential energy part, and check the inter-atomic distance at every timestep. If an atom in the reactive region of MOLE_A meets one atom in the reactive region of MOLE_B and no other atom is involved before their collision, we believe such trajectory is effective and should be re-run with the further BOMD simulation.
Otherwise, the trajectory is assigned as non-reactive trajectory that is not taken into account in the next BOMD simulation. The idea of the virtual dynamics shares some similarity with the discrete molecular dynamics (a kind of high-energy limit of BOMD) 83 widely used in biological system simulation. The computational setups are given in Section 3 and the additional implementation details are given in SI.

BOMD
After initial sampling and virtual dynamics, the selected initial conditions are taken to run the on-the-fly BOMD simulations. Because different reactive regions, initial collision energies and initial orientations are considered here, the BOMD simulation step still requires a large number of trajectories. To save the computational cost, the BOMD should be run at a less computationally demanding electronic structure method. This setup is acceptable, because the BOMD step is only for the initial assessment of reactions. More additional discussions on the BOMD step are given in the next section and the SI.

Identification of the reaction events with modified connectivity matrix and HMM
To find elementary reaction steps in BOMD trajectories, the bond cleavage and formation are detected based on the modified inter-atomic connectivity matrix and HMM.
The element of the atomic connectivity matrix A is defined based on the interatomic distances as the following: (1) where Rij is the distance between two atoms, rij is the standard covalent bond length of the two atoms, which is taken from reference. 84 Here the upper limit defines the threshold distance to divide the bound and unbound statuses, and the value 1.4 is taken from previous works. 53,54,63 In the current BOMD simulation, the high initial collision energy is given, which indicates that the inter-atomic distance may experience a large amplitude motion even without reaction events. Thus, the employment of the connectivity matrix may provide many reaction events. If not filtered, a large number of structures should be taken into account in the next intermediate and TS optimization, largely reducing the efficiency of next steps. To remove the 'fake' reaction events in this step, the approach based on a two-state HMM is taken to analyze the atomic connectivity.
HMM is widely used in many research fields 53 where T = P(Xt|Xt-1) and 1-α (1-β) defines the probability of transition between different states. At each time step, the probability of output state observed is given by (3) where O = P(Vt|Xt), λ (μ) represents the probability to assign State 1 (State 2) as its own, while 1-λ represents the probability of the wrong assignment. Figure 1 (b) shows the relationship between 'hidden' states, observation states, transition probability and output probability.
In the reaction detection step, the connectivity matrix is calculated to identify whether a bond is broken or not. In principle, the hidden states can be viewed as the identification on whether the bond status changes or not for any pair of atoms. However, many chemical bond statues may vary at the same time and thus it is necessary to monitor the overall result. In implementation, when no bond breaks or forms, we assume that the hidden state does not change. Otherwise, the hidden state changes. The observable states now are whether the connectivity matrix changes or not. More additional discussions on the HMM model are given in SI.
The initial probabilities of two hidden states can be set as 0.5. The transition matrix is given by (4) and the output probability matrix is (5) Based on the HMM analysis of the atomic connectivity matrix, the reaction events in trajectories are identified. The values in the transition matrix T and output matrix O are not the same as the ones used in previous works that ran reactive dynamics at high temperature. 53,54,63 The current work only puts the high kinetic energy in the direction of collision between two reactive regions initially and the inner energy of each molecule is still low. In this way, the molecular geometry is not far from the ground-state minima and no strange structures appear before collision. However, the excessive collision energy induces the large-amplitude bond stretching motions during collisions. In this situation, two originally unbound atoms may quickly become very close to each other, and then return to their unbound status again immediately. Such features appear quite often in the collision dynamics, and it is necessary to get rid of them in the identification of reaction events. For this purpose, a modified connectivity matrix is proposed as (6) The elements of the B matrix, Bij, displays the ultrafast oscillation patterns, when such high-energy non-reactive collision happens. If we can detect such high-oscillation feature in the HMM analyses, such 'fake' events can be filtered out. For this purpose, we use the B matrix instead of the original A matrix to conduct the HMM analyses. Also for the same purpose, we adjust the values of two matrices (T and O) in HMM. This explains why the parameters used in the T and O matrices are different to previous works. 53,54,63 At the end, the current setup makes sure that only one reaction event is detected for most trajectories. gives the transition progress of Chain X and the output progress of visible sequence V.
T is the transition matrix and O is the output matrix.

Transition state search
When the reaction time tr is determined, a time interval [ts, tf] centered at tr is defined. Then a few of geometries equally distributed in this time interval are chosen as starting guesses to perform the intermediate optimization. In this situation, two timeadjacent selected geometries may give the same or different optimization results. In principle, an elementary step should take place in the latter case. After optimization, the latter case can be identified by checking the atomic connectivity matrix and this sub-time interval (labelled as [ti, ti+1]) containing the elementary reaction step is known.
Next, several geometries in the trajectory propagation within this sub-time interval are taken to perform the TS optimization. Alternatively, the intermediates optimized from the geometries at ti and ti+1 are chosen to perform the CI-NEB task.
After the TS optimization, the atomic connectivity matrix is used here to divide all optimized structures into several groups. To avoid repeated calculations, the frequency analysis is conducted for the representative structures of each obtained TS groups. The real TS structures are identified, when only one imaginary frequency is observed and the vibrational motion of this mode is localized in the reactive region.

IRC and reaction mechanism construction
IRC calculations are performed for all real TS geometries. The end points of the IRC path are optimized to compare with the reactant/product for an elementary reaction.
This step makes sure that the obtained reaction pathways are dynamically accessible.
At the end, the connection of all reaction species gives the multi-reaction mechanism.

Computational Details
We mainly considered the bimolecular reactions of two systems, namely a penicillin G anion + a H2O molecule and a penicillin G anion + a OH radical reactions.
As a complimentary study, we also explored two simple models, NH3 + CH3CHO reactions and NH3 + HCHO reactions, which was studied in previous works. 26,77 Initially, all electronic structure calculations were performed using the Gaussian In the virtual dynamics, we originally sampled 20000 initial conditions with different initial orientations for each pair of reactive regions. The additional discussion on the orientation sampling in the virtual dynamics is given in SI. The initial distance between them was set as 7 Å. For each pair of reactive regions, about 5000 initial conditions may give effective collisions.
After the virtual dynamics simulation, 700 initial conditions randomly chosen from the above effective samples were taken in the BOMD simulations. In BOMD run, the nuclear motion was integrated with the velocity-Verlet algorithm. The BOMD simulation was run up to 100 fs with the time step of 0.5 fs. The BOMD dynamics was run at the HF level with the 3-21G 102-107 basis set using our home-made code interfaced with Gaussian 16 package.
To avoid the missing of low-energy reactions, different collision energies were chosen in the BOMD simulation. The threshold between low-energy and high-energy reactions was set to be a rather high value 40 kcal/mol for reaction barrier. The reaction with such high barrier is very difficult to take place at the room temperature. In this implementation, we firstly ran the preliminary BOMD with one collision energy to check the reaction ratio. Next, we adjusted the initial collision energies to re-run the BOMD. Here, we first ran some test trajectories (~20) with some initial energies. If no reaction happened, we just enlarged the collision energy by a factor of 1.2 to run the test collision dynamics again. The collision energy range in these test calculations was 0.37-0.66 Hartree. For each selected collision energy, we firstly ran 500 trajectories to obtain the main results and another 200 trajectories were calculated to check convergence. When no new low-energy reaction was found in these additional 200 BOMD trajectories, we believed the calculation results are converged at this collision energy. Next, we needed to examine the low-energy reactions under different initial collision energies. Once no new low-energy reaction was found, the calculation was thought to be converged and we did not increase the collision energies any more.
The modified atomic connectivity matrix was calculated along the trajectory propagation. The HMM was employed to select the reaction events. At each reaction time, we defined a time interval ~30 fs, and chose six geometries equally distributed within such interval for the intermediate optimization, see S1-S6 in Figure 2. If two timeadjacent structures (for instance Si and Si+1) give different optimization results, we may assume that this sub-time interval contains an elementary reaction step. In the current work, we fixed the time interval to 30 fs. Although the optimization of this value is not considered in the current work, this value is acceptable because we got enough number of the TS structures (see below). Certainly, the dependence of the results on this time interval may be an interesting topic in the future research.
Next, three geometries were chosen within the time interval from Si to Si+1 to perform the direct TS optimization at the B3LYP/6-31G(d,p) level using Gaussian 16.
If the TS optimization was not converged, we tried to employ the CI-NEB approach to redo the TS search again, in which two intermediate geometries optimized from Si and Si+1 were taken as the starting and ending points. The CI-NEB work was conducted at the B3LYP/6-31G(d,p) level using Orca. 108 Here eight images were used along the guess reaction pathways and all other details in the CI-NEB calculations followed the standard setups in Orca.
After the collection of all optimized TS geometries, the frequency analyses and IRC calculations were performed to build the reaction mechanism, which were all conducted at the same level of theory using Gaussian 16. To make sure that the IRC results correctly reflect the reactions given by the BOMD, we also examined the geometry similarity between the optimized endpoints of IRC and two intermediate geometries optimized from Si and Si+1, by comparing their connectivity matrices. At the end, all species were connected to build the reaction network.

Results
Here we examine the feasibility and applicability of our proposed theoretical approach by studying two typical penicillin G anion-involved reactions that are environmentally significant. The first one is penicillin G anion + H2O reaction. And the second one is the penicillin G anion + OH radical reaction. For illustration, we label penicillin G anion as MOLE_A and the second specie (H2O or OH radical) as MOLE_B.
As a complimentary study, we also tested the performance of the current approach by using another two simple models, NH3 + CH3CHO reactions and NH3 + HCHO reactions, which were studied by the work of Zimmerman. 26,77 Since these two systems are rather small and only used for additional test, we put their results in the SI and mainly focus on two medium-sized compound reactions. All energies shown in this work are given in the electronic energies, instead of Gibbs energies. This gives the rather qualitative understanding of the reaction mechanism.

Reactions between penicillin G anion and H2O molecule
The optimized geometry of penicillin G anion is given in SI Figure S1. The NBO charges and CFF values are shown in SI Table S1. As shown in Figure 3 In principle, it is also possible to include more atoms into the reactive regions.
However, the current calculations show that all obtained elementary reactions between the reactive regions of MOLE_A and H2O already display rather high energy barriers.
At the same time, the reactivity of other non-selected atoms is much weaker according to their NBO charges and CFF values. This indicates that the reaction probability between the non-selected atoms and H2O should be extremely small. Thus we believe the above selection should already cover the primary reactions. Following the automatic reaction-mechanism discovery procedure described above, totally nine elementary reactions are obtained, as shown in Figure 4.    Figure 4). Alternatively, it is also possible that the O atom of H2O connects to the C atom of the C-N bond firstly, and then the further reaction takes place, resulting in the two-step reactions (Pathway 4 to Pathway 5, and Pathway 7 to Pathway 8 in Figure 4).
Among all hydrolysis reactions, the one-step hydrolysis reaction of the lactam region (Pathway 3 in Figure 4) shows the lowest reaction energy barrier of 35.8 kcal/mol, and thus we take this reaction as a typical example for illustration. All critical geometries (R3, TS3 and P3) in this elementary reaction step are given in Figure 5. All penicillin G anion + H2O reactions display rather high energy barriers, indicating that these reactions are not easy to take place at the room temperature. Given the high stability of the penicillin G anion in aqueous environments, we may explain why such compounds can be detected in wastewater and become a group of typical pollutants. Thus, some additional efforts must be taken for their management, for instance by adding some chemicals or using photocatalysis reactions.

Reactions between penicillin G anion and OH radical
In the second example, we study the reactions between penicillin G anion and OH radical. Because of the high reactivity, the OH radical should easily react with more regions of penicillin G anion. Therefore, we only use CFF values here to choose the reactive regions of the MOLE_A. As shown in Figure 6, five reactive regions are selected for MOLE_A. Reactive region 1 includes the six-member aromatic benzene ring. Reactive region 2 contains the peptide bond. Reactive region 3 contains the amide bond of the lactam moiety of the four-membered ring. Reactive region 4 contains the COOgroup and Reactive region 5 refers to the C-S moiety in the five-member ring.
The whole OH radical is selected as reactive region for MOLE_B. Totally nineteen different reaction pathways are found for penicillin G anion + OH radical reactions, as shown in Figure 7. These reaction pathways are roughly divided into four types. In Type-Ⅰ reactions, the OH radical simply attaches to one atom of MOLE_A. In Type-II reactions, after the OH radical connects to one atom of MOLE_A, the additional bond cleavage or formation happens, or the second-step reaction occurs.
In Type-Ⅲ reactions, the OH radical takes one H atom away from penicillin G anion and H2O is formed. Type-Ⅳ reactions lead to the bond cleavage of OH radical. Several critical geometries are shown in Figures 8.  All four types of reactions are found for the benzene moiety (Reactive region 1). Figures 6-8, Pathways 1, 2, and 3 are all Type-Ⅰ reactions with low-energy barriers, in which the OH radical simply connects to one C atom in the six-membered ring. Pathway 1 may be the barrierless reaction with E ~0.4 kcal/mol. As Type-Ⅱ reactions, both Pathways 4 and 5 are composed of two steps. The OH radical connects to one C atom in the ring, and then the H atom originally connecting with this C atom flies away. The only Type-Ⅲ reaction detected in this region follows Pathway 6, in which the OH radical captures the H atom of six-member ring and the H2O molecule is formed. The Type-Ⅳ reaction here follows Pathway 7, which includes two elementary steps. After the OH radical connecting with the C atom in the ring, the H atom of this OH radical simply goes away, while the very high barrier (E ~70.2 kcal/mol) seems to prevent this pathway.

As shown in
Only one two-step reaction exists in the peptide bond region (Reactive region 2), which is assigned as the Type-II reaction (Pathway 8). The first-step reaction denotes that the OH radical connects with the C atom in the peptide bond. In the second step, the C-N bond becomes longer (1.866 Å at TS8' in Figure 8) and the C-N bond breaks (P8' in Figure 8). The first step is the rate-deterministic step with the reaction barrier Overall, several low-energy reaction pathways are found in the penicillin G anion + OH radical reactions, due to the high reactivity of the OH radical. Thus, it is very easy to induce the chemical transformation of the penicillin G anion under the existence of OH radicals.

Discussion
Here we discuss several aspects of the current theoretical approach for the automatic exploration of reaction mechanism.
One great advantage of this theoretical method is that it can naturally be implemented in the parallel way, because different trajectories and TS finding tasks can be run independently. The price is that more convergence tests are necessary for the determination of some simulation setups, such as selecting reactive regions, choosing initial collision energies and so on.
The reactive-region selection and the virtual-dynamics step are very essential, because they allow us to run much smaller number of trajectories in the BOMD step.
These tricks largely reduce the number of non-reactive trajectories in the BOMD step and save a huge amount of simulation time. When the reaction involves the large-sized compounds, such approaches become extremely useful to reduce computational cost.
Certainly, great caution should be paid to avoid the loss of important reactions. Thus, the preliminary test calculations on how to set up these prescreening calculations are helpful.
In the selection of the reactive regions, the NBO charges and CFF values are employed and more details are given in the SI. We also noticed that the previous work by Reiher and co-workers also chose to use the similar molecular descriptors to define the reactive sites. 29 Certainly, the current selection standard may not always be optimal.
For instance, we may directly use some loose thresholds to include more atoms or try different criteria. These setups can always be examined by checking whether the observed low-energy reactions do not change with such parameters. At the same time, the involvement of chemical intuition may be highly helpful.
In a few situations, some reactions involve some atoms out of the reactive region, for instance, Pathway 6, 8 and 9 in the penicillin G anion + H2O reactions and Pathway 13 in the penicillin G anion + OH reactions. Such results may be attributed by a few reasons. For instance, the virtual dynamics and the real BOMD were run with different setups. In the BOMD dynamics the potential energy was considered and the vibrational level was sampled initially. In addition, for a two-step reaction, the second step may be fully determined by the product of the first step. Due to these reasons, a few reactions found by the current approach may involve some atoms out of the pre-defined reactive regions. However, such situations happen very occasionally.
In this work, the HF/3-21G level is used in BOMD part to generate the initial guesses for the next TS searching tasks. In fact, many BOMD-based approaches 46,53 used the low-cost electronic-structure methods, such as semiempirical method or HF, to speed up the trajectory propagation. This method in principle works here, because this step is only used to generate the initial guess structure for the successive reaction pathway construction. As a complimentary test of the employment of the HF/3-21G level in BOMD, we also studied the bimolecular reactions of small systems shown in the works by Zimmerman, 26,77 see SI for details. Our results on the low-energy reactions are more-or-less consistent with his previous findings, although some differences exist in the high-energy reactions. This indicates that our selection of the HF/3-21G level in the BOMD is reasonable. Certainly, the limitation of the employment of the HF/3-21G level is also clear, including that it may give the improper description on the reaction barriers, provide the inaccurate description on the bond dissociation limit and miss some important reactions.
Next, we employ the HMM analysis of the modified atomic connectivity matrix in the determination of the reaction events, avoiding the selection of a large number of non-reactive 'fake' reactions caused by high collision energies. This provides more reasonable initial guess for the intermediate and TS searching, and significantly improve the efficiency of the reaction pathway construction.
One BOMD trajectory can only follow a single pathway and does not give us two different reaction pathways. In principle, it can experience the reactions involving successive steps, while the connection between the multi-step reaction and BOMD trajectory is not trivial. Taking a two-step reaction as an example, it is certainly true that a single BOMD trajectory may experience both two steps. However, when we used the HMM as a filter, it is possible that only one reactive event is identified. Since we took a time interval centered at the located reactive event and chose several initial configurations for the successive intermediate and TS search, these geometries may already cover two reaction steps. As the consequence, we may find two successive reaction steps. In other situations, maybe only the first step or the second step was obtained from the reaction pathway construction. However, at the end, we need to connect all involved species to build a reaction network. In this step, we compared the similarity of all obtained species based on the inter-atomic connectivity matrix. When two species show the similar connectivity matrix, we assign them as the same intermediate. Thus, the same intermediate may be composed of different isomers for such large system. In the current work, we only put them together and did not pay attention on the relation between these isomers. Such limitation will be the research topic of our further work.
In the current work, the selection of collision energies is fully determined by the ratio of reactive trajectories in all trajectories. That's why we recommended to take the iterative procedure to determine the suitable collision energies. When no more lowenergy reaction (energy barrier < 40 kcal/mol) is found with the increasing the collision energy, we believe that the calculations are converged. Normally different reaction products may be found under different collision energies. For example, with the increasing of the collision energy, we found that more trajectories gave high-energy reactions. For different reactive regions, different optimal collision energies may be chosen according to the tested results. The current iterative procedure normally gives the rather large collision energy ( > 0.37 Hartree). We also noticed that the very high kinetic energies were employed in other works of the BOMD-based reaction mechanism discovery. 46,54 For instance, the ab initio nanoreactor used ~ 2000K to run BOMD and the instantaneous temperature imparted by the piston reached to ~ 10000K. 54 And in the TSSCDS 46 method, the more than 20 kcal/mol (per normal mode) energy was initialized to run the BOMD. They obtained reactions with barrier height ~10-100 kcal/mol. In fact, it may be very difficult to give the direct connection between reactive barrier height and such high kinetic energies energy. The reason may be that the vibrational energy redistribution is generally extremely fast in the classical dynamics simulations. Even if we add the very high initial collision energies, such energy is quickly dissipated into the whole system due to the large molecular size. Thus, it is not easy to characterize the energy flow during the collision. The characterization of such fast energy flux is beyond the scope of the current work. We wish to emphasize that the reasonable reactive network is obtained with such collision energies. Thus, for practice purpose we believe such selection of the collision energy is reasonable.
We try to address the convergence of our calculations by examining whether the low-energy reactions can be found as many as possible in the current approach. In practice, with the increase of the collision energy, we found more trajectories gave highenergy reactions. It is also important to known whether there are enough trajectories used in the BOMD step. In practice, the iterative procedure was employed to know how many trajectories are required. For instance, we ran 500 trajectories first and then 200 trajectories additionally to check whether new reactions were found. This makes sure that the results are converged. Taking the reactions between a penicillin G anion + a OH radical as example, we try to analyze all reactions in Reactive region 1, see Table 1 In principle the searching of possible reactions in the high-dimensional potential surface of a complicated system is a kind of so-called "NP-hard" problem that is basically computationally unsolvable. We can never make sure that 100% of lowenergy reaction pathways can be found, no matter which theoretical approach is employed. All approaches, including the current one, only try to find relevant reactions as many as possible. For instance, in the test example of NH3 + CH3CHO reaction provided by the previous work by Zimmerman, 77 we found the same two low-energy reactions, as well as several high-energy reactions, see SI for details. For the second example NH3 + HCHO reaction, we also found the same lowest-energy reaction (32 kcal/mol). For the second lowest-energy reactions, although we also clarified it as a H2 elimination reaction, the reactive barrier and products are different: 63 kcal/mol in this work and 80 kcal/mol in the work by Zimmerman. 26 In addition, the different products are found here. The products are H2, CO and NH3 in the work by Zimmerman, while they are H2 and NH2CO in the current work. Thus, the low-energy reactions (< 40 kcal/mol) are more-or-less addressed by our approaches. When much larger systems are involved, such as a penicillin G anion + a H2O molecule and a penicillin G anion + a OH radical, the current approach finds many low-energy reactions in a rather automatic way, therefore we believe this protocol is useful.
In the current work, we only use the B3LYP/6-31G(d,p) level to optimize the TS geometries and build the reaction pathways. In principle, this electronic level is enough to get qualitative answers for our chosen reactions. If quantitative results are necessary, more accurate electronic-structure methods, such as DFT with more accurate functionals, the second-order Møller-Plesset and even couple-cluster approaches, become suitable choices in the TS searching. This may improve the accuracy in the description of each reaction mechanism, by giving more precise energies for all species.
It is also possible to describe solvent effects by using the continuum solvation model. 105 All these additional improvements can be introduced, because any electronic structure method can in principle be applied in our current approach.

ASSOCIATED CONTENT Supporting Information
This information is available free of charge via the Internet at http://pubs.acs.org

Code Availability
The code in this work is available from the website:

Supporting Information for
An automatic approach to explore multireaction mechanism for medium-sized bimolecular reactions via collision dynamics simulations and transition state searches  Figure S1. Optimized structure of penicillin G anion with atomic labels.

Selection of reactive regions
In the selection of reactive regions, it should be possible to allow that the positive charged atom only interacts with the negative charged atom. However, this may result in too many chosen reactive pairs.
In practice, the selection was performed according to the following procedure.
In the calculations of the penicillin G anion + water reactions, the whole water was taken as a reactant region because of its small size. The reactive region of the penicillin G anion was defined by the below procedure: (1) The NBO charge and CFF value for each atom were obtained by electronic-structure calculations.
(2) The NBO charges were first considered. We first defined a fragment composed of three adjacent atoms. All H atoms were not included in these fragments because they generally display very small NBO charges. The "effective NBO value" of this fragment was given by summing the absolute NBO charges of these three atoms.
Please notice that such effective NBO value is not the overall NBO charge of this In the calculations of the penicillin G anion + OH radical reactions, the whole OH radical was taken as the whole reactant region. Here the OH radical in principle shows the very high reactivity and it can react with many parts of penicillin G anion. Therefore, we selected five reactive regions for penicillin G anion only according to the CFF values.
At the end, in the selected reactive region, the NBO charges and CFF values of the included atoms are given in Table S2. In the treatment of reactions involve only small compounds (NH3, HCHO, CH3CHO), we did not perform such pre-screening step by using NBO charges and CFF values.
Certainly, the above selection standard may not always be optimal. We may directly use some loose threshold values to include more atoms. We even can try to use different criteria to examine the convergence to check whether the observed low-energy reactions do not change with such parameters. At the same time, the involvement of chemical intuition may be highly helpful.

Initial intermolecular velocity in the collision dynamics
For two species, we tried to connect the two centers of reactive regions by a line.
Then we added the inter-molecular velocities along the direction of this line and the initial velocities were determined by the chosen collision energies. This makes sure that the direct collision takes place between two reactive regions. In the initial sampling, the total linear momentum was set to zero. Figure S2 shows the motion of a water molecule towards the selected reactive center of a penicillin G anion.

Orientation samplings
The sampling of different orientations of two reactive species is explained by considering a typical system: a NH3 molecule and a HCHO molecule. Here the NH3 molecule was located at the center of a sphere with radius 7 Å. Next a HCHO molecule was generated at a randomly chosen point on the surface of this sphere. The orientation of this HCHO molecule was also randomly sampled.
Similar to the approach explained in Figure S2, we performed the collision dynamics with proper sampling of the initial intermolecular velocity. Figure S3 shows the sampling of all possible orientations.
Different initial intermolecular velocities were taken in account in real BOMD collision dynamics. We also considered the sampling of the vibrational motion by using Wigner sampling of the lowest vibrational state.
In the virtual dynamics, the velocity of the larger compound was set to be zero and the velocity of the smaller compound was taken according to the low boundary of the collision energy. No sampling based on the vibrational motion was considered.

Distances checking in virtual dynamics
In the virtual dynamics, we monitored all pairs of inter-atomic distances between two species, as shown in Figure S4. Suppose that Atom_a, Atom_b and Atom_c belong to MOLE_B (for instance H or O), the reaction region of MOLE_A and non-reaction region of MOLE_A, respectively. When the distance between Atom_a and Atom_b is lower than the 1.4 times of their covalent bond length, we assume that the possible collision may take place between two reaction regions. Such collision is assumed to be effective, only when the distance between Atom_a and Atom_c is always larger than 1.4 times of their corresponding covalent bond length before the collision between Atom_a and Atom_b. Figure S4. All inter-atomic distances between two reaction centers are monitored in the virtual dynamics.

Additional discussions on HMM
In the HMM model, we define a hidden state (not visible) that evolves over time being and an observable state that is visible at any given time. We wish to derive the "hidden state" status at a given time step by minoring "observable state". Here a hidden state does not correspond to only a single observable state. Instead, under such hidden state, the appearance of the visible observable is determined by a probability.
In the current work, the breaking and formation of a chemical bond were characterized by the connectivity matrix. Next, we defined the binary hidden states based on the changing of connectivity matrix. When any element of the connectivity matrix changes (bond breaking or formation), the observable is set to be 1. If all elements remain unchanged (no bond breaking or formation) between two adjacent steps, the observable is set to be 0. Such binary observable states give us the direct indicator on the changing of inter-atomic connectivity. However, the changing of the connectivity matrix may reflect the bond breaking or formation, while it may only indicate the large vibrational motion. At the same time, in principle, the breaking of an old chemical bond or the formation of a new bond may prefer to give the observable value of 1 but not always. Instead, it gives either 0 or 1 according to some possibility.
The HMM model simply takes this possibility into account.
In the current BOMD, the energy is very high and thus the very strong stretching motion may take place. This gives too many oscillations of the binary observable states between 0 and 1. If we take all of them into the next optimization step, the computational cost becomes extremely high. At the same time, the selection of all state jumps is also not necessary because many of them just correspond to the strong bond stretching motion. We used HMM to filter out many oscillations and only chose a few representative ones. In one word, the original trajectory gives a large number of oscillations of binary state, while the employment of the HMM model largely reduces the number of the selection events in the BOMD. Similar approach was also used by previous works to reduce the number of the chosen reaction events. 1,2,3

Additional examples of bimolecular reactions
We tried to test this theoretical method by examining bimolecular reactions of small systems given by the works by Zimmerman. 4,5 Here we chose two test systems: (1) a NH3 molecule + a CH3CHO molecule; (2) a NH3 molecule + a HCHO molecule.
Since all these compounds are very small, the procedure would reckon the whole molecules as reactive regions and it is not necessary to perform the prescreening of the     Figure S7. All NH3 + HCHO reactions found by the current approach. E is the energy barrier of the elementary step. Figure S8. Two lowest-energy NH3 + HCHO reactions and the relevant transition state structures.
As shown in Figure S7 and Figure S8, several reaction pathways are found by using the current approach.
Among them, the lowest reaction pathway ( Figure S8 (a)) is the same as the one shown in the work of Zimmerman 5 .
The second lowest-energy reaction ( Figure S8 (b)) is a little different to their result.
Although both works clarified that this is a H2 elimination reaction, the details are different. The reactive barrier is 63 kcal/mol in our work and 80 kcal/mol in the work by Zimmerman. The detailed study on the origin of such difference is beyond the scope of our current investigation, because the barriers of both reactions are larger than 40 kcal/mol. The products are H2, CO and NH3 in the work by Zimmerman, while in our work they are H2 and NH2CO.
Except them, we also found several high-energy reactions.
To summarize our observations on these two examples, we more-or-less got the consistent results as the findings by the work of Zimmerman for the low-energy (< 40 kcal/mol) reactions. For the higher-energy reactions (for example, the second lowest NH3 + HCHO reaction with the barrier height larger than 40 kcal/mol), some deviations exist between the work of Zimmerman and the current one. For the very high-energy reactions, two approaches may give different answers. However, we only focus on the low-energy reactions, and thus we believe the current approach is trustworthy.
The finding of all reactions is in fact a NP-hard problem. In this sense, all approaches just tried to find the reactions as many as possible and no available method can make sure that all reactions are found. Figure S9. All species involved in the penicillin G anion + H2O reaction network.

All species in penicillin G anion + H2O reactions
Here we label each structure according to the following rule: the capital letters, R, TS and P that denote reactants, TS, and products respectively. The numbers after these capital letters correspond to the reaction indices given in Figure 4. The dashed line and the blue arrows show the key atomic distances with the unit of Å. Figure S10. All species involved in the penicillin G anion + OH reaction network.

All species in penicillin G anion + OH reactions
Here we label each structure according to the following rule: the capital letters, R, TS and P that denote reactants, TS, and products respectively. The numbers after these capital letters correspond to the reaction indices given in Figure 7.