最終更新：2025-05-18
概要
本記事では、自作モジュール(MultiOptPy)で、Pd錯体が関わる還元的脱離の遷移状態構造を算出してみる。電子状態計算ソフトウェアにPySCFを使用する。 MultiOptPyは電子状態計算ソフトウェアを用いた分子構造最適化手法の勉強を目的として作成したpythonモジュールである。
MultiOptPyのレポジトリ:https://github.com/ss0832/MultiOptPy
使用した自作モジュールMultiOptPyのバージョン
1.9.5b
環境
WSL2 (Ubuntu-22.04)
手順
1. 初期構造の準備
モデル反応系として、還元的脱離によりクロロメタンが生じるPd錯体をxyzファイル形式で用意した。今回は名前をreductive_elimination.xyzとした。
以下の初期構造を使用した。（求めたい計算レベルで構造最適化した構造を使用。）
10
Pd 0.964650052418 0.891771433775 -0.009519082792
P -1.546431504098 0.909920579639 0.060646416449
H -2.240595368384 2.021899687438 0.584472129325
H -2.211050355998 0.752916176074 -1.174497966093
H -2.132439261334 -0.135998111616 0.803269457853
Cl 3.347143445207 1.041725005112 -0.071422638318
C 0.948661078845 -1.129406998296 -0.049088240278
H 0.736374718077 -1.456905031744 0.959297763472
H 0.197156870856 -1.472928521361 -0.746569121481
H 1.936530324410 -1.422994219020 -0.356588718137
2. 遷移状態構造最適化のための初期構造の算出
※あらかじめpythonで利用できる電子状態計算ソフトウェアのPySCF(導入方法はpipを使えば可能である。)が使える環境を用意しておく。
a. 以下のURLにアクセスし、Source codeをダウンロードする。
wget https://github.com/ss0832/MultiOptPy/archive/refs/tags/v1.9.5b.zip
unzip v1.9.5b.zip
cd MultiOptPy-1.9.5b
b. 初期構造をカレントディレクトリにreductive_elimination.xyzとして保存し、以下を実行する。
python optmain.py reductive_elimination.xyz -opt rsirfo_fsb -pyscf -spin 0 -bs 6-31g -sub_bs Pd LanL2DZ -func hf -ka 1.0 30 6,1,7 -ns 100
-opt rsirfo_fsbは準ニュートン法であるRS-I-RFO法を構造最適化に使用することを示す。初期のへシアンに関しては、特にオプションで指定しない限り、単位行列が使われる。-bsは基底関数の指定（特に指定しなければ6-31G(d)が使われる。）今回はPople系の基底関数である6-31Gを使用。-sub_bsは二種類以上の基底関数とそれに対応するECPを指定するために使用する。今回は、6-31GにPdが対応していないため、代わりにLanL2DZを割り当てるために使用した。割り当てたい元素の後に、割り当てたい基底関数名を指定する。-funcは計算手法の指定（Hartree-Fock(HF), post-HF(MP2等), DFTの汎関数）が可能。（特に指定しなければB3LYPが使われる。）今回はHartree-Fock法を使用。-pyscfは使用する電子状態計算ソフトウェアをPySCFにすることを示す。-spinはスピン多重度の指定である。PySCFを使用するときは目的とするスピン多重度に1を引いた値を指定する。-kaは指定した３つの原子がなす角度に対して人工的な調和ポテンシャルを付与することを指定する。初めに力の定数kを原子単位系で指定する。次に、平衡となる角度Θをdegree単位で指定する。最後に、ポテンシャルを付与する対象a,b,cを指定する。これにより指定した3津の原子がなす角度∠a-b-cに対して力の定数k、平衡値Θの調和ポテンシャルが構造最適化時に付与される。-ns 100は、構造最適化時の電子状態計算ソフトウェアの計算の反復回数の上限を100回とすることを示す。（デフォルトでは1000回）今回は初期パスを求めることが目的である。そのため、構造最適化で収束させるまでの時間コストを削減するためにデフォルトよりも少なくした。
これを実行すると、計算レベルHF/6-31G (LanL2DZ for Pd)で指定した人工力ポテンシャルをかけた上で初期構造を構造最適化する。
結果はyyyy_mm_dd（今日の年月日）ディレクトリの中に存在するディレクトリを開いて確認できる。
正常終了していれば、このディレクトリ中に、reductive_elimination_traj.xyzが存在するので、これをコピーして、MultiOptPy-1.9.5bディレクトリに置く。
reductive_elimination_traj.xyzは構造最適化の過程をAvogadro（公式ページ:https://avogadro.cc/ ）等で可視化して確認できるようにしている。このreductive_elimination_traj.xyzは次のNEB計算に使用する。
※reductive_elimination_traj.xyzをアニメーションとして表示したい場合は、[https://github.com/ss0832/molecule_movie] を使うと良い。
（この人工力ポテンシャルを加えて行った構造最適化の結果はreductive_elimination_optimized.xyzで確認できる。構造を可視化して、生成系になっているか確認する。反応系のままであれば、-maの設定を見直してb.をやり直す。）
c. reductive_elimination_traj.xyzを初期パスとして、NEB法で経路の緩和を行う。
NEB法を用いることで、先ほど得られたreductive_elimination_traj.xyz全体のエネルギーを下げることができる。これにより、パスのエネルギー極大値を持つ構造を遷移状態構造に近づける。（この時点ではまだ正確な遷移状態構造は求められていない。）
LUP法と呼ばれるNEB法とよく似た方法があるが、本モジュールではNEB法を用いている。
python nebmain.py reductive_elimination_traj.xyz -pyscf -spin 0 -bs 6-31g -sub_bs Pd LanL2DZ -func hf -ns 10 -nd 0.2
-nd 0.2はノード間の距離を0.20Åとして初期パスを作成することを示す。-ns 10は10回分NEB法による経路の緩和を行うことを示す。-bsは基底関数の指定（特に指定しなければ6-31G(d)が使われる。）今回はPople系の基底関数である6-31Gを使用。-sub_bsは二種類以上の基底関数とそれに対応するECPを指定するために使用する。今回は、6-31GにPdが対応していないため、代わりにLanL2DZを割り当てるために使用した。割り当てたい元素の後に、割り当てたい基底関数名を指定する。-funcは計算手法の指定（Hartree-Fock(HF), post-HF(MP2等), DFTの汎関数）が可能。（特に指定しなければB3LYPが使われる。）今回はHartree-Fock法を使用。-pyscfは使用する電子状態計算ソフトウェアをPySCFにすることを示す。-spinはスピン多重度の指定である。PySCFを使用するときは目的とするスピン多重度に1を引いた値を指定する。
(参考)
d. 初期構造の決定及び遷移状態構造の計算
MultiOptPy-1.9.5bと同じディレクトリ内に、NEBという名前を含むディレクトリが生成されている。 そのディレクトリ内のenergy_plot.csvを確認し、緩和後のパスのエネルギー極大値を示す構造を確認する。
<パスの緩和後の各ノードのエネルギー一覧(単位
-967.3533776682,-967.3533766816091,-967.3468252258808,-967.3301079444882,-967.310561093444,-967.2950880241922,-967.2828633901217,-967.2656894837386,-967.2504592121413,-967.2497427001297,-967.263114953713,-967.2788523768011,-967.2973752348205,-967.3101814527679,-967.3167743726967,-967.3233900059749,-967.3276802894632,-967.3302307445492,-967.3320222190255,-967.3332698723179,-967.334021212773,-967.3348051446567,-967.335271914242
私が実行した環境では、9番である。 path_ITR_10_reductive_eliminationディレクトリ内のreductive_elimination_traj_9.xyzを可視化する。
※こちら[https://ss0832.github.io/molecule_viewer/] を使うことでも可視化は可能である。
reductive_elimination_traj_9.xyz
10
0 0
Pd 0.668724799811 1.392434289376 -0.013679976681
P -1.691458516629 0.936133077128 0.061664340658
H -2.378186295672 2.034839070313 0.598599403130
H -2.311523643417 0.724289078021 -1.197082267894
H -2.197021118777 -0.165166673408 0.819068156917
Cl 2.927512660866 0.173697087633 -0.016038619634
C 1.298535396011 -1.144009914200 -0.058259918459
H 1.050957732908 -1.244971799254 0.968326378081
H 0.471028078303 -1.186225729080 -0.740685036524
H 2.161430906596 -1.521018486527 -0.421912459595
構造が壊れていないので、これを遷移状態を求めるための初期構造とする。
reductive_elimination_traj_9.xyzをMultiOptPy-1.9.5bと同じディレクトリ内にコピーする。
そして、以下を実行する。
python optmain.py reductive_elimination_traj_9.xyz -opt rsirfo_bofill -order 1 -fc 5 -pyscf -spin 0 -bs 6-31g -sub_bs Pd LanL2DZ -func hf -tcc -freq
-opt rsirfo_bofillは遷移状態構造の最適化向けのoptimizerを指定することを意味する。準ニュートン法であるRS-I-RFO法を使用する。今回は-fcで正確なへシアンを計算するようにしているので、初期へシアンは正確なへシアンを使用するようになっている。-order 1は一次の鞍点を求めることを指定する。（デフォルトだと極小値を求めるようになっている。）-fc 5は5回の反復回数当たり1回正確なへシアンを計算することを指定する。-freqは収束条件を満たした後に基準振動解析を行うことを示す。（自前で実装しているため、あくまで目安として使用することを推奨する。各振動モードをvibration_animation内のxyzファイルで可視化できる。）-tccはデフォルトよりも厳しい収束条件を課すことを示す。（電子状態計算ソフトウェアGaussianの収束条件のtight相当である。）
実行して得られた正確な遷移状態構造を以下に示す。
(実行して得られた正確な遷移状態構造はreductive_elimination_traj_9_optimized.xyzとして保存されている。)
10
OptimizedStructure
Pd 0.527694499687 0.859219291911 -0.261726306988
P -1.948272534942 1.071495592831 -0.005663592215
H -2.443896999758 1.902077003986 1.030374777658
H -2.689904853830 1.617715972478 -1.082952808137
H -2.758731266829 -0.066255220226 0.241767247715
Cl 2.952275870334 0.542827257817 -0.425173342926
C 1.618875230443 -1.322361830385 0.084638230109
H 1.521927126939 -1.233111697004 1.146365593187
H 0.727334738594 -1.543829899838 -0.468517584314
H 2.492698189361 -1.827776471571 -0.259112214089
70回程度の反復計算で遷移状態構造が得られた。-freqオプションにより生成されたnormal_modes.txtやvibration_animationディレクトリ内の振動モードのアニメーションを確認した。
以下に-freqオプションで生成されたnormal_modes.txtの一部を示す。
Mode 0 1 2
Freq [cm^-1] -486.7707 18.4676 48.4120
Reduced mass [au] 3.6509 1.0824 1.6216
Force const [Dyne/A] -0.5097 0.0002 0.0022
Char temp [K] 0.0000 26.5708 69.6541
Normal mode x y z x y z x y z
Pd -0.02547 0.01754 -0.00027 0.00059 0.00060 0.00256 0.01071 0.03152 0.00740
P 0.00430 0.00558 -0.00123 -0.00049 -0.00311 -0.00855 0.00692 -0.06068 -0.00158
H -0.02143 -0.00301 -0.00204 -0.03165 -0.43307 0.32156 -0.04271 -0.05538 -0.02970
H -0.01932 -0.00142 0.00714 0.09218 0.49762 0.18190 -0.02218 -0.13877 -0.02103
H -0.00212 0.00460 -0.00079 -0.06487 -0.07626 -0.55579 0.07288 -0.09984 0.03416
Cl -0.00678 -0.06864 0.01114 -0.00232 -0.00446 -0.02550 -0.00079 -0.03588 -0.02328
C 0.20164 0.02092 -0.02038 0.00248 0.01296 0.06093 -0.08619 0.00590 0.00679
H 0.05266 0.06310 -0.03859 -0.01830 0.05310 0.05584 -0.50690 0.05785 -0.03430
H 0.23477 0.25224 -0.16447 0.01493 -0.01263 0.05134 0.14565 -0.05006 -0.34226
H 0.13380 -0.19775 0.12057 0.01204 0.00425 0.09826 0.06881 0.01408 0.39087
(...snip...)
その結果、虚振動が１つであることが確認できた。
次に、vibration_animation内のmode_1_487i_wave_number.xyzをAvogadroで確認すると、還元的脱離で想定されるような反応系と生成系をつなぐ方向に振動していることを確認できた。
終わりに
自作モジュールで、遷移状態構造を算出する手順を説明した。
参考
- 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.