最終更新：2025-07-19
概要
本記事では、自作モジュール(MultiOptPy)で、Vilsmeier-Haack反応のある１つの素過程の遷移状態構造を算出してみる。電子状態計算ソフトウェアにPySCFを使用する。 MultiOptPyは電子状態計算ソフトウェアを用いた分子構造最適化手法の勉強を目的として作成したpythonモジュールである。
MultiOptPyのレポジトリ:https://github.com/ss0832/MultiOptPy
Vilsmeier-Haack反応:
- https://www.chem-station.com/odos/2009/06/vilsmeier-haack-vilsmeier-haac.html
- https://en.wikipedia.org/wiki/Vilsmeier%E2%80%93Haack_reaction
使用した自作モジュール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. 初期構造の準備
モデル反応系として、以下の構造を用意した。今回は名前をvilsmeier_haack_rxn.xyzとした。 初期構造は以下のものを使用した。
24
1 0
C 0.86572620 -1.74333752 0.00002948
C 2.26712605 -1.74340509 0.00067179
C 2.96788488 -0.52979116 0.00014884
C 2.26724385 0.68389033 -0.00101454
C 0.86584400 0.68395789 -0.00165557
C 0.16508518 -0.52965603 -0.00113594
H 0.33068129 -2.66995868 0.00042970
H 2.80208102 -2.67007784 0.00156093
H 4.03788477 -0.52984275 0.00063752
H 2.80228877 1.61051149 -0.00141271
H 0.33088903 1.61063065 -0.00254450
H -0.90491471 -0.52960444 -0.00162599
C 0.72357821 -0.98050517 4.02566970
H 1.26227261 -1.90124608 4.10900341
N 1.36228736 0.12697164 3.82832049
C 0.62221163 1.39191433 3.71383078
H 1.23230768 2.19557710 4.06993324
H 0.36739756 1.56456421 2.68905675
H -0.27176493 1.33634228 4.29916629
C 2.82816895 0.12052531 3.71854974
H 3.14131532 0.91770003 3.07718260
H 3.25908558 0.25322489 4.68891137
H 3.15110820 -0.81404122 3.30965410
Cl -1.03149107 -0.97278741 4.15709372
2. 遷移状態構造最適化のための初期構造の算出
※あらかじめpythonで利用できる電子状態計算ソフトウェアのPySCF(導入方法はpipを使えば可能である。)が使える環境を用意しておく。
a. 初期構造をカレントディレクトリにvilsmeier_haack_rxn.xyzとして保存し、以下を実行する。
python optmain.py vilsmeier_haack_rxn.xyz -pyscf -spin 0 -elec 1 -func hf -bs 3-21g -opt rsirfo_fsb -ma 150 1 13 -ns 75-opt rsirfo_fsbは準ニュートン法であるRS-I-RFO法を構造最適化に使用することを示す。初期のへシアンに関しては、特にオプションで指定しない限り、単位行列が使われる。-pyscfは使用する電子状態計算ソフトウェアをPySCFにすることを示す。-spinはスピン多重度の指定である。PySCFを使用するときは目的とするスピン多重度に1を引いた値を指定する。-elec Nは形式電荷をNとすることを示す。-ma yyy a bはyyykJ/molの活性化障壁を超えうるペア同士を近づける力を原子のラベル番号aとbのペアに構造最適化時に加えることを示す。-func hfはHartree-Fock法を用いることを示す。-bs 3-21gは基底関数にPople系である3-21Gを用いることを示す。-ns 75は電子状態計算の反復回数を最大75回までとすることを示す。収束条件を満たすまで計算すると時間がかかるため、反復回数を制限した。
これを実行すると、計算レベルHF/3-21Gで指定した人工力ポテンシャルをかけた上で初期構造を構造最適化することができる。
結果はyyyy_mm_dd（今日の年月日）ディレクトリの中に存在するディレクトリを開いて確認できる。
正常終了していれば、このディレクトリ中に、vilsmeier_haack_rxn_traj.xyzが存在するので、これをコピーして、MultiOptPy-1.9.12cディレクトリに置く。
vilsmeier_haack_rxn_traj.xyzは構造最適化の過程をAvogadro（公式ページ:https://avogadro.cc/ ）等で可視化して確認できるようにしている。このvilsmeier_haack_rxn_traj.xyzは次のNEB計算に使用する。
※vilsmeier_haack_rxn_traj.xyzをアニメーションとして表示したい場合は、[https://github.com/ss0832/molecule_movie] を使うと良い。
（この人工力ポテンシャルを加えて行った構造最適化の結果はvilsmeier_haack_rxn_optimized.xyzで確認できる。構造を可視化して、生成系になっているか確認する。反応系のままであれば、-maの設定を見直してb.をやり直す。今回の場合以下の構造が得られた。生成系の安定構造に近いものが得られていると判断できるため、次のNEB計算を行うことが可能である。）
b. vilsmeier_haack_rxn_traj.xyzを初期パスとして、NEB法で経路の緩和を行う。
NEB法を用いることで、先ほど得られたvilsmeier_haack_rxn_traj.xyz全体のエネルギーを下げることができる。これにより、パスのエネルギー極大値を持つ構造を遷移状態構造に近づける。（この時点ではまだ正確な遷移状態構造は求められていない。）
python nebmain.py vilsmeier_haack_rxn_traj.xyz -ns 10 -fc 10 -func hf -bs 3-21g -pyscf -spin 0 -elec 1 -spng -nd 0.45-nd 0.45はノード間の距離を0.45Åとして初期パスを作成することを示す。-ns nはn回分NEB法による経路の緩和を行うことを示す。-fc 10は正確なへシアンを経路を緩和する反復計算10回ごとに計算することを示す。これを指定すると経路緩和のアルゴリズムがへシアンを使わない場合とは別のものになる。ノードごとにへシアンを計算するため、計算コストが急増する。-pyscfは使用する電子状態計算ソフトウェアをPySCFにすることを示す。-spngは緩和中のパスのエネルギープロファイルや各ノードの勾配のRMS値をmatplotlibで可視化するオプションである。-spinはスピン多重度の指定である。PySCFを使用するときは目的とするスピン多重度に1を引いた値を指定する。
c. 初期構造の決定及び遷移状態構造の計算
MultiOptPy-1.9.12cと同じディレクトリ内に、NEBという名前を含むディレクトリが生成されている。 そのディレクトリ内のenergy_plot.csvを確認し、緩和後のパスのエネルギー極大値を示す構造を確認する。
パスの緩和後の各ノードのエネルギー一覧(単位
※bias_force_rms.csvにて、各Iterationごとのすべてのノードの勾配のRMS値を確認できる。
私が実行した環境では、12番のノード(グラフでは13番)がエネルギー極大値を示していた。 path_ITR_10_vilsmeier_haack_rxnディレクトリ内のvilsmeier_haack_rxn_traj_12.xyzから遷移状態構造を求めようとしたがうまくいかなかった。そこで、エネルギー極大値を示す構造に近い勾配のRMS値が極小値を示す構造vilsmeier_haack_rxn_traj_9.xyzを使用した。しかしながら、こちらに関してもうまくいかなかった。
エネルギー極大値を示す構造でも、勾配のRMS値が極小値を示す構造でも目的とする遷移状態が求められなかったので、エネルギー極大値を示す構造の反応物側の隣のノードの構造(vilsmeier_haack_rxn_traj_11.xyz)を遷移状態構造の初期構造とした。（結合の組み換えに関わるC-C結合長が共有結合長よりも長い構造を使用した。）
※こちら[https://ss0832.github.io/molecule_viewer/] を使うことでも可視化は可能である。
vilsmeier_haack_rxn_traj_11.xyz
24
1 0
C -0.665535750741 -1.261969132048 -0.955253628666
C 0.723417379569 -1.383968592920 -1.204083178714
C 1.427876985302 -0.328995229630 -1.711857018642
C 0.748487323966 0.839239765699 -2.061681295642
C -0.629643362577 0.947878601723 -1.924805025325
C -1.340126436891 -0.101779064732 -1.395270596073
H -1.235377215659 -2.169306785266 -0.859458396852
H 1.219387801469 -2.304594345029 -0.959108336652
H 2.483514397399 -0.404954707846 -1.876930650880
H 1.298762675297 1.658840097994 -2.477831994844
H -1.133841614231 1.832146281414 -2.257984227883
H -2.405154725423 -0.041987364816 -1.287737228986
C -0.596607615947 -1.014689951705 1.012702168264
H -0.058559272505 -1.921129099822 1.173095231685
N -0.061150503367 0.106896965024 1.421853044968
C -0.860045755446 1.333832456218 1.634698740032
H -0.195216164773 2.108694798246 1.980347882501
H -1.325458033894 1.639707744169 0.708424710000
H -1.621161460265 1.147947350032 2.378004606053
C 1.387935404571 0.149489995731 1.757336503760
H 1.865965586407 0.934788904305 1.191049377619
H 1.495314270328 0.337182118564 2.817078798445
H 1.849631525087 -0.795617111757 1.513226132372
Cl -2.372415437677 -1.307653693548 1.384184383461
構造が壊れていないので、これを遷移状態を求めるための初期構造とする。
vilsmeier_haack_rxn_traj_17.xyzをMultiOptPy-1.9.12cと同じディレクトリ内にコピーする。
そして、以下を実行する。
python optmain.py vilsmeier_haack_rxn_traj_11.xyz -opt rsirfo_bofill -func hf -bs 3-21g -pyscf -spin 0 -elec 1 -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と同等）
実行して得られた正確な遷移状態構造を以下に示す。
(実行して得られた正確な遷移状態構造はvilsmeier_haack_rxn_traj_11_optimized.xyzとして保存されている。)
24
OptimizedStructure
C -0.661562604711 -1.265878831712 -0.902670691991
C 0.734086159028 -1.383939256473 -1.275476097479
C 1.405929407198 -0.329973813031 -1.798223385997
C 0.698427056896 0.847194886875 -2.091927851627
C -0.684412483938 0.947347275300 -1.929052630391
C -1.364696927261 -0.100494155794 -1.387040554704
H -1.222498641906 -2.186361681738 -0.919227488198
H 1.239236468544 -2.308003578725 -1.070076962051
H 2.451993277690 -0.388259575700 -2.014888415979
H 1.234898928541 1.685585304298 -2.492022995727
H -1.194473491808 1.835718859963 -2.237800033819
H -2.426748425375 -0.056911657151 -1.250267078804
C -0.538242611472 -1.011649243264 0.858622974356
H -0.005836105196 -1.892252479482 1.157403154962
N 0.027087901198 0.148347630522 1.271717066784
C -0.763764649201 1.370698400609 1.499693755132
H -0.124094843807 2.228975180508 1.355842787041
H -1.587799398781 1.426966181345 0.807781067137
H -1.157832862092 1.382985207949 2.508820858391
C 1.319320220225 0.093211707531 1.990860177855
H 1.886575693902 0.982260845117 1.756365653736
H 1.161642195941 0.044223837919 3.061977120644
H 1.886464606250 -0.771059296175 1.674582039266
Cl -2.313698869866 -1.298731748691 1.425007531463
40回程度の反復計算で遷移状態構造が得られた。-freqオプションにより生成されたnormal_modes.txtやvibration_animationディレクトリ内の振動モードのアニメーションを確認した。
以下に-freqオプションで生成されたnormal_modes.txtの一部を示す。
Mode 0 1 2
Freq [cm^-1] -260.9727 72.0127 88.3133
Reduced mass [au] 6.1150 3.7080 4.1690
Force const [Dyne/A] -0.2454 0.0113 0.0192
Char temp [K] 0.0000 103.6102 127.0631
Normal mode x y z x y z x y z
C -0.01226 -0.04670 -0.15571 0.04655 -0.00420 0.01738 0.01507 -0.05416 0.00113
C -0.02693 -0.02258 -0.03923 0.04961 0.05055 0.02152 0.02933 0.00492 0.02317
C -0.00631 0.01446 0.02305 -0.00362 0.05846 -0.03420 0.00117 0.05159 0.08564
C 0.00809 0.01084 0.02212 -0.05914 0.01088 -0.08988 -0.03958 0.03514 0.11250
C 0.01692 -0.00191 0.02278 -0.05787 -0.03801 -0.05996 -0.05362 -0.03594 0.04404
C 0.00196 -0.04484 -0.04115 -0.00407 -0.04470 -0.00123 -0.02567 -0.08013 -0.01247
H -0.01526 -0.03994 0.03412 0.08491 -0.02867 0.02816 0.05162 -0.07637 -0.00937
H -0.02713 -0.02838 -0.06596 0.09642 0.08362 0.05514 0.05895 0.01802 0.00879
H 0.00265 0.01691 0.06235 -0.00308 0.09494 -0.04153 0.01092 0.10659 0.11823
H 0.00933 0.01436 0.03380 -0.10330 0.01135 -0.14795 -0.05862 0.07865 0.17793
H 0.00957 0.00694 0.05798 -0.10298 -0.07294 -0.08621 -0.08059 -0.05058 0.04673
H -0.00148 -0.05242 -0.06921 -0.00279 -0.07978 0.01834 -0.03201 -0.13180 -0.04727
C -0.00602 0.03115 0.13137 0.00188 -0.00947 0.03765 -0.00429 -0.01131 -0.03216
H 0.01215 -0.00189 -0.01853 -0.02862 -0.02921 0.03322 -0.02580 -0.01044 0.01149
N -0.03826 -0.02020 0.10343 0.01503 -0.03642 0.07274 0.02785 -0.00193 -0.07187
C -0.00507 0.02075 0.01896 0.06122 0.00195 0.01888 0.04135 0.00183 -0.03762
H 0.01875 -0.00853 -0.03542 0.09887 -0.02752 0.00806 0.04679 -0.00264 -0.03970
H -0.00039 -0.01420 0.01257 0.07836 0.01413 -0.00046 0.02825 0.01930 -0.02006
H 0.00363 0.08275 0.02023 0.04152 0.04668 0.01057 0.06206 -0.00630 -0.02935
C 0.05396 0.03153 -0.04468 0.01038 -0.05984 0.07892 0.01797 -0.04110 -0.05658
H 0.02404 0.01994 -0.15070 0.08071 -0.12704 -0.00669 -0.00822 0.00188 0.04451
H 0.13910 0.09085 -0.03138 -0.00145 0.05920 0.08241 0.01017 -0.16054 -0.06282
H 0.01328 0.00892 -0.05021 -0.05556 -0.13728 0.17248 0.05163 0.01324 -0.14594
Cl 0.00154 0.00785 -0.01420 -0.02671 0.03193 -0.02900 -0.00819 0.05088 -0.01656
(...snip...)その結果、虚振動が１つであることが確認できた。
次に、vibration_animation内のmode_1_261i_wave_number.xyzをAvogadroで確認すると、想定される反応系と生成系をつなぐ方向に振動していることを確認できた。
終わりに
自作モジュールで、電子状態計算モジュールであるPySCFを用いて、計算レベルHF/3-21GでVilsmeier-Haack反応の１つの素過程の遷移状態構造を算出する手順を説明した。
参考
- 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.