最終更新：2025-08-19

概要#

本記事では、自作モジュール(MultiOptPy)で、SN2反応（基質-methyl trifluoromethanesulfonate, 求核剤-fluoride 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.12c

環境#

Windows 11

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

Source codeのダウンロード#

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

環境構築手順#

今回は、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.12cディレクトリ内に存在するsoftware_path.confに対して、uma-s-1.ptの絶対パスを用いて以下を追記する。

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

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

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

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

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

※自作モジュールでの具体的な使用の仕方に関しては、ase_calculation_tools.py を参照

手順#

1. 初期構造の準備#

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

13
OptimizedStructure
C     -1.397191095072     -0.503357183150      0.812700352725
H     -1.275480258399     -1.573355838428      0.669452231922
H     -1.074937104241     -0.192948081206      1.803741774349
H     -2.428050081873     -0.212143918950      0.652264634483
O     -0.639231977547      0.209550759498     -0.197220664975
F     -2.224768956906     -1.991056691866      3.353635027375
S      0.921855415126      0.127148446561     -0.061265529786
O      1.387530637631      0.989143892904      0.970578168528
O      1.364685860616     -1.223721479809     -0.143594728573
C      1.266598487840      0.945578341758     -1.675913284471
F      2.574721006165      1.161272565667     -1.733141545725
F      0.622030392065      2.100172002919     -1.769294657694
F      0.902237674595      0.163717184100     -2.681941778158

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

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

python .\optmain.py sn2_TfO_F.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_TfO_F_traj.xyzが存在するので、これをコピーして、MultiOptPy-v1.12cディレクトリに置く。

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

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

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

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

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

python .\nebmain.py sn2_TfO_F_traj.xyz -os uma-s-1 -ns 20 -modelhess -spng -nd 0.12 -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.12cと同じディレクトリ内に、NEBという名前を含むディレクトリが生成されている。 そのディレクトリ内のenergy_plot.csvを確認し、緩和後のパスのエネルギー極大値を示す構造を確認する。

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

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

経路緩和の結果、経路のエネルギー極大値を示す構造付近の勾配のRMS値の極小値を示す構造の中から目視で、ITR. 5の経路の6番(グラフ上では7番)の構造を遷移状態構造を求める初期構造として採用して遷移状態構造を求めようとした。

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

sn2_TfO_F_traj_6_itr5.xyz

13
-1 1
C      -1.511561824733     -0.650405908653      1.059011640615
H      -1.277774648086     -1.644658308534      0.744044920760
H      -1.061812058260     -0.218725100586      1.922250957206
H      -2.459040967070     -0.255533119587      0.726864313635
O      -0.623742469011      0.202616392612     -0.198204112761
F      -2.133182708110     -1.719876830408      2.874123820794
S       0.867675848268      0.125010746229     -0.057197606459
O       1.401601477930      0.991673376997      0.951692701267
O       1.385886912357     -1.207186879224     -0.156587548376
C       1.272470916866      0.935556352985     -1.653949590552
F       2.586370894415      1.162979650694     -1.750239038706
F       0.642358126611      2.109612552660     -1.773279459566
F       0.910750498823      0.168937074816     -2.688530997858

しかしながら、うまくいかなかったため、エネルギー極大値をとる構造付近の勾配のRMS値の最小値をとる構造の生成物側で1つ隣のノードを使用したが、これもうまくいかなかった。

sn2_TfO_F_traj_7_itr5.xyz

13
-1 1
C      -1.575460305272     -0.714185348066      1.158231034818
H      -1.269353581918     -1.651900118057      0.744578804961
H      -1.043644517785     -0.208711249385      1.930561658658
H      -2.458981422539     -0.255020465798      0.726210884846
O      -0.615041927582      0.212337926513     -0.210027635310
F      -2.119304922788     -1.677300748695      2.795380247976
S       0.859396123617      0.127858674309     -0.061838792554
O       1.403964116075      0.991756277847      0.948840207173
O       1.385155896589     -1.205188636698     -0.159900305408
C       1.279187943066      0.934154650101     -1.651154105068
F       2.591312674040      1.164924178761     -1.753679551305
F       0.647486150382      2.110899949870     -1.776153011157
F       0.915283774116      0.170374909299     -2.691049437631

そこで以下のような操作を行った。 a, sn2_TfO_F_optimized.xyzをsn2_TfO_F_reverse.xyzに改名して以下のようなコマンドを実行した。

python .\optmain.py .\sn2_TfO_F_reverse.xyz -opt fire -os uma-s-1 -ma 465 1 5 -elec -1 -ns 120

• -opt fireはFIREと呼ばれる最急降下法と似た連続最適化アルゴリズムを用いて構造最適化を行うことを指定する。準ニュートン法だと、経験上構造最適化途中にエネルギーが上がることがあるので、時間をかけてもエネルギーを確実に下げたいときに使用する。

sn2_TfO_F_reverse.xyz

13
OptimizedStructure
C     -1.897507649579     -1.095752431252      1.777177253589
H     -1.359856126639     -1.800472537204      1.143734781875
H     -1.162595454619     -0.421467690576      2.222699471455
H     -2.578216498150     -0.518712844367      1.147996749121
O     -0.583031904297      0.405267695015     -0.486886653381
F     -2.570533037711     -1.743202183440      2.719345066672
S      0.847662941267      0.291905918664     -0.289621221682
O      1.397541724792      1.142135344544      0.742882304827
O      1.370091193190     -1.055964554320     -0.338915021399
C      1.503511478737      1.045206379481     -1.842246655769
F      2.845906745211      1.128164363572     -1.841413327401
F      1.032768955558      2.291556077671     -2.027398692353
F      1.154257632240      0.331336462213     -2.927354055555

b, 前のコマンドで得られたITR. 105の構造を記録したファイルsn2_TfO_F_reverse_0.xyzをsn2_TfO_F_reverse_itr105.xyzに改名して以下のコマンドを実行した。

python .\optmain.py .\sn2_TfO_F_reverse_itr105.xyz -opt fire -os uma-s-1  -elec -1 -ns 120

sn2_TfO_F_reverse_itr105.xyz

13
-1 1
C      -1.300518756379     -0.586193055337      0.967823990509
H      -1.181623012881     -1.668956322503      0.907693854227
H      -1.026471075608     -0.243925991362      1.968108697102
H      -2.341693141130     -0.333932169545      0.788543365695
O      -0.588372240255      0.053021970705     -0.013398668219
F      -2.614029408724     -1.592946480837      2.526333197407
S       0.944547869816      0.118681919332     -0.025751689109
O       1.445596062743      1.054286872172      0.927771220091
O       1.542034242555     -1.169330392616     -0.166304689580
C       1.200020090335      0.937806751520     -1.662431855159
F       2.522749288990      1.085294340393     -1.806109726621
F       0.641897202276      2.145188555490     -1.733083730600
F       0.755862878262      0.201004002587     -2.679193965743

結果を以下に示す。

構造最適化の結果

13
OptimizedStructure
C     -1.927086981613     -1.073086342294      1.723136702839
H     -1.366831081267     -1.819625340378      1.167028798524
H     -1.245166080738     -0.393938041283      2.227576119285
H     -2.576359300860     -0.517676251288      1.052706538680
O     -0.503691597723      0.335393773851     -0.364533794898
F     -2.719690800809     -1.721748597992      2.689001387479
S      0.941060182542      0.276226969763     -0.252120356174
O      1.519826538813      1.152811273399      0.741980328853
O      1.510948423435     -1.051336917438     -0.328422482243
C      1.482853998376      1.046706277557     -1.841405982691
F      2.822531514075      1.128608088478     -1.929874074085
F      1.001867888419      2.295479802532     -1.978719240123
F      1.059737297351      0.342185305093     -2.906353945444

生成物の構造をSN2反応の遷移状態構造を経由して反応物の構造で収束した。本記事の操作では遷移状態構造を見つけることはできなかった。

終わりに#

　 　自作モジュールで、UMAモデルのニューラルネットワークポテンシャル(NNP)を用いて、SN2反応（基質-methyl Trifluoromethanesulfonate, 求核剤-fluoride 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.

個人的な技術的補足#

-modelhessを使わない場合のアルゴリズムを用いたnebmain.pyによる経路緩和の結果

コマンド#

python .\nebmain.py sn2_TfO_F_traj.xyz -os uma-s-1 -ns 20 -spng -nd 0.12 -elec -1

結果#

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

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

緩和後の経路のエネルギー極大値を示すノードの分子構造 sn2_TfO_F_traj_5.xyz

13
-1 1
C      -1.518212863552     -0.624951600833      1.009912595544
H      -1.289674642173     -1.619897872367      0.693916033805
H      -1.091575720215     -0.213920992966      1.889330544710
H      -2.460328061626     -0.218329695223      0.659838834166
O      -0.623955489214      0.198175266960     -0.173633386275
F      -2.147236570262     -1.804623206232      3.037126405412
S       0.883970047969      0.125016136106     -0.051529909165
O       1.412484219300      0.995407443482      0.952854416990
O       1.389682170676     -1.208693371716     -0.158440340161
C       1.279482572028      0.936547730253     -1.655338139082
F       2.596335520730      1.161240428213     -1.748808394956
F       0.649551329339      2.108319746800     -1.771434572635
F       0.919477487001      0.165709987523     -2.683794088352