Supporting Information for “Discovering Surface Reaction Pathways Using Accelerated Molecular Dynamics and Network Analysis Tools”

We present an automated method that maps surface reaction pathways with no experimental data and with minimal human interventions. In this method, bias potentials promoting surface reactions are applied to enable statistical samplings of the surface reaction within the timescale of ab initio molecular dynamics (MD) simulations, and elementary reactions are extracted automatically using an extension of the method constructed for gas- or liquid-phase reactions. By converting the extracted elementary reaction data to directed graph data, MD trajectories can be efficiently mapped onto reaction pathways using a network analysis tool. To demonstrate the power of the method, it was applied to the steam reforming of methane on the Rh(111) surface and to propane reforming on the Pt(111) and Pt3Sn(111) surfaces. We discover new energetically favorable pathways for both reactions and reproduce the experimentally-observed materials-dependence of the surface reaction activity and the selectivity for the propane reforming reactions.


I. SAMPLING ELEMENTARY REACTIONS
Here, we describe a method to extract elementary reactions from the trajectories of molecular dynamics (MD) simulations.

A. Identification of chemical species
At first, we describe a method to extract the chemical species contained in a system from the atomic coordinates obtained by MD simulations. For each atomic pair, bonding or no bonding was described using a "bond matrix", where i and j are the labels specifying atoms (here, atoms composing the slab are excluded).
This matrix was updated every 5 steps of the MD simulation. In this study, bonding or no bonding was determined from the inter-atomic distance; if the distance is shorter than the cutoff radius, the pair of atoms is judged to form a bond. Different cutoff parameters were used to detect bond creations (M [i, j] : 0 → 1) and bond breakings (M [i, j] : 1 → 0). Shorter cutoff parameters were used for bond creations, and longer cutoff parameters were used for bond breakings to avoid wrong assignments of molecular vibrations (bond stretching and contraction) and non-reactive collisions (temporary reduction of inter-atomic distance) Isolated atomic groups can be identified by the block-diagonalization of the bond matrix.
The block-diagonalization can be executed by permutations of rows and columns of the bond matrix. This process is schematically shown in FIG. 1. An isolated atomic group can be found in each block in the block-diagonalized matrix, and the identified atomic group corresponds to the isolated molecule. In this study, a brute-force method was used for the block-diagonalization. The composition formula (C x H y O z ) of each molecule can be constructed from the atoms that constitute each block. However, if isomers present, the composition formula is not sufficient. In such a case, we can use SMILES format to represent molecules in text lines. The molecular structure file like XYZ format can be constructed for each molecule from the MD trajectory data and atomic IDs. The molecular structure files can be converted to SMILES format using OpenBabel software, for example.
In this manner, the type and number of chemical species present in the system at every MD step can be determined.

B. Identification of chemical reactions
Here, we describe a method to extract elementary reactions from the trajectories of molecular dynamics simulations. Changes in the bond matrix during the MD simulations mean changes in the bonds between atoms caused by chemical reactions. Therefore, chemical reactions can be extracted by detecting the time variation of the bond matrix.
Let M t be the bond matrix at t-th MD step, and assume that M t is the matrix, whose rows and columns are arranged in the atomic ID order before the block diagonalization as shown in FIG. 1. Here, we compare M t and M t−1 . If there are differences in the matrix elements between M t and M t−1 , atoms corresponding to these matrix elements are judged to be involved in chemical reactions.
It should be, however, noted that strong molecular vibrations and inter-molecular collisions at high temperatures can result in temporary changes in the bond matrix, and those changes can be wrongly assigned to chemical reactions even though we try to prevent them by using two different cutoff values. To avoid the wrong assignments, not only M t and M t−1 but also M t+1 and M t−2 were used to judge whether chemical reactions happen or not. Even if there is a difference between M t and M t−1 , if M t−2 = M t or M t−1 = M t+1 , the change from M t−1 to M t is judged to just a temporal change and is not assigned to a reaction. The use of the history of the bond matrix in this way can drastically reduce the probability of the wrong assignments.
When the change in the bond matrix is assigned to a chemical reaction, the molecules involved to the reaction can be determined from the atomic IDs that exhibit the changes in the matrix elements. The molecules involved in the reaction obtained from M t−1 are assigned to the reactants, for example, A and B, and those obtained from M t are assigned to the products, C and D. By using the identified species, the reaction equation can be output as A + B → C + D. As in this example, multiple molecules can be involved in the identified chemical reaction. Without any sorting of the chemical species, the reaction, A + B → C + D and B + A → C + D (or A + B → D + C and so on) can be judged to be different reaction although they are actually identical. To prevent this, the extracted molecules are sorted according to a predetermined procedure. In this study, the molecules were sorted according to their molecular names (composition formula).