DFT計算をエンジンとしたクラスター展開・相分離シミュレーション: ATATとの連携#

Advance/PHASE Logo

新材料の設計において、第一原理計算（DFT）から得られる絶対零度でのミクロな電子状態やエネルギーを、有限温度におけるマクロな熱力学特性（相図、固溶限、相転移など）へとスケールアップする手法が求められています。本解析では、第一原理計算ソフトウェアAdvance/PHASEをエネルギー評価エンジンとして活用し、合金理論自動化ツールキット「ATAT (Alloy Theoretic Automated Toolkit)」と連携させることで、パラジウム－水素（Pd-H）系の複雑な相分離挙動シミュレーションを行いました。

Keywords: 第一原理計算 (DFT), ATAT, クラスター展開法 (Cluster Expansion), モンテカルロ法 (GCMC), Pd-H系, 水素吸蔵合金, 相分離, ヒステリシス

理論背景と計算手順#

クラスター展開法（Cluster Expansion: CE）#

第一原理計算を用いて任意の温度・濃度の相図を描くためには、膨大な数の原子配置パターンに対するエネルギー評価が必要となります。これを現実的な計算コストで実現するのがクラスター展開法 [1] です。Pd-H系のような侵入型合金にCE法を適用する場合、水素の侵入位置である「八面体サイト」を一つの格子とみなします。そして、各サイトに水素が存在すれば+1、空孔（水素がない状態）であれば-1という「スピン変数」 $\sigma_i$ を割り当て、イジングモデルのように系全体の配置エネルギー $E(\sigma)$ を定式化します。

$E(\sigma) = J_0 + \sum_{i} J_i \sigma_i + \sum_{i,j} J_{ij} \sigma_i \sigma_j + \sum_{i,j,k} J_{ijk} \sigma_i \sigma_j \sigma_k + \dots$

ここで、 $J$ は有効クラスター相互作用（ECI: Effective Cluster Interactions）と呼ばれるフィッティングパラメータです。

本解析では、ATAT [1] が自動生成した百数十パターン（数点除外）のスーパーセル構造に対して、Advance/PHASEによる第一原理構造最適化・全エネルギー計算を実行しました。そのDFT計算結果をATATにフィードバックし、高精度なECIを抽出しています。

ATATとの連携の工夫1：格子膨張の考慮#

Pd-H系では、水素の吸蔵に伴って格子体積が約10%も膨張する（ $\alpha$ 相から $\beta$ 相への転移）ことが知られています [2]。剛体格子（Rigid Lattice）を前提とする単純なCEモデルでは、この体積膨張に伴う弾性ひずみエネルギーを捉えきれません。本解析では、Advance/PHASEの構造最適化において、セル体積と内部原子座標のフルリラクゼーションを許容しました。最適化した各構造のDFTエネルギーをATATに学習させることで、局所的な格子体積の変動をCEモデルに反映させています。なお、マクロな相分離に伴う長距離の弾性ひずみエネルギー（Constituent Strain）の厳密な評価には追加の定式化が必要となりますが、本事例では計算パイプラインの構築を優先し、短距離相互作用に基づく基本モデルを採用しています。

ATATとの連携の工夫2：計算自動化とジョブ管理#

ATATを用いたクラスター展開では、数十から数百に及ぶスーパーセル構造に対して第一原理計算を連続的に実行する必要があります。本解析では、PythonおよびShellスクリプトを活用し、Advance/PHASEと計算サーバーをシームレスに繋ぐ以下の自動化・ジョブ管理機構を構築しました。

構造の自動生成と状態管理： ATATが提示する構造候補に対して、Advance/PHASE用の入力ファイルを自動生成し、計算準備完了（ready）、計算実行中（wait）、計算完了・データ抽出済み（done）といったディレクトリステータスをスクリプトで自動管理するシステムを構築しました。これにより、膨大な計算の抜け漏れや重複を防いでいます。また、計算が完了した構造から、波動関数や電荷密度などの巨大な不要バイナリファイルを自動削除してディスク容量を節約するスクリプトを導入しました。
ジョブ管理システム（PBS等）との連携・監視： 計算サーバーのジョブ管理システム（PBSなど）と連携し、投入中のジョブステータス（submitted など）を常時監視する仕組みを導入しました。サーバーの最大投入可能ジョブ数を超えないように、キューの空き状況に応じて自動的に次のDFT計算ジョブを投入する制御を行い、限られた計算リソースを最大限に活用してハイスループットな解析を実現しています。

計算結果と考察#

1. ECIの収束性とCEモデルの予測精度#

Advance/PHASEで得られたエネルギー群からATATによって抽出されたECIのクラスター半径依存性を図1に、各構造の生成熱に対する残差（フィッティングエラー）を図2に示します。

ECI vs Cluster radius
図1. 有効クラスター相互作用（ECI）とクラスター半径の関係。長距離領域で若干のばらつきが見られるものの、絶対値としては数 meV/atom の精度範囲内に十分に収まっており、物理モデルとして妥当な相互作用の減衰が確認できます。

Residual error of the fit
図2. Pd-Hの各構造パターンに対するフィッティングの残差エラー。誤差は最大でも数 meV/atom のオーダーに抑えられており、Advance/PHASEのDFTエネルギーが精度を損なうことなくマクロな統計力学モデルへマッピングされていることが分かります。

2. 基底状態の探索（Convex Hull）#

図3は、横軸に水素濃度 $x_H$ 、縦軸に生成熱をプロットしたConvex Hull（凸包）です。赤点がAdvance/PHASEによるDFT計算値、黒線が絶対零度（0 K）における熱力学的な安定状態の結線を示しています。両端だけでなく、 $x_H \approx 0.57$ や $0.71$ などの中間濃度にも深い谷（安定な規則相：Ordered Phase）が存在する結果が得られました。なお、どのような規則相が最安定な基底状態として現れるかは、採用する計算セルサイズの上限や、長距離弾性ひずみ考慮の有無といったモデルの制約に依存することが知られています。本解析のモデル設定においては、この中間相の強固な安定性が、後の有限温度での相転移挙動（プラトーの形成など）に大きな影響を与えています。

Calculated energies and Convex Hull
図3. Pd-H系 (PdH_x, x：0～1) 基底状態のConvex Hull（凸包）。中間の水素濃度領域にPd-H特有の安定な規則相が存在することを示しています。

3. グランドカノニカル・モンテカルロ法による相分離とヒステリシスの再現#

構築した高精度CEモデルを用い、ATATの emc2 モジュールにて温度300 Kでのグランドカノニカル・モンテカルロ（GCMC）シミュレーションを実施しました。化学ポテンシャル $\mu$ （外部圧力に対応）を連続的に増減させ、平衡となる水素濃度 $x_H$ を評価した結果（等温線：Isotherm）を図4に示します。

Hysteresis Loop at 300 K
図4. 300 Kにおける化学ポテンシャル $\mu$ vs 水素濃度 $x_H$ のヒステリシスループ。青線が水素吸蔵過程（Sweep UP）、赤線が放出過程（Sweep DOWN）を示します。

図4から以下の特徴が読み取れます。 * ヒステリシスの発現（二相共存）： 行き（青線）と帰り（赤線）で経路が異なり、中央に大きな隙間（ヒステリシスループ）が形成されています。これは、水素がまばらに溶け込んだ $\alpha$ 相から、水素化パラジウムである $\beta$ 相へと不連続に濃度がジャンプする「第一種相転移」がシミュレートされている証拠です。 * 中間相（プラトー）の存在： 0から1へ一気に濃度が飛ぶのではなく、 $x_H \approx 0.5$ 付近で階段状の平坦部（プラトー）を形成しています。絶対零度（図3）では $x_H = 0.5$ の状態は凸包のわずかに上に位置していますが、有限温度（300 K）においてはエントロピーの効果等によってこの濃度での相が安定化し、シミュレーション上で明確なプラトーとして現れています。現実のPd-H系では、後述する格子振動や四面体サイトへの占有効果などによってこのプラトーは消失し、より滑らかな転移になることが示唆されていますが、シミュレーション手法としてはモデルのエネルギー地形をGCMCで再現できていることが確認できます。

4. モンテカルロ・スナップショットによる相分離状態の可視化#

図5は、300 KにおけるGCMCシミュレーションのスナップショットです。計算には1458原子（Pd原子：729個、H/空孔サイト：729個）の大きなスーパーセルを用いています。図中の白球は水素原子（H）、茶球は空孔（X: Vacancy）を表しています（ホストであるPd格子は省略して表示しています）。

Monte Carlo Snapshot showing Phase Separation
図5. 300 KにおけるGCMCシミュレーションのスナップショット。スーパーセル内において、水素が密集した領域（H: 白球）と空孔が密集した領域（X: 茶球）が分離し、二相共存状態を形成している様子。

周期境界条件が適用された有限サイズのシミュレーションセル内において、系が界面エネルギーを最小化しようと再配列した結果、このように水素領域（ $\beta$ 相）と空孔領域（ $\alpha$ 相）が層状（スラブ状）に分かれた相分離構造として観察されます。第一原理DFT計算から抽出されたミクロな相互作用パラメータが、モンテカルロ法の枠組みにおいて相分離のダイナミクスを適切に駆動していることが視覚的にも確認できました。

5. 本モデルの物理的限界とより高度な解析への展望#

本事例では第一原理計算とCE・GCMC法をシームレスに連携させるパイプラインを構築・実証しましたが、実材料のPd-H系の相図を定量的に完全再現するには、以下の物理的要因をモデルに組み込む拡張が必要になることが最新の研究で指摘されています。

長距離弾性ひずみ（Constituent Strain）の明示的考慮： 短距離のクラスター展開では切り捨てられてしまう長距離の弾性ひずみは、コヒーレントな相分離におけるヒステリシスや臨界温度に決定的な影響を与えます [2]。
四面体サイト（T-site）占有モデルの導入： 本解析では水素が八面体サイト（O-site）のみを占有する格子モデルを用いましたが、Pd-H系の高濃度・高温域の正確な記述には、T-site占有を許容したモデル化が有効とされています [3]。
格子振動（振動エントロピー）の付加： 静的モデルの限界として原子の熱的な揺らぎが欠落していますが、有限温度での振動エントロピーの差は合金の相安定性に無視できない補正を与えます [4, 5]。

Advance/PHASEとATATの連携システムは、これらの高度な拡張（Constituent Strain項の追加やフォノン計算による振動エントロピーの加算など）を組み込むための柔軟な基盤を備えています。

まとめ#

本事例では、第一原理計算ソフトウェアAdvance/PHASEのDFT計算とATATを連携させることで、パラジウム中の水素吸蔵プロセスをマルチスケールに解析しました。Advance/PHASEによるエネルギー・体積評価から高精度なCEモデルを構築し、ATATによるモンテカルロシミュレーションやPythonスクリプトなどによるジョブ自動化を通じて、「絶対零度のミクロな原子配列」から「有限温度でのマクロな相分離・ヒステリシス挙動」へのハイスループットなスケールアップを自動実行する計算パイプラインを構築しました。このような計算アプローチは、水素吸蔵合金のみならず、二次電池の電極材料（リチウム挿入・脱離反応）や、各種合金の析出・時効熱処理プロセスなど、様々な材料開発において強力な計算手法の基盤となります。

本解析の詳細や、研究への適用可能性に関するご相談はこちら

お問い合わせ

参考文献#

A. van de Walle, M. Asta, and G. Ceder, "The alloy theoretic automated toolkit: A user guide," CALPHAD 26, 539 (2002).
J. M. Rahm, J. Löfgren, and P. Erhart, "Quantitative predictions of thermodynamic hysteresis: Temperature-dependent character of the phase transition in Pd-H", Acta Materialia 227, 117697 (2022).
P. L. Hagelstein, "Statistical Mechanics Models for PdDx and PdHx Phase Diagrams with both O-site and T-site Occupation", J. Condensed Matter Nucl. Sci. 29, 401 (2018).
A. van de Walle and G. Ceder, "The effect of lattice vibrations on substitutional alloy thermodynamics", Rev. Mod. Phys. 74, 11 (2002).
B. Fultz, "Vibrational thermodynamics of materials", Prog. Mater. Sci. 55, 247 (2010).