5903 words
30 minutes
【計算化学】自作pythonライブラリで遷移状態構造を求めてみる(PySCF使用, フルオロメタンと水酸化物イオンが関わるSN2反応, 基底関数の選択ミスによる失敗例)
最終更新：2026-01-27
概要
本記事では、自作ライブラリ(MultiOptPy)で、フルオロメタンと水酸化物イオンが関わるSN2反応の素過程の遷移状態構造を算出してみる。計算レベルは、ωB97X-D/STO-3Gとした。
低廉な基底関数を使った場合、SN2反応の遷移状態構造に見られるWalden反転が見られず、化学的に妥当な遷移状態構造が得られない失敗例に関して示す。
MultiOptPyは電子状態計算ソフトウェアを用いた分子構造最適化手法の勉強を目的として作成したpythonライブラリである。
MultiOptPyのレポジトリ:https://github.com/ss0832/MultiOptPy
使用した自作ライブラリMultiOptPyのバージョン
v1.20.5Video Demo
https://www.youtube.com/watch?v=AE61iY2HZ8Y
環境
Ubuntu (WSL2) 環境構築手順(Linux)
すでにAnacondaの導入が終わっている前提で、環境の構築手段について述べる。以下を実行すると可能である。
## 1. Download and install MultiOptPy:
wget https://github.com/ss0832/MultiOptPy/archive/refs/tags/v1.20.5.zip
unzip v1.20.5.zip
cd MultiOptPy-1.20.5
## 2. Create and activate a conda environment:
conda create -n <env-name> python=3.12 pip
conda activate <env-name>
pip install git+https://github.com/ss0832/MultiOptPy.git@v1.20.5
# or pip install multioptpy
## 3. Copy the test configuration file and run the AutoTS workflow（任意で行う）:
cp test/config_autots_run_xtb_test.json .
python run_autots.py aldol_rxn.xyz -cfg config_autots_run_xtb_test.json
手順
1. 初期構造の準備
モデル反応系として、以下の構造を用意した。今回はファイルの名前をSN2_3.xyzとした。
初期構造は以下のものを使用した。
7
-1 1
C 0.000000 0.000000 0.000000
F 1.330000 0.000000 0.000000
H -0.358125 1.040070 0.000000
H -0.358125 -0.520035 -0.900727
H -0.358125 -0.520035 0.900727
O -2.878752 -0.129174 -0.021944
H -3.668163 -0.589323 -0.3475122. 遷移状態構造最適化
run_autots.pyを適切に使用することで、自動的に遷移状態構造が得られる。以下にその手順を説明していく。
初期構造をMultiOptPy-1.20.5ディレクトリ内にSN2_3.xyzとして保存する。その後、同じディレクトリ内で、config_SN2_3.jsonを作成し、以下のように記述する。
config_SN2_3.json
{
"work_dir": "SN2_3",
"top_n_candidates": 3,
"step1_settings": {
"usextb": "GFN2-xTB",
"opt_method": ["rsirfo_block_fsb"],
"use_model_hessian": "fischerd3",
"spin_multiplicity": 0,
"electronic_charge": -1,
"manual_AFIR": ["130", "1", "6"]
},
"step2_settings": {
"NSTEP": 20,
"basisset": "STO-3G",
"functional": "wB97XD",
"use_model_hessian": "fischerd3",
"save_pict": true,
"node_distance": 0.50,
"align_distances_energy_predicted": 2,
"spin_multiplicity": 0,
"electronic_charge": -1,
"pyscf": true,
"dft_grid": 6
},
"step3_settings": {
"opt_method": ["rsirfo_block_bofill"],
"basisset": "STO-3G",
"functional": "wB97XD",
"calc_exact_hess": 5,
"tight_convergence_criteria": true,
"max_trust_radius": 0.2,
"frequency_analysis": true,
"spin_multiplicity": 0,
"electronic_charge": -1,
"pyscf": true,
"dft_grid": 6
},
"step4_settings": {
"opt_method": ["rsirfo_block_bofill"],
"basisset": "STO-3G",
"functional": "wB97XD",
"spin_multiplicity": 0,
"electronic_charge": -1,
"calc_exact_hess": 10,
"tight_convergence_criteria": true,
"frequency_analysis": true,
"pyscf": true,
"intrinsic_reaction_coordinates": ["0.5", "200", "lqa"],
"dft_grid": 6,
"step4b_opt_method": ["rsirfo_block_fsb"]
}
}その後、以下のコマンドを実行する。
python run_autots.py SN2_3.xyz -cfg config_SN2_3.jsonこれにより、これまでの似た内容の記事で行ってきたコマンドの操作をまとめ、遷移状態構造を求める処理を自動的に行う。
具体的な処理の流れは、
Step1. バイアスポテンシャルによるNEB法のための初期経路の作成
Step2. NEB法等による経路の緩和
Step3. NEB法等により得られた経路のエネルギー極大値を示す構造のうち、エネルギー値が上位の最大で3個
(`run_autots.py`にて、`--top_n X`で最大値を変更可能)の構造を初期構造とした遷移状態構造の算出
(Step4.得られた遷移状態構造に対するIRC計算とIRC経路の末端に存在する構造に対する構造最適化。
こちらは`--run_step4`をコマンドで追記しなければ行わない。)となっている。
run_autots.pyのオプションの説明：
-cfg YYY.jsonは、workflowを実行するためのオプションが記されたJSONファイルの読み込み先を指定する。
これらの一連の結果は、(jsonファイルの"work_dir"にて指定した名前)のディレクトリの中に存在するファイルを開いて確認できる。
以下にすべてのstepで共通のオプションに関する説明を載せる。
"opt_method": ["rsirfo_block_fsb"]は準ニュートン法であるRS-I-RFO法を構造最適化に使用することを示す。初期のへシアンに関しては、特にオプションで指定しない限り、単位行列が使われる。(以前のHessian更新法とは細かな点で異なる方法を使用している。具体的には、複数の座標変位や勾配変位を用いてHessianの更新を行う。)"spin_multiplicity": Zはスピン多重度の指定である。PySCFを使用するときは目的とするスピン多重度に1を引いた値を指定する。（デフォルトでは1が指定される。）"electronic_charge": 0は形式電荷をMとすることを示す。（デフォルトでは0が指定される。）"othersoft": "uma-s-1p1"は今回使用するNNPを指定している。これを使用する際にASEライブラリが必要である。"use_model_hessian": "fischerd3"は、計算コストが非常に低い数式を使用して、近似したHessianを生成する機能を呼び出すオプションである。デフォルトではこの機能は使用されない。"grid": 6はDFTの計算における数値積分の格子の細かさを指定するオプションである。デフォルトでは3だが、これでは積分精度としては不十分なので、6を指定した。(これでも同じ条件でSTO-3Gを使った際に勾配が収束せずに求められない遷移状態構造が存在した。, 整数値と実際の数値積分グリッドの細かさの対応はPySCFのドキュメントを参照)"pyscf": trueは量子化学計算エンジンにPySCFを使用することを指定する。"functional": "wB97XD"は計算手法の指定を行う。自作モジュールでのPySCFの使用において、DFTの汎関数の指定か、HF法の指定が可能である。"basisset": "XXX"は基底関数を指定可能である。
※オプションの説明はMultiOptPy-v1.20.5/OPTION_README.mdにて示されている。
Step 1
Step1では、omolのデータセットを使用したuma-s-1p1モデルのNNPで得たエネルギーに対して、指定した人工力ポテンシャルを加えた上で初期構造を構造最適化を行っている。
以下のJSON内で記述したバイアスポテンシャルで、次の経路緩和アルゴリズムの初期経路として用いるトラジェクトリーを生成する。
"manual_AFIR": ["yyy", "a", "b]：yyykJ/molの活性化障壁を超えうるペア同士を近づける力を原子のラベル番号aとbのペアに構造最適化時に加えることを示す。"shape_conditions": ["yyy", "lt", "a,b"]: 構造最適化中に、ラベル番号aの原子とラベル番号bの原子の間の距離yyy(Å)よりも大きくなった時に構造最適化を途中で打ち切る。“lt”を”gt”にすると、yyy(Å)よりも小さくなった時に途中で打ち切る。"shape_conditions": ["yyy", "lt", "a,b,c"]: 構造最適化中に、角度の中心がラベル番号bの原子で、ラベル番号aの原子とラベル番号cで作る角度yyy(degrees)よりも大きくなった時に構造最適化を途中で打ち切る。“lt”を”gt”にすると、yyy(degrees)よりも小さくなった時に途中で打ち切る。"shape_conditions": ["yyy", "lt", "a,b,c,d"]: 構造最適化中に、ラベル番号bとラベル番号cを軸とした、ラベル番号aの原子とラベル番号dの原子がなす二面角yyy(degrees)よりも大きくなった時に構造最適化を途中で打ち切る。“lt”を”gt”にすると、yyy(degrees)よりも小さくなった時に途中で打ち切る。"keep_pot": ["c", "d", "a,b"]:力の定数ca.u.で、平行距離dÅの調和ポテンシャルをa番とb番の原子ペアにかける。初期経路生成時に開裂を防ぎたい結合を保持するとき等に使用する。"usextb": "GFN2-xTB"はGFN2-xTBを使用することを指定する。QSM法(NEB法の類似法)にて初期構造を手早く生成するときに使用可能である。
Step1が正常終了していれば作成されたwork_dirディレクトリ中に、SN2_3_step1_traj.xyzが存在する。必要に応じて確認し、目的に沿った初期経路が得られているか確認する。もし想定とは異なる場合は、プロセスをkillして再度設定を見直してやり直す。
SN2_3_step1_traj.xyzは構造最適化の過程をAvogadro（公式ページ:https://avogadro.cc/ ）等で可視化して確認できるようにしている。このSN2_3_step1_traj.xyzはStep2の経路緩和計算に使用している。
※SN2_3_step1_traj.xyzをアニメーションとして表示したい場合は、[https://github.com/ss0832/molecule_movie] を使うと良い。
Step 2
Step2では、NEB法を用いることで、先ほど得られたSN2_3_step1_traj.xyz全体のエネルギーを下げることができる。これにより、パスのエネルギー極大値を持つ構造を遷移状態構造に近づける。（この時点ではまだ正確な遷移状態構造は求められていない。）
Step2固有のオプションについて以下に示す。
-
"NSTEP": nはn回分NEB法による経路の緩和を行うことを示す。 -
"save_pict": trueは緩和中のパスのエネルギープロファイルや各ノードの勾配のRMS値をmatplotlibで可視化するオプションである。 -
"node_distance": yyy: 入力された経路を経路座標上でyyy(Å)間隔で線形補間により構造を再配置して初期経路とする。 -
"align_distances_energy_predicted": a:経路緩和a回に1回、経路座標上で等間隔に線形補間により構造を置きなおす。エネルギー極大値を示す構造に対しては前後のノードの情報を使って経路座標を変数とした多項式を作り、真のエネルギー極大値の場所を連続最適化により推定し、再配置する。
MultiOptPy-v1.20.5/"work_dir"と同じディレクトリ内に、NEBという名前を含むディレクトリが生成されている。
そのディレクトリ内のenergy_plot.csvを確認し、緩和後のパスのエネルギー極大値を示す構造を確認する。
経路の緩和後の各ノードのエネルギー一覧(単位
※bias_force_rms.csvにて、各Iterationごとのすべてのノードの勾配のRMS値を確認できる。
経路緩和の結果、以下の構造がstep3の初期構造として自動的に用いられた。“work_dir”内のSN2_3_step3_TS_Opt_Inputs内に保存されたSN2_3_ts_guess_X.xyzにて確認が可能である。ts_guessの番号が小さい順にエネルギー値が高い構造を示すようになっている。
※こちら[https://ss0832.github.io/molecule_viewer/] を使うことでも可視化は可能である。
SN2_3_ts_guess_1.xyz
7
-1 0
C 0.104191995868 -0.015591503387 -0.013616728667
F 2.072907669223 0.137958697783 -0.217064533011
H 0.129799862604 1.075534348909 -0.214074663690
H 0.188331498640 -0.595581016301 -0.956581942548
H 0.915626455978 -0.311584678653 0.703617842902
O -1.346777705740 -0.364383877632 0.654834501443
H -2.064079776573 0.073648029281 0.042885523571SN2_3_ts_guess_2.xyz
7
-1 0
C 0.634484803738 0.086703015381 -0.121256394883
F 2.126226533018 0.118808615208 -0.294437194405
H 0.304444351097 1.138189970802 -0.046597012938
H 0.200056683408 -0.408404430271 -1.009723321019
H 0.361819290422 -0.464375997892 0.797412728026
O -1.441402143943 0.071928650820 0.591909876953
H -2.185629517740 -0.542849824049 0.082691318267※この初期構造から、初めに示したjsonファイルの条件では遷移状態構造を求められなかった。
SN2_3_ts_guess_3.xyz
7
-1 0
C -0.127617172233 0.030433251439 -0.576856471881
F 1.528822514988 -0.020840699996 1.216379555276
H -0.224258560207 1.017159998075 -1.078593119989
H -0.133204898446 -0.748463666205 -1.364227995699
H 0.917672265455 0.021461398211 0.238597198223
O -1.249033914179 -0.160485108305 0.341415703822
H -0.712380235378 -0.139265173220 1.223285130249Step 3
step3のオプションで、追加での説明を要するものを以下に示す。
"opt_method": ["rsirfo_block_bofill"]は遷移状態構造の最適化向けのoptimizerを指定することを意味する。準ニュートン法であるRS-I-RFO法を使用する。今回は-fcで正確なHessianを計算するようにしているので、初期Hessianは正確なHessianを使用するようになっている。(Bofill法によるHessianの更新法を細かい点で変更している。具体的には、複数の座標変位や勾配変位を用いてHessianの更新を行う。)"saddle_order": 1は一次の鞍点を求めることを指定する。（step3のデフォルトでは一次の鞍点を指定する。それ以外の値の指定は、プログラムの使用目的上想定していないので、行わないことを勧める。）"calc_exact_hess": 5は5回の反復回数当たり1回正確なHessianを計算することを指定する。"frequency_analysis": trueは収束条件を満たした後に基準振動解析を行うことを示す。（自前で実装しているため、あくまで目安として使用することを推奨する。各振動モードをvibration_animation内のxyzファイルで可視化できる。）UMAモデルから算出されるHessianは数値微分により求めているため、原子数Zが多いとZの二乗オーダーで計算コストが急増する。"tight_convergence_criteria": trueは収束条件を厳しくすることを示す。（Gaussianのtightと同等）"max_trust_radius": Dは一回の反復計算ごとの計算されるステップ幅の最大値をDÅ以下にすることを示す。デフォルトでは、"saddle_order": 1を指定すると0.1Åが指定される。"detect_negative_eigenvalues": trueは、初めの計算時(ITR. 0)に、任意の次数の鞍点(遷移状態構造等)を求める際に、正確なへシアンから算出した固有値に１つも負の固有値がない場合、計算を打ち切るオプションである。
実行して得られた正確な遷移状態構造と思われる構造を以下に示す。
(実行して得られた正確な遷移状態構造は計算開始時に、MultiOptPy-1.20.5/"work_dir"ディレクトリ内に生成された新規ディレクトリ内のSN2_3_ts_final_X.xyzとして保存されている。)
SN2_3_ts_final_1.xyz
7
OptimizedStructure
C 0.152331745794 -0.027690262394 0.060609300070
F 2.020335931437 0.163729510192 -0.332913285506
H 0.075493419620 1.069177254486 -0.073355486903
H 0.082408927229 -0.629177915651 -0.866905420953
H 1.021748785663 -0.323157445636 0.700511422550
O -1.391417215014 -0.309951930365 0.651223303676
H -1.960901594728 0.057070789367 -0.13916983293411回の反復計算により、停留点に収束した構造が得られた。"frequency_analysis": trueオプションにより生成されたnormal_modes.txtやvibration_animationディレクトリ内の振動モードのアニメーションを確認した。
以下に"frequency_analysis": trueオプションで生成されたnormal_modes.txtの一部を示す。
Mode 0 1 2
Freq [cm^-1] -1786.1665 233.0252 360.3721
Reduced mass [au] 1.2488 1.1151 3.7278
Force const [Dyne/A] -2.3474 0.0357 0.2852
Char temp [K] 0.0000 335.2713 518.4950
Normal mode x y z x y z x y z
C -0.11456 -0.02066 0.04322 -0.00024 0.07326 0.03445 -0.06682 0.08094 -0.17458
F 0.03068 -0.00621 0.01356 0.00013 -0.03055 -0.01437 0.00476 -0.03493 0.07509
H 0.02944 -0.02285 -0.05800 0.06015 0.10384 0.23849 -0.09546 0.07873 -0.19780
H 0.02914 0.05936 -0.01956 -0.06150 0.24890 -0.07314 -0.09610 0.10104 -0.18852
H 0.63702 0.25884 -0.54853 0.00051 -0.09911 -0.04619 -0.13831 0.04248 -0.09449
O 0.00737 0.00456 -0.00970 -0.00013 0.01425 0.00659 0.07516 -0.02743 0.05991
H -0.02684 -0.00467 0.00973 0.00344 -0.77631 -0.36308 -0.15719 -0.09230 0.19332
(...snip...)その結果、虚振動が１つであることが確認できた。つまりこの構造は遷移状態構造である。
SN2_3_ts_final_3.xyz
7
OptimizedStructure
C -0.093029583234 0.065928373308 -0.623030526026
F 1.566714260982 -0.045551286493 1.299654585562
H -0.334930026356 1.003496671813 -1.191915085041
H -0.218032746998 -0.741124740919 -1.393890647724
H 1.002114584940 0.024177370584 0.371155029282
O -1.226120713253 -0.120275915807 0.329555346651
H -0.696715776082 -0.186650472485 1.2084712972967回の反復計算により、停留点に収束した構造が得られた。"frequency_analysis": trueオプションにより生成されたnormal_modes.txtやvibration_animationディレクトリ内の振動モードのアニメーションを確認した。
以下に"frequency_analysis": trueオプションで生成されたnormal_modes.txtの一部を示す。
Mode 0 1 2
Freq [cm^-1] -195.6074 392.6833 452.7282
Reduced mass [au] 5.7323 1.1241 1.1015
Force const [Dyne/A] -0.1292 0.1021 0.1330
Char temp [K] 0.0000 564.9836 651.3748
Normal mode x y z x y z x y z
C -0.00920 -0.00856 0.06861 0.00573 -0.08554 -0.00990 -0.00226 0.03450 0.00402
F 0.13788 0.00725 0.01707 -0.00233 0.03419 0.00395 0.00119 -0.01973 -0.00235
H 0.05695 -0.00393 0.04528 0.39753 0.28166 0.41766 -0.02252 0.03187 -0.01257
H 0.05648 0.00120 0.04567 -0.43748 0.31442 -0.34864 0.01830 0.03083 0.01995
H 0.00817 -0.01278 0.11420 0.01536 -0.23060 -0.02667 0.00061 -0.00876 -0.00115
O -0.14483 0.00009 -0.08493 0.00025 -0.00302 -0.00034 0.00448 -0.06505 -0.00749
H -0.31263 -0.02072 0.00397 -0.00371 0.05642 0.00655 -0.06294 0.93952 0.10901
(...snip...)その結果、虚振動が１つであることが確認できた。つまりこの構造は遷移状態構造である。
次に、vibration_animation内の虚振動を示す分子振動が示されたxyzファイル(mode_1_XXXi_wave_number.xyz)をAvogadroで確認すると、求められた遷移状態構造の中に、想定される反応系と生成系をつなぐ方向に振動している構造が存在しないことが確認された。
NEB法等により緩和した経路から得られた遷移状態構造を求めるための初期構造 (No.2)はDFTの数値積分のグリッドが粗かったのか、収束せずに遷移状態構造が求められない結果となった。
求められた遷移状態構造No.1を以下のコマンドでIRC計算を行った結果、2つの安定構造に近い末端構造が得られた。
optmain SN2_3_ts_guess_1_optimized.xyz -opt rsirfo_block_bofill -func wb97xd -bs sto-3g -freq -tcc -irc lqa 0.1 200 -order 1 -fc 10 -pyscf -spin 0 -elec -1 -grid 6SN2_3_ts_guess_1_optimized_0_irc_endpoint_1.xyz
7
Terminal structure
C -0.7142618496 0.1972242792 -0.4283007243
F 1.9004852825 0.1220282322 -0.2445608213
H -1.0918801331 1.2436041415 -0.5993258017
H -1.0854561752 -0.3453445621 -1.3418610715
H 0.8281307221 0.1383114016 -0.2887868061
O -1.4804709104 -0.3440977408 0.7234410842
H -2.4759004336 -0.2241471414 0.4583322283
SN2_3_ts_guess_1_optimized_0_irc_endpoint_2.xyz
7
Terminal structure
C 0.0001566281 0.0010429701 -0.0022239162
F 1.5426329718 0.2047098585 -0.4247241214
H -0.2113493316 1.1045376655 -0.0719813084
H -0.2042167216 -0.6548211852 -0.8940621885
H 0.2834552177 -0.4593687645 0.9855785421
O -1.6846416455 -0.2514879560 0.5235929426
H -2.2132797914 0.1295662728 -0.2964256937SN2_3_ts_guess_1_optimized_0_irc_endpoint_2.xyzを以下のコマンドの条件により構造最適化を行った。
optmain SN2_3_ts_guess_1_optimized_0_irc_endpoint_2.xyz -opt rsirfo_block_fsb -func wb97xd -bs sto-3g -freq -tcc -fc 50 -pyscf -spin 0 -elec -1 -grid 6構造最適化の結果が以下のとおりである。
7
OptimizedStructure
C 0.354160984490 -0.010510188263 0.025550903159
F 1.892658332905 0.196375377346 -0.403852212329
H 0.142989576311 1.093813439003 -0.041965930474
H 0.150149365799 -0.668631998001 -0.865480332509
H 0.640088628859 -0.470721671632 1.012975538726
O -1.327937948402 -0.261177304889 0.547360047711
H -1.852108939962 0.120852346437 -0.274588014284基準振動解析の結果が以下のとおりである。
Mode 0 1 2
Freq [cm^-1] 294.2041 438.2861 472.6989
Reduced mass [au] 1.0594 3.1362 3.4793
Force const [Dyne/A] 0.0540 0.3550 0.4581
Char temp [K] 423.2940 630.5959 680.1083
Normal mode x y z x y z x y z
C 0.00017 -0.03867 -0.01810 -0.06782 0.07384 -0.15912 -0.00080 0.17842 0.08319
F -0.00003 0.01987 0.00928 -0.00310 -0.03854 0.08267 0.00020 -0.07760 -0.03613
H -0.07918 -0.06820 -0.27732 -0.07115 0.07104 -0.26362 0.00832 0.17293 0.09067
H 0.08048 -0.25591 0.12584 -0.07145 0.15578 -0.22404 -0.00986 0.18070 0.07398
H -0.00067 0.15901 0.07427 -0.07617 0.06527 -0.14138 -0.00144 0.34474 0.16094
O 0.00008 -0.03308 -0.01543 0.08010 -0.02253 0.04918 0.00058 -0.09139 -0.04269
H -0.00327 0.77607 0.36276 -0.18659 -0.08728 0.18482 -0.00043 0.09043 0.04247
標準的な有機化学の教科書でよくみられるSN2反応の遷移状態構造に近い構造が安定構造として得られた。これは化学的に妥当ではない。
終わりに
自作ライブラリでフルオロメタンと水酸化物イオンが関わるSN2反応の化学的に妥当ではない遷移状態構造を算出する手順を説明した。
参考
- https://github.com/ss0832/MultiOptPy (自作ライブラリMultiOptPyのレポジトリ)
- https://avogadro.cc/ (Avogadro、分子構造可視化ツール)
- https://ai.meta.com/blog/meta-fair-science-new-open-source-releases/ (UMAの公開に関する記事)
- https://github.com/facebookresearch/fairchem (FAIR Chemistryの提供するGitHubのレポジトリ)
- https://fair-chem.github.io/ (同上のレポジトリの内容に関して説明したサイト)
- https://huggingface.co/facebook/UMA (NNPの配布サイト, Hugging Faceへのアカウント登録と配布元の使用許諾が必要である。)
- arXiv preprint arXiv:2505.08762 (2025). (プレプリント)
- 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.
- Recent developments in the PySCF program package, Q. Sun, X. Zhang, S. Banerjee, P. Bao, M. Barbry, N. S. Blunt, N. A. Bogdanov, G. H. Booth, J. Chen, Z.-H. Cui, J. J. Eriksen, Y. Gao, S. Guo, J. Hermann, M. R. Hermes, K. Koh, P. Koval, S. Lehtola, Z. Li, J. Liu, N. Mardirossian, J. D. McClain, M. Motta, B. Mussard, H. Q. Pham, A. Pulkin, W. Purwanto, P. J. Robinson, E. Ronca, E. R. Sayfutyarova, M. Scheurer, H. F. Schurkus, J. E. T. Smith, C. Sun, S.-N. Sun, S. Upadhyay, L. K. Wagner, X. Wang, A. White, J. Daniel Whitfield, M. J. Williamson, S. Wouters, J. Yang, J. M. Yu, T. Zhu, T. C. Berkelbach, S. Sharma, A. Yu. Sokolov, and G. K.-L. Chan, J. Chem. Phys. 153, 024109 (2020)
- PySCF: the Python-based simulations of chemistry framework, Q. Sun, T. C. Berkelbach, N. S. Blunt, G. H. Booth, S. Guo, Z. Li, J. Liu, J. McClain, S. Sharma, S. Wouters, and G. K.-L. Chan, WIREs Comput. Mol. Sci. 8, e1340 (2018)
- Libcint: An efficient general integral library for Gaussian basis functions, Q. Sun, J. Comp. Chem. 36, 1664 (2015)
- J. Chem. Theory Comput. 2019, 15, 3, 1652–1671 https://doi.org/10.1021/acs.jctc.8b01176
補足: SN2_3_ts_guess_2.xyzからDFTの数値積分の格子をより細かくして計算した結果(-grid 6→8にした。)
コマンド：
optmain SN2_3_ts_guess_2.xyz -opt rsirfo_block_bofill -func wb97xd -bs sto-3g -freq -tcc -order 1 -fc 5 -pyscf -spin 0 -elec -1 -grid 8結果を表したログ
frequencies:
[ 0.000000000000+1316.666704746570j 302.256827512645 +0.000000000000j
380.886456007748 +0.000000000000j 392.600992206038 +0.000000000000j
447.167665767012 +0.000000000000j 767.827393136594 +0.000000000000j
891.375135820368 +0.000000000000j 1241.777064021346 +0.000000000000j
1288.674166455666 +0.000000000000j 1498.201512227227 +0.000000000000j
1585.107688772386 +0.000000000000j 3194.269810330159 +0.000000000000j
3272.995358649674 +0.000000000000j 3392.250493271638 +0.000000000000j
3409.150284430916 +0.000000000000j]
Iteration: 10
Energy ratio (actual/predicted): [[-0.151912473645]]
Adaptive factor: 1.505016106616674
Current trust radius: 0.22578948509973876
Decrease trust radius
user_difined_min_trust_radius: 0.01
user_difined_max_trust_radius: 0.1
==================================================
RS-I-RFO Iteration 10
==================================================
Block Bofill update method (Mean Weight)
Gradient norm: 0.000029
Converged: Gradient norm 0.000029 below threshold 0.000100
Converged: Energy change 0.000000 below threshold 0.000001
Condition number of Hessian: 176.97, Ill-conditioned: False
Gradient norm (0.000029) < threshold (0.010000). Using ADAPTIVE trust radius.
Trust radius adjusted: 0.100000 → 0.025000
Found 1 negative eigenvalues (target for saddle order: 1)
Safeguard Solver: Converged in 3 iterations
Initial step with alpha=1.000000 has norm=0.000091
Initial step is within trust radius (0.025000), using it directly
Predicted energy change: -0.000000
Step quality assessment: poor (avg ratio: 6.032)
==================================================================================
trust radii (unit. ang.): 0.15002463037244232
step radii (unit. ang.): 9.093147351666438e-05
Optimizer instance 0: <multioptpy.Optimizer.rsirfo.RSIRFO object at 0x7af9a4ae7a10>
EIGENVALUES (NORMAL COORDINATE, NUMBER OF VALUES: 15):
-0.15104490 0.00590499 0.01501166 0.03689555 0.04080889 0.07791045
0.08967299 0.10981702 0.14008097 0.14439681 0.34070427 0.47451875
0.79411398 1.02108768 1.04501840
==================================================================================
caluculation results (unit a.u.):
Value Threshold
ENERGY : -212.226053833864
BIAS ENERGY : -212.226053833864
Maximum Force : 0.000012607051 0.000015000000
RMS Force : 0.000006415373 0.000010000000
Maximum Displacement : 0.000049140187 0.000062392949
RMS Displacement : 0.000019842874 0.000043584627
ENERGY SHIFT : 0.000000002516
BIAS ENERGY SHIFT : 0.000000002516
=====================================================
converged!!!
=====================================================
C 0.546814628839 0.109998777279 -0.085543963331
F 2.051053351342 -0.010846123021 -0.366170590968
H 0.415605465513 1.205445997160 0.003996530452
H 0.104636895203 -0.328417296548 -0.999394172348
H 0.225298223759 -0.417348433912 0.830930953932
O -1.476869479291 0.151854016198 0.550964704504
H -1.866539085364 -0.710686937156 0.065216537760
====================================================
Performing vibrational analysis...
====================================================
Is Exact Hessian calculated? : True
Point group: C1
Temperature 298.1500 [K]
Pressure 101325.00 [Pa]
Rotational constants [GHz] 584.03369 24.32345 24.01078
Symmetry number 1
Zero-point energy (ZPE) 0.0502925459 [Eh]
tot elec trans rot vib
Entropy [Eh/K] 0.0002612196 0.0000000000 0.0000600984 0.0000278802 0.0001732409
Cv [Eh/K] 0.0000312850 0.0000000000 0.0000047502 0.0000047502 0.0000217846
Cp [Eh/K] 0.0000344518 0.0000000000 0.0000079170 0.0000047502 0.0000217846
Internal energy [Eh] -212.1685916480 -212.2260538339 0.0014162773 0.0014162773 0.0546296312
Total internal energy [Eh] -212.1685916480
Electronic energy [Eh] -212.2260538339
Enthalpy [Eh] -212.1676474632 -212.2260538339 0.0023604622 0.0014162773 0.0546296312
Total enthalpy [Eh] -212.1676474632
Gibbs free energy [Eh] -212.2455300740 -212.2260538339 -0.0155578851 -0.0068962098 0.0029778547
Total Gibbs free energy [Eh] -212.2455300740
Mode 0 1 2
Freq [cm^-1] -1319.3994 304.1942 380.5508
Reduced mass [au] 1.2268 1.3758 5.1129
Force const [Dyne/A] -1.2583 0.0750 0.4363
Char temp [K] 0.0000 437.6676 547.5277
Normal mode x y z x y z x y z
C 0.11357 -0.00309 -0.02996 -0.01369 -0.10012 -0.05813 -0.07173 0.11111 -0.09535
F -0.00171 0.00004 0.00045 0.01214 0.05080 0.02573 -0.06939 -0.04695 0.06858
H -0.03789 -0.01867 -0.04803 -0.12057 -0.10265 -0.21968 0.00463 0.13423 -0.19378
H -0.02628 0.06348 0.00477 0.04558 -0.27455 -0.00490 0.00326 0.06343 -0.11549
H -0.84225 0.25522 -0.12011 0.01490 -0.02826 -0.00552 -0.04879 0.15352 -0.06214
O -0.02392 -0.01631 0.03073 0.02115 0.01079 0.05697 0.12998 -0.05169 0.02092
H -0.03377 -0.00513 0.02388 -0.34145 0.46870 -0.46681 0.14006 0.03115 -0.11796
Mode 3 4 5
Freq [cm^-1] 392.0658 447.4030 770.1983
Reduced mass [au] 1.8139 5.0881 1.1987
Force const [Dyne/A] 0.1643 0.6001 0.4190
Char temp [K] 564.0952 643.7131 1108.1434
Normal mode x y z x y z x y z
C -0.04187 -0.10831 -0.10423 0.01901 0.07632 -0.13181 -0.03989 -0.06191 -0.02809
F -0.00658 0.04847 0.05444 0.08985 -0.03735 0.03422 0.01791 0.01627 0.00510
H -0.03056 -0.12674 0.14801 0.01099 0.07043 -0.13527 0.04054 -0.05609 -0.02694
H 0.01928 0.11099 -0.24520 0.01041 0.10662 -0.13696 -0.03074 -0.03632 -0.04346
H -0.09070 -0.29197 -0.22746 0.22840 0.00616 -0.11829 -0.07682 -0.06469 -0.04525
O 0.05027 0.05312 0.00387 -0.12997 -0.02180 0.07804 -0.03883 0.04382 0.04368
H -0.07332 -0.15939 0.47804 -0.10723 -0.04200 0.07621 0.82068 -0.10781 -0.33928
Mode 6 7 8
Freq [cm^-1] 890.2187 1242.2961 1290.2312
Reduced mass [au] 5.1552 1.0817 1.0566
Force const [Dyne/A] 2.4071 0.9836 1.0364
Char temp [K] 1280.8260 1787.3870 1856.3548
Normal mode x y z x y z x y z
C -0.19681 0.01743 0.04953 0.00334 0.04835 0.04166 0.03475 0.02738 -0.01209
F 0.13369 -0.01073 -0.02911 -0.01308 -0.02024 -0.01970 0.01700 -0.01944 0.00601
H -0.07363 0.02840 -0.00946 -0.39467 -0.01067 0.02670 -0.84240 -0.11715 0.09135
H -0.04312 0.02325 -0.01878 0.78493 -0.18960 -0.26094 -0.38233 0.07657 0.18889
H -0.32686 0.06551 0.02906 -0.06390 -0.15710 -0.10150 0.07328 -0.08739 -0.06467
O 0.02245 -0.01000 -0.00589 -0.01103 0.01072 0.01270 0.02169 0.01132 -0.00982
H -0.08934 0.03629 0.05167 0.05544 -0.00688 0.00958 0.07293 -0.01133 -0.02904
Mode 9 10 11
Freq [cm^-1] 1498.7370 1585.2045 3194.7998
Reduced mass [au] 1.1225 1.0784 1.0211
Force const [Dyne/A] 1.4856 1.5966 6.1403
Char temp [K] 2156.3482 2280.7556 4596.6041
Normal mode x y z x y z x y z
C 0.01397 -0.05640 0.07589 0.01042 0.06125 0.04567 0.03304 -0.00547 -0.00487
F -0.00691 0.00467 -0.00509 -0.00043 0.00034 -0.00035 0.00393 -0.00076 -0.00005
H -0.11405 -0.01249 -0.63257 0.14847 0.12066 -0.39741 -0.08733 0.57120 0.04967
H 0.02867 0.62173 -0.25947 -0.14160 -0.22159 0.26185 -0.24229 -0.21997 -0.45556
H 0.05404 -0.03469 0.09747 -0.13723 -0.63405 -0.40843 -0.12330 -0.28620 0.47884
O -0.00009 0.00060 -0.00084 0.00137 0.00119 -0.00154 -0.00175 -0.00054 -0.00157
H -0.00322 -0.00059 0.00026 -0.00733 -0.01962 0.03118 0.01316 0.02295 0.01069
Mode 12 13 14
Freq [cm^-1] 3276.7075 3393.2213 3410.0477
Reduced mass [au] 1.0697 1.1124 1.1148
Force const [Dyne/A] 6.7669 7.5462 7.6376
Char temp [K] 4714.4510 4882.0884 4906.2978
Normal mode x y z x y z x y z
C 0.00570 0.00611 0.00073 0.00030 -0.07030 0.06000 0.01672 0.05969 0.06982
F -0.00187 0.00014 0.00067 0.00038 0.00007 -0.00053 0.00033 0.00048 0.00013
H 0.00764 -0.06354 -0.00701 -0.08602 0.54248 0.04915 0.07681 -0.51536 -0.04505
H -0.01689 -0.00501 -0.01030 -0.08314 -0.07135 -0.14617 -0.32696 -0.29571 -0.61277
H 0.01395 -0.00306 -0.00493 0.16556 0.37064 -0.62222 0.04545 0.10393 -0.17086
O -0.02374 -0.04953 -0.02806 -0.00112 -0.00206 0.00008 0.00130 0.00192 0.00119
H 0.33929 0.78218 0.44635 0.01060 0.02682 0.01339 -0.02113 -0.04313 -0.02413
C 0.546814628839 0.109998777279 -0.085543963331
F 2.051053351342 -0.010846123021 -0.366170590968
H 0.415605465513 1.205445997160 0.003996530452
H 0.104636895203 -0.328417296548 -0.999394172348
H 0.225298223759 -0.417348433912 0.830930953932
O -1.476869479291 0.151854016198 0.550964704504
H -1.866539085364 -0.710686937156 0.065216537760
7
OptimizedStructure
C 0.546814628839 0.109998777279 -0.085543963331
F 2.051053351342 -0.010846123021 -0.366170590968
H 0.415605465513 1.205445997160 0.003996530452
H 0.104636895203 -0.328417296548 -0.999394172348
H 0.225298223759 -0.417348433912 0.830930953932
O -1.476869479291 0.151854016198 0.550964704504
H -1.866539085364 -0.710686937156 0.065216537760注意：DFTの数値積分に使用するグリッドの細かさに関する条件を変えてしまっているので、本文の結果とそのまま比較は不可能である。本文中のjsonファイルの条件を補足と合わせる必要がある。
【計算化学】自作pythonライブラリで遷移状態構造を求めてみる(PySCF使用, フルオロメタンと水酸化物イオンが関わるSN2反応, 基底関数の選択ミスによる失敗例)
https://ss0832.github.io/posts/20260127_mop_usage_sn2_3/ Related Posts
【計算化学】自作pythonライブラリで遷移状態構造を求めてみる(Pd錯体からの還元的脱離によるトルエン生成反応, NNP(UMA)使用)
2026-02-21
【計算化学】自作pythonライブラリで遷移状態構造を求めてみる(アゾジカルボン酸エステルによるN-保護ジアルキルピロール誘導体のC(sp3)–H官能基化(Triphenylborane触媒), NNP(UMA)使用)
2026-02-21
【計算化学】自作pythonライブラリで遷移状態構造を求めてみる(アゾジカルボン酸エステルによるN-保護ジアルキルピロール誘導体のC(sp3)–H官能基化(Tris(biphenyl-4-yl)borane触媒), NNP(UMA)使用)
2026-02-20
【計算化学】自作pythonライブラリで遷移状態構造を求めてみる(アゾジカルボン酸エステルによるN-保護ジアルキルピロール誘導体のC(sp3)–H官能基化(Tris(pentachlorophenyl)borane触媒), NNP(UMA)使用)
2026-02-20
【計算化学】自作pythonライブラリで遷移状態構造を求めてみる(アゾジカルボン酸エステルによるN-保護ジアルキルピロール誘導体のC(sp3)–H官能基化(Tris(pentafluorophenyl)borane触媒), NNP(UMA)使用)
2026-02-20
【計算化学】自作pythonライブラリで遷移状態構造を求めてみる(アゾジカルボン酸エステルによるN-保護ジアルキルピロール誘導体のC(sp3)–H官能基化(Tris(4-tert-butylphenyl)borane触媒), NNP(UMA)使用)
2026-02-20