最終更新：2025-09-26
概要
本記事では、自作ライブラリ(MultiOptPy)で、BH9データセット 1. Radical rearrangement and addition No. 10の素過程の遷移状態構造を算出してみる。計算レベルは、Meta社のFAIR Chemistryが開発したニューラルネットワークポテンシャル(NNP)であるUMA(Meta’s Universal Model for Atoms)とした。
MultiOptPyは電子状態計算ソフトウェアを用いた分子構造最適化手法の勉強を目的として作成したpythonライブラリである。
MultiOptPyのレポジトリ:https://github.com/ss0832/MultiOptPy
BH9のデータセットについて:
- J. Chem. Theory Comput. 2022, 18, 1, 151–166 https://doi.org/10.1021/acs.jctc.1c00694
今回使用したニューラルネットワークポテンシャルについて：
- 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.17.2環境
Windows 11※Windows 11環境下でAnaconda PowerShell Promptを使用した。
Source codeのダウンロード
wget https://github.com/ss0832/MultiOptPy/archive/refs/tags/v1.17.2.zip
unzip v1.17.2.zip
cd MultiOptPy-v1.17.2環境構築手順
今回は、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の自作ライブラリへの導入方法について説明する。
1, 以下のサイトにアクセスして、uma-s-1p1.ptをダウンロードする。（使用許諾が下りていれば可能である。）
https://huggingface.co/facebook/UMA2, ダウンロード後、MultiOptPy-v1.17.2ディレクトリ内に存在するsoftware_path.confに対して、uma-s-1p1.ptの絶対パスを用いて以下を追記する。
uma-s-1p1::(uma-s-1p1.ptの絶対パス)これで、MultiOptPy-v1.17.2がNNPuma-s-1p1を使用できるようになる。
使用するNNPに関する具体的な説明
今回使用するNNPについて具体的に説明する。
- UMAのModel Checkpointは
uma-s-1p1を使用した。 - 小分子系のトレーニングセットである
Omol25(omol)を使用して学習したニューラルネットワークポテンシャルを使用する。
※自作ライブラリでの具体的な使用の仕方に関しては、ase_calculation_tools.py を参照
手順
1. 初期構造の準備
モデル反応系として、以下の構造を用意した。今回はファイルの名前をbh9_1_10.xyzとした。 初期構造は以下のものを使用した。（uma-s-1を使って構造最適化したものを使用している。）
18
OptimizedStructure
N -0.324494291275 -0.305349020975 1.831207404966
C -1.309136612073 0.566576325643 2.394976894185
H -1.183825864889 1.615676729071 2.124680803875
H -2.294354043368 0.198357322999 2.095272309350
H -1.276265399255 0.469613487376 3.485824016604
C 0.770093009018 0.301191854382 1.231467923390
O 1.212525633288 1.352582881146 1.654035834525
C 1.420632345878 -0.496844491130 0.128795512231
H 2.162572290389 0.144610322484 -0.344840218691
H 1.943450194226 -1.329423684290 0.606561203957
C 0.415028340091 -1.036287333176 -0.894963793976
H 0.942299351682 -1.640284683723 -1.634898769856
H -0.289724855004 -1.695253907812 -0.380624892816
C -0.334036826543 0.065352788951 -1.582366968136
H -0.917790296972 0.724059934286 -0.942173646456
C -0.317597370902 0.292624850722 -2.886420379511
H 0.252525759303 -0.337945641538 -3.559542357397
H -0.871901363595 1.110742265587 -3.326990876245
2. 遷移状態構造最適化のための初期構造の算出
a. 初期構造をカレントディレクトリにbh9_1_10.xyzとして保存し、以下を実行する。
python optmain.py .\bh9_1_10.xyz -ma 150 1 14 -os uma-s-1p1 -fc 50 -opt rsirfo_fsb -spin 2 -elec 0 -opt rsirfo_fsbは準ニュートン法であるRS-I-RFO法を構造最適化に使用することを示す。初期のへシアンに関しては、特にオプションで指定しない限り、単位行列が使われる。-spinはスピン多重度の指定である。PySCFを使用するときは目的とするスピン多重度に1を引いた値を指定する。（デフォルトでは1が指定される。）今回はUMAモデルのNNPを使用し、ラジカル反応を扱うためスピンだ重度を2とした。-elec Mは形式電荷をMとすることを示す。（デフォルトでは0が指定される。）-ma yyy a bはyyykJ/molの活性化障壁を超えうるペア同士を近づける力を原子のラベル番号aとbのペアに構造最適化時に加えることを示す。-os uma-s-1p1は今回使用するNNPを指定している。これを使用する際にASEライブラリが必要である。-fc XはX回の反復計算ごとに正確なHessianを計算することを指定する。Hessianの計算コストは勾配の場合と比べると大きいが、構造最適化の収束までの時間を短縮できる。
これを実行すると、omolのデータセットを使用したuma-s-1p1モデルのNNPで得たエネルギーに対して、指定した人工力ポテンシャルを加えた上で初期構造を構造最適化することができる。
結果はyyyy_mm_dd（今日の年月日）ディレクトリの中に存在するディレクトリを開いて確認できる。
正常終了していれば、このディレクトリ中に、bh9_1_10_traj.xyzが存在するので、これをコピーして、MultiOptPy-v1.17.2ディレクトリに置く。
bh9_1_10_traj.xyzは構造最適化の過程をAvogadro（公式ページ:https://avogadro.cc/ ）等で可視化して確認できるようにしている。このbh9_1_10_traj.xyzは次のNEB計算に使用する。
※bh9_1_10_traj.xyzをアニメーションとして表示したい場合は、[https://github.com/ss0832/molecule_movie] を使うと良い。
（この人工力ポテンシャルを加えて行った構造最適化の結果はbh9_1_10_optimized.xyzで確認できる。構造を可視化して、生成系になっているか確認する。反応系のままであれば、-maの設定を見直してb.をやり直す。今回の場合以下の構造が得られた。生成系の安定構造に近いものが得られていると判断できるため、次のNEB計算を行うことが可能である。）
b. bh9_1_10_traj.xyzを初期パスとして、NEB法で経路の緩和を行う。
NEB法を用いることで、先ほど得られたbh9_1_10_traj.xyz全体のエネルギーを下げることができる。これにより、パスのエネルギー極大値を持つ構造を遷移状態構造に近づける。（この時点ではまだ正確な遷移状態構造は求められていない。）
python nebmain.py bh9_1_10_traj.xyz -os uma-s-1p1 -ns 15 -aneb 3 5 -modelhess -spng -ndb 0.5 -ad 9999 -spin 2 -elec 0-ad Xは線形補間で、各ノード間の距離を全て等しくするための処理である。X回の反復計算ごとに本処理を行う。Xを-nsよりも大きな数値を指定することで、初期経路に対してのみ処理を行うことが出来る。-ndb Nはノード間の距離をN Åとして初期経路を作成することを示す。経路作成時に元のノードをベルンシュタイン多項式を用いてがたついた経路を滑らかにする。
→プログラムの仕様上-adの処理を行った後に、-ndbの処理を行うようになっている。
-ns nはn回分NEB法による経路の緩和を行うことを示す。-fc MはM回あたりの経路緩和回数に対して1回だけ正確なHessianを計算し、経路緩和に使用する。これを使用すると、Hessianを使用しない場合の経路緩和アルゴリズムとは別のものを使用して、経路緩和を行う。-spngは緩和中のパスのエネルギープロファイルや各ノードの勾配のRMS値をmatplotlibで可視化するオプションである。-os uma-s-1p1は今回使用するNNPを指定している。これを使用する際にASEライブラリが必要である。-aneb A Bこれを指定すると、(B+1)回の緩和ごとに、エネルギー極大値を示すノードと前後のノードの間に線形補間でA個の新規ノードを内挿するようにできる。デフォルトではこのような操作は行われない。このオプションを使用するとノードの数が徐々に増えるため、計算コストが使用しない場合と比べて増加する。一方で、エネルギー極大値を示すノード周辺にノードを増加させるため、緩和している経路中のノードが遷移状態構造付近に存在する可能性が高くすることが出来る。
c. 初期構造の決定及び遷移状態構造の計算
MultiOptPy-v1.17.2と同じディレクトリ内に、NEBという名前を含むディレクトリが生成されている。 そのディレクトリ内のenergy_plot.csvを確認し、緩和後のパスのエネルギー極大値を示す構造を確認する。
パスの緩和後の各ノードのエネルギー一覧(単位
※bias_force_rms.csvにて、各Iterationごとのすべてのノードの勾配のRMS値を確認できる。
経路緩和の結果、経路のエネルギー極大値を示す構造の中から目視で、緩和後に得られた経路の30番(グラフ上では31番)の構造を遷移状態構造を求める初期構造として採用した。
※こちら[https://ss0832.github.io/molecule_viewer/] を使うことでも可視化は可能である。
bh9_1_10_traj_30.xyz
18
0 2
N -0.473105008466 0.085781788392 0.956731280021
C -1.325776895760 0.970167497690 1.699740594087
H -1.037724157016 2.016391746837 1.604265945758
H -2.360848970054 0.821641123161 1.394911704657
H -1.248938278821 0.703376959776 2.759847076298
C 0.872950435714 0.168440441047 1.126795516588
O 1.399437722807 0.714111995582 2.081293797113
C 1.601745741735 -0.582131764521 0.040746645898
H 2.154034220292 0.152649878266 -0.551178681528
H 2.339404205848 -1.246240015822 0.492811359273
C 0.557427777414 -1.313862437039 -0.807000859864
H 0.989772709819 -1.602389962725 -1.765001376723
H 0.228673726958 -2.214507790658 -0.292143513533
C -0.622976602327 -0.383889789934 -1.028848570107
H -1.626982645217 -0.670085055427 -0.744603668897
C -0.512246706370 0.556363684516 -2.022628099557
H 0.412307506170 0.671055058436 -2.581201076462
H -1.347154782725 1.153126642424 -2.364538073022
構造が壊れていないので、これを遷移状態を求めるための初期構造とする。
遷移状態構造を求めるための初期構造を含むxyzファイルをMultiOptPy-v1.17.2と同じディレクトリ内にコピーする。
そして、以下を実行する。
python optmain.py .\bh9_1_10_traj_30.xyz -spin 2 -elec 0 -os uma-s-1p1 -fc 5 -tcc -freq -opt rsirfo_bofill -order 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と同等）
実行して得られた正確な遷移状態構造と思われる構造を以下に示す。
(実行して得られた正確な遷移状態構造は計算開始時に、yyyy_mm_ddディレクトリ内に生成された新規ディレクトリ内のYYY_optimized.xyzとして保存されている。)
18
OptimizedStructure
N -0.482333631312 -0.009759047398 0.997336739450
C -1.313213525327 1.026812142082 1.549332519037
H -0.927877386847 2.031926173844 1.368442763547
H -2.319925259540 0.927515605825 1.144178008544
H -1.374942042728 0.880052201422 2.631526771138
C 0.877572550300 0.176768962284 1.035006096511
O 1.436884384126 0.883931140353 1.852259995937
C 1.607709980765 -0.641673342897 -0.011411412316
H 2.176533388369 0.058332933066 -0.627587842154
H 2.336717927696 -1.285894630247 0.479752389530
C 0.584129864598 -1.421041606048 -0.830604425683
H 1.003896025335 -1.749665047852 -1.783410014437
H 0.274099937551 -2.310216765019 -0.283754416350
C -0.630618587655 -0.549350822429 -1.067820197731
H -1.601806520927 -1.014848562406 -0.953366623271
C -0.566523151100 0.655797737076 -1.678273168939
H 0.381296712733 1.129080609921 -1.903931050985
H -1.461600666038 1.212232318421 -1.917676131829
停留点に収束した分子構造が得られた。-freqオプションにより生成されたnormal_modes.txtやvibration_animationディレクトリ内の振動モードのアニメーションを確認した。
以下に-freqオプションで生成されたnormal_modes.txtの一部を示す。
Mode 0 1 2
Freq [cm^-1] -369.8127 80.8194 103.7522
Reduced mass [au] 8.2892 4.4685 1.0351
Force const [Dyne/A] -0.6679 0.0172 0.0066
Char temp [K] 0.0000 116.2810 149.2762
Normal mode x y z x y z x y z
N -0.00647 -0.01301 0.15634 -0.01009 -0.05057 0.01873 -0.00654 -0.01060 0.00286
C 0.00451 0.02994 0.04123 -0.02036 -0.06701 0.03449 0.00615 0.01201 -0.02116
H 0.03157 0.01405 0.01329 -0.05704 -0.05782 0.00759 -0.25467 0.02817 -0.49622
H 0.02103 0.01878 0.00241 -0.03356 -0.10430 0.07646 -0.13270 -0.30629 0.40867
H -0.05701 0.07965 0.04766 0.02758 -0.05392 0.03920 0.48449 0.36645 0.05454
C -0.00072 0.01333 0.01970 -0.01291 0.00547 -0.02850 -0.00466 -0.01010 0.00934
O 0.02109 0.01080 0.00264 -0.01498 0.10770 -0.11678 -0.00877 -0.01844 0.01857
C -0.02506 0.01314 0.00022 -0.01275 -0.07106 0.03023 0.00030 0.00222 0.00317
H -0.01994 0.01372 0.00581 0.07238 -0.13094 0.04238 -0.00163 0.01112 0.01143
H -0.02346 0.01440 0.00225 -0.08655 -0.11172 0.08758 0.00226 0.00058 -0.00180
C -0.00163 -0.00360 -0.01708 -0.02928 -0.02492 0.00466 0.00514 0.00733 -0.00718
H 0.07500 -0.07056 0.03887 -0.02742 -0.03246 0.00796 0.01040 0.01464 -0.00745
H -0.03736 0.02751 0.01403 -0.07891 -0.01650 -0.01061 0.00463 0.00260 -0.01508
C 0.00289 -0.04671 -0.20667 0.01456 0.03301 -0.00309 0.00315 0.00484 -0.00817
H 0.00213 -0.02202 -0.11397 -0.00923 0.06638 -0.06706 0.00453 -0.00037 -0.01689
C 0.00051 -0.01292 -0.02515 0.08865 0.06414 0.06457 -0.00027 0.00902 0.00044
H -0.00370 0.02612 0.02981 0.11708 0.04076 0.13498 -0.00166 0.01704 0.01079
H -0.00094 -0.01105 -0.01885 0.12121 0.11278 0.05596 -0.00244 0.00448 -0.00199
(...snip...)その結果、虚振動が１つであることが確認できた。つまりこの構造は遷移状態構造である。
次に、vibration_animation内の虚振動を示す分子振動が示されたxyzファイル(mode_1_XXXi_wave_number.xyz)をAvogadroで確認すると、想定される反応系と生成系をつなぐ方向に振動していることを確認できた。
終わりに
自作ライブラリで、UMAモデルのニューラルネットワークポテンシャル(NNP)を用いて、BH9データセット 1. Radical rearrangement and addition No. 10の素過程のある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.
- J. Chem. Theory Comput. 2022, 18, 1, 151–166 https://doi.org/10.1021/acs.jctc.1c00694 (BH9のデータセットについて)