最終更新：2026-02-11
概要
本記事では、自作ライブラリ(MultiOptPy)で、1-ブロモペンタンへのリチウムアセチリドのSN2反応によるアルキニル化の素過程の遷移状態構造を算出してみる。計算レベルは、Meta社のFAIR Chemistryが開発したニューラルネットワークポテンシャル(NNP)であるUMA(Meta’s Universal Model for Atoms)とした。
MultiOptPyは電子状態計算ソフトウェアを用いた分子構造最適化手法の勉強を目的として作成したpythonライブラリである。
MultiOptPyのレポジトリ:https://github.com/ss0832/MultiOptPy
今回使用したニューラルネットワークポテンシャルについて：
- 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). (プレプリント)
使用した自作ライブラリMultiOptPyのバージョン
v1.20.8Video Demo
https://www.youtube.com/watch?v=AE61iY2HZ8Y
環境
Windows 11※Windows 11環境下でAnaconda PowerShell Promptを使用した。
Source codeのダウンロード(Unixコマンド)
wget https://github.com/ss0832/MultiOptPy/archive/refs/tags/v1.20.8.zip
unzip v1.20.8.zip
cd MultiOptPy-1.20.8https://github.com/ss0832/MultiOptPy/releases/tag/v1.20.8 にアクセスしてzipファイルをダウンロードする。Unixコマンドの場合とはディレクトリ名が異なるので都度読み替えていただけると良い。
移動先のディレクトリでrequirements.txtを参照することで、本ソースコードで必要なモジュールを把握することが出来る。導入方法は各自の状況に合わせて適宜LLMとの対話などで調べると良い。
次に述べる環境構築手順を使用する場合は、環境構築が終わった後、pip install -r requirements.txtで本自作モジュールが動作させるために最低限必要なモジュールを導入することが可能である。
環境構築手順
今回は、Windows 11のPower Shellを使用した。初めに、NNPを使用できる環境が整ったAnaconda PowerShell Promptを用意する手順を説明する。
1, https://repo.anaconda.com/archive/ より、Anaconda3-2025.06-1-Windows-x86_64.exeでAnacondaをインストールする。
2, 検索機能を使い、スタートからAnaconda PowerShell Promptを開く。
3, 以下のコマンドを実行し、仮想環境を作成する。
conda create -n (任意の仮想環境名) python=3.12.74, 先ほど作成した仮想環境をconda activate (仮想環境名)で起動させる。
5, 以下のコマンドを実行し、必要なライブラリを導入する。
pip install ase==3.26.0 fairchem-core==2.7.1 torch==2.6.0- fairchem-coreは、FAIR Chemistryが管理しているNNPを動作させるために必要なライブラリである。
- aseはNNPに電子エネルギーを算出したい分子構造を渡すために必要なインターフェイスの役割を果たすために必要なライブラリである。
- torchはPyTorchライブラリを指す。これはニューラルネットワークなどの機械学習を行ったり、学習結果を扱ったりするために必須なライブラリである。
これで、Anaconda PowerShell Promptから仮想環境を立ち上げることで、NNPを使用する準備が整えることが出来る。
次に、NNPを使用するために必要なModelの情報が保存されている.ptファイルのダウンロードおよびNNPの自作ライブラリへの導入方法について説明する。
UMAを使用可能にするための手順
1, 以下のサイトにアクセスして、uma-s-1p1.ptをダウンロードする。（使用許諾が下りていれば可能である。）
https://huggingface.co/facebook/UMA2, ダウンロード後、MultiOptPy-1.20.8ディレクトリ内に存在するsoftware_path.confに対して、uma-s-1p1.ptの絶対パスを用いて以下を追記する。
uma-s-1p1::(uma-s-1p1.ptの絶対パス)これで、MultiOptPy-1.20.8がNNPuma-s-1p1を使用できるようになる。
※Linuxの場合
すでにAnacondaの導入が終わっている前提で、環境の構築手段について述べる。以下を実行すると可能である。
## 1. Download and install MultiOptPy:
wget https://github.com/ss0832/MultiOptPy/archive/refs/tags/v1.20.8.zip
unzip v1.20.8.zip
cd MultiOptPy-1.20.8
## 2. Create and activate a conda environment:
conda env create -f environment.yml
conda activate test_mop
## 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
#### Installation via pip (Linux)
conda create -n <env-name> python=3.12 pip
conda activate <env-name>
pip install git+https://github.com/ss0832/MultiOptPy.git@v1.20.8
# or pip install multioptpy
wget https://github.com/ss0832/MultiOptPy/archive/refs/tags/v1.20.8.zip
unzip v1.20.8.zip
cd MultiOptPy-1.20.8
あとは、「UMAを使用可能にするための手順」と同じように進めれば問題なく導入可能である。
使用するNNPに関する具体的な説明
今回使用するNNPについて具体的に説明する。
- UMAのModel Checkpointは
uma-s-1p1を使用した。 - 小分子系のトレーニングセットである
Omol25(omol)を使用して学習したニューラルネットワークポテンシャルを使用する。
※自作ライブラリでの具体的な使用の仕方に関しては、ase_calculation_tools.py を参照すると良い。omol以外のモデルを使用したい場合は、現バージョンでは、multioptpy/Calculator/ase_tools/firechem.py内の、self.task_nameを編集することで対応可能である。
手順
1. 初期構造の準備
モデル反応系として、以下の構造を用意した。今回はファイルの名前をsn2_alkyne_n_PeBr.xyzとした。
初期構造は以下のものを使用した。
21
OptimizedStructure
C 2.007074729955 -0.421348750868 -1.254989266941
H 2.576829348174 -0.108762354535 -0.377438320312
H 2.372004587418 0.095424293179 -2.141679070556
C -0.198676639097 3.352987944251 -0.105903475087
H -1.243555115783 3.523029476170 -0.055357464165
C 0.977226434178 3.078407768215 -0.191086748779
Li 2.776463074377 3.446586677022 -0.068915041332
Br 2.438310662611 -2.309432204379 -1.527775892114
C 0.525422610568 -0.208115512535 -1.043970737448
H 0.375025084441 0.870476776620 -0.923921726390
H -0.031097677595 -0.541945854097 -1.923056999654
C 0.012903666113 -0.928622494646 0.197141720301
H 0.160080799039 -2.006358937264 0.088090851587
H 0.588236929193 -0.595386795849 1.065459902929
C -1.464343083134 -0.649060470289 0.444822236252
H -2.042026618809 -0.968094522970 -0.425549555518
H -1.604940616355 0.426942247068 0.562918939379
C -1.964754943317 -1.376830805860 1.686062641155
H -2.974945890218 -1.063575430608 1.938909961849
H -1.968421232978 -2.452522321058 1.522137318635
H -1.316816108779 -1.163798727568 2.5341007262072. 遷移状態構造最適化
run_autots.pyを適切に使用することで、自動的に遷移状態構造が得られる。以下にその手順を説明していく。
初期構造をMultiOptPy-1.20.8ディレクトリ内にsn2_alkyne_n_PeBr.xyzとして保存する。その後、同じディレクトリ内で、config_sn2_alkyne_n_PeBr.jsonを作成し、以下のように記述する。
config_sn2_alkyne_n_PeBr.json
{
"work_dir": "sn2_alkyne_n_PeBr",
"top_n_candidates": 3,
"step1_settings": {
"othersoft": "uma-s-1p1",
"opt_method": ["rsirfo_block_fsb"],
"use_model_hessian": "fischerd3old",
"spin_multiplicity": 1,
"electronic_charge": 0,
"manual_AFIR": ["310", "1", "6"]
},
"step2_settings": {
"othersoft": "uma-s-1p1",
"NSTEP": 20,
"use_model_hessian": "fischerd3old",
"save_pict": true,
"node_distance": 0.40,
"align_distances_energy_predicted": 2,
"spin_multiplicity": 1,
"electronic_charge": 0
},
"step3_settings": {
"othersoft": "uma-s-1p1",
"opt_method": ["rsirfo_block_bofill"],
"calc_exact_hess": 5,
"tight_convergence_criteria": true,
"max_trust_radius": 0.2,
"frequency_analysis": true,
"spin_multiplicity": 1,
"electronic_charge": 0
},
"step4_settings": {
"othersoft": "uma-s-1p1",
"opt_method": ["rsirfo_block_bofill"],
"spin_multiplicity": 1,
"electronic_charge": 0,
"calc_exact_hess": 10,
"tight_convergence_criteria": true,
"frequency_analysis": true,
"intrinsic_reaction_coordinates": ["0.5", "200", "lqa"],
"step4b_opt_method": ["rsirfo_block_fsb"]
}
}その後、以下のコマンドを実行する。
python run_autots.py sn2_alkyne_n_PeBr.xyz -cfg config_sn2_alkyne_n_PeBr.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を生成する機能を呼び出すオプションである。デフォルトではこの機能は使用されない。
※オプションの説明はMultiOptPy-v1.20.8/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番の原子ペアにかける。初期経路生成時に開裂を防ぎたい結合を保持するとき等に使用する。
Step1が正常終了していれば作成されたwork_dirディレクトリ中に、sn2_alkyne_n_PeBr_step1_traj.xyzが存在する。必要に応じて確認し、目的に沿った初期経路が得られているか確認する。もし想定とは異なる場合は、プロセスをkillして再度設定を見直してやり直す。
sn2_alkyne_n_PeBr_step1_traj.xyzは構造最適化の過程をAvogadro（公式ページ:https://avogadro.cc/ ）等で可視化して確認できるようにしている。このsn2_alkyne_n_PeBr_step1_traj.xyzはStep2のNEB計算に使用している。
※sn2_alkyne_n_PeBr_step1_traj.xyzをアニメーションとして表示したい場合は、[https://github.com/ss0832/molecule_movie] を使うと良い。
Step 2
Step2では、NEB法を用いることで、先ほど得られたsn2_alkyne_n_PeBr_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.8/"work_dir"と同じディレクトリ内に、NEBという名前を含むディレクトリが生成されている。
そのディレクトリ内のenergy_plot.csvを確認し、緩和後のパスのエネルギー極大値を示す構造を確認する。
経路の緩和後の各ノードのエネルギー一覧(単位
※bias_force_rms.csvにて、各Iterationごとのすべてのノードの勾配のRMS値を確認できる。
経路緩和の結果、以下の構造がstep3の初期構造として自動的に用いられた。“work_dir”内のsn2_alkyne_n_PeBr_step3_TS_Opt_Inputs内に保存されたsn2_alkyne_n_PeBr_ts_guess_X.xyzにて確認が可能である。ts_guessの番号が小さい順にエネルギー値が高い構造を示すようになっている。
※こちら[https://ss0832.github.io/molecule_viewer/] を使うことでも可視化は可能である。
sn2_alkyne_n_PeBr_ts_guess_1.xyz
21
0 1
C 1.717803758003 0.363282044879 -1.078475327069
H 2.351675276755 0.170176765556 -0.234574881869
H 2.190860620880 0.659513120003 -1.995573473482
C 0.407575823079 3.025363509218 -0.137714921978
H -0.474398781547 3.580537850356 0.100560714223
C 1.295878493826 2.252214155821 -0.461156667381
Li 3.038404840188 3.221853329122 -0.229093231050
Br 2.898569810329 -1.902073440219 -1.689140758547
C 0.301129276146 -0.160218470600 -1.098612633770
H -0.387576990753 0.655507139863 -1.317385236007
H 0.227172282984 -0.856991478366 -1.934284717016
C -0.107624347167 -0.854854347444 0.201304253693
H 0.171134205002 -1.910010795838 0.148100589885
H 0.457283386923 -0.431131960649 1.038361735093
C -1.595200283769 -0.705799834631 0.494851490986
H -2.179092614932 -1.049514432640 -0.365123974289
H -1.820928556056 0.359884108731 0.612094631379
C -2.038620169678 -1.457727381787 1.743092731593
H -3.088132009828 -1.266640199350 1.971318821629
H -1.916746636566 -2.535080794201 1.618872810260
H -1.449167383820 -1.158288887824 2.612578043717Step 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.8/"work_dir"ディレクトリ内に生成された新規ディレクトリ内のsn2_alkyne_n_PeBr_ts_final_X.xyzとして保存されている。)
sn2_alkyne_n_PeBr_ts_final_1.xyz
21
OptimizedStructure
C 1.321442282241 0.580523839678 -1.028745114234
H 1.968819530664 0.525625205852 -0.178405951038
H 1.617866441426 1.090723538386 -1.921586373672
C 1.335311311870 3.086926056077 -0.122815701965
H 0.407109529624 3.449197802349 0.257656634818
C 2.470330689201 2.806280901033 -0.508422207934
Li 3.967249582461 1.589004561000 -1.208642311653
Br 3.521226634947 -0.723715069100 -1.854130716598
C -0.014300584461 -0.020377753266 -0.963987101711
H -0.731969461690 0.701504257927 -1.371551589559
H -0.004677924757 -0.830167547815 -1.709370226206
C -0.428410797336 -0.535390184940 0.405427848936
H 0.313788559879 -1.261693132898 0.750547633175
H -0.407675603387 0.294186061078 1.118890392144
C -1.810060470696 -1.176333780418 0.403700398377
H -1.822934176672 -2.003063140905 -0.313069561835
H -2.545971040160 -0.448766816489 0.047554473001
C -2.213746486382 -1.683757367161 1.781478860962
H -3.204487028033 -2.137678056035 1.765153708926
H -1.507944471277 -2.434105787743 2.142534822479
H -2.230966517461 -0.868923586609 2.50778208359076回の反復計算により、停留点に収束した構造が得られた。"frequency_analysis": trueオプションにより生成されたnormal_modes.txtやvibration_animationディレクトリ内の振動モードのアニメーションを確認した。
以下に"frequency_analysis": trueオプションで生成されたnormal_modes.txtの一部を示す。
Mode 0 1 2
Freq [cm^-1] -300.9333 31.3525 59.9584
Reduced mass [au] 6.5617 3.4929 1.6814
Force const [Dyne/A] -0.3501 0.0020 0.0036
Char temp [K] 0.0000 45.1093 86.2667
Normal mode x y z x y z x y z
C -0.06279 0.13870 0.05549 0.00743 -0.06008 0.00106 0.00355 -0.02368 -0.00998
H 0.05099 -0.13771 -0.05015 0.00416 -0.05425 0.00416 0.00136 -0.01533 -0.00789
H 0.02242 -0.07427 -0.03727 -0.00350 -0.03124 0.01410 -0.00067 -0.02182 -0.01039
C -0.05998 -0.09922 -0.02086 -0.04155 -0.06083 0.05471 -0.04401 -0.02732 -0.00401
H -0.05314 -0.08412 -0.01719 -0.04719 -0.08307 0.06210 -0.05038 -0.04622 -0.00141
C -0.04128 -0.06618 -0.01012 -0.03816 -0.03642 0.04771 -0.04019 -0.00634 -0.00752
Li 0.01446 -0.05247 -0.02090 -0.01498 0.01157 0.02135 -0.02022 0.01612 -0.00406
Br 0.03593 -0.01012 -0.01160 0.03156 0.01581 -0.03055 0.01834 0.00452 0.00126
C -0.04784 0.10970 0.04212 0.01805 -0.08432 -0.00797 0.00032 -0.01800 -0.00643
H -0.11279 0.02795 0.02107 0.01375 -0.13621 -0.08975 0.00483 0.03529 0.07605
H 0.07964 0.11131 0.03394 0.07293 -0.14249 0.05505 -0.05856 0.03868 -0.06636
C -0.02082 0.03243 0.01771 -0.02270 0.00637 0.01468 0.03893 -0.09120 -0.02258
H 0.00286 0.03860 -0.02067 -0.02992 0.03566 0.09177 -0.03714 -0.22999 -0.15169
H -0.03872 0.00052 0.05579 -0.05163 0.05726 -0.04380 0.21656 -0.16075 0.05361
C -0.00514 0.00137 0.00382 -0.01982 0.00012 0.01412 -0.05676 0.11520 0.03220
H 0.01390 0.00442 0.00014 0.01270 -0.06365 0.08716 -0.25300 0.20387 -0.06662
H -0.01694 -0.01062 0.00353 -0.00741 -0.03400 -0.08116 0.02023 0.27109 0.19147
C -0.00218 -0.00706 0.00101 -0.07611 0.12323 0.04300 -0.00098 0.01026 0.00996
H 0.00219 -0.01647 -0.00072 -0.07276 0.11597 0.04158 -0.08026 0.18163 0.05539
H 0.00444 -0.00063 0.00123 -0.08750 0.16037 0.14251 -0.09175 -0.15565 -0.15761
H -0.01085 -0.00797 0.00201 -0.11243 0.18965 -0.03240 0.21345 -0.07825 0.11442
(...snip...)その結果、虚振動が１つであることが確認できた。つまりこの構造は遷移状態構造である。
IRC計算により、求められた遷移状態構造に対応する生成系と始原系に対応する安定構造に近い構造を求めた。 コマンド：
python .\optmain.py .\sn2_alkyne_n_PeBr_ts_final_1.xyz -os uma-s-1p1 -opt rsprfo_block_bofill -fc 10 -order 1 -tcc -freq -irc 0.1 200 lqa これにより求められたIRC経路の末端構造（生成系と始原系に対応する安定構造に近い構造）を以下に示す。
sn2_alkyne_n_PeBr_ts_final_1_0_irc_endpoint_1.xyz
21
Terminal structure
C -1.0533916447 1.4347395420 -0.5253247909
H -0.2240921280 0.7264360306 -0.6764652729
H -1.1636151374 1.9547789720 -1.4763392362
C -0.6966959876 2.5380838813 0.5338468676
H -1.3742604174 3.3079552081 0.8605240367
C 0.4889809007 2.1161097129 0.7071283254
Li 2.2368130936 1.0467623025 0.3013119500
Br 2.1302494736 -0.7100820090 -1.0873036012
C -2.3073922862 0.6617691290 -0.1321136910
H -3.1570329988 1.3509511866 -0.0902853058
H -2.5085707431 -0.0357642848 -0.9500222554
C -2.2206147937 -0.1118926733 1.1769904462
H -1.3403960401 -0.7628446474 1.1618247008
H -2.0773737795 0.5826740264 2.0110449078
C -3.4793022714 -0.9406051267 1.4147972202
H -3.5641381569 -1.6922090490 0.6246870289
H -4.3566573008 -0.2931907419 1.3147338968
C -3.5047929759 -1.6199555073 2.7771083784
H -4.4082840692 -2.2163281445 2.9070441274
H -2.6464004203 -2.2827020608 2.9012730191
H -3.4728565830 -0.8824411801 3.5817454745sn2_alkyne_n_PeBr_ts_final_1_0_irc_endpoint_2.xyz
21
Terminal structure
C -0.3385212601 0.3321473809 -0.9645089973
H -0.2654132199 1.2270019495 -0.3487062847
H -0.3486696692 0.6075863019 -2.0136088452
C 0.2932654820 3.8242950575 0.7842852274
H -0.4636044571 4.5053643897 1.0959575462
C 1.1810555873 3.0639569747 0.4347270411
Li 2.4326338614 1.6800040455 -0.0810564135
Br 1.4367190716 -0.5441332427 -0.7626275685
C -1.4897862634 -0.5751595801 -0.5987611716
H -2.3680258586 -0.1059659293 -1.0587849274
H -1.3807083979 -1.5478972455 -1.0855701880
C -1.7688219546 -0.7419260533 0.8908613081
H -0.9675907954 -1.3143561117 1.3680647858
H -1.7732350321 0.2443240501 1.3658818954
C -3.1143390495 -1.4264012169 1.1242355065
H -3.0915751817 -2.4278021164 0.6832012603
H -3.8863825711 -0.8725355460 0.5808334927
C -3.5042261929 -1.5189554951 2.5934630731
H -4.4745836050 -2.0010208944 2.7158493170
H -2.7711978982 -2.0950024426 3.1614088035
H -3.5649666069 -0.5263050750 3.0434411307
構造最適化の結果得られた生成物側の幾何構造に基づくと、sn2_alkyne_n_PeBr_ts_final_1.xyz は 反応によるアルキニル化の遷移状態ではなく、アルキリデンカルベン様化学種の生成に対応する構造であると推測される。
最適化のアルゴリズムを変えて、遷移状態構造を求め直すこととした。
python .\optmain.py .\sn2_alkyne_n_PeBr_ts_guess_1.xyz -os uma-s-1p1 -opt rsprfo_block_bofill -fc 5 -order 1 -tcc -freq -tr 0.2遷移状態の探索アルゴリズムとして、ポテンシャルエネルギー曲面を変換した像関数に対する標準的なRFO法の適用に代わり、ヘシアンの固有ベクトル空間を最大化モードと最小化モードに分割して処理する分割RFO（P-RFO）法を採用した。(-opt rsprfo_block_bofillで使用可能である。)
312回の反復計算により、以下の構造が得られた。
sn2_alkyne_n_PeBr_ts_guess_1_optimized.xyz
21
OptimizedStructure
C 1.142595177162 0.724326943664 -0.868681288601
H 1.493002371258 1.035087913107 0.094032648443
H 1.553504695107 1.056458629563 -1.799743014349
C 2.179718739664 3.856491313880 -0.596173072845
H 2.843872282683 4.673809896051 -0.417380187429
C 1.216485570493 3.108437150765 -0.786130055111
Li 3.196313126064 1.871223372993 -0.748948137627
Br 3.263344188204 -0.628079059832 -0.791170689696
C -0.118309724828 -0.040487896099 -0.911880154681
H -0.885075715166 0.685794496090 -1.211690089922
H -0.063648814336 -0.762256545690 -1.731740599762
C -0.489574267723 -0.707453176719 0.405254879352
H 0.315722185000 -1.388047244541 0.696478105238
H -0.556383532173 0.056951707345 1.186156535626
C -1.805165603798 -1.470312503632 0.323638344276
H -1.734638708270 -2.225943135289 -0.464627779287
H -2.602999710967 -0.784379257475 0.022902925922
C -2.169161890317 -2.137368769595 1.643312410212
H -3.112280788019 -2.678882395023 1.570203377048
H -1.397926963717 -2.847903491393 1.946467123571
H -2.269392616321 -1.397467948169 2.439718719623
以下に"-freqオプションで生成されたnormal_modes.txtの一部を示す。
Mode 0 1 2
Freq [cm^-1] -394.2816 25.2815 48.2404
Reduced mass [au] 6.2493 3.7154 3.7847
Force const [Dyne/A] -0.5724 0.0014 0.0052
Char temp [K] 0.0000 36.3745 69.4071
Normal mode x y z x y z x y z
C -0.09887 0.17365 0.00984 0.00269 -0.05016 -0.01765 -0.00774 0.01484 -0.07819
H 0.09583 -0.12162 0.02911 -0.00064 -0.04677 -0.01777 -0.01613 0.04161 -0.08405
H 0.07722 -0.09030 -0.00194 -0.00574 -0.03849 -0.01729 0.01259 -0.02783 -0.08444
C -0.05610 -0.05016 -0.01402 -0.06791 -0.02113 0.10071 -0.02404 -0.01631 0.14314
H -0.08625 -0.01918 -0.03859 -0.09638 -0.01443 0.17560 -0.04992 -0.03974 0.34626
C -0.01362 -0.10712 0.01389 -0.03425 -0.04134 0.01185 0.00915 0.00657 -0.11273
Li 0.04547 -0.06583 0.00136 -0.02014 0.00610 -0.00129 -0.00575 -0.00477 -0.04511
Br 0.02717 -0.01062 0.00060 0.04294 0.00692 -0.02425 -0.02289 -0.01213 0.01871
C -0.04921 0.07677 0.00068 0.01363 -0.06759 -0.01806 -0.01103 0.01903 -0.05480
H -0.14508 -0.04344 -0.02826 0.01678 -0.09820 -0.10025 -0.01805 0.01882 -0.03711
H 0.07014 0.07136 0.01098 0.05817 -0.11823 0.02871 -0.03126 0.01363 -0.05094
C -0.00292 0.02424 -0.00472 -0.03281 0.01259 0.00848 0.02697 0.02834 -0.03862
H 0.01909 0.03763 -0.03723 -0.04241 0.03196 0.08118 0.04388 0.04271 -0.05212
H -0.01043 0.00247 0.01627 -0.06270 0.06194 -0.04242 0.03835 0.03636 -0.04534
C 0.00546 0.00374 -0.00246 -0.02908 0.00599 0.00955 0.03621 0.00879 0.00817
H 0.01168 0.00439 -0.00255 0.00344 -0.05193 0.06799 0.02159 0.00467 0.01083
H -0.00193 -0.00412 -0.00176 -0.01610 -0.01561 -0.07419 0.01455 -0.00619 0.03159
C 0.00607 0.00079 -0.00269 -0.08543 0.10419 0.04363 0.09514 0.01120 0.02582
H 0.00756 -0.00164 -0.00582 -0.08183 0.09797 0.04339 0.10128 -0.00469 0.06455
H 0.00676 0.00263 -0.00082 -0.09809 0.12739 0.13034 0.11823 0.02596 0.00159
H 0.00176 0.00014 -0.00248 -0.12046 0.16382 -0.01620 0.11287 0.01472 0.02482
(...snip...)その結果、虚振動が１つであることが確認できた。つまりこの構造(sn2_alkyne_n_PeBr_ts_guess_1_optimized.xyz)は遷移状態構造である。
IRC計算により、求められた遷移状態構造に対応する生成系と始原系に対応する安定構造に近い構造を求めた。
コマンド：
python .\optmain.py .\sn2_alkyne_n_PeBr_ts_guess_1_optimized.xyz -os uma-s-1p1 -irc 0.5 300 lqa -opt rsirfo_block_bofill -fc 10 -order 1 -tcc -freq -tr 0.2これにより求められたIRC経路の末端構造（生成系と始原系に対応する安定構造に近い構造）を以下に示す。
sn2_alkyne_n_PeBr_ts_guess_1_optimized_0_irc_endpoint_1.xyz
21
Terminal structure
C -1.0978490309 1.2700718175 -1.2273642271
H -0.3556194499 0.4955919161 -1.4484021050
H -1.2926770183 1.8036818381 -2.1600406741
C -0.0043115745 3.0039054882 0.4831356227
H 0.3526658639 3.7526217479 1.1521026190
C -0.4972009717 2.2223315454 -0.2895977483
Li 1.5412005686 1.2908371516 0.3382129144
Br 2.1613967525 -0.7085787243 -0.4019890423
C -2.3821950239 0.6271773790 -0.6903262806
H -3.1106844644 1.4074931686 -0.4527881465
H -2.8109385608 0.0322284502 -1.5005576758
C -2.1540622915 -0.2651077245 0.5227487100
H -1.3866370816 -1.0090396056 0.2853496573
H -1.7575332081 0.3312591673 1.3518290765
C -3.4315434503 -0.9645695256 0.9723779818
H -3.8337752693 -1.5449461535 0.1360769956
H -4.1888426417 -0.2122527132 1.2152125594
C -3.2108489778 -1.8802080141 2.1697217286
H -4.1353644694 -2.3732460351 2.4720419078
H -2.4778837549 -2.6548355545 1.9371365004
H -2.8356872393 -1.3183381479 3.0280185885
sn2_alkyne_n_PeBr_ts_guess_1_optimized_0_irc_endpoint_2.xyz
21
Terminal structure
C -0.3020032650 0.2991806662 -0.9928496257
H -0.3964637794 1.2442964933 -0.4628765094
H -0.0902088045 0.4923635470 -2.0388136999
C 0.4693019013 4.3841078956 -0.2251319763
H -0.1375363223 5.2567140018 -0.2953746487
C 1.1757423269 3.3924720765 -0.1384000173
Li 2.3124426631 1.8388130238 0.0213250916
Br 1.3911731498 -0.4825955683 -0.3084562556
C -1.4558535119 -0.6519690869 -0.7886848168
H -2.3080382923 -0.1942760788 -1.3047557342
H -1.2558120187 -1.5964279924 -1.3019853762
C -1.8299710688 -0.9074073702 0.6647263732
H -0.9883148838 -1.3714066263 1.1871061213
H -2.0151125986 0.0495776756 1.1637252596
C -3.0590026034 -1.7995273251 0.7965997631
H -2.8706069143 -2.7484481547 0.2850237745
H -3.9004202435 -1.3324521725 0.2758890666
C -3.4368924393 -2.0645024695 2.2479646910
H -4.3193851236 -2.7003375380 2.3216243643
H -2.6223411472 -2.5605066565 2.7793646268
H -3.6526355728 -1.1312438794 2.7715159211
構造最適化の結果得られた生成物側の幾何構造に基づくと、sn2_alkyne_n_PeBr_ts_guess_1_optimized.xyz は標準的な 反応によるアルキニル化の遷移状態であると考えられる。
終わりに
自作ライブラリで、UMAモデルのニューラルネットワークポテンシャル(NNP, uma-s-1p1)を用いて、1-ブロモペンタンへのリチウムアセチリドの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.
- The Journal of Physical Chemistry, 1985, 89, 1, 52–57
補足
1. 像関数（Image Potential）に基づくアプローチ
遷移状態（一次の鞍点）の探索を、数学的な極小点探索問題へと変換して処理する手法である。
- 原理: 元のポテンシャルエネルギー曲面（PES）において、特定の方向（想定される反応座標など）に沿った勾配の符号を反転させる操作を行い、仮想的な関数（像関数）を再定義する。
- 特徴: この数学的変換により、元のPES上において鞍点として存在する遷移状態が、像関数上では極小点として表現される構造的対応を持つ。したがって、構築された像関数に対しては、標準的な有理関数最適化（RFO: Rational Function Optimization）等の極小点探索アルゴリズムをそのまま適用することが可能となる。
2. 分割RFO法（Partitioned RFO; P-RFO法）
像関数のような仮想的なポテンシャル面を構築せず、元のPES上においてヘシアン行列の固有ベクトル空間を直接操作する手法である（Banerjee et al., 1985）。
- 原理: ヘシアン行列の固有ベクトル空間を2つのブロックに明示的に分割（Partition）する。遷移状態探索においては、系全体のエネルギーを最大化すべき方向（通常、1つの負の固有値に対応するモード）と、最小化すべき方向（残りの正の固有値に対応する全モード）に分離される。
- 特徴: 分割された各ブロックに対して独立した有理関数近似を適用し、それぞれ固有のシフトパラメータ（）の最適化を実行する。これにより、特定の1モードに沿ったエネルギーの最大化と、それに直交する部分空間でのエネルギーの最小化を並列して進行させる。