最終更新：2025-06-29
概要
本記事では、自作モジュール(MultiOptPy)で、Cope脱離反応の遷移状態構造を算出してみる。電子状態計算ソフトウェアにPySCFを使用する。 MultiOptPyは電子状態計算ソフトウェアを用いた分子構造最適化手法の勉強を目的として作成したpythonモジュールである。
MultiOptPyのレポジトリ:https://github.com/ss0832/MultiOptPy
Cope脱離反応:
- https://www.chem-station.com/odos/2009/08/cope-cope-elimination.html
- https://en.wikipedia.org/wiki/Cope_reaction
使用した自作モジュールMultiOptPyのバージョン
1.9.12b環境
WSL2 (Ubuntu-22.04)Source codeのダウンロード
wget https://github.com/ss0832/MultiOptPy/archive/refs/tags/v1.9.12b.zip
unzip v1.9.12b.zip
cd MultiOptPy-1.9.12b手順
1. 初期構造の準備
モデル反応系として、以下の構造を用意した。今回は名前をcope_elimination.xyzとした。 初期構造は以下のものを使用した。
20
0 1
C -1.37702991 1.37788020 0.22017445
H -1.42697748 1.89877175 -0.71314074
H -2.36790550 1.16460025 0.56306980
C -0.65184375 2.25375222 1.25870862
H -1.36918338 2.68244532 1.92694965
H 0.03753287 1.65222596 1.81352936
C 0.11358826 3.37837475 0.53694828
H -0.57578835 3.97990101 -0.01787246
H 0.61745138 3.98693518 1.25852722
H 0.83092789 2.94968165 -0.13129275
N -0.63933713 0.11921395 0.03998389
O -0.57585237 -0.54285382 1.22625366
C 0.72195924 0.41222472 -0.43109664
H 1.25891929 -0.50394732 -0.56225575
H 0.67201167 0.93311627 -1.36441182
H 1.22582236 1.02078515 0.29048230
C -1.33156029 -0.71684571 -0.95134418
H -1.38150786 -0.19595416 -1.88465936
H -0.79460024 -1.63301774 -1.08250329
H -2.32243588 -0.93012566 -0.60844883
2. 遷移状態構造最適化のための初期構造の算出
※あらかじめpythonで利用できる電子状態計算ソフトウェアのPySCF(導入方法はpipを使えば可能である。)が使える環境を用意しておく。
a. 初期構造をカレントディレクトリにcope_elimination.xyzとして保存し、以下を実行する。
python optmain.py cope_elimination.xyz -pyscf -spin 0 -elec 0 -func hf -bs 3-21g -opt rsirfo_fsb -ma 300 5 12 -ns 50-opt rsirfo_fsbは準ニュートン法であるRS-I-RFO法を構造最適化に使用することを示す。初期のへシアンに関しては、特にオプションで指定しない限り、単位行列が使われる。-pyscfは使用する電子状態計算ソフトウェアをPySCFにすることを示す。-spinはスピン多重度の指定である。PySCFを使用するときは目的とするスピン多重度に1を引いた値を指定する。-elec 1は形式電荷を1とすることを示す。-ma 300 5 12は300kJ/molの活性化障壁を超えうるペア同士を近づける力を原子のラベル番号12と5のペアに構造最適化時に加えることを示す。-func hfはHartree-Fock法を用いることを示す。-bs 3-21gは基底関数にPople系である3-21Gを用いることを示す。-ns 50は50回電子状態計算を反復計算を行う。（デフォルトは1000回）脱離反応で1分子から2分子にわかれる反応で構造最適化の収束性が単分子の場合よりも悪くなるので、時間短縮のために指定している。
これを実行すると、計算レベルHF/3-21Gで指定した人工力ポテンシャルをかけた上で初期構造を構造最適化することができる。
結果はyyyy_mm_dd（今日の年月日）ディレクトリの中に存在するディレクトリを開いて確認できる。
正常終了していれば、このディレクトリ中に、cope_elimination_traj.xyzが存在するので、これをコピーして、MultiOptPy-1.9.12bディレクトリに置く。
cope_elimination_traj.xyzは構造最適化の過程をAvogadro（公式ページ:https://avogadro.cc/ ）等で可視化して確認できるようにしている。このcope_elimination_traj.xyzは次のNEB計算に使用する。
※cope_elimination_traj.xyzをアニメーションとして表示したい場合は、[https://github.com/ss0832/molecule_movie] を使うと良い。
（この人工力ポテンシャルを加えて行った構造最適化の結果はcope_elimination_optimized.xyzで確認できる。構造を可視化して、生成系になっているか確認する。反応系のままであれば、-maの設定を見直してb.をやり直す。今回の場合以下の構造が得られた。生成系の安定構造に近いものが得られていると判断できるため、次のNEB計算を行うことが可能である。）
b. cope_elimination_traj.xyzを初期パスとして、NEB法で経路の緩和を行う。
NEB法を用いることで、先ほど得られたcope_elimination_traj.xyz全体のエネルギーを下げることができる。これにより、パスのエネルギー極大値を持つ構造を遷移状態構造に近づける。（この時点ではまだ正確な遷移状態構造は求められていない。）
python nebmain.py cope_elimination_traj.xyz -ns 10 -modelhess -func hf -bs 3-21g -pyscf -spin 0 -elec 0 -spng -nd 0.5-nd 0.5はノード間の距離を0.5Åとして初期パスを作成することを示す。-ns 10は10回分NEB法による経路の緩和を行うことを示す。-modelhessは少ない計算コストで算出された近似へシアンを用いることを示す。これを指定すると経路緩和のアルゴリズムが-fcを1以上に指定した場合と同様のものになる。-pyscfは使用する電子状態計算ソフトウェアをPySCFにすることを示す。-spngは緩和中のパスのエネルギープロファイルや各ノードの勾配のRMS値をmatplotlibで可視化するオプションである。-spinはスピン多重度の指定である。PySCFを使用するときは目的とするスピン多重度に1を引いた値を指定する。
c. 初期構造の決定及び遷移状態構造の計算
MultiOptPy-1.9.12bと同じディレクトリ内に、NEBという名前を含むディレクトリが生成されている。 そのディレクトリ内のenergy_plot.csvを確認し、緩和後のパスのエネルギー極大値を示す構造を確認する。
パスの緩和後の各ノードのエネルギー一覧(単位
私が実行した環境では、13番のノード(グラフでは14番)がエネルギー極大値を示していた。 path_ITR_10_cope_eliminationディレクトリ内のcope_elimination_traj_13.xyzを可視化する。
※こちら[https://ss0832.github.io/molecule_viewer/] を使うことでも可視化は可能である。
cope_elimination_traj_13.xyz
20
0 0
C -0.425427318617 0.595487772791 -0.457251709162
H -0.099122412924 0.677867968222 -1.487258758089
H -1.494223269952 0.775330044323 -0.394718605513
C 0.334233693381 1.236595240428 0.685950829733
H -0.074491209951 -0.036458408699 1.478243458150
H 1.402252911232 1.127173538496 0.516089778125
C 0.095326589304 2.739169820372 0.991263001823
H -0.944367940651 3.052960458560 0.952147368511
H 0.455469339962 2.974797174888 1.990917679823
H 0.652863694284 3.396765229670 0.318045790504
N -0.205659360435 -0.994118856352 -0.169698424423
O -0.279663536868 -1.125913073833 1.263722630518
C 1.202634039629 -1.328320887432 -0.555941850209
H 1.342045688603 -2.394877933736 -0.507860304103
H 1.393911991420 -0.966067553681 -1.555189991286
H 1.836571085058 -0.844569423810 0.160599533413
C -1.142154887977 -2.023463568660 -0.720627508178
H -1.222976287389 -1.905479730984 -1.790839865745
H -0.727321408627 -2.993347297932 -0.483239414216
H -2.099901399483 -1.963530512631 -0.234353639676
構造が壊れていないので、これを遷移状態を求めるための初期構造とする。
cope_elimination_traj_4.xyzをMultiOptPy-1.9.12bと同じディレクトリ内にコピーする。
そして、以下を実行する。
python optmain.py cope_elimination_traj_13.xyz -opt rsirfo_bofill -func hf -bs 3-21g -pyscf -spin 0 -elec 0 -fc 5 -order 1 -tcc -freq-opt rsirfo_bofillは遷移状態構造の最適化向けのoptimizerを指定することを意味する。準ニュートン法であるRS-I-RFO法を使用する。今回は-fcで正確なへシアンを計算するようにしているので、初期へシアンは正確なへシアンを使用するようになっている。-order 1は一次の鞍点を求めることを指定する。（デフォルトだと極小値を求めるようになっている。）-fc 5は5回の反復回数当たり1回正確なへシアンを計算することを指定する。-freqは収束条件を満たした後に基準振動解析を行うことを示す。（自前で実装しているため、あくまで目安として使用することを推奨する。各振動モードをvibration_animation内のxyzファイルで可視化できる。）-tccは収束条件を厳しくすることを示す。（Gaussianのtightと同等）
実行して得られた正確な遷移状態構造を以下に示す。
(実行して得られた正確な遷移状態構造はcope_elimination_traj_4_optimized.xyzとして保存されている。)
20
OptimizedStructure
C 0.273263644284 0.803744456009 -0.707484269949
H 1.064266120134 0.652350556758 -1.420084088604
H -0.698976568121 0.888968057247 -1.161951286700
C 0.547195483821 1.505879852956 0.486236882943
H 0.547589478355 0.238050864355 1.241261880491
H 1.548506250215 1.893298574928 0.584396200498
C -0.549356268609 2.375767804839 1.105208123177
H -1.508676563180 1.865084365999 1.080409500305
H -0.329192208084 2.588433319209 2.146405668013
H -0.669547036659 3.330350762607 0.595979242329
N 0.110724991583 -0.987145991583 -0.165154519238
O 0.377149841524 -0.914850051243 1.227017952873
C 1.144818882774 -1.837067490899 -0.781847952530
H 1.084750563412 -2.836702003552 -0.374031111223
H 1.005845261305 -1.862756073969 -1.854687664335
H 2.105126800113 -1.411223068399 -0.539416284706
C -1.264833011985 -1.486954803263 -0.340193747339
H -1.521353174227 -1.491988824587 -1.391410743271
H -1.343888951411 -2.485348478014 0.067608236263
H -1.923413535246 -0.827891829396 0.201737981001
30回程度の反復計算で遷移状態構造が得られた。-freqオプションにより生成されたnormal_modes.txtやvibration_animationディレクトリ内の振動モードのアニメーションを確認した。
以下に-freqオプションで生成されたnormal_modes.txtの一部を示す。
Mode 0 1 2
Freq [cm^-1] -1254.6297 78.0603 131.9957
Reduced mass [au] 1.4895 2.7896 1.2699
Force const [Dyne/A] -1.3814 0.0100 0.0130
Char temp [K] 0.0000 112.3114 189.9123
Normal mode x y z x y z x y z
C -0.02193 -0.10139 -0.02682 -0.09446 -0.00269 -0.00795 -0.01621 -0.00880 0.02453
H -0.00983 0.05444 -0.04575 -0.17072 -0.00009 -0.09288 -0.00885 -0.00998 0.03323
H 0.00322 0.05029 -0.04854 -0.14280 0.00095 0.09603 -0.00987 0.00092 0.01335
C 0.00974 -0.06784 0.08781 0.01979 -0.00807 -0.03053 -0.03256 -0.01423 0.03364
H 0.07559 0.74098 -0.21659 -0.02572 0.00599 -0.01595 -0.05344 -0.01775 0.02467
H -0.03032 0.09847 -0.10091 0.04459 -0.05298 -0.10365 -0.01370 -0.06703 0.06494
C -0.00151 0.00585 -0.00793 0.10784 0.04950 0.04462 0.00841 0.08889 -0.03789
H 0.00562 -0.00524 0.00036 0.08212 0.09443 0.12177 -0.11363 0.33943 -0.46857
H 0.00461 0.00981 -0.01214 0.19884 0.04727 0.02595 -0.25438 -0.24152 0.08582
H 0.00686 -0.00885 -0.02854 0.11599 0.05135 0.04636 0.40236 0.25356 0.17854
N 0.00832 0.06343 -0.00397 -0.02703 -0.00903 -0.00275 0.00217 -0.01457 0.00279
O -0.00147 -0.01086 0.00157 -0.07906 0.01172 0.00604 0.02393 -0.03504 -0.00059
C 0.00017 0.01169 -0.00833 0.05606 0.05757 0.04273 0.00191 0.01411 -0.03693
H -0.04298 0.03392 0.02875 0.12112 0.06258 0.06397 0.01806 0.00029 -0.06856
H 0.01556 -0.01289 -0.00984 0.08981 0.02393 0.03936 -0.01349 0.04692 -0.03587
H 0.00158 -0.00771 0.01053 0.01628 0.13749 0.06119 0.00100 0.01603 -0.03712
C 0.00018 0.01167 -0.00819 0.01420 -0.10130 -0.05642 0.00601 -0.03354 0.02407
H -0.02078 -0.00764 -0.00314 0.05204 -0.10470 -0.06586 -0.01459 -0.03234 0.02921
H 0.05726 0.01945 0.01063 0.06739 -0.11289 -0.07391 0.03004 -0.03651 0.02146
H 0.00023 -0.00680 0.01039 -0.05017 -0.15442 -0.07004 0.00690 -0.04626 0.04115
(...snip...)その結果、虚振動が１つであることが確認できた。
次に、vibration_animation内のmode_1_1255i_wave_number.xyzをAvogadroで確認すると、想定される反応系と生成系をつなぐ方向に振動していることを確認できた。
終わりに
自作モジュールで、電子状態計算モジュールであるPySCFを用いて、計算レベルHF/3-21GでCope脱離反応の遷移状態構造を算出する手順を説明した。
参考
- https://github.com/pyscf/pyscf (PySCFのgithubのレポジトリ)
- https://github.com/ss0832/MultiOptPy (自作モジュールMultiOptPyのレポジトリ)
- https://avogadro.cc/ (Avogadro、分子構造可視化ツール)
- The Journal of Chemical Physics 2010, 132, 241102.
- The Journal of Chemical Physics 1991, 94, 751–760.
- In Classical and Quantum Dynamics in Condensed Phase Simulations; WORLD SCIENTIFIC: LERICI, Villa Marigola, 1998; pp 385–404.
- The Journal of Chemical Physics, 2020, 153, 024109.
- The Journal of Chemical Physics, 2022, 144, 214108.