反応経路解析 - 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 14:43:15 -41.288286 3.372026
BFGS: 1 14:43:15 -41.050344 9.227204
BFGS: 2 14:43:15 -41.352642 0.937929
BFGS: 3 14:43:15 -41.360438 0.846813
BFGS: 4 14:43:15 -41.367057 0.414904
BFGS: 5 14:43:16 -41.370927 0.222537
BFGS: 6 14:43:16 -41.373629 0.187670
BFGS: 7 14:43:16 -41.374702 0.132663
BFGS: 8 14:43:16 -41.376004 0.143219
BFGS: 9 14:43:16 -41.377716 0.194960
BFGS: 10 14:43:16 -41.379438 0.192537
BFGS: 11 14:43:16 -41.380632 0.108121
BFGS: 12 14:43:16 -41.381465 0.084829
BFGS: 13 14:43:16 -41.382037 0.099253
BFGS: 14 14:43:16 -41.382358 0.063755
BFGS: 15 14:43:16 -41.382590 0.051461
BFGS: 16 14:43:16 -41.382839 0.060071
BFGS: 17 14:43:16 -41.383148 0.073773
BFGS: 18 14:43:16 -41.383478 0.065854
BFGS: 19 14:43:16 -41.383675 0.038820
BFGS: 20 14:43:16 -41.383762 0.025375
BFGS: 21 14:43:16 -41.383797 0.025825
BFGS: 22 14:43:16 -41.383857 0.029769
BFGS: 23 14:43:17 -41.383902 0.027067
BFGS: 24 14:43:17 -41.383944 0.014089
BFGS: 25 14:43:17 -41.383969 0.009115
Step Time Energy fmax
BFGS: 0 14:43:17 -42.925393 0.730117
BFGS: 1 14:43:17 -42.911635 1.663660
BFGS: 2 14:43:17 -42.928937 0.251647
BFGS: 3 14:43:17 -42.928764 0.549965
BFGS: 4 14:43:17 -42.931020 0.205362
BFGS: 5 14:43:17 -42.933635 0.252208
BFGS: 6 14:43:17 -42.937590 0.417864
BFGS: 7 14:43:17 -42.942083 0.419333
BFGS: 8 14:43:17 -42.946968 0.288300
BFGS: 9 14:43:17 -42.950906 0.179478
BFGS: 10 14:43:17 -42.954262 0.270988
BFGS: 11 14:43:17 -42.956673 0.231530
BFGS: 12 14:43:17 -42.958089 0.118511
BFGS: 13 14:43:17 -42.959289 0.146021
BFGS: 14 14:43:17 -42.960235 0.184809
BFGS: 15 14:43:17 -42.960773 0.116383
BFGS: 16 14:43:18 -42.961079 0.068664
BFGS: 17 14:43:18 -42.961431 0.119403
BFGS: 18 14:43:18 -42.962067 0.201879
BFGS: 19 14:43:18 -42.963184 0.261755
BFGS: 20 14:43:18 -42.964901 0.241058
BFGS: 21 14:43:18 -42.966974 0.131766
BFGS: 22 14:43:18 -42.969068 0.111116
BFGS: 23 14:43:18 -42.970628 0.153895
BFGS: 24 14:43:18 -42.971380 0.111868
BFGS: 25 14:43:18 -42.971700 0.061562
BFGS: 26 14:43:18 -42.971887 0.035948
BFGS: 27 14:43:18 -42.971988 0.045980
BFGS: 28 14:43:18 -42.972005 0.024488
BFGS: 29 14:43:18 -42.972040 0.013546
BFGS: 30 14:43:18 -42.972020 0.021100
BFGS: 31 14:43:18 -42.972060 0.040577
BFGS: 32 14:43:18 -42.972098 0.064734
BFGS: 33 14:43:18 -42.972179 0.091686
BFGS: 34 14:43:19 -42.972329 0.095636
BFGS: 35 14:43:19 -42.972412 0.052399
BFGS: 36 14:43:19 -42.972583 0.078331
BFGS: 37 14:43:19 -42.972768 0.136407
BFGS: 38 14:43:19 -42.972924 0.122819
BFGS: 39 14:43:19 -42.972981 0.193897
BFGS: 40 14:43:19 -42.973057 0.047662
BFGS: 41 14:43:19 -42.973107 0.047561
BFGS: 42 14:43:19 -42.973201 0.088515
BFGS: 43 14:43:19 -42.973237 0.042605
BFGS: 44 14:43:19 -42.973302 0.058457
BFGS: 45 14:43:19 -42.973299 0.032923
BFGS: 46 14:43:19 -42.973306 0.007837
可視化を行い、反応前後の構造を確認します。
[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
FIRE: 0 14:43:27 -38.199809 23.490956
FIRE: 1 14:43:27 -39.254621 6.505189
/tmp/ipykernel_25517/2722927025.py:8: FutureWarning: Please import NEB from ase.mep, not ase.neb.
neb = NEB(images, k=0.1, climb=True, allow_shared_calculator=False, parallel=True)
FIRE: 2 14:43:27 -39.008617 11.163567
FIRE: 3 14:43:27 -39.405844 7.057152
FIRE: 4 14:43:27 -39.450020 6.904893
FIRE: 5 14:43:27 -39.462435 5.353684
FIRE: 6 14:43:27 -39.478987 3.718861
FIRE: 7 14:43:27 -39.489786 3.016063
FIRE: 8 14:43:27 -39.491198 2.897458
FIRE: 9 14:43:27 -39.488352 3.915679
FIRE: 10 14:43:27 -39.488606 4.039161
FIRE: 11 14:43:28 -39.495308 3.512413
FIRE: 12 14:43:28 -39.507829 2.810146
FIRE: 13 14:43:28 -39.519259 1.817654
FIRE: 14 14:43:28 -39.525174 3.265853
FIRE: 15 14:43:28 -39.539539 2.892644
FIRE: 16 14:43:28 -39.572475 1.963444
FIRE: 17 14:43:28 -39.604964 2.303778
FIRE: 18 14:43:28 -39.606562 2.028061
FIRE: 19 14:43:28 -39.609304 1.513654
FIRE: 20 14:43:28 -39.612562 0.948627
FIRE: 21 14:43:28 -39.615873 0.909884
FIRE: 22 14:43:28 -39.619310 1.198012
FIRE: 23 14:43:29 -39.623196 1.533371
FIRE: 24 14:43:29 -39.627874 1.530900
FIRE: 25 14:43:29 -39.633955 1.165640
FIRE: 26 14:43:29 -39.641255 0.737797
FIRE: 27 14:43:29 -39.649404 0.800478
FIRE: 28 14:43:29 -39.658483 1.200173
FIRE: 29 14:43:29 -39.669464 1.085851
FIRE: 30 14:43:29 -39.682563 0.602644
FIRE: 31 14:43:29 -39.696479 0.907790
FIRE: 32 14:43:29 -39.710647 1.034686
FIRE: 33 14:43:29 -39.725190 0.620371
FIRE: 34 14:43:29 -39.737168 1.120402
FIRE: 35 14:43:29 -39.747340 0.514008
FIRE: 36 14:43:30 -39.752792 1.085785
FIRE: 37 14:43:30 -39.755755 0.841935
FIRE: 38 14:43:30 -39.753593 0.537654
FIRE: 39 14:43:30 -39.748983 1.064867
FIRE: 40 14:43:30 -39.752546 0.434176
FIRE: 41 14:43:30 -39.753732 0.642694
FIRE: 42 14:43:30 -39.754038 0.579657
FIRE: 43 14:43:30 -39.754581 0.468322
FIRE: 44 14:43:30 -39.755184 0.342757
FIRE: 45 14:43:30 -39.755784 0.333293
FIRE: 46 14:43:30 -39.756348 0.378653
FIRE: 47 14:43:30 -39.757035 0.433039
FIRE: 48 14:43:31 -39.757874 0.411865
FIRE: 49 14:43:31 -39.758987 0.319442
FIRE: 50 14:43:31 -39.760280 0.299208
FIRE: 51 14:43:31 -39.761602 0.399914
FIRE: 52 14:43:31 -39.762984 0.443207
FIRE: 53 14:43:31 -39.764582 0.367706
FIRE: 54 14:43:31 -39.766270 0.265688
FIRE: 55 14:43:31 -39.767651 0.306568
FIRE: 56 14:43:31 -39.768935 0.271192
FIRE: 57 14:43:31 -39.769791 0.382343
FIRE: 58 14:43:31 -39.770105 0.282708
FIRE: 59 14:43:31 -39.769909 0.320444
FIRE: 60 14:43:32 -39.769471 0.269958
FIRE: 61 14:43:32 -39.768488 0.312202
FIRE: 62 14:43:32 -39.768029 0.365170
FIRE: 63 14:43:32 -39.768429 0.204112
FIRE: 64 14:43:32 -39.768648 0.220764
FIRE: 65 14:43:32 -39.768724 0.337588
FIRE: 66 14:43:32 -39.769099 0.279886
FIRE: 67 14:43:32 -39.769610 0.161808
FIRE: 68 14:43:32 -39.769876 0.233502
FIRE: 69 14:43:32 -39.770145 0.223181
FIRE: 70 14:43:32 -39.770639 0.153745
FIRE: 71 14:43:33 -39.770829 0.265734
FIRE: 72 14:43:33 -39.771117 0.161783
FIRE: 73 14:43:33 -39.771184 0.234144
FIRE: 74 14:43:33 -39.771296 0.141023
FIRE: 75 14:43:33 -39.770925 0.323859
FIRE: 76 14:43:33 -39.771124 0.230935
FIRE: 77 14:43:33 -39.771124 0.242028
FIRE: 78 14:43:33 -39.770795 0.364439
FIRE: 79 14:43:33 -39.771262 0.119569
FIRE: 80 14:43:33 -39.771155 0.249417
FIRE: 81 14:43:33 -39.771227 0.204118
FIRE: 82 14:43:33 -39.771324 0.144965
FIRE: 83 14:43:33 -39.771389 0.111566
FIRE: 84 14:43:34 -39.771388 0.133739
FIRE: 85 14:43:34 -39.771367 0.164100
FIRE: 86 14:43:34 -39.771395 0.168005
FIRE: 87 14:43:34 -39.771444 0.144330
FIRE: 88 14:43:34 -39.771531 0.111891
FIRE: 89 14:43:34 -39.771593 0.123567
FIRE: 90 14:43:34 -39.771598 0.152820
FIRE: 91 14:43:34 -39.771657 0.138157
FIRE: 92 14:43:34 -39.771756 0.104620
FIRE: 93 14:43:34 -39.771798 0.132303
FIRE: 94 14:43:34 -39.771839 0.114626
FIRE: 95 14:43:34 -39.771912 0.116256
FIRE: 96 14:43:35 -39.771909 0.119960
FIRE: 97 14:43:35 -39.771974 0.113715
FIRE: 98 14:43:35 -39.771948 0.106299
FIRE: 99 14:43:35 -39.771956 0.120135
FIRE: 100 14:43:35 -39.772026 0.129044
FIRE: 101 14:43:35 -39.772069 0.153474
FIRE: 102 14:43:35 -39.772048 0.258043
FIRE: 103 14:43:35 -39.772141 0.064037
FIRE: 104 14:43:35 -39.772104 0.207688
FIRE: 105 14:43:35 -39.772121 0.166595
FIRE: 106 14:43:35 -39.772154 0.094351
FIRE: 107 14:43:35 -39.772153 0.059529
FIRE: 108 14:43:36 -39.772165 0.118646
FIRE: 109 14:43:36 -39.772160 0.113806
FIRE: 110 14:43:36 -39.772162 0.104494
FIRE: 111 14:43:36 -39.772173 0.091514
FIRE: 112 14:43:36 -39.772172 0.076141
FIRE: 113 14:43:36 -39.772176 0.063003
FIRE: 114 14:43:36 -39.772170 0.055623
FIRE: 115 14:43:36 -39.772164 0.055830
FIRE: 116 14:43:36 -39.772183 0.061785
FIRE: 117 14:43:36 -39.772179 0.068822
FIRE: 118 14:43:36 -39.772177 0.072309
FIRE: 119 14:43:36 -39.772180 0.069552
FIRE: 120 14:43:37 -39.772191 0.061102
FIRE: 121 14:43:37 -39.772186 0.053883
FIRE: 122 14:43:37 -39.772210 0.058257
FIRE: 123 14:43:37 -39.772206 0.067317
FIRE: 124 14:43:37 -39.772221 0.065830
FIRE: 125 14:43:37 -39.772235 0.053900
FIRE: 126 14:43:37 -39.772238 0.054710
FIRE: 127 14:43:37 -39.772259 0.058955
FIRE: 128 14:43:37 -39.772262 0.049006
もしログを抑制したい場合は、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
/tmp/ipykernel_25517/4187733795.py:8: FutureWarning: Please import NEB from ase.mep, not ase.neb.
neb = NEB(_images, k=0.1, climb=True, allow_shared_calculator=False, parallel=True)
[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など |
参考文献¶
反応経路探索手法についてさらに知りたい方は、以下の文献も参考になるでしょう。