不均一系触媒上の反応解析（NEB法）¶

前節では有機分子を用いて、反応経路探索を行い活性化エネルギーを算出するNEB法の使い方を学びました。本節では、より実用的な応用例として、不均一系触媒Rh上でのNO解離の反応解析を行ってみましょう。

流れとしては以下の目次のようになります。触媒となるRh Slab部分の構造を作成し、その上に反応させたい対象となるNOを吸着させます。その後、反応前後の状態を作成し、NEB法を用いて活性化エネルギーを算出します。

※本チュートリアルの内容はmatlantis-contribで公開されている計算事例を基にしています。

matlantis-contrib/matlantis_contrib_examples/NEB_Solid_catalyst

セットアップ¶

[1]:

import numpy as np
import pandas as pd
from IPython.display import display_png
from IPython.display import Image as ImageWidget
import ipywidgets as widgets
import matplotlib as mpl
import matplotlib.pyplot as plt
import os
import glob
from pathlib import Path
from PIL import Image, ImageDraw

import ase
from ase import Atoms, units
from ase.units import Bohr,Rydberg,kJ,kB,fs,Hartree,mol,kcal
from ase.io import read, write
from ase.build import surface, molecule, add_adsorbate
from ase.constraints import FixAtoms, FixedPlane, FixBondLength, ExpCellFilter
from ase.neb import NEB
from ase.vibrations import Vibrations
from ase.thermochemistry import IdealGasThermo
from ase.visualize import view
from ase.build.rotate import minimize_rotation_and_translation
from ase.optimize import BFGS, LBFGS, FIRE
from ase.md import MDLogger
from ase.io import read, write, Trajectory
from ase.build import sort

from pfp_api_client.pfp.calculators.ase_calculator import ASECalculator
from pfp_api_client.pfp.estimator import Estimator, EstimatorCalcMode


estimator = Estimator(calc_mode=EstimatorCalcMode.CRYSTAL_PLUS_D3, model_version="v3.0.0")
calculator = ASECalculator(estimator)

1. BulkからSlab作成¶

触媒として、今回はRhを用います。 Slab構造は、Bulk構造の構造ファイルから作成ができます。

1-1 Bulk構造を読み込みから作成まで¶

今回はMaterials Projectからダウンロードしたcifファイルをinputフォルダに入れて読み込み

Input cif file is from
A. Jain, S.P. Ong, G. Hautier, W. Chen, W.D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder, K.A. Persson (*=equal contributions)
The Materials Project: A materials genome approach to accelerating materials innovation APL Materials, 2013, 1(1), 011002.
doi:10.1063/1.4812323
[bibtex]
Licensed under CC BY 4.0

[2]:

bulk = read("../input/Rh_mp-74_conventional_standard.cif")
print("Number of atoms =", len(bulk))
print("Initial lattice constant =", bulk.cell.cellpar())


bulk.calc = calculator
opt = LBFGS(ExpCellFilter(bulk))
opt.run()
print ("Optimized lattice constant =", bulk.cell.cellpar())

Number of atoms = 4
Initial lattice constant = [ 3.843898  3.843898  3.843898 90.       90.       90.      ]
       Step     Time          Energy         fmax
*Force-consistent energies used in optimization.
LBFGS:    0 07:14:41      -25.513042*       3.6306
LBFGS:    1 07:14:42      -24.835918*      13.6274
LBFGS:    2 07:14:42      -25.584866*       0.6765
LBFGS:    3 07:14:42      -25.587189*       0.1175
LBFGS:    4 07:14:43      -25.587259*       0.0013
Optimized lattice constant = [ 3.793282    3.79328225  3.79328165 90.00000028 90.00000252 90.00001485]

sort関数は、元素番号に応じて順序をソートする関数です。

[3]:

from pfcc_extras.visualize.view import view_ngl

bulk = bulk.repeat([2, 2, 2])
bulk = sort(bulk)
# Shift positions a bit to prevent surface to be cut at wrong place when `makesurface` is called
bulk.positions += [0.01, 0, 0]

view_ngl(bulk, representations=["ball+stick"])

1-2 Slab構造を作成まで¶

bulk構造から任意のミラー指数でslab構造を作成。 miller_indices=(x,y,z)で指定できます。 makesurface は中で surface 関数を使用して表面構造を作成しています。

[5]:

from ase import Atoms
from ase.build import surface


def makesurface(
    atoms: Atoms, miller_indices=(1, 1, 1), layers=4, rep=[4, 4, 1]
) -> Atoms:
    s1 = surface(atoms, miller_indices, layers)
    s1.center(vacuum=10.0, axis=2)
    s1 = s1.repeat(rep)
    s1.set_positions(s1.get_positions() - [0, 0, min(s1.get_positions()[:, 2])])
    s1.pbc = True
    return s1

[6]:

slab = makesurface(bulk, miller_indices=(1, 1, 1), layers=2, rep=[1, 1, 1])
slab = sort(slab)
# adjust `positions` before `wrap`
slab.positions += [1, 1, 0]
slab.wrap()

view_ngl(slab, representations=["ball+stick"])

1-3 作成したslabのz座標を確認¶

slabの最も高い座標確認（吸着構造を作成するときに必要） slabの層ごとの座標を確認（何層目までを固定するのか決めるのに必要）

[7]:

import plotly.express as px

# Check z_position of atoms
atoms = slab
df = pd.DataFrame({
    "x": atoms.positions[:, 0],
    "y": atoms.positions[:, 1],
    "z": atoms.positions[:, 2],
    "symbol": atoms.symbols,
})
display(df)

coord = "z"
df_sorted = df.sort_values(coord).reset_index().rename({"index": "atom_index"}, axis=1)
fig = px.scatter(df_sorted, x=df_sorted.index, y=coord, color="symbol", hover_data=["x", "y", "z", "atom_index"])
fig.show()

	x	y	z	symbol
0	1.007071	0.995918	4.458544e-08	Rh
1	10.394965	1.770218	2.190052e+00	Rh
2	6.371580	8.738923	2.190052e+00	Rh
3	2.348199	1.770218	2.190052e+00	Rh
4	14.418346	8.738923	2.190052e+00	Rh
...	...	...	...	...
59	5.030453	3.318820	6.570156e+00	Rh
60	1.007071	0.995918	6.570156e+00	Rh
61	5.030453	4.867420	4.380104e+00	Rh
62	6.371581	2.544519	4.380104e+00	Rh
63	7.712709	4.867420	4.380104e+00	Rh

64 rows × 4 columns