【計算化学】自作pythonライブラリで遷移状態構造を求めてみる(BH9データセット, 1. Radical rearrangement and additionのNo. 8の素過程, NNP使用)

最終更新：2025-09-25

概要#

本記事では、自作ライブラリ(MultiOptPy)で、BH9データセット 1. Radical rearrangement and addition No. 8の素過程の遷移状態構造を算出してみる。計算レベルは、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.16

環境#

Windows 11

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

Source codeのダウンロード#

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

環境構築手順#

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

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

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

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

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

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

※自作ライブラリでの具体的な使用の仕方に関しては、ase_calculation_tools.py を参照

手順#

1. 初期構造の準備#

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

14
OptimizedStructure
C     -0.904432576392      0.899025254761      1.593301981759
H     -1.671818540601      1.142311887901      0.843263502027
N     -0.823283325120      1.599944495125      2.611557650057
C     -0.004335739072     -0.283911135123      1.300319080644
H     -0.642671524114     -1.156422204024      1.141839858113
H      0.619513522265     -0.470738628441      2.173660066160
C      0.872031481803     -0.062987228996      0.058135069510
H      1.535682122449      0.785322721267      0.244292134272
H      1.503067253894     -0.943712233563     -0.075515905991
C      0.081570545438      0.175633765857     -1.194211233682
H     -0.440558693494      1.125416755870     -1.271165083990
C     -0.016559684248     -0.690350926250     -2.191345103599
H      0.497476218104     -1.644805712350     -2.161169416336
H     -0.605681060912     -0.474726812034     -3.072962598943

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

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

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

python optmain.py .\bh9_1_8.xyz -ma 150 3 10 -os uma-s-1 -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-1は今回使用するNNPを指定している。これを使用する際にASEライブラリが必要である。
-fc XはX回の反復計算ごとに正確なHessianを計算することを指定する。Hessianの計算コストは勾配の場合と比べると大きいが、構造最適化の収束までの時間を短縮できる。

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

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

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

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

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

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

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

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

python nebmain.py bh9_1_8_traj.xyz -os uma-s-1 -ns 15 -aneb 3 5 -modelhess -spng -ndb 0.2 -spin 2 -elec 0

-ndb Nはノード間の距離をN Åとして初期経路を作成することを示す。経路作成時に元のノードをベルンシュタイン多項式を用いてがたついた経路を滑らかにする。
-ns nはn回分NEB法による経路の緩和を行うことを示す。
-fc MはM回あたりの経路緩和回数に対して1回だけ正確なHessianを計算し、経路緩和に使用する。これを使用すると、Hessianを使用しない場合の経路緩和アルゴリズムとは別のものを使用して、経路緩和を行う。
-spngは緩和中のパスのエネルギープロファイルや各ノードの勾配のRMS値をmatplotlibで可視化するオプションである。
-os uma-s-1は今回使用するNNPを指定している。これを使用する際にASEライブラリが必要である。
-aneb A Bこれを指定すると、(B+1)回の緩和ごとに、エネルギー極大値を示すノードと前後のノードの間に線形補間でA個の新規ノードを内挿するようにできる。デフォルトではこのような操作は行われない。このオプションを使用するとノードの数が徐々に増えるため、計算コストが使用しない場合と比べて増加する。一方で、エネルギー極大値を示すノード周辺にノードを増加させるため、緩和している経路中のノードが遷移状態構造付近に存在する可能性が高くすることが出来る。

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

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

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

NEB計算の結果の可視化

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

経路緩和の結果、経路のエネルギー極大値を示す構造の中から目視で、緩和後に得られた経路の32番(グラフ上では33番)の構造を遷移状態構造を求める初期構造として採用した。

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

bh9_1_8_traj_32.xyz

14
0 2
C      -1.089019112139      0.528115273224      1.154924361569
H      -1.824238309385      0.594949218350      1.966268054270
N      -1.140874272559      1.312018343469      0.182791792649
C      -0.011538065312     -0.528250756629      1.214060332766
H      -0.433699193421     -1.490498033214      0.905357308915
H       0.362916062115     -0.644090223861      2.229637739026
C       1.048805051347     -0.020410158048      0.242944198943
H       1.606917316519      0.782361498182      0.723038022684
H       1.751881911339     -0.807947730024     -0.032654971893
C       0.320471650742      0.529892020130     -0.961839753208
H       0.575746059740      1.518002071281     -1.318867084417
C      -0.237671507259     -0.346500911595     -1.859519559242
H      -0.288155188934     -1.406366463343     -1.638524525969
H      -0.641542402792     -0.021274147921     -2.807615916096

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

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

遷移状態構造を求めるための初期構造を含むxyzファイルをMultiOptPy-v1.16と同じディレクトリ内にコピーする。

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

python optmain.py .\bh9_1_8_traj_32.xyz -spin 2 -elec 0 -os uma-s-1 -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として保存されている。)

14
OptimizedStructure
C     -1.115168329763      0.499887954754      1.131464431378
H     -1.841386511898      0.591172332862      1.946655576183
N     -1.228673727394      1.230104805610      0.124993168624
C      0.006938286158     -0.508523913072      1.246613725119
H     -0.388192423974     -1.496367599387      0.991098184425
H      0.386391062450     -0.571366679530      2.267550299244
C      1.057024981106     -0.034128978911      0.246534548851
H      1.658472105013      0.751756332501      0.703643996203
H      1.729733849178     -0.844972804092     -0.035738377475
C      0.343944502074      0.516814805722     -0.967760393852
H      0.609233121280      1.514958920230     -1.290851691000
C     -0.223963350315     -0.323704246164     -1.884474475849
H     -0.343864057847     -1.380413441481     -1.677681844658
H     -0.650489506069      0.054782510958     -2.802047147193

遷移状態構造

停留点に収束した分子構造が得られた。-freqオプションにより生成されたnormal_modes.txtやvibration_animationディレクトリ内の振動モードのアニメーションを確認した。

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

Mode                                 0                   1                   2
Freq [cm^-1]                     -613.4949            139.1904            232.6260
Reduced mass [au]                  10.2475              2.4067              2.3491
Force const [Dyne/A]               -2.2724              0.0275              0.0749
Char temp [K]                       0.0000            200.2639            334.6969
Normal mode                   x         y         z            x         y         z            x         y         z     
       C                -0.02556    0.00622    0.02790   -0.00842   -0.03360   -0.04358   -0.04484   -0.06682   -0.08192
       H                 0.00197   -0.05199    0.06564   -0.03085   -0.07994   -0.05945   -0.14581   -0.17657   -0.15820
       N                -0.12130    0.06344    0.10022    0.05819    0.06031    0.01561   -0.00360   -0.01996   -0.05208
       C                -0.00540    0.01290   -0.02033   -0.02537   -0.06355   -0.09733    0.03884    0.04089    0.04165
       H                -0.00987    0.01584   -0.02020   -0.00648   -0.02814   -0.26275    0.12600   -0.03185    0.19008
       H                -0.01228    0.01469   -0.01757   -0.08878   -0.20186   -0.08270    0.01893    0.21986    0.06136
       C                 0.01015   -0.00473   -0.01000    0.01836    0.08722    0.02592    0.04825   -0.02693    0.01266
       H                -0.01248    0.01550   -0.01448   -0.13685    0.16191    0.10180    0.10494   -0.05491   -0.01488
       H                 0.03369   -0.01482    0.07233    0.15313    0.18494    0.06496   -0.00105   -0.07270    0.02496
       C                 0.14484   -0.09169   -0.12557    0.02275    0.01175   -0.00629    0.08179   -0.00477    0.00351
       H                -0.06667    0.00610    0.01331    0.04663   -0.02779   -0.10756    0.15161   -0.02557   -0.00350
       C                 0.02280    0.00377    0.00267   -0.05491   -0.05880    0.10606   -0.10613    0.07280    0.06038
       H                 0.00504   -0.00071   -0.01342   -0.10694   -0.02864    0.22926   -0.28879    0.10021    0.09373
       H                -0.00235    0.00943    0.01429   -0.07167   -0.13987    0.08083   -0.12890    0.13819    0.09809

(...snip...)

その結果、虚振動が１つであることが確認できた。つまりこの構造は遷移状態構造である。

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

終わりに#

　　自作ライブラリで、UMAモデルのニューラルネットワークポテンシャル(NNP)を用いて、BH9データセット 1. Radical rearrangement and addition No. 8の素過程のある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のデータセットについて)

概要#