Home
2196 words
11 minutes
【計算化学】自作pythonモジュールで遷移状態構造を求めてみる（分子内求核置換反応）
2025-06-02
Computational Chemistry
MultiOptPy
/
python

最終更新：2025-06-03

概要#

本記事では、自作モジュール(MultiOptPy)で、分子内求核置換反応(SNi反応)の遷移状態構造を算出してみる。電子状態計算ソフトウェアにPySCFを使用する。 MultiOptPyは電子状態計算ソフトウェアを用いた分子構造最適化手法の勉強を目的として作成したpythonモジュールである。

MultiOptPyのレポジトリ:https://github.com/ss0832/MultiOptPy

求核置換反応：

使用した自作モジュールMultiOptPyのバージョン#

1.9.9b

環境#

WSL2 (Ubuntu-22.04)

Source codeのダウンロード#

wget https://github.com/ss0832/MultiOptPy/archive/refs/tags/v1.9.9b.zip
unzip v1.9.9b.zip
cd MultiOptPy-1.9.9b

手順#

1. 初期構造の準備#

モデル反応系として、以下の構造を用意した。今回は名前をsni.xyzとした。

8

C     -0.238288037000     -1.130438310795     -0.268060261498
H      0.174608309256     -2.079640497297      0.023307938416
H      0.081491224063     -0.857235383546     -1.259869249380
H     -1.317151018947     -1.152064207442     -0.204398150975
O      0.319110662208     -0.183627321277      0.707651047767
S     -0.224535437451      1.381033139551      0.710269397007
O      0.654412001897      2.081982919584      1.761228418728
Cl     0.550352295973      1.939989661222     -1.470129140063
初期パスを求めるための初期構造

2. 遷移状態構造最適化のための初期構造の算出#

※あらかじめpythonで利用できる電子状態計算ソフトウェアのPySCF(導入方法はpipを使えば可能である。)が使える環境を用意しておく。

a. 初期構造をカレントディレクトリにsni.xyzとして保存し、以下を実行する。#

python optmain.py sni.xyz -opt rsirfo_fsb -func hf -bs 3-21g -pyscf -spin 0 -ma 200 1 8 -400 1 5 -ns 40
  • -opt rsirfo_fsbは準ニュートン法であるRS-I-RFO法を構造最適化に使用することを示す。初期のへシアンに関しては、特にオプションで指定しない限り、単位行列が使われる。
  • -pyscfは使用する電子状態計算ソフトウェアをPySCFにすることを示す。
  • -spinはスピン多重度の指定である。PySCFを使用するときは目的とするスピン多重度に1を引いた値を指定する。
  • -ma 200 1 8 -400 1 5は200kJ/molの活性化障壁を超えうるペア同士を近づける力を原子のラベル番号1と8のペアに、400kJ/molの活性化障壁を超えうるペア同士を遠ざける力を原子のラベル番号1と5のペアに、構造最適化時にそれぞれ加えることを示す。
  • -ns 40反復計算回数を40回に制限している。これは生成系中の２つの分子が人工力によって離れていき構造最適化が終わらないのでつけている。

これを実行すると、計算レベルHF/3-21Gで指定した人工力ポテンシャルをかけた上で初期構造を構造最適化することができる。

結果はyyyy_mm_dd（今日の年月日）ディレクトリの中に存在するディレクトリを開いて確認できる。

正常終了していれば、このディレクトリ中に、sni_traj.xyzが存在するので、これをコピーして、MultiOptPy-1.9.9bディレクトリに置く。

sni_traj.xyzは構造最適化の過程をAvogadro（公式ページ:https://avogadro.cc/ ）等で可視化して確認できるようにしている。このsni_traj.xyzは次のNEB計算に使用する。

※sni_traj.xyzをアニメーションとして表示したい場合は、[https://github.com/ss0832/molecule_movie] を使うと良い。

（この人工力ポテンシャルを加えて行った構造最適化の結果はsni_optimized.xyzで確認できる。構造を可視化して、生成系になっているか確認する。反応系のままであれば、-maの設定を見直してb.をやり直す。今回の場合以下の構造が得られた。中間体の安定構造に近いものが得られていると判断できるため、次のNEB計算を行うことが可能である。）

初期パスを求めるための構造最適化の結果（安定構造ではない）

b. sni_traj.xyzを初期パスとして、NEB法で経路の緩和を行う。#

NEB法を用いることで、先ほど得られたsni_traj.xyz全体のエネルギーを下げることができる。これにより、パスのエネルギー極大値を持つ構造を遷移状態構造に近づける。（この時点ではまだ正確な遷移状態構造は求められていない。）

LUP法と呼ばれるNEB法とよく似た方法があるが、本モジュールではNEB法を用いている。

python nebmain.py sni_traj.xyz -nd 0.4 -func hf -bs 3-21g -pyscf -spin 0 -spng -ns 20
  • -nd 0.4はノード間の距離を0.4Åとして初期パスを作成することを示す。
  • -ns 20は20回分NEB法による経路の緩和を行うことを示す。今回は200kJ/mol以上の強めの人工力を使用して初期パスを求めたため、10回では緩和が足りないと判断した。
  • -pyscfは使用する電子状態計算ソフトウェアをPySCFにすることを示す。
  • -spngは緩和中のパスのエネルギープロファイルや各ノードの勾配のRMS値をmatplotlibで可視化するオプションである。

(参考)

c. 初期構造の決定及び遷移状態構造の計算#

MultiOptPy-1.9.9bと同じディレクトリ内に、NEBという名前を含むディレクトリが生成されている。 そのディレクトリ内のenergy_plot.csvを確認し、緩和後のパスのエネルギー極大値を示す構造を確認する。

パスの緩和後の各ノードのエネルギー一覧(単位) (energy_plot.csvに保存されている。)

NEB計算の結果の可視化

私が実行した環境では、11番のノード(グラフでは12番)がエネルギー極大値を示していた。 path_ITR_20_sniディレクトリ内のsni_traj_11.xyzを可視化する。

※こちら[https://ss0832.github.io/molecule_viewer/] を使うことでも可視化は可能である。

sni_traj_11.xyz

8
0 0
C      -0.067423055121     -1.289336532977     -0.897093357450
H      -0.298744132598     -1.897623133986     -0.036128672547
H       0.956042041575     -1.005617800824     -1.013275152469
H      -0.850650002294     -0.744166928793     -1.360051491609
O       0.020418351335      0.029220616370      1.697169363068
S      -0.099142911765      1.535720821363      1.403767749258
O       1.263630276727      2.224240271391      1.400331252170
Cl     -0.924130567858      1.147562687455     -1.194719690421
NEB法により緩和したパスのエネルギー極大値を示す構造

構造が壊れていないので、これを遷移状態を求めるための初期構造とする。

sni_traj_11.xyzをMultiOptPy-1.9.9bと同じディレクトリ内にコピーする。

そして、以下を実行する。

python optmain.py sni_traj_11.xyz -opt rsirfo_bofill -func hf -bs 3-21g -pyscf -spin 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と同等）

実行して得られた正確な遷移状態構造を以下に示す。

(実行して得られた正確な遷移状態構造はsni_traj_11_optimized.xyzとして保存されている。)

8
OptimizedStructure
C     -0.097210835249     -1.115175616977     -0.686250494553
H     -0.320101743542     -1.757387960431      0.132320984320
H      0.917169448251     -0.832609948796     -0.852686767133
H     -0.772074932489     -1.078749669191     -1.509711666343
O      0.027243921035     -0.086881809201      1.225614196670
S     -0.076240607349      1.463840507703      1.402014822207
O      1.308679049962      2.099986510162      1.408413910538
Cl    -0.987464300618      1.306977986731     -1.119714985705
遷移状態構造

20回程度の反復計算で遷移状態構造が得られた。-freqオプションにより生成されたnormal_modes.txtやvibration_animationディレクトリ内の振動モードのアニメーションを確認した。

以下に-freqオプションで生成されたnormal_modes.txtの一部を示す。

Mode                                 0                   1                   2
Freq [cm^-1]                     -657.3659             83.9694            129.8648
Reduced mass [au]                   2.5282              4.8892              1.2782
Force const [Dyne/A]               -0.6437              0.0203              0.0127
Char temp [K]                       0.0000            120.8133            186.8465
Normal mode                   x         y         z            x         y         z            x         y         z     
       C                 0.05796   -0.04782    0.15397    0.14055    0.05558    0.00183   -0.01180   -0.03686    0.04035
       H                 0.24639   -0.32534   -0.01041    0.04816    0.01249   -0.05676   -0.50992   -0.00836   -0.07083
       H                 0.02509   -0.03663   -0.02534    0.15300    0.10590    0.16190    0.12919   -0.29758    0.46337
       H                -0.14313    0.22623    0.33308    0.26081    0.05564   -0.09648    0.32410    0.17245   -0.22439
       O                -0.01003   -0.04531   -0.09188   -0.11668    0.00089    0.04571    0.04665   -0.02070   -0.02380
       S                -0.00934   -0.00445   -0.02781   -0.00516    0.01427   -0.00177   -0.00802   -0.02077   -0.01718
       O                -0.00542   -0.01175   -0.00090    0.04419   -0.08993   -0.08814   -0.03263    0.03376   -0.03678
      Cl                -0.00798    0.05049    0.00646   -0.02367    0.00359    0.02015    0.00660    0.02951    0.02472

(...snip...)

その結果、虚振動が１つであることが確認できた。

次に、vibration_animation内のmode_1_657i_wave_number.xyzをAvogadroで確認すると、想定される反応系と生成系をつなぐ方向に振動していることを確認できた。

終わりに#

　 　自作モジュールで、電子状態計算モジュールであるPySCFを用いて、計算レベルHF/3-21Gで分子内求核置換反応(SNi反応)の遷移状態構造を算出する手順を説明した。

参考#

  • 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.
【計算化学】自作pythonモジュールで遷移状態構造を求めてみる（分子内求核置換反応）
https://ss0832.github.io/posts/20250602_mop_usage_12/
Author
ss0832
Published at
2025-06-02
【計算化学】自作pythonモジュールで遷移状態構造を求めてみる（Williamsonエーテル合成）
【計算化学】自作pythonモジュールで二次の鞍点構造を求めてみる（ホスフィン1分子）