GRRM20 with Matlantis Application Examples

Overview

We introduce an application case provided by ENEOS Corporation as part of the joint development of GRRM20 with Matlantis. To demonstrate the accelerated calculation speed of GRRM20 with Matlantis and its applicability to large-scale systems, we performed reaction pathway searches in periodic boundary systems for 1) the oxidation reaction of carbon monoxide and 2) the dehydrogenation of methylcyclohexane (MCH) on a platinum catalyst.
In 1), we were able to obtain reaction pathway search results equivalent to those of the previous study [1] using conventional GRRM calculations with DFT. In 2), we succeeded in automatically searching the dehydrogenation process of MCH, which is expected to be a hydrogen carrier in the construction of a CO₂-free hydrogen supply chain developed by ENEOS. Also, we discovered reaction pathways that are difficult to investigate using conventional NEB and String methods, such as a reaction pathway in which two hydrogen atoms are removed at once, demonstrating the usefulness of GRRM20 with Matlantis.

CO oxidation reaction

The reaction pathway search using GRRM for the following two reactions on a platinum catalyst surface was performed. Both reactions are multi-step reactions, and it is known that analysis using the NEB method and the String method, which are widely known as reaction pathway search methods, is difficult.
1) Oxidation reaction of carbon monoxide (CO) (68-atom system)
2) Dehydrogenation of methylcyclohexane (MCH) (85-atom system)

Background

The oxidation reaction of CO on a platinum catalyst is shown in Fig. 1. Since there are few possible reaction combinations involving CO and O₂, which comprise four atoms in total, and the search space is limited, this is considered to be a relatively easy calculation to perform when searching for reaction pathways. There has been previous research on this system using GRRM and SIESTA [1].

Fig.1 Carbon monoxide (CO) oxidation reaction on a platinum slab

In the previous research, the computational details is as follows.
・Calculation level: PBE/DZP
・Vacuum layer: 15 Å
・Gamma parameter for SC-AFIR: 250

Platinum slab modeling

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

Tab.1 Platinum slab model creation conditions

item	Details
Calculation target (number of atoms)	Platinum slab (Pt 64 atoms) 1) CO+O₂ (4 atoms) 2) MCH (21 atoms)
Platinum crystal structure	mp-126 (Material Project [2]) (Miller index 111)
vacuum layer	20 Å
Structural optimization calculation conditions	・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 the global reaction pathway search, we applied kinetic navigation using the rate constant matrix contraction (RCMC) method [3][4], and performed SC-AFIR (Single component Artificial Force Induced Reaction) [5] for the CO oxidation. Then, RePATH was performed, and all reaction paths obtained by SC-AFIR were relaxed by LUP (locally updated planes) [6][7]. The initial structure is shown in Fig. 3. All slab atoms (Pt) and cell vectors (TV) were fixed in the SC-AFIR and RePATH calculations.

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

Options in the GRRM input file

In order to preferentially search for reactions in which reactants are adsorbed on the slab surface, we set a weak force (Gamma=25) between the four platinum atoms (numbers 6, 10, 13, and 17, colored light blue in Fig. 4) located in the center of the slab and the reactants (CO and O₂).
Four platinum atoms on the slab surface (numbers 6, 9, 16, and 18 in Fig. 4) were designated as AFIR targets (this was changed from the previous study's settings (numbers 6, 10, 13, and 17) to make it easier to obtain bond dissociation structures).
・The AFIR threshold was set to Gamma=300.
・The options in the GRRM input file are shown in Fig 5.

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

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

Results

The number of structures discovered by SC-AFIR and RePATH calculations is shown in Tab. 2. The number of reaction pathways discovered is the sum of the number of structures at the energy maxima of the reaction pathways (path tops, PTs) and the number of transition structures (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 by SC-AFIR were classified, and the most stable structures for each composition were extracted. In Tab. 3, the structures are compared with those in previous studies, and the structures are shown together with the relative energies based on CO ₂ (g) + O* (Tab.3 (a)).
・As shown in Fig. 6, the relative energy ordering of the structures found by GRRM20 with Matlantis is consistent with previous studies.
・In this study, the OC-OO intermediate reported in the previous study was not discovered, and we were not able to completely reproduce the previous study. However, we did discover the C* + 3O* intermediate, which was not discovered in the previous study (also confirmed in the reaction pathway network in Fig. 7).
・There are several possible reasons for this, such as different calculation conditions from the previous studies.
・In order to confirm that the C*+3O* intermediate is a new reaction pathway, it is necessary to perform a reaction pathway search with a lower threshold and increase the number of search pathways.
・As shown in Fig. 8, the multi-step reaction CO* + O ₂ * → CO ₂ (g) + O* was also extracted.

Tab. 3 Comparison of the platinum slab model with previous studies [1]
Relative energy of CO oxidation reaction pathway CO ₂ (g) + O* (previous research: ΔG _DFT, GRRM20 with Matlantis results: ΔG _PFP)
and its structure ((a)-(e) and (g) are the results of GRRM20 with Matlantis, (f) is taken from a previous study [1])

Fig. 6 Comparison of relative energies [kJ/mol]
from GRRM20 with Matlantis (ΔG _PFP)
with those from the previous study (ΔG _DFT)  as presented in Tab. 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

Background

MCH (methylcyclohexane) has been attracting attention as a hydrogen carrier for efficiently transporting CO₂-free hydrogen produced from renewable energy, and ENEOS is promoting development including its proprietary technology Direct MCH ^® [8]. The dehydrogenation reaction of MCH is the reaction for extracting hydrogen from MCH, and as shown in Fig. 9, it is possible to extract three hydrogen molecules (six hydrogen atoms H*) in the process of MCH becoming toluene. Understanding the elementary reactions is important when examining catalysts for dehydrogenation reactions.

Fig.9 Dehydrogenation reaction of methylcyclohexane (MCH)

About Direct MCH?​
・It is ENEOS's proprietary technology 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 a global reaction pathway search, we applied kinetic navigation using the rate constant matrix contraction method (RCMC) [3][4] and performed SC-AFIR2 (SC-AFIRx [9] is a method to systematically search higher energy TS) for the MCH dehydrogenation reaction. We then performed RePATH and relaxed all the obtained reaction pathways using locally updated planes (LUP) [6][7]. The initial structure is shown in Fig. 11. All slab atoms (Pt) and cell vectors (TV) were fixed in the SC-AFIR2 and RePATH calculations.

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

Options in the GRRM input file

・In order to preferentially search for reactions in which the reactants were adsorbed on the slab surface, we set a weak force (Gamma=15) between the nine platinum atoms located in the center of the slab (atoms 22, 25, 29-33, 36, and 37 in Fig. 12) and some of the C and H atoms that make up the reactant molecule MCH (colored light blue in Fig. 12).
・To improve the efficiency of the dehydrogenation process search, atoms that match the above (atoms colored light blue in the slab and MCH in Fig. 12) were specified as targets for reaction pathway search.
・The options in the GRRM input file are shown in Fig.13.

Fig.12 MCH dehydrogenation reaction
Platinum slab atom numbering (left)
and MCH (right)

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 discovered by SC-AFIR2 and RePATH calculations is shown in Tab.4. The number of reaction pathways discovered is the sum of the number of structures at the energy maxima of the reaction pathways (path tops, PTs) and the number of transition structures (transition states, TSs).

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
	SC-AFIR	RePATH
EQs	414	1060
PTs	245	725
TSs	163	481

・The equilibrium structures obtained by SC-AFIR2 were classified, and the most stable structures for each composition were extracted. Since there were a total of 37 structures extracted, Tab. 5 shows the main structures from which hydrogen H* is eliminated, along with the relative energy based on C ₇ H ₇ *+7H*.
・As shown in Tab. 5, we obtained reaction pathways in which up to seven hydrogen atoms (H*) were removed from a structure in which one hydrogen (H*) was removed.
・The dehydrogenation reaction pathway (Tab.5 (g)) in which six H* are eliminated from MCH to form toluene (TOL) was automatically explored.
・We also discovered a reaction pathway in which two hydrogen atoms are lost at once, such as C ₇ H ₁₃ * + H* → C ₇ H ₁₁ *+ 3H*, which is difficult to explore using 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

1) CO oxidation reaction and 2) MCH dehydrogenation reaction on a platinum slab were analyzed using GRRM20 with Matlantis, with automated reaction pathway search and reaction pathway relaxation calculations (RePATH) using SC-AFIR and SC-AFIR2.

・Comparing the CO oxidation reaction on a platinum catalyst surface with a previous study using GRRM calculations with SIESTA, we were able to obtain an equivalent reaction pathway and demonstrated a similar level of accuracy.
・In this study, the OC-OO intermediate from the previous study was not discovered, and the results were not a perfect reproduction of the previous study. However, we did discover the C* + 3O* intermediate, which was not discovered in the previous study.

・We performed a reaction pathway search for the dehydrogenation process of methylcyclohexane (MCH) on a platinum catalyst surface, and obtained one to seven hydrogen desorption reaction pathways, including a reaction in which six hydrogen atoms are removed from MCH to form toluene. This demonstrated that GRRM20 with Matlantis is an effective tool for investigating dehydrogenation catalysts for MCH, which is expected to be used as a hydrogen carrier.
・We also discovered a reaction pathway in which two hydrogen atoms are lost at once, which is difficult to explore using the conventional 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, HB Schlegel, J. Chem. Phys., 107, 375(1997)
[8] https://www.eneos-rd.com/research/carbon-neutral/dmch.html
[9] S. Maeda, T. Taketsugu, K. Morokuma, J. Comput. Chem. 35, 166 (2014)