反応経路解析 - Reaction Path Analysis¶
前章まで安定構造付近の挙動を見ていましたが、いよいよ本章では、物質の反応を解析していきます。
X→Yへの反応に伴う、活性化エネルギー \(E_a\) と生成エンタルピー \(\Delta H\) の関係。 Wikipediaより。
ある反応物Xから生成物Yができる反応過程において、その反応経路中でエネルギーが最大となる点を 遷移状態 (TS: Transition State) と呼び、反応物・生成物から遷移状態までのエネルギー差を活性化エネルギー (Activation energy) \(E_a\) と呼びます。
アレニウスの式によると、反応の速度定数 \(k\) は、
で決まります。ここで、
\(A\): 頻度因子(温度に無関係な定数)
\(E_a\): 活性化エネルギー
\(k_B\): ボルツマン定数
\(T\): 絶対温度
です。つまり、化学反応の起こりやすさは活性化エネルギー \(E_a\) によって決まります。
反応物・生成物のエネルギーが変わることはありませんが、活性化エネルギーは周辺環境、つまり触媒の存在によって反応経路が変わることで変わります。 この活性化エネルギーが低くなるような触媒を見つけられれば、対象とする反応をより起こしやすくすることができます。
参考文献
「分子レベルで見た触媒の働き 反応はなぜ速く進むのか」 松本吉泰
「熱力学で理解する化学反応のしくみ」 平山令明
反応物X、生成物Yがわかっていたとしても、遷移状態を見つけることは自明ではありません。 本節では、反応前後の構造がわかっているときにその反応途中の経路を探索するNEB (Nudged Elastic Band)法を用いて、活性化エネルギーの算出を行ってみます。
NEB¶
NEB法のイメージ図: 初期配置として反応経路候補(赤線)を生成後、構造緩和を行い遷移状態を通るような反応経路(青線)を探索する
※ 以下の内容はMatlantis製品内に含まれるNEB Tutorialとほぼ同じ内容です。
今回は、クルチウス転位と呼ばれる有機化学反応の一例を題材にします。
NEB法では、ある注目したい反応についてその前後の構造からスタートしますが、今回は直接構造を書き下して構造を設定します。
[1]:
import numpy as np
from IPython.display import Image
import matplotlib.pyplot as plt
import ase
from ase import Atoms
from ase.visualize import view
from ase.optimize import BFGS
from ase.optimize import FIRE
from ase.neb import NEB
from ase.build.rotate import minimize_rotation_and_translation
from pfp_api_client.pfp.calculators.ase_calculator import ASECalculator
from pfp_api_client.pfp.estimator import Estimator
from pfp_api_client.pfp.estimator import Estimator, EstimatorCalcMode
estimator = Estimator(calc_mode=EstimatorCalcMode.CRYSTAL, model_version="v2.0.0")
calculator = ASECalculator(estimator)
[2]:
react = ase.Atoms(
symbols="C2NON2H3",
positions = [
[ 2.12, -0.48, 0.00],
[ 0.76, 0.16, -0.00],
[-0.28, -0.79, -0.00],
[ 0.57, 1.35, 0.00],
[-1.42, -0.28, 0.00],
[-2.48, 0.07, 0.00],
[ 2.67, -0.14, -0.88],
[ 2.06, -1.57, -0.00],
[ 2.67, -0.14, 0.88],
],
)
prod = ase.Atoms(
symbols="C2NON2H3",
positions=[
[ 2.30, -0.87, 0.00],
[ 0.68, 0.92, -0.00],
[ 0.99, -0.25, 0.00],
[ 0.25, 2.03, 0.00],
[-2.34, 0.11, 0.00],
[-3.13, -0.67, -0.00],
[ 2.87, -0.57, -0.89],
[ 2.19, -1.96, 0.00],
[ 2.87, -0.57, 0.89],
],
)
まずは構造最適化計算を行います。構造最適化計算のやり方は今までと同様です。
また、後の作業のため、minimize_rotation_and_translation methodを用いて、反応前後の2つの構造がなるべく近くなるように位置を移動させておきます。
[3]:
react.calc = calculator
opt = BFGS(react)
opt.run(fmax=0.01)
prod.calc = calculator
opt = BFGS(prod)
opt.run(fmax=0.01)
minimize_rotation_and_translation(react, prod)
Step Time Energy fmax
BFGS: 0 07:11:23 -41.288279 3.3720
BFGS: 1 07:11:23 -41.050347 9.2272
BFGS: 2 07:11:24 -41.352632 0.9379
BFGS: 3 07:11:24 -41.360422 0.8468
BFGS: 4 07:11:24 -41.367067 0.4149
BFGS: 5 07:11:24 -41.370915 0.2225
BFGS: 6 07:11:25 -41.373628 0.1877
BFGS: 7 07:11:25 -41.374702 0.1327
BFGS: 8 07:11:25 -41.375998 0.1432
BFGS: 9 07:11:25 -41.377709 0.1950
BFGS: 10 07:11:26 -41.379433 0.1925
BFGS: 11 07:11:26 -41.380617 0.1081
BFGS: 12 07:11:26 -41.381467 0.0848
BFGS: 13 07:11:27 -41.382043 0.0992
BFGS: 14 07:11:27 -41.382362 0.0638
BFGS: 15 07:11:27 -41.382599 0.0515
BFGS: 16 07:11:27 -41.382840 0.0601
BFGS: 17 07:11:28 -41.383157 0.0738
BFGS: 18 07:11:28 -41.383481 0.0658
BFGS: 19 07:11:28 -41.383667 0.0388
BFGS: 20 07:11:28 -41.383759 0.0254
BFGS: 21 07:11:29 -41.383802 0.0258
BFGS: 22 07:11:29 -41.383865 0.0298
BFGS: 23 07:11:29 -41.383911 0.0271
BFGS: 24 07:11:29 -41.383936 0.0141
BFGS: 25 07:11:30 -41.383955 0.0091
Step Time Energy fmax
BFGS: 0 07:11:30 -42.925378 0.7301
BFGS: 1 07:11:30 -42.911607 1.6637
BFGS: 2 07:11:30 -42.928921 0.2516
BFGS: 3 07:11:31 -42.928744 0.5499
BFGS: 4 07:11:31 -42.931010 0.2054
BFGS: 5 07:11:31 -42.933623 0.2522
BFGS: 6 07:11:32 -42.937579 0.4179
BFGS: 7 07:11:32 -42.942069 0.4193
BFGS: 8 07:11:32 -42.946957 0.2883
BFGS: 9 07:11:32 -42.950890 0.1795
BFGS: 10 07:11:33 -42.954245 0.2710
BFGS: 11 07:11:33 -42.956652 0.2315
BFGS: 12 07:11:33 -42.958064 0.1185
BFGS: 13 07:11:34 -42.959282 0.1460
BFGS: 14 07:11:34 -42.960218 0.1848
BFGS: 15 07:11:34 -42.960760 0.1164
BFGS: 16 07:11:34 -42.961065 0.0687
BFGS: 17 07:11:34 -42.961415 0.1194
BFGS: 18 07:11:35 -42.962055 0.2019
BFGS: 19 07:11:35 -42.963172 0.2617
BFGS: 20 07:11:35 -42.964888 0.2410
BFGS: 21 07:11:36 -42.966962 0.1318
BFGS: 22 07:11:36 -42.969055 0.1111
BFGS: 23 07:11:36 -42.970605 0.1539
BFGS: 24 07:11:36 -42.971380 0.1119
BFGS: 25 07:11:36 -42.971683 0.0615
BFGS: 26 07:11:37 -42.971868 0.0360
BFGS: 27 07:11:37 -42.971954 0.0460
BFGS: 28 07:11:37 -42.972000 0.0245
BFGS: 29 07:11:37 -42.972008 0.0135
BFGS: 30 07:11:38 -42.972006 0.0211
BFGS: 31 07:11:38 -42.972034 0.0406
BFGS: 32 07:11:38 -42.972068 0.0647
BFGS: 33 07:11:38 -42.972157 0.0917
BFGS: 34 07:11:39 -42.972297 0.0956
BFGS: 35 07:11:39 -42.972386 0.0525
BFGS: 36 07:11:39 -42.972556 0.0783
BFGS: 37 07:11:39 -42.972755 0.1364
BFGS: 38 07:11:40 -42.972918 0.1228
BFGS: 39 07:11:40 -42.972963 0.1925
BFGS: 40 07:11:40 -42.973029 0.0464
BFGS: 41 07:11:41 -42.973089 0.0477
BFGS: 42 07:11:41 -42.973177 0.0895
BFGS: 43 07:11:41 -42.973234 0.0430
BFGS: 44 07:11:41 -42.973260 0.0586
BFGS: 45 07:11:42 -42.973297 0.0315
BFGS: 46 07:11:42 -42.973299 0.0083
可視化を行い、反応前後の構造を確認します。
[4]:
from ase.io import write
from IPython.display import Image
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
write("output/curtius_react.png", react, rotation="0x,0y,0z")
write("output/curtius_prod.png", prod, rotation="0x,0y,0z")
fig, axes = plt.subplots(1, 2, figsize=(6, 3))
ax0, ax1 = axes
ax0.imshow(mpimg.imread("output/curtius_react.png"))
ax0.set_axis_off()
ax0.set_title("Reactant")
ax1.imshow(mpimg.imread("output/curtius_prod.png"))
ax1.set_axis_off()
ax1.set_title("Product")
fig.show()
[5]:
from pfcc_extras.visualize.view import view_ngl
view_ngl([react, prod], representations=["ball+stick"], replace_structure=True)
ここからがNEB法の出番です。NEB法は、反応の前後の構造を離散的に補完した構造を多数用意し、それら複数の構造を位相空間上で互いに結合させつつ構造最適化を行うことで反応経路を探す手法です。
中間点の補完やNEBの実行そのもののはASEの組み込みのものを使うことができます。 以下の例のように、allow_shared_calculator=Falseを指定し、各 image に対して個別のCalculatorをセットすることで高速化が可能です。
NEBは複数の構造を内部で使う手法のため、MD計算などと比較するとステップあたりの計算時間は大きくなります。 従来、DFTなどの量子計算手法を用いる場合は数日から数ヶ月かかることもありますが、Matlantisでは本例では数分とかからずに終わらせることができます。
計算は以下のようなステップで行われています。
まず
imagesとしてASE Atomsのリストを作成します。最初が反応物
reactで、最後が生成物prodになるようにします。中間
imageは後のneb.interpolate()で座標値が変更されるので、ここではコピーを作成しておきます。今回は、中間
imageを7つ、合計9つでimagesを作成しました。
それぞれの Atoms に対して
calculatorを作成し、セットします。NEBクラスを作成します。
kが各imagesをつなげるバネの強さを表すバネ定数ですclimb=Trueとすることで、Climbing Image NEB法と呼ばれる、遷移状態を見つけるためにエネルギー勾配を上る手法を用いています。
neb.interpolate()ここで
imagesの線形補間をおこない、生成物から反応物へ座標が徐々に変わるようなリストを作成しています。この時点で生成された反応経路の初期配置(候補)は、コメントアウトしてある
view_nglの可視化を行うことで、確認できます。
FIREでopt
得られた反応経路候補をFIRE法を用いて最適化することで、適切な反応経路に修正していきます。
[6]:
images = [react.copy()]
images += [react.copy() for i in range (7)]
images += [prod.copy()]
for image in images:
estimator = Estimator(calc_mode=EstimatorCalcMode.CRYSTAL, model_version="v2.0.0")
calculator = ASECalculator(estimator)
image.calc = calculator
neb = NEB(images, k=0.1, climb=True, allow_shared_calculator=False, parallel=True)
neb.interpolate()
# Check interpolated images
# view_ngl(images, representations=["ball+stick"], replace_structure=True)
opt = FIRE(neb)
status = opt.run(fmax=0.05, steps=500)
Step Time Energy fmax
*Force-consistent energies used in optimization.
FIRE: 0 07:13:17 -38.199767* 23.4912
FIRE: 1 07:13:17 -39.254614* 6.5052
FIRE: 2 07:13:18 -39.008607* 11.1636
FIRE: 3 07:13:18 -39.405843* 7.0572
FIRE: 4 07:13:18 -39.450036* 6.9048
FIRE: 5 07:13:18 -39.462446* 5.3536
FIRE: 6 07:13:19 -39.478978* 3.7188
FIRE: 7 07:13:19 -39.489785* 3.0161
FIRE: 8 07:13:19 -39.491196* 2.8974
FIRE: 9 07:13:20 -39.488352* 3.9156
FIRE: 10 07:13:20 -39.488611* 4.0391
FIRE: 11 07:13:20 -39.495302* 3.5124
FIRE: 12 07:13:21 -39.507837* 2.8101
FIRE: 13 07:13:21 -39.519245* 1.8177
FIRE: 14 07:13:21 -39.525173* 3.2658
FIRE: 15 07:13:21 -39.539541* 2.8926
FIRE: 16 07:13:22 -39.572465* 1.9634
FIRE: 17 07:13:22 -39.604970* 2.3038
FIRE: 18 07:13:22 -39.606554* 2.0281
FIRE: 19 07:13:23 -39.609300* 1.5137
FIRE: 20 07:13:23 -39.612558* 0.9486
FIRE: 21 07:13:23 -39.615882* 0.9099
FIRE: 22 07:13:24 -39.619306* 1.1980
FIRE: 23 07:13:24 -39.623187* 1.5334
FIRE: 24 07:13:24 -39.627865* 1.5309
FIRE: 25 07:13:24 -39.633946* 1.1657
FIRE: 26 07:13:25 -39.641256* 0.7378
FIRE: 27 07:13:25 -39.649407* 0.8005
FIRE: 28 07:13:25 -39.658478* 1.2002
FIRE: 29 07:13:26 -39.669462* 1.0858
FIRE: 30 07:13:26 -39.682556* 0.6026
FIRE: 31 07:13:26 -39.696465* 0.9078
FIRE: 32 07:13:27 -39.710635* 1.0347
FIRE: 33 07:13:27 -39.725172* 0.6204
FIRE: 34 07:13:27 -39.737158* 1.1204
FIRE: 35 07:13:27 -39.747343* 0.5140
FIRE: 36 07:13:28 -39.752792* 1.0859
FIRE: 37 07:13:28 -39.755747* 0.8419
FIRE: 38 07:13:28 -39.753584* 0.5378
FIRE: 39 07:13:29 -39.748985* 1.0645
FIRE: 40 07:13:29 -39.752538* 0.4342
FIRE: 41 07:13:29 -39.753741* 0.6426
FIRE: 42 07:13:29 -39.754036* 0.5796
FIRE: 43 07:13:30 -39.754571* 0.4683
FIRE: 44 07:13:30 -39.755176* 0.3427
FIRE: 45 07:13:30 -39.755772* 0.3333
FIRE: 46 07:13:31 -39.756359* 0.3786
FIRE: 47 07:13:31 -39.757009* 0.4330
FIRE: 48 07:13:31 -39.757876* 0.4118
FIRE: 49 07:13:32 -39.758968* 0.3195
FIRE: 50 07:13:32 -39.760281* 0.2992
FIRE: 51 07:13:32 -39.761587* 0.3999
FIRE: 52 07:13:33 -39.762970* 0.4432
FIRE: 53 07:13:33 -39.764575* 0.3677
FIRE: 54 07:13:34 -39.766266* 0.2657
FIRE: 55 07:13:34 -39.767656* 0.3065
FIRE: 56 07:13:35 -39.768928* 0.2712
FIRE: 57 07:13:35 -39.769798* 0.3823
FIRE: 58 07:13:35 -39.770114* 0.2827
FIRE: 59 07:13:36 -39.769901* 0.3204
FIRE: 60 07:13:36 -39.769478* 0.2699
FIRE: 61 07:13:37 -39.768488* 0.3123
FIRE: 62 07:13:37 -39.768018* 0.3649
FIRE: 63 07:13:37 -39.768416* 0.2041
FIRE: 64 07:13:38 -39.768642* 0.2207
FIRE: 65 07:13:38 -39.768710* 0.3374
FIRE: 66 07:13:39 -39.769102* 0.2798
FIRE: 67 07:13:39 -39.769614* 0.1618
FIRE: 68 07:13:39 -39.769868* 0.2333
FIRE: 69 07:13:40 -39.770150* 0.2230
FIRE: 70 07:13:40 -39.770635* 0.1537
FIRE: 71 07:13:41 -39.770827* 0.2656
FIRE: 72 07:13:41 -39.771113* 0.1617
FIRE: 73 07:13:42 -39.771178* 0.2339
FIRE: 74 07:13:42 -39.771294* 0.1410
FIRE: 75 07:13:42 -39.770918* 0.3239
FIRE: 76 07:13:43 -39.771114* 0.2307
FIRE: 77 07:13:43 -39.771130* 0.2419
FIRE: 78 07:13:44 -39.770793* 0.3642
FIRE: 79 07:13:44 -39.771259* 0.1196
FIRE: 80 07:13:44 -39.771153* 0.2493
FIRE: 81 07:13:45 -39.771224* 0.2040
FIRE: 82 07:13:45 -39.771311* 0.1449
FIRE: 83 07:13:46 -39.771377* 0.1116
FIRE: 84 07:13:46 -39.771394* 0.1337
FIRE: 85 07:13:46 -39.771365* 0.1640
FIRE: 86 07:13:47 -39.771384* 0.1680
FIRE: 87 07:13:47 -39.771443* 0.1443
FIRE: 88 07:13:48 -39.771532* 0.1119
FIRE: 89 07:13:48 -39.771601* 0.1236
FIRE: 90 07:13:49 -39.771593* 0.1528
FIRE: 91 07:13:49 -39.771657* 0.1381
FIRE: 92 07:13:49 -39.771759* 0.1046
FIRE: 93 07:13:50 -39.771788* 0.1323
FIRE: 94 07:13:50 -39.771831* 0.1146
FIRE: 95 07:13:51 -39.771917* 0.1162
FIRE: 96 07:13:51 -39.771895* 0.1199
FIRE: 97 07:13:51 -39.771969* 0.1137
FIRE: 98 07:13:52 -39.771947* 0.1062
FIRE: 99 07:13:52 -39.771940* 0.1201
FIRE: 100 07:13:53 -39.772026* 0.1290
FIRE: 101 07:13:53 -39.772065* 0.1533
FIRE: 102 07:13:53 -39.772045* 0.2577
FIRE: 103 07:13:54 -39.772143* 0.0641
FIRE: 104 07:13:54 -39.772100* 0.2081
FIRE: 105 07:13:55 -39.772121* 0.1670
FIRE: 106 07:13:55 -39.772161* 0.0945
FIRE: 107 07:13:55 -39.772160* 0.0595
FIRE: 108 07:13:56 -39.772143* 0.1189
FIRE: 109 07:13:56 -39.772150* 0.1140
FIRE: 110 07:13:57 -39.772159* 0.1047
FIRE: 111 07:13:57 -39.772169* 0.0916
FIRE: 112 07:13:57 -39.772171* 0.0762
FIRE: 113 07:13:58 -39.772176* 0.0630
FIRE: 114 07:13:58 -39.772174* 0.0556
FIRE: 115 07:13:59 -39.772165* 0.0558
FIRE: 116 07:13:59 -39.772176* 0.0618
FIRE: 117 07:13:59 -39.772177* 0.0688
FIRE: 118 07:14:00 -39.772181* 0.0723
FIRE: 119 07:14:00 -39.772179* 0.0695
FIRE: 120 07:14:01 -39.772188* 0.0611
FIRE: 121 07:14:01 -39.772200* 0.0539
FIRE: 122 07:14:01 -39.772202* 0.0583
FIRE: 123 07:14:02 -39.772204* 0.0673
FIRE: 124 07:14:02 -39.772207* 0.0658
FIRE: 125 07:14:03 -39.772225* 0.0539
FIRE: 126 07:14:03 -39.772234* 0.0547
FIRE: 127 07:14:04 -39.772249* 0.0589
FIRE: 128 07:14:04 -39.772255* 0.0490
もしログを抑制したい場合は、FIRE(neb, logfile=None)のような指定をするとデフォルトのoptのログが表示されなくなります。
結果を見てみましょう。以下はmatplotlibを用いてエネルギーの軌跡を可視化したものです。
エネルギーが一度上がった後で下がる反応が見えていると思います。 このエネルギーが最大の点(遷移状態、transition state)と左右の安定点それぞれとのエネルギー差が活性化エネルギーに相当します。
[7]:
energies = np.array([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")
fig.show()
活性化エネルギーは以下のように計算できます。
ここでは、順反応X→Yにおける活性化エネルギーを E_act_forward、逆反応Y→Xにおける活性化エネルギーを E_act_backwardとしています。
[8]:
# Transition state takes maximum energy in the reaction path
ts_index = np.argmax(energies)
E_act_forward = energies[ts_index] - energies[0]
E_act_backward = energies[ts_index] - energies[-1]
print(f"ts_index = {ts_index}")
print(f"E_act_forward = {E_act_forward:.2f} eV")
print(f"E_act_backward = {E_act_backward:.2f} eV")
ts_index = 2
E_act_forward = 1.61 eV
E_act_backward = 3.20 eV
結果を見てみましょう。構造についても同様にNotebookに可視化してもいいですが、画像をファイルとして保存することもできます。今回は試しにそちらを実行してみましょう。
今回はNEBのため、結果は構造1つではなく一組の軌跡として出力されています。ASEからアニメーションGIFとして保存してみましょう。
[9]:
fig = plt.figure(facecolor="white")
ax = fig.add_subplot()
ase.io.write(
"output/curtius_NEB.gif",
images,
format="gif",
ax=ax
)
上には最終構造の静止画が表示されていると思いますが、左のファイルビューアから“curtius_NEB.gif”をダブルクリックして開くと直接アニメーションGIFを見ることができます。
以下のようにしてファイルの中身を表示することができます。
[10]:
Image("output/curtius_NEB.gif", format="gif")
[10]:
<IPython.core.display.Image object>
最後に追加で、NEB法により反応経路の最適化前後がどう変わったのかを比較してみましょう。
neb.interpolate()関数を呼び、FIRE optを行う前の反応経路の候補構造 interpolated_images とNEB法で最適化された後の images のエネルギーを比べてみます。
[11]:
_images = [react.copy()]
_images += [react.copy() for i in range (7)]
_images += [prod.copy()]
for image in _images:
estimator = Estimator()
calculator = ASECalculator(estimator)
image.calc = calculator
neb = NEB(_images, k=0.1, climb=True, allow_shared_calculator=False, parallel=True)
neb.interpolate()
interpolated_images = _images
[12]:
initial_energies = [image.get_total_energy() for image in interpolated_images]
opt_energies = [image.get_total_energy() for image in images]
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.plot(initial_energies, label="neb interpolated", color="red")
ax.plot(energies, label="after neb opt")
ax.set_xlabel("replica")
ax.set_ylabel("energy")
ax.legend()
fig.show()
上記図からもわかるように、NEBの反応経路探索により、反応経路の候補構造(赤線)よりも活性化エネルギーエネルギーが低くなる経路=遷移状態のエネルギーが低い経路(青線) が見つかったことがわかります。
[コラム] 反応経路解析手法と分子動力学法の比較¶
本節では、反応を見るために始状態・終状態を指定したNEB法により反応経路の探索を行いました。
次章で説明するような分子動力学法では、原子の時間発展を追うことができますが、そのような形でシミュレーションを行い自然と反応が起きるまで待つというアプローチはどうでしょうか？ この方法でも活性化エネルギーの解析ができる例はあり、以下のLi ionの拡散現象はそのような例となります。
しかし、一般的には反応時間が課題となります。 反応というのはレアイベントです。 アレニウスの式を見ても分かる通り、活性化エネルギー \(E_a\) に対して指数関数的に発生頻度が変わる現象です。
そのため、活性化エネルギーがとても低い場合はMDでもその反応を見ることができますが、 すこしでも高くなってくると、MDで扱える時間スケール (〜ns程度)ではその反応が起きず、効率的に見つけることはできません。 さらに活性化エネルギーが高くなれば日常の世界でもおこらないような現象となります。
参考までに、具体的な値を通して活性化エネルギーと現象の起こりやすさの関係を見てみましょう。ほぼ室温である300 Kでは、アレニウスの式の指数関数内部の分母にあたる\(k_B T\)の値はおおまかに0.026 eVとなります。これは、エネルギーが1 eV高い現象の起きやすさが\(\exp\left(-1/0.026\right)\)倍、つまり\(2\times10^{-17}\)倍に下がるということを意味します。これはマクロスケールと原子スケールの差を考慮してもなお大きな値であり、レアイベントとみなすことができます。 同様に計算すると、温度が1000 Kおよび2000 Kの場合は\(k_B T\)はそれぞれ0.086 eV, 0.172 eVとなり、エネルギーが1 eV高い現象の起きやすさはそれぞれ約\(9\times10^{-6}\)倍、 \(0.003\)倍となります。これは300 Kのときの値と比べると大幅に大きくなっており、マクロスケールでは頻繁に起きる現象とみなすことができます。このように、一般に現象の起きやすさは活性化エネルギーと温度に対して非常に鋭敏に変化するものであることが知られています。 NEB法のような反応経路解析を使うことによって、このような指数関数スケールにまたがるレアイベントを効率的に分析することが可能となります。
[13]:
from ase.units import kB
from math import exp
for T in [300, 1000, 2000]:
kBT = kB * T
ratio = exp(-1 / (kBT))
print(f"----- T = {T} K -----")
print(f"kB T : {kBT:.3f}")
print(f"ratio: {ratio:.2e}")
----- T = 300 K -----
kB T : 0.026
ratio: 1.59e-17
----- T = 1000 K -----
kB T : 0.086
ratio: 9.12e-06
----- T = 2000 K -----
kB T : 0.172
ratio: 3.02e-03
このように、MDではみることが難しいような反応も、NEB法などを使うことで効率的に探索することができる可能性があります。
この他、反応経路探索に関するいくつかのアプローチ方法を以下に記載します。
アプローチ | 特徴 | 手法例 |
|---|---|---|
特定の反応経路を対象とする反応経路探索 | ある2つの安定構造をつなぐ経路を出発点として、近傍の最小エネルギー経路(MEP: Minimum Energy Path)を求める手法。比較的少ない計算コストでMEPおよび対応する遷移状態を見つけられるが、いわゆる局所探索のため、結果は探索前に想定した反応経路に強く依存する。 | NEB法・String法など |
反応経路を事前に想定しない反応経路探索 | 反応経路や反応後の構造を事前に仮定せずに、反応経路を1つ以上見つける方法。標準的なMDを行う方法以外に、効率的にエネルギー曲面を探索することを目的としたメタダイナミクス、局所的なエネルギー曲面の情報を使って探索する手法(ADDF)などがある。 | MD、メタダイナミクス、ADDFなど |
参考文献¶
反応経路探索手法についてさらに知りたい方は、以下の文献も参考になるでしょう。