最終更新：2025-06-28
概要
本記事では、自作モジュール(MultiOptPy)で、分子内Mannich反応の遷移状態構造を算出してみる。電子状態計算ソフトウェアにPySCFを使用する。 MultiOptPyは電子状態計算ソフトウェアを用いた分子構造最適化手法の勉強を目的として作成したpythonモジュールである。
MultiOptPyのレポジトリ:https://github.com/ss0832/MultiOptPy
分子内Mannich反応:
- https://www.chem-station.com/odos/2009/06/mannich.html
- https://en.wikipedia.org/wiki/Mannich_reaction
使用した自作モジュールMultiOptPyのバージョン
1.9.12b環境
WSL2 (Ubuntu-22.04)Source codeのダウンロード
wget https://github.com/ss0832/MultiOptPy/archive/refs/tags/v1.9.12b.zip
unzip v1.9.12b.zip
cd MultiOptPy-1.9.12b手順
1. 初期構造の準備
モデル反応系として、以下の構造を用意した。今回は名前をintramolecular_mannich_rxn.xyzとした。 初期構造は以下のものを使用した。
24
1 0
C 1.32598549 0.40416980 0.08373468
N 0.33091890 -0.42450882 0.35066769
H 1.28333659 1.47128067 0.14973063
C 2.26719821 -1.69636356 0.43014509
H 2.85795254 -2.54060013 0.14174619
H 2.41660449 -1.53099737 1.47667831
C 2.56048760 -0.42239174 -0.33560373
H 2.55309571 -0.58237946 -1.39354951
H 3.50220926 0.01347840 -0.07468700
C -1.01286371 -0.02140769 0.78959570
H -0.94801858 0.89110033 1.34459499
H -1.63851135 0.12518766 -0.06595957
H -1.43018866 -0.78709732 1.40964351
C 0.77881534 -1.82478807 0.14120789
H 0.30048654 -2.53461957 0.78326770
C 0.51724409 -2.24825705 -1.31613713
H 0.91498524 -3.22802284 -1.47972495
H 0.99150745 -1.55545893 -1.97946587
C -0.99936171 -2.25870714 -1.58334054
C -1.61526941 -1.13986850 -2.03657870
H -1.04781445 -0.24922738 -2.20877995
H -2.66901487 -1.14712948 -2.22223368
O -1.75773591 -3.44900353 -1.35320344
H -2.50274792 -3.47571998 -1.95805328
2. 遷移状態構造最適化のための初期構造の算出
※あらかじめpythonで利用できる電子状態計算ソフトウェアのPySCF(導入方法はpipを使えば可能である。)が使える環境を用意しておく。
a. 初期構造をカレントディレクトリにintramolecular_mannich_rxn.xyzとして保存し、以下を実行する。
python optmain.py intramolecular_mannich_rxn.xyz -pyscf -spin 0 -elec 1 -func hf -bs 3-21g -opt rsirfo_fsb -ma 150 1 20-opt rsirfo_fsbは準ニュートン法であるRS-I-RFO法を構造最適化に使用することを示す。初期のへシアンに関しては、特にオプションで指定しない限り、単位行列が使われる。-pyscfは使用する電子状態計算ソフトウェアをPySCFにすることを示す。-spinはスピン多重度の指定である。PySCFを使用するときは目的とするスピン多重度に1を引いた値を指定する。-elec 1は形式電荷を1とすることを示す。-ma 150 1 20は150kJ/molの活性化障壁を超えうるペア同士を近づける力を原子のラベル番号1と20のペアに、構造最適化時に加えることを示す。-func hfはHartree-Fock法を用いることを示す。-bs 3-21gは基底関数にPople系である3-21Gを用いることを示す。
これを実行すると、計算レベルHF/3-21Gで指定した人工力ポテンシャルをかけた上で初期構造を構造最適化することができる。
結果はyyyy_mm_dd（今日の年月日）ディレクトリの中に存在するディレクトリを開いて確認できる。
正常終了していれば、このディレクトリ中に、intramolecular_mannich_rxn_traj.xyzが存在するので、これをコピーして、MultiOptPy-1.9.12bディレクトリに置く。
intramolecular_mannich_rxn_traj.xyzは構造最適化の過程をAvogadro（公式ページ:https://avogadro.cc/ ）等で可視化して確認できるようにしている。このintramolecular_mannich_rxn_traj.xyzは次のNEB計算に使用する。
※intramolecular_mannich_rxn_traj.xyzをアニメーションとして表示したい場合は、[https://github.com/ss0832/molecule_movie] を使うと良い。
（この人工力ポテンシャルを加えて行った構造最適化の結果はintramolecular_mannich_rxn_optimized.xyzで確認できる。構造を可視化して、生成系になっているか確認する。反応系のままであれば、-maの設定を見直してb.をやり直す。今回の場合以下の構造が得られた。生成系の安定構造に近いものが得られていると判断できるため、次のNEB計算を行うことが可能である。）
b. intramolecular_mannich_rxn_traj.xyzを初期パスとして、NEB法で経路の緩和を行う。
NEB法を用いることで、先ほど得られたintramolecular_mannich_rxn_traj.xyz全体のエネルギーを下げることができる。これにより、パスのエネルギー極大値を持つ構造を遷移状態構造に近づける。（この時点ではまだ正確な遷移状態構造は求められていない。）
python nebmain.py intramolecular_mannich_rxn_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 10は10回分NEB法による経路の緩和を行うことを示す。-modelhessは少ない計算コストで算出された近似へシアンを用いることを示す。これを指定すると経路緩和のアルゴリズムが-fcを1以上に指定した場合と同様のものになる。-pyscfは使用する電子状態計算ソフトウェアをPySCFにすることを示す。-spngは緩和中のパスのエネルギープロファイルや各ノードの勾配のRMS値をmatplotlibで可視化するオプションである。-spinはスピン多重度の指定である。PySCFを使用するときは目的とするスピン多重度に1を引いた値を指定する。
c. 初期構造の決定及び遷移状態構造の計算
MultiOptPy-1.9.12bと同じディレクトリ内に、NEBという名前を含むディレクトリが生成されている。 そのディレクトリ内のenergy_plot.csvを確認し、緩和後のパスのエネルギー極大値を示す構造を確認する。
パスの緩和後の各ノードのエネルギー一覧(単位
私が実行した環境では、15番のノード(グラフでは16番)がエネルギー極大値を示していた。 path_ITR_10_intramolecular_mannich_rxnディレクトリ内のintramolecular_mannich_rxn_traj_15.xyzを可視化する。
※こちら[https://ss0832.github.io/molecule_viewer/] を使うことでも可視化は可能である。
intramolecular_mannich_rxn_traj_15.xyz
24
1 0
C 0.228718620580 1.266586237390 0.010714317344
N -0.169485523231 0.454065507566 0.971806521574
H -0.243212246893 2.212493661962 -0.157005879612
C 2.021558733749 -0.247482414576 0.590774099815
H 2.652203147268 -0.949774464621 0.066117330660
H 2.502475236588 -0.005878320964 1.527609697426
C 1.722801756942 1.042159202905 -0.247388529834
H 1.974482253655 0.910918839310 -1.283597336859
H 2.255396763100 1.916180128990 0.114465306677
C -1.484278162726 0.536615701515 1.651836809113
H -1.450418263490 1.289795271618 2.422959391121
H -2.257492092761 0.778712713020 0.940890573975
H -1.710612668960 -0.418507347207 2.099538343241
C 0.603294037818 -0.800888713286 0.844689232947
H 0.485590155687 -1.407982158143 1.723817499523
C 0.077389954638 -1.588900242403 -0.465489737254
H -0.411284544854 -2.502582917494 -0.174120027499
H 0.930237705561 -1.828732534149 -1.083004148989
C -0.930364202809 -0.750286612410 -1.179113913607
C -0.624440520687 0.386245760240 -1.841703618892
H 0.319858007788 0.488867006302 -2.320399052407
H -1.392583502716 1.054676853756 -2.186022650765
O -2.169988076613 -1.094496288207 -0.770953617949
H -2.929846567635 -0.741804871116 -1.256420609749
構造が壊れていないので、これを遷移状態を求めるための初期構造とする。
intramolecular_mannich_rxn_traj_15.xyzをMultiOptPy-1.9.12bと同じディレクトリ内にコピーする。
そして、以下を実行する。
python optmain.py intramolecular_mannich_rxn_traj_15.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と同等）
実行して得られた正確な遷移状態構造を以下に示す。
(実行して得られた正確な遷移状態構造はintramolecular_mannich_rxn_traj_15_optimized.xyzとして保存されている。)
24
OptimizedStructure
C 0.094178841794 1.264453245849 -0.110620372291
N -0.182459075294 0.497853311189 0.979960685002
H -0.483297096659 2.151325212908 -0.276220735598
C 1.983962688094 -0.150012025152 0.320837371686
H 2.410417644166 -0.793278347102 -0.433708852295
H 2.695794114507 -0.053453012342 1.126374720338
C 1.611944112602 1.258544955356 -0.245464464980
H 1.974751093085 1.449603483004 -1.240307855554
H 1.982856249286 2.034772012256 0.413069654068
C -1.448475275713 0.510987517257 1.728713911548
H -1.244654138341 0.218817959834 2.748419631302
H -1.843431966138 1.516871658901 1.735728604547
H -2.180026647187 -0.162336262925 1.303095497761
C 0.634394740064 -0.708271359015 0.863848055356
H 0.730110147814 -1.219959322652 1.806029696028
C -0.067604467813 -1.684096995429 -0.190812778885
H -0.694409612048 -2.409076538516 0.305250059672
H 0.706319672208 -2.200113836890 -0.744564019547
C -0.905468735451 -0.877786621427 -1.116195757644
C -0.381181042427 0.222139682836 -1.769430655312
H 0.617852052845 0.116063022569 -2.137664869135
H -1.029272810047 0.817706924455 -2.388835738042
O -2.197794122539 -1.142337538988 -1.037011032805
H -2.784506366808 -0.658417125974 -1.640490755219
90回程度の反復計算で遷移状態構造が得られた。-freqオプションにより生成されたnormal_modes.txtやvibration_animationディレクトリ内の振動モードのアニメーションを確認した。
以下に-freqオプションで生成されたnormal_modes.txtの一部を示す。
Mode 0 1 2
Freq [cm^-1] -466.6009 100.4336 135.2121
Reduced mass [au] 7.7823 2.6971 2.6981
Force const [Dyne/A] -0.9983 0.0160 0.0291
Char temp [K] 0.0000 144.5015 194.5401
Normal mode x y z x y z x y z
C -0.03800 -0.07350 -0.17082 -0.02053 -0.03605 -0.04250 -0.04263 -0.01963 -0.02106
N 0.00541 -0.00886 0.00243 -0.03979 -0.03335 -0.03947 -0.02125 -0.01849 -0.00626
H 0.03482 0.00496 -0.03305 -0.01196 -0.03015 -0.04199 -0.08201 -0.04777 -0.03640
C -0.00570 -0.00071 -0.00189 -0.00100 0.01115 0.07811 -0.02305 0.01135 -0.09050
H 0.00380 -0.00401 0.00419 0.09028 -0.00156 0.14127 -0.10897 0.04578 -0.16941
H -0.01175 -0.00084 0.00386 -0.07457 0.07640 0.13572 0.05740 -0.05604 -0.15396
C -0.00460 0.00057 -0.00440 -0.01409 -0.02417 -0.00581 -0.04262 0.04720 0.01859
H 0.06965 0.00516 0.02002 0.00749 -0.07146 -0.00675 -0.02209 0.13823 0.04310
H -0.07084 -0.00460 0.03175 -0.03945 0.01621 -0.03927 -0.08827 0.00263 0.09680
C 0.01001 0.00360 -0.00393 -0.01354 0.00188 0.00761 -0.00001 -0.04421 0.03201
H 0.00116 0.00297 -0.00129 0.04900 0.20883 0.05472 0.05498 0.06344 0.05220
H 0.02573 0.01077 0.00411 -0.09244 -0.02744 -0.16401 -0.07317 -0.07225 -0.04503
H -0.02899 0.02129 0.02395 0.02429 -0.13256 0.15982 0.03311 -0.13292 0.11635
C -0.00947 0.00716 0.00099 -0.01895 -0.01729 -0.00680 0.00720 0.00188 -0.00525
H -0.00707 -0.01286 -0.00957 -0.06925 -0.03481 -0.01151 0.06328 0.02771 0.00322
C -0.00288 -0.02269 0.00083 0.06264 -0.01675 -0.08171 -0.03523 -0.02634 0.05793
H 0.02103 -0.04262 0.00171 0.10794 -0.08544 -0.12628 -0.09293 0.05014 0.09727
H -0.00521 0.01156 -0.03666 0.10424 0.07848 -0.10826 -0.06568 -0.11453 0.09738
C -0.01843 -0.01244 0.05143 0.02174 0.00610 -0.01614 0.03134 -0.00936 0.01066
C 0.03430 0.11396 0.15160 0.00466 -0.00160 -0.04493 0.11188 -0.06135 -0.01374
H -0.03698 -0.03965 -0.02354 -0.01307 -0.00580 -0.09144 0.12860 -0.11328 0.04773
H 0.01750 0.07117 0.12047 -0.02005 0.01169 -0.00648 0.17006 -0.06097 -0.07382
O 0.02054 -0.00969 -0.02931 0.01640 0.07942 0.11416 0.01162 0.10060 0.00999
H 0.00017 0.06396 0.04183 -0.01927 0.11405 0.17655 0.05463 0.12632 -0.01164
(...snip...)その結果、虚振動が１つであることが確認できた。
次に、vibration_animation内のmode_1_467i_wave_number.xyzをAvogadroで確認すると、想定される反応系と生成系をつなぐ方向に振動していることを確認できた。
終わりに
自作モジュールで、電子状態計算モジュールであるPySCFを用いて、計算レベルHF/3-21Gで分子内Mannich反応の遷移状態構造を算出する手順を説明した。
参考
- 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.