【計算化学】自作pythonモジュールで遷移状態構造を求めてみる(Kolbe-Schmitt反応1）

最終更新：2025-07-02

概要#

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

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

Kolbe-Schmitt反応:

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

1.9.12c

環境#

WSL2 (Ubuntu-22.04)

Source codeのダウンロード#

wget https://github.com/ss0832/MultiOptPy/archive/refs/tags/v1.9.12c.zip
unzip v1.9.12c.zip
cd MultiOptPy-1.9.12c

手順#

1. 初期構造の準備#

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

16
0 1
 C                 -1.26394156   -1.26459405    0.00008088
 C                  0.13745828   -1.26476121    0.00073857
 C                  0.83830337   -0.05119711    0.00019378
 C                  0.13774862    1.16253415   -0.00100884
 C                 -1.26365122    1.16270130   -0.00166709
 C                 -1.96449631   -0.05086280   -0.00112112
 H                 -1.79905234   -2.19117717    0.00049642
 H                  0.67234737   -2.19147196    0.00165655
 H                  1.90830324   -0.05132473    0.00069749
 H                  0.67285940    2.08911727   -0.00142434
 H                 -3.03449618   -0.05073517   -0.00162437
 O                 -1.97850300    2.40120258   -0.00289588
 Na                -2.32085306    2.99451365   -1.94567184
 C                 -1.40601036    0.68682475   -3.79029766
 O                 -0.14761036    0.68682475   -3.79029766
 O                 -2.66441036    0.68682475   -3.79029766

初期パスを求めるための初期構造

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

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

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

python optmain.py kolbe_schmitt_rxn.xyz -pyscf -spin 0 -elec 0 -func hf -bs 3-21g -opt rsirfo_fsb -ma 200 13 16 200 6 14

-opt rsirfo_fsbは準ニュートン法であるRS-I-RFO法を構造最適化に使用することを示す。初期のへシアンに関しては、特にオプションで指定しない限り、単位行列が使われる。
-pyscfは使用する電子状態計算ソフトウェアをPySCFにすることを示す。
-spinはスピン多重度の指定である。PySCFを使用するときは目的とするスピン多重度に1を引いた値を指定する。
-elec 1は形式電荷を1とすることを示す。
-ma yyy a bはyyykJ/molの活性化障壁を超えうるペア同士を近づける力を原子のラベル番号aとbのペアに構造最適化時に加えることを示す。
-func hfはHartree-Fock法を用いることを示す。
-bs 3-21gは基底関数にPople系である3-21Gを用いることを示す。

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

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

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

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

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

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

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

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

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

python nebmain.py kolbe_schmitt_rxn_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 nはn回分NEB法による経路の緩和を行うことを示す。
-modelhessは少ない計算コストで算出された近似へシアンを用いることを示す。これを指定すると経路緩和のアルゴリズムが-fcを1以上に指定した場合と同様のものになる。
-pyscfは使用する電子状態計算ソフトウェアをPySCFにすることを示す。
-spngは緩和中のパスのエネルギープロファイルや各ノードの勾配のRMS値をmatplotlibで可視化するオプションである。
-spinはスピン多重度の指定である。PySCFを使用するときは目的とするスピン多重度に1を引いた値を指定する。

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

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

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

NEB計算の結果の可視化

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

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

kolbe_schmitt_rxn_traj_8.xyz

16
0 0
C      -0.332351743236     -1.508435596519      0.465903485973
C       0.973532901636     -1.545587682967      0.720228743669
C       1.701828209580     -0.308226145510      0.739348795384
C       1.088533018211      0.878994294703      0.673726999682
C      -0.329927905533      0.979042820335      0.467121356091
C      -0.965209756808     -0.275949018185      0.059453862972
H      -0.947752483635     -2.375752195039      0.546269612642
H       1.483686884797     -2.453247841409      0.953024958359
H       2.748591573213     -0.347977961928      0.944253863990
H       1.593838776493      1.799419225791      0.834649752380
H      -2.038243304772     -0.258842439849      0.097280996659
O      -0.935952110095      2.069653266643      0.460897553869
Na     -2.389954202773      2.535546989347     -1.087383311777
C      -0.594552821323      0.130596988789     -1.668748859124
O       0.319375223653     -0.379748036715     -2.259944158053
O      -1.375442259407      1.060513332515     -1.946083652716

NEB法により緩和したパスのエネルギー極大値を示す構造

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

kolbe_schmitt_rxn_traj_8.xyzをMultiOptPy-1.9.12cと同じディレクトリ内にコピーする。

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

python optmain.py kolbe_schmitt_rxn_traj_8.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と同等）

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

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

16
OptimizedStructure
C     -0.212263660485     -1.470440765436      0.346991513291
C      1.080078549428     -1.484103074672      0.722372074566
C      1.733804620475     -0.240158083678      0.976165575225
C      1.088856290735      0.939715619999      0.875680597397
C     -0.305012048486      1.003042686100      0.517989509837
C     -0.915168664583     -0.241869895379      0.129438346863
H     -0.727169247008     -2.388424328032      0.142698073132
H      1.622602269622     -2.400526137530      0.830123559755
H      2.767082626127     -0.254234520804      1.266721369127
H      1.576431283112      1.865169183510      1.105041608659
H     -1.991301160154     -0.269197976501      0.177149664245
O     -0.929939482223      2.108323319616      0.487357045699
Na    -1.987099582063      2.790793983439     -1.138364617510
C     -0.907814822271     -0.022342865672     -1.905222291481
O     -0.401938042587     -0.959041439607     -2.414531842469
O     -1.491148929637      1.023294294648     -2.119610186337

遷移状態構造

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

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

Mode                                 0                   1                   2
Freq [cm^-1]                     -403.0109             47.2205             88.4808
Reduced mass [au]                  11.0091              8.4497              7.7605
Force const [Dyne/A]               -1.0535              0.0111              0.0358
Char temp [K]                       0.0000             67.9398            127.3042
Normal mode                   x         y         z            x         y         z            x         y         z
       C                -0.01531    0.01534    0.04710    0.06431   -0.00359   -0.03759    0.01197    0.01883   -0.02288
       C                 0.02529    0.00670   -0.00321    0.05700    0.02938   -0.00821   -0.01006   -0.01234    0.05170
       C                 0.00815   -0.01804   -0.00597    0.00814    0.04379    0.04394    0.00057   -0.02867    0.10698
       C                -0.01198   -0.00936   -0.00596   -0.02945    0.02444    0.05645    0.03474   -0.01336    0.07022
       C                -0.00268   -0.01731    0.02846   -0.02204   -0.01282    0.02118    0.05854    0.02129   -0.01728
       C                -0.00328   -0.01438    0.18245    0.02376   -0.02225   -0.01258    0.04372    0.03462   -0.03955
       H                -0.01971    0.01250    0.06123    0.10168   -0.01471   -0.08024    0.00384    0.03168   -0.06023
       H                 0.02304    0.00001   -0.04442    0.08531    0.04432   -0.02332   -0.03583   -0.02512    0.07336
       H                 0.01198   -0.00412   -0.01677    0.00095    0.07110    0.07108   -0.01934   -0.05447    0.17670
       H                -0.00671   -0.00856   -0.01625   -0.06656    0.03532    0.09178    0.04081   -0.02518    0.10418
       H                -0.00982   -0.00486   -0.04211    0.02362   -0.05977   -0.01648    0.04205    0.05184   -0.04764
       O                -0.01339    0.01361   -0.01796   -0.05304   -0.02946    0.02775    0.08665    0.03640   -0.04706
      Na                 0.00221    0.01000    0.02104    0.05522    0.02859   -0.01209   -0.08457   -0.04351    0.03193
       C                 0.00239    0.02624   -0.19808   -0.03273   -0.01916   -0.02606   -0.01466    0.00949   -0.04987
       O                 0.01488   -0.01701   -0.01941   -0.14384   -0.07696   -0.02920   -0.05273   -0.00935   -0.05138
       O                -0.00652   -0.00256   -0.02280    0.05660    0.03068   -0.01173   -0.00800    0.01442   -0.03748

(...snip...)

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

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

終わりに#

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

参考#

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.

概要#