2701 words
14 minutes
【計算化学】自作pythonモジュールで遷移状態構造を求めてみる（Diels-Alder反応）
最終更新：2025-06-29
概要
本記事では、自作モジュール(MultiOptPy)で、Diels-Alder反応の遷移状態構造を算出してみる。電子状態計算ソフトウェアにPySCFを使用する。 MultiOptPyは電子状態計算ソフトウェアを用いた分子構造最適化手法の勉強を目的として作成したpythonモジュールである。
MultiOptPyのレポジトリ:https://github.com/ss0832/MultiOptPy
Diels-Alder反応:
- https://www.chem-station.com/odos/2009/06/diels-alder.html
- https://en.wikipedia.org/wiki/Diels%E2%80%93Alder_reaction
使用した自作モジュールMultiOptPyのバージョン
1.9.12b環境
WSL2 (Ubuntu-22.04)Source codeのダウンロード
wget https://github.com/ss0832/MultiOptPy/archive/refs/tags/v1.9.12b.zip
unzip v1.9.12b.zip
cd MultiOptPy-1.9.12b手順
1. 初期構造の準備
モデル反応系として、プソラレンを含む以下の構造を用意した。今回は名前をdiels_alder_rxn.xyzとした。 初期構造は以下のものを使用した。
30
0 1
C 2.72104784 1.44878733 0.90937246
C 2.38025690 2.31720676 -0.14573910
C 1.81076227 0.43912806 1.32982341
C 1.10976392 2.19186944 -0.69255357
H 3.06933485 3.05309519 -0.50426359
C 4.06997061 1.60936951 1.61595997
C 0.52905337 0.35886019 0.73962342
O 2.16522375 -0.53085231 2.33535677
C 0.22379198 1.23346740 -0.30201506
C 0.42131350 3.01724437 -1.74326493
C 4.30748158 0.88190071 2.72081521
H 4.80070089 2.29919603 1.24843771
H -0.18516018 -0.36333060 1.07613415
C 3.26910513 -0.15624757 3.17764046
O -0.99179577 1.26786804 -1.09691379
C -0.82840962 2.50771178 -1.85916448
H 0.83010261 3.84829879 -2.27911940
H 5.20605766 1.03053198 3.28239313
O 3.39441433 -0.69806372 4.30649094
H -1.60653048 2.96014346 -2.43772606
C 5.31866628 -0.48991089 -2.02256963
H 5.87678882 -1.37594024 -2.24245878
H 4.24946995 -0.50807851 -2.05984083
C 5.96596287 0.65529381 -1.69686985
H 5.40784035 1.54132340 -1.47698162
O 7.39488892 0.67957487 -1.64706246
C 7.86129074 1.98950986 -1.98088828
H 7.46794690 2.69719140 -1.28134024
H 8.93048715 2.00767821 -1.94361981
H 7.53442414 2.24382111 -2.96749054
2. 遷移状態構造最適化のための初期構造の算出
※あらかじめpythonで利用できる電子状態計算ソフトウェアのPySCF(導入方法はpipを使えば可能である。)が使える環境を用意しておく。
a. 初期構造をカレントディレクトリにdiels_alder_rxn.xyzとして保存し、以下を実行する。
python optmain.py diels_alder_rxn.xyz -pyscf -spin 0 -elec 0 -func hf -bs 3-21g -opt rsirfo_fsb -ma 200 21 11 200 3 24 -modelhess-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を用いることを示す。-modelhessは準ニュートン法による構造最適化時に、入力した分子構造を基に、計算コストの低い近似へシアンを生成する。
これを実行すると、計算レベルHF/3-21Gで指定した人工力ポテンシャルをかけた上で初期構造を構造最適化することができる。
結果はyyyy_mm_dd（今日の年月日）ディレクトリの中に存在するディレクトリを開いて確認できる。
正常終了していれば、このディレクトリ中に、diels_alder_rxn_traj.xyzが存在するので、これをコピーして、MultiOptPy-1.9.12bディレクトリに置く。
diels_alder_rxn_traj.xyzは構造最適化の過程をAvogadro（公式ページ:https://avogadro.cc/ ）等で可視化して確認できるようにしている。このdiels_alder_rxn_traj.xyzは次のNEB計算に使用する。
※diels_alder_rxn_traj.xyzをアニメーションとして表示したい場合は、[https://github.com/ss0832/molecule_movie] を使うと良い。
（この人工力ポテンシャルを加えて行った構造最適化の結果はdiels_alder_rxn_optimized.xyzで確認できる。構造を可視化して、生成系になっているか確認する。反応系のままであれば、-maの設定を見直してb.をやり直す。今回の場合以下の構造が得られた。生成系の安定構造に近いものが得られていると判断できるため、次のNEB計算を行うことが可能である。）
b. diels_alder_rxn_traj.xyzを初期パスとして、NEB法で経路の緩和を行う。
NEB法を用いることで、先ほど得られたdiels_alder_rxn_traj.xyz全体のエネルギーを下げることができる。これにより、パスのエネルギー極大値を持つ構造を遷移状態構造に近づける。（この時点ではまだ正確な遷移状態構造は求められていない。）
python nebmain.py diels_alder_rxn_traj.xyz -ns 3 -modelhess -func hf -bs 3-21g -pyscf -spin 0 -elec 0 -spng -nd 0.5-nd 0.5はノード間の距離を0.5Åとして初期パスを作成することを示す。-ns 3は3回分NEB法による経路の緩和を行うことを示す。-modelhessは少ない計算コストで算出された近似へシアンを用いることを示す。これを指定すると経路緩和のアルゴリズムが-fcを1以上に指定した場合と同様のものになる。-pyscfは使用する電子状態計算ソフトウェアをPySCFにすることを示す。-spngは緩和中のパスのエネルギープロファイルや各ノードの勾配のRMS値をmatplotlibで可視化するオプションである。-spinはスピン多重度の指定である。PySCFを使用するときは目的とするスピン多重度に1を引いた値を指定する。
c. 初期構造の決定及び遷移状態構造の計算
MultiOptPy-1.9.12bと同じディレクトリ内に、NEBという名前を含むディレクトリが生成されている。 そのディレクトリ内のenergy_plot.csvを確認し、緩和後のパスのエネルギー極大値を示す構造を確認する。
パスの緩和後の各ノードのエネルギー一覧(単位
私が実行した環境では、30番のノード(グラフでは31番)がエネルギー極大値を示していた。 path_ITR_3_diels_alder_rxnディレクトリ内のdiels_alder_rxn_traj_30.xyzを可視化する。
※こちら[https://ss0832.github.io/molecule_viewer/] を使うことでも可視化は可能である。
diels_alder_rxn_traj_30.xyz
30
0 0
C -0.287825425533 0.497323306171 1.163923703707
C -0.827109173531 1.650641120130 0.552772727488
C -0.704178490232 -0.775414306544 0.721106805037
C -1.763507562282 1.508448001059 -0.416185138182
H -0.485824614680 2.613015825546 0.872874130386
C 0.786689382387 0.509550469847 2.068121801026
C -1.733499056133 -0.918701881440 -0.233380681136
O -0.504920631303 -1.885275557056 1.533549020840
C -2.195673113160 0.210048438767 -0.793794419642
C -2.548795530936 2.433677037783 -1.236133019190
C 1.344550444059 -0.694770479572 2.407989113337
H 1.277563735430 1.432207435297 2.315297529492
H -2.118765805081 -1.885877824297 -0.463976155024
C 0.515819576341 -1.900269407031 2.440821137728
O -3.181966989421 0.322493142774 -1.750544590880
C -3.357375830693 1.692582293645 -1.987481325197
H -2.483496492627 3.497065380960 -1.229162902521
H 2.236779518123 -0.749154426089 2.994853912955
O 0.696694140613 -2.866990762394 3.131079530198
H -4.099604477931 1.931362823257 -2.709570744060
C 2.055199135320 -1.422196287982 0.526456020322
H 2.866606575777 -0.745874771678 0.642722616776
H 2.210747532813 -2.428958950720 0.873744588637
C 1.101863611676 -1.177168914844 -0.410617907957
H 0.526228500339 -1.956535012659 -0.855383933890
O 1.137110236411 0.041951688604 -1.104592121592
C 2.191492915805 0.221556816452 -2.079116061206
H 2.272348686305 1.285692229614 -2.243343932583
H 3.142026989044 -0.166084289781 -1.726104633761
H 1.930822213101 -0.274343137819 -3.005925071108
構造が壊れていないので、これを遷移状態を求めるための初期構造とする。
diels_alder_rxn_traj_30.xyzをMultiOptPy-1.9.12bと同じディレクトリ内にコピーする。
そして、以下を実行する。
python optmain.py diels_alder_rxn_traj_30.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と同等）
実行して得られた正確な遷移状態構造を以下に示す。
(実行して得られた正確な遷移状態構造はdiels_alder_rxn_traj_30_optimized.xyzとして保存されている。)
30
OptimizedStructure
C -0.315801467495 0.478999674366 1.275338372244
C -0.892611378064 1.661454115712 0.720237587268
C -0.599230988464 -0.791730647390 0.678730444299
C -1.712714450187 1.546953823221 -0.339529095731
H -0.661680924248 2.610125731421 1.162323087725
C 0.640458315893 0.461131348240 2.264906817284
C -1.568023536287 -0.896627079063 -0.367755673251
O -0.487011415033 -1.932445001828 1.481727115593
C -2.027992827690 0.251009591450 -0.873170724280
C -2.471206199506 2.494747720001 -1.159726922138
C 1.274812860367 -0.749088559426 2.526290364506
H 1.022764035464 1.373783272351 2.676297976643
H -1.869441130312 -1.857437869461 -0.721715538994
C 0.490758726131 -1.989044966779 2.436162626818
O -2.909131892265 0.407843829547 -1.927573792523
C -3.143595591446 1.783947018769 -2.057055526904
H -2.479861913394 3.554335822396 -1.052330397703
H 2.081201041344 -0.796449201727 3.228786233164
O 0.662732723115 -2.994563710192 3.069879043527
H -3.811658654516 2.054991803328 -2.837771110429
C 2.173579135685 -1.175000505819 0.696701036706
H 2.900851354153 -0.395063649092 0.792618438038
H 2.503265827536 -2.159583185163 0.967754733598
C 1.212759675516 -1.045537126373 -0.302212298446
H 0.855128131396 -1.933496754579 -0.788520887490
O 1.212557495629 0.116571778105 -1.089356361744
C 1.862364317512 -0.007406580694 -2.371930220140
H 1.743116337906 0.944919869423 -2.863708181314
H 2.915867046340 -0.233260717330 -2.256392318485
H 1.397745344918 -0.784079843416 -2.969004827839
40回程度の反復計算で遷移状態構造が得られた。-freqオプションにより生成されたnormal_modes.txtやvibration_animationディレクトリ内の振動モードのアニメーションを確認した。
以下に-freqオプションで生成されたnormal_modes.txtの一部を示す。
Mode 0 1 2
Freq [cm^-1] -909.9287 42.9005 80.7197
Reduced mass [au] 8.7151 3.0396 4.6246
Force const [Dyne/A] -4.2515 0.0033 0.0178
Char temp [K] 0.0000 61.7242 116.1376
Normal mode x y z x y z x y z
C -0.02713 -0.01187 -0.02018 0.01081 0.00594 -0.02155 -0.00439 -0.00816 0.01534
C 0.01441 -0.01831 0.00279 0.00975 0.00573 -0.02103 0.01835 -0.00742 -0.00683
C -0.12507 0.02496 0.06868 -0.00582 0.00347 -0.00816 -0.02822 -0.01240 0.03465
C -0.00866 -0.00983 -0.01082 -0.01277 0.00296 -0.00335 0.02834 -0.01125 -0.01375
H 0.02442 -0.01579 -0.00628 0.02375 0.00747 -0.03206 0.04145 -0.00267 -0.02913
C 0.02254 -0.02154 0.02247 0.02013 0.00585 -0.03042 0.01878 -0.00738 -0.00735
C 0.00963 -0.00016 0.01765 -0.02351 0.00229 0.00775 -0.04220 -0.01598 0.04807
O 0.01149 0.00551 -0.01456 -0.00900 0.00925 0.00115 -0.01382 -0.00532 0.04009
C 0.00009 0.02001 -0.00245 -0.03053 0.00092 0.01175 -0.00548 -0.01586 0.01654
C 0.00009 -0.00226 0.00031 -0.02658 0.00120 0.00741 0.08074 -0.00890 -0.05955
C -0.06799 0.04181 0.10452 0.00787 0.00222 -0.02017 0.02286 -0.00867 -0.01932
H 0.04690 -0.00026 -0.04513 0.03211 0.00615 -0.04210 0.03745 -0.00651 -0.02657
H 0.00080 0.00201 0.01866 -0.03369 0.00159 0.01855 -0.06134 -0.01725 0.06799
C 0.00805 0.00865 -0.00258 -0.00353 0.00878 -0.00496 0.01899 -0.00737 0.00780
O 0.00482 0.00227 0.00322 -0.05387 -0.00172 0.03112 0.01804 -0.01647 -0.00305
C -0.00272 -0.00004 -0.00370 -0.05053 -0.00170 0.02760 0.07502 -0.01169 -0.05310
H 0.00114 -0.00227 0.00250 -0.01954 0.00192 0.00081 0.11995 -0.00503 -0.09467
H -0.00141 0.01879 0.02213 0.01258 0.00115 -0.02571 0.04331 -0.00743 -0.04270
O 0.00136 0.00327 -0.00568 -0.00825 0.01276 0.00291 0.03874 -0.00771 0.00162
H -0.00319 -0.00397 -0.00355 -0.06721 -0.00381 0.04113 0.10498 -0.01076 -0.07840
C 0.03778 -0.02628 -0.11893 -0.00967 -0.03020 -0.02596 -0.03153 -0.01304 -0.04475
H -0.03653 0.01402 0.06924 0.00086 -0.03926 -0.03020 -0.02524 -0.01777 -0.05390
H -0.05868 -0.00321 0.06356 -0.02120 -0.03272 -0.02109 -0.02706 -0.01796 -0.06774
C 0.13319 -0.02522 -0.04763 -0.01631 -0.02369 -0.01844 -0.06089 0.00097 -0.01522
H -0.10224 0.02096 0.05720 -0.02664 -0.02193 -0.01435 -0.06611 0.00949 -0.02694
O -0.00377 -0.00203 0.00073 -0.01103 -0.02656 -0.02457 -0.09844 0.02146 0.01475
C -0.00142 0.00475 -0.00334 0.18838 0.02263 0.07052 -0.02786 0.11067 0.04320
H -0.00312 0.00444 -0.00335 0.18951 0.01365 0.05276 -0.06982 0.12003 0.07107
H -0.00411 0.00219 0.00569 0.18851 0.10263 0.22763 -0.01778 0.18026 0.08689
H 0.00130 0.00251 -0.00278 0.33777 -0.01344 0.00150 0.05624 0.09943 -0.00739
(...snip...)その結果、虚振動が１つであることが確認できた。
次に、vibration_animation内のmode_1_910i_wave_number.xyzをAvogadroで確認すると、想定される反応系と生成系をつなぐ方向に振動していることを確認できた。
終わりに
自作モジュールで、電子状態計算モジュールであるPySCFを用いて、計算レベルHF/3-21GでDiels-Alder反応の遷移状態構造を算出する手順を説明した。
参考
- 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.
- https://en.wikipedia.org/wiki/Psoralen (ソラレン)
- https://ja.wikipedia.org/wiki/%E3%83%95%E3%83%A9%E3%83%8E%E3%82%AF%E3%83%9E%E3%83%AA%E3%83%B3 (フラノクマリン)
- https://ja.wikipedia.org/wiki/%E3%82%B8%E3%83%A3%E3%82%A4%E3%82%A2%E3%83%B3%E3%83%88%E3%83%BB%E3%83%9B%E3%82%B0%E3%82%A6%E3%82%A3%E3%83%BC%E3%83%89(バイカルハナウド)
【計算化学】自作pythonモジュールで遷移状態構造を求めてみる（Diels-Alder反応）
https://ss0832.github.io/posts/20250629_mop_usage_22/