【計算化学】自作pythonモジュールで遷移状態構造を求めてみる（プロトンリレーありのホルムアルデヒドの水和反応）

最終更新：2025-04-28

概要#

本記事では、自作モジュール(MultiOptPy)で、ホルムアルデヒドの水和反応（プロトンリレーあり）の遷移状態構造を算出してみる。 MultiOptPyは電子状態計算ソフトウェアを用いた分子構造最適化手法の勉強を目的として作成したpythonモジュールである。プロトンリレーがある遷移状態構造は、環の歪みが4員環よりも少ない6員環遷移状態をとるため、プロトンリレーなしよりも活性化障壁が低くなる。したがって、似た素過程の遷移状態を求めるとき（例えばアセタール化等）、基質以外の反応にかかわる分子を2分子として、求めるとより現実的（実際にその形をとる可能性の高い）遷移状態構造になると考えられる。

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

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

1.8.1c

環境#

WSL2 (Ubuntu-22.04)

手順#

1. 初期構造の準備#

モデル反応系として、以下の構造をxyzファイル形式で用意した。今回は名前をhydration_2.xyzとした。この構造はホルムアルデヒドと水の２分子で構成されている。

(事前に構造最適化したものを使用すると失敗しにくい。)

10

C     -0.484976158346     -0.083814372337      1.479580729356
O      0.029569254001      1.003296662380      1.296465079297
O     -0.845814366047     -0.938175060839     -0.828994131148
H     -0.059072720380     -0.438531280996     -1.136993448093
H     -1.604922643547     -0.678084770966     -1.366065497740
H      0.069770304261     -1.005746444793      1.455746792313
H     -1.530798919018     -0.190793874043      1.723207701785
O      1.231153754004      0.725739427946     -1.148745688119
H      2.169231659043      0.508328915606     -1.202301945550
H      1.025859836029      1.097780798042     -0.271899592102

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

※あらかじめpythonで利用できる電子状態計算ソフトウェアのpsi4(導入方法は、psi4公式ページ [https://psicode.org/installs/v191/] を参照)が使える環境を用意しておく。

a. 以下のURLにアクセスし、Source codeをダウンロードする。#

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

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

python optmain.py hydration_2.xyz -ma 100 1 3 100 4 8 100 2 10 -bs 3-21g -func hf

-ma 100 1 3 100 4 8 100 2 10は原子のラベル番号1と3のペアと原子のラベル番号4と8のペア,原子のラベル番号2と10のペアに、100kJ/molの活性化障壁を超えうる人工力ポテンシャルをそれぞれ付加することを示す。(詳しくはAFIR法 [https://doi.org/10.1063/1.3457903] を参照)
-bsは基底関数の指定（特に指定しなければ6-31G(d)が使われる。）今回は3-21Gを使用した。
-funcは計算手法の指定（Hartree-Fock(HF), post-HF(MP2等), DFTの汎関数）が可能。（特に指定しなければB3LYPが使われる。）今回はHartree-Fock法を使用した

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

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

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

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

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

（この人工力ポテンシャルを加えて行った構造最適化の結果はhydration_2_optimized.xyzで確認できる。構造を可視化して、生成系になっているか確認する。反応系のままであれば、-maの設定を見直してb.をやり直す。）

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

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

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

python nebmain.py hydration_2_traj.xyz -ns 10 -spng -nd 0.15 -bs 3-21g -func hf

-spngは経路の緩和途中の各ノードのエネルギーや勾配のRMSを可視化する。
-nd 0.15はノード間の距離を0.15Åとして初期パスを作成することを示す。
-ns 10は10回分NEB法による経路の緩和を行うことを示す。
-bsは基底関数の指定（特に指定しなければ6-31G(d)が使われる。）今回は3-21Gを使用。
-funcは計算手法の指定（Hartree-Fock(HF), post-HF(MP2等), DFTの汎関数）が可能。（特に指定しなければB3LYPが使われる。）今回はHartree-Fock法を使用

(参考)

NEB法：https://doi.org/10.1142/9789812839664_0016
LUP法：https://doi.org/10.1063/1.460343

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

MultiOptPy-1.8.1cと同じディレクトリ内に、NEBという名前を含むディレクトリが生成されている。このディレクトリ内のplot_energy_9.pngを開き、エネルギーが極大値を示すノード番号を確認する。

私が実行した環境では、8番である。これらの２つの構造を確認する。 path_ITR_10_hydration_2ディレクトリ内のhydration_2_traj_8.xyzを可視化する。

hydration_2_traj_8.xyz

10
0 1
C      -0.507903462688     -0.173298880670      1.072003594524
O       0.125203423393      0.984638546479      0.991784288467
O      -0.766774658094     -0.699454522878     -0.488081413531
H       0.185276240295     -0.240224726482     -0.990809415614
H      -1.497781478560     -0.693637430024     -1.151542245045
H       0.057688828511     -1.053094613433      1.349791213539
H      -1.479014829751     -0.145113611972      1.526834643861
O       1.082757517610      0.563036362224     -1.076782433065
H       2.029887883307      0.506957830403     -1.224439528585
H       0.770660535977      0.950191046352     -0.008758704551

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

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

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

python optmain.py hydration_2_traj_8.xyz -opt rfo3_bofill -order 1 -fc 5 -bs 3-21g -func hf

-opt rfo3_bofillは遷移状態構造の最適化向けのoptimizerを指定する。
-order 1は一次の鞍点を求めることを指定する。（デフォルトだと極小値を求めるようになっている。）
-fc 5は5回の反復回数当たり1回正確なへシアンを計算することを指定する。

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

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

10
OptimizedStructure
C     -0.510624370151     -0.194103202553      1.054376574120
O      0.011615984299      1.002385072475      0.988427800453
O     -0.701263765283     -0.775604450822     -0.516779318578
H      0.199163840143     -0.273242580436     -1.016667712591
H     -1.537769530820     -0.501234286376     -0.918068395572
H      0.096709429933     -1.020256550766      1.394829388901
H     -1.528331275886     -0.279796341311      1.406367790981
O      1.112752004123      0.604620861423     -1.076917491175
H      2.046530015219      0.419489824122     -1.211032161109
H      0.811217668423      1.017741654243     -0.104536475430

参考#

https://github.com/psi4/psi4 (Psi4のgithubのレポジトリ)
https://github.com/ss0832/MultiOptPy (自作モジュールMultiOptPyのレポジトリ)
https://avogadro.cc/ (Avogadro、分子構造可視化ツール)
The Journal of Chemical Physics 2010, 132 (24), 241102.
The Journal of Chemical Physics 1991, 94 (1), 751–760.
In Classical and Quantum Dynamics in Condensed Phase Simulations; WORLD SCIENTIFIC: LERICI, Villa Marigola, 1998; pp 385–404.

概要#