最終更新：2025-06-12
概要
本記事では、自作モジュール(MultiOptPy)で、Cannizzaro反応のある1つの素過程の遷移状態構造を算出してみる。電子状態計算ソフトウェアにPySCFを使用する。 MultiOptPyは電子状態計算ソフトウェアを用いた分子構造最適化手法の勉強を目的として作成したpythonモジュールである。
MultiOptPyのレポジトリ:https://github.com/ss0832/MultiOptPy
Cannizzaro反応:
- https://www.chem-station.com/odos/2009/06/cannizzaro-cannizzaro-reaction.html
- https://en.wikipedia.org/wiki/Cannizzaro_reaction
使用した自作モジュールMultiOptPyのバージョン
1.9.12環境
WSL2 (Ubuntu-22.04)Source codeのダウンロード
wget https://github.com/ss0832/MultiOptPy/archive/refs/tags/v1.9.12.zip
unzip v1.9.12.zip
cd MultiOptPy-1.9.12手順
1. 初期構造の準備
モデル反応系として、以下の構造を用意した。今回は名前をcannizzaro_rxn.xyzとした。 以下のコマンドを実行して初期構造を作成した。
python optmain.py cannizzaro_rxn.xyz -opt rsirfo_fsb -func hf -bs 3-21g -pyscf -spin 0 -elec 0 -kp 2.0 3.0 11,3 2.0 3.0 11,4 2.0 3.0 2,11-kp k x a,bは、構造最適化時に原子のラベル番号a,bの原子ペアに対して、力の定数k[a.u.]、平行距離x[angstroms]の調和ポテンシャルを加えることを意味する。-func hfはHartree-Fock法を用いることを示す。-bs 3-21Gは基底関数にPople系である3-21Gを用いることを示す。-opt rsirfo_fsbは準ニュートン法であるRS-I-RFO法を構造最適化に使用することを示す。初期のへシアンに関しては、特にオプションで指定しない限り、単位行列が使われる。-pyscfは使用する電子状態計算ソフトウェアをPySCFにすることを示す。-spinはスピン多重度の指定である。PySCFを使用するときは目的とするスピン多重度に1を引いた値を指定する。-elecは形式電荷の指定が可能である。
cannizzaro_rxn.xyz
14
0 1
C -1.60395977 1.21172897 0.00039014
H -1.24700511 1.71574069 0.87414958
O -3.03395971 1.21178626 0.00082841
O -1.12734732 -0.13650693 -0.00033738
C -1.09098263 1.93821359 -1.25685866
H -1.44760873 2.94703346 -1.25631429
H -0.02098268 1.93817072 -1.25718659
H -1.44793730 1.43420187 -2.13061810
K -3.93052438 3.74797827 0.00219697
K 1.56265255 -0.13661470 -0.00116181
C -1.30373697 0.73089654 3.39961136
H -0.76856835 -0.19565327 3.39961136
H -2.37373695 0.73070183 3.39961136
O -0.67473530 1.82081739 3.39961136
最終的に、以下の構造を初期構造とした。
14
OptimizedStructure
C 0.026279416161 0.367560704256 -0.486434562435
H 1.024895765957 0.635476929396 -0.029914722940
O -0.987599132498 1.051238689133 0.177129401473
O -0.116066607519 -1.044781523307 -0.463901481489
C 0.142770526156 0.863594865051 -1.961301629867
H 0.273389290879 1.938930404855 -1.984482969668
H 0.970469300222 0.382654481837 -2.479163596355
H -0.772117378708 0.618535530326 -2.501894198151
K -2.586647812030 -0.661863766583 -0.566455695504
K 2.100716912603 -1.751666341885 0.233902208028
C -0.029727707131 -0.583516322765 2.500812741353
H -0.731284543384 -0.762503636827 3.318488379562
H -0.360235941816 0.117693100216 1.731535805260
O 1.045157911108 -1.171353113704 2.511680320731
2. 遷移状態構造最適化のための初期構造の算出
※あらかじめpythonで利用できる電子状態計算ソフトウェアのPySCF(導入方法はpipを使えば可能である。)が使える環境を用意しておく。
a. 初期構造をカレントディレクトリにcannizzaro_rxn.xyzとして保存し、以下を実行する。
python optmain.py cannizzaro_rxn.xyz -opt rsirfo_fsb -func hf -bs 3-21g -pyscf -spin 0 -elec 0 -ma 100 2 11 100 10 14 -ns 100-opt rsirfo_fsbは準ニュートン法であるRS-I-RFO法を構造最適化に使用することを示す。初期のへシアンに関しては、特にオプションで指定しない限り、単位行列が使われる。-pyscfは使用する電子状態計算ソフトウェアをPySCFにすることを示す。-spinはスピン多重度の指定である。PySCFを使用するときは目的とするスピン多重度に1を引いた値を指定する。-ma 100 2 11 100 10 14は100kJ/molの活性化障壁を超えうるペア同士を近づける力を原子のラベル番号2と11のペアとラベル番号10と14のペアに、構造最適化時にそれぞれ加えることを示す。-func hfはHartree-Fock法を用いることを示す。-bs 3-21Gは基底関数にPople系である3-21Gを用いることを示す。-ns 100は構造最適化時の電子状態計算の反復回数を最大100回とすることを示す。(デフォルトは1000回である。)
これを実行すると、計算レベルHF/3-21Gで指定した人工力ポテンシャルをかけた上で初期構造を構造最適化することができる。
結果はyyyy_mm_dd（今日の年月日）ディレクトリの中に存在するディレクトリを開いて確認できる。
正常終了していれば、このディレクトリ中に、cannizzaro_rxn_traj.xyzが存在するので、これをコピーして、MultiOptPy-1.9.12ディレクトリに置く。
cannizzaro_rxn_traj.xyzは構造最適化の過程をAvogadro（公式ページ:https://avogadro.cc/ ）等で可視化して確認できるようにしている。このcannizzaro_rxn_traj.xyzは次のNEB計算に使用する。
※cannizzaro_rxn_traj.xyzをアニメーションとして表示したい場合は、[https://github.com/ss0832/molecule_movie] を使うと良い。
（この人工力ポテンシャルを加えて行った構造最適化の結果はcannizzaro_rxn_optimized.xyzで確認できる。構造を可視化して、生成系になっているか確認する。反応系のままであれば、-maの設定を見直してb.をやり直す。今回の場合以下の構造が得られた。生成系の安定構造に近いものが得られていると判断できるため、次のNEB計算を行うことが可能である。）
b. cannizzaro_rxn_traj.xyzを初期パスとして、NEB法で経路の緩和を行う。
NEB法を用いることで、先ほど得られたcannizzaro_rxn_traj.xyz全体のエネルギーを下げることができる。これにより、パスのエネルギー極大値を持つ構造を遷移状態構造に近づける。（この時点ではまだ正確な遷移状態構造は求められていない。）
python nebmain.py cannizzaro_rxn_traj.xyz -ns 10 -modelhess -func hf -bs 3-21g -pyscf -spin 0 -elec 0 -spng -nd 0.5-nd 0.5はノード間の距離を0.5Åとして初期パスを作成することを示す。-ns 10は10回分NEB法による経路の緩和を行うことを示す。-modelhessは少ない計算コストで算出された近似へシアンを用いることを示す。これを指定すると経路緩和のアルゴリズムが-fcを1以上に指定した場合と同様のものになる。-pyscfは使用する電子状態計算ソフトウェアをPySCFにすることを示す。-spngは緩和中のパスのエネルギープロファイルや各ノードの勾配のRMS値をmatplotlibで可視化するオプションである。-spinはスピン多重度の指定である。PySCFを使用するときは目的とするスピン多重度に1を引いた値を指定する。
c. 初期構造の決定及び遷移状態構造の計算
MultiOptPy-1.9.12と同じディレクトリ内に、NEBという名前を含むディレクトリが生成されている。 そのディレクトリ内のenergy_plot.csvを確認し、緩和後のパスのエネルギー極大値を示す構造を確認する。
パスの緩和後の各ノードのエネルギー一覧(単位
私が実行した環境では、14番のノード(グラフでは15番)がエネルギー極大値を示していた。 path_ITR_10_cannizzaro_rxnディレクトリ内のcannizzaro_rxn_traj_14.xyzを可視化する。
※こちら[https://ss0832.github.io/molecule_viewer/] を使うことでも可視化は可能である。
cannizzaro_rxn_traj_14.xyz
14
0 0
C -0.288541323143 -0.049079848011 -0.505020589300
H 0.491414120902 0.057349638018 0.348638147216
O -1.543882192816 0.170318663545 -0.009778230838
O -0.059438411463 -1.297631345637 -1.087752122013
C 0.185672765982 1.118763806634 -1.398287286611
H 0.155517437909 2.056560527600 -0.853993036896
H 1.182635442376 0.941298193164 -1.792804684745
H -0.491484472386 1.195633238647 -2.248188309962
K -2.541370688546 -1.404588356781 -1.536780090118
K 1.206451071703 -2.623383637827 0.610269244800
C 0.421189045388 -0.022712918632 2.072415292355
H 1.136540165271 0.790398558381 2.134655127539
H -0.599257412459 0.230584915352 1.811043117503
O 0.744554451281 -1.163511434452 2.455583421071
構造が壊れていないので、これを遷移状態を求めるための初期構造とする。
cannizzaro_rxn_traj_14.xyzをMultiOptPy-1.9.12と同じディレクトリ内にコピーする。
そして、以下を実行する。
python optmain.py cannizzaro_rxn_traj_14.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と同等）
実行して得られた正確な遷移状態構造を以下に示す。
(実行して得られた正確な遷移状態構造はcannizzaro_rxn_traj_14_optimized.xyzとして保存されている。)
14
OptimizedStructure
C -0.397254934533 -0.003094695702 -0.578761265595
H 0.421595723176 0.033960496597 0.472813955017
O -1.576310793480 0.378661974508 -0.109577722219
O -0.276254384584 -1.263674386922 -1.095832166096
C 0.318544659059 1.061261684906 -1.441551922820
H 0.295560058900 2.018462736294 -0.939139994682
H 1.342688855756 0.771007967975 -1.648092732498
H -0.210397746951 1.147780435339 -2.388504379272
K -2.801204623399 -1.448012208315 -1.182541513821
K 1.529626482394 -2.278451818761 0.285728231418
C 0.343819757078 -0.099586407985 2.060868935726
H 0.789066421622 0.882172253117 2.205675547131
H -0.722133693146 -0.065040526159 1.883474542284
O 0.942654218108 -1.135447504892 2.475440485426
30回程度の反復計算で遷移状態構造が得られた。-freqオプションにより生成されたnormal_modes.txtやvibration_animationディレクトリ内の振動モードのアニメーションを確認した。
以下に-freqオプションで生成されたnormal_modes.txtの一部を示す。
Mode 0 1 2
Freq [cm^-1] -1046.1007 36.9421 70.9369
Reduced mass [au] 1.7831 11.1831 9.7594
Force const [Dyne/A] -1.1497 0.0090 0.0289
Char temp [K] 0.0000 53.1514 102.0624
Normal mode x y z x y z x y z
C 0.10122 0.01230 0.08050 -0.02367 -0.02337 0.03408 -0.01516 -0.04710 0.01937
H -0.20597 -0.02154 -0.62826 -0.00798 -0.01902 0.01119 -0.03414 -0.03590 0.02291
O -0.03247 0.00940 0.01914 -0.02254 -0.05278 0.06133 -0.01058 0.00274 -0.01101
O -0.00345 -0.03958 -0.02973 -0.00933 -0.01357 0.01181 -0.04291 -0.06507 0.05614
C -0.01511 0.00737 -0.03557 -0.06247 0.00209 0.03206 0.02898 -0.09559 -0.00734
H -0.00215 0.01258 -0.02710 -0.07084 -0.00415 0.04387 0.04559 -0.07655 -0.04313
H -0.00139 0.01904 -0.03678 -0.06157 0.02426 0.00685 0.02619 -0.12870 0.02098
H -0.01634 0.01355 -0.03215 -0.08375 0.00423 0.04374 0.04968 -0.11702 -0.02083
K 0.00343 0.00290 -0.00065 -0.00903 0.01838 -0.08317 -0.05672 0.04855 -0.02957
K -0.00101 -0.00213 -0.00719 -0.02264 -0.02541 0.02379 0.07937 0.04846 -0.00560
C -0.01543 -0.02429 0.10991 0.07630 0.03782 0.00349 -0.03250 -0.03752 0.01782
H 0.07555 -0.04229 -0.15501 0.05960 0.05422 -0.05501 -0.07628 -0.01694 0.01277
H 0.04549 -0.04663 -0.21472 0.06464 0.00599 0.06423 -0.03494 -0.08091 0.02345
O -0.01640 0.03586 -0.01755 0.12272 0.06695 0.01203 0.01383 -0.01005 0.01712
(...snip...)その結果、虚振動が１つであることが確認できた。
次に、vibration_animation内のmode_1_1046i_wave_number.xyzをAvogadroで確認すると、想定される反応系と生成系をつなぐ方向に振動していることを確認できた。
終わりに
自作モジュールで、電子状態計算モジュールであるPySCFを用いて、計算レベルHF/3-21GでCannizzaro反応のある1つの素過程の遷移状態構造を算出する手順を説明した。
参考
- 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.