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

最終更新：2025-10-11

概要#

本記事では、自作ライブラリ(MultiOptPy)で、MOBH35データセットのNo. 15の素過程の遷移状態構造を算出してみる。計算レベルは、Meta社のFAIR Chemistryが開発したニューラルネットワークポテンシャル(NNP)であるUMA(Meta’s Universal Model for Atoms)とした。

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

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

MOBH35のデータセットについて:

J. Phys. Chem. A 2019, 123, 3761−3781 http://dx.doi.org/10.1021/acs.jpca.9b01546

この文献のSupporting Informationから、データセットの詳細を確認できる。

有機金属錯体が関わる化学反応の活性化障壁及び反応熱の値のDFTの汎関数のエネルギーの精度の比較のベンチマークに使用される。

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

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.18c

環境#

Windows 11

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

Source codeのダウンロード(Unixコマンド)#

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

https://github.com/ss0832/MultiOptPy/releases/tag/v1.18c にアクセスしてzipファイルをダウンロードする。Unixコマンドの場合とはディレクトリ名が異なるので都度読み替えていただけると良い。

環境構築手順#

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

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

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

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

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

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

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

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

手順#

1. 初期構造の準備#

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

19
OptimizedStructure
C     -0.354616855909     -1.492565868471      0.512974918089
C     -1.673392644636     -1.895175472601      0.444041823583
C     -2.676708297246     -1.099701909003     -0.105014420808
C     -2.487055436244      0.185935733398     -0.606910206222
C     -1.264932052912      0.827480059941     -0.575452944994
Rh     0.500069774222      0.205609346509      0.020033706530
P      1.271897829687      2.408491842413      0.061107349639
P      2.596550463794     -0.819395098835      0.090125384636
H      0.271853689573     -2.225992834541      1.032808057153
H     -1.966434948895     -2.862874725269      0.842279049700
H     -3.681463912662     -1.512602720149     -0.144573613670
H     -1.300893947678      1.803741494731     -1.071671769558
H     -3.347355393438      0.679848237941     -1.050572015383
H      3.156547439718     -1.115941334171     -1.166302589848
H      2.754353408791     -2.104364990168      0.645816830375
H      3.775506342338     -0.269865999316      0.634879186649
H      2.528880840293      2.860330034956     -0.391067178678
H      1.379687793595      2.977957961503      1.343362309072
H      0.517505907608      3.449086241131     -0.515863876266

初期経路を求めるための初期構造

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

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

python optmain.py .\mobh35_15_Rh.xyz -ma 250 1 5 -os uma-s-1p1 -modelhess -opt rsirfo_block_fsb -spin 1 -elec 0

-opt rsirfo_block_fsbは準ニュートン法であるRS-I-RFO法を構造最適化に使用することを示す。初期のへシアンに関しては、特にオプションで指定しない限り、単位行列が使われる。(以前のHessian更新法とは細かな点で異なる方法を使用している。具体的には、複数の座標変位や勾配変位を用いてHessianの更新を行う。)
-spinはスピン多重度の指定である。PySCFを使用するときは目的とするスピン多重度に1を引いた値を指定する。（デフォルトでは1が指定される。）
-elec Mは形式電荷をMとすることを示す。（デフォルトでは0が指定される。）
-ma yyy a bはyyykJ/molの活性化障壁を超えうるペア同士を近づける力を原子のラベル番号aとbのペアに構造最適化時に加えることを示す。
-os uma-s-1p1は今回使用するNNPを指定している。これを使用する際にASEライブラリが必要である。
-modelhessは、計算コストが非常に低い数式を使用して、近似したHessianを生成するオプションである。

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

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

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

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

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

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

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

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

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

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

NEB計算の結果の可視化

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

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

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

mobh35_15_Rh_traj_12.xyz

19
0 1
C      -0.642010995915     -1.022551789211      0.576716351690
C      -1.842778671957     -1.680398039351      0.682746623756
C      -2.862332494680     -1.179449401377     -0.117294264179
C      -2.439072826326     -0.083691574598     -0.858031535921
C      -1.133314552901      0.290521097091     -0.651497312572
Rh      0.770616631660      0.319822066132      0.029068639890
P       1.273724384959      2.497164402875      0.436915696567
P       2.688470948016     -0.860876418223     -0.276520054347
H       0.046901911789     -1.135333764938      1.433256301021
H      -1.978370279003     -2.453991969673      1.426146354460
H      -3.861408180897     -1.591706324923     -0.158300583229
H      -0.667582419861      0.874602095621     -1.466082987522
H      -3.027436177939      0.390734970733     -1.631492675555
H       2.870916465158     -1.333955798798     -1.587727626226
H       2.779392989101     -2.108302336872      0.369423216911
H       4.025119910172     -0.474386932580     -0.067271835704
H       2.489862200931      3.184850531899      0.269064724814
H       1.053406344993      2.900845264073      1.765302872021
H       0.455894812701      3.466103922119     -0.174421905875

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

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

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

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

python optmain.py .\mobh35_15_Rh_traj_12.xyz -spin 1 -elec 0 -os uma-s-1p1  -tcc -freq -opt rsirfo_block_bofill -order 1 -fc 5 -tr 0.1

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

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

(実行して得られた正確な遷移状態構造は計算開始時に、yyyy_mm_ddディレクトリ内に生成された新規ディレクトリ内のYYY_optimized.xyzとして保存されている。)

19
OptimizedStructure
C     -0.433003390728     -1.520537130479     -0.232050168075
C     -1.604444893995     -1.912065386639      0.290389206494
C     -2.647587705118     -0.924151720688      0.571695699089
C     -2.478757780942      0.252967441008     -0.044212592846
C     -1.275699772575      0.399605380175     -0.856310210550
Rh     0.511083750022      0.394789379859     -0.247576189469
P      1.359396975612      2.485350439414     -0.274404681148
P      2.449810761989     -0.613229173442      0.475060649133
H      0.219683444987     -2.300440206385     -0.621575306868
H     -1.830219034521     -2.965778172485      0.427293808469
H     -3.518018647466     -1.151638008494      1.172628729527
H     -1.451234984497      0.243782441304     -1.927887645073
H     -3.199121882646      1.062250413580     -0.013781765734
H      3.714905424701      0.006175234077      0.490294903895
H      2.851202117165     -1.784912119895     -0.189932554359
H      2.521537804065     -1.130407915012      1.781566570484
H      2.708785008509      2.753073662460      0.030799744386
H      0.799954740814      3.428994314428      0.604473332844
H      1.301728064624      3.276171127212     -1.436471530198

遷移状態構造

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

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

Mode                                 0                   1                   2
Freq [cm^-1]                     -179.3350             49.3622             60.6815
Reduced mass [au]                   3.5525              1.1048              2.5180
Force const [Dyne/A]               -0.0673              0.0016              0.0055
Char temp [K]                       0.0000             71.0212             87.3072
Normal mode                   x         y         z            x         y         z            x         y         z     
       C                 0.10014   -0.05497    0.08233   -0.00983    0.00131    0.02380   -0.02303    0.00203    0.07979
       C                 0.05888   -0.00056    0.02759   -0.01889    0.00390    0.00639   -0.04074    0.02297    0.05752
       C                 0.04277    0.00300   -0.04616   -0.02960    0.00231   -0.03062   -0.06266    0.02801   -0.04624
       C                 0.00665    0.02189   -0.02412   -0.01529   -0.00350   -0.03749   -0.02943   -0.00041   -0.09073
       C                -0.00254    0.16318    0.01045    0.00531   -0.01198   -0.00724    0.01812   -0.03823   -0.02392
      Rh                -0.02172   -0.01269   -0.00327    0.00144   -0.00294    0.00711    0.00487   -0.00585    0.02042
       P                -0.00917    0.01566   -0.00190    0.00303   -0.00490    0.01107   -0.01051   -0.00522    0.01717
       P                -0.00000   -0.02882   -0.00752    0.01754    0.01579   -0.01607    0.04536    0.01771   -0.06582
       H                 0.11548   -0.11108    0.21618    0.00312   -0.00251    0.05402   -0.00490   -0.01989    0.15423
       H                 0.02594    0.00866    0.04909   -0.01442    0.00488    0.02089   -0.03285    0.02889    0.11509
       H                -0.00480   -0.02730   -0.12427   -0.04442    0.00556   -0.05095   -0.09486    0.05074   -0.08458
       H                -0.00490    0.33838   -0.01741    0.03106   -0.01636   -0.01109    0.08012   -0.06626   -0.03018
       H                -0.07867   -0.05052   -0.06963   -0.01664   -0.00372   -0.06542   -0.03031    0.00171   -0.17054
       H                -0.00210   -0.02332   -0.01876    0.06787   -0.07853   -0.38184   -0.00575    0.11764    0.09374
       H                -0.00556   -0.03531    0.00149   -0.18885   -0.18255    0.21027    0.16089    0.18363   -0.28926
       H                -0.00372   -0.01196   -0.00037    0.20310    0.34278    0.10388    0.04915   -0.24309   -0.16910
       H                -0.00469   -0.00879   -0.00219    0.08820   -0.03166   -0.34080   -0.01350    0.00330    0.02383
       H                 0.00752    0.02783   -0.00443    0.24841   -0.12386    0.29550   -0.02299   -0.01880    0.02391
       H                 0.00350    0.01559   -0.00279   -0.34359    0.14595    0.13124   -0.01338    0.00441    0.02371
       
(...snip...)

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

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

終わりに#

　　自作ライブラリで、UMAモデルのニューラルネットワークポテンシャル(NNP, uma-s-1p1)を用いて、MOBH35データセットのNo. 15の反応のある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. Phys. Chem. A 2019, 123, 3761−3781 http://dx.doi.org/10.1021/acs.jpca.9b01546 (MOBH35のデータセットについて)

概要#