メインコンテンツへスキップ

IRC計算実行方法

  • 更新

IRCとは内部反応座標 (Intrinsic Reaction Coordinate)のことで、NEBなどの手法によって見つけた反応の遷移状態（Transition state）を起点として、反応物(Reactant)と生成物(Product)の方向にそれぞれエネルギーが降るように構造を最適化して最小のエネルギの反応経路を生成する手法です。

この記事ではZádor groupが開発したSellaを使ったIRC計算の実行方法について説明します。

https://github.com/zadorlab/sella

Sellaのinstallは 

pip install Sella

で実行でききます。https://pypi.org/project/Sella/

記事作成時点（2024年7月）では、

  • python 3.9
  • ase 3.23.0
  • Sella 2.3.4

を使って動作確認をしております。

 

1. setup

必要なライブラリをimportします。

import numpy as np
from IPython.display import Image
import matplotlib.pyplot as plt

import ase
from ase.constraints import FixAtoms
from ase import Atoms

from sella import IRC
from sella import Sella, Constraints

from pfcc_extras import view_ngl

from pfp_api_client import ASECalculator, Estimator
estimator = Estimator()
calculator = ASECalculator(estimator)

 

2. TS構造の読み込み

NEB計算などを実施してTS構造を作成してください。

その計算結果から最もエネルギーが高い点かつ、その構造のforceが下がっている構造がTS構造になります。

（参考）Atomistic Simulation Tutorial 5-2, 不均一系触媒上の反応解析（NEB法） , 4. NEB計算結果の確認と遷移状態構造取得

images = ase.io.read("NEB_images.xyz", index="::")
for i in images:
    i.calc = calculator
    energies = [image.get_total_energy() for image in images]

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.plot(energies)
ax.set_xlabel("replica")
ax.set_ylabel("energy")
ax.grid(True)
ax.xaxis.set_ticks(np.array(range(0,20,2)))

fig.show()
def calc_max_force(atoms):
    return ((atoms.get_forces() ** 2).sum(axis=1).max()) ** 0.5

mforces = [calc_max_force(image) for image in images]
plt.plot(range(len(mforces)), mforces)
plt.xlabel("replica")
plt.ylabel("max force [eV]")
plt.xticks(np.arange(0, len(mforces), 2))
plt.grid(True)
plt.show()

上図の場合だと、replica No = 8がTS構造になります。

 

3. TS optimization

NEBで得られたTS構造は、厳密な鞍点まで収束させる操作が入っていませんので、このTS構造から振動解析やIRC計算を実施すると、虚振動が１つにならなかったり、おかしな経路が生成したりすることがあります。ここでは、sella というライブラリを用いて、TS構造を収束させます。

TS = images[8].copy()

c = FixAtoms(indices=[atom.index for atom in TS if atom.position[2] <= 2])
TS.set_constraint(c)

TS.calc = calculator
print((((TS.get_forces()**2).sum(1))**0.5).max())

TSopt = Sella(TS)
TSopt.run(fmax=0.03)
print((((TS.get_forces()**2).sum(1))**0.5).max())

 

4. IRC計算

TS opt後の構造を用いてIRC計算を実行します。反応の前後の方向それぞれに対して実行する必要がありますので、direction='forward'とdirection='reverse'をそれぞれ実行しております。

TS_f = TS.copy()
TS_f.calc = calculator
TS_r = TS.copy()
TS_r.calc = calculator

irc_f = IRC(TS_f, trajectory='irc_f.traj', dx=0.1, eta=1e-4, gamma=0.4)
irc_f.run(fmax=0.01, steps=9999, direction='forward')
irc_r = IRC(TS_r, trajectory='irc_r.traj', dx=0.1, eta=1e-4, gamma=0.4)
irc_r.run(fmax=0.01, steps=9999, direction='reverse')

 

生成したそれぞれのtrajを繋ぎます。

IRC_f_images = ase.io.read("irc_f.traj", index="::")
IRC_f_images.reverse()
IRC_r_images = ase.io.read("irc_r.traj", index="::")
IRC_images = IRC_f_images + IRC_r_images

 

それをviewerで表示することも可能です。

v = view_ngl(IRC_images[87], representations=["ball+stick"], replace_structure=True)
v.view.control.spin([1,0,0], np.deg2rad(-90))
v

 

gifやapng動画として出力することも可能です。

from pfcc_extras.visualize.povray import traj_to_apng
traj_to_apng(IRC_images[::4], rotation="0x,0y,0z", clean=True)

anim.png

 

energyとforceの履歴を確認しましょう。

energies = [image.get_total_energy() for image in IRC_images]

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.plot(energies)
ax.set_xlabel("replica")
ax.set_ylabel("energy")
ax.grid(True)
fig.show()
def calc_max_force(atoms):
    return ((atoms.get_forces() ** 2).sum(axis=1).max()) ** 0.5

mforces = [calc_max_force(image) for image in IRC_images]
plt.plot(range(len(mforces)), mforces)
plt.xlabel("replica")
plt.ylabel("max force [eV]")
plt.xticks(np.arange(0, len(mforces), 20))
plt.grid(True)
plt.show()

 

5. Tips

TS optがうまくいかない場合やIRC計算がうまくいかない場合は、NEBで生成したTS構造が妥当なのかを見直してください。

 

以上です。

 

 

 

 

 

 

 

 

 

 

関連記事