最終更新：2025-05-29
概要
本記事では、自作モジュール(MultiOptPy)で、Swern酸化の反応機構の１段階目の遷移状態構造を算出してみる。電子状態計算ソフトウェアにPySCFを使用する。 MultiOptPyは電子状態計算ソフトウェアを用いた分子構造最適化手法の勉強を目的として作成したpythonモジュールである。
MultiOptPyのレポジトリ:https://github.com/ss0832/MultiOptPy
Swern酸化:
- https://www.chem-station.com/odos/2009/06/swern-swern-oxidation.html
- https://en.wikipedia.org/wiki/Swern_oxidation
使用した自作モジュールMultiOptPyのバージョン
1.9.7c環境
WSL2 (Ubuntu-22.04)Source codeのダウンロード
wget https://github.com/ss0832/MultiOptPy/archive/refs/tags/v1.9.7c.zip
unzip v1.9.7c.zip
cd MultiOptPy-1.9.7c手順
1. 初期構造の準備
モデル反応系として、以下の構造を用意した。今回は名前をswern_oxid_stage_1.xyzとした。 以下のコマンドを実行して構造最適化して得られた構造を使用した。
python optmain.py swern_oxid_stage_1.gjf -kp 2.0 3.0 2,11 -opt rsirfo_fsb -pyscf -spin 0-kp k x a,bは、構造最適化時に原子のラベル番号a,bの原子ペアに対して、力の定数k[a.u.]、平行距離x[angstroms]の調和ポテンシャルを加えることを意味する。
16
OptimizedStructure
S 0.546379587477 -1.002002566876 0.688052941288
O 0.776662285723 0.064125459690 -0.369692119643
C 0.818793354591 -2.606477145280 -0.152255450331
H 0.500207746095 -3.430076725096 0.493591791581
H 0.265663313207 -2.610192763063 -1.096164315979
H 1.889984270316 -2.683511866603 -0.353939710658
C -1.267311919755 -1.172185381210 0.897178760738
H -1.736345367934 -1.299117872221 -0.083004143002
H -1.491121440468 -2.019018676927 1.554036476395
H -1.620533241104 -0.244986166053 1.355288706515
C -0.511821078786 2.769111389464 -0.240230578458
C 0.958548387328 2.868290457416 -0.733175737744
O -0.866107828051 2.877869203332 0.888511930904
O 1.293010256418 2.822414481989 -1.869701530863
Cl -1.639000627131 2.504766991218 -1.590015756800
Cl 2.082992302076 3.160991180223 0.611518736056
2. 遷移状態構造最適化のための初期構造の算出
※あらかじめpythonで利用できる電子状態計算ソフトウェアのPySCF(導入方法はpipを使えば可能である。)が使える環境を用意しておく。
a. 初期構造をカレントディレクトリにswern_oxid_stage_1.xyzとして保存し、以下を実行する。
python optmain.py swern_oxid_stage_1.xyz -ma 100 2 12 50 1 16 -opt rsirfo_fsb -pyscf -spin 0-opt rsirfo_fsbは準ニュートン法であるRS-I-RFO法を構造最適化に使用することを示す。初期のへシアンに関しては、特にオプションで指定しない限り、単位行列が使われる。-pyscfは使用する電子状態計算ソフトウェアをPySCFにすることを示す。-spinはスピン多重度の指定である。PySCFを使用するときは目的とするスピン多重度に1を引いた値を指定する。-ma 100 2 12 50 1 16は100kJ/molの活性化障壁を超えうる力を、原子のラベル番号2と12のペアに構造最適化時に加える。そして、50kJ/molの活性化障壁を超えうる力を原子のラベル番号1と16のペアに、構造最適化時に加えることを示す。
これを実行すると、-bs（基底関数の指定）と-func（計算手法の指定）がないので、デフォルトの設定である計算レベルB3LYP/6-31G(d)で指定した人工力ポテンシャルをかけた上で初期構造を構造最適化することができる。
結果はyyyy_mm_dd（今日の年月日）ディレクトリの中に存在するディレクトリを開いて確認できる。
正常終了していれば、このディレクトリ中に、swern_oxid_stage_1_traj.xyzが存在するので、これをコピーして、MultiOptPy-1.9.7cディレクトリに置く。
swern_oxid_stage_1_traj.xyzは構造最適化の過程をAvogadro（公式ページ:https://avogadro.cc/ ）等で可視化して確認できるようにしている。このswern_oxid_stage_1_traj.xyzは次のNEB計算に使用する。
※swern_oxid_stage_1_traj.xyzをアニメーションとして表示したい場合は、[https://github.com/ss0832/molecule_movie] を使うと良い。
（この人工力ポテンシャルを加えて行った構造最適化の結果はswern_oxid_stage_1_optimized.xyzで確認できる。構造を可視化して、生成系になっているか確認する。反応系のままであれば、-maの設定を見直してb.をやり直す。今回の場合以下の構造が得られた。中間体の安定構造に近いものが得られていると判断できるため、次のNEB計算を行うことが可能である。）
b. swern_oxid_stage_1_traj.xyzを初期パスとして、NEB法で経路の緩和を行う。
NEB法を用いることで、先ほど得られたswern_oxid_stage_1_traj.xyz全体のエネルギーを下げることができる。これにより、パスのエネルギー極大値を持つ構造を遷移状態構造に近づける。（この時点ではまだ正確な遷移状態構造は求められていない。）
LUP法と呼ばれるNEB法とよく似た方法があるが、本モジュールではNEB法を用いている。
python nebmain.py swern_oxid_stage_1_traj.xyz -ns 10 -spng -nd 0.4 -pyscf -spin 0-nd 0.4はノード間の距離を0.4Åとして初期パスを作成することを示す。-ns 10は10回分NEB法による経路の緩和を行うことを示す。-pyscfは使用する電子状態計算ソフトウェアをPySCFにすることを示す。-spngは緩和中のパスのエネルギープロファイルや各ノードの勾配のRMS値をmatplotlibで可視化するオプションである。
(参考)
c. 初期構造の決定及び遷移状態構造の計算
MultiOptPy-1.9.7cと同じディレクトリ内に、NEBという名前を含むディレクトリが生成されている。 そのディレクトリ内のenergy_plot.csvを確認し、緩和後のパスのエネルギー極大値を示す構造を確認する。
パスの緩和後の各ノードのエネルギー一覧(単位
私が実行した環境では、10番のノードがエネルギー極大値を示していた。 path_ITR_10_swern_oxid_stage_1ディレクトリ内のswern_oxid_stage_1_traj_10.xyzを可視化する。
※こちら[https://ss0832.github.io/molecule_viewer/] を使うことでも可視化は可能である。
swern_oxid_stage_1_traj_10.xyz
16
0 0
S 0.556326436138 -0.612926030472 0.615284206072
O 0.467592030298 0.247866259600 -0.774120084406
C 0.986200982390 -2.153189255544 -0.229352007119
H 0.835504888774 -2.973027435033 0.475452105858
H 0.379657356242 -2.268147731421 -1.129107887489
H 2.042991261176 -2.067859994514 -0.490008272449
C -1.219145344112 -0.858131044006 0.923539060760
H -1.725398562064 -1.095398454770 -0.015229627223
H -1.332196926329 -1.658771323747 1.658355055595
H -1.601523334886 0.078057153434 1.337569060992
C -0.614982790620 2.392667687189 -0.246254805320
C 0.693931747175 1.734653375248 -0.765824580076
O -0.952583313330 2.655943583974 0.858742047942
O 1.289825616994 2.142933215626 -1.718001722302
Cl -1.721790907173 2.699201919162 -1.653756295210
Cl 1.915590859329 1.736128075272 1.152713744376
構造が壊れていないので、これを遷移状態を求めるための初期構造とする。
swern_oxid_stage_1_traj_10.xyzをMultiOptPy-1.9.7cと同じディレクトリ内にコピーする。
そして、以下を実行する。
python optmain.py swern_oxid_stage_1_traj_10.xyz -opt rsirfo_bofill -pyscf -spin 0 -order 1 -fc 5 -freq -grid 5-opt rsirfo_bofillは遷移状態構造の最適化向けのoptimizerを指定することを意味する。準ニュートン法であるRS-I-RFO法を使用する。今回は-fcで正確なへシアンを計算するようにしているので、初期へシアンは正確なへシアンを使用するようになっている。-order 1は一次の鞍点を求めることを指定する。（デフォルトだと極小値を求めるようになっている。）-fc 5は5回の反復回数当たり1回正確なへシアンを計算することを指定する。-freqは収束条件を満たした後に基準振動解析を行うことを示す。（自前で実装しているため、あくまで目安として使用することを推奨する。各振動モードをvibration_animation内のxyzファイルで可視化できる。）-gridはDFT計算の数値積分計算時の格子の細かさを指定できる。1~10の間で指定可能であり、数値が大きいほど格子が細かくなる。今回は塩化物イオンが遊離する初期構造を用いたため、DFT計算時の格子が粗い場合、うまくPESを描けない場合がある。そこで、デフォルトの3よりも細かくした。
実行して得られた正確な遷移状態構造を以下に示す。
(実行して得られた正確な遷移状態構造はswern_oxid_stage_1_traj_18_optimized.xyzとして保存されている。)
16
OptimizedStructure
S 0.551900409702 -0.683127948792 0.481254908658
O -0.117489512701 0.356894425565 -0.514687975853
C 0.853850131375 -2.082013082258 -0.633651956401
H 1.167371096977 -2.944830910382 -0.038729786162
H -0.050986561439 -2.303145779360 -1.204332810277
H 1.661000926403 -1.777919055379 -1.303493534917
C -0.878223062869 -1.300524327450 1.409601331626
H -1.665334188222 -1.599877057695 0.714001796410
H -0.556065621724 -2.138060360165 2.035290963414
H -1.206486889589 -0.463067603899 2.028266612393
C -0.313651322228 2.838373508025 -1.060596961407
C 0.634351604251 1.855879818700 -0.340823273343
O 0.091951762491 3.858008746370 -1.512648794976
O 1.813494823755 1.818417110479 -0.532181530548
Cl -2.046004231938 2.376731482295 -1.186052639834
Cl 0.060320635757 2.188261033946 1.658783651218
140回程度の反復計算で遷移状態構造が得られた。-freqオプションにより生成されたnormal_modes.txtやvibration_animationディレクトリ内の振動モードのアニメーションを確認した。
以下に-freqオプションで生成されたnormal_modes.txtの一部を示す。
Mode 0 1 2
Freq [cm^-1] -150.8228 44.8310 50.6461
Reduced mass [au] 10.3422 6.8775 6.2218
Force const [Dyne/A] -0.1386 0.0081 0.0094
Char temp [K] 0.0000 64.5017 72.8685
Normal mode x y z x y z x y z
S -0.02827 -0.03466 -0.02588 0.01012 -0.00969 -0.01334 0.01978 -0.00085 -0.00135
O -0.04223 -0.06544 0.00463 -0.03712 -0.00112 0.03094 -0.07226 0.00458 0.05911
C 0.00345 -0.05407 0.00787 -0.07705 -0.01315 -0.03255 0.04103 0.04943 -0.05750
H 0.02702 -0.04025 0.01538 -0.05094 -0.01812 -0.05359 0.10396 0.04422 -0.09820
H 0.00451 -0.07991 0.01563 -0.11674 0.00024 0.02513 0.03019 0.02642 -0.03140
H -0.00984 -0.04590 -0.00461 -0.11604 -0.02350 -0.08416 0.00234 0.11074 -0.07642
C 0.01179 -0.07136 0.01721 0.06611 0.00108 0.08097 0.07900 -0.08041 0.03329
H 0.00167 -0.08540 0.03400 0.01734 0.00708 0.13353 0.06366 -0.09498 0.05711
H 0.03179 -0.06310 0.01767 0.10153 -0.00158 0.05924 0.13610 -0.08231 0.00131
H -0.00019 -0.07975 0.02195 0.11286 0.00296 0.10319 0.06943 -0.10815 0.06630
C -0.00316 0.03173 -0.00130 -0.00814 0.01852 0.02760 0.01168 -0.01392 -0.01649
C 0.01661 0.13148 0.11772 -0.01125 0.00657 0.01323 -0.03564 -0.03690 0.01341
O -0.00908 0.01030 -0.05745 -0.01794 0.05869 0.10967 0.05892 -0.04705 -0.04893
O -0.00062 0.06901 0.06136 -0.01568 -0.01380 -0.00960 -0.03235 -0.07859 0.02087
Cl 0.01268 0.00221 0.00775 0.01304 -0.03014 -0.09137 -0.00920 0.05956 -0.00653
Cl 0.02549 0.02153 -0.03943 0.02198 0.01545 0.00773 -0.03264 0.02766 0.00527
(...snip...)その結果、虚振動が１つであることが確認できた。
次に、vibration_animation内のmode_1_151i_wave_number.xyzをAvogadroで確認すると、想定されるような反応系と生成系をつなぐ方向に振動していることを確認できた。
終わりに
自作モジュールで、電子状態計算モジュールであるPySCFを用いて、計算レベルB3LYP/6-31G(d) (grid dense level = 5)でSwern酸化の反応機構の一段階目の遷移状態構造を算出する手順を説明した。
参考
- 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.