最終更新：2025-05-19
概要
本記事では、自作モジュール(MultiOptPy)で、Pd錯体が関わる還元的脱離の遷移状態構造を算出してみる。今回はPd錯体にホスフィン配位子ではなく、二座ホスフィン配位子である1,2-ビス(ジフェニルホスフィノ)エチレン (DPPV)の置換基部位が水素の二座ホスフィン配位子を配位させた場合について求める。電子状態計算ソフトウェアにPySCFを使用する。 MultiOptPyは電子状態計算ソフトウェアを用いた分子構造最適化手法の勉強を目的として作成したpythonモジュールである。
MultiOptPyのレポジトリ:https://github.com/ss0832/MultiOptPy
使用した自作モジュールMultiOptPyのバージョン
1.9.6
環境
WSL2 (Ubuntu-22.04)
手順
1. 初期構造の準備
モデル反応系として、還元的脱離によりクロロメタンが生じるPd錯体をxyzファイル形式で用意した。今回は名前をdppv_reductive_elimination.xyzとした。
以下の初期構造を使用した。（求めたい計算レベルで構造最適化した構造を使用。）
16
Pd 0.690574651117 -0.467474892242 -0.929686552917
P -0.563925518225 -0.662787059995 1.051316678351
H -1.725450702919 -1.448471508237 1.064812481125
H 0.034554453275 -1.133611852994 2.229292119826
P 0.090183759522 1.991385284605 -0.688064145578
H 1.049370193987 2.989129819866 -0.455964827041
H -0.724755586321 2.679446000401 -1.599855346143
C -1.191390172412 0.984915529424 1.553494521126
H -1.796075740785 1.058580794020 2.451710634718
C -0.921726126253 2.078242344180 0.846904809564
H -1.310029427184 3.033657806583 1.184440819179
C 1.046967157040 -2.473394285785 -0.909972632465
H 0.107836736017 -3.030261349280 -0.969710188398
H 1.661283328052 -2.714136060905 -1.771794505942
H 1.575098859544 -2.760875902045 0.003433834399
Br 1.977484135545 -0.124344667596 -3.060357699806
2. 遷移状態構造最適化のための初期構造の算出
※あらかじめpythonで利用できる電子状態計算ソフトウェアのPySCF(導入方法はpipを使えば可能である。)が使える環境を用意しておく。
a. 以下のURLにアクセスし、Source codeをダウンロードする。
wget https://github.com/ss0832/MultiOptPy/archive/refs/tags/v1.9.6.zip
unzip v1.9.6.zip
cd MultiOptPy-1.9.6
b. 初期構造をカレントディレクトリにdppv_reductive_elimination.xyzとして保存し、以下を実行する。
python optmain.py dppv_reductive_elimination.xyz -opt rsirfo_fsb -pyscf -spin 0 -bs def2svp -func hf -ka 1.0 30 12,1,16 -ns 50 -ecp default def2svp
-opt rsirfo_fsbは準ニュートン法であるRS-I-RFO法を構造最適化に使用することを示す。初期のへシアンに関しては、特にオプションで指定しない限り、単位行列が使われる。-bsは基底関数の指定（特に指定しなければ6-31G(d)が使われる。）今回はKarlsruhe系の基底関数であるdef2-SVPを使用。-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 50は、構造最適化時の電子状態計算ソフトウェアの計算の反復回数の上限を50回とすることを示す。（デフォルトでは1000回）今回は初期パスを求めることが目的である。そのため、構造最適化で収束させるまでの時間コストを削減するためにデフォルトよりも少なくした。-ecpは有効内核ポテンシャルの指定を行える。今回は、Pdに対してdef2-SVPの有効内核ポテンシャルを割り当てている。
これを実行すると、計算レベルHF/def2-SVPで指定した人工力ポテンシャルをかけた上で初期構造を構造最適化する。
結果はyyyy_mm_dd（今日の年月日）ディレクトリの中に存在するディレクトリを開いて確認できる。
正常終了していれば、このディレクトリ中に、dppv_reductive_elimination_traj.xyzが存在するので、これをコピーして、MultiOptPy-1.9.6ディレクトリに置く。
dppv_reductive_elimination_traj.xyzは構造最適化の過程をAvogadro（公式ページ:https://avogadro.cc/ ）等で可視化して確認できるようにしている。このdppv_reductive_elimination_traj.xyzは次のNEB計算に使用する。
※dppv_reductive_elimination_traj.xyzをアニメーションとして表示したい場合は、[https://github.com/ss0832/molecule_movie] を使うと良い。
（この人工力ポテンシャルを加えて行った構造最適化の結果はdppv_reductive_elimination_optimized.xyzで確認できる。構造を可視化して、生成系になっているか確認する。反応系のままであれば、-maの設定を見直してb.をやり直す。）
c. dppv_reductive_elimination_traj.xyzを初期パスとして、NEB法で経路の緩和を行う。
NEB法を用いることで、先ほど得られたdppv_reductive_elimination_traj.xyz全体のエネルギーを下げることができる。これにより、パスのエネルギー極大値を持つ構造を遷移状態構造に近づける。（この時点ではまだ正確な遷移状態構造は求められていない。）
LUP法と呼ばれるNEB法とよく似た方法があるが、本モジュールではNEB法を用いている。
python nebmain.py dppv_reductive_elimination_traj.xyz -pyscf -spin 0 -bs def2svp -func hf -ns 10 -nd 0.25 -ecp default def2svp
-nd 0.25はノード間の距離を0.25Åとして初期パスを作成することを示す。-ns 10は10回分NEB法による経路の緩和を行うことを示す。-bsは基底関数の指定（特に指定しなければ6-31G(d)が使われる。）今回はKarlsruhe系の基底関数であるdef2-SVPを使用。-funcは計算手法の指定（Hartree-Fock(HF), post-HF(MP2等), DFTの汎関数）が可能。（特に指定しなければB3LYPが使われる。）今回はHartree-Fock法を使用。-pyscfは使用する電子状態計算ソフトウェアをPySCFにすることを示す。-spinはスピン多重度の指定である。PySCFを使用するときは目的とするスピン多重度に1を引いた値を指定する。-ecpは有効内核ポテンシャルの指定を行える。今回は、Pdに対してdef2-SVPの有効内核ポテンシャルを割り当てている。
(参考)
d. 初期構造の決定及び遷移状態構造の計算
MultiOptPy-1.9.6と同じディレクトリ内に、NEBという名前を含むディレクトリが生成されている。 そのディレクトリ内のenergy_plot.csvを確認し、緩和後のパスのエネルギー極大値を示す構造を確認する。
<パスの緩和後の各ノードのエネルギー一覧(単位
-3499.1277851525 -3499.1277847142 -3499.1151714093 -3499.0977944898 -3499.0851151491 -3499.0671648915 -3499.0323867988 -3499.032528 -3499.0236940477 -3499.0212092458 -3499.0391832949 -3499.053215144 -3499.0596235272 -3499.0750396029 -3499.0819618187 -3499.0863139933 -3499.0892642792 -3499.0919380625 -3499.0931769624 -3499.0948141609 -3499.0960983343 -3499.0970924158 -3499.097682317
私が実行した環境では、9番である。 path_ITR_10_dppv_reductive_eliminationディレクトリ内のdppv_reductive_elimination_traj_9.xyzを可視化する。
※こちら[https://ss0832.github.io/molecule_viewer/] を使うことでも可視化は可能である。
dppv_reductive_elimination_traj_9.xyz
16
0 0
Pd 0.613539493170 -0.257093487413 -0.858533053142
P -0.634851818421 -0.638062080771 1.150298002099
H -1.803575045235 -1.430938186419 1.200897358306
H -0.061493681792 -1.118499647730 2.351642591331
P 0.042676436629 2.061149286563 -0.629451576586
H 1.000933862108 3.075508321148 -0.403062539749
H -0.768758182546 2.758962704857 -1.552328206857
C -1.281385504362 1.052463549911 1.675936088917
H -1.898837633636 1.193315634792 2.581368797866
C -0.998386708658 2.139264156975 0.946755889015
H -1.394731166617 3.098040055897 1.291534921095
C 1.281893362878 -2.494845669700 -1.271523365625
H 0.224476010861 -2.712382412471 -1.337830424992
H 1.838172710215 -3.067981255572 -1.934286466212
H 1.765215042469 -2.424773291711 -0.305323975676
Br 2.075112822937 -1.234127678355 -2.906094039789
構造が壊れていないので、これを遷移状態を求めるための初期構造とする。
dppv_reductive_elimination_traj_9.xyzをMultiOptPy-1.9.6と同じディレクトリ内にコピーする。
そして、以下を実行する。
python optmain.py dppv_reductive_elimination_traj_9.xyz -opt rsirfo_bofill -order 1 -fc 5 -pyscf -spin 0 -bs def2svp -func hf -tcc -freq -ecp default def2svp
-opt rsirfo_bofillは遷移状態構造の最適化向けのoptimizerを指定することを意味する。準ニュートン法であるRS-I-RFO法を使用する。今回は-fcで正確なへシアンを計算するようにしているので、初期へシアンは正確なへシアンを使用するようになっている。-order 1は一次の鞍点を求めることを指定する。（デフォルトだと極小値を求めるようになっている。）-fc 5は5回の反復回数当たり1回正確なへシアンを計算することを指定する。-freqは収束条件を満たした後に基準振動解析を行うことを示す。（自前で実装しているため、あくまで目安として使用することを推奨する。各振動モードをvibration_animation内のxyzファイルで可視化できる。）-tccはデフォルトよりも厳しい収束条件を課すことを示す。（電子状態計算ソフトウェアGaussianの収束条件のtight相当である。）
実行して得られた正確な遷移状態構造を以下に示す。
二座ホスフィン配位子の両方のドナー原子が金属に配位した構造を想定していたが、今回の計算レベルHF/def2-SVPの計算では、片方のドナー原子が配位した状態の遷移状態構造が算出された。
(実行して得られた正確な遷移状態構造はdppv_reductive_elimination_traj_9_optimized.xyzとして保存されている。)
16
OptimizedStructure
Pd 1.586795920124 -1.033121787256 0.153685895829
P 0.016844447105 -0.332973970769 1.819487467936
H -0.738798401675 -1.345057357063 2.439007059604
H 0.562250971480 0.204274661913 3.001388580846
P -0.380042846739 1.833041034426 -1.027522928531
H 0.852559538146 2.173000656795 -0.439784838840
H -0.754780955292 3.154035405135 -1.351262325139
C -1.327509355342 0.873553976761 1.517624755097
H -2.101242955727 0.919847680524 2.277973732011
C -1.440646078545 1.686293067106 0.472908379297
H -2.318845312987 2.325510893376 0.447850343145
C 0.858329614683 -2.118054047422 -1.851886033803
H 0.459474093032 -1.189444194519 -2.233816296619
H 1.231741140972 -2.787214498127 -2.606754329414
H 0.260185744843 -2.618280196043 -1.101250222649
Br 3.233684435925 -1.745411324835 -1.517649238770
200回程度の反復計算で遷移状態構造が得られた。-freqオプションにより生成されたnormal_modes.txtやvibration_animationディレクトリ内の振動モードのアニメーションを確認した。
以下に-freqオプションで生成されたnormal_modes.txtの一部を示す。
Mode 0 1 2
Freq [cm^-1] -443.7014 23.5802 29.1995
Reduced mass [au] 4.0743 3.6549 3.7355
Force const [Dyne/A] -0.4726 0.0012 0.0019
Char temp [K] 0.0000 33.9266 42.0115
Normal mode x y z x y z x y z
Pd 0.00834 -0.01496 -0.02893 -0.01540 -0.00661 0.00027 -0.00239 -0.02213 0.02777
P -0.01262 -0.00133 -0.00026 -0.03407 -0.03965 -0.00620 0.01052 0.01850 0.02464
H 0.00650 -0.00519 0.01101 -0.12416 -0.03066 -0.09979 0.03459 0.03821 0.08554
H 0.03934 -0.00505 -0.02404 -0.03626 -0.15960 0.04885 0.00110 0.10093 -0.00821
P -0.00099 0.00140 -0.00087 0.01701 -0.05094 -0.02680 0.07441 -0.01794 -0.01582
H -0.00026 -0.00807 0.00460 0.05701 -0.07084 -0.09923 0.04932 -0.02900 0.04271
H 0.00656 0.00486 0.00852 0.05590 -0.05197 -0.07623 0.08761 -0.01136 -0.00544
C -0.00238 0.00152 0.00189 0.06138 0.07755 0.03937 -0.01753 -0.03217 -0.05255
H -0.00142 -0.00027 0.00231 0.10361 0.15938 0.07758 -0.05781 -0.05179 -0.09247
C 0.00093 0.00084 0.00068 0.08012 0.06996 0.03138 0.00785 -0.04592 -0.06620
H -0.00012 -0.00078 0.00185 0.13466 0.14616 0.06697 -0.01302 -0.07646 -0.11689
C -0.17656 0.04491 0.12689 0.01379 -0.09366 0.03840 -0.01231 -0.01172 0.01817
H -0.14201 0.00878 0.00184 -0.07107 -0.12032 0.06171 0.12058 -0.03599 -0.17651
H 0.02687 0.04608 0.23006 0.04163 -0.05401 0.01707 -0.03043 -0.19829 0.17401
H -0.33797 0.02732 -0.01227 0.07883 -0.15045 0.05222 -0.13184 0.19609 0.06254
Br 0.02636 0.01200 0.01673 0.00068 0.04047 -0.00464 -0.02755 0.04399 -0.02500
(...snip...)
その結果、虚振動が１つであることが確認できた。
次に、vibration_animation内のmode_1_444i_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.