最終更新：2025-05-24
概要
本記事では、自作モジュール(MultiOptPy)で、1,3,5-ヘキサトリエンの電子環状反応の遷移状態を算出してみる。電子状態計算ソフトウェアにPySCFを使用する。 MultiOptPyは電子状態計算ソフトウェアを用いた分子構造最適化手法の勉強を目的として作成したpythonモジュールである。
MultiOptPyのレポジトリ:https://github.com/ss0832/MultiOptPy
使用した自作モジュールMultiOptPyのバージョン
1.9.7c環境
WSL2 (Ubuntu-22.04)手順
1. 初期構造の準備
モデル反応系として、以下の構造を用意した。今回は名前をelectrocyclic_rxn_2.xyzとした。
14
0 1
C -2.51090453 1.81267715 2.56969022
H -1.86098900 1.15583690 3.10919990
H -3.56428389 1.79558817 2.75677281
C -1.99990041 2.66623817 1.64943188
H -0.94652083 2.68332851 1.46235066
C -2.93529270 3.61159802 0.87294232
H -3.98867207 3.59450896 1.06002483
C -2.42428828 4.46516082 -0.04731421
H -3.07420359 5.12200245 -0.58682248
C -0.90820933 4.48975769 -0.31657173
H -0.47152067 3.80503071 -1.01324379
C -0.08945414 5.24524305 0.45508160
H -0.50246928 5.88323965 1.20826075
H 0.97001027 5.20374093 0.31116239
2. 遷移状態構造最適化のための初期構造の算出
※あらかじめpythonで利用できる電子状態計算ソフトウェアのPySCF(導入方法はpipを使えば可能である。)が使える環境を用意しておく。
a. 以下のURLにアクセスし、Source codeをダウンロードする。
wget https://github.com/ss0832/MultiOptPy/archive/refs/tags/v1.9.7c.zip
unzip v1.9.7c.zip
cd MultiOptPy-1.9.7cb. 初期構造をカレントディレクトリにelectrocyclic_rxn_2.xyzとして保存し、以下を実行する。
python optmain.py electrocyclic_rxn_2.xyz -ma 200 1 12 -opt rsirfo_fsb -pyscf -spin 0-opt rsirfo_fsbは準ニュートン法であるRS-I-RFO法を構造最適化に使用することを示す。初期のへシアンに関しては、特にオプションで指定しない限り、単位行列が使われる。-pyscfは使用する電子状態計算ソフトウェアをPySCFにすることを示す。-spinはスピン多重度の指定である。PySCFを使用するときは目的とするスピン多重度に1を引いた値を指定する。-ma 6\200 1 12は200kJ/molの活性化障壁を超えうる力を原子のラベル番号1と7のペアに、構造最適化時に加えることを示す。
これを実行すると、-bs（基底関数の指定）と-func（計算手法の指定）がないので、デフォルトの設定である計算レベルB3LYP/6-31G(d)で指定した人工力ポテンシャルをかけた上で初期構造を構造最適化することができる。
結果はyyyy_mm_dd（今日の年月日）ディレクトリの中に存在するディレクトリを開いて確認できる。
正常終了していれば、このディレクトリ中に、electrocyclic_rxn_2_traj.xyzが存在するので、これをコピーして、MultiOptPy-1.9.7cディレクトリに置く。
electrocyclic_rxn_2_traj.xyzは構造最適化の過程をAvogadro（公式ページ:https://avogadro.cc/ ）等で可視化して確認できるようにしている。このelectrocyclic_rxn_2_traj.xyzは次のNEB計算に使用する。
※electrocyclic_rxn_2_traj.xyzをアニメーションとして表示したい場合は、[https://github.com/ss0832/molecule_movie] を使うと良い。
（この人工力ポテンシャルを加えて行った構造最適化の結果はelectrocyclic_rxn_2_optimized.xyzで確認できる。構造を可視化して、生成系になっているか確認する。反応系のままであれば、-maの設定を見直してb.をやり直す。）
c. electrocyclic_rxn_2_traj.xyzを初期パスとして、NEB法で経路の緩和を行う。
NEB法を用いることで、先ほど得られたelectrocyclic_rxn_2_traj.xyz全体のエネルギーを下げることができる。これにより、パスのエネルギー極大値を持つ構造を遷移状態構造に近づける。（この時点ではまだ正確な遷移状態構造は求められていない。）
LUP法と呼ばれるNEB法とよく似た方法があるが、本モジュールではNEB法を用いている。
python nebmain.py electrocyclic_rxn_2_traj.xyz -pyscf -spin 0 -ns 10 -spng -nd 0.5-nd 0.5はノード間の距離を0.5Åとして初期パスを作成することを示す。-ns 10は10回分NEB法による経路の緩和を行うことを示す。-pyscfは使用する電子状態計算ソフトウェアをPySCFにすることを示す。-spngは緩和中のパスのエネルギープロファイルや各ノードの勾配のRMS値をmatplotlibで可視化するオプションである。
(参考)
d. 初期構造の決定及び遷移状態構造の計算
MultiOptPy-1.9.7cと同じディレクトリ内に、NEBという名前を含むディレクトリが生成されている。 そのディレクトリ内のenergy_plot.csvを確認し、緩和後のパスのエネルギー極大値を示す構造を確認する。
パスの緩和後の各ノードのエネルギー一覧(単位
-233.3752611012588,-233.3734617038813,-233.36328233538123,-233.35796242878564,-233.35461311519904,-233.34842711976867,-233.3442292621993,-233.34664897118557,-233.34694618480293,-233.34865751578272,-233.3451143622679,-233.34299598479438,-233.34341705318926,-233.34568245475012,-233.34768356994982,-233.34889707016288,-233.34552776318847,-233.34150753791792,-233.33487737747234,-233.35163114775045,-233.366630820221,-233.38817625386773,-233.39918304653793,-233.40663796149838,-233.4115647233333,-233.41350504482764,-233.41410231149803-233.3754660277664,-233.37380849526085,-233.3633807465015,-233.35814418419753,-233.35482074180723,-233.3487472975622,-233.34458287809704,-233.34685802804333,-233.34722397035924,-233.3489695988609,-233.34546873517448,-233.34342454003553,-233.34378414609867,-233.3459835085816,-233.34812147771518,-233.3494379705279,-233.34562491438035,-233.34160400519963,-233.33506754146592,-233.3526670492169,-233.36761178642166,-233.3887815136188,-233.39958102566604,-233.40674443412908,-233.41165327001082,-233.41353937054484,-233.414133968936-233.37565702024102,-233.3740969985464,-233.36347183207442,-233.35831093967752,-233.35502227856597,-233.34905642658558,-233.34494705735653,-233.3470690851741,-233.34750034163687,-233.34931435582388,-233.34580171017492,-233.34385500604364,-233.34415731264824,-233.3462858371738,-233.34856259242363,-233.34999576863157,-233.34572222282216,-233.34170089938274,-233.33525237713334,-233.35370329007478,-233.36861695542297,-233.38937610776256,-233.3999321971194,-233.40685264624216,-233.4117344394628,-233.41356887611005,-233.41416271411595私が実行した環境では、18番のノードがエネルギー極大値を示していた。 path_ITR_10_electrocyclic_rxn_2ディレクトリ内のelectrocyclic_rxn_2_traj_18.xyzを可視化する。
※こちら[https://ss0832.github.io/molecule_viewer/] を使うことでも可視化は可能である。
electrocyclic_rxn_2_traj_18.xyz
14
0 0
C 0.180108310834 -0.860431729532 1.151941612224
H 0.760118246714 -1.515880709921 1.791284304904
H -0.228614573546 0.020148805179 1.602065783489
C -0.726230606420 -1.376438493423 0.214266918972
H -0.945955901197 -2.438929862966 0.122752088608
C -1.362830515699 -0.455198191341 -0.594919795098
H -2.251171223538 -0.729020458335 -1.151069626094
C -0.756096018100 0.764469772783 -0.961191757330
H -1.261809653207 1.313207328521 -1.748876887691
C 0.567577612749 1.120884719979 -0.725024134694
H 1.019479048062 1.579737670948 -1.613449437645
C 1.380316680151 0.717296975741 0.367107734631
H 1.164178312705 1.053613828940 1.370069393977
H 2.460930280494 0.806540343426 0.175043801746
構造が壊れていないので、これを遷移状態を求めるための初期構造とする。
electrocyclic_rxn_2_traj_18.xyzをMultiOptPy-1.9.7cと同じディレクトリ内にコピーする。
そして、以下を実行する。
python optmain.py electrocyclic_rxn_2_traj_18.xyz -opt rsirfo_bofill -fc 5 -order 1 -pyscf -spin 0 -freq-opt rsirfo_bofillは遷移状態構造の最適化向けのoptimizerを指定することを意味する。準ニュートン法であるRS-I-RFO法を使用する。今回は-fcで正確なへシアンを計算するようにしているので、初期へシアンは正確なへシアンを使用するようになっている。-order 1は一次の鞍点を求めることを指定する。（デフォルトだと極小値を求めるようになっている。）-fc 5は5回の反復回数当たり1回正確なへシアンを計算することを指定する。-freqは収束条件を満たした後に基準振動解析を行うことを示す。（自前で実装しているため、あくまで目安として使用することを推奨する。各振動モードをvibration_animation内のxyzファイルで可視化できる。）
実行して得られた正確な遷移状態構造を以下に示す。
(実行して得られた正確な遷移状態構造はelectrocyclic_rxn_2_traj_18_optimized.xyzとして保存されている。)
14
OptimizedStructure
C 0.229907745449 -0.968431525688 1.231812534917
H 0.860818399766 -1.784659671327 1.585337021061
H 0.024550860672 -0.234105291326 1.994053481174
C -0.770490637219 -1.264433024523 0.300296027688
H -0.979261650004 -2.314994571220 0.086200727884
C -1.392111491422 -0.337943968698 -0.543145616189
H -2.293915583187 -0.658079888874 -1.063770687221
C -0.746749956890 0.790318646715 -1.079717541001
H -1.214980582862 1.227512204616 -1.960858988545
C 0.590406651356 1.115787226110 -0.830516077579
H 1.169112404153 1.440918030107 -1.697960989544
C 1.275969710911 0.864712899782 0.362342105138
H 0.880634573507 1.270528599801 1.279674423098
H 2.366109555769 0.852870334523 0.336253579122
20回程度の反復計算で遷移状態構造が得られた。-freqオプションにより生成されたnormal_modes.txtやvibration_animationディレクトリ内の振動モードのアニメーションを確認した。
以下に-freqオプションで生成されたnormal_modes.txtの一部を示す。
Mode 0 1 2
Freq [cm^-1] -571.4990 101.9795 348.1487
Reduced mass [au] 3.3881 2.2559 2.1995
Force const [Dyne/A] -0.6520 0.0138 0.1571
Char temp [K] 0.0000 146.7258 500.9082
Normal mode x y z x y z x y z
C 0.09700 0.12759 -0.06512 0.08223 -0.08824 -0.03111 -0.01540 0.06199 -0.00969
H 0.13474 0.12985 -0.11393 0.06970 -0.16821 -0.18719 0.07514 0.15079 0.02837
H -0.09122 -0.07247 0.09031 0.22183 -0.16823 0.08181 -0.16487 0.04872 -0.03507
C -0.02417 0.02383 -0.00488 -0.02461 0.03947 0.03831 0.09480 -0.01507 -0.08038
H -0.11937 0.01036 0.16235 -0.12635 0.05694 0.04441 0.26306 -0.01674 -0.23059
C 0.00098 -0.02953 0.00347 -0.01113 0.06110 0.04624 -0.02458 0.04196 0.04964
H 0.03560 -0.04647 -0.04538 -0.05701 0.08508 0.11083 -0.06594 0.08015 0.09721
C 0.02354 0.00988 -0.01517 0.03536 -0.01889 -0.06609 -0.02793 0.03676 0.05204
H 0.03903 -0.04071 -0.04805 0.05679 -0.08526 -0.11032 -0.06380 0.08535 0.09487
C -0.03322 0.00787 0.00271 0.03005 -0.03030 -0.04267 0.03766 -0.11419 -0.03326
H -0.02006 0.18409 0.07963 0.09985 -0.10405 -0.02195 0.08392 -0.32926 -0.08159
C -0.06024 -0.14809 0.06564 -0.10215 0.05391 0.04725 -0.05832 -0.01318 0.02515
H 0.03204 0.14211 -0.01186 -0.28144 0.06395 -0.03300 -0.14365 0.08469 -0.05319
H -0.05707 -0.20610 0.04582 -0.09952 0.11673 0.21152 -0.05813 -0.08305 0.13826
(...snip...)その結果、虚振動が１つであることが確認できた。
次に、vibration_animation内のmode_1_571i_wave_number.xyzをAvogadroで確認すると、還元的脱離で想定されるような反応系と生成系をつなぐ方向に振動していることを確認できた。
終わりに
自作モジュールで、計算レベルB3LYP/6-31G(d)で1,3,5-ヘキサトリエンの電子環状反応の遷移状態構造を算出する手順を説明した。
参考
- 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.