【機械学習ポテンシャル】分子構造最適化の環境構築

【機械学習ポテンシャル】分子構造最適化の環境構築

SHARE

機械学習ポテンシャル（MLポテンシャル）を使って、DFT 並みの精度を保ちつつ 100～10,000 倍高速に分子構造最適化を行う方法を解説する記事です。Google Colab 上で水分子やエタノールの最適化結果をエネルギーや結合長の収束グラフで可視化します。学ぶと、ブラウザだけで手軽に高精度・高速な分子シミュレーションを試せるようになり、試行錯誤のサイクルが飛躍的に短縮することができます！

動作検証済み環境

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
【PySCF】Pythonで始める量子化学計算

PySCFを使った無料で始められる量子化学計算の技術書を販売中

PySCFを使い構造最適化や振動数計算、物性解析など高度な量子化学計算ができます！

技術書ページへ

【初心者向け】Gaussianを使った量子化学計算入門

Gaussianを使った量子化学計算の初心者向け技術書を販売中

しばしば出くわすエラーへの対処法をはじめ
Gaussianと無料ソフトウェア Avogadro を組み合わせた物性解析手法が学べます！

技術書ページへ

機械学習 (ML) ポテンシャルとは？

MLポテンシャル (Machine-Learning Interatomic Potential) とは、

「量子化学計算（DFT / ab initio）が返すエネルギーと力を、ニューラルネットなどの統計モデルでほぼ同精度・桁違い高速に近似するもの」 です[1]。

以下でもう少し詳しく見ていきます。

  1. なぜ速いのか
    • DFT は電子波動関数を自洽的に解くため計算量が O(N³) 以上。
    • ML ポテンシャルは 原子座標 → ニューラルネットへの順伝播 だけなので O(N) に近い。
    • その結果、10²〜10⁴ 倍 速い評価が可能になります。
  2. 精度は大丈夫？
    • 学習データを高精度 DFT / CCSD(T) で網羅的にサンプリング。
    • 多くの有機分子での誤差 (MAE) は < 1 kcal mol⁻¹[2]。
    • ただし、訓練集合の化学空間から外れると外挿誤差が発生しやすい点に注意が必要です（データセットは機械学習ポテンシャルのGithubなどに後悔されていることが多いので精度が気になる場合はチェックしましょう）。

使用ツールの概要

本記事では主に以下のツール／ライブラリを使います。いずれも Python エコシステムで広く用いられています。Pythonで結果が得られるので後のデータ分析に対してすぐに接続できることも特徴ですね。

  1. TorchANI
    • PyTorch ベースで開発された ML ポテンシャルライブラリの一つ。
    • 事前学習済みの ANI 系列モデル（ANI-1x, ANI-2x）を簡単に利用可能。
    • 軽元素（H, C, N, O, F, S, Cl）向けに幅広くトレーニングされており、分子構造のエネルギーと力を高速に評価[3]。
  2. ASE (Atomic Simulation Environment)
    • Python で書かれた原子スケールのシミュレーション用フレームワーク[4]。
    • 特徴
      • Atoms クラス：原子システム（原子種・座標・セル情報など）を一元管理
      • Calculator インターフェース：DFT コードや ML ポテンシャルを「黒箱」として組み込める
      • 最適化 (Optimizer) 機能：BFGS, FIRE, MD, NEB など多数のアルゴリズムを提供
      • IO 機能：XYZ, CIF, POSCAR, VASP 形式などさまざまなファイル読み書きに対応
      • 可視化支援：トラジェクトリファイル (.traj) に結果を保存し、OVITO や ASE GUI で再生可能
      • ASEで動く様々な機械学習ポテンシャルが開発されている。
    • なぜ使うのか？
      • 分子構造や計算バックエンドを一貫した API で操作できるため、異なる計算コード（DFT も ML も）を切り替えて動かせる
      • Python スクリプト内で完結するので、ワークフローの自動化や可視化が容易

MLポテンシャルの実装方法

ここからは、Google Colab 上で TorchANI＋ASE を用い、水分子とエタノールの構造最適化を行うコードを順に解説します。まずはASEとTorchANIのパッケージをインストールしましょう！

# Colab で最初に一行だけ
!pip install -q ase torchani

1) モデルと Calculator の準備

import torchani
from ase import Atoms
from ase.optimize import BFGS
# ① 事前学習済みモデル（ANI-2x）をロード
model = torchani.models.ANI2x()
# ② ASE 流儀の Calculator に変換(必要な原子はここで定義してください)
ani_calc = torchani.ase.Calculator(['H', 'O', 'C'], model=model)

ここで使用する原子が不足していると後の計算でエラーになるので注意して下さい！

2) 水分子の構造最適化

# ③ 水分子の Atoms オブジェクト作成
water = Atoms('H2O', positions=[( 0.0000, 0.757, 0.586), ( 0.0000, -0.757, 0.586), ( 0.0000, 0.000, 0.000)], calculator=ani_calc)
# ④ 初期エネルギー計算
print(f'Initial Energy = {water.get_potential_energy():.6f} eV')
# ⑤ BFGS で構造最適化
opt = BFGS(water, trajectory='water_opt.traj')
opt.run(fmax=1e-3)
print(f'Optimized Energy = {water.get_potential_energy():.6f} eV')
print(f'Converged in {opt.nsteps} steps')

3) オープンソース構造（エタノール）の計算

それではエタノールの構造を計算します。構造情報が入っているXYZファイルはは IQmol リポジトリから取得します[5]。

from ase.io import read
import urllib.request
import io # ioモジュールをインポートします
from ase.optimize import BFGS # Import BFGS again as it's used later
# ⑥ GitHub 上の XYZ ファイルを読み込み
url = '<https://raw.githubusercontent.com/nutjunkie/IQmol/master/share/fragments/Molecules/Alcohols/Ethanol.xyz>'
# URLからデータをダウンロードします
with urllib.request.urlopen(url) as response: xyz_data = response.read().decode('utf-8')
# ダウンロードしたデータを文字列IOオブジェクトとして扱います
xyz_file_object = io.StringIO(xyz_data)
# Read the data from the file-like object
ethanol = read(xyz_file_object, format='xyz')
ethanol.set_calculator(ani_calc)
# ⑦ 初期エネルギー計算
print(f'Ethanol Initial Energy = {ethanol.get_potential_energy():.6f} eV')
# ⑧ 最適化
opt_eth = BFGS(ethanol, trajectory='ethanol_opt.traj')
opt_eth.run(fmax=1e-3)
print(f'Ethanol Optimized Energy = {ethanol.get_potential_energy():.6f} eV')
print(f'Converged in {opt_eth.nsteps} steps')

構造最適化収束までのエネルギーとO-H結合長の可視化

次に、上で生成した water_opt.traj に含まれる各最適化ステップの情報を使って、以下を可視化します。

  1. エネルギー収束：ステップごとに get_potential_energy() で得られるエネルギーの変化をプロット
  2. 幾何収束：例として O–H 結合長 (get_distance(2,0)) のステップ変化をプロット

これにより、最適化がどのように進行し、エネルギーと構造がどのように収束したかを直感的に把握できます。

1) Unicode マイナス & フォント設定

import matplotlib as mpl
# Unicode マイナスを通常のマイナスに置き換え
mpl.rcParams['axes.unicode_minus'] = False

2) プロットコード（英語軸名）

import matplotlib.pyplot as plt
from ase.io import read
# トラジェクトリ読み込み
traj = read('water_opt.traj', index=':')
# 各ステップのエネルギーとステップ番号
energies = [a.get_potential_energy() for a in traj]
steps = list(range(len(energies)))
# 各ステップの O–H 結合長 (O: idx=2, H: idx=0)
oh_dist = [a.get_distance(2, 0) for a in traj]
# Energy Convergence
plt.figure()
plt.plot(steps, energies, marker='o')
plt.xlabel('Step')
plt.ylabel('Energy (eV)')
plt.title('Energy Convergence')
plt.tight_layout()
plt.show()
# Geometry Convergence
plt.figure()
plt.plot(steps, oh_dist, marker='o')
plt.xlabel('Step')
plt.ylabel('O–H Bond Length (Å)')
plt.title('Geometry Convergence')
plt.tight_layout()
plt.show()

以上の図を見ると、構造最適化が順調に進んでいることが分かります。また、O-H結合長も安定したことが分かりました。

可視化コードの解説

ここでは上記のコードを解説します。

# Unicode マイナス設定
import matplotlib as mpl
mpl.rcParams['axes.unicode_minus'] = False
  • フォント未対応の Unicode “−” を通常の “-” に置き換え、マイナス記号が欠落しないようにします。
traj = read('water_opt.traj', index=':')
  • ASE の read で、.traj に保存された全ステップを読み込みます。
energies = [a.get_potential_energy() for a in traj]
oh_dist = [a.get_distance(2, 0) for a in traj]
  • get_potential_energy() と get_distance() をループで呼び出し、リスト化。
  • これを Matplotlib で折れ線グラフとしてプロットします。

最後に

本記事では、ML ポテンシャルの基礎から、Colab 上での環境構築 、水分子 & エタノールの最適化 、収束挙動の可視化 までを一貫して実践しました。いかがでしょうか？機械学習ポテンシャルは近年のマテリアルインフォマティクスを助ける非常に強力なツールなので是非この機会に学んでいきましょう!

参考文献

Q. Ge, D. B. Kaffes, M. O. De La Rosa-Padilla, and A. R. Canning, “IQmol: an open-source molecular editor and quantum chemistry viewer,” GitHub repository, https://github.com/nutjunkie/IQmol, accessed May 2025.

J. S. Smith, O. Isayev, and A. E. Roitberg, “ANI-1: an extensible neural network potential with DFT accuracy at force field computational cost,” Chemical Science, vol. 8, no. 4, pp. 3192–3203, 2017.

J. S. Smith, R. Zubatyuk, B. Nebgen, N. Lubbers, C. Barros, A. S. Roitberg, O. Isayev, and S. Tretiak, “The ANI-2x dataset of chalcogen-containing molecules for ∼kcal mol⁻¹ accuracy,” Journal of Chemical Theory and Computation, vol. 15, no. 11, pp. 8139–8147, 2019.

TorchANI Development Team, “TorchANI: PyTorch Implementation of ANI Neural Network Potentials,” GitHub repository, https://github.com/aiqm/torchani, accessed May 2025.

A. H. Larsen, J. J. Mortensen, J. Blomqvist, I. E. Castelli, R. Christensen, M. Dułak, J. Friis, M. N. Groves, B. Hammer, C. Hargus, E. D. Hermes, P. C. Jennings, P. B. Jensen, J. Kermode, J. R. Kitchin, E. Larsen, K. Le Page, C. J. Langreth, M. Lundqvist, D. Lu, T. Mao, A. Molina, and K. W. Jacobsen, “The Atomic Simulation Environment—a Python library for working with atoms,” Journal of Physics: Condensed Matter, vol. 29, no. 27, 273002, 2017.

【PySCF】Pythonで始める量子化学計算

PySCFを使った無料で始められる量子化学計算の技術書を販売中

PySCFを使い構造最適化や振動数計算、物性解析など高度な量子化学計算ができます！

技術書ページへ

【初心者向け】Gaussianを使った量子化学計算入門

Gaussianを使った量子化学計算の初心者向け技術書を販売中

しばしば出くわすエラーへの対処法をはじめ
Gaussianと無料ソフトウェア Avogadro を組み合わせた物性解析手法が学べます！

技術書ページへ

コメントを残す

メールアドレスが公開されることはありません。 ※ が付いている欄は必須項目です

前の記事

【RNA-seq】GSEAで遺伝子セットエンリッチメント解析を行う方法【R】
【RNA-seq】GSEAで遺伝子セットエンリッチメント解析を…

次の記事

Alphafold Analysisによるタンパク質の分析【タ…
【タンパク質デザイン】Alphafold Analysisによるタンパク質の分析