最終更新：2025-06-03
概要
本記事では、自作モジュール(MultiOptPy)で、Williamsonエーテル合成の遷移状態構造を算出してみる。電子状態計算ソフトウェアにPySCFを使用する。 MultiOptPyは電子状態計算ソフトウェアを用いた分子構造最適化手法の勉強を目的として作成したpythonモジュールである。
MultiOptPyのレポジトリ:https://github.com/ss0832/MultiOptPy
Williamsonエーテル合成：
- https://www.chem-station.com/odos/2009/06/williamson-williamson-ether-sy.html
- https://en.wikipedia.org/wiki/Williamson_ether_synthesis
使用した自作モジュールMultiOptPyのバージョン
1.9.11環境
WSL2 (Ubuntu-22.04)Source codeのダウンロード
wget https://github.com/ss0832/MultiOptPy/archive/refs/tags/v1.9.11.zip
unzip v1.9.11.zip
cd MultiOptPy-1.9.11手順
1. 初期構造の準備
モデル反応系として、以下の構造を用意した。今回は名前をwilliamson_ether_synthesis.xyzとした。
10
-1 0
C -1.95556172 0.11554152 -0.36617882
H -1.59924313 -0.89338715 -0.36607468
H -1.59889175 0.61967080 -1.23998670
H -3.02556164 0.11591071 -0.36596983
O -1.47843441 0.78968654 0.80120201
C -0.74133075 1.80811209 4.59850205
H -0.38467652 0.79930203 4.59850195
H -0.38465784 2.31251063 3.72485078
H -1.81133075 1.80812528 4.59850202
Br -0.10465304 2.70848671 6.15801064
2. 遷移状態構造最適化のための初期構造の算出
※あらかじめpythonで利用できる電子状態計算ソフトウェアのPySCF(導入方法はpipを使えば可能である。)が使える環境を用意しておく。
a. 初期構造をカレントディレクトリにwilliamson_ether_synthesis.xyzとして保存し、以下を実行する。
python optmain.py williamson_ether_synthesis.xyz -opt rsirfo_fsb -func hf -bs 3-21g -pyscf -elec -1 -spin 0 -ma 150 5 6-opt rsirfo_fsbは準ニュートン法であるRS-I-RFO法を構造最適化に使用することを示す。初期のへシアンに関しては、特にオプションで指定しない限り、単位行列が使われる。-pyscfは使用する電子状態計算ソフトウェアをPySCFにすることを示す。-spinはスピン多重度の指定である。PySCFを使用するときは目的とするスピン多重度に1を引いた値を指定する。-elec -1は形式電荷を-1に指定することを示す。-ma 150 5 6は150kJ/molの活性化障壁を超えうるペア同士を近づける力を原子のラベル番号5と6のペアに、構造最適化時にそれぞれ加えることを示す。-func hfはHartree-Fock法を用いることを示す。-bs 3-21Gは基底関数にPople系である3-21Gを用いることを示す。
これを実行すると、計算レベルHF/3-21Gで指定した人工力ポテンシャルをかけた上で初期構造を構造最適化することができる。
結果はyyyy_mm_dd（今日の年月日）ディレクトリの中に存在するディレクトリを開いて確認できる。
正常終了していれば、このディレクトリ中に、williamson_ether_synthesis_traj.xyzが存在するので、これをコピーして、MultiOptPy-1.9.11ディレクトリに置く。
williamson_ether_synthesis_traj.xyzは構造最適化の過程をAvogadro（公式ページ:https://avogadro.cc/ ）等で可視化して確認できるようにしている。このwilliamson_ether_synthesis_traj.xyzは次のNEB計算に使用する。
※williamson_ether_synthesis_traj.xyzをアニメーションとして表示したい場合は、[https://github.com/ss0832/molecule_movie] を使うと良い。
（この人工力ポテンシャルを加えて行った構造最適化の結果はwilliamson_ether_synthesis_optimized.xyzで確認できる。構造を可視化して、生成系になっているか確認する。反応系のままであれば、-maの設定を見直してb.をやり直す。今回の場合以下の構造が得られた。安定構造に近いものが得られていると判断できるため、次のNEB計算を行うことが可能である。）
b. williamson_ether_synthesis_traj.xyzを初期パスとして、NEB法で経路の緩和を行う。
NEB法を用いることで、先ほど得られたwilliamson_ether_synthesis_traj.xyz全体のエネルギーを下げることができる。これにより、パスのエネルギー極大値を持つ構造を遷移状態構造に近づける。（この時点ではまだ正確な遷移状態構造は求められていない。）
LUP法と呼ばれるNEB法とよく似た方法があるが、本モジュールではNEB法を用いている。
python nebmain.py williamson_ether_synthesis_traj.xyz -pyscf -ns 10 -nd 0.3 -elec -1 -spin 0 -spng -func hf -bs 3-21g -fc 10-nd 0.3はノード間の距離を0.3Åとして初期パスを作成することを示す。-ns 10は10回分NEB法による経路の緩和を行うことを示す。-fc 10は10回に一回正確なへシアンを計算させることを指示する。これに伴い、経路の緩和に使うアルゴリズムも変わる。-pyscfは使用する電子状態計算ソフトウェアをPySCFにすることを示す。-spngは緩和中のパスのエネルギープロファイルや各ノードの勾配のRMS値をmatplotlibで可視化するオプションである。-spinはスピン多重度の指定である。PySCFを使用するときは目的とするスピン多重度に1を引いた値を指定する。-elec -1は形式電荷を-1に指定することを示す。
(参考)
c. 初期構造の決定及び遷移状態構造の計算
MultiOptPy-1.9.11と同じディレクトリ内に、NEBという名前を含むディレクトリが生成されている。 そのディレクトリ内のenergy_plot.csvを確認し、緩和後のパスのエネルギー極大値を示す構造を確認する。
パスの緩和後の各ノードのエネルギー一覧(単位
私が実行した環境では、12番のノード(グラフでは13番)がエネルギー極大値を示していた。 path_ITR_10_williamson_ether_synthesisディレクトリ内のwilliamson_ether_synthesis_traj_12.xyzを可視化する。
※こちら[https://ss0832.github.io/molecule_viewer/] を使うことでも可視化は可能である。
williamson_ether_synthesis_traj_12.xyz
10
-1 0
C -0.430116079830 -0.619884815036 -1.631312392463
H -0.128885447215 -1.606892649081 -2.050081255813
H 0.030886997161 0.096535869591 -2.359175235168
H -1.514482658255 -0.568464153221 -1.917196339847
O -0.105685414640 -0.523900343032 -0.306895430117
C 0.320767278300 0.504978126343 1.358859847620
H 0.754170510959 -0.447156051591 1.535549593048
H 0.752669119350 1.091545579627 0.575969053362
H -0.723580549889 0.600280012940 1.536130091877
Br 1.044256244059 1.472958423461 3.258152067501
構造が壊れていないので、これを遷移状態を求めるための初期構造とする。
williamson_ether_synthesis_traj_12.xyzをMultiOptPy-1.9.11と同じディレクトリ内にコピーする。
そして、以下を実行する。
python optmain.py williamson_ether_synthesis_traj_12.xyz -opt rsirfo_bofill -fc 5 -order 1 -func hf -bs 3-21g -pyscf -elec -1 -spin 0 -freq -tcc-opt rsirfo_bofillは遷移状態構造の最適化向けのoptimizerを指定することを意味する。準ニュートン法であるRS-I-RFO法を使用する。今回は-fcで正確なへシアンを計算するようにしているので、初期へシアンは正確なへシアンを使用するようになっている。-order 1は一次の鞍点を求めることを指定する。（デフォルトだと極小値を求めるようになっている。）-fc 5は5回の反復回数当たり1回正確なへシアンを計算することを指定する。-freqは収束条件を満たした後に基準振動解析を行うことを示す。（自前で実装しているため、あくまで目安として使用することを推奨する。各振動モードをvibration_animation内のxyzファイルで可視化できる。）-tccは収束条件を厳しくすることを示す。（Gaussianのtightと同等）
実行して得られた正確な遷移状態構造を以下に示す。
(実行して得られた正確な遷移状態構造はwilliamson_ether_synthesis_traj_12_optimized.xyzとして保存されている。)
10
OptimizedStructure
C -0.328278402043 -0.517255700742 -1.641755786629
H -0.825581204485 -1.504114896147 -1.795185971978
H -0.132239056486 -0.171786777333 -2.683673109537
H -1.177087967778 0.136798576650 -1.324017227574
O 0.765130233220 -0.530946280204 -0.805718668850
C 0.340396175846 0.361341989768 1.461098457777
H 0.914308301611 -0.451288307501 1.854895030467
H 0.764709256729 0.750209065111 0.556067235013
H -0.701140323625 0.135269481081 1.362045072574
Br 0.379782987011 1.791772849318 3.016244968737
60回程度の反復計算で遷移状態構造が得られた。-freqオプションにより生成されたnormal_modes.txtやvibration_animationディレクトリ内の振動モードのアニメーションを確認した。
以下に-freqオプションで生成されたnormal_modes.txtの一部を示す。
Mode 0 1 2
Freq [cm^-1] -197.6935 71.6138 113.4820
Reduced mass [au] 2.7387 2.9029 1.1357
Force const [Dyne/A] -0.0631 0.0088 0.0086
Char temp [K] 0.0000 103.0363 163.2753
Normal mode x y z x y z x y z
C -0.00686 -0.02110 0.01357 0.03953 -0.12612 -0.11446 0.01222 0.01269 -0.04254
H -0.06627 0.00478 0.03497 0.16028 -0.19295 -0.08102 0.30469 -0.06112 -0.51302
H 0.05251 -0.07641 0.00683 0.09625 -0.19124 -0.12648 -0.01845 0.53080 0.12311
H 0.02224 0.04644 -0.04472 -0.06198 -0.19304 -0.24794 -0.19704 -0.36995 0.18924
O -0.03618 -0.05131 0.05026 -0.03480 0.05858 -0.01649 -0.04244 -0.05113 0.02724
C 0.07817 0.09683 -0.17001 -0.10273 0.05455 -0.01014 0.02521 -0.00178 0.00072
H 0.07438 0.07612 -0.22386 -0.13484 0.01182 -0.05297 0.04366 0.01239 0.00160
H 0.14948 0.34431 0.00266 -0.10213 0.05120 -0.01246 0.01437 0.01987 0.00885
H 0.08445 0.13508 -0.25776 -0.11792 0.10917 0.01082 0.02987 -0.01990 0.00119
Br -0.00756 -0.00789 0.01975 0.01871 0.00418 0.02880 0.00065 0.00727 0.00325
(...snip...)その結果、虚振動が１つであることが確認できた。
次に、vibration_animation内のmode_1_198i_wave_number.xyzをAvogadroで確認すると、想定される反応系と生成系をつなぐ方向に振動していることを確認できた。IRC計算すると、想定される反応系と生成系をつないでいることを確認できる。
終わりに
自作モジュールで、電子状態計算モジュールであるPySCFを用いて、計算レベルHF/3-21Gでwilliamsonエーテル合成の遷移状態構造を算出する手順を説明した。
参考
- 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.