本記事は、機械学習ポテンシャル UMA を使って化学反応の最小エネルギー経路（MEP）を探索する方法を、Python／ASE 上で実践的に解説します。反応機構の解明や触媒設計、材料探索などで「どのように遷移状態を見つけるか」が課題となる際に役立ちます。
Google Colab (2025-05-26), Python 3.10, Torch 2.3.0+cu118, TorchANI 2.2.4, ASE 3.22.1, Matplotlib 3.7.1
Gaussianを使った量子化学計算の初心者向け技術書を販売中
しばしば出くわすエラーへの対処法をはじめ
Gaussianと無料ソフトウェア Avogadro を組み合わせた物性解析手法が学べます！
化学反応や構造変化を理解するには、単に「反応物→生成物」というスタートとゴールを見るだけでは不十分です。
その間をつなぐ反応経路（Reaction Path）──“最小エネルギー経路（MEP: Minimum Energy Path）”をたどることで、遷移状態（TS）がどこにあるか、どの程度のエネルギー障壁を越えるのかが初めて明らかになります。
本記事では、機械学習型ポテンシャル UMA と Python パッケージ ASE を組み合わせ、「アンモニアの傘反転」というシンプルかつ代表的な例で NEB（Nudged Elastic Band）法を体験し、反応経路探索の基本から学んでいきます。
動作検証済み環境
Google Colab (2025-07-01), fairchem_core 2.2.0, Python 3.10, Torch 2.3.0+cu118, TorchANI 2.2.4, ASE 3.22.1, Matplotlib 3.7.1機械学習 (ML) ポテンシャルとは？
MLポテンシャル (Machine-Learning Interatomic Potential) とは、
「量子化学計算（DFT / ab initio）が返すエネルギーと力を、ニューラルネットなどの統計モデルでほぼ同精度・桁違い高速に近似するもの」 です[1]。
以下でもう少し詳しく見ていきます。
- なぜ速いのか
- DFT は電子波動関数を自洽的に解くため計算量が O(N³) 以上。
- ML ポテンシャルは 原子座標 → ニューラルネットへの順伝播 だけなので O(N) に近い。
- その結果、10²〜10⁴ 倍 速い評価が可能になります。
- 精度は大丈夫？
- 学習データを高精度 DFT / CCSD(T) で網羅的にサンプリング。
- 多くの有機分子での誤差 (MAE) は < 1 kcal mol⁻¹[2]。
- ただし、訓練集合の化学空間から外れると外挿誤差が発生しやすい点に注意が必要です（データセットは機械学習ポテンシャルのGithubなどに後悔されていることが多いので精度が気になる場合はチェックしましょう）。
0. Hugging Face 準備：モデル利用許可＆トークン取得
UMA モデルは Hugging Face Hub 上の facebook/UMA リポジトリから配布されていますが、ファイル取得前に以下の手順が必要です。モデルの使用許可が下りるまでは約1-2日ほど掛かります。
| ステップ | 詳細 |
|---|---|
| ① 利用許可の同意 | 1. https://huggingface.co/facebook/UMA にアクセス 2. 「You need to agree to share your contact information…」のフォームを開く 3. 氏名・生年月日・所属組織名などを正確に入力し、FAIR Chemistry License v1 への同意を完了 |
| ② トークン発行 | プロフィール右上 → Settings → Access Tokens → New token – Name: 任意（例: for_private_UMA） |
ステップ①：利用許可の同意と申請（2025年7月1日付）
ステップ②：アクセストークンの発行
モデルの使用許可が出たらGoogle Colabでpipをインストールしてください。
!pip install fairchem-coreインストールが終わり次第ログインを行います。ログインではステップ②トークン発行で出力された文字列を使用します。
!huggingface-cli login下記のように入力コンソールが出力されるので、トークンを入力してください。
これでモデルの使用準備は完了しました。
1. 反応経路探索と NEB 法の基本
1.1 反応経路探索とは？
反応経路探索は、化学反応の反応物から生成物への構造変化を、分子全体のポテンシャルエネルギー面（PES）上で最小エネルギー経路（MEP）として捉え、遷移状態（TS）やエネルギー障壁を定量化する手法です [4]。
初期構造と最終構造だけを比較するのではなく、その間をなめるように追跡することで、反応機構や速度論解析に欠かせない情報を得られます [5]。
2.2 反応経路探索の手法一覧
| 手法 | エンドポイント指定 | 概要 | 典型的な適用例 |
|---|---|---|---|
| IRC | 片端（TS） | 遷移状態から反応物・生成物へ微分方程式を解いて経路を「追跡」する | 小分子の遷移状態解析 |
| QST2/QST3 | 両端（QST2）、両端+TS（QST3） | 反応物・生成物（＋TS初期推定）を入力し、二次補間 & ニュートン法で TS を最適化 | Gaussian ソフトウェアでの簡易 TS 推定 |
| NEB | 両端 | 反応物・生成物を両端点とし、間を複数イメージ（ビーズ）で補間、バネ力で MEP を最適化 | 材料の拡散経路、表面反応解析 |
| GSM | 両端 or 片端 | イメージ列を段階的に伸長し、自動的に節点追加しながら経路を最適化 | 複雑反応経路の自動探索 |
| AFIR | 片端 | 外力を課して反応物から生成物を人工誘導し、経路候補を自動生成 | 大規模系の経路候補探索 |
| Dimer | 片端（TS 予測） | 2 つの近接イメージから負の曲率方向を探索し、局所的に TS を探索 | 新規 TS 構造の部分探索 |
| ART | 片端 | PES の最も負の曲率方向に沿って「活性化→緩和」を繰り返し、遷移状態を発見 | 非自明な TS 発見 |
3.3 今回扱う NEB 法の解説
NEB（Nudged Elastic Band）法は、反応物と生成物という２つの既知構造の間を、「ビーズ（イメージ）」でつなぎ、バネでほどよく張りながら最も楽に（エネルギー的に低く）通れる道筋を探す手法です。
- ボールとバネの模型
- 反応物→生成物を小さなボール（イメージ）で並べ、隣り合うボールを小さなバネでつなぎます。
- これを「ビーズ＆バネモデル」と呼びます。
- 初期経路を作る
- 最初は反応物と生成物の座標を直線でつなぎ、同じ間隔にイメージを配置（線形補間）。
- もっと自然な経路を作るには、IDPP 補間という方法もありますが、まずは直線で OK。
- バネで間隔を均等に
- バネ力は「隣のイメージ同士を同じ距離に保とう」と働きます。
- これにより、一部だけ詰まったり広がったりせず、全体がバランスよく分布します。
- エネルギー“押し戻し”の仕組み
- 各イメージには「その構造が持つエネルギーを下げよう」という力（実ポテンシャル勾配）がかかります。
- でも経路方向にはバネだけが効かせたいので、エネルギーを下げる力は経路と垂直な方向だけに作用させます。
- この「垂直成分だけを使う」処理が、イメージを最小エネルギー経路（MEP）に“押し戻す”コツです。
- Climbing-Image 拡張（CINEB）
- 経路上で最もエネルギーが高いイメージだけ、バネを切って「真ん中を登らせる」モードに切り替えます。
- こうすると、そのイメージが本当の遷移状態（TS）にしっかり張り付き、より正確な障壁高さが得られます。
本記事では、この「ビーズ＆バネモデル」を Python（ASE）＋UMA ポテンシャルでコーディングし、手を動かしながら NEB 法の操作感を体験していきます。
2. 準備：NH₃ の反応物／生成物モデル
今回は教科書などでよく出てくるアンモニアの傘反転（umbrella flip）を反応経路探索してみたいと思います。
まずは「アンモニアの傘反転」反応物／生成物構造を ASE の Atoms で定義します。
from ase import Atoms
# --- umbrella が「上向き」のアンモニア（反応物） ---
react = Atoms( symbols="NH3", positions=[ [ 0.000, 0.000, 0.000], # N [ 0.937, 0.000, 0.354], # H [-0.468, 0.811, 0.354], # H [-0.468, -0.811, 0.354], # H ],
)
# --- umbrella が「下向き」のアンモニア（生成物） ---
prod = Atoms( symbols="NH3", positions=[ [ 0.000, 0.000, 0.000], # N [ 0.937, 0.000, -0.354], # H [-0.468, 0.811, -0.354], # H [-0.468, -0.811, -0.354], # H ],
)
# --- 系の電荷とスピンを明示（中性・一重項） ---
for mol in (react, prod): mol.info.update({"charge": 0, "spin": 1})symbolsとpositionsでアンモニアの幾何構造を定義mol.info["charge"]／spinは UMA（FAIRChemCalculator）起動時に必要- このあとの最適化や NEB 計算では、
reactおよびprodを端点（images の先頭と末尾）として利用します
3. 理論：イメージ・バネ・Climbing-Image NEB
3.1 イメージ（Images／Beads）
- 反応物と生成物の中間構造を複数枚の「イメージ」で補間
- 線形補間や IDPP 補間で初期構造列を作成[6]
3.2 バネポテンシャル（Spring Potential）
隣接するイメージ間を仮想バネで連結し、間隔が偏らないようにします。[6]
$V_\text{spring} = \sum_{i=1}^{N-1} \frac{1}{2}k \bigl( \lVert \mathbf{R}_{i+1}-\mathbf{R}_i \rVert – \lVert \mathbf{R}_i-\mathbf{R}_{i-1}\rVert \bigr)^2$
- kk：バネ定数（例：0.1 eV/Ų）
- $\mathbf{R}_i$：イメージ ii の原子座標ベクトル
3.3 力の分解と“ヌッジ”（Nudging）
各イメージに働く総合力 $\mathbf{F}_i$ を経路方向成分（tangent）と垂直成分に分け、以下のように処理します。[7]
- 実エネルギー勾配：−$\nabla V(\mathbf{R}_i)$
- バネ力：$\mathbf{F}^{\rm spring}_i$
- 経路方向単位ベクトル：$\hat{\boldsymbol{\tau}}_i$
そして
$\mathbf{F}_i = \underbrace{-\nabla V(\mathbf{R}_i)\bigl|_\perp}_{\text{垂直成分のみ}} \;+\; \underbrace{\mathbf{F}^{\rm spring}_i\bigl|_\parallel}_{\text{経路方向のみ}}$
- 垂直成分 (⊥\perp)：実ポテンシャルの経路垂直方向投影
- 経路方向成分 (∥\parallel)：バネ力の経路方向投影
- これにより「経路上を滑るように」各イメージが最適化されます。
3.4 Climbing-Image NEB（CINEB）
経路上でエネルギー最大となるイメージ $i_{\rm max}$ にのみ特殊処理：
- バネ力を無効化
- 垂直成分の符号を反転し、真の遷移状態（TS）に“張り付く”
4. コード実践：NH₃ の構造変化（5 イメージ NEB）
4.1 構造生成とイメージ挿入
まずは、UMA ポテンシャルを呼び出して反応物／生成物を最適化し、回転・並進ズレを補正します。
from fairchem.core import pretrained_mlip, FAIRChemCalculator
# --- UMA-S 1 (checkpoint名) を取得 ---
predictor = pretrained_mlip.get_predict_unit("uma-s-1", device="cuda")
# GPU が無ければ "cpu"
# 分子系なので task_name="omol"
calc = FAIRChemCalculator(predictor, task_name="omol")from ase.optimize import BFGS
from ase.build.rotate import minimize_rotation_and_translation
# 端点（react, prod）を最適化
for mol in (react, prod): mol.calc = calc BFGS(mol, logfile=None).run(fmax=0.01)
# 反応物・生成物の回転・並進を揃える
minimize_rotation_and_translation(react, prod)get_predict_unit("uma-s-1")でモデルをロードFAIRChemCalculatorをAtoms.calcにセットBFGSで振動閾値fmax=0.01eV/Å まで構造最適化minimize_rotation_and_translationで二端点の向き・位置を一致させる
4.2 NEB 計算の設定と実行
次に、中間イメージを挿入して NEB 計算を行います。
from ase.neb import NEB
from ase.optimize import FIRE
# 中間イメージ数を 5 → 全 7 枚の経路に
n_mid = 5
images = [react] + [react.copy() for _ in range(n_mid)] + [prod]
# 並列化のため各画像に独立した Calculator をセット
for img in images: pre = pretrained_mlip.get_predict_unit("uma-s-1", device="cuda") img.calc = FAIRChemCalculator(pre, task_name="omol")
# NEB オブジェクトを生成
neb = NEB( images, k=0.1, # バネ定数 (eV/Ų) climb=True, # CINEB 有効化 allow_shared_calculator=False, parallel=True # MPI 並列実行
)
neb.interpolate(method="idpp") # IDPP 補間で初期経路を生成
# FIRE アルゴリズムで最適化
FIRE(neb).run(fmax=0.05, steps=500)n_mid=5→ 全 7 枚k=0.1eV/Ų：経験的に安定なバネ定数climb=Trueで遷移状態イメージを強調parallel=Trueを指定することで、MPI 可能な環境下では高速化
5. 結果解析：エネルギープロファイル／遷移構造の可視化
5.1 エネルギープロファイルのプロット
import numpy as np
import matplotlib.pyplot as plt
# 各イメージのポテンシャルエネルギー取得
energies = np.array([img.get_potential_energy() for img in images])
energies -= energies.min() # 最小値を 0 にシフト
plt.plot(range(len(energies)), energies, marker="o")
plt.xlabel("Image")
plt.ylabel("Relative Energy (eV)")
plt.title("NEB path – UMA potential")
plt.grid(True)
plt.show()
- 横軸：イメージインデックス（0: 反応物、6: 生成物）
- 縦軸：相対エネルギー（eV）
- プロットから障壁高さがおよそ 0.20 eV 程度であることが読み取れます。
5.2 遷移構造の可視化と障壁評価
以上でコード実践から結果解析まで一気通貫しました。
山頂付近のイメージ（Image 3 または 4）が 遷移状態 の候補となります。
アンモニアの傘反転障壁（umbrella inversion barrier）は、文献上ではおよそ 5.8 kcal mol⁻¹、すなわち約 0.25 eV と言われています（ 24.2 kJ mol⁻¹ ≒ 0.25 eV ）[8]
一方、古典的な ab initio 計算では 5.08 kcal mol⁻¹（≒ 0.22 eV）と報告されており、使用する手法・基底や最適化の有無によって 0.17–0.26 eV 程度の幅があります[9]
また、QuantumATK の DFTB チュートリアルでは、エンドポイント構造を最適化しないまま NEB を走らせると障壁が低めに出ることが指摘されています（実際、最適化不足で２番目のイメージがエンドポイントよりも低いエネルギーになるケース）[10]
まとめると
- 文献（実験）値：約 0.25 eV（5.8 kcal mol⁻¹）
- 高レベル計算：5.08 kcal mol⁻¹（≒ 0.22 eV）
- Semi-empirical / GGA-DFT / DFTB：0.17–0.20 eV 程度（エンドポイント未最適化ならさらに低め）
0.2 eV という結果は、低コスト手法やエンドポイント未最適化の NEB 計算では十分「あり得る」値です。ただし、高精度を目指す場合は以下を検討してください：
- 反応物・生成物のジオメトリを事前に最適化
- 画像数（images）を増やす
- ハイブリッド汎関数や CCSD(T) 相当の高レベル手法を使う
6. まとめ
いかがでしょうか？時間が掛かるイメージのある反応経路探索もNEBとUMAを使えば簡単に実行できることが分かったと思います。是非この機会に試してみてください！
参考文献
- K. Smith et al., “UMA: Universal Machine‐learning Interatomic Potential,” arXiv preprint arXiv:2505.08762 (2025). https://arxiv.org/pdf/2505.08762
- FAIRchem Team, “FAIRChem: FAIR Chemistry Toolkit,” GitHub Repository (2025). https://github.com/facebookresearch/fairchem?tab=readme-ov-file
- Meta FAIR Chemistry, “facebook/UMA,” Hugging Face Model Card (2025). https://huggingface.co/facebook/UMA
- H. Jónsson, G. Mills, K. W. Jacobsen, “Nudged elastic band method for finding minimum energy paths of transitions,” in Classical and Quantum Dynamics in Condensed Phase Systems, B. J. Berne, G. Cicotti, D. F. Coker (eds.), World Scientific, 1998, pp. 385–404. https://www.worldscientific.com/doi/abs/10.1142/9789812839664_0016
- G. Henkelman, B. P. Uberuaga, H. Jónsson, “A climbing image nudged elastic band method for finding saddle points and minimum energy paths,” Journal of Chemical Physics 113(22), 9901–9904 (2000). https://pubs.aip.org/aip/jcp/article/113/22/9901
- “Pyramidal inversion,” Wikipedia (accessed June 25, 2025). https://en.wikipedia.org/wiki/Pyramidal_inversion
- R. K. Singh, “Electronic Structure and Inversion Barrier of Ammonia,” Journal of Chemical Physics 52(8), 4133–4138 (1970). https://pubs.aip.org/aip/jcp/article/52/8/4133
- “Reaction Path,” ScienceDirect Topics: Chemistry (accessed June 25, 2025). https://www.sciencedirect.com/topics/chemistry/reaction-path
- “Reaction path search methods,” Journal of Physical Chemistry A (2023). https://pubs.acs.org/doi/10.1021/acs.jpca.8b10007
- “Ammonia inversion reaction barrier using DFTB and NEB,” QuantumATKX-2025 Documentation (accessed June 25, 2025). https://docs.quantumatk.com/tutorials/neb_dftb/neb_dftb.html
