MatterSim応用編：MatterSimとPythonでマテリアルズインフォマティクス(MI)に挑戦！

Posted at 2025-04-30

近年、マテリアルズインフォマティクス（MI）の重要性が急速に高まっています。
MI とは、実験やシミュレーションで得られる多様なデータを活用し、データ駆動型で新規材料を探索するアプローチです。

本記事では Python と MatterSim を組み合わせ、実際に小規模な MI タスクへ取り組む手順を解説します。
シミュレーションデータを自動生成し、そこから機械学習モデルを構築する一連の流れを体験いただけます。

補足
MatterSim のインストール方法や基本操作は、下記の記事にまとめています。必要に応じてご参照ください。
MatterSimを使って材料物性予測！物性予測・フォノン計算までやってみた【Python/ASE】

2. MI で MatterSim が活躍する具体課題：高エントロピー合金（HEA）設計

本稿で扱うテーマは 高エントロピー合金（HEA）の固溶体設計 です。

HEA とは
5 元素以上を等モル比で混合した新規合金であり、優れた機械特性・耐食性を示すことから次世代材料として注目されています。
課題
組成空間が膨大なため、実験のみで網羅的に探索するのは現実的ではありません。
MI と高速シミュレーションを組み合わせ、有望な組成を効率よく絞り込む方法が求められています。

近年報告された MI × HEA の代表例は次のとおりです。

耐食 HEA の機械学習探索
Zeng et al. は物理特徴量と機械学習を組み合わせ、AlCrFeCoNi 系 HEA の耐食性を高精度に予測しました

Zeng, C. et al. (2023) Machine learning accelerated discovery of corrosion-resistant high-entropy alloys. arXiv:2307.06384
Transformer による HEA 特性予測
Kamnis et al. は Transformer ベースの言語モデルで元素間相互作用を捉え、従来手法を上回る性能を達成しました

Kamnis, S. & Delibasis, K. (2024) High-entropy alloy property predictions using a Transformer-based language model. arXiv:2411.04861

これらの研究は、データ駆動型手法が HEA 探索に有効であることを示しています。
本記事でも MatterSim により物性データを作成し、機械学習モデルで 固溶体になりやすい組成 を直接予測する流れを示します。

3. HEA データセットの作成

3.1 方針

格子：FCC を仮定
元素候補：Ni, Co, Cr, Fe, Mn, Cu, Al, Ti, V, Mo
生成数：200 構造
予測項目：エネルギー（eV/atom）を中心に、体積・密度も取得

3.2 構造生成

Ni-FCC（格子定数 3.52 Å）を 5×5×5 に拡張し、毎回ランダムに選んだ 5 元素で全原子を置換します。
これにより多様な組成・原子配置を持つ 200 構造を作成し、hea_structures フォルダに保存します。
10 構造ごとに進捗を表示するため、大量生成でも処理状況を把握できます。

import random
import numpy as np
from ase.build import bulk, make_supercell
from ase import Atoms
from ase.io import write
import os

# 出力ディレクトリ作成
output_dir = "hea_structures"
os.makedirs(output_dir, exist_ok=True)

# 元素リスト（10種類程度）
elements_master = ["Ni", "Co", "Cr", "Fe", "Mn", "Cu", "Al", "Ti", "V", "Mo"]

# 元構造（NiのFCC構造）を作成
base = bulk('Ni', 'fcc', a=3.52)  # 格子定数 a=3.52 Å（NiのFCC標準値）
supercell = base.repeat((5, 5, 5))

# 置換して構造を作成
n_structures = 200
for i in range(n_structures):
    atoms = supercell.copy()

    # 毎回ランダムに5元素を選択
    selected_elements = random.sample(elements_master, 5)

    for atom in atoms:
        atom.symbol = random.choice(selected_elements)  # 5元素からランダムに置換

    filename = os.path.join(output_dir, f"hea_{i:03d}.xyz")
    write(filename, atoms)

    # 進捗表示（10個ごと）
    if (i + 1) % 10 == 0:
        print(f"{i + 1}/{n_structures} 個生成完了")

print(f"全ての{n_structures}個の擬似HEA構造を{output_dir}フォルダに保存しました。")

3.3 構造最適化と物性予測

生成した各構造に対し、MatterSim を ASE の Calculator として設定し BFGS で構造最適化を行います。
最適化後の

エネルギー（eV/atom）
体積（Å³）
密度（g/cm³）

を取得し、hea_dataset_relaxed.csv にまとめます。

import glob, os, torch, pandas as pd
from ase.io import read, write
from ase.optimize import BFGS
from mattersim.forcefield import MatterSimCalculator  # 公式 ASE-Calculator

# 事前設定
device = "cuda" if torch.cuda.is_available() else "cpu"
calc   = MatterSimCalculator(device=device,
                             load_path="MatterSim-v1.0.0-1M.pth")  # 適宜パスを調整

# 入出力ディレクトリ
in_dir  = "hea_structures"
out_dir = "hea_structures_relaxed"

# 結果フォルダが無ければ作成
os.makedirs(out_dir, exist_ok=True)

# 構造ファイルを昇順（hea_000.xyz → hea_001.xyz → …）で取得
files = sorted(glob.glob(os.path.join(in_dir, "hea_*.xyz")))

#メインループ
data = []
for i, f in enumerate(files):
    atoms = read(f)
    atoms.calc = calc

    # 構造最適化（BFGS、収束条件 fmax=0.05 eV/Å）
    BFGS(atoms).run(fmax=0.05)

    # 物性値の取得
    e_pa  = atoms.get_potential_energy() / len(atoms)        # eV/atom
    vol   = atoms.get_volume()                               # Å³
    dens  = atoms.get_masses().sum() / (vol * 1.66054)       # g/cm³
    comp  = atoms.get_chemical_formula(mode="hill")

    data.append({"file": os.path.basename(f),
                 "formula": comp,
                 "E (eV/atom)": e_pa,
                 "V (Å³)": vol,
                 "ρ (g/cm³)": dens})

    # 最適化後構造を保存
    out_path = os.path.join(out_dir, os.path.basename(f))
    write(out_path, atoms)

    # 進捗表示（10構造ごと）
    if (i + 1) % 10 == 0:
        print(f"{i + 1}/{len(files)} 構造を最適化完了")

# 結果の出力
pd.DataFrame(data).to_csv("hea_dataset_relaxed.csv", index=False)
print("→ hea_dataset_relaxed.csv を生成しました。")

10分ほど時間がかかると思いますが、200構造に対するエネルギーなどの物性値が求められると思います。

3.4 純元素基準エネルギー

ΔE_f（形成エネルギー）を計算するため、上記 10 元素の純元素結晶エネルギーを同様に評価し、element_ref.csv として保存します。
結晶相は簡便に FCC とし、単位セルを 3×3×3 に拡張してリラックス後のエネルギーを用います。

import os, torch, pandas as pd
from ase.build import bulk
from ase.optimize import BFGS
from mattersim.forcefield import MatterSimCalculator

device = "cuda" if torch.cuda.is_available() else "cpu"
calc   = MatterSimCalculator(device=device,
                             load_path="MatterSim-v1.0.0-1M.pth")

elements = ["Ni", "Co", "Cr", "Fe", "Mn", "Cu", "Al", "Ti", "V", "Mo"]
rows = []

for el in elements:
    atoms = bulk(el, "fcc", a=3.52).repeat((3, 3, 3))
    atoms.calc = calc
    BFGS(atoms).run(fmax=0.05)
    e_pa = atoms.get_potential_energy() / len(atoms)
    rows.append({"element": el, "E_ref (eV/atom)": e_pa})
    print(f"{el}: {e_pa:.3f} eV/atom")

pd.DataFrame(rows).to_csv("element_ref.csv", index=False)
print("→ element_ref.csv を保存しました。")

4. ΔE_f を用いた機械学習モデル

4.1 ΔE_f の算出

合金 200 件のエネルギーと、純元素エネルギーの濃度加重平均との差から ΔE_f を求めました。
ΔE_f < 0 なら固溶体が熱力学的に安定、ΔE_f > 0 なら相分離傾向を示唆します。

import pandas as pd
from pymatgen.core import Composition

# --- データ読み込み ---------------------------------------------------
df_alloy = pd.read_csv("hea_dataset_relaxed.csv")
df_ref   = pd.read_csv("element_ref.csv")

# 参照辞書 { 'Co': -3.99, 'Ni': -4.11, ... }
ref_dict = df_ref.set_index("element")["E_ref (eV/atom)"].to_dict()

# --- ΔE_f 計算関数 ----------------------------------------------------
def formation_energy(row):
    comp = Composition(row["formula"])
    e_mix = row["E (eV/atom)"]
    e_ref_weighted = sum(
        ref_dict[el.symbol] * comp.get_atomic_fraction(el)
        for el in comp
    )
    return e_mix - e_ref_weighted  # 負なら混合が安定

# ΔE_f 列を追加
df_alloy["ΔE_f (eV/atom)"] = df_alloy.apply(formation_energy, axis=1)
df_alloy.to_csv("hea_dataset_with_dEf.csv", index=False)

print("→ ΔE_f を付加したデータセットを保存しました")

4.2 ランダムフォレスト回帰

特徴量：10 元素の組成フラクション（10 次元）
学習／テスト：8:2 に分割（160 件／40 件）
モデル：ランダムフォレスト回帰（決定木 300 本）

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from pymatgen.core import Composition
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_absolute_error, r2_score

# データ読み込み
df_alloy = pd.read_csv("hea_dataset_relaxed.csv")
df_ref   = pd.read_csv("element_ref.csv")

# 純元素エネルギー辞書 { 'Ni': -4.11, ... }
ref_dict = df_ref.set_index("element")["E_ref (eV/atom)"].to_dict()

# ΔE_f を計算して列を追加
def formation_energy(row):
    comp = Composition(row["formula"])
    e_mix = row["E (eV/atom)"]
    e_ref = sum(ref_dict[el.symbol] * comp.get_atomic_fraction(el) for el in comp)
    return e_mix - e_ref

df_alloy["ΔE_f (eV/atom)"] = df_alloy.apply(formation_energy, axis=1)
df_alloy.to_csv("hea_dataset_with_dEf.csv", index=False)
print("ΔE_f を計算し hea_dataset_with_dEf.csv に保存しました")

# 元素フラクションを特徴量に変換
elem_list = ["Ni", "Co", "Cr", "Fe", "Mn", "Cu", "Al", "Ti", "V", "Mo"]
def comp_to_vec(formula):
    comp = Composition(formula)
    return np.array([comp.get_atomic_fraction(el) if el in comp else 0.0
                     for el in elem_list])

X = np.vstack(df_alloy["formula"].apply(comp_to_vec).values)
y = df_alloy["ΔE_f (eV/atom)"].values

# データ分割
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.2, random_state=42)

# ランダムフォレスト回帰
rf = RandomForestRegressor(n_estimators=300, random_state=42, n_jobs=-1)
rf.fit(X_tr, y_tr)
y_pred = rf.predict(X_te)

mae = mean_absolute_error(y_te, y_pred)
r2  = r2_score(y_te, y_pred)
print(f"テスト MAE : {mae:.3f} eV/atom")
print(f"テスト R²  : {r2:.3f}")

# 散布図（真値 vs 予測）
plt.figure(figsize=(4,4))
plt.scatter(y_te, y_pred, alpha=0.7)
plt.plot([y_te.min(), y_te.max()], [y_te.min(), y_te.max()], ls="--")
plt.xlabel("True ΔE_f (eV/atom)")
plt.ylabel("Predicted ΔE_f (eV/atom)")
plt.title("ΔE_f prediction")
plt.tight_layout()
plt.savefig("deltaEf_pred_vs_true.png", dpi=300)
print("散布図を deltaEf_pred_vs_true.png に保存しました")

# 特徴量重要度バーグラフ
importances = rf.feature_importances_
order = np.argsort(importances)[::-1]

plt.figure(figsize=(6,3))
plt.bar(range(len(elem_list)), importances[order])
plt.xticks(range(len(elem_list)), np.array(elem_list)[order])
plt.ylabel("Feature importance")
plt.title("Element fraction importance")
plt.tight_layout()
plt.savefig("feature_importance.png", dpi=300)
print("特徴量重要度グラフを feature_importance.png に保存しました")

手元の環境で実施すると以下のようになります。

ΔE_f を計算し hea_dataset_with_dEf.csv に保存しました
テスト MAE : 0.046 eV/atom
テスト R²  : 0.883
散布図を deltaEf_pred_vs_true.png に保存しました
特徴量重要度グラフを feature_importance.png に保存しました

実際に生成されたグラフを確認してみましょう。まずはエネルギーの予測値です。R² が0.88なので散布図（真値 vs 予測）は対角線上に良く沿い、モデルが ΔE_f を高精度で再現できていることが確認できました。

また、特徴量重要度は以下のようになりました。

1. Mo : 0.562
2. Al : 0.125
3. Ti : 0.094
4. Cu : 0.092
5. Cr : 0.033

4.3 固溶体になりやすい組成の判定

文献値より ΔE_f < –0.05 eV/atom で単相形成が期待されます。
モデル誤差（MAE 0.046 eV/atom）を考慮し、ΔE_f_pred < –0.07 eV/atom を固溶体有望域と設定しました。

Mo / Ti を多く含み、Cu / Al を抑えた組成 がこの閾値を満たしやすく、逆に Mo / Ti を欠き Cu / Al が高濃度の組成は不安定域に入る傾向が見られました。

5. まとめと今後の展望

ΔE_f 回帰モデル（MAE 0.046 eV/atom、R² 0.883） により、組成フラクションだけから固溶体安定性を迅速に予測できることを示しました。
モデルは Mo・Ti による安定化、Cu・Al による不安定化 を定量的に捉えています。
実務では
- ΔE_f_pred < –0.07 eV/atom を候補抽出し詳細計算・実験へ
- ΔE_f_pred > +0.05 eV/atom を一次スクリーニングで除外
  することで、計算／実験コストを大幅に削減できます。

今後の拡張

体積・密度・原子サイズ差などの物理特徴量を追加し、モデル精度を向上
不確実性推定付きアクティブラーニングでサンプル効率を最大化
相分離／単相分類モデルの導入により、固溶体判定をさらに明確化

おわりに

MatterSim は第一原理計算に比べて軽量でありながら物性を高精度に予測できるため、MI や AI を活用した材料開発における 高速スクリーニング や 初期探索 に最適です。
今回示したワークフローをベースに、より大規模な組成空間や別物性へ拡張し、固溶体 HEA の探索をさらに加速させてみてください。

参考資料

🙌 もし参考になったら「いいね」や「ストック」お願いします！
X（Twitter）: @ap123hao （フォロー歓迎です）

You get articles that match your needs
You can efficiently read back useful information
You can use dark theme

What you can do with signing up