【計算化学】自作pythonモジュールで遷移状態構造を求めてみる(ケテンの[4+2]付加環化反応）

最終更新：2025-07-28

概要#

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

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

ケテンの[4+2]付加環化反応について：

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

モデル反応系として、以下の構造を用意した。今回はファイルの名前をketene_42_cycloaddition.xyzとした。初期構造は以下のものを使用した。

15

C     -1.103526951971     -0.581207884523     -1.730463328277
C      0.010283551945     -0.205084029153     -2.265536480479
H      0.298980006538      0.823020384479     -2.236674291882
H      0.674473246357     -0.935957267595     -2.676606622964
O     -2.105189191051     -0.919846293171     -1.240822899703
C      0.506169644883     -1.376039720567      1.074997884669
H      1.070997008867     -2.285737915590      1.146217526445
H     -0.552128861637     -1.486159712423      0.939290081985
C      1.100782427869     -0.199790746468      1.150003486619
H      2.167220882144     -0.166071381726      1.280824476959
C      0.446043198443      1.120599531691      1.048466166971
H      1.106512360336      1.962802686665      1.150828500573
C     -0.838427848492      1.345273105712      0.838044792993
H     -1.223448266178      2.345682189580      0.784364908581
H     -1.558741208054      0.558517053088      0.737065797510

初期パスを求めるための初期構造

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

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

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

python optmain.py ketene_42_cycloaddition.xyz -pyscf  -spin 0 -elec 0 -func hf -bs 3-21g -ma 200 2 13 200 1 6 -opt rsirfo_fsb

-opt rsirfo_fsbは準ニュートン法であるRS-I-RFO法を構造最適化に使用することを示す。初期のへシアンに関しては、特にオプションで指定しない限り、単位行列が使われる。
-pyscfは使用する電子状態計算ソフトウェアをPySCFにすることを示す。
-spinはスピン多重度の指定である。PySCFを使用するときは目的とするスピン多重度に1を引いた値を指定する。
-elec Mは形式電荷をMとすることを示す。
-ma yyy a bはyyykJ/molの活性化障壁を超えうるペア同士を近づける力を原子のラベル番号aとbのペアに構造最適化時に加えることを示す。
-func hfはHartree-Fock法を用いることを示す。
-bs 3-21gは基底関数にPople系である3-21Gを用いることを示す。
-modelhessは-optでニュートン法ファミリーのoptimzierを使用した際に、電子状態計算でへシアンを算出する場合よりも大幅に低い計算コストで、精度の粗い初期へシアンを計算することを指定するオプションである。

これを実行すると、計算レベルHF/3-21Gで指定した人工力ポテンシャルをかけた上で初期構造を構造最適化することができる。

結果はyyyy_mm_dd（今日の年月日）ディレクトリの中に存在するディレクトリを開いて確認できる。

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

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

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

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

15
OptimizedStructure
C     -0.731299116983     -1.014244837014     -0.565422444696
C     -0.453776508919      0.333969125679     -1.210982822213
H     -1.249974454591      0.535667300098     -1.914660608627
H      0.478841544372      0.247080049209     -1.764470635218
O     -1.697573132565     -1.687306291557     -0.860580358028
C      0.263607298029     -1.444340577474      0.395752228376
H      1.052118051702     -1.975682987322     -0.139531132121
H     -0.182415914180     -2.162120492911      1.074713475862
C      0.883050304453     -0.297337541293      1.170340231542
H      1.556970976543     -0.568998367251      1.960058355330
C      0.613684198404      0.967301400310      0.922558734899
H      1.056230540842      1.742769973518      1.518359032475
C     -0.327725610219      1.390124201958     -0.190026171538
H      0.036948475544      2.308960269341     -0.637259468642
H     -1.298686652433      1.624158774710      0.241151582598

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

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

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

python nebmain.py ketene_42_cycloaddition_traj.xyz -ns 10 -modelhess -nd 0.3 -func hf -bs 3-21g -pyscf -elec 0 -spin 0 -spng

-nd Nはノード間の距離をN Åとして初期パスを作成することを示す。
-ns nはn回分NEB法による経路の緩和を行うことを示す。
-modelhessは少ない計算コストで算出された近似へシアンを用いることを示す。これを指定すると経路緩和のアルゴリズムが-fcを1以上に指定した場合と同様のものになる。
-pyscfは使用する電子状態計算ソフトウェアをPySCFにすることを示す。
-spngは緩和中のパスのエネルギープロファイルや各ノードの勾配のRMS値をmatplotlibで可視化するオプションである。
-spinはスピン多重度の指定である。PySCFを使用するときは目的とするスピン多重度に1を引いた値を指定する。

c. 初期構造の決定及び遷移状態構造の計算#

MultiOptPy-1.9.12cと同じディレクトリ内に、NEBという名前を含むディレクトリが生成されている。そのディレクトリ内のenergy_plot.csvを確認し、緩和後のパスのエネルギー極大値を示す構造を確認する。

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

NEB計算の結果の可視化

※bias_force_rms.csvにて、各Iterationごとのすべてのノードの勾配のRMS値を確認できる。

私が実行した環境では、21番のノード(グラフでは22番)がエネルギー極大値を示していた。

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

ketene_42_cycloaddition_traj_21.xyz

15
0 0
C      -0.564537716268     -1.107516665902     -1.118106519244
C      -0.362954594268      0.227882468767     -1.448568409746
H      -1.088623018125      0.719031364136     -2.054315647367
H       0.659960704142      0.504367234090     -1.566644083900
O      -1.005469096874     -2.190424861021     -1.341487769867
C       0.072689913954     -1.285938187942      0.756538927320
H       0.418758547026     -2.296102177155      0.760269935267
H      -0.956708771579     -1.152663299442      1.029487556434
C       1.009057767995     -0.297974248249      1.066413692695
H       2.022907971444     -0.607670782306      1.211177902179
C       0.732122606203      1.029216388086      0.905779244310
H       1.518569012649      1.755402681857      0.936585825124
C      -0.479440183758      1.398755075756      0.314354971156
H      -0.593865519036      2.426350640282      0.016620064644
H      -1.382467623504      0.877284369041      0.531894310994

NEB法により緩和したパスのエネルギー極大値を示す構造

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

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

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

python optmain.py ketene_42_cycloaddition_traj_21.xyz -opt rsirfo_bofill -fc 5 -order 1 -func hf -bs 3-21g -pyscf -elec -spin 0 -tcc -freq

-opt rsirfo_bofillは遷移状態構造の最適化向けのoptimizerを指定することを意味する。準ニュートン法であるRS-I-RFO法を使用する。今回は-fcで正確なへシアンを計算するようにしているので、初期へシアンは正確なへシアンを使用するようになっている。
-order 1は一次の鞍点を求めることを指定する。（デフォルトだと極小値を求めるようになっている。）
-fc 5は5回の反復回数当たり1回正確なへシアンを計算することを指定する。
-freqは収束条件を満たした後に基準振動解析を行うことを示す。（自前で実装しているため、あくまで目安として使用することを推奨する。各振動モードをvibration_animation内のxyzファイルで可視化できる。）
-tccは収束条件を厳しくすることを示す。（Gaussianのtightと同等）

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

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

15
OptimizedStructure
C     -0.434987363913     -1.080590767312     -1.253216296583
C     -0.450651404460      0.256716303863     -1.511629042882
H     -1.365703213775      0.657493701779     -1.904548056095
H      0.456114369383      0.713947200024     -1.845719575031
O     -0.587949943075     -2.213282213925     -1.548596677439
C      0.025541836499     -1.299364598826      0.912983504339
H      0.263700194470     -2.343282908643      0.973810587623
H     -1.016000314299     -1.069951304724      1.042009090257
C      0.990112198584     -0.354811421758      1.142291314242
H      1.995147101141     -0.681254211190      1.338580338118
C      0.776156880690      0.999132742296      0.888880990926
H      1.618142779215      1.664207110983      0.872159877414
C     -0.414040256695      1.403840537398      0.330826965246
H     -0.515861868288      2.405651200115     -0.043532153336
H     -1.339720995479      0.941548629919      0.605699133202

遷移状態構造

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

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

Mode                                 0                   1                   2
Freq [cm^-1]                     -954.8974             79.9983            233.8784
Reduced mass [au]                   6.6752              3.4927              4.4277
Force const [Dyne/A]               -3.5861              0.0132              0.1427
Char temp [K]                       0.0000            115.0997            336.4989
Normal mode                   x         y         z            x         y         z            x         y         z
       C                 0.03604   -0.03088    0.16517    0.00642   -0.00348   -0.00583    0.02999   -0.01857    0.06037
       C                -0.00213    0.07027    0.08594    0.12831   -0.00812   -0.02679    0.02412   -0.03620   -0.03114
       H                 0.01640   -0.04123   -0.08585    0.20617    0.05421   -0.14353    0.04798   -0.04838   -0.09446
       H                -0.01009   -0.04156   -0.11997    0.20391   -0.06499    0.10250    0.04749    0.00078    0.09106
       O                -0.00703    0.01252   -0.03758   -0.13494    0.00772    0.02282    0.02516   -0.03234    0.07640
       C                -0.03546    0.00369   -0.10292    0.09157   -0.02174   -0.03761   -0.05646    0.02937   -0.14099
       H                -0.01659    0.00901   -0.02736    0.15756   -0.00832   -0.06278   -0.08647    0.01929   -0.21529
       H                -0.01224   -0.00865    0.15088    0.07656   -0.09484   -0.02683   -0.03952    0.05169   -0.03547
       C                 0.01416    0.04258   -0.00180    0.02655    0.03495    0.00756   -0.01440   -0.01852   -0.10017
       H                -0.00661    0.01542    0.05973    0.04405    0.09584    0.02009   -0.01322   -0.06107   -0.18194
       C                 0.02835   -0.03657    0.00760   -0.05362    0.02766    0.03837   -0.00244    0.00500    0.01913
       H                -0.00958    0.01121    0.05661   -0.09603    0.08259    0.08480    0.00809   -0.00722    0.04858
       C                -0.02538   -0.06312   -0.11163   -0.05704   -0.03389   -0.00056   -0.01024    0.07525    0.10766
       H                -0.00576   -0.03212   -0.03977   -0.10304   -0.03200    0.01700   -0.00474    0.09282    0.15529
       H                -0.02943    0.05606    0.09786   -0.04042   -0.10013   -0.05746   -0.00866    0.03285    0.03342

(...snip...)

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

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

終わりに#

　　自作モジュールで、電子状態計算モジュールであるPySCFを用いて、計算レベルHF/3-21Gでケテンの[4+2]付加環化反応の遷移状態構造を算出する手順を説明した。

参考#

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.

概要#