GRRM20 with Matlantis Application Examples

Abstract

We introduce an application case provided by ENEOS Corporation as part of the joint development among three companies of GRRM20 with Matlantis. To demonstrate the accelerated computational speed of GRRM20 with Matlantis and its potential for large-scale systems, we conducted reaction pathway searches in periodic boundary systems for the following reactions: 1) oxidation of carbon monoxide on a platinum catalyst, 2) dehydrogenation of methylcyclohexane (MCH). For the former, we obtained reaction pathways comparable with those found in the preceding study [1] using GRRM calculations with DFT. For the latter, we successfully automated the search for the dehydrogenation process of MCH, which is expected to be a hydrogen carrier in the CO₂-free hydrogen supply chain developed by ENEOS. We discovered reaction pathways, such as the simultaneous removal of two hydrogen atoms, which are difficult to examine using traditional Nudged Elastic Band (NEB) and String methods, thereby demonstrating the usefulness of GRRM20 with Matlantis.

CO oxidation reaction

We have conducted reaction pathway searches using GRRM for the following two reactions on a platinum (Pt) catalyst surface, both of which involve multi-step processes. These reactions are known to be difficult to analyze using widely recognized reaction pathway search methods such as NEB and String methods.

1) Oxidation reaction of carbon monoxide (CO) (68-atoms system)
2) Dehydrogenation of methylcyclohexane (MCH) (85-atoms system)

Background

The oxidation of CO on a Pt catalyst is depicted in Fig. 1. Owing to the limited number of potential reaction combinations involving CO and O₂, which comprise four atoms in total, the search space is inherently constrained. This makes the calculation for reaction pathway search within this system relatively straightforward. There is a preceding study[1] which has explored this system employing GRRM and SIESTA.

Fig.1 Oxidation reaction of carbon monoxide (CO) on a platinum slab surface

In the preceding study[1], the computational method is as follows.
Calculation level: PBE/DZP
Vacuum layer thickness: 15 Å
Gamma parameter for SC-AFIR : 250

Platinum slab modeling

The platinum slab model was created using the conditions outlined in Tab. 1, as illustrated in Fig. 2.

Tab.1 Conditions for creating a platinum slab model

Items	Details
Calculation system (Number of atoms)	Platinum slab (Pt 64 atoms) 1) CO+O₂ (4 atoms) 2) MCH (21 atoms)
Platinum crystal structure	mp-126 (Materials Project [2]) (Miller index 111)
Vacuum layer	20 Å
Calculation Conditions for Structural Optimization	・Using ExpCellFilter ・LBFGS ・fmax= 0.005
PFP version	v5.0.0
PFP calculation mode	1) CRYSTAL_U0 2) CRYSTAL_U0_PLUS_D3

Fig. 2 Platinum slab model (Pt₆₄)

Reaction pathway search

To perform an global search for reaction pathways, we applied kinetic navigation via the rate constant matrix contraction (RCMC) method [3][4].  For the oxidation of CO, we performed SC-AFIR method (Single Component Artificial Force Induced Reaction) [5]. Subsequently, all reaction pathways identified by SC-AFIR were relaxed through LUP (Locally Updated Planes) [6][7] using RePATH. The initial structure is shown in Fig. 3. All slab atoms (Pt) and translation vectors (TV) were fixed during both SC-AFIR and RePATH calculations.

Fig. 3 Initial structure in the reaction pathway search (CO oxidation reaction)

Options in GRRM input file

・To preferentially search for reactions with molecules adsorbed on the slab surface, a weak force (Gamma=25) was set between four platinum atoms (Nos. 6, 10, 13, and 17, colored in light blue in Fig. 4) and the reacting molecules (CO and O₂) located at the center of the slab.
・Four Pt atoms on the slab surface (No. 6, 9, 16, and 18 in Fig. 4) were set as AFIR targets (changed from the settings in the previous study where Nos. 6, 10, 13, and 17 were used) to facilitate finding bond dissociation structures.
・The AFIR threshold was set to Gamma=300.
・The options in the GRRM input file are depicted in Fig. 5.

Fig.4 Numbering of slab Pt atoms (left) and reacting molecules (CO and O₂) in the CO oxidation reaction

Fig. 5 Part of the GRRM input files for SC-AFIR (left) and RePATH (right) for the CO oxidation reaction on a platinum slab

Results

The number of structures identified through SC-AFIR and RePATH calculations is presented in Tab. 2. The total number of discovered pathways corresponds to the sum of the path tops (PTs) and transition states (TSs).

Tab.2 Equilibrium structures (EQs), path tops (PTs) and transition states (TSs) identified in the calculation of CO oxidation on platinum slabs using SC-AFIR and RePATH.

	CO oxidation reaction
	SC-AFIR	RePATH
EQs	103	96
PTs	134	168
TSs	163	98

The equilibrium structures obtained with SC-AFIR were classified, and the most stable structures for each composition were extracted. Tab. 3 compares these structures with  those from the preceding study, along with their relative energies with reference to CO₂(g)+O* (see Tab. 3 (a)).
・As shown in Fig. 6, the relative energy order of the structures identified in GRRM20 with Matlantis consistent with the findings from the preceding study.
・Although the OC-OO intermediate reported in the preceding study was not found in this study, preventing a complete reproduction of the preceding results, the C* + 3O* intermediate not reported in the preceding study was discovered (this intermediate was also confirmed in the reaction pathway network depicted in Fig. 7).
・There could be several reasons for discrepancies, including differences in calculation conditions between this and the preceding study.
・ To confirm the C*+3O* intermediate as a novel reaction pathway, it is necessary to undertake more extensive search of reaction pathways by lowering thresholds, thereby increasing the number of pathways examined.
・Additionally, a multi-step reaction CO* + O₂* → CO₂(g) + O* has been identified, as illustrated in Fig. 8.

Tab. 3 Relative energies and structures of the CO oxidation reaction pathway on platinum slab model, referenced to CO₂(g) + O* compared to the preceding study [1].  The relative energies resulting from GRRM20 with Matlantis are denoted on ΔG_PFP, while the preceding study’s findings are represented as ΔG_DFT. Structures (a) to (e) and (g) are derived from GRRM 20 with Matlantis, while structure (f) is cited from the preceding study [1].

Fig. 6 Comparison of relative energies [kJ/mol] from GRRM20 with Matlantis (ΔG_PFP) with those from the preceding study (ΔG_DFT), as presented in Table 3.

Fig.7 Reaction pathway network (color coded according to relative energy ΔG, ranging from blue to red representing 0 to 350 [kJ/mol])

Fig. 8 Reaction pathway for CO* + O₂* → CO₂(g) + O*.

MCH dehydrogenation reaction

Background

MCH (methylcyclohexane) is attracting attention as a hydrogen carrier for the efficient transport of carbon dioxide-free hydrogen produced from renewable energy sources. ENEOS is developing its original technologies, Direct MCH® [8], among others. The dehydrogenation reaction of MCH is a process to extract hydrogen from MCH, and as illustrated in Fig. 9, it allows for the release of three hydrogen molecules (a total of six hydrogen atoms, denoted H*) in the conversion of MCH into toluene. A thorough understanding of the fundamental reactions is essential for investigating catalysts for dehydrogenation reactions.

Fig.9 Dehydrogenation reaction of methylcyclohexane (MCH)

About Direct MCH🄬
・It is a proprietary technology owned by ENEOS for the production of MCH, which serves as a hydrogen carrier.
・ENEOS is developing technologies to create a CO₂-free hydrogen supply chain that employs MCH as a hydrogen carrier.

Fig.10 Conceptual diagram illustrating the CO₂-free hydrogen supply chain utilizing MCH (Source: ENEOS Corporation)

Platinum slab modeling

The same platinum slab model was used for the MCH reaction as for the CO oxidation reaction.

Reaction pathway search

 To perform an extensive search for global reaction pathways, we employed kinetic navigation using rate constant matrix contraction (RCMC) method [3][4]. Subsequently, we utilized SC-AFIR2 (SC-AFIRx [9] is a method designed for systematically searching transition states of higher energy) for the MCH dehydrogenation reaction. Following this, RePATH was executed, and all resulting reaction pathways were relaxed using LUP (locally updated planes) [6][7]. The initial structure is depicted in Fig. 11, with all slab atoms (Pt) and translation vectors (TV) being fixed during the SC-AFIR2 and RePATH calculations.

Fig. 11 Initial structure in the reaction pathway search (MCH dehydrogenation reaction)

Options settings in the GRRM input file

・To preferentially search for the reactions where the reactive molecule is adsorbed onto the slab surface, we set a weak force (Gamma=15) between the nine central Pt atoms of the slab (atoms numbered 22, 25, 29 to 33, 36, 37 in Fig. 12) and selected carbon (C) and hydrogen (H) atoms of the reactive molecule MCH (indicated in light blue in Fig. 12).
・In order to improve the efficiency of the search for the dehydrogenation pathways, we targeted the same atoms that are colored in light blue (associated with both the slab and MCH in Fig. 12) in our reaction pathway search.
・The settings for the options in the GRRM input file are presented in Fig. 13.

Fig.12 Numbering of platinum slab atoms (left) and MCH (right) for the MCH dehydrogenation reaction.

Fig.13 Part of GRRM input files for SC-AFIR2 (left) and RePATH (right) in the MCH dehydrogenation reaction on platinum slabs.

Results

The number of structures identified by the SC-AFIR2 and RePATH calculations is presented in Tab. 4. The number of reaction pathways discovered is given by the total number of path tops (PTs) and transition states (TSs) structures.

Tab.4 Number of equilibrium structures (EQs), structures at energy maxima of reaction pathways (path tops, PTs) and transition states (TSs) identified in SC-AFIR2 and RePATH calculations for the dehydrogenation of MCH on platinum slabs.

	MCH dehydrogenation reaction
	SC-AFIR	RePATH
EQs	414	1060
PTs	245	725
TSs	163	481

・The equilibrium structures derived from SC-AFIR2 were classified, and the most stable structures for each composition were identified. Due to the extraction of a total of 37 compositions, Tab. 5 selectively presents the principal structures involving the release of hydrogen (H*), along with their relative energies with reference to C₇H₇*+7H*.
・As depicted in Tab. 5, reaction pathways encompassing the dissociation of one H* atom up to seven H* atoms were determined.
・The dehydrogenation pathway where six H* atoms are removed from MCH resulting in toluene (TOL) as demonstrated in Tab. 5 (g) was automatically searched.
・Moreover, we uncovered pathways like C₇H₁₃* + H* → C₇H₁₁*+ 3H*, which involve the simultaneous dissociation of two hydrogen atoms, which are difficult to explore using conventional approaches such as NEB calculations.

Tab.5 Key structures observed in the dehydrogenation pathway of MCH on platinum catalyst surfaces, with the dissociated hydrogen atoms (H*) highlighted in red, and their respective relative energies [kJ/mol] (with the structure having seven dissociated H* atoms (h) as the reference)

Summary

We used GRRM20 with Matlantis to perform automated reaction pathway searches using SC-AFIR and SC-AFIR2 methods and relaxation calculations (RePATH) for 1) the CO oxidation reaction and 2) the MCH dehydrogenation reaction on a platinum slab.

・When we compared our CO oxidation reaction pathways on the platinum catalyst surface to those from the preceding study that used GRRM with SIESTA, we found that we could achieve similar reaction pathways with comparable accuracy.
・In our results, the previously reported OC-OO intermediate was not observed, indicating that our replication of the preceding study was not complete; however, we discovered an unreported C* + 3O* intermediate.
・Regarding the MCH dehydrogenation on a platinum catalyst surface, our pathway search identified several processes for removing one to seven hydrogen atoms, including a significant reaction that removes six hydrogens from MCH to form toluene. This supports the use of MCH as a potential hydrogen carrier and confirms the usefulness of GRRM20 with Matlantis as a tool for investigating dehydrogenation catalysts.
・Additionally, we identified complex reaction pathways involving the simultaneous removal of two hydrogen atoms, which are difficult to find using conventional methods such as NEB and String methods.

References

[1] K. Sugiyama, Y. Sumiya, M. Takagi, K. Saita, S. Maeda, Phys.Chem.Chem.Phys., 21, 14366(2019)
[2] Commentary: The Materials Project: A materials genome approach to accelerating materials innovation Anubhav Jain, Shyue Ping Ong, Geoffroy Hautier, Wei Chen, William Davidson Richards, Stephen Dacek, Shreyas Cholia, Dan Gunter, David Skinner, Gerbrand Ceder, and Kristin A. Persson
[3] Y. Sumiya, Y. Nagahata, T. Komatsuzaki, T. Taketsugu, S. Maeda, J. Phys. Chem. A, 119, 11641(2015)
[4] Y. Sumiya, S. Maeda, Chem. Lett., 49, 553(2020)
[5] S. Maeda, Y. Harabuchi, Comput. Mol. Sci., 11, e1538(2021)
[6] C. Choi, R. Elber, J. Chem. Phys., 94, 751(1991)
[7] P. Y. Ayala, H. B. Schlegel, J. Chem. Phys., 107, 375(1997)
[8] https://www.eneos.co.jp/english/rd/research/carbon-neutral/dmch.html
[9] S. Maeda, T. Taketsugu, K. Morokuma, J. Comput. Chem. 35, 166 (2014)