概要
Microsoft 社が開発した機械学習汎関数 Skala の使用方法
Skala は機械学習ポテンシャルではなく、DFT 計算用のディープラーニングベースの交換相関汎関数。
$O(N^3$) と書かれているので、機械学習ポテンシャルのように一瞬で計算が完了するわけではない。
PySCF を通して動かすことができる。
Azure AI Foundry カタログと github を通して公開されている。
github：https://github.com/microsoft/skala
X：https://x.com/MSFTResearch/status/1976318700216013069
論文：https://arxiv.org/abs/2506.14665
公式HP：https://labs.ai.azure.com/projects/skala/
検証環境
とあるスパコン
PBS 系ジョブスケジューラー
環境構築方法
conda 環境を設定する
conda create -n skala-cu118 python=3.11 -y
source /home/xxxx/miniconda3/etc/profile.d/conda.sh
conda activate skala-cu118
conda config --env --set channel_priority strict
conda config --env --set solver classic
pip install --upgrade pip
pip install torch==2.3.1+cu118 --index-url https://download.pytorch.org/whl/cu118
# まずビルド系を使わない前提でpipを最新化
pip install -U pip wheel setuptools
# NumPyはwheelのみ、かつ 2.0未満（1.26.4が多くの環境で安定）
pip install --only-binary=:all: "numpy<2"
# 残りも基本 wheel のみで（ビルドに落ちないように）
pip install --only-binary=:all: microsoft-skala pyscf e3nn qcelemental opt_einsum_fx \
huggingface_hub ase azure-core azure-identity
CONDA_SOLVER=classic conda install -n skala-cu118 -c conda-forge dftd3-python
# 構造最適化用に以下をインストール
pip install pyberny
pip install geometric
以下で動作チェック
python - <<'PY'
import torch, numpy as np, pyscf
print("torch:", torch.__version__, "CUDA:", getattr(torch.version, "cuda", None))
print("CUDA available:", torch.cuda.is_available())
print("numpy:", np.__version__)
print("pyscf:", pyscf.__version__)
PY
#動作チェックの結果は、以下の通り
torch: 2.3.1+cu118 CUDA: 11.8
CUDA available: True
numpy: 1.26.4
pyscf: 2.9.0
今回インストールしたものは、github で提供されている environment.yaml を参考に作成した。pytorch は、CPU 版ではなく GPU 版を入れる。
mamba を使わない方が簡単にインストールできる。
name: skala
channels:
- pytorch
- nvidia
- conda-forge
dependencies:
- python
- numpy
- ase
- azure-core
- azure-identity
- dftd3-python
- e3nn
- opt_einsum_fx
- pyscf <2.10.0
- qcelemental
- pytorch # ← GPU ビルド（CPU指定を削除）
- pytorch-cuda=12.1 # ← NVIDIA GPU 用
# Testing and development
- pre-commit
- pytest
- pytest-cov
- pip
- pip:
- huggingface_hub
実行方法
以下のように run_skala.py を作成する
gaussian と pySCF では、スピンの指定方法が異なるので注意が必要。
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Skala + PySCF:
- エネルギー計算
- 構造最適化（geomeTRIC）
使用例:
# エネルギーだけ
python run_skala.py --xyz mol.xyz --charge 1 --mult 1
# 最適化だけ
python run_skala.py --xyz mol.xyz --charge 1 --mult 1 --opt
"""
from __future__ import annotations
import argparse
from pathlib import Path
def read_xyz_to_atom_string(xyz_path: Path) -> str:
lines = xyz_path.read_text().strip().splitlines()
if len(lines) < 3:
raise ValueError("XYZが不適（原子数+コメント+座標が必要）")
atoms = []
for line in lines[2:]:
if not line.strip():
continue
sym, x, y, z = line.split()[:4]
atoms.append(f"{sym} {float(x):.10f} {float(y):.10f} {float(z):.10f}")
if not atoms:
raise ValueError("XYZから座標を読み取れませんでした。")
return "; ".join(atoms)
def write_xyz(path: Path, mol) -> None:
"""PySCF Mole から簡易 XYZ を書き出し（Å）"""
import numpy as np
coords = mol.atom_coords() # Bohr → to Angstrom
symbols = [mol.atom_symbol(i) for i in range(mol.natm)]
ang = coords * 0.529177210903
with path.open("w") as f:
f.write(f"{mol.natm}\noptimized by Skala/PySCF\n")
for s, (x, y, z) in zip(symbols, ang):
f.write(f"{s:2s} {x:16.10f} {y:16.10f} {z:16.10f}\n")
def main():
ap = argparse.ArgumentParser()
ap.add_argument("--xyz", required=True, type=Path, help="入力 .xyz（Å）")
ap.add_argument("--charge", type=int, default=0, help="電荷（例：+1 → 1, -1 → -1）")
ap.add_argument("--mult", type=int, default=1, help="多重度（Gaussianと同じ指定：1=一重項, 2=二重項, ...）")
ap.add_argument("--basis", type=str, default="6-31+G**", help="基底関数（既定: 6-31+G**）")
ap.add_argument("--opt", action="store_true", help="構造最適化を実行")
ap.add_argument("--freq", action="store_true", help="振動数解析を実行")
ap.add_argument("--grid-level", type=int, default=None, help="数値積分グリッドレベル（任意）")
args = ap.parse_args()
# Gaussianの多重度 → PySCFのspin(=2S) に変換
spin = args.mult - 1
atom_spec = read_xyz_to_atom_string(args.xyz)
# ===== PySCF + Skala =====
from pyscf import gto
from skala.pyscf import SkalaKS
mol = gto.M(
atom=atom_spec,
basis=args.basis,
charge=args.charge,
spin=spin, # 例：mult=1 → spin=0（閉殻RKS）
verbose=4,
)
mf = SkalaKS(mol, xc="skala")
if args.grid_level is not None:
mf.grids.level = args.grid_level
print(">>> SCF (Skala) start")
e0 = mf.kernel()
print(f"[Skala] SCF energy = {e0:.10f} Ha")
# ===== 構造最適化 =====
mol_eq = mol
if args.opt:
print(">>> Geometry optimization (geomeTRIC) start")
# geomeTRIC optimizer（PySCF同梱ラッパ）
from pyscf.geomopt.geometric_solver import optimize
# optimize は MF Scanner を使って内部で繰り返しSCFと勾配評価を行い、最終分子(Mole)を返す
mol_eq = optimize(mf)
# 最適化後、等価な MF を作り直してエネルギーを一度評価（キャッシュ目的）
mf = SkalaKS(mol_eq, xc="skala")
if args.grid_level is not None:
mf.grids.level = args.grid_level
eopt = mf.kernel()
print(f"[Skala] Optimized energy = {eopt:.10f} Ha")
# 最適化構造をXYZで保存
out_xyz = args.xyz.with_suffix(".opt.xyz")
write_xyz(out_xyz, mol_eq)
print(f"[Skala] Optimized geometry written to: {out_xyz}")
# ===== 振動数解析 =====
if args.freq:
print(">>> Harmonic frequency analysis start")
# 解析は対象Moleに対して行う（最適化したなら mol_eq）
mf_freq = SkalaKS(mol_eq, xc="skala")
if args.grid_level is not None:
mf_freq.grids.level = args.grid_level
mf_freq.kernel()
# まずは標準の Hessian を試す（実装が無い場合は例外 → 数値Hessianへ）
try:
hess = mf_freq.Hessian().kernel()
except Exception as err:
print("[WARN] Analytic Hessian not available; falling back to numerical Hessian.")
# 数値Hessian（PySCFの数値差分）にフォールバック
# 注：環境によりモジュール名が異なる場合があります
try:
from pyscf.hessian import numerical_hessian as numh
hess = numh.Hessian(mf_freq).kernel()
except Exception:
# 旧版互換のフォールバック
from pyscf.hessian import thermo as _thermo # ここは後の解析にも使う
# 数値Hessianを自前生成: geomeTRICやpyscf.toolsが必要な場合あり
raise RuntimeError(
"Numerical Hessian fallback failed. PySCFのhessianモジュールのバージョンを確認してください。"
) from err
# 振動計算解析（cm^-1 などを取得）
from pyscf.hessian import thermo
vibinfo = thermo.harmonic_analysis(mol_eq, hess, exclude_trans=True, exclude_rot=True)
# vibinfo = (freqs, modes, IRintens, Ramanintens, thermo_info, ...)
freqs_cm1 = vibinfo[0]
print("Frequencies (cm^-1):")
for i, w in enumerate(freqs_cm1, 1):
print(f" {i:3d}: {w:12.4f}")
# 想定外の虚数振動がある場合の注意書き
if any(w < 0.0 for w in freqs_cm1):
print("[NOTE] 虚振動が検出されました。")
if __name__ == "__main__":
main()
スパコンで実行する場合は、以下のようにジョブスクリプトを作成する。
#!/bin/bash
#PJM -L "elapse=72:00:00"
#PJM -S
#PJM -j
#PJM "--norestart"
#PJM --mail-list "sample@gmail.com"
#PJM -m "e"
set -euo pipefail
## SET YOUR DIRECTORY
USERHOME="/home/xxxx"
# ===== Optional: Modules =====
. /usr/share/Modules/init/sh
module load gcc/11.3.0
# pipのcu118を使うならCUDAモジュールは通常不要（競合の可能性あり）
# module load cuda/12.4.1
# ===== Conda =====
export PATH="${USERHOME}/miniconda3/condabin/:${PATH}"
# conda初期化
source "${USERHOME}/miniconda3/etc/profile.d/conda.sh"
conda activate skala-cu118
# ===== Hugging Face: 認証を無効化（401回避）=====
unset HF_TOKEN
unset HUGGINGFACE_HUB_TOKEN
unset HF_HUB_OFFLINE
rm -f ~/.cache/huggingface/token ~/.huggingface/token || true
export HF_HOME="${SCRATCH:-$HOME}/.cache/huggingface"
# ===== libgomp 競合回避 =====
export LD_LIBRARY_PATH="${CONDA_PREFIX}/lib:${LD_LIBRARY_PATH:-}"
export LD_PRELOAD="${CONDA_PREFIX}/lib/libgomp.so.1"
# ===== 実行 =====
python run_skala.py --xyz ./file_path/file_name.xyz --charge 1 --mult 1 --opt
opt の収束判定
log ファイルの最初の方に閾値が書いてある。
> ===== Optimization Info: ====
> Job type: Energy minimization
> Maximum number of optimization cycles: 300
> Initial / maximum trust radius (Angstrom): 0.100 / 0.300
> Convergence Criteria:
> Will converge when all 5 criteria are reached:
> |Delta-E| < 1.00e-06
> RMS-Grad < 3.00e-04
> Max-Grad < 4.50e-04
> RMS-Disp < 1.20e-03
> Max-Disp < 1.80e-03
Gaussian では、RMS Force, Max Force, RMS Displacement, Max Displacement で判定を行うが、PySCF では エネルギー変化も判定に用いている。PySCF の log ファイル内で Grad と書いてあるのが、Gaussian の Force に対応している。
各 step で SCF 計算が収束すると以下のように書かれる。
この部分が収束判定に使われる。
converged SCF energy = -976.466780030509
Step 8 : Displace = [92m7.082e-04[0m/[92m1.608e-03[0m (rms/max) Trust = 7.084e-04 ([1;91mx[0m) Grad = [0m3.725e-04[0m/[0m1.675e-03[0m (rms/max) E (change) = -976.4667800305 ([91m+2.522e-05[0m) Quality = [91m-3.029[0m
（注意）geomeTRIC の進捗サマリは、色付き表示のANSIエスケープコードが混ざっているため読みづらい。
- [92m = 緑にする
- [91m = 赤にする
- [1;91m = 太字＋赤
- [0m = 色をリセット（「0 m」という数値ではない）
この部分を Gaussian 形式で分かりやすく書くと以下になる。
Convergence criteria at Step 8
-------------------------------------------------------------
Item Value Threshold Converged?
Maximum Force 0.0016750 0.000450 No
RMS Force 0.0003725 0.000300 No
Maximum Displacement 0.0016080 0.001800 Yes
RMS Displacement 0.0007082 0.001200 Yes
Energy Change 0.0000252 0.000001 No
-------------------------------------------------------------
注意点
現時点では、典型元素用に開発されており、将来的には遷移金属にも適用される見込み。
解析的 Hessian が提供されていないため、現時点で用途は限られる。
構造最適化と一点計算に使うのが良いと思われる。
遷移状態構造の最適化、IRC 計算、振動計算等は厳しい。
上記のスクリプト内に振動数計算の部分を残してあるが、使うことはお勧めしない。