【計算化学】自作pythonモジュールで遷移状態構造を求めてみる(SN2反応（基質-chloromethane, 求核剤-bromide ion）, NNP使用）

最終更新：2025-08-14

概要#

本記事では、自作モジュール(MultiOptPy)で、SN2反応（基質-chloromethane, 求核剤-bromide ion）の遷移状態構造を算出してみる。使用する計算レベルは、Meta社のFAIR Chemistryが開発したニューラルネットワークポテンシャル(NNP)であるUMA(Meta’s Universal Model for Atoms)を使用する。

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

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

SN2反応について：

https://en.wikipedia.org/wiki/SN2_reaction

今回使用したニューラルネットワークポテンシャルについて：

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へのアカウント登録と配布元の使用許諾が必要である。)

使用した自作モジュールMultiOptPyのバージョン#

v1.12b

環境#

Windows 11

※Windows 11環境下でAnaconda PowerShell Promptを使用した。

Source codeのダウンロード#

wget https://github.com/ss0832/MultiOptPy/archive/refs/tags/v1.12b.zip
unzip v1.12b.zip
cd MultiOptPy-v1.12b

環境構築手順#

今回は、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.7

4, 先ほど作成した仮想環境をconda activate (仮想環境名)で起動させる。

5, 以下のコマンドを実行し、必要なライブラリを導入する。

pip install fairchem-core==2.2.0 ase==3.25.0 torch==2.6.0

fairchem-coreは、FAIR Chemistryが管理しているNNPを動作させるために必要なライブラリである。
aseはNNPに電子エネルギーを算出したい分子構造を渡すために必要なインターフェイスの役割を果たすために必要なライブラリである。
torchはPyTorchライブラリを指す。これはニューラルネットワークなどの機械学習を行ったり、学習結果を扱ったりするために必須なライブラリである。

これで、Anaconda PowerShell Promptから仮想環境を立ち上げることで、NNPを使用する準備が整えることが出来る。

次に、NNPを使用するために必要なModelの情報が保存されている.ptファイルのダウンロードおよびNNPの自作モジュールへの導入方法について説明する。

1, 以下のサイトにアクセスして、uma-s-1.ptをダウンロードする。（使用許諾が下りていれば可能である。）

https://huggingface.co/facebook/UMA

2, ダウンロード後、MultiOptPy-v1.12bディレクトリ内に存在するsoftware_path.confに対して、uma-s-1.ptの絶対パスを用いて以下を追記する。

uma-s-1::(uma-s-1.ptの絶対パス)

これで、MultiOptPy-v1.12bがNNPuma-s-1を使用できるようになる。

使用するNNPに関する具体的な説明#

今回使用するNNPについて具体的に説明する。

UMAのModel Checkpointはuma-s-1を使用した。
小分子系のトレーニングセットであるOmol25(omol)を使用して学習したニューラルネットワークポテンシャルを使用する。

※自作モジュールでの具体的な使用の仕方に関しては、https://github.com/ss0832/MultiOptPy/blob/main/multioptpy/ase_calculation_tools.py を参照

手順#

1. 初期構造の準備#

モデル反応系として、以下の構造を用意した。今回はファイルの名前をsn2_Cl_Br.xyzとした。初期構造は以下のものを使用した。

6

C      0.134497662922      0.262821827927     -0.436820313373
H      0.496066659825     -0.760059576032     -0.442386154475
H      0.498156268875      0.784387618101      0.442190929474
H     -0.949778075901      0.279794049359     -0.463902893826
Br    -0.927075879399     -1.661560687040      2.789829085360
Cl     0.748133363678      1.094616767686     -1.888910653160

初期パスを求めるための初期構造

2. 遷移状態構造最適化のための初期構造の算出#

a. 初期構造をカレントディレクトリに`sn2_Cl_Br.xyz`として保存し、以下を実行する。#

python .\optmain.py sn2_Cl_Br.xyz -opt rsirfo_fsb -os uma-s-1 -ma 300 1 6 -elec -1

-opt rsirfo_fsbは準ニュートン法であるRS-I-RFO法を構造最適化に使用することを示す。初期のへシアンに関しては、特にオプションで指定しない限り、単位行列が使われる。
-spinはスピン多重度の指定である。PySCFを使用するときは目的とするスピン多重度に1を引いた値を指定する。（デフォルトでは1が指定される。）
-elec Mは形式電荷をMとすることを示す。（デフォルトでは0が指定される。）
-ma yyy a bはyyykJ/molの活性化障壁を超えうるペア同士を近づける力を原子のラベル番号aとbのペアに構造最適化時に加えることを示す。
-os uma-s-1は今回使用するNNPを指定している。これを使用する際にASEライブラリが必要である。

これを実行すると、omolのデータセットを使用したuma-s-1モデルのNNPで得たエネルギーに対して、指定した人工力ポテンシャルを加えた上で初期構造を構造最適化することができる。

結果はyyyy_mm_dd（今日の年月日）ディレクトリの中に存在するディレクトリを開いて確認できる。

正常終了していれば、このディレクトリ中に、sn2_Cl_Br_traj.xyzが存在するので、これをコピーして、MultiOptPy-v1.12bディレクトリに置く。

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

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

（この人工力ポテンシャルを加えて行った構造最適化の結果はsn2_Cl_Br_optimized.xyzで確認できる。構造を可視化して、生成系になっているか確認する。反応系のままであれば、-maの設定を見直してb.をやり直す。今回の場合以下の構造が得られた。生成系の安定構造に近いものが得られていると判断できるため、次のNEB計算を行うことが可能である。）

6
OptimizedStructure
C     -0.143628752597     -0.204053732939      0.353636691009
H      0.456449829219     -0.843657824379     -0.282905564139
H      0.451667759437      0.665642198150      0.605934036264
H     -0.998191428987      0.131812632834     -0.222039729670
Cl     0.912109528675      1.338234788200     -2.308136112783
Br    -0.678406935747     -1.087978061866      1.853510679319

初期パスを求めるための構造最適化の結果（安定構造ではない）

b. `sn2_Cl_Br_traj.xyz`を初期パスとして、NEB法で経路の緩和を行う。#

NEB法を用いることで、先ほど得られたsn2_Cl_Br_traj.xyz全体のエネルギーを下げることができる。これにより、パスのエネルギー極大値を持つ構造を遷移状態構造に近づける。（この時点ではまだ正確な遷移状態構造は求められていない。）

python .\nebmain.py sn2_Cl_Br_traj.xyz -os uma-s-1 -ns 10 -modelhess -spng -nd 0.15 -elec -1

-nd Nはノード間の距離をN Åとして初期パスを作成することを示す。
-ns nはn回分NEB法による経路の緩和を行うことを示す。
-fc MはM回あたりの経路緩和回数に対して1回だけ正確なHessianを計算し、経路緩和に使用する。これを使用すると、Hessianを使用しない場合の経路緩和アルゴリズムとは別のものを使用して、経路緩和を行う。
-spngは緩和中のパスのエネルギープロファイルや各ノードの勾配のRMS値をmatplotlibで可視化するオプションである。
-os uma-s-1は今回使用するNNPを指定している。これを使用する際にASEライブラリが必要である。

c. 初期構造の決定及び遷移状態構造の計算#

MultiOptPy-v1.12bと同じディレクトリ内に、NEBという名前を含むディレクトリが生成されている。そのディレクトリ内のenergy_plot.csvを確認し、緩和後のパスのエネルギー極大値を示す構造を確認する。

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

NEB計算の結果の可視化

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

私が実行した環境では、9番のノード(グラフでは10番)がエネルギー極大値を示していた。また、Walden反転の図解のような反応に関わる分子構造の部位が平面となっている。

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

sn2_Cl_Br_traj_9.xyz

6
-1 1
C      -0.091519824954     -0.068085308824      0.117753702843
H       0.461421640983     -0.861741032507     -0.351215704705
H       0.455971182183      0.732718630561      0.579717987265
H      -1.049672710724      0.183771996110     -0.309874066256
Cl      0.836913767008      1.153931804428     -1.979133557690
Br     -0.613114054497     -1.140596089768      1.942751638543

NEB法により緩和したパスのエネルギー極大値を示す構造

構造が壊れていないので、これを遷移状態を求めるための初期構造とする。

sn2_Cl_Br_traj_9.xyzをMultiOptPy-v1.12bと同じディレクトリ内にコピーする。

そして、以下を実行する。

python .\optmain.py sn2_Cl_Br_traj_9.xyz -freq -tcc -opt rsirfo_bofill -fc 5 -order 1 -os uma-s-1

-opt rsirfo_bofillは遷移状態構造の最適化向けのoptimizerを指定することを意味する。準ニュートン法であるRS-I-RFO法を使用する。今回は-fcで正確なへシアンを計算するようにしているので、初期へシアンは正確なへシアンを使用するようになっている。
-order 1は一次の鞍点を求めることを指定する。（デフォルトだと極小値を求めるようになっている。）
-fc 5は5回の反復回数当たり1回正確なへシアンを計算することを指定する。
-freqは収束条件を満たした後に基準振動解析を行うことを示す。（自前で実装しているため、あくまで目安として使用することを推奨する。各振動モードをvibration_animation内のxyzファイルで可視化できる。）UMAモデルから算出されるHessianは数値微分により求めているため、原子数Zが多いとZの二乗オーダーで計算コストが急増する。
-tccは収束条件を厳しくすることを示す。（Gaussianのtightと同等）

実行して得られた正確な遷移状態構造を以下に示す。

(実行して得られた正確な遷移状態構造はsn2_Cl_Br_traj_9_optimized.xyzとして保存されている。)

6
OptimizedStructure
C     -0.003519468323     -0.005572729696      0.009523736698
H      0.518409543744     -0.870644233864     -0.346225048682
H      0.515626947067      0.730218023971      0.589536197801
H     -1.013377918261      0.173314911007     -0.299512782405
Cl     0.714870710541      1.135708676895     -1.940800181315
Br    -0.732009814768     -1.163024648313      1.987478077903

遷移状態構造

10回程度の反復計算で遷移状態構造が得られた。-freqオプションにより生成されたnormal_modes.txtやvibration_animationディレクトリ内の振動モードのアニメーションを確認した。

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

Mode                                 0                   1                   2
Freq [cm^-1]                     -397.8345            179.8777            190.7259
Reduced mass [au]                  12.8454             36.1442              4.1846
Force const [Dyne/A]               -1.1979              0.6890              0.0897
Char temp [K]                       0.0000            258.8039            274.4120
Normal mode                   x         y         z            x         y         z            x         y         z
       C                -0.08195   -0.13021    0.22251   -0.01239   -0.01954    0.03338    0.17740    0.10071    0.12429
       H                -0.01016   -0.00550    0.02092   -0.01142   -0.01566    0.02926    0.18199    0.10307    0.12718
       H                -0.01014   -0.01555    0.01504   -0.01141   -0.01785    0.02795    0.17936    0.09953    0.13116
       H                -0.00055   -0.01207    0.02064   -0.00930   -0.01709    0.02919    0.18231    0.10601    0.12252
      Cl                 0.01537    0.02442   -0.04173   -0.04048   -0.06434    0.10996   -0.03872   -0.02199   -0.02707
      Br                 0.00592    0.00940   -0.01606    0.02023    0.03213   -0.05490   -0.01676   -0.00951   -0.01177

(...snip...)

その結果、虚振動が１つであることが確認できた。

次に、vibration_animation内のmode_1_398i_wave_number.xyzをAvogadroで確認すると、想定される反応系と生成系をつなぐ方向に振動していることを確認できた。

終わりに#

　　自作モジュールで、UMAモデルのニューラルネットワークポテンシャル(NNP)を用いて、SN2反応（基質-chloromethane, 求核剤-bromide ion）のある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へのアカウント登録と配布元の使用許諾が必要である。)
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.

概要#