最終更新：2025-07-19
概要
本記事では、自作モジュール(MultiOptPy)で、Criegee反応のある１つの素過程の遷移状態構造を算出してみる。電子状態計算ソフトウェアにPySCFを使用する。 MultiOptPyは電子状態計算ソフトウェアを用いた分子構造最適化手法の勉強を目的として作成したpythonモジュールである。
MultiOptPyのレポジトリ:https://github.com/ss0832/MultiOptPy
Criegee(酸化)反応:
- https://www.chem-station.com/odos/2009/08/criegee-criegee-glycol-oxidati.html
- https://en.wikipedia.org/wiki/Criegee_oxidation
使用した自作モジュール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. 初期構造の準備
モデル反応系として、以下の構造を用意した。今回は名前をcriegee_rxn.xyzとした。 初期構造は以下のものを使用した。
23
OptimizedStructure
C -1.193395288935 -2.899578087464 0.380451722914
C 0.185885236032 -3.123260274641 -0.267647880205
H -1.891013068565 -3.650983690794 0.010612662115
H -1.111867525854 -3.010927777500 1.466143413091
H 0.059642289344 -3.355398975267 -1.329601096961
H 0.675807423354 -3.971409045563 0.210916719775
O -1.710209161064 -1.629628671882 0.074854678854
O 1.011351476212 -1.996255685911 -0.126409109187
O -0.138225832914 0.860846595699 -1.785427719548
O 0.180261872548 0.997611117577 1.386578647899
C 1.296553791476 1.487968662769 1.939644175612
C -0.769806421477 2.000966021315 -2.102505195788
C 2.589887983325 1.075584586549 1.279844178044
H 2.693981959249 -0.011726313319 1.272084823603
H 2.613263238416 1.431090156217 0.246549497057
H 3.426118399674 1.509792496418 1.819435903911
C -1.929246608483 2.374352556342 -1.213045929382
H -2.669058651096 1.570570866888 -1.197103805155
H -1.583239522975 2.548793895340 -0.191988750945
H -2.396299523449 3.280573253272 -1.586139721160
O 1.234377011788 2.200851889853 2.872089107632
O -0.417019726045 2.634459509253 -3.027015659741
Pb -0.157749350561 -0.324293085152 -0.132320662438
2. 遷移状態構造最適化のための初期構造の算出
※あらかじめpythonで利用できる電子状態計算ソフトウェアのPySCF(導入方法はpipを使えば可能である。)が使える環境を用意しておく。
a. 初期構造をカレントディレクトリにcriegee_rxn.xyzとして保存し、以下を実行する。
python optmain.py criegee_rxn.xyz -opt rsirfo_fsb -func hf -bs def2svp -ecp default def2svp -pyscf -spin 0 -elec 0 -ma -300 7,8 23 -300 1 2 -ns 40-opt rsirfo_fsbは準ニュートン法であるRS-I-RFO法を構造最適化に使用することを示す。初期のへシアンに関しては、特にオプションで指定しない限り、単位行列が使われる。-pyscfは使用する電子状態計算ソフトウェアをPySCFにすることを示す。-spinはスピン多重度の指定である。PySCFを使用するときは目的とするスピン多重度に1を引いた値を指定する。-elec Nは形式電荷をNとすることを示す。-ma yyy a b -zzz c dはyyykJ/molの活性化障壁を超えうるペア同士を近づける力を原子のラベル番号aとbのペアに構造最適化時に加えることを示す。また、zzz kJ/molの活性化障壁を超えうるペア同士を遠ざける力を原子のラベル番号cとdに構造最適化時に加える。-func hfはHartree-Fock法を用いることを示す。-bs def2svpは基底関数にKarlsruhe系であるdef2-SVPを用いることを示す。-ecp default def2svpは有効内核ポテンシャル(ECP)を基底関数で、ECPが割り当てられているすべての原子に対して用いることを示す。PySCFを使用する場合は、必ず-ecpの指定をしないとECPが割り当てられないので、不可解な計算結果を得ることになる。-ns 40は電子状態計算の反復回数を最大40回までとすることを示す。収束条件を満たすまで計算すると時間がかかるため、反復回数を制限した。
これを実行すると、計算レベルHF/def2-SVPで指定した人工力ポテンシャルをかけた上で初期構造を構造最適化することができる。
結果はyyyy_mm_dd（今日の年月日）ディレクトリの中に存在するディレクトリを開いて確認できる。
正常終了していれば、このディレクトリ中に、criegee_rxn_traj.xyzが存在するので、これをコピーして、MultiOptPy-1.9.12cディレクトリに置く。
criegee_rxn_traj.xyzは構造最適化の過程をAvogadro（公式ページ:https://avogadro.cc/ ）等で可視化して確認できるようにしている。このcriegee_rxn_traj.xyzは次のNEB計算に使用する。
※criegee_rxn_traj.xyzをアニメーションとして表示したい場合は、[https://github.com/ss0832/molecule_movie] を使うと良い。
（この人工力ポテンシャルを加えて行った構造最適化の結果はcriegee_rxn_optimized.xyzで確認できる。構造を可視化して、生成系になっているか確認する。反応系のままであれば、-maの設定を見直してb.をやり直す。今回の場合以下の構造が得られた。生成系の安定構造に近いものが得られていると判断できるため、次のNEB計算を行うことが可能である。）
b. criegee_rxn_traj.xyzを初期パスとして、NEB法で経路の緩和を行う。
NEB法を用いることで、先ほど得られたcriegee_rxn_traj.xyz全体のエネルギーを下げることができる。これにより、パスのエネルギー極大値を持つ構造を遷移状態構造に近づける。（この時点ではまだ正確な遷移状態構造は求められていない。）
python nebmain.py criegee_rxn_traj.xyz -ns 10 -modelhess -func hf -bs def2svp -ecp default def2svp -pyscf -spin 0 -elec 0 -spng -nd 0.5 -nd Mはノード間の距離を M Åとして初期パスを作成することを示す。-ns nはn回分NEB法による経路の緩和を行うことを示す。-fc 10は正確なへシアンを経路を緩和する反復計算10回ごとに計算することを示す。これを指定すると経路緩和のアルゴリズムがへシアンを使わない場合とは別のものになる。ノードごとにへシアンを計算するため、計算コストが急増する。-pyscfは使用する電子状態計算ソフトウェアをPySCFにすることを示す。-spngは緩和中のパスのエネルギープロファイルや各ノードの勾配のRMS値をmatplotlibで可視化するオプションである。-bs def2svpは基底関数にKarlsruhe系であるdef2-SVPを用いることを示す。-ecp default def2svpは有効内核ポテンシャル(ECP)を基底関数で、ECPが割り当てられているすべての原子に対して用いることを示す。PySCFを使用する場合は、必ず-ecpの指定をしないとECPが割り当てられないので、不可解な計算結果を得ることになる。
c. 初期構造の決定及び遷移状態構造の計算
MultiOptPy-1.9.12cと同じディレクトリ内に、NEBという名前を含むディレクトリが生成されている。 そのディレクトリ内のenergy_plot.csvを確認し、緩和後のパスのエネルギー極大値を示す構造を確認する。
パスの緩和後の各ノードのエネルギー一覧(単位
※bias_force_rms.csvにて、各Iterationごとのすべてのノードの勾配のRMS値を確認できる。
私が実行した環境では、4番のノード(グラフでは5番)がエネルギー極大値を示していた。 path_ITR_10_criegee_rxnディレクトリ内のcriegee_rxn_traj_4.xyzから遷移状態構造を求めようとしたがうまくいかなかった。10回の反復計算後の4番のノードの勾配のRMS値が極小値を示していなかったため、目的とする正確な遷移状態から遠ざかった構造が得られたと仮定した。そこで、5回の反復計算計算後の4番のノードの勾配のRMS値を確認した。こちらに関しては極小値を示していたため、目的とする正確な遷移状態構造に近いものが得られていると考え、path_ITR_5_criegee_rxnディレクトリ内のcriegee_rxn_traj_4.xyzをcriegee_rxn_traj_4_itr5.xyzと改名して遷移状態構造を求めた。
※こちら[https://ss0832.github.io/molecule_viewer/] を使うことでも可視化は可能である。
criegee_rxn_traj_4_itr5.xyz
23
0 0
C -1.402126570970 -2.946845139516 0.389226730751
C 0.367813404551 -3.241873589925 -0.268366556946
H -1.982775542774 -3.522629401300 -0.334689097151
H -1.387800149249 -3.413243316101 1.384301170163
H 0.215618667218 -3.825657221916 -1.187662208916
H 0.780978968212 -3.842284008745 0.545758652629
O -1.735647397626 -1.665930160143 0.445499734043
O 1.036026658339 -2.119914167006 -0.497602133623
O -0.149057624349 0.993491926652 -1.764278487410
O 0.268380545556 1.133795753264 1.327595779701
C 1.368486851238 1.549076579955 1.950287648432
C -0.820046449814 2.091598482655 -2.120804637707
C 2.662669902014 1.113027186404 1.306304132637
H 2.726027952686 0.021259571077 1.290748341062
H 2.708153111791 1.475469885447 0.276640741420
H 3.507554116296 1.506291273821 1.862176606950
C -1.990849267375 2.450231247468 -1.239639356672
H -2.712898269992 1.629926292218 -1.232559973762
H -1.664554923911 2.637117165317 -0.214424666530
H -2.475068021712 3.340794500557 -1.626709930320
O 1.306620425225 2.206445531854 2.925265088399
O -0.504655872328 2.703444039302 -3.075418425981
Pb -0.122850513025 -0.273592431340 -0.141649151170
構造が壊れていないので、これを遷移状態を求めるための初期構造とする。
criegee_rxn_traj_17.xyzをMultiOptPy-1.9.12cと同じディレクトリ内にコピーする。
そして、以下を実行する。
python optmain.py criegee_rxn_traj_4_itr5.xyz -opt rsirfo_bofill -func hf -bs def2svp -ecp default def2svp -pyscf -spin 0 -elec 0 -order 1 -fc 5 -tcc -freq-opt rsirfo_bofillは遷移状態構造の最適化向けのoptimizerを指定することを意味する。準ニュートン法であるRS-I-RFO法を使用する。今回は-fcで正確なへシアンを計算するようにしているので、初期へシアンは正確なへシアンを使用するようになっている。-order 1は一次の鞍点を求めることを指定する。（デフォルトだと極小値を求めるようになっている。）-fc 5は5回の反復回数当たり1回正確なへシアンを計算することを指定する。-freqは収束条件を満たした後に基準振動解析を行うことを示す。（自前で実装しているため、あくまで目安として使用することを推奨する。各振動モードをvibration_animation内のxyzファイルで可視化できる。）-tccは収束条件を厳しくすることを示す。（Gaussianのtightと同等）-bs def2svpは基底関数にKarlsruhe系であるdef2-SVPを用いることを示す。-ecp default def2svpは有効内核ポテンシャル(ECP)を基底関数で、ECPが割り当てられているすべての原子に対して用いることを示す。PySCFを使用する場合は、必ず-ecpの指定をしないとECPが割り当てられないので、不可解な計算結果を得ることになる。
実行して得られた正確な遷移状態構造を以下に示す。
(実行して得られた正確な遷移状態構造はcriegee_rxn_traj_4_itr5_optimized.xyzとして保存されている。)
23
OptimizedStructure
C -1.427858648995 -3.093143521562 0.436858303300
C 0.176750033905 -3.372551442824 -0.625949155008
H -2.121227756357 -3.564362192683 -0.259223059959
H -1.146461218280 -3.697149839172 1.299693797976
H -0.300074910181 -3.778508088605 -1.518322622978
H 0.636888894578 -4.106658618633 0.034880560266
O -1.634333742377 -1.856130682478 0.671504041606
O 0.831957142357 -2.290532297258 -0.792771049638
O -1.316483793884 1.537640200497 0.067972895088
O 1.795392550607 0.937352713472 0.004411076217
C 1.906366834923 0.952259531410 1.252733312209
C -1.400292489155 1.663128949217 -1.176243402767
C 3.011847666280 1.748437865021 1.881352780105
H 3.942361075367 1.563695699610 1.348365571835
H 2.768798729985 2.806091559612 1.773086019630
H 3.109829116527 1.503036927613 2.934505677767
C -2.122434922654 2.848024880139 -1.746929052761
H -3.060488637715 2.993737543399 -1.215062007503
H -1.503862321516 3.730959723332 -1.581944741212
H -2.293366352426 2.716556505542 -2.811096073707
O 1.092928376495 0.324284809001 1.941654954116
O -0.872073247322 0.817098504011 -1.908271220860
Pb -0.074162380163 -0.383268728662 -0.011206603722230回程度の反復計算で遷移状態構造が得られた。-freqオプションにより生成されたnormal_modes.txtやvibration_animationディレクトリ内の振動モードのアニメーションを確認した。
以下に-freqオプションで生成されたnormal_modes.txtの一部を示す。
Mode 0 1 2
Freq [cm^-1] -994.1451 33.8752 37.2990
Reduced mass [au] 7.4947 4.0354 1.2596
Force const [Dyne/A] -4.3642 0.0027 0.0010
Char temp [K] 0.0000 48.7389 53.6649
Normal mode x y z x y z x y z
C -0.13748 0.03398 0.09087 0.07511 -0.02087 0.09352 0.01058 -0.01056 0.01787
C 0.14122 -0.01465 -0.09030 0.07645 -0.01394 0.09373 -0.01362 -0.00475 -0.01870
H 0.09087 -0.01646 -0.09493 0.08694 -0.08425 0.12476 -0.00227 -0.01515 0.03386
H 0.10876 -0.01402 -0.02913 0.10698 0.03523 0.12252 0.03200 -0.00751 0.01314
H -0.10646 0.02591 0.02948 0.08480 -0.07942 0.11916 -0.03240 0.00602 -0.01367
H -0.09149 0.01328 0.09483 0.11071 0.03859 0.12835 -0.00330 -0.01305 -0.03520
O -0.03764 -0.07316 0.01252 0.02884 -0.01658 0.03046 0.00996 -0.01093 0.01989
O 0.00780 -0.08108 -0.01703 0.03260 0.00289 0.03102 -0.01309 -0.00521 -0.02048
O -0.01165 0.01765 -0.00229 -0.00995 0.00810 0.00901 0.01094 0.00461 0.01013
O 0.01736 0.01188 0.00315 -0.01236 -0.00430 0.00863 -0.00839 0.00896 -0.00970
C 0.00189 0.00230 0.00016 -0.01073 -0.03391 0.00942 -0.01196 0.01629 -0.01078
C -0.00090 0.00285 -0.00001 0.00259 0.03488 0.01144 0.01695 0.00990 0.01158
C 0.00241 0.00165 0.00147 0.00410 -0.06494 0.02288 -0.04954 0.05993 -0.00027
H 0.00431 0.00473 0.00213 0.01588 -0.14378 0.06967 -0.09457 0.36041 -0.17904
H 0.00430 0.00262 0.00244 0.07251 -0.05592 -0.03819 -0.27133 0.04057 0.29396
H 0.00440 0.00389 0.00286 -0.05153 -0.02581 0.03703 0.14254 -0.15838 -0.06839
C -0.00160 0.00251 -0.00135 0.02780 0.05735 0.02646 0.06841 0.03693 0.00293
H -0.00221 0.00611 -0.00181 0.06773 0.12355 0.07764 0.22173 0.29071 0.19977
H -0.00299 0.00417 -0.00224 0.08873 0.02688 -0.03600 0.27197 -0.04777 -0.29737
H -0.00260 0.00540 -0.00259 -0.03877 0.04095 0.03903 -0.19448 -0.09901 0.06134
O 0.01158 0.00857 0.01562 -0.01540 -0.04153 -0.00436 -0.00275 0.00348 -0.00961
O -0.00738 0.01314 -0.01499 0.00123 0.04445 -0.00183 0.00367 0.00147 0.00986
Pb 0.00118 0.00610 0.00018 -0.01467 0.00353 -0.02358 -0.00157 -0.00813 -0.00020
(...snip...)その結果、虚振動が１つであることが確認できた。
次に、vibration_animation内のmode_1_994i_wave_number.xyzをAvogadroで確認すると、想定される反応系と生成系をつなぐ方向に振動していることを確認できた。
終わりに
自作モジュールで、電子状態計算モジュールであるPySCFを用いて、計算レベルHF/def2-SVPでCriegee反応の１つの素過程の遷移状態構造を算出する手順を説明した。
参考
- 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.