Home
2796 words
14 minutes
【計算化学】自作pythonモジュールで遷移状態構造を求めてみる（分子内Diels-Alder反応）
2025-06-29
Computational Chemistry
MultiOptPy
/
python

最終更新：2025-06-29

概要#

本記事では、自作モジュール(MultiOptPy)で、分子内Diels-Alder反応の遷移状態構造を算出してみる。電子状態計算ソフトウェアにPySCFを使用する。 MultiOptPyは電子状態計算ソフトウェアを用いた分子構造最適化手法の勉強を目的として作成したpythonモジュールである。

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

Diels-Alder反応:

使用した自作モジュール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. 初期構造の準備#

モデル反応系として、Imperatorinを含む以下の構造を用意した。今回は名前をintramolecular_diels_alder_rxn.xyzとした。 初期構造は以下のものを使用した。

34
0 1
 C                  2.71530609    1.49977504    0.71124648
 C                  2.10392950    2.39648435   -0.18635200
 C                  1.97916266    0.40491854    1.24434421
 C                  0.75541804    2.20803809   -0.45998169
 H                  2.65423782    3.19914006   -0.63110426
 C                  4.17251123    1.72048420    1.12656020
 C                  0.60769952    0.26281553    0.93377730
 O                  2.59928775   -0.58922703    2.08382160
 C                  0.03345773    1.16761659    0.04231018
 C                 -0.19131080    3.03619743   -1.28281723
 C                  4.68544074    0.95391000    2.10410540
 H                  4.76425102    2.48080894    0.66110310
 C                  3.83397152   -0.18064707    2.69777814
 O                 -1.32452690    1.15566983   -0.47386485
 C                 -1.40377641    2.44344849   -1.16688654
 H                  0.04191367    3.92265645   -1.83481066
 H                  5.67319992    1.13752306    2.47223317
 O                  4.23153520   -0.77184672    3.73508073
 H                 -2.31390633    2.86818234   -1.53589155
 O                 -0.16798268   -0.79241567    1.50797609
 C                 -0.10172129   -1.94571602    0.66511921
 H                 -0.68212685   -2.73529461    1.09476445
 H                 -0.48985126   -1.70338258   -0.30210860
 C                  1.36360416   -2.40010794    0.53114414
 H                  2.15417469   -1.68157446    0.59120987
 C                  1.65179974   -3.71002568    0.33717031
 C                  3.11712530   -4.16441791    0.20319752
 H                  3.21546270   -5.16447256    0.57078485
 H                  3.40774087   -4.13256007   -0.82608751
 H                  3.74628813   -3.51193517    0.77181022
 C                  0.51396915   -4.74417635    0.25071830
 H                  0.29920198   -5.12078330    1.22895256
 H                 -0.36114280   -4.27876159   -0.15235511
 H                  0.81327767   -5.55151742   -0.38450986
初期パスを求めるための初期構造

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

※あらかじめpythonで利用できる電子状態計算ソフトウェアのPySCF(導入方法はpipを使えば可能である。)が使える環境を用意しておく。

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

python optmain.py intramolecular_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（今日の年月日）ディレクトリの中に存在するディレクトリを開いて確認できる。

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

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

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

（この人工力ポテンシャルを加えて行った構造最適化の結果はintramolecular_diels_alder_rxn_optimized.xyzで確認できる。構造を可視化して、生成系になっているか確認する。反応系のままであれば、-maの設定を見直してb.をやり直す。今回の場合以下の構造が得られた。生成系の安定構造に近いものが得られていると判断できるため、次のNEB計算を行うことが可能である。）

初期パスを求めるための構造最適化の結果（安定構造ではない）

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

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

python nebmain.py intramolecular_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を確認し、緩和後のパスのエネルギー極大値を示す構造を確認する。

パスの緩和後の各ノードのエネルギー一覧(単位) (energy_plot.csvに保存されている。)

NEB計算の結果の可視化

私が実行した環境では、29番のノード(グラフでは30番)がエネルギー極大値を示していた。 path_ITR_3_intramolecular_diels_alder_rxnディレクトリ内のintramolecular_diels_alder_rxn_traj_29.xyzを可視化する。

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

intramolecular_diels_alder_rxn_traj_29.xyz

34
0 0
C       0.839305925060      1.313619671902      0.357891080948
C       0.589571906090      2.419557215578     -0.505677237850
C      -0.290045132477      0.590272943085      0.875578091260
C      -0.682248945879      2.776088601086     -0.753074732134
H       1.406191939877      2.929613251024     -0.978739439158
C       2.019887696091      0.646740216760      0.517542597434
C      -1.639266279232      0.980804025468      0.572051851952
O      -0.119113878177     -0.042110336748      2.112073901441
C      -1.787472433271      2.024595842220     -0.235004842825
C      -1.287886449189      3.756048027656     -1.656105605096
C       1.993090658534     -0.562195633712      1.218439816179
H       2.915789262494      0.936272259539      0.001238274842
C       1.074156005478     -0.686320954103      2.336330705225
O      -2.954684699327      2.510719128397     -0.774298843176
C      -2.606893708916      3.548919096267     -1.642484620909
H      -0.725084046353      4.436253553231     -2.252167361609
H       2.861171905864     -1.175901776425      1.280056118754
O       1.221484862742     -1.327801472726      3.343734861180
H      -3.451607768331      3.958135450610     -2.145650975811
O      -2.518176410139     -0.075001379739      0.671000097751
C      -1.842023415700     -1.258985221445      0.067000731933
H      -1.875648772552     -2.076032169645      0.763413911504
H      -2.442910374257     -1.492053968254     -0.799237249731
C      -0.402877038768     -0.885560367177     -0.307361347367
H      -0.321779776758     -0.328806827736     -1.232254175334
C       0.684381812753     -1.754450594658      0.005657732552
C       1.727374278709     -2.019263917741     -1.061609743275
H       2.674751300829     -2.238149005611     -0.590058590973
H       1.419873948724     -2.864464927787     -1.672361462854
H       1.843735456315     -1.172716704281     -1.716634689038
C       0.460277192245     -2.971800215761      0.896497559708
H       0.282071276647     -2.734171640579      1.933885979922
H      -0.388786837151     -3.539038535740      0.537587838433
H       1.323390538022     -3.622813632954      0.832739766123
NEB法により緩和したパスのエネルギー極大値を示す構造

構造が壊れていないので、これを遷移状態を求めるための初期構造とする。

intramolecular_diels_alder_rxn_traj_29.xyzをMultiOptPy-1.9.12bと同じディレクトリ内にコピーする。

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

python optmain.py intramolecular_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と同等）

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

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

34
OptimizedStructure
C      0.819008649999      1.298693873919      0.282380662743
C      0.621019265710      2.448770278770     -0.542995983153
C     -0.320298387597      0.603458983032      0.774971366491
C     -0.644808050833      2.828254994973     -0.816830736788
H      1.468992203569      2.978903795526     -0.927705293133
C      2.007014727486      0.647715710746      0.548136156774
C     -1.654443055021      0.971938156784      0.401886065450
O     -0.206850765499     -0.008786988204      2.029218013595
C     -1.776715471256      2.041481350967     -0.386166671976
C     -1.222881357770      3.883522917776     -1.652694515923
C      1.952294095451     -0.547835582255      1.256110891200
H      2.924957542983      0.972930956813      0.099917295076
C      0.959963202306     -0.674171630389      2.316342045009
O     -2.924113842731      2.555181164800     -0.963786054831
C     -2.534234019289      3.678678380989     -1.698072876654
H     -0.684626437521      4.668085604260     -2.131113645685
H      2.838521222850     -1.129033594614      1.402093871576
O      1.058618902332     -1.298653573281      3.338775591208
H     -3.326724174038      4.198670696365     -2.178639023020
O     -2.554090080064     -0.047855454591      0.567920060516
C     -1.859416872845     -1.348208416976      0.341853349328
H     -1.794644985091     -1.841506287522      1.292073274844
H     -2.487901347796     -1.903965523175     -0.335987089486
C     -0.481078148765     -1.070308237503     -0.261307654897
H     -0.533234322837     -0.633541121118     -1.242806883929
C      0.659958328477     -1.853726593633     -0.013513406108
C      1.644224067235     -1.977552301675     -1.169958004079
H      2.600734690394     -2.359253459454     -0.837358164354
H      1.243450541878     -2.670307891695     -1.906419355319
H      1.796295231393     -1.021723444122     -1.653849384454
C      0.607346346055     -3.100266325950      0.864384727246
H      0.003318605670     -2.985943942929      1.750191985569
H      0.198585432113     -3.916130313366      0.272699545405
H      1.601758263051     -3.387516183270      1.180249841762
遷移状態構造

40回程度の反復計算で遷移状態構造が得られた。-freqオプションにより生成されたnormal_modes.txtやvibration_animationディレクトリ内の振動モードのアニメーションを確認した。

以下に-freqオプションで生成されたnormal_modes.txtの一部を示す。

Mode                                 0                   1                   2
Freq [cm^-1]                     -801.1929             68.7542            124.2088
Reduced mass [au]                   9.4124              4.8673              3.9188
Force const [Dyne/A]               -3.5598              0.0136              0.0356
Char temp [K]                       0.0000             98.9219            178.7088
Normal mode                   x         y         z            x         y         z            x         y         z
       C                 0.03149   -0.02217   -0.01708   -0.00402   -0.01197   -0.03948    0.00386   -0.02914   -0.04132
       C                -0.00205    0.01900   -0.00322   -0.00714    0.00693   -0.01261    0.00933   -0.04133   -0.05926
       C                 0.00107   -0.12819   -0.08008   -0.00071   -0.02169   -0.04607   -0.00162   -0.01004   -0.03392
       C                 0.01477    0.00342    0.00324   -0.00850    0.01941    0.01083    0.00980   -0.02269   -0.03653
       H                -0.00201    0.02966    0.01270   -0.00861    0.01906    0.00098    0.01021   -0.05173   -0.07152
       C                -0.01657    0.02390   -0.01883    0.00020   -0.00209   -0.03097   -0.00051   -0.02064    0.00198
       C                -0.02264    0.00574   -0.00772   -0.00064   -0.02372   -0.04989   -0.00158   -0.01388   -0.03362
       O                -0.00043    0.01327    0.01742   -0.00574   -0.00621   -0.03633   -0.01888    0.03074   -0.01269
       C                -0.00931   -0.00853    0.01546   -0.00650    0.00453   -0.01073   -0.00119   -0.00157   -0.01885
       C                 0.00118    0.00201    0.00084   -0.01500    0.06085    0.06783    0.01488    0.00963    0.00175
       C                -0.06678   -0.08699   -0.06166    0.00947    0.01246   -0.00774   -0.01192   -0.00031    0.03434
       H                 0.01739    0.05260    0.07329   -0.00021    0.00183   -0.02916    0.00831   -0.03304    0.01109
       C                -0.00456    0.00371    0.00947    0.00562    0.02562   -0.01025   -0.03836    0.01177    0.01162
       O                -0.00705    0.00102   -0.00198   -0.01179    0.03492    0.02624   -0.00405    0.04812    0.02960
       C                 0.00396   -0.00201    0.00147   -0.01675    0.06974    0.07689    0.00646    0.05504    0.04547
       H                -0.00106    0.00289    0.00013   -0.01830    0.08256    0.09971    0.02347    0.00475    0.00329
       H                -0.03143   -0.02567   -0.00903    0.01395    0.02315    0.00707   -0.01481    0.00257    0.06330
       O                 0.00105   -0.00073    0.00572    0.01313    0.06020    0.01037   -0.07696    0.00685    0.01310
       H                 0.00734    0.00020   -0.00349   -0.02128    0.09796    0.11482    0.00594    0.09503    0.08961
       O                -0.01226   -0.00477   -0.00208    0.00122   -0.02597   -0.06636    0.00258   -0.01332   -0.00458
       C                -0.00188    0.00144    0.00058    0.00596   -0.02670   -0.03645    0.01201   -0.02777    0.09561
       H                -0.00195    0.00170    0.00031   -0.00895   -0.01264   -0.02847    0.05078    0.05519    0.13677
       H                 0.02203   -0.00656   -0.01279    0.01668   -0.03692   -0.03813   -0.00389   -0.09254    0.16400
       C                 0.01149    0.12182    0.06222    0.01695   -0.03959   -0.01847   -0.01346   -0.04705    0.03086
       H                -0.01783   -0.04418   -0.01459    0.03672   -0.06387   -0.03034   -0.06259   -0.07361    0.02166
       C                 0.09395    0.05780    0.07586    0.00690   -0.04008    0.02479    0.01723   -0.01091   -0.00044
       C                -0.00184   -0.00182   -0.00094    0.02700   -0.09709    0.04845   -0.00050    0.07484   -0.02464
       H                -0.00009   -0.02851   -0.03690    0.01808   -0.08943    0.08294    0.02490    0.10967   -0.05794
       H                -0.05446    0.01024    0.01543    0.03592   -0.12764    0.07212    0.01726    0.08369   -0.04244
       H                 0.00768   -0.00472   -0.00361    0.04308   -0.12000    0.00845   -0.05708    0.10048    0.00801
       C                 0.00006    0.00069   -0.00372   -0.01310   -0.00286    0.07723    0.09591   -0.03061   -0.02604
       H                 0.00806    0.00196    0.00197   -0.03749    0.03161    0.05596    0.12855   -0.06704    0.00149
       H                -0.03155    0.03510   -0.02898    0.00514   -0.02878    0.10019    0.09903   -0.03244   -0.02551
       H                -0.01052   -0.04334   -0.01021   -0.02110    0.01314    0.11654    0.11937   -0.00354   -0.07472

(...snip...)

その結果、虚振動が１つであることが確認できた。

次に、vibration_animation内のmode_1_801i_wave_number.xyzをAvogadroで確認すると、想定される反応系と生成系をつなぐ方向に振動していることを確認できた。

終わりに#

　 　自作モジュールで、電子状態計算モジュールであるPySCFを用いて、計算レベルHF/3-21Gで分子内Diels-Alder反応の遷移状態構造を算出する手順を説明した。

参考#

【計算化学】自作pythonモジュールで遷移状態構造を求めてみる（分子内Diels-Alder反応）
https://ss0832.github.io/posts/20250629_mop_usage_23/
Author
ss0832
Published at
2025-06-29
【計算化学】自作pythonモジュールで遷移状態構造を求めてみる（Diels-Alder反応）
【計算化学】自作pythonモジュールで遷移状態構造を求めてみる（分子内Mannich反応）