mt-property-analysis

name

mt-property-analysis

description

物性計算とMD後処理を扱うスキルです。 ElasticTensorFeature, PostElasticPropertiesFeature, 弾性定数, elastic, Young率, バルク弾性率, ForceConstantFeature, PostPhononBandFeature, PostPhononDOSFeature, PostPhononThermochemistryFeature, フォノン, phonon, 力定数, QHA, 準調和近似, VibrationFeature, 振動解析, PostMDDiffusionFeature, 拡散係数, MSD, 平均二乗変位, PostEMDViscosityFeature, PostNEMDViscosityFeature, 粘度, Green-Kubo, PostNEMDThermalConductivityFeature, 熱伝導率, PostMDSpecificHeatFeature, 比熱, PostMDThermalExpansionFeature, 熱膨張, GasStandardFormationEnthalpyFeature, PostVibrationGasThermoFeature, ガス熱化学, RDF, 動径分布関数, 振動状態密度 に関するコード生成時に使用してください。

物性計算

概要

Matlantis を用いた物性計算と MD 後処理のガイドです。PFP（Preferred Potential）の高速な原子間ポテンシャル計算を活用して、弾性定数、フォノン、準調和近似（QHA）、拡散係数、粘度、熱伝導率、動径分布関数（RDF）、振動状態密度など、幅広い物性を評価できます。

matlantis-features ライブラリは、これらの物性計算を高レベル API として提供しています。各 Feature は callable なオブジェクトであり、result = feature(atoms) のパターンで統一的に利用できます。多くの物性計算は 2 段階パイプライン（計算 Feature + 後処理 PostFeature）で構成されています。

本ガイドでは、理論的な背景と実装を分離し、MD trajectory の生成コードと後処理コードを明確に区別して説明します。

ワークフロー

物性計算の標準的なワークフローは以下の通りです。

1. 構造の準備: 対象系の構造を作成し、構造最適化を完了する
2. 計算手法の選択: 目的の物性に応じた Feature クラスを選択する
3. 一次計算の実行: ForceConstantFeature、ElasticTensorFeature、MDFeature 等を実行する
4. 後処理の実行: Post-Feature で結果を解析・変換する
5. 結果の解析: 数値の妥当性を検証し、可視化する

Feature 選択の判断基準

計算対象	一次 Feature	後処理 Feature
弾性定数	ElasticTensorFeature	PostElasticPropertiesFeature
フォノンバンド	ForceConstantFeature	PostPhononBandFeature
フォノン DOS	ForceConstantFeature	PostPhononDOSFeature
熱力学量	ForceConstantFeature	PostPhononThermochemistryFeature
QHA	ForceConstantFeature	PostPhononQHAGibbsFeature
拡散係数	MDFeature	PostMDDiffusionFeature
比熱	MDFeature	PostMDSpecificHeatFeature
熱膨張係数	MDFeature	PostMDThermalExpansionFeature
粘度（EMD）	MDFeature	PostEMDViscosityFeature
粘度（NEMD）	MDFeature	PostNEMDViscosityFeature
熱伝導率	MDFeature	PostNEMDThermalConductivityFeature
振動解析	VibrationFeature	---
気体熱力学	VibrationFeature	PostVibrationGasThermoFeature
標準生成エンタルピー	---	GasStandardFormationEnthalpyFeature

実装パターン

パターン A: 弾性定数

最適化済み結晶に小ひずみを付加し、応力応答から弾性テンソルを算出します。Young 率、剪断弾性率、体積弾性率、Poisson 比を得ることができます。

理論

弾性定数は、結晶に微小ひずみを付加した際の応力応答から求めます。

1. 原子位置とセルの両方を最適化した結晶構造を用意する
2. 対角ひずみ（伸張・圧縮）と非対角ひずみ（剪断）を複数の大きさで付加する
3. 各ひずみ下で原子位置のみを再緩和し（セル形状は固定）、応力を計算する
4. ひずみ-応力の線形フィッティングから 6×6 の弾性定数テンソルを求める

実装

import numpy as np
from matlantis_features.features.elasticity import (
    ElasticTensorFeature,
    PostElasticPropertiesFeature,
)
from matlantis_features.features.common import FireASEOptFeature
from matlantis_features.utils.calculators import pfp_estimator_fn

estimator_fn = pfp_estimator_fn(model_version="v9.0.0", calc_mode="R2SCAN")

# Step 0: 原子位置とセルの両方を最適化（弾性テンソル計算の前提）
#   セルも緩和するため filter=True の optimizer を使う
cell_opt = FireASEOptFeature(filter=True, estimator_fn=estimator_fn)
opt_result = cell_opt(atoms)

# Step 1: 弾性テンソルの計算
#   diagonal_strains: 伸張・圧縮ひずみ / off_diagonal_strains: 剪断ひずみ（いずれも必須引数）
diagonal_strains = [-0.015, -0.01, -0.005, 0.005, 0.01, 0.015]
off_diagonal_strains = [-0.03, -0.02, -0.01, 0.01, 0.02, 0.03]
elastic = ElasticTensorFeature(
    diagonal_strains=diagonal_strains,
    off_diagonal_strains=off_diagonal_strains,
    # ひずみ印加後の「原子位置のみ」再緩和用。filter を付けない（=セル固定）新しいインスタンスを渡す。
    # Step 0 の filter=True の optimizer は使い回さない（後述の注意参照）。
    optimizer=FireASEOptFeature(estimator_fn=estimator_fn),
    check_symmetry=True,        # 格子対称性に基づく既約テンソルの検証
    estimator_fn=estimator_fn,
)
elastic_result = elastic(opt_result.atoms)

# 6x6 弾性テンソル（GPa 単位）。ElasticTensor オブジェクトの elastic_tensor 属性で取得
print(np.round(elastic_result.elastic_tensor))

# Step 2: 物性値の算出
props = PostElasticPropertiesFeature()
results = props(elastic_result)

# 結果の取得（属性名は単数形: young_modulus / poisson_ratio）
print(f"Bulk Modulus    {results.bulk_modulus:3.2f} GPa")
print(f"Young's Modulus {results.young_modulus:3.2f} GPa")
print(f"Shear Modulus   {results.shear_modulus:3.2f} GPa")
print(f"Poisson ratio   {results.poisson_ratio:3.2f}")

注意:

• ElasticTensorFeature は diagonal_strains・off_diagonal_strains・optimizer が必須引数です。estimator_fn だけでは初期化できません。
• optimizer は使い回さない: optimizer にはひずみ印加後の原子位置のみを緩和するもの（filter=False、デフォルト）を渡します。Step 0 のセル最適化用（filter=True）optimizer をそのまま渡すと、各ひずみ状態でセルまで緩んで応力が消え、非現実的な弾性定数になります（ElasticTensorFeature は filter 有効時に警告を出します）。事前最適化用とは別のインスタンスを用意してください。
• 物性値の属性名は単数形です: young_modulus・poisson_ratio（youngs_modulus / poissons_ratio は存在しません）。ほかに bulk_modulus・shear_modulus・anisotropy・vicker_hardness・compliance_tensor などが取得できます。方向依存の値は results.get_prop_results(dir1=[1,1,1], dir2=[1,-1,0]) で取得します。

パターン B: フォノン計算

フォノンバンド構造、状態密度（DOS）、および熱力学量を計算します。力定数の計算を起点として、段階的に各種物性を算出します。

理論

フォノン計算は有限変位法（Finite Displacement Method）に基づきます。

1. primitive cell を使用する（計算効率とバンド構造の正確性のため）
2. 十分なサイズのスーパーセルを構築する
3. 各原子を微小変位させてフォルダー定数を計算する
4. フォノンバンド、DOS、熱力学量を段階的に算出する

実装

from matlantis_features.features.phonon import (
    ForceConstantFeature,
    PostPhononBandFeature,
    PostPhononDOSFeature,
    PostPhononThermochemistryFeature,
)

estimator_fn = pfp_estimator_fn(model_version="v9.0.0", calc_mode="R2SCAN")

# Step 1: 力定数の計算
# supercell: スーパーセルのサイズ（大きいほど正確だが遅い。4-6で開始推奨）
# delta: 変位量 (Angstrom)
fc = ForceConstantFeature(
    supercell=6,
    delta=0.01,
    estimator_fn=estimator_fn,
)
fc_result = fc(primitive_atoms)

# Step 2a: バンド構造
band = PostPhononBandFeature()
band_result = band(fc_result)

# Step 2b: 状態密度
dos = PostPhononDOSFeature(mesh=[20, 20, 20])
dos_result = dos(fc_result)

# Step 2c: 熱力学量（自由エネルギー、エントロピー、比熱など）
thermo = PostPhononThermochemistryFeature()
thermo_result = thermo(fc_result)

パターン C: 準調和近似（QHA）

QHA は、異なる体積でのフォノン計算結果を統合し、有限温度での熱力学量を精密に評価します。

from matlantis_features.features.phonon import (
    ForceConstantFeature,
    PostPhononQHAGibbsFeature,
)

# 複数体積でのForceConstant計算結果を用いてQHA Gibbsエネルギーを算出
fc = ForceConstantFeature(supercell=6, delta=0.01, estimator_fn=estimator_fn)

# 体積を変えた複数構造でfc_resultを取得
# ...

qha = PostPhononQHAGibbsFeature()
qha_result = qha(fc_results)

パターン D: 拡散係数

MD trajectory から平均二乗変位（MSD）を計算し、拡散係数を算出します。

理論

1. 十分な長さの MD trajectory を生成する
2. MSD を計算する（周期境界を考慮した unwrapping が必要）
3. MSD の線形領域を特定してフィッティングする
4. Einstein の関係式から拡散係数を求める

実装

from matlantis_features.features.md import (
    MDFeature,
    PostMDDiffusionFeature,
)

estimator_fn = pfp_estimator_fn(model_version="v9.0.0", calc_mode="R2SCAN")

# Step 1: MD trajectoryの生成
md = MDFeature(
    estimator_fn=estimator_fn,
    # MD パラメータ（温度、ステップ数、アンサンブル等）
)
md_result = md(atoms)

# Step 2: 拡散係数の計算
diffusion = PostMDDiffusionFeature()
diffusion_result = diffusion(md_result)

MSD 解析の注意点

• Unwrapping: 周期境界条件下では原子がセル境界を越えると座標が折り返されるため、MSD 計算前に座標を unwrap する必要があります
• 線形領域の特定: MSD が時間に対して線形になる領域（バリスティック領域の後、統計誤差が大きくなる前）を選択する必要があります
• 統計精度: 複数の時間起点から MSD を平均化し、統計誤差を低減させます

パターン E: 粘度

平衡 MD（EMD）の Green-Kubo 法、または非平衡 MD（NEMD）による粘度計算です。

理論（Green-Kubo 法）

• 平衡 MD で応力テンソルの自己相関関数を計算する
• Green-Kubo 積分により粘度を求める
• 統計誤差の評価と平均化が重要

実装

from matlantis_features.features.md import (
    MDFeature,
    PostEMDViscosityFeature,
    PostNEMDViscosityFeature,
)

# EMD（Green-Kubo法）
md_result = md(atoms)
emd_viscosity = PostEMDViscosityFeature()
emd_result = emd_viscosity(md_result)

# NEMD
nemd_viscosity = PostNEMDViscosityFeature()
nemd_result = nemd_viscosity(md_result)

パターン F: 熱伝導率

from matlantis_features.features.md import (
    MDFeature,
    PostNEMDThermalConductivityFeature,
)

md_result = md(atoms)
thermal = PostNEMDThermalConductivityFeature()
thermal_result = thermal(md_result)

統計誤差や平均化の必要性に注意してください。単一の MD run からの結果は統計的に不十分な場合があります。

パターン G: 比熱・熱膨張係数

from matlantis_features.features.md import (
    PostMDSpecificHeatFeature,
    PostMDThermalExpansionFeature,
)

# 比熱
specific_heat = PostMDSpecificHeatFeature()
sh_result = specific_heat(md_result)

# 熱膨張係数
thermal_expansion = PostMDThermalExpansionFeature()
te_result = thermal_expansion(md_result)

パターン H: 動径分布関数（RDF）

MD trajectory から RDF を計算し、構造の特徴を解析します。

# MD trajectoryからRDFを計算
# ase.geometry.analysis または matlantis-features の後処理機能を使用
from ase.geometry.analysis import Analysis

analysis = Analysis(trajectory)
rdf = analysis.get_rdf(rmax=10.0, nbins=200)

パターン I: 振動解析と気体熱力学

分子の振動モードと、気体状態の熱力学量を計算します。

from matlantis_features.features.common.vibration import VibrationFeature
from matlantis_features.features.common.gas_thermo import PostVibrationGasThermoFeature
from matlantis_features.features.common.gas_formation_enthalpy import (
    GasStandardFormationEnthalpyFeature,
)

estimator_fn = pfp_estimator_fn(model_version="v9.0.0", calc_mode="R2SCAN")

# 振動解析
vib = VibrationFeature(delta=0.01, estimator_fn=estimator_fn)
vib_result = vib(molecule_atoms)

# 気体熱力学量（理想気体近似）
gas_thermo = PostVibrationGasThermoFeature()
gas_result = gas_thermo(vib_result)

# 標準生成エンタルピー
formation = GasStandardFormationEnthalpyFeature(estimator_fn=estimator_fn)
formation_result = formation(molecule_atoms)

パターン J: 振動状態密度（VDOS from MD）

MD trajectory から振動状態密度を計算します。速度自己相関関数のフーリエ変換に基づきます。

# MD trajectoryから速度自己相関関数を計算し、
# フーリエ変換で振動状態密度を得る
# 実装は解析ライブラリに依存

ベストプラクティス

構造最適化の徹底

物性計算の入力構造は必ず事前に最適化してください。結晶系では filter=True または UnitCellASEFilter を使ってセル形状も最適化する必要があります。

Primitive Cell の使用（フォノン計算）

フォノン計算では primitive cell を使用してください。Conventional cell を使用すると、不要なバンドフォールディングが発生し、結果の解釈が困難になります。

スーパーセルサイズの選択（フォノン計算）

ForceConstantFeature(supercell=N) の N は大きいほど正確ですが計算コストが増大します。N=4-6 から開始し、結果が収束しない場合に増やしてください。

理論と実装の分離

物性計算のコードを書く際は、理論的な説明（なぜこの手法を使うか）と実装（具体的なコード）を明確に分離してください。

MD 生成と後処理の分離

MD trajectory の生成コードと後処理（MSD 計算、Green-Kubo 積分等）のコードは分離してください。これにより、同じ trajectory に対して異なる解析を適用したり、trajectory の生成条件を変更したりすることが容易になります。

統計誤差への配慮

粘度、熱伝導率などの輸送特性の計算では、統計誤差が大きくなりやすいです。以下の点に注意してください。

• 十分な長さの MD trajectory を生成する
• 複数の独立した MD run から結果を平均化する
• 相関時間を考慮したブロック平均を使用する

estimator_fn の使用

matlantis-features の Feature クラスには、Estimator インスタンスではなく estimator_fn（ファクトリ関数）を渡してください。Feature 内部で必要に応じて新しい Estimator が生成され、インスタンス共有による問題を防ぎます。

よくあるエラーと対処

エラー・症状	原因	対処法
フォノンバンドに虚の周波数（負の値）が出る	構造が十分に最適化されていない	より厳密な fmax（例: 0.001 eV/A）で再最適化する
フォノン計算が非常に遅い	supercell が大きすぎる	supercell を小さくして結果の収束を確認する
MSD が線形にならない	MD が平衡化していない、または trajectory が短すぎる	平衡化フェーズを長くする、trajectory を延長する
粘度の値が不安定	統計サンプル不足	MD run を延長する、複数 run を平均化する
MultiCalculatorUseDetected	Estimator を共有した	estimator_fn パターンを使用する
弾性定数が対称でない	構造が十分に最適化されていない	セル最適化（filter=True）を含む最適化を実行する
Conventional cell で不要なバンドが出る	primitive cell を使っていない	spglib 等で primitive cell に変換する

関連ガイド

• 構造最適化: 物性計算の前段階として必須
• 反応経路探索 (reaction/SKILL.md): 振動解析を遷移状態検証に使用
• 可視化 (visualization/SKILL.md): 計算結果の構造確認
• Group Drive ストレージ (storage/SKILL.md): 大規模計算結果の保存