【計算化学】自作pythonライブラリで遷移状態構造を求めてみる(PySCF使用, フルオロメタンと水酸化物イオンが関わるSN2反応)

4022 words

20 minutes

【計算化学】自作pythonライブラリで遷移状態構造を求めてみる(PySCF使用, フルオロメタンと水酸化物イオンが関わるSN2反応)

2026-01-27

Computational Chemistry

MultiOptPy

python

最終更新：2026-01-27

概要#

本記事では、自作ライブラリ(MultiOptPy)で、フルオロメタンと水酸化物イオンが関わるSN2反応の素過程の遷移状態構造を算出してみる。計算レベルは、ωB97X-D/6-31Gとした。

STO-3Gから6-31Gに高価な基底関数に変えたことで、モデルの記述性能が向上し、SN2反応の遷移状態構造として妥当なWalden反転が見られたことを示す。

MultiOptPyは電子状態計算ソフトウェアを用いた分子構造最適化手法の勉強を目的として作成したpythonライブラリである。

MultiOptPyのレポジトリ:https://github.com/ss0832/MultiOptPy

使用した自作ライブラリMultiOptPyのバージョン#

v1.20.5

Video 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_4.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.347512

初期経路を求めるための初期構造

2. 遷移状態構造最適化#

　run_autots.pyを適切に使用することで、自動的に遷移状態構造が得られる。以下にその手順を説明していく。

初期構造をMultiOptPy-1.20.5ディレクトリ内にSN2_4.xyzとして保存する。その後、同じディレクトリ内で、config_SN2_4.jsonを作成し、以下のように記述する。

config_SN2_4.json

{
  "work_dir": "SN2_4",
  "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": "6-31G",
	"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": "6-31G",
	"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": "6-31G",
	"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_4.xyz -cfg config_SN2_4.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"]:力の定数c a.u.で、平行距離dÅの調和ポテンシャルをa番とb番の原子ペアにかける。初期経路生成時に開裂を防ぎたい結合を保持するとき等に使用する。
"usextb": "GFN2-xTB"はGFN2-xTBを使用することを指定する。QSM法(NEB法の類似法)にて初期構造を手早く生成するときに使用可能である。

Step1が正常終了していれば作成されたwork_dirディレクトリ中に、SN2_4_step1_traj.xyzが存在する。必要に応じて確認し、目的に沿った初期経路が得られているか確認する。もし想定とは異なる場合は、プロセスをkillして再度設定を見直してやり直す。

SN2_4_step1_traj.xyzは構造最適化の過程をAvogadro（公式ページ:https://avogadro.cc/ ）等で可視化して確認できるようにしている。このSN2_4_step1_traj.xyzはStep2のNEB計算に使用している。

※SN2_4_step1_traj.xyzをアニメーションとして表示したい場合は、[https://github.com/ss0832/molecule_movie] を使うと良い。

Step 2#

Step2では、NEB法を用いることで、先ほど得られたSN2_4_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を確認し、緩和後のパスのエネルギー極大値を示す構造を確認する。

経路の緩和後の各ノードのエネルギー一覧(単位) (energy_plot.csvに保存されている。)

NEB計算の結果の可視化

※bias_force_rms.csvにて、各Iterationごとのすべてのノードの勾配のRMS値を確認できる。

経路緩和の結果、以下の構造がstep3の初期構造として自動的に用いられた。“work_dir”内のSN2_4_step3_TS_Opt_Inputs内に保存されたSN2_4_ts_guess_X.xyzにて確認が可能である。ts_guessの番号が小さい順にエネルギー値が高い構造を示すようになっている。

※こちら[https://ss0832.github.io/molecule_viewer/] を使うことでも可視化は可能である。

SN2_4_ts_guess_1.xyz

7
-1 0
C       0.336243477453      0.055828440371      0.020836910971
F       2.214003755989      0.158497384462     -0.164062978899
H       0.296243614799      1.116593555996      0.117412586145
H       0.270033265084     -0.378188131565     -0.961697666562
H       0.473301182622     -0.563486243061      0.877335058541
O      -1.529574516698      0.073485798528      0.370745209888
H      -2.060250779249     -0.462730804731     -0.260569120084

NEB法により緩和した経路から得られた遷移状態構造を求めるための初期構造 (No.1)

Step 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_4_ts_final_X.xyzとして保存されている。)

SN2_4_ts_final_1.xyz

7
OptimizedStructure
C      0.388912739426      0.107810823429     -0.005031694316
F      2.146664331901      0.030159839170     -0.139419416520
H      0.356744940624      1.163286910272      0.206088451094
H      0.234335266143     -0.227049673602     -1.017499986129
H      0.358472753062     -0.589920264410      0.815556238916
O     -1.602006345203      0.222419836586      0.152460880191
H     -1.883123685953     -0.706707471445     -0.012154473236

遷移状態構造 (No.1)

16回の反復計算により、停留点に収束した構造が得られた。"frequency_analysis": trueオプションにより生成されたnormal_modes.txtやvibration_animationディレクトリ内の振動モードのアニメーションを確認した。

以下に"frequency_analysis": trueオプションで生成されたnormal_modes.txtの一部を示す。

Mode                                 0                   1                   2
Freq [cm^-1]                     -449.1335            184.5076            295.5290
Reduced mass [au]                  10.0413              1.0569              3.1229
Force const [Dyne/A]               -1.1934              0.0212              0.1607
Char temp [K]                       0.0000            265.4652            425.2003
Normal mode                   x         y         z            x         y         z            x         y         z     
       C                -0.24971    0.00975    0.01883    0.00074   -0.00214    0.01133   -0.00602   -0.19630   -0.03855
       F                 0.06087   -0.00296   -0.00470   -0.00046    0.00145   -0.00756   -0.00004    0.08707    0.01728
       H                 0.04501    0.01956    0.00084   -0.01668    0.04912   -0.24755   -0.05980   -0.20428   -0.03663
       H                 0.08753   -0.01429   -0.02431    0.02015   -0.24385    0.08558    0.01451   -0.22477   -0.03608
       H                 0.08952   -0.02022    0.00544   -0.00260    0.19281    0.17358    0.01320   -0.22101   -0.05414
       O                 0.09702   -0.00366   -0.00729    0.00354   -0.01051    0.05297    0.02506    0.07453    0.01313
       H                 0.06394    0.01270   -0.00174   -0.05722    0.16687   -0.84466   -0.29320    0.16324    0.05175
       
(...snip...)

その結果、虚振動が１つであることが確認できた。つまりこの構造は遷移状態構造である。STO-3Gを基底関数として用いた場合と異なりこちらはWalden反転が見られる遷移状態構造であることがわかる。

次に、vibration_animation内の虚振動を示す分子振動が示されたxyzファイル(mode_1_XXXi_wave_number.xyz)をAvogadroで確認すると、求められた遷移状態構造の中に、想定される反応系と生成系をつなぐ方向に振動している構造が存在することを確認できた。

終わりに#

　　自作ライブラリを用いて、フルオロメタンと水酸化物イオンが関わるSN2反応のある1つの遷移状態構造を算出する手順を説明した。

参考#

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