最終更新：2025-07-07
概要
本記事では、自作モジュール(MultiOptPy)で、Darzens縮合反応のある１つの素過程の遷移状態構造を算出してみる。電子状態計算ソフトウェアにPySCFを使用する。 MultiOptPyは電子状態計算ソフトウェアを用いた分子構造最適化手法の勉強を目的として作成したpythonモジュールである。
MultiOptPyのレポジトリ:https://github.com/ss0832/MultiOptPy
Darzens縮合反応:
- https://www.chem-station.com/odos/2009/07/darzens-darzens-condensation.html
- https://en.wikipedia.org/wiki/Darzens_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. 初期構造の準備
モデル反応系として、以下の構造を用意した。今回は名前をdarzens_condensation.xyzとした。 初期構造は以下のものを使用した。
23
-1 0
C -2.82118419 2.03205358 0.17807529
O -3.35800107 1.11027829 -0.48955374
C -1.32037284 1.96408752 0.51652160
H -0.75798531 2.44141918 -0.25857500
H -1.14140484 2.46301931 1.44600399
H -1.01895686 0.94060090 0.59729005
C -3.66505182 3.22806633 0.65665751
H -4.07529005 3.01199223 1.62097933
H -3.04635568 4.09881409 0.71921496
H -4.45983333 3.40438854 -0.03770035
C -3.81618886 -0.01477644 5.12430622
Cl -3.00300293 1.40092105 5.78167161
C -5.29667857 0.07114131 4.70913354
H -5.51992821 1.06722306 4.38838326
H -5.48231789 -0.61197222 3.90676583
H -5.91644156 -0.18213080 5.54378741
C -3.04723684 -1.33942951 4.96428419
O -2.43395438 -1.58554020 3.89335715
O -3.03012184 -2.28979325 6.03265483
C -2.91327435 -3.60960947 5.49475912
H -2.00534693 -3.68708979 4.93389853
H -2.90046803 -4.32072080 6.29416932
H -3.74657675 -3.80857262 4.85372810
2. 遷移状態構造最適化のための初期構造の算出
※あらかじめpythonで利用できる電子状態計算ソフトウェアのPySCF(導入方法はpipを使えば可能である。)が使える環境を用意しておく。
a. 初期構造をカレントディレクトリにdarzens_condensation.xyzとして保存し、以下を実行する。
python optmain.py darzens_condensation.xyz -pyscf -spin 0 -elec -1 -func hf -bs 3-21g -opt rsirfo_fsb -ma 150 1 11 -ns 55-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を用いることを示す。-ns 55は電子状態計算の反復回数を最大55回までとすることを示す。これは、塩化物イオンが解離したため、計算効率を上げるために指定した
これを実行すると、計算レベルHF/3-21Gで指定した人工力ポテンシャルをかけた上で初期構造を構造最適化することができる。
結果はyyyy_mm_dd（今日の年月日）ディレクトリの中に存在するディレクトリを開いて確認できる。
正常終了していれば、このディレクトリ中に、darzens_condensation_traj.xyzが存在するので、これをコピーして、MultiOptPy-1.9.12cディレクトリに置く。
darzens_condensation_traj.xyzは構造最適化の過程をAvogadro（公式ページ:https://avogadro.cc/ ）等で可視化して確認できるようにしている。このdarzens_condensation_traj.xyzは次のNEB計算に使用する。
※darzens_condensation_traj.xyzをアニメーションとして表示したい場合は、[https://github.com/ss0832/molecule_movie] を使うと良い。
（この人工力ポテンシャルを加えて行った構造最適化の結果はdarzens_condensation_optimized.xyzで確認できる。構造を可視化して、生成系になっているか確認する。反応系のままであれば、-maの設定を見直してb.をやり直す。今回の場合以下の構造が得られた。生成系の安定構造に近いものが得られていると判断できるため、次のNEB計算を行うことが可能である。）
b. darzens_condensation_traj.xyzを初期パスとして、NEB法で経路の緩和を行う。
NEB法を用いることで、先ほど得られたdarzens_condensation_traj.xyz全体のエネルギーを下げることができる。これにより、パスのエネルギー極大値を持つ構造を遷移状態構造に近づける。（この時点ではまだ正確な遷移状態構造は求められていない。）
python nebmain.py darzens_condensation_traj.xyz -ns 10 -modelhess -func hf -bs 3-21g -pyscf -spin 0 -elec -1 -spng -nd 0.5-nd 0.5はノード間の距離を0.5Åとして初期パスを作成することを示す。-ns nはn回分NEB法による経路の緩和を行うことを示す。-modelhessは少ない計算コストで算出された近似へシアンを用いることを示す。これを指定すると経路緩和のアルゴリズムが-fcを1以上に指定した場合と同様のものになる。-pyscfは使用する電子状態計算ソフトウェアをPySCFにすることを示す。-spngは緩和中のパスのエネルギープロファイルや各ノードの勾配のRMS値をmatplotlibで可視化するオプションである。-spinはスピン多重度の指定である。PySCFを使用するときは目的とするスピン多重度に1を引いた値を指定する。
c. 初期構造の決定及び遷移状態構造の計算
MultiOptPy-1.9.12cと同じディレクトリ内に、NEBという名前を含むディレクトリが生成されている。 そのディレクトリ内のenergy_plot.csvを確認し、緩和後のパスのエネルギー極大値を示す構造を確認する。
パスの緩和後の各ノードのエネルギー一覧(単位
※bias_force_rms.csvにて、各Iterationごとのすべてのノードの勾配のRMS値を確認できる。
私が実行した環境では、15番のノード(グラフでは16番)がエネルギー極大値を示していた。 path_ITR_10_darzens_condensationディレクトリ内のdarzens_condensation_traj_15.xyzから遷移状態構造を求めようとしたがうまくいかなかったため、エネルギー極大値を示す構造に近い勾配のRMS値が極小値を示す構造darzens_condensation_traj_17.xyzを使用した。
※こちら[https://ss0832.github.io/molecule_viewer/] を使うことでも可視化は可能である。
darzens_condensation_traj_17.xyz
23
-1 0
C 0.055758079459 0.945928957797 -1.547089412206
O -0.517445717911 0.074084895000 -2.262018852897
C 1.580968529070 1.092854153662 -1.599524642599
H 1.808947097102 1.664713674955 -2.501081136257
H 1.966328740741 1.628563305108 -0.745284757005
H 2.049471544965 0.126184606453 -1.667485273928
C -0.614683517509 2.357889942462 -1.465304327182
H -1.007123676301 2.612362382048 -0.492600060698
H 0.097943870615 3.128768606438 -1.755224734910
H -1.411432689978 2.356706865846 -2.203156474340
C -0.154807161280 0.196756073080 0.321848285481
Cl 0.896713108982 1.012537473129 1.892737659788
C -1.594806870248 0.168177402436 0.817307200349
H -2.052292972813 1.150607885359 0.865541657199
H -2.094654138322 -0.441857799269 0.072110087109
H -1.714037572663 -0.332478211657 1.762220985011
C 0.446008194061 -1.209532302881 0.423100231432
O 1.533183141745 -1.547871567567 0.042236829999
O -0.382679823960 -2.142941605150 1.053774884865
C 0.175867246096 -3.036922419240 2.046353737169
H 1.254593418172 -3.059979970233 2.017250902562
H -0.143626297689 -2.686424978406 3.019682342288
H -0.178192532333 -4.058127369370 1.904604868769
構造が壊れていないので、これを遷移状態を求めるための初期構造とする。
darzens_condensation_traj_17.xyzをMultiOptPy-1.9.12cと同じディレクトリ内にコピーする。
そして、以下を実行する。
python optmain.py darzens_condensation_traj_17.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と同等）
実行して得られた正確な遷移状態構造を以下に示す。
(実行して得られた正確な遷移状態構造はdarzens_condensation_traj_17_optimized.xyzとして保存されている。)
23
OptimizedStructure
C 0.179778772189 0.416807657107 -1.137467178081
O -0.238966648809 -0.725384004841 -1.474081770529
C 1.685307416057 0.685302152405 -1.074381407062
H 2.067827951342 0.644551192227 -2.093466314573
H 1.897841744887 1.652988280169 -0.647809356003
H 2.181349349383 -0.071212682676 -0.486401638428
C -0.593246162497 1.654124726529 -1.605498859013
H -0.369758993686 2.519673502757 -1.000411260906
H -0.304927506417 1.856147405881 -2.636965419179
H -1.653953655908 1.446194151326 -1.592797796916
C -0.242303080129 0.592053584680 0.905978907804
Cl 0.484167373803 2.174337148239 1.743754412167
C -1.744716340104 0.605024765613 1.053432833235
H -2.146132393348 1.559547135339 0.741516984456
H -2.122667316852 -0.170641181318 0.400097249793
H -2.081226297868 0.405639290368 2.067290836337
C 0.452158427448 -0.571510024830 1.403439376639
O 1.604943787377 -0.705271270736 1.764470587881
O -0.367463279002 -1.675713503838 1.271516635081
C 0.234919023387 -2.864140347898 0.699606563069
H 0.386919821079 -2.675237911542 -0.349648656372
H 1.166272877333 -3.090354988499 1.196829073746
H -0.476124869667 -3.662925076462 0.850996196853
160回程度の反復計算で遷移状態構造が得られた。-freqオプションにより生成されたnormal_modes.txtやvibration_animationディレクトリ内の振動モードのアニメーションを確認した。
以下に-freqオプションで生成されたnormal_modes.txtの一部を示す。
Mode 0 1 2
Freq [cm^-1] -289.9622 58.1750 116.1395
Reduced mass [au] 6.9217 3.5713 4.9138
Force const [Dyne/A] -0.3429 0.0071 0.0391
Char temp [K] 0.0000 83.7008 167.0988
Normal mode x y z x y z x y z
C -0.02892 0.03621 0.18277 -0.01823 -0.00928 -0.00254 -0.01032 -0.01784 0.01937
O -0.00935 0.00537 0.07394 -0.11320 0.03088 -0.01699 0.02161 -0.02468 0.00227
C -0.00229 0.02438 0.01403 0.00268 -0.12516 -0.00843 -0.01615 0.01228 0.03557
H -0.11608 0.00285 -0.02536 -0.01423 -0.08569 -0.01645 0.00306 -0.02892 0.04475
H 0.06590 0.02572 -0.01785 0.07678 -0.16755 0.05138 -0.03806 0.03366 -0.00406
H 0.04945 0.01503 -0.03982 -0.03855 -0.19600 -0.06488 -0.01246 0.04654 0.07838
C 0.00345 0.01459 0.01803 0.07648 0.06150 0.02916 -0.04089 -0.04117 0.01071
H -0.00798 0.08187 -0.06906 0.14459 0.02928 0.04897 -0.07402 -0.03191 0.00934
H 0.03163 -0.10255 0.00532 0.08860 0.06315 0.03279 -0.03720 -0.03610 0.01283
H -0.00057 0.02765 0.00770 0.06059 0.14291 0.02747 -0.03422 -0.07644 0.00083
C 0.02979 0.00303 -0.14929 0.00184 -0.01044 0.00979 -0.01165 -0.00566 0.01753
Cl -0.00360 -0.01873 -0.02650 -0.02629 0.01296 -0.01291 0.04138 0.00587 -0.04844
C 0.02771 -0.01687 -0.04919 -0.00096 -0.03490 -0.02259 -0.00867 0.03639 0.05344
H 0.00248 -0.01533 -0.00658 -0.00687 -0.03275 -0.00830 0.00260 0.02272 -0.00234
H 0.01034 -0.01333 -0.04174 0.02970 -0.02332 -0.05442 -0.04557 -0.00147 0.12036
H 0.11135 -0.04057 -0.03015 -0.02267 -0.06598 -0.03615 0.01488 0.11857 0.07811
C 0.00851 -0.02164 -0.03331 0.01538 0.00890 0.02820 -0.02704 -0.00568 0.03808
O -0.01330 0.01259 0.00893 0.00840 0.03014 0.05796 -0.05214 0.00341 0.11985
O 0.00718 -0.00251 0.01844 0.03714 -0.00953 0.00978 -0.00226 -0.01504 -0.05172
C -0.01343 -0.00604 -0.01713 0.04930 0.02336 -0.04735 0.04325 0.03633 -0.12178
H -0.03495 0.04578 -0.01046 -0.06536 0.03238 -0.06243 0.04622 0.11056 -0.10644
H -0.01195 -0.02716 -0.02930 0.11205 0.08295 -0.13719 0.04301 0.03659 -0.12067
H -0.02449 0.00407 -0.02530 0.11533 -0.02168 0.02554 0.06712 0.00438 -0.17779
(...snip...)その結果、虚振動が１つであることが確認できた。
次に、vibration_animation内のmode_1_290i_wave_number.xyzをAvogadroで確認すると、想定される反応系と生成系をつなぐ方向に振動していることを確認できた。
終わりに
自作モジュールで、電子状態計算モジュールであるPySCFを用いて、計算レベルHF/3-21GでDarzens縮合反応の１つの素過程の遷移状態構造を算出する手順を説明した。
参考
- 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.