AIと量子計算で「夢の素材」を設計する——マテリアルズ・インフォマティクス完全実装ガイド
元記事: 新素材の発見に「200年」かかっていた時代が終わる——AIと量子計算が「夢の物質」を設計する
本記事では上記 note 記事のコンセプトをベースに、実際に手元で動かせる実装コードを追加しています。
対象読者
- 材料科学 × AI に興味がある研究者・エンジニア
- GNN（グラフニューラルネットワーク）を素材設計に応用したい方
- 量子化学シミュレーション（VQE）を試したい方
- Materials Project / pymatgen でデータ駆動型研究をはじめたい方
なぜ素材発見にAIが必要なのか
118元素 × 無限の結晶構造の組み合わせ空間を人手で探索すれば 200年 かかる。GNoME（Google DeepMind）はGNNで 220万種以上 の安定結晶構造を予測し、従来比10倍以上のペースで候補を提示した。
従来のアプローチ:
試行錯誤実験 → DFT計算 → 合成 → 評価（数年/素材）
マテリアルズ・インフォマティクス（MI）:
データベース → ML予測 → 逆設計 → DFT検証 → 合成（数週間/候補群）
システム全体アーキテクチャ
┌─────────────────────────────────────────────────────────────────┐
│ Materials Informatics Pipeline │
├──────────────┬──────────────┬──────────────┬───────────────────┤
│ Data Layer │ ML Layer │ Design Layer│ Validation Layer │
│ │ │ │ │
│ Materials │ GNN Property │ Inverse │ DFT (VASP/ │
│ Project API │ Predictor │ Design Gen │ Quantum ESPRESSO) │
│ │ │ │ │
│ ICSD/MatNavi │ Formation │ Diffusion- │ VQE (Qiskit) │
│ (NIMS) │ Energy / Eg │ based Gen │ Quantum Chem │
│ │ │ │ │
│ AFLOW / OQMD │ Stability │ Matlantis │ ASE Relaxation │
│ │ Classifier │ MLFF │ │
└──────────────┴──────────────┴──────────────┴───────────────────┘
↓ ↓ ↓
pymatgen PyG / DGL Crystal Diffusion VAE
環境構築
1. 基本パッケージのインストール
# Python 3.10+ 推奨
pip install pymatgen mp-api
# グラフニューラルネットワーク
pip install torch torch-geometric
pip install dgl dglgo -f https://data.dgl.ai/wheels/repo.html
# 量子化学シミュレーション
pip install qiskit qiskit-nature qiskit-aer
pip install pyscf # ab initio計算ライブラリ
# 原子シミュレーション
pip install ase # Atomic Simulation Environment
# 可視化・分析
pip install matminer # 特徴量エンジニアリング
pip install nglview # Jupyter上での結晶構造可視化
pip install matplotlib scikit-learn pandas numpy
# (オプション) NVIDIA GPU用 PyTorch Geometric 高速化
# pip install pyg-lib torch-scatter torch-sparse -f https://data.pyg.org/whl/torch-2.1.0+cu121.html
2. Materials Project API キーの取得
# https://materialsproject.org でアカウント作成後、API Keyを取得
export MP_API_KEY="your_api_key_here"
# または .env ファイルに記載
echo "MP_API_KEY=your_api_key_here" >> .env
Step 1: Materials Project からデータ取得
1-1. pymatgen で結晶構造データベースにアクセス
"""
materials_project_fetcher.py
Materials Project API から結晶データを取得・解析するモジュール
"""
import os
from mp_api.client import MPRester
import pandas as pd
import numpy as np
from pymatgen.core import Structure, Composition
from pymatgen.analysis.phase_diagram import PhaseDiagram, PDEntry
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
class MaterialsDataFetcher:
"""Materials Project APIからデータを取得・整形するクラス"""
def __init__(self, api_key: str = None):
self.api_key = api_key or os.environ.get("MP_API_KEY")
if not self.api_key:
raise ValueError("MP_API_KEY environment variable not set")
def fetch_materials_by_elements(
self,
elements: list[str],
max_sites: int = 20,
properties: list[str] = None
) -> pd.DataFrame:
"""
指定元素を含む材料を取得
Args:
elements: 元素リスト (例: ["Li", "Fe", "O"])
max_sites: 最大サイト数（構造の複雑さ制限）
properties: 取得するプロパティ一覧
Returns:
材料データのDataFrame
"""
default_properties = [
"material_id", "formula_pretty", "energy_above_hull",
"band_gap", "formation_energy_per_atom", "density",
"volume", "nsites", "spacegroup", "is_stable",
"theoretical", "structure"
]
props = properties or default_properties
with MPRester(self.api_key) as mpr:
docs = mpr.materials.summary.search(
elements=elements,
num_sites=(1, max_sites),
fields=props
)
records = []
for doc in docs:
record = {
"material_id": doc.material_id,
"formula": doc.formula_pretty,
"energy_above_hull": doc.energy_above_hull,
"band_gap": doc.band_gap,
"formation_energy": doc.formation_energy_per_atom,
"density": doc.density,
"nsites": doc.nsites,
"spacegroup_number": doc.spacegroup.number if doc.spacegroup else None,
"is_stable": doc.is_stable,
}
records.append(record)
df = pd.DataFrame(records)
print(f"取得した材料数: {len(df)}")
return df
def get_structure(self, material_id: str) -> Structure:
"""特定材料の結晶構造を取得"""
with MPRester(self.api_key) as mpr:
structure = mpr.get_structure_by_material_id(material_id)
return structure
def analyze_structure(self, structure: Structure) -> dict:
"""結晶構造を解析して特徴量を抽出"""
analyzer = SpacegroupAnalyzer(structure)
return {
"space_group_symbol": analyzer.get_space_group_symbol(),
"space_group_number": analyzer.get_space_group_number(),
"crystal_system": analyzer.get_crystal_system(),
"lattice_type": analyzer.get_lattice_type(),
"num_atoms": len(structure),
"volume": structure.volume,
"density": structure.density,
"a": structure.lattice.a,
"b": structure.lattice.b,
"c": structure.lattice.c,
"alpha": structure.lattice.alpha,
"beta": structure.lattice.beta,
"gamma": structure.lattice.gamma,
}
def fetch_battery_candidates(self) -> pd.DataFrame:
"""全固体電池候補材料（Li系酸化物・硫化物）を取得"""
print("全固体電池候補を検索中...")
# Li含有安定材料でバンドギャップが小さいもの
with MPRester(self.api_key) as mpr:
docs = mpr.materials.summary.search(
elements=["Li"],
energy_above_hull=(0, 0.05), # 安定性閾値 50meV/atom
fields=["material_id", "formula_pretty", "band_gap",
"formation_energy_per_atom", "energy_above_hull",
"is_stable", "density"]
)
records = [
{
"material_id": d.material_id,
"formula": d.formula_pretty,
"band_gap_eV": d.band_gap,
"formation_energy": d.formation_energy_per_atom,
"energy_above_hull": d.energy_above_hull,
"density": d.density,
"is_stable": d.is_stable,
}
for d in docs
]
df = pd.DataFrame(records)
print(f"候補材料: {len(df)} 件")
return df.sort_values("energy_above_hull")
# 使用例
if __name__ == "__main__":
fetcher = MaterialsDataFetcher()
# Li-Fe-O系の材料を検索
df = fetcher.fetch_materials_by_elements(["Li", "Fe", "O"])
print("\n=== Li-Fe-O系材料（バンドギャップ上位10件）===")
print(df.nlargest(10, "band_gap")[["formula", "band_gap", "energy_above_hull", "is_stable"]])
# 全固体電池候補
battery_df = fetcher.fetch_battery_candidates()
print("\n=== 全固体電池候補（安定性上位10件）===")
print(battery_df.head(10)[["formula", "band_gap_eV", "formation_energy", "density"]])
1-2. matminer で特徴量エンジニアリング
"""
feature_engineering.py
pymatgen 構造 → ML用特徴量ベクトルへの変換
"""
from matminer.featurizers.composition import (
ElementProperty, Stoichiometry, ValenceOrbital, IonProperty
)
from matminer.featurizers.structure import (
SiteStatsFingerprint, CoulombMatrix, RadialDistributionFunction
)
import pandas as pd
from pymatgen.core import Structure
class MaterialFeatureExtractor:
"""結晶構造から ML 用特徴量を抽出"""
def __init__(self):
# 組成ベース特徴量
self.comp_featurizers = [
ElementProperty.from_preset("magpie"), # Magpie特徴量（279次元）
Stoichiometry(),
ValenceOrbital(props=["frac"]),
]
# 構造ベース特徴量
self.struct_featurizers = [
SiteStatsFingerprint.from_preset("CrystalNNFingerprint_ops"),
]
def extract_composition_features(self, df: pd.DataFrame) -> pd.DataFrame:
"""
組成式から特徴量を抽出
df には 'composition' 列が必要（pymatgen.Composition オブジェクト）
"""
from pymatgen.core import Composition
if "composition" not in df.columns:
df["composition"] = df["formula"].apply(Composition)
for featurizer in self.comp_featurizers:
featurizer.set_n_jobs(4)
try:
df = featurizer.featurize_dataframe(df, "composition", ignore_errors=True)
print(f" {featurizer.__class__.__name__}: OK")
except Exception as e:
print(f" {featurizer.__class__.__name__}: SKIP ({e})")
return df
def extract_structure_features(self, structures: list[Structure]) -> pd.DataFrame:
"""結晶構造から構造的特徴量を抽出"""
fingerprints = []
featurizer = SiteStatsFingerprint.from_preset("CrystalNNFingerprint_ops")
for i, struct in enumerate(structures):
try:
fp = featurizer.featurize(struct)
fingerprints.append(fp)
if i % 10 == 0:
print(f" 構造 {i+1}/{len(structures)} 処理中...")
except Exception as e:
fingerprints.append([None] * featurizer.feature_labels().__len__())
print(f" 構造 {i+1} スキップ: {e}")
cols = featurizer.feature_labels()
return pd.DataFrame(fingerprints, columns=cols)
Step 2: GNN による結晶特性予測
2-1. 結晶構造をグラフに変換
"""
crystal_graph.py
結晶構造 → PyTorch Geometric グラフ変換
"""
import torch
import numpy as np
from torch_geometric.data import Data
from pymatgen.core import Structure
from pymatgen.analysis.graphs import StructureGraph
from pymatgen.analysis.local_env import CrystalNN
# 元素の原子番号 → one-hot または埋め込み用インデックス
ELEMENT_FEATURES = {
"H": 1, "Li": 3, "Be": 4, "B": 5, "C": 6, "N": 7, "O": 8, "F": 9,
"Na": 11, "Mg": 12, "Al": 13, "Si": 14, "P": 15, "S": 16, "Cl": 17,
"K": 19, "Ca": 20, "Ti": 22, "V": 23, "Cr": 24, "Mn": 25, "Fe": 26,
"Co": 27, "Ni": 28, "Cu": 29, "Zn": 30, "Ga": 31, "Ge": 32, "As": 33,
# ... 必要に応じて追加
}
NUM_ELEMENTS = 118
def structure_to_graph(structure: Structure, cutoff: float = 6.0) -> Data:
"""
pymatgen Structure → PyTorch Geometric Data オブジェクト
ノード特徴量: 原子番号のone-hot + 電気陰性度 + 原子半径 + ...
エッジ特徴量: 距離 + 方向ベクトル (Gaussian basis expansion)
Args:
structure: pymatgen の結晶構造
cutoff: 結合を考慮する最大距離 [Å]
Returns:
PyTorch Geometric の Data オブジェクト
"""
from pymatgen.core.periodic_table import Element
# ノード特徴量（原子ごと）
node_features = []
for site in structure:
elem = site.specie
z = elem.Z # 原子番号
# 基本特徴量: 原子番号, 電気陰性度, 共有結合半径, イオン化エネルギー
features = [
z / 118.0, # 正規化原子番号
float(elem.X) if elem.X else 0.0, # 電気陰性度 (Pauling scale)
float(elem.atomic_radius or 1.5) / 3.0, # 原子半径
float(elem.row) / 7.0, # 周期
float(elem.group) / 18.0, # 族
]
# one-hot encoding for atomic number (簡略版: 上位30元素)
z_onehot = [0.0] * 30
if z <= 30:
z_onehot[z - 1] = 1.0
features.extend(z_onehot)
node_features.append(features)
x = torch.tensor(node_features, dtype=torch.float)
# エッジ（原子間結合）の構築
edge_src, edge_dst, edge_attrs = [], [], []
all_neighbors = structure.get_all_neighbors(cutoff)
for i, neighbors in enumerate(all_neighbors):
for neighbor in neighbors:
j = neighbor.index
dist = neighbor.nn_distance
# Gaussian Basis Expansion で距離を多次元特徴量に
gauss_centers = np.linspace(0, cutoff, 64)
gauss_width = 0.5
gauss_features = np.exp(-((dist - gauss_centers) ** 2) / (2 * gauss_width ** 2))
edge_src.append(i)
edge_dst.append(j)
edge_attrs.append(gauss_features.tolist())
if len(edge_src) == 0:
# エッジなし（孤立原子）のフォールバック
edge_index = torch.zeros((2, 0), dtype=torch.long)
edge_attr = torch.zeros((0, 64), dtype=torch.float)
else:
edge_index = torch.tensor([edge_src, edge_dst], dtype=torch.long)
edge_attr = torch.tensor(edge_attrs, dtype=torch.float)
return Data(x=x, edge_index=edge_index, edge_attr=edge_attr)
2-2. GNN モデル定義（Crystal Graph Convolutional Neural Network）
"""
cgcnn_model.py
結晶グラフ畳み込みニューラルネットワーク（CGCNN）
バンドギャップ・形成エネルギーの予測に使用
"""
import torch
import torch.nn as nn
import torch.nn.functional as F
from torch_geometric.nn import (
CGConv, global_mean_pool, global_add_pool,
NNConv, GATConv
)
from torch_geometric.data import Data, DataLoader
import numpy as np
class CrystalGNN(nn.Module):
"""
Crystal Graph Neural Network
参考: Xie & Grossman (2018) "Crystal Graph Convolutional Neural Networks
for an Accurate and Interpretable Prediction of Material Properties"
"""
def __init__(
self,
node_features: int = 35, # ノード特徴量次元
edge_features: int = 64, # エッジ特徴量次元（Gaussian basis）
hidden_dim: int = 128,
num_conv_layers: int = 4,
output_dim: int = 1, # 予測値の次元（回帰: 1）
dropout: float = 0.1,
):
super().__init__()
# 初期埋め込み層
self.node_embedding = nn.Linear(node_features, hidden_dim)
# Crystal Graph Convolution 層（CGConv は PyG に実装済み）
self.conv_layers = nn.ModuleList([
CGConv(
channels=hidden_dim,
dim=edge_features,
batch_norm=True,
)
for _ in range(num_conv_layers)
])
# グローバルプーリング後の全結合層
self.fc_layers = nn.Sequential(
nn.Linear(hidden_dim, hidden_dim),
nn.BatchNorm1d(hidden_dim),
nn.SiLU(), # Swish activation
nn.Dropout(dropout),
nn.Linear(hidden_dim, hidden_dim // 2),
nn.SiLU(),
nn.Linear(hidden_dim // 2, output_dim),
)
def forward(self, data: Data) -> torch.Tensor:
x, edge_index, edge_attr = data.x, data.edge_index, data.edge_attr
# ノード特徴量を埋め込み空間へ
x = self.node_embedding(x)
x = F.silu(x)
# グラフ畳み込み
for conv in self.conv_layers:
x_new = conv(x, edge_index, edge_attr)
x = x + x_new # Residual connection
x = F.silu(x)
# グラフ全体を1ベクトルに集約（平均プーリング）
x = global_mean_pool(x, data.batch)
# 予測
out = self.fc_layers(x)
return out.squeeze(-1)
class MaterialsGNNTrainer:
"""GNNモデルの学習・評価を管理するクラス"""
def __init__(
self,
model: CrystalGNN,
target: str = "band_gap", # 予測ターゲット
device: str = "auto",
):
self.model = model
self.target = target
if device == "auto":
self.device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
else:
self.device = torch.device(device)
self.model.to(self.device)
print(f"使用デバイス: {self.device}")
print(f"モデルパラメータ数: {sum(p.numel() for p in model.parameters()):,}")
def train_epoch(self, loader: DataLoader, optimizer, criterion) -> float:
self.model.train()
total_loss = 0.0
for batch in loader:
batch = batch.to(self.device)
optimizer.zero_grad()
pred = self.model(batch)
target = batch.y.float()
loss = criterion(pred, target)
loss.backward()
# Gradient clipping for stability
torch.nn.utils.clip_grad_norm_(self.model.parameters(), max_norm=1.0)
optimizer.step()
total_loss += loss.item()
return total_loss / len(loader)
@torch.no_grad()
def evaluate(self, loader: DataLoader) -> dict:
self.model.eval()
preds, targets = [], []
for batch in loader:
batch = batch.to(self.device)
pred = self.model(batch)
preds.extend(pred.cpu().numpy())
targets.extend(batch.y.cpu().numpy())
preds = np.array(preds)
targets = np.array(targets)
mae = np.mean(np.abs(preds - targets))
rmse = np.sqrt(np.mean((preds - targets) ** 2))
return {"MAE": mae, "RMSE": rmse, "predictions": preds, "targets": targets}
def train(
self,
train_loader: DataLoader,
val_loader: DataLoader,
epochs: int = 100,
lr: float = 1e-3,
patience: int = 20,
) -> dict:
"""学習ループ（Early Stopping付き）"""
optimizer = torch.optim.Adam(self.model.parameters(), lr=lr)
scheduler = torch.optim.lr_scheduler.ReduceLROnPlateau(
optimizer, patience=10, factor=0.5, min_lr=1e-5
)
criterion = nn.MSELoss()
history = {"train_loss": [], "val_mae": [], "val_rmse": []}
best_val_mae = float("inf")
patience_counter = 0
for epoch in range(epochs):
train_loss = self.train_epoch(train_loader, optimizer, criterion)
val_metrics = self.evaluate(val_loader)
scheduler.step(val_metrics["MAE"])
history["train_loss"].append(train_loss)
history["val_mae"].append(val_metrics["MAE"])
history["val_rmse"].append(val_metrics["RMSE"])
if epoch % 10 == 0:
print(
f"Epoch {epoch:3d} | "
f"Train Loss: {train_loss:.4f} | "
f"Val MAE: {val_metrics['MAE']:.4f} eV | "
f"Val RMSE: {val_metrics['RMSE']:.4f} eV"
)
# Early Stopping
if val_metrics["MAE"] < best_val_mae:
best_val_mae = val_metrics["MAE"]
torch.save(self.model.state_dict(), "best_gnn_model.pth")
patience_counter = 0
else:
patience_counter += 1
if patience_counter >= patience:
print(f"Early Stopping at epoch {epoch}")
break
print(f"\n最良 Val MAE: {best_val_mae:.4f} eV")
return history
# 学習実行例
if __name__ == "__main__":
from torch_geometric.data import DataLoader
from crystal_graph import structure_to_graph
# ダミーデータで動作確認
import random
def create_dummy_dataset(n: int = 100):
"""動作確認用ダミーデータセット生成"""
from torch_geometric.data import Data
dataset = []
for _ in range(n):
n_atoms = random.randint(2, 10)
x = torch.randn(n_atoms, 35)
# ランダムエッジ
n_edges = n_atoms * 3
edge_index = torch.randint(0, n_atoms, (2, n_edges))
edge_attr = torch.randn(n_edges, 64)
y = torch.tensor([random.uniform(0, 5)]) # バンドギャップ (eV)
dataset.append(Data(x=x, edge_index=edge_index, edge_attr=edge_attr, y=y))
return dataset
dataset = create_dummy_dataset(200)
train_data = dataset[:160]
val_data = dataset[160:]
train_loader = DataLoader(train_data, batch_size=32, shuffle=True)
val_loader = DataLoader(val_data, batch_size=32)
model = CrystalGNN(node_features=35, edge_features=64, hidden_dim=128, num_conv_layers=4)
trainer = MaterialsGNNTrainer(model, target="band_gap")
history = trainer.train(train_loader, val_loader, epochs=50)
Step 3: 逆設計（Inverse Design）— 目標特性から構造を生成
"""
inverse_design.py
目標特性（バンドギャップ・形成エネルギー）から結晶構造を逆設計
Crystal Diffusion Variational AutoEncoder (CDVAE) の簡易実装
"""
import torch
import torch.nn as nn
import torch.nn.functional as F
import numpy as np
from typing import Optional
class CrystalVAE(nn.Module):
"""
結晶構造の変分オートエンコーダー
潜在空間から新しい結晶構造を生成
参考: CDVAE (Xie et al., 2021), DiffCSP (Jiao et al., 2023)
"""
def __init__(
self,
input_dim: int = 256, # 入力特徴量次元（GNNの出力）
latent_dim: int = 64, # 潜在空間次元
condition_dim: int = 4, # 条件（バンドギャップ, 形成E, etc.）
):
super().__init__()
self.latent_dim = latent_dim
# エンコーダー
self.encoder = nn.Sequential(
nn.Linear(input_dim + condition_dim, 512),
nn.LayerNorm(512),
nn.SiLU(),
nn.Linear(512, 256),
nn.SiLU(),
)
self.fc_mu = nn.Linear(256, latent_dim)
self.fc_logvar = nn.Linear(256, latent_dim)
# デコーダー（結晶構造パラメータを再構成）
self.decoder = nn.Sequential(
nn.Linear(latent_dim + condition_dim, 256),
nn.LayerNorm(256),
nn.SiLU(),
nn.Linear(256, 512),
nn.SiLU(),
nn.Linear(512, input_dim), # 構造特徴量の再構成
)
# 格子定数予測ヘッド (a, b, c, α, β, γ)
self.lattice_head = nn.Sequential(
nn.Linear(latent_dim, 64),
nn.SiLU(),
nn.Linear(64, 6),
nn.Softplus(), # 正値を保証
)
def encode(self, x: torch.Tensor, condition: torch.Tensor):
h = self.encoder(torch.cat([x, condition], dim=-1))
return self.fc_mu(h), self.fc_logvar(h)
def reparameterize(self, mu: torch.Tensor, logvar: torch.Tensor) -> torch.Tensor:
"""Reparameterization trick"""
std = torch.exp(0.5 * logvar)
eps = torch.randn_like(std)
return mu + eps * std
def decode(self, z: torch.Tensor, condition: torch.Tensor) -> dict:
zc = torch.cat([z, condition], dim=-1)
recon = self.decoder(zc)
lattice = self.lattice_head(z) # [a, b, c, α, β, γ]
return {"reconstruction": recon, "lattice": lattice}
def forward(self, x: torch.Tensor, condition: torch.Tensor) -> dict:
mu, logvar = self.encode(x, condition)
z = self.reparameterize(mu, logvar)
decoded = self.decode(z, condition)
return {**decoded, "mu": mu, "logvar": logvar, "z": z}
def vae_loss(self, recon: torch.Tensor, x: torch.Tensor,
mu: torch.Tensor, logvar: torch.Tensor,
beta: float = 1.0) -> torch.Tensor:
"""VAE損失 = 再構成誤差 + β * KLダイバージェンス"""
recon_loss = F.mse_loss(recon, x, reduction="sum")
kl_loss = -0.5 * torch.sum(1 + logvar - mu.pow(2) - logvar.exp())
return recon_loss + beta * kl_loss
class PropertyGuidedSampler:
"""
目標特性を指定して結晶構造を生成するサンプラー
Bayesian Optimization + VAE latent space exploration
"""
def __init__(self, vae: CrystalVAE, predictor: nn.Module):
"""
Args:
vae: 学習済み CrystalVAE
predictor: バンドギャップ等を予測する GNN
"""
self.vae = vae
self.predictor = predictor
self.device = next(vae.parameters()).device
def sample_with_target(
self,
target_band_gap: float,
target_formation_energy: float = -2.0,
n_samples: int = 100,
n_iterations: int = 500,
temperature: float = 1.0,
) -> list[dict]:
"""
目標バンドギャップを持つ構造を逆設計
Args:
target_band_gap: 目標バンドギャップ [eV]
target_formation_energy: 目標形成エネルギー [eV/atom]
n_samples: 生成する候補数
n_iterations: 最適化イテレーション数
temperature: サンプリング温度（高いほど多様性UP）
Returns:
候補構造のリスト（予測特性付き）
"""
self.vae.eval()
condition = torch.tensor(
[[target_band_gap, target_formation_energy, 0, 0]],
device=self.device, dtype=torch.float
).expand(n_samples, -1)
# ランダム初期化から最適化
z = torch.randn(n_samples, self.vae.latent_dim,
device=self.device, requires_grad=True)
optimizer = torch.optim.Adam([z], lr=0.01)
for i in range(n_iterations):
optimizer.zero_grad()
# デコード
decoded = self.vae.decode(z * temperature, condition)
lattice_params = decoded["lattice"]
# 目標特性との差を損失に（プレースホルダー）
# 実際は GNN predictor で予測して目標との差を最小化
lattice_loss = torch.sum(
(lattice_params[:, :3] - 5.0) ** 2 # 格子定数 ~5Å を目標
)
# 潜在空間の正則化
reg_loss = 0.01 * torch.sum(z ** 2)
loss = lattice_loss + reg_loss
loss.backward()
optimizer.step()
if i % 100 == 0:
print(f" Iter {i}: loss={loss.item():.4f}, "
f"avg_a={lattice_params[:, 0].mean().item():.2f}Å")
# 生成結果をまとめる
results = []
with torch.no_grad():
decoded = self.vae.decode(z, condition)
lattice = decoded["lattice"].cpu().numpy()
for idx in range(n_samples):
results.append({
"idx": idx,
"target_band_gap": target_band_gap,
"lattice_a": float(lattice[idx, 0]),
"lattice_b": float(lattice[idx, 1]),
"lattice_c": float(lattice[idx, 2]),
"lattice_alpha": float(lattice[idx, 3]),
"latent_vector": z[idx].detach().cpu().numpy().tolist(),
})
print(f"\n{n_samples} 件の候補構造を生成しました")
return results
# 逆設計の実行例
if __name__ == "__main__":
print("=== 逆設計デモ ===")
print("目標: バンドギャップ 2.0 eV の半導体材料")
vae = CrystalVAE(input_dim=128, latent_dim=64, condition_dim=4)
# ダミー予測器（実際は学習済みGNNを使用）
dummy_predictor = nn.Linear(128, 1)
sampler = PropertyGuidedSampler(vae, dummy_predictor)
candidates = sampler.sample_with_target(
target_band_gap=2.0,
n_samples=10,
n_iterations=200,
)
print("\n生成された候補構造（上位3件）:")
for c in candidates[:3]:
print(f" 候補 {c['idx']}: a={c['lattice_a']:.2f}Å, "
f"b={c['lattice_b']:.2f}Å, c={c['lattice_c']:.2f}Å")
Step 4: ASE で構造の最適化・DFTとの連携
"""
structure_optimizer.py
ASE（Atomic Simulation Environment）で構造緩和・DFT計算の前処理
"""
from ase import Atoms
from ase.optimize import BFGS, LBFGS
from ase.io import write, read
from ase.build import bulk, surface, add_adsorbate
import numpy as np
from pymatgen.core import Structure
from pymatgen.io.ase import AseAtomsAdaptor
class StructureOptimizer:
"""ASE を使った結晶構造の前処理・最適化クラス"""
def __init__(self):
self.adaptor = AseAtomsAdaptor()
def pymatgen_to_ase(self, structure: Structure) -> Atoms:
"""pymatgen Structure → ASE Atoms 変換"""
return self.adaptor.get_atoms(structure)
def ase_to_pymatgen(self, atoms: Atoms) -> Structure:
"""ASE Atoms → pymatgen Structure 変換"""
return self.adaptor.get_structure(atoms)
def create_battery_anode(self) -> Atoms:
"""
シリコンアノード構造の生成例
全固体電池の負極材料（Si は理論容量 3579 mAh/g）
"""
# ダイヤモンド構造のシリコン
si = bulk("Si", "diamond", a=5.43) # 格子定数 5.43 Å
si = si.repeat([2, 2, 2]) # 2×2×2 スーパーセル
print(f"Si アノード構造: {len(si)} 原子, 体積={si.get_volume():.1f} ų")
return si
def create_perovskite(self, A: str = "Ba", B: str = "Ti", X: str = "O") -> Atoms:
"""
ペロブスカイト構造 ABO3 を生成
太陽電池・強誘電体材料として重要
"""
from ase.build import bulk
# BaTiO3 ペロブスカイト（立方晶）
a = 4.01 # 格子定数 [Å]
atoms = Atoms(
symbols=[A, B, X, X, X],
positions=[
[0, 0, 0], # A サイト（Ba）
[a/2, a/2, a/2], # B サイト（Ti）
[a/2, a/2, 0], # O1
[a/2, 0, a/2], # O2
[0, a/2, a/2], # O3
],
cell=[a, a, a],
pbc=True,
)
print(f"ペロブスカイト {A}{B}{X}3: {len(atoms)} 原子")
return atoms
def create_catalyst_surface(
self,
element: str = "Pt",
facet: tuple = (1, 1, 1),
n_layers: int = 4,
) -> Atoms:
"""
触媒表面モデルの生成
燃料電池・CO2還元触媒の設計に使用
"""
slab = surface(element, facet, n_layers, vacuum=10.0)
slab = slab.repeat([3, 3, 1]) # 3×3 スーパーセル
# CO2 吸着分子を配置
add_adsorbate(slab, "C", height=2.0, position=(slab.cell[0, 0]/2, slab.cell[1, 1]/2))
print(f"{element}({facet[0]}{facet[1]}{facet[2]}) 触媒表面: {len(slab)} 原子")
return slab
def relax_structure_with_emt(self, atoms: Atoms, fmax: float = 0.05) -> Atoms:
"""
EMT（Effective Medium Theory）ポテンシャルで構造緩和
DFT の前処理・動作確認用（金属のみ対応）
本番では VASP/Quantum ESPRESSO/CP2K を使用する
"""
from ase.calculators.emt import EMT
atoms.calc = EMT()
# エネルギーと力の初期値
e_initial = atoms.get_potential_energy()
f_initial = np.max(np.abs(atoms.get_forces()))
print(f"初期エネルギー: {e_initial:.4f} eV")
print(f"初期最大力: {f_initial:.4f} eV/Å")
# BFGS 法で構造最適化
optimizer = BFGS(atoms, trajectory="relaxation.traj", logfile="relax.log")
optimizer.run(fmax=fmax)
e_final = atoms.get_potential_energy()
print(f"最終エネルギー: {e_final:.4f} eV")
print(f"エネルギー差: {e_final - e_initial:.4f} eV")
return atoms
def export_for_dft(self, atoms: Atoms, format: str = "vasp") -> str:
"""
DFT計算ソフト向けに構造ファイルを出力
対応形式:
- "vasp": VASP の POSCAR ファイル
- "espresso-in": Quantum ESPRESSO の入力ファイル
- "cif": CIF 形式（汎用）
"""
filename_map = {
"vasp": "POSCAR",
"espresso-in": "structure.in",
"cif": "structure.cif",
}
filename = filename_map.get(format, "structure.cif")
write(filename, atoms, format=format)
print(f"構造ファイルを出力: {filename}")
return filename
# 使用例
if __name__ == "__main__":
optimizer = StructureOptimizer()
# 全固体電池アノード（Si）の構造生成
si_anode = optimizer.create_battery_anode()
# ペロブスカイト太陽電池材料（BaTiO3）
batio3 = optimizer.create_perovskite("Ba", "Ti", "O")
# Pt(111) 触媒表面
pt_surface = optimizer.create_catalyst_surface("Pt", (1, 1, 1))
# DFT用ファイルに出力
optimizer.export_for_dft(batio3, format="cif")
# 簡易緩和（EMTポテンシャル、金属系のみ動作確認用）
from ase.build import bulk
au_bulk = bulk("Au", "fcc", a=4.08).repeat([2, 2, 2])
relaxed_au = optimizer.relax_structure_with_emt(au_bulk)
Step 5: VQE による量子化学シミュレーション
"""
vqe_materials.py
Variational Quantum Eigensolver (VQE) で電子構造を量子計算
小分子・活性サイトのエネルギーを精密計算
"""
from qiskit import QuantumCircuit
from qiskit.primitives import Estimator
from qiskit_nature.second_q.drivers import PySCFDriver
from qiskit_nature.second_q.mappers import JordanWignerMapper, ParityMapper
from qiskit_nature.second_q.circuit.library import (
UCCSD, HartreeFock, EvolvedOperatorAnsatz
)
from qiskit.algorithms.optimizers import COBYLA, L_BFGS_B, SPSA
from qiskit.algorithms.minimum_eigensolvers import VQE
import numpy as np
class QuantumChemistrySimulator:
"""
VQE を使った量子化学計算
触媒活性サイト・電解質分子のエネルギー精密計算に使用
"""
def __init__(self):
self.mapper = JordanWignerMapper() # フェルミオン → クービット変換
def compute_h2_ground_state(self, bond_length: float = 0.74) -> dict:
"""
水素分子 H2 の基底状態エネルギーを VQE で計算
最も基本的な量子化学ベンチマーク
Args:
bond_length: H-H 結合距離 [Å]
Returns:
{"energy": float, "hartree_fock_energy": float, "correlation_energy": float}
"""
print(f"\n=== H2 VQE 計算 (bond_length={bond_length}Å) ===")
# PySCF で分子軌道を計算
driver = PySCFDriver(
atom=f"H 0 0 0; H 0 0 {bond_length}",
basis="sto3g", # STO-3G 基底関数
)
problem = driver.run()
# フェルミオン演算子 → クービット演算子
fermionic_op = problem.hamiltonian.second_q_op()
qubit_op = self.mapper.map(fermionic_op)
print(f"量子ビット数: {qubit_op.num_qubits}")
print(f"パウリ項数: {len(qubit_op)}")
# Hartree-Fock 初期状態
hf_state = HartreeFock(
num_spatial_orbitals=problem.num_spatial_orbitals,
num_particles=problem.num_particles,
qubit_mapper=self.mapper,
)
# UCCSD アンザッツ（化学精度の VQE 標準）
ansatz = UCCSD(
num_spatial_orbitals=problem.num_spatial_orbitals,
num_particles=problem.num_particles,
qubit_mapper=self.mapper,
initial_state=hf_state,
)
print(f"回路深さ: {ansatz.decompose().depth()}")
print(f"パラメータ数: {ansatz.num_parameters}")
# VQE 実行
estimator = Estimator()
optimizer = COBYLA(maxiter=500, tol=1e-6)
vqe = VQE(estimator=estimator, ansatz=ansatz, optimizer=optimizer)
result = vqe.compute_minimum_eigenvalue(qubit_op)
# 核間反発エネルギーを加算して全エネルギーを求める
nuclear_repulsion = problem.nuclear_repulsion_energy
total_energy = result.eigenvalue.real + nuclear_repulsion
# ハートリー-フォック エネルギー（参照値）
hf_energy = problem.reference_energy
correlation_energy = total_energy - hf_energy
output = {
"bond_length_angstrom": bond_length,
"vqe_total_energy_hartree": total_energy,
"hartree_fock_energy_hartree": hf_energy,
"correlation_energy_hartree": correlation_energy,
"n_qubits": qubit_op.num_qubits,
"n_parameters": ansatz.num_parameters,
"optimizer_iterations": result.optimizer_evals,
}
print(f"\nVQE 全エネルギー: {total_energy:.6f} Hartree")
print(f"HF エネルギー: {hf_energy:.6f} Hartree")
print(f"相関エネルギー: {correlation_energy:.6f} Hartree")
return output
def potential_energy_scan(
self,
bond_lengths: list[float] = None,
) -> list[dict]:
"""
H2 分子の結合距離 vs エネルギーカーブを計算
解離エネルギーの評価・ポテンシャルエネルギー曲面の生成
"""
if bond_lengths is None:
bond_lengths = np.linspace(0.5, 2.5, 10).tolist()
print(f"ポテンシャルエネルギーカーブを計算中 ({len(bond_lengths)} 点)")
results = []
for bl in bond_lengths:
try:
result = self.compute_h2_ground_state(bond_length=bl)
results.append(result)
print(f" d={bl:.2f}Å: E={result['vqe_total_energy_hartree']:.6f} Ha")
except Exception as e:
print(f" d={bl:.2f}Å: 計算失敗 ({e})")
return results
def simulate_lithium_hydride(self) -> dict:
"""
LiH (リチウムハイドライド) の基底状態エネルギー計算
全固体電池の電解質研究に関連
"""
print("\n=== LiH VQE 計算 ===")
driver = PySCFDriver(
atom="Li 0 0 0; H 0 0 1.595", # LiH 平衡結合距離
basis="sto3g",
)
problem = driver.run()
fermionic_op = problem.hamiltonian.second_q_op()
# パリティマッパー（2量子ビット削減）を使用
parity_mapper = ParityMapper(num_particles=problem.num_particles)
qubit_op = parity_mapper.map(fermionic_op)
print(f"量子ビット数（パリティマッパー）: {qubit_op.num_qubits}")
# UCCSD アンザッツ
ansatz = UCCSD(
num_spatial_orbitals=problem.num_spatial_orbitals,
num_particles=problem.num_particles,
qubit_mapper=parity_mapper,
initial_state=HartreeFock(
num_spatial_orbitals=problem.num_spatial_orbitals,
num_particles=problem.num_particles,
qubit_mapper=parity_mapper,
),
)
estimator = Estimator()
vqe = VQE(
estimator=estimator,
ansatz=ansatz,
optimizer=L_BFGS_B(maxiter=1000),
)
result = vqe.compute_minimum_eigenvalue(qubit_op)
total_energy = result.eigenvalue.real + problem.nuclear_repulsion_energy
print(f"LiH 基底状態エネルギー: {total_energy:.6f} Hartree")
return {"molecule": "LiH", "ground_state_energy_hartree": total_energy}
# 量子化学デモ
if __name__ == "__main__":
simulator = QuantumChemistrySimulator()
# H2 の基底状態
h2_result = simulator.compute_h2_ground_state(bond_length=0.74)
# ポテンシャルエネルギーカーブ（5点）
pes = simulator.potential_energy_scan(
bond_lengths=[0.5, 0.7, 0.74, 1.0, 1.5]
)
print("\n=== ポテンシャルエネルギーカーブ ===")
for r in pes:
print(f" d={r['bond_length_angstrom']:.2f}Å: "
f"E={r['vqe_total_energy_hartree']:.6f} Ha")
Step 6: Matlantis API 連携（機械学習力場）
"""
matlantis_client.py
Preferred Networks の Matlantis（汎用機械学習力場）を使った高速構造緩和
DFT と同等精度で 1000倍以上高速
※ Matlantis の利用には有料サブスクリプションが必要
"""
import numpy as np
class MatlantisSimulator:
"""
Matlantis MLFF（Machine Learning Force Field）ラッパー
DFT比:
精度: ほぼ同等（MAE ~0.1 eV/Å）
速度: 1000-10000倍高速
対応元素: 72 元素
実際の使用には pfcc-extras パッケージとAPIキーが必要:
pip install pfcc-extras
export MATLANTIS_TOKEN=your_token
"""
def relax_structure(self, atoms) -> dict:
"""
構造緩和のデモ（Matlantis実際の呼び出し例）
実際のコード（APIキー取得後):
from pfcc_extras.structure.ase_interface import ASEInterface
from ase.optimize import FIRE
calc = ASEInterface(model_type="universal-fast") # または "universal-accurate"
atoms.calc = calc
optimizer = FIRE(atoms, logfile="matlantis_relax.log")
optimizer.run(fmax=0.03) # 最大力 0.03 eV/Å で収束
energy = atoms.get_potential_energy()
forces = atoms.get_forces()
stress = atoms.get_stress()
return {
"energy": energy,
"forces": forces,
"stress": stress,
"n_steps": optimizer.get_number_of_steps(),
}
"""
# デモ用のダミー返値
print("=== Matlantis 構造緩和（デモ） ===")
n_atoms = 10
return {
"energy": -50.0 + np.random.normal(0, 0.1),
"forces": np.random.randn(n_atoms, 3) * 0.01,
"max_force": 0.025,
"n_steps": 87,
"converged": True,
}
def compute_phonons(self, atoms, supercell_matrix=None) -> dict:
"""
フォノン分散計算（熱力学安定性の評価）
実際のコード:
from pfcc_extras.structure.ase_interface import ASEInterface
from phonopy import Phonopy
from phonopy.structure.atoms import PhonopyAtoms
calc = ASEInterface(model_type="universal-accurate")
# phonopy と連携してフォノン計算...
"""
print("フォノン分散計算中... (デモ)")
# 架空のフォノン分散（実際は Phonopy で計算）
q_points = np.linspace(0, 1, 50)
n_bands = 6 # 2原子セルのフォノン
frequencies = np.abs(np.random.randn(len(q_points), n_bands) * 5 + 10)
has_imaginary = np.any(frequencies < 0)
return {
"q_points": q_points.tolist(),
"frequencies_THz": frequencies.tolist(),
"has_imaginary_modes": has_imaginary,
"dynamically_stable": not has_imaginary,
}
def screen_candidates(
self,
structures: list,
target_property: str = "band_gap",
target_value: float = 2.0,
) -> list[dict]:
"""
大量候補構造のハイスループットスクリーニング
GNoME が生成した 220万構造を Matlantis でスクリーニングする流れ:
1. Matlantis で全構造を高速緩和（DFTより1000倍速い）
2. 安定性フィルタ（energy_above_hull < 50meV/atom）
3. 目標特性フィルタ（バンドギャップ・弾性率等）
4. フォノン計算で動的安定性確認
5. DFT で最終確認（上位候補のみ）
"""
print(f"\n=== ハイスループットスクリーニング ===")
print(f"候補数: {len(structures)} 件")
print(f"目標特性: {target_property} ≈ {target_value}")
results = []
for i, struct in enumerate(structures):
# デモ用ダミー計算
relax_result = self.relax_structure(struct)
# 架空の予測バンドギャップ
predicted_property = target_value + np.random.normal(0, 0.5)
results.append({
"idx": i,
"energy": relax_result["energy"],
"converged": relax_result["converged"],
f"predicted_{target_property}": predicted_property,
"property_error": abs(predicted_property - target_value),
"pass_filter": abs(predicted_property - target_value) < 0.3,
})
passed = [r for r in results if r["pass_filter"]]
print(f"フィルター通過: {len(passed)}/{len(results)} 件")
return sorted(results, key=lambda x: x["property_error"])
Step 7: 統合パイプライン
"""
mi_pipeline.py
マテリアルズ・インフォマティクス統合パイプライン
データ取得 → 特徴量抽出 → GNN学習 → 逆設計 → スクリーニング → DFT検証
"""
import json
import pandas as pd
import numpy as np
from pathlib import Path
class MaterialsInformaticsPipeline:
"""
MI の完全パイプライン統合管理
フェーズ:
1. データ収集（Materials Project / ICSD）
2. 特徴量エンジニアリング（matminer）
3. GNN 学習（CGCNN）
4. 逆設計（Crystal VAE）
5. 高速スクリーニング（Matlantis MLFF）
6. DFT 検証（VASP / Quantum ESPRESSO）
"""
def __init__(self, target_property: str = "band_gap", output_dir: str = "./mi_results"):
self.target_property = target_property
self.output_dir = Path(output_dir)
self.output_dir.mkdir(parents=True, exist_ok=True)
print(f"MI パイプライン初期化")
print(f" 目標特性: {target_property}")
print(f" 出力ディレクトリ: {self.output_dir}")
def phase1_data_collection(self, elements: list[str], n_max: int = 5000) -> pd.DataFrame:
"""フェーズ1: データ収集"""
print(f"\n{'='*50}")
print(f"フェーズ1: データ収集")
print(f" 対象元素: {elements}")
# 実際は MaterialsDataFetcher を使用
# fetcher = MaterialsDataFetcher()
# df = fetcher.fetch_materials_by_elements(elements)
# デモ用ダミーデータ
np.random.seed(42)
n = min(n_max, 1000)
df = pd.DataFrame({
"material_id": [f"mp-{i:05d}" for i in range(n)],
"formula": [f"{''.join(np.random.choice(elements, np.random.randint(2, 4)))}" for _ in range(n)],
"band_gap": np.abs(np.random.randn(n) * 2.0),
"formation_energy": -np.abs(np.random.randn(n) * 2.0),
"energy_above_hull": np.abs(np.random.randn(n) * 0.1),
"is_stable": np.random.random(n) > 0.3,
})
stable_df = df[df["is_stable"]]
print(f" 取得件数: {len(df)} 件")
print(f" 安定材料: {len(stable_df)} 件")
df.to_csv(self.output_dir / "raw_materials.csv", index=False)
return df
def phase2_featurize(self, df: pd.DataFrame) -> pd.DataFrame:
"""フェーズ2: 特徴量エンジニアリング"""
print(f"\n{'='*50}")
print(f"フェーズ2: 特徴量エンジニアリング")
# 実際は MaterialFeatureExtractor を使用
# extractor = MaterialFeatureExtractor()
# df = extractor.extract_composition_features(df)
# デモ: ダミー特徴量
n_features = 50
feature_cols = [f"feature_{i}" for i in range(n_features)]
feature_matrix = np.random.randn(len(df), n_features)
feature_df = pd.DataFrame(feature_matrix, columns=feature_cols)
df = pd.concat([df.reset_index(drop=True), feature_df], axis=1)
print(f" 特徴量次元: {n_features}")
df.to_csv(self.output_dir / "featurized_materials.csv", index=False)
return df
def phase3_train_gnn(self, df: pd.DataFrame) -> dict:
"""フェーズ3: GNN 学習"""
print(f"\n{'='*50}")
print(f"フェーズ3: GNN 学習 (target={self.target_property})")
# 実際は CrystalGNN + MaterialsGNNTrainer を使用
# trainer = MaterialsGNNTrainer(model, target=self.target_property)
# history = trainer.train(train_loader, val_loader)
# デモ: モック結果
metrics = {
"best_val_mae_eV": 0.12,
"best_val_rmse_eV": 0.18,
"epochs_trained": 87,
"model_path": str(self.output_dir / "best_gnn_model.pth"),
}
print(f" 最良 Val MAE: {metrics['best_val_mae_eV']:.3f} eV")
print(f" 最良 Val RMSE: {metrics['best_val_rmse_eV']:.3f} eV")
with open(self.output_dir / "gnn_metrics.json", "w") as f:
json.dump(metrics, f, indent=2)
return metrics
def phase4_inverse_design(
self,
target_value: float,
n_candidates: int = 1000,
) -> list[dict]:
"""フェーズ4: 逆設計"""
print(f"\n{'='*50}")
print(f"フェーズ4: 逆設計")
print(f" 目標 {self.target_property}: {target_value}")
# 実際は PropertyGuidedSampler を使用
# sampler = PropertyGuidedSampler(vae, predictor)
# candidates = sampler.sample_with_target(target_band_gap=target_value)
# デモ: ダミー候補
candidates = [
{
"idx": i,
"predicted_band_gap": target_value + np.random.normal(0, 0.3),
"predicted_formation_energy": -2.0 + np.random.normal(0, 0.5),
"lattice_a": 4.0 + np.random.normal(0, 0.5),
}
for i in range(n_candidates)
]
# 目標値に近い候補を選択
candidates.sort(key=lambda x: abs(x["predicted_band_gap"] - target_value))
top_candidates = candidates[:100] # 上位100件
print(f" 生成候補: {n_candidates} 件")
print(f" 目標値付近（±0.3eV）: {len([c for c in candidates if abs(c['predicted_band_gap'] - target_value) < 0.3])} 件")
return top_candidates
def phase5_screening(self, candidates: list[dict]) -> list[dict]:
"""フェーズ5: Matlantis でハイスループットスクリーニング"""
print(f"\n{'='*50}")
print(f"フェーズ5: ハイスループットスクリーニング")
print(f" 候補数: {len(candidates)} 件")
# 安定性・フォノン安定性でフィルタリング
screened = []
for c in candidates:
# Matlantis で緩和・安定性チェック（デモ）
energy_above_hull = np.abs(np.random.randn()) * 0.05
dynamically_stable = np.random.random() > 0.3
if energy_above_hull < 0.05 and dynamically_stable:
c["energy_above_hull"] = energy_above_hull
c["dynamically_stable"] = dynamically_stable
screened.append(c)
print(f" スクリーニング通過: {len(screened)}/{len(candidates)} 件")
return screened
def phase6_dft_validation(self, candidates: list[dict], top_n: int = 10) -> list[dict]:
"""フェーズ6: DFT最終検証"""
print(f"\n{'='*50}")
print(f"フェーズ6: DFT 最終検証（上位 {top_n} 件）")
# 実際は VASP/Quantum ESPRESSO ジョブを投入
# ASEInterface + VASPCalculator 等を使用
validated = []
for c in candidates[:top_n]:
# DFT 計算結果（デモ）
c["dft_band_gap"] = c["predicted_band_gap"] + np.random.normal(0, 0.1)
c["dft_formation_energy"] = c["predicted_formation_energy"] + np.random.normal(0, 0.05)
c["dft_validated"] = True
validated.append(c)
print(f" 候補 {c['idx']}: "
f"予測Eg={c['predicted_band_gap']:.2f}eV, "
f"DFT Eg={c['dft_band_gap']:.2f}eV")
return validated
def run_full_pipeline(
self,
elements: list[str] = ["Li", "Na", "O", "S", "P"],
target_value: float = 2.0,
) -> dict:
"""全フェーズを順次実行"""
print("\n" + "="*60)
print("マテリアルズ・インフォマティクス パイプライン 開始")
print("="*60)
df = self.phase1_data_collection(elements)
df = self.phase2_featurize(df)
metrics = self.phase3_train_gnn(df)
candidates = self.phase4_inverse_design(target_value, n_candidates=500)
screened = self.phase5_screening(candidates)
validated = self.phase6_dft_validation(screened, top_n=5)
summary = {
"total_database_size": len(df),
"gnn_val_mae": metrics["best_val_mae_eV"],
"candidates_generated": len(candidates),
"after_screening": len(screened),
"dft_validated": len(validated),
"top_candidates": validated,
}
print("\n" + "="*60)
print("=== パイプライン完了 ===")
print(f" データベース: {summary['total_database_size']} 件")
print(f" GNN MAE: {summary['gnn_val_mae']:.3f} eV")
print(f" 逆設計候補: {summary['candidates_generated']} 件")
print(f" スクリーニング後: {summary['after_screening']} 件")
print(f" DFT検証済み: {summary['dft_validated']} 件")
print("="*60)
with open(self.output_dir / "pipeline_summary.json", "w", encoding="utf-8") as f:
json.dump(summary, f, ensure_ascii=False, indent=2)
return summary
# パイプライン実行
if __name__ == "__main__":
pipeline = MaterialsInformaticsPipeline(
target_property="band_gap",
output_dir="./mi_battery_results"
)
# 全固体電池材料の探索
result = pipeline.run_full_pipeline(
elements=["Li", "Na", "O", "S", "P", "Si"],
target_value=2.5, # 目標バンドギャップ 2.5 eV（固体電解質向け）
)
実装フェーズロードマップ
| フェーズ | 期間 | 内容 | 成果物 |
|---|---|---|---|
| Phase 1 | 1〜2週 | pymatgen + Materials Project API でデータ収集・可視化 | 材料DB CSV, 可視化ダッシュボード |
| Phase 2 | 2〜4週 | matminer で特徴量エンジニアリング → ランダムフォレスト/XGBoost でベースライン構築 | ベースライン予測モデル（MAE < 0.5eV） |
| Phase 3 | 1〜2ヶ月 | CGCNN/MEGNET で GNN 実装 → Materials Project 全データで学習 | GNN モデル（MAE < 0.2eV） |
| Phase 4 | 2〜3ヶ月 | Crystal VAE / DiffCSP で逆設計パイプライン構築 | 新規候補構造の自動生成 |
| Phase 5 | 3〜6ヶ月 | Matlantis API でハイスループットスクリーニング | 上位候補リスト（安定性・目標特性） |
| Phase 6 | 継続 | DFT (VASP/QE) → 合成実験 → フィードバックループ | 論文・特許 |
主要ライブラリ比較
| ライブラリ | 用途 | 特徴 |
|---|---|---|
| pymatgen | 結晶構造操作・解析 | Materials Project 標準、豊富な解析機能 |
| matminer | 特徴量エンジニアリング | Magpie/SOAP/ElementProperty 等 280+ 特徴量 |
| PyTorch Geometric | GNN 実装 | CGConv/SchNet/DimeNet++ 内蔵 |
| DGL | GNN 実装（別系統） | HeteroGraph 対応、MPNN/AttentiveFP |
| ASE | 原子シミュレーション | DFT コード（VASP, QE, CP2K）との I/O |
| Qiskit Nature | 量子化学 | VQE/QEOM、PySCF/GAMESS ドライバー |
| Matlantis | ML力場（商用） | 72元素対応、DFT比 1000x 高速 |
| CHGNET | ML力場（OSS） | 磁気モーメント対応、無料 |
| ALIGNN | GNN（線グラフ） | 角度情報を含む高精度モデル |
まとめ
従来の材料探索（200年かかる）
└── 試行錯誤 → DFT → 合成 → 評価 → 繰り返し
マテリアルズ・インフォマティクス（AIで加速）
├── GNN でバンドギャップ・形成エネルギーを予測（MAE < 0.1eV）
├── 逆設計で目標特性から構造を自動生成
├── Matlantis で DFT 比 1000倍速スクリーニング
└── 量子コンピュータ（VQE）で難しい電子構造を精密計算
応用分野:
全固体電池（Li固体電解質） → Eg > 3eV、高イオン伝導率
室温超伝導体 → 高転移温度 Tc
太陽電池 → Eg ≈ 1.3eV（Shockley-Queisser 最適）
パワー半導体（GaN/Ga₂O₃） → 広バンドギャップ、高耐圧
CO₂還元触媒 → 吸着エネルギー最適化
参考リンク
- Materials Project — 結晶構造・特性データベース（無料）
- matminer ドキュメント — 特徴量エンジニアリング
- PyTorch Geometric — GNN ライブラリ
- CGCNN 論文 (2018) — Crystal Graph CNN
- GNoME 論文 (2023) — Google DeepMind 結晶生成
- Matlantis — 汎用 ML 力場（Preferred Networks）
- CDVAE 論文 (2022) — Crystal Diffusion VAE
- Qiskit Nature — 量子化学シミュレーション
- ASE ドキュメント — 原子シミュレーション環境
- NIMS MatNavi — 日本の材料データベース
元記事: 新素材の発見に「200年」かかっていた時代が終わる——AIと量子計算が「夢の物質」を設計する by @airbticsquantum