最終更新：2025-06-05
概要
本記事では、自作モジュール(MultiOptPy)で、SN2反応の遷移状態構造を算出してみる。電子状態計算ソフトウェアにPySCFを使用する。 MultiOptPyは電子状態計算ソフトウェアを用いた分子構造最適化手法の勉強を目的として作成したpythonモジュールである。
MultiOptPyのレポジトリ:https://github.com/ss0832/MultiOptPy
使用した自作モジュールMultiOptPyのバージョン
1.9.11環境
WSL2 (Ubuntu-22.04)Source codeのダウンロード
wget https://github.com/ss0832/MultiOptPy/archive/refs/tags/v1.9.11.zip
unzip v1.9.11.zip
cd MultiOptPy-1.9.11手順
1. 初期構造の準備
モデル反応系として、以下の構造を用意した。今回は名前をsn2_2.xyzとした。
6
C -0.338898486699 -0.076862210703 0.014675421482
H -0.032501363200 -1.103144177303 0.014254107896
H -0.030440655315 0.438971927475 0.899746706206
H -0.029411803977 0.438402705961 -0.870064342390
Br -2.393713509620 -0.088847370762 0.017081849967
Cl 2.824965818811 0.391479125331 -0.075693743161
2. 遷移状態構造最適化のための初期構造の算出
※あらかじめpythonで利用できる電子状態計算ソフトウェアのPySCF(導入方法はpipを使えば可能である。)が使える環境を用意しておく。
a. 初期構造をカレントディレクトリにsn2_2.xyzとして保存し、以下を実行する。
python optmain.py sn2_2.xyz -opt rsirfo_fsb -func hf -bs 3-21g -pyscf -elec -1 -spin 0 -ma 100 1 6-opt rsirfo_fsbは準ニュートン法であるRS-I-RFO法を構造最適化に使用することを示す。初期のへシアンに関しては、特にオプションで指定しない限り、単位行列が使われる。-pyscfは使用する電子状態計算ソフトウェアをPySCFにすることを示す。-spinはスピン多重度の指定である。PySCFを使用するときは目的とするスピン多重度に1を引いた値を指定する。-elec -1は形式電荷を-1に指定することを示す。-ma 100 1 6は100kJ/molの活性化障壁を超えうるペア同士を近づける力を原子のラベル番号1と6のペアに、構造最適化時にそれぞれ加えることを示す。-func hfはHartree-Fock法を用いることを示す。-bs 3-21Gは基底関数にPople系である3-21Gを用いることを示す。
これを実行すると、計算レベルHF/3-21Gで指定した人工力ポテンシャルをかけた上で初期構造を構造最適化することができる。
結果はyyyy_mm_dd（今日の年月日）ディレクトリの中に存在するディレクトリを開いて確認できる。
正常終了していれば、このディレクトリ中に、sn2_2_traj.xyzが存在するので、これをコピーして、MultiOptPy-1.9.11ディレクトリに置く。
sn2_2_traj.xyzは構造最適化の過程をAvogadro（公式ページ:https://avogadro.cc/ ）等で可視化して確認できるようにしている。このsn2_2_traj.xyzは次のNEB計算に使用する。
※sn2_2_traj.xyzをアニメーションとして表示したい場合は、[https://github.com/ss0832/molecule_movie] を使うと良い。
（この人工力ポテンシャルを加えて行った構造最適化の結果はsn2_2_optimized.xyzで確認できる。構造を可視化して、生成系になっているか確認する。反応系のままであれば、-maの設定を見直してb.をやり直す。今回の場合以下の構造が得られた。安定構造に近いものが得られていると判断できるため、次のNEB計算を行うことが可能である。）
b. sn2_2_traj.xyzを初期パスとして、NEB法で経路の緩和を行う。
NEB法を用いることで、先ほど得られたsn2_2_traj.xyz全体のエネルギーを下げることができる。これにより、パスのエネルギー極大値を持つ構造を遷移状態構造に近づける。（この時点ではまだ正確な遷移状態構造は求められていない。）
python nebmain.py sn2_2_traj.xyz -ns 15 -func hf -bs 3-21g -pyscf -elec -1 -spin 0 -nd 0.2 -spng -modelhess-nd 0.2はノード間の距離を0.2Åとして初期パスを作成することを示す。-ns 15は15回分NEB法による経路の緩和を行うことを示す。原子数が少ないので、時間をかけてじっくり経路を緩和してもよいと考えて回数を増やした。-modelhessは少ない計算コストで算出された近似へシアンを用いることを示す。これを指定すると経路緩和のアルゴリズムが-fcを1以上に指定した場合と同様のものになる。-pyscfは使用する電子状態計算ソフトウェアをPySCFにすることを示す。-spngは緩和中のパスのエネルギープロファイルや各ノードの勾配のRMS値をmatplotlibで可視化するオプションである。-spinはスピン多重度の指定である。PySCFを使用するときは目的とするスピン多重度に1を引いた値を指定する。-elec -1は形式電荷を-1に指定することを示す。
c. 初期構造の決定及び遷移状態構造の計算
MultiOptPy-1.9.11と同じディレクトリ内に、NEBという名前を含むディレクトリが生成されている。 そのディレクトリ内のenergy_plot.csvを確認し、緩和後のパスのエネルギー極大値を示す構造を確認する。
パスの緩和後の各ノードのエネルギー一覧(単位
私が実行した環境では、5番のノード(グラフでは6番)がエネルギー極大値を示していた。 path_ITR_15_sn2_2ディレクトリ内のsn2_2_traj_5.xyzを可視化する。
※こちら[https://ss0832.github.io/molecule_viewer/] を使うことでも可視化は可能である。
sn2_2_traj_5.xyz
6
-1 0
C -0.082498123650 -0.078452999944 -0.000082129706
H 0.065181277197 -1.142512258015 0.000367685514
H 0.023333197617 0.443658310940 0.922024092403
H -0.000386218042 0.443609220720 -0.924458519148
Br -2.449665838055 -0.036681896521 0.030663125724
Cl 2.444035704932 0.370379622820 -0.028514254787
構造が壊れていないので、これを遷移状態を求めるための初期構造とする。
sn2_2_traj_5.xyzをMultiOptPy-1.9.11と同じディレクトリ内にコピーする。
そして、以下を実行する。
python optmain.py sn2_2_traj_5.xyz -opt rsirfo_bofill -func hf -bs 3-21g -pyscf -elec -1 -spin 0 -order 1 -fc 5 -freq -tcc-opt rsirfo_bofillは遷移状態構造の最適化向けのoptimizerを指定することを意味する。準ニュートン法であるRS-I-RFO法を使用する。今回は-fcで正確なへシアンを計算するようにしているので、初期へシアンは正確なへシアンを使用するようになっている。-order 1は一次の鞍点を求めることを指定する。（デフォルトだと極小値を求めるようになっている。）-fc 5は5回の反復回数当たり1回正確なへシアンを計算することを指定する。-freqは収束条件を満たした後に基準振動解析を行うことを示す。（自前で実装しているため、あくまで目安として使用することを推奨する。各振動モードをvibration_animation内のxyzファイルで可視化できる。）-tccは収束条件を厳しくすることを示す。（Gaussianのtightと同等）
実行して得られた正確な遷移状態構造を以下に示す。
(実行して得られた正確な遷移状態構造はsn2_2_traj_5_optimized.xyzとして保存されている。)
6
OptimizedStructure
C 0.077406676190 0.005698078622 -0.000944261168
H 0.101480321444 -1.054598865785 -0.000526542082
H -0.004086741355 0.530523831813 0.917057469010
H -0.026379715087 0.529295070126 -0.917396616746
Br -2.526507180197 -0.186094710268 0.030826558248
Cl 2.378086639004 0.175176595492 -0.029016607261
10回程度の反復計算で遷移状態構造が得られた。-freqオプションにより生成されたnormal_modes.txtやvibration_animationディレクトリ内の振動モードのアニメーションを確認した。
以下に-freqオプションで生成されたnormal_modes.txtの一部を示す。
Mode 0 1 2
Freq [cm^-1] -327.4811 191.9089 191.9095
Reduced mass [au] 12.1303 4.1596 4.1596
Force const [Dyne/A] -0.7665 0.0903 0.0903
Char temp [K] 0.0000 276.1141 276.1150
Normal mode x y z x y z x y z
C -0.27121 -0.01998 0.00331 0.01734 -0.23721 -0.00847 -0.00350 0.00823 -0.23782
H -0.04388 -0.01269 0.00054 0.06049 -0.24045 -0.00914 -0.00508 0.00834 -0.24206
H -0.04482 0.00143 0.00871 -0.00505 -0.24370 -0.00933 -0.03977 0.00481 -0.24323
H -0.04502 0.00142 -0.00762 -0.00241 -0.24345 -0.00750 0.03412 0.01210 -0.24420
Br 0.02551 0.00188 -0.00031 -0.00145 0.02157 0.00077 0.00030 -0.00075 0.02161
Cl 0.03936 0.00290 -0.00048 -0.00419 0.05370 0.00192 0.00083 -0.00186 0.05386
(...snip...)その結果、虚振動が１つであることが確認できた。
次に、vibration_animation内のmode_1_327i_wave_number.xyzをAvogadroで確認すると、想定される反応系と生成系をつなぐ方向に振動していることを確認できた。
終わりに
自作モジュールで、電子状態計算モジュールであるPySCFを用いて、計算レベルHF/3-21Gでブロモメタンと塩化物イオンが関わるSN2の遷移状態構造を算出する手順を説明した。
参考
- 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.