材料計算学
本書は、材料物理・材料工学で用いられる計算手法を、数式に基づいて一貫して学ぶための教科書として設計するものである。連続体から原子スケールまでの記述を往復しながら、計算結果の信頼性を担保するための理論的根拠と手法選択の判断軸を身につける構成である。
第1章 数値モデル化と支配方程式
本章では、材料現象を計算可能な形にするために、どの量を変数として選び、どの方程式を立てるのかを数式の形で整理する。連続体としての記述と原子・粒子としての記述を往復しながら、モデル化の前提、支配方程式の導出、条件設定の意味を一つずつ確認する。
1.1 連続体モデルと粒子モデル
材料計算では、状態を場として表す連続体モデルと、原子や粒子の集まりとして表す粒子モデルの二つが基本である。連続体モデルでは、位置
両者の違いは、何を「最小単位」とみなすかにある。連続体モデルは、材料内部の微細構造を直接には追わず、ある体積要素の平均量で状態を表す考え方である。これは、観測したい現象の空間スケールが原子間距離より十分大きく、平均化した量の変化で現象を説明できる場合に有効である。例えば、厚さ数mmの試験片の温度分布を議論するなら、原子の個数を数えるよりも、温度場
この平均化の考え方を数式で表すと、連続体の場は粒子の情報から粗視化して得られる量として定義できる。例えば、密度場
のように書ける。ここで
粒子モデルは、局所構造や熱揺らぎ、欠陥の原子レベル機構を直接追える点が強い。例えば拡散は、原子のランダムなジャンプの積み重ねであり、粒子モデルではその起点をそのまま記述できる。一方で粒子数
連続体モデルは、長時間・大域の挙動を扱いやすい点が強い。例えば拡散方程式や弾性方程式は、空間に分布する量の滑らかな変化として現象を表すため、メッシュ幅を現象の代表長さに合わせて選べる。ただし、材料定数(弾性定数、拡散係数、界面エネルギーなど)が入力として必要になり、それらがどの条件で有効な値かを理解しないと結果の解釈を誤りやすい。連続体モデルは「何を捨てたか」を明示しながら使う姿勢が重要である。
材料計算の多くは、この二つを使い分けるだけではなく、接続する。粒子モデルで得た局所物性を連続体モデルに渡し、連続体モデルで得た場を粒子モデルの外場として与える。以降の章で扱うDFT、MD、MC、フェーズフィールド、有限要素法は、この「どの自由度で書くか」という選択の上に配置される。
1.2 保存則と生成消滅項
材料現象の支配方程式は、保存則から作るのが最も堅牢である。保存則は、ある領域に含まれる量の増減が、境界を通る流れと内部での生成消滅で決まるという主張である。これをまず積分形で書くと、領域
となる。ここで
次に、発散定理
を用いると、積分形は
となる。これは任意の領域
が得られる。この導出は、材料計算で繰り返し現れる基本操作であり、方程式の形がなぜこの形になるのかを理解する最短経路である。
この一般形の強みは、物理を
を採用すると、
が得られる。
同じ枠組みで熱伝導も書ける。保存量を単位体積あたりの内部エネルギー
が得られる。
運動量の保存則では、
のように書ける。弾性体では
生成消滅項
以降の章で扱う多くの方程式は、見かけは異なっても、この保存則の枠組みか、次節の変分原理の枠組みに分類できる。どちらの出発点から来ているかを意識すると、式変形や境界条件の意味が理解しやすくなる。
1.3 変分原理と自由エネルギー
材料の平衡状態は、自由エネルギーが最小となる状態として決まる場合が多い。ここで言う自由エネルギーは、温度一定ならヘルムホルツ自由エネルギー、圧力一定ならギブズ自由エネルギーなどであり、条件に応じて適切な熱力学ポテンシャルを選ぶ必要がある。連続体の場
のように書ける。
平衡条件は、
となる。第二項に部分積分(多次元では発散定理に相当)を適用すると
が得られる。境界での変分の扱いが適切に定まっているとすると、領域内部での平衡条件は
となる。これがオイラー・ラグランジュ方程式であり、変分原理から支配方程式が導かれる基本形である。
材料で特に重要なのは、勾配エネルギー項を含む自由エネルギーである。例えば
を考える。ここで
なので、平衡条件は
となる。すなわち、勾配エネルギー項はラプラシアン
変分原理は平衡条件だけでなく、時間発展のモデル化にも使える。平衡へ向かって自由エネルギーが減少する緩和を表したいなら、非保存量
と置くのが自然である。ここで
であり、結果として
が得られる。これは非保存秩序変数の緩和方程式の基本形である。
一方、保存される量(濃度
と置き、保存則
とする。ここで
この節で導いた「自由エネルギーから駆動力が出る」という構造は、DFTでは電子密度に対する変分、フェーズフィールドでは秩序変数に対する変分、有限要素法では弱形式(エネルギー原理)として再登場する。式の見た目が変わっても、根にある考え方は共通である。
1.4 初期条件と境界条件
支配方程式は、方程式だけでは解が決まらず、初期条件と境界条件が必要である。初期条件は
代表的な境界条件は次の三つである。ディリクレ条件は場の値を固定する条件であり、例えば温度固定なら
である。ノイマン条件は法線方向の勾配や流束を固定する条件であり、熱流束
である。ロビン条件は値と流束の線形結合を与え、対流境界
のように、外部との熱交換を表すのに使われる。
周期境界条件は、無限に繰り返す材料やバルクの一部を表したいときに用いる。例えば一次元なら
のように、値と勾配の一致を課す。これは境界の影響を消すための道具であり、欠陥や界面がある場合は、周期像との相互作用が生じないようセルサイズを十分大きく取る必要がある。
吸収境界条件は、波動や弾性波、電磁波などで計算領域の端から反射が戻ってくるのを抑える目的で用いる。例えば波動方程式で、領域端で外向きに進む波を想定し
のような形を置くと、端での反射が減る場合がある。吸収の良さは周波数や入射角に依存するため、問題に応じて設計が必要である。
初期条件は、時間依存問題では必須である。拡散方程式なら
境界条件は数値計算の実装とも深く結びつく。有限要素法の弱形式では、ディリクレ条件は本質境界条件として試験関数空間を制約する形で入る。一方、ノイマン条件は自然境界条件として境界積分項に現れ、弱形式に組み込まれる。これは、方程式の導出(部分積分で境界項が出る)と境界条件の種類が対応していることを意味する。
材料計算における条件設定は、物理の再現だけでなく、問題が数学的に解ける形になっているかにも関わる。例えば拡散方程式で全境界をノイマン条件(流束指定)にすると、総量保存のため定数分だけ解が不定になる場合がある。そのときは平均値を別途固定するなど、解の自由度と条件数の関係を意識して設計する必要がある。
1.5 無次元化とスケールの切り出し
無次元化は、式の中に隠れている重要な比を取り出し、何が支配的かを見える形にする操作である。さらに、数値計算では適切な格子幅
とおき、方程式を
例として拡散方程式
を無次元化する。
であり、
となる。
で見積もれる。
次に移流拡散方程式
を考える。速度の代表値を
となる。ここで
となる。
このような無次元数は、数値条件にも直結する。例えば陽解法で拡散方程式を解く場合、空間を格子幅
のような上限が現れる(
が現れる。したがって、現象の支配項が何であるかを無次元化で把握すると、どの制約が支配的になるかも見通せる。
さらに、無次元化はモデルの単純化にも使える。ある無次元数が十分小さいなら、その項を無視した近似が成立しうる。例えば
材料計算では、同じ名前の方程式でも材料定数が桁違いになることがある。無次元化は、材料ごとに式の意味を揃える作法であり、異なる系の結果を比較するための共通のものさしになる。
1.6 連続体と離散体の接続
連続体モデルと粒子モデルは、どちらか一方を選べば終わる関係ではない。材料物性は階層的であり、原子スケールの相互作用がメソスケールの組織を作り、それがマクロな応答(強度、導電率、磁性など)を決める。したがって、計算でも階層をまたぐ情報の受け渡しが必要になる。
第一の方向は、離散体から連続体への粗視化である。代表例は、原子計算から材料定数を抽出して連続体モデルに渡す操作である。例えば弾性定数
と展開できることを用い、いくつかの独立なひずみ状態を計算して二次係数をフィットする。ここで重要なのは、微小ひずみ領域で二次近似が妥当であることを確認し、数値ノイズや応力の収束誤差が係数に混ざらないようにする点である。
拡散係数
で定義できる。ここで
第二の方向は、連続体から離散体への情報供給である。例えば、連続体側で得られた応力場や温度場を、原子計算の境界条件や外場として与え、局所構造変化を調べることがある。温度勾配下の拡散、応力誘起の相変態、電場下のイオン移動などでは、マクロな場がミクロな遷移率を変えるため、この方向の接続が重要になる。
連続体と離散体の接続では、同じ名前の量でも定義が異なる場合がある点に注意が必要である。例えば応力は、連続体では面に作用する力の密度として定義されるが、原子系ではビリアル応力として
のように表される。ここで
さらに、界面や粒界のように連続体の仮定が破れやすい領域では、接続の設計が計算結果の解釈を左右する。例えばフェーズフィールドで界面エネルギー
連続体と離散体を接続する作業は、材料計算の中心である。支配方程式の形を理解し、保存則と自由エネルギーという二つの出発点を押さえたうえで、どの物理量をどの階層で定義し、どう平均し、どう渡すのかを一貫させることが、計算結果を材料科学の議論へつなげる基本動作である。
第2章 離散化と偏微分方程式の数値解法
本章では、偏微分方程式をコンピュータで解くために、連続な式を有限個の未知数へ落とし込む手順を、式変形から順に整理する。空間離散と時間離散の設計が、解の安定性・精度・計算量を同時に決めることを、最小の例から一般形へ広げて理解する。
2.1 空間離散の設計
偏微分方程式は、未知量が連続な位置
まず、一次元で直感を作る。例えば
である。これらを引き算すると
となり、中心差分
が得られる。誤差が
が得られ、拡散方程式の空間離散の基礎になる。
差分法は、このような導関数近似を多次元に拡張して方程式を点ごとに満たす。例えば二次元のラプラシアンは
となる。格子が直交で、境界形状が単純な場合に実装が簡潔で高速に動きやすい。一方で、材料の実形状が複雑である場合や、界面が曲線・曲面である場合には、格子に形状を合わせることが難しくなる。
有限体積法は、保存則の積分形から出発し、セル平均を未知数として扱う。保存則
をセル
である。発散定理により
となるので、セル内総量の変化はセル境界を通る流束の和で決まる形になる。これを離散化すると、一次元では
のようになり、数値的に「出入りの差で増減する」構造がそのまま残る。材料の拡散・熱伝導・流体輸送などで、保存が重要な量を扱うときに強い。ただし、高次精度を狙う場合や界面近傍での流束評価(再構成)を丁寧に設計しないと、数値拡散や振動が出やすい。
有限要素法は、微分方程式を弱形式に変換してから近似する。例としてポアソン方程式
を考える。試験関数
とする。左辺に部分積分(発散定理)を適用すると
が得られる。境界項はノイマン境界条件として自然に入り、ディリクレ境界条件は
と展開する。
となる。ここで
スペクトル法は、未知関数をフーリエ級数や直交多項式で展開する。周期境界の一次元で
と書けるとき、微分は係数空間で
となり、高い精度で微分が計算できる。解が十分滑らかである場合、誤差が非常に速く減る性質を持つが、界面の不連続や局所的な急峻勾配があるとギブス現象が出やすい。材料では界面や欠陥が本質の場合が多いため、スペクトル法は適用条件を見極めて使う必要がある。
材料現象では、界面近傍の急峻な勾配、結晶方位に依存する異方性、複雑な形状、複数場(温度・応力・濃度など)の連成が同時に現れる。そのため、どの離散化が何を得意とし、どの条件で不利になるのかを、式の構造と幾何の扱いから説明できることが重要である。
| 手法 | 強み | 注意点 |
|---|---|---|
| 差分法 | 構造が単純な領域で高速化しやすく、導関数近似の意味が追いやすい | 形状が複雑だと格子生成や境界近傍の精度が難しくなりやすい |
| 有限体積法 | 保存則が離散形でも保たれ、輸送問題で物理の整合性を保ちやすい | 流束再構成や界面条件の設計が結果を左右しやすい |
| 有限要素法 | 複雑形状・異方性・連成を弱形式で統一でき、境界条件も整理しやすい | 要素品質と数値積分、拘束条件の実装が精度と安定性に直結する |
| スペクトル法 | 解が滑らかな場合に非常に高精度になりやすい | 不連続・界面を含む問題では振動が出やすく工夫が要る |
2.2 時間積分と安定性
空間離散により、偏微分方程式は有限次元の常微分方程式系に落ちる。未知ベクトルを
の形になる。時間積分は、この系を離散時間
最も基本の陽解法は前進オイラー法で、
である。実装が簡潔で並列化もしやすい一方、安定性が時間刻みを強く制約する。安定性を理解するために、線形の試験方程式
を考える。前進オイラー法では
であり、増幅率
となり、
具体例として一次元拡散方程式
を中心差分で離散化すると
となる。前進オイラー法で時間離散すると
である。フォン・ノイマン解析により、安定には
が得られる。格子を細かくして界面を解像しようとすると、時間刻みが二乗で小さくなるため、計算回数が急増する。
波動や移流では別の制約が出る。例えば移流方程式
は情報が速度
が現れる(
これらの制約を緩めるのが陰解法である。後退オイラー法は
であり、右辺が未知の
となり、
精度の観点では、時間積分の次数も意識する必要がある。前進オイラー法と後退オイラー法は一次精度であり、時間刻みを半分にすると誤差はおよそ半分になる。二次精度の代表としてクランク・ニコルソン法があり、線形問題では
である。拡散型では安定である一方、強い非線形や不連続があると振動的な誤差が残る場合があるため、安定性と減衰性の両方を見て選ぶ必要がある。
材料計算では、拡散のように硬い項と、比較的緩やかな項が混ざることが多い。このとき、硬い項だけ陰に、その他を陽にするIMEX(Implicit-Explicit)法や、演算子分割が有効になる。例えば
で
のようにして、安定性と計算量を両立させる設計ができる。時間積分は、単に式を進める操作ではなく、安定性と計算コストの配分を決める設計問題である。
2.3 疎行列と線形ソルバ
空間離散と時間離散を行うと、多くの場合で線形方程式
を解くことになる。有限要素法の剛性行列、陰解法の時間発展行列、ポアソン方程式の離散化行列などが代表である。ここで重要なのは、
直接法(ガウス消去、LU分解、Cholesky分解)は頑健であり、中規模問題では扱いやすい。しかし三次元で自由度
反復法の理解のために、線形方程式を残差で見る。近似解
を残差と呼び、
という部分空間上で最良近似を探す方法である。行列の全要素を扱うのではなく、疎行列の行列ベクトル積だけを中心に計算が進むため、大規模問題で有利になる。
非対称の場合(移流を含む、非自己共役な離散化、あるいは連成によりブロック構造になる場合)はGMRESなどが用いられる。GMRESは残差ノルムを最小化するが、反復回数とともに記憶量が増えるため、再直交化やリスタートなどの工夫を用いる。材料計算では、物性の強い異方性や不連続があると、行列が悪条件になりやすく、ソルバの選択が計算時間を支配する。
ここで決定的に重要なのが前処理である。前処理は、元の問題
を解く。理想は
線形ソルバの停止判定は、単に反復回数を固定すればよいわけではない。残差が十分小さいかを基準にするが、残差の大きさは右辺
を用いることが多い。さらに、物理量として意味のある誤差(例えばエネルギー誤差や保存量誤差)を併せて監視しないと、数値的には収束していても物理的には不適切な解になる場合がある。線形ソルバは、PDEの解法の中心部であり、離散化と同じくらい「設計」の対象である。
2.4 非線形ソルバと変分問題
材料の相変態、塑性、接触、化学反応、非線形拡散などでは、離散化後に非線形方程式
を解く必要がある。ここで
の形になる。非線形性は、構成式が非線形であること、自由エネルギーが多井戸形であること、あるいは制約条件が入ることから生じる。
基本となるのはニュートン法である。現在の近似
とする。ここで
を繰り返す。解の近くでは二次収束するため非常に速いが、初期値が悪いと発散することがある。また、各反復で線形方程式を解く必要があり、その線形ソルバの性能が全体を支配する。
ニュートン法を安定に使うために、線形探索や減速を入れる。例えば更新を
とし、残差が確実に減るように
ニュートン法が重い場合、固定点反復(ピカール反復)を使うこともある。例えば
の形に書けるとき
で更新する。収束は遅いが、ヤコビアンを作らずに済み、実装が簡単になることがある。塑性の内部変数更新や、弱い非線形の拡散係数更新などで用いられる。ただし、収束条件は写像の縮小性に依存するため、緩和係数を入れて
のようにする場合も多い。
変分問題として書ける場合は、非線形ソルバの見通しがよくなる。例えばエネルギー汎関数
は、最適性条件
例えば、時間発展がエネルギーを減らすべき現象では、離散化後も
が成立するように時間積分を設計すると、数値的不安定や非物理的振動を抑えやすい。この考え方は、フェーズフィールドの安定化スキームや、非線形拡散の陰解法設計で特に重要である。非線形ソルバは、単に収束させるだけでなく、物理として守るべき性質を壊さない形で解くことが目標になる。
2.5 マルチグリッドと前処理
大規模PDEの計算時間の多くは線形ソルバが占めるため、収束を速くする工夫が中心課題になる。前処理はその中核であり、マルチグリッドは前処理の最も強力な考え方の一つである。発想は、誤差には波長の短い成分と長い成分があり、反復法は短波長誤差を減らすのが得意だが、長波長誤差を減らすのが苦手である、という観察に基づく。
例えば一次元ポアソン方程式を考えると、Jacobi法やGauss-Seidel法は局所的な平均化なので、格子1つ2つで振動する誤差はすぐ減る。一方、領域全体でゆっくり変化する誤差は、局所平均化ではなかなか消えない。ここで、粗い格子へ移ると、その「ゆっくり変化する誤差」が粗格子では短波長に見えるため、少ない反復で減らせる。これがマルチグリッドの基本原理である。
マルチグリッドの構成は、滑らか化(smoothing)、制限(restriction)、粗格子補正(coarse-grid correction)、補間(prolongation)からなる。線形方程式
を解き(あるいは近似的に解き)、誤差補正
幾何学的マルチグリッドは、格子やメッシュが階層構造(2倍粗くするなど)を持つ場合に適用しやすい。材料計算では複雑形状の有限要素メッシュが多く、規則的な粗格子が作りにくい場合がある。そのときは代数的マルチグリッド(AMG)が有効であり、行列
マルチグリッドは、それ自体をソルバとして使うだけでなく、Krylov法の前処理として使うことが多い。例えばCGやGMRESの前処理にAMGを入れると、反復回数が大きく減る場合がある。ここでの狙いは、
材料計算では、メッシュの局所細分化、異方性、物性の桁違いの変化がしばしば現れ、これらがソルバを難しくする。マルチグリッドと前処理は、単なる計算高速化ではなく、計算自体を成立させるための中心技術である。離散化と同じく、どの誤差をどの階層で減らすのか、という視点で理解すると、手法の選択理由が明確になる。
2.6 並列計算と再現性
材料計算は自由度が大きくなりやすく、単一計算機では時間もメモリも足りないことが多い。そこで分散メモリ並列(MPI)による領域分割が標準となる。領域
疎行列計算では、計算の中心は行列ベクトル積
並列計算で重要なのは、再現性が揺れる可能性である。浮動小数点の加算は結合法則を満たさず、
が一般に成り立つ。並列では、総和(内積、ノルム、残差評価など)を複数プロセスで分担し、最後に集約する。この集約の順序が実行ごとに変わると丸め誤差の入り方が変わり、最終結果がわずかに変化する場合がある。線形・非線形反復は、停止判定や分岐によりこの微小差を増幅することがあるため、数値の揺れを完全にゼロにするのは難しい。
この問題への対応は、目的に応じて複数ある。第一に、停止判定を物理量に基づいて頑健にし、丸め誤差程度の揺れで分岐しないようにする。例えば相対残差の閾値を過度に厳しくしすぎず、必要精度に合わせる。第二に、可能なら決定的な還元(同じ順序で総和を取る)や高精度和(補償和)を用いて揺れを減らす。第三に、検証では一度の実行結果だけでなく、メッシュ・時間刻み・ソルバ設定を変えたときに物理的結論が不変であることを示す。
疎行列ライブラリ(PETScなど)は、並列疎行列、Krylovソルバ、前処理、非線形ソルバを統一的に提供し、現代の材料計算の基盤となっている。重要なのは、ライブラリを使うこと自体ではなく、数値解法の構造(残差、安定性、前処理、スケーリング)を理解したうえで、設定が何を意味しているかを説明できることである。並列化は計算を速くする手段であると同時に、数値誤差の扱い方を再点検する契機でもある。
小まとめ
離散化は、連続な材料方程式を有限次元の代数方程式へ変換する操作であり、その選択が安定性・精度・計算量を決める。今後の章で扱うDFT、モンテカルロ、分子動力学、フェーズフィールド、有限要素法のいずれも、内部では本章の考え方に沿って数値的に解かれているため、方程式から解法までを一つの流れとして捉えることが材料計算学の基礎になる。
第3章 密度汎関数理論の理論基盤
本章では、量子多体系としての電子問題を、電子密度という場の変数へ落とし込む論理を、定理と変分計算として整理する。計算で現れるKohn–Sham方程式が、どのエネルギー汎関数の最小化問題として導かれるのかを、式のつながりが追える形で示す。
3.1 電子密度と基礎定理
電子の基底状態は、多体波動関数
を主変数に据えることで回避する理論である。密度
Hohenberg–Kohn(HK)の第一定理は、外場ポテンシャル
を定める。
HKの第二定理は、密度に関する変分原理を与える。すなわち、ある外場
を最小にする密度が基底状態密度
は普遍汎関数であり、外場に依存しない形で定義される。
変分原理の計算としては、粒子数保存
を拘束条件として、ラグランジュ未定乗数(化学ポテンシャル)
を満たす密度が平衡密度になる。ここでの汎関数微分は、後にKohn–Sham方程式の形に直結する。HK定理は、密度という場を用いた変分問題として電子基底状態が定まることを保証し、計算に向けた数学的な入口を与える。
ただし、HK定理は存在定理であり、
3.2 Kohn–Sham方程式
Kohn–Sham(KS)形式の要点は、相互作用する電子系の基底状態密度を、ある非相互作用電子系の密度として再現できるようにすることである。非相互作用系では、
と表される(スピンを明示する場合は後述する)。
KSのエネルギー汎関数は、普遍汎関数を
と分解して定義する。
である。残りの差分が交換相関エネルギー
としてまとめられる。ここで
次に、軌道
と書く。ここで
は原子単位系での表現である。軌道は正規直交条件
を満たす必要があるため、ラグランジュ未定乗数行列
をとる。
が得られる。ここで
である。
と書ける。これがKohn–Sham方程式である。見た目は一電子シュレディンガー方程式に似ているが、有効ポテンシャルが密度
自己無撞着(SCF)の構造は次のように理解できる。密度
を解いて軌道
を作り、
KS方程式は、量子力学の基礎式と変分原理から導かれるため、式の意味を追えることが強みである。以降の節では、この枠組みの中で未知である交換相関汎関数と、計算上の工夫(擬ポテンシャル、相対論、精度管理)を体系的に整理する。
3.3 交換相関汎関数の近似
交換相関汎関数
LDA(局所密度近似)は、各点で電子密度が一様電子気体であるとみなし、
と置く。
で与えられる(原子単位系)。相関は量子モンテカルロなどの参照データをパラメータ化して与える。LDAは、密度がゆっくり変化する金属的結合に強い一方、分子や界面など密度変化が大きい領域では誤差が増えやすい。
GGA(一般化勾配近似)は密度勾配を導入し、
の形を取る。交換部分では、無次元勾配
を用いて、LDA交換に増強因子
meta-GGAは、密度と勾配に加えて、運動エネルギー密度
などの軌道由来の量を使うことで、結合様式(共有結合、金属結合、弱結合など)の識別力を上げる。SCANは多くの厳密条件を満たすように設計され、meta-GGAの中で広く参照されてきた。数値安定性を改善する工夫を加えたr2SCANのような改良も提案されており、実装と収束性の観点から選択されることがある。meta-GGAは表面・界面、遷移金属化合物などで改善する例がある一方で、ポテンシャルの非線形性が増すためSCFが不安定になりやすい場合もあり、混合法や初期密度の工夫が効きやすい。
ハイブリッド汎関数は、厳密交換(Fock交換)を混合する。厳密交換エネルギーは占有軌道から
で与えられ、局所・半局所近似より非局所的な交換を取り込む。一般形として
のように書ける。バンドギャップや局在電子の記述が改善する例がある一方、計算量が増え、金属では収束が難しくなる場合がある。スクリーン型(短距離だけ厳密交換を入れる)などの工夫は、この計算負担と収束性を緩和するために導入される。
弱い分散相互作用(ファンデルワールス力)を扱うには、半局所汎関数だけでは不十分な場合がある。経験的な分散補正(
| 近似の階層 | 参照量 | 取り込む情報 | 計算上の性質 |
|---|---|---|---|
| LDA | 局所一様電子気体 | 収束しやすいことが多い | |
| GGA | 密度の空間変化 | 多くの材料で基準として使われる | |
| meta-GGA | 局所環境の識別力向上 | 非線形性が増し収束が難しくなる場合がある | |
| ハイブリッド | 厳密交換混合 | 非局所交換 | 計算量増大、金属で難しい場合がある |
交換相関汎関数は、物性予測の中心誤差源になりやすい。したがって、同じ構造でも汎関数を変えて結果の変化を見て、議論の頑健性を確かめる姿勢が、材料物性へ接続するうえで不可欠である(第4章で相互検証の観点を整理する)。
3.4 擬ポテンシャルとPAW
平面波基底で周期系を扱う場合、原子核近傍の内殻電子は波動関数が急峻に振動するため、これをそのまま展開すると非常に高い平面波カットオフが必要になる。材料物性の多くは価電子が担うため、内殻電子を凍結し、価電子だけを有効ポテンシャルで扱う考え方が擬ポテンシャルである。目的は、価電子の散乱特性を保ったまま、波動関数を滑らかにして基底数を減らすことである。
擬ポテンシャルの基本要求は、あるカットオフ半径
のように、角運動量
ultrasoft擬ポテンシャルはノルム保存条件を緩めることで、より低いカットオフで計算できるようにする。その代わり、重なり演算子が単位ではなくなり、一般化固有値問題
として解く必要がある。ここで
PAW(projector augmented-wave)法は、全電子波動関数と擬波動関数を線形変換で結ぶ枠組みである。概念的には
という変換演算子
擬ポテンシャルやPAWでは、どの電子を価電子に含めるか(半芯状態を含めるか)が物性に影響することがある。例えば磁性や結合に寄与する準位が核近傍にある場合、半芯を凍結すると誤差が増えることがある。逆に価電子を増やすと計算量が増えるため、対象物性(格子定数、磁気異方性、電場勾配など)に応じて選ぶ必要がある。
また、擬ポテンシャルは生成条件(参照配置、交換相関汎関数、相対論取り扱い)に依存するため、DFT計算では汎関数と擬ポテンシャルの整合が重要である。汎関数だけを変えて擬ポテンシャルを据え置くと、厳密には理論整合が崩れる場合がある。材料計算ではこの点をどの程度厳密に扱うかが現実的な判断になるが、少なくとも差の出やすい物性では慎重に相互検証を行うべきである。
3.5 スピン、相対論、磁性
材料機能においてスピン自由度は中心的である。スピン分極DFTでは、電子密度をスピンごとに分けて
とし、交換相関エネルギーを
となる。ここで
非共線磁性では、スピン量子化軸が空間的に変化するため、スピン密度はベクトル場
スピン軌道相互作用(SOC)は、磁気異方性、Dzyaloshinskii–Moriya相互作用、異常ホール効果などに直接つながる。相対論的な出発点はDirac方程式であるが、固体計算では多くの場合、Pauli近似としてSOC項が有効ハミルトニアンに現れる。中心ポテンシャルを仮定すると、原子内SOCは
の形で書け、
相対論効果はSOCだけではない。重元素ではスカラー相対論(質量増大や収縮効果)が結合長やバンド位置に影響する。擬ポテンシャルやPAWの生成時に相対論を含めたものが提供されるのは、この効果を価電子有効ポテンシャルに取り込むためである。SOCを後から摂動的に入れる方法(第二変分法)を用いる実装も多く、計算負担と精度のバランスを取る設計になっている。
磁性材料を扱うときは、スピン分極の初期条件(初期磁気モーメント、反強磁性配置など)が収束先を変えることがある。これはエネルギー地形が多谷性を持つためであり、どの磁気秩序が安定かを比較するには、複数の初期条件とセル設定で相互に確認する必要がある。SOCを含めるとさらに自由度が増え、SCFが難しくなる場合があるため、段階的に計算設定を強めていくのが合理的である。
3.6 精度管理と相互検証
DFTは原理計算であるが、数値計算である以上、近似と離散化の影響を管理しなければならない。誤差の源は大きく分けて、交換相関汎関数の近似、擬ポテンシャルやPAWの近似、基底と離散化(平面波カットオフ、実空間グリッド、数値積分)、ブリルアンゾーン積分(
平面波DFTでは、波動関数の展開をエネルギーカットオフ
スーパーセル計算では、周期境界により欠陥像どうしが相互作用する。電荷欠陥では、長距離クーロン相互作用により収束が遅く、電荷補正やポテンシャル整列が必要になることがある。界面や表面でも、真空層厚さ、スラブ厚さ、双極子補正などの設計が結果を左右する。これらは数値条件の選択に見えるが、根は周期境界という数学的仮定にあるため、相互作用の長さスケールを見積もることが本質である。
相互検証は、誤差の切り分けに有効である。例えば、汎関数を変える、擬ポテンシャル(あるいはPAWデータセット)を変える、コードを変えるといった比較により、物理由来の差と数値由来の差を分離できる。全エネルギーの絶対値はコード間で直接比較できない場合が多いが、同一手順で評価したエネルギー差(形成エネルギー、相安定差など)は比較しやすい。比較は単に一致を確認する作業ではなく、どの仮定が結論に効いているかを理解する作業である。
本章の結論は、DFTが定理と変分原理で立ち上がる厳密な枠組みを持つ一方で、実用では交換相関近似と数値設定が物性値を左右するという点である。次章では、DFT計算を材料物性へ結びつけるために、スケーリング、ブリルアンゾーン積分、構造最適化、格子力学、欠陥・界面、有限温度といった具体的課題を、式と計算の観点から整理する。
第4章 DFT計算から材料物性へ
本章では、Kohn–Sham方程式を解いた結果を、材料科学で議論する物性(構造安定性、弾性、フォノン、電子状態、欠陥形成、界面エネルギー、有限温度自由エネルギーなど)へ変換する手順を、数学的定義と数値上の要点として整理する。DFTの出力をそのまま読むのではなく、何を計算し、どの差を物性として解釈するのかを、式の形で明確にする。
4.1 基底関数と計算スケーリング
Kohn–Sham方程式は連続な偏微分方程式であるため、何らかの基底に展開して有限次元の固有値問題として解く。平面波基底は周期境界と相性が良く、基底をエネルギーカットオフで系統的に増やせるため、収束の扱いが比較的分かりやすい。周期セル体積を
に比例して増える。固有値問題の計算は概ね行列ベクトル積と直交化が中心であり、バンド数が増えるほど計算負担が増大する。一般に、バンド法のコストは系サイズに対して
局在基底(原子軌道基底、数値原子軌道、ガウス基底など)は、原子の周りに局在した関数で波動関数を展開する。基底数は原子数に比例しやすく、疎行列構造が得られるため、大規模系で有利になりやすい。数値局在基底を用いると、行列要素の評価が実空間積分として整理でき、線形スケーリング法(密度行列の局在性を利用して
が
ただし、局在基底は基底の選び方が収束性に影響しやすい。平面波では
計算スケーリングの理解は、材料物性の対象に応じて手法を選ぶ判断につながる。例えば、フォノンや欠陥でスーパーセルが必要な場合、平面波で
4.2 ブリルアンゾーン積分と金属の収束
周期系の電子状態はブロッホ波で記述され、波数
であり、エネルギーも同様に
ここで
金属の難しさは、フェルミ面近傍で占有数が急に変わり、積分の被積分関数が滑らかでなくなる点にある。絶縁体では占有がギャップで分離されるため、
で滑らかにし、数値積分を安定化する方法が用いられる。これはスメアリングと呼ばれ、実際には
になる。ここでエントロピー項
で与えられる。
テトラへドロン法は、BZを小さな四面体に分割し、エネルギー分散を線形補間して積分する方法である。絶縁体や半導体では高精度に積分でき、金属でも改良版(占有の扱いを工夫した方法)が用いられることがある。ただし、磁性金属や準二次元系など、フェルミ面が複雑で状態密度が鋭く変化する場合には、
低次元系(薄膜、表面、界面)では、周期方向が限られるため、
4.3 構造最適化と格子力学
DFT計算で最初に行うことが多いのは、原子位置
で定義され、平衡条件は
構造最適化は、数値的には多変数最小化問題であり、共役勾配法、準ニュートン法(BFGSなど)、減衰付き分子動力学などが用いられる。これらは、力と(場合により)近似ヘッセ行列から更新量を決め、エネルギーを下げる方向へ進む。重要なのは、SCF収束が不十分だと力がノイズを含み、最適化が不安定になる点である。したがって、力収束閾値とSCF収束閾値は連動して設定する必要がある。
格子力学は、平衡構造の周りで原子を微小変位させたときの二次の応答からフォノンを求める枠組みである。原子
と展開できる。
で定義される。質量
が得られ、固有値問題
を解くことでフォノン分散
力定数は有限変位法(原子を少し動かして力差分から求める)でも、密度汎関数摂動論(DFPT)でも求められる。有限変位法では、変位量
のように数値微分する。
フォノンから熱力学へ接続するには、調和近似の振動自由エネルギー
を用いる。これにより、温度依存の自由エネルギー
4.4 励起状態と強相関の拡張
標準DFTは基底状態理論であり、Kohn–Sham固有値
GW法は、多体摂動論に基づき、自己エネルギー
で近似する(
を解くことで得られる。これによりバンドギャップが改善する例が多いが、計算量が大きく、収束(空虚準位数、誘電関数のカットオフ、
BSE(Bethe–Salpeter方程式)は、電子と正孔の相互作用を取り込み、励起子効果を含む光学応答を記述する。線形応答として、電子正孔対の有効ハミルトニアンを作り、その固有値として励起子準位を得る枠組みである。TDDFT(時間依存DFT)は、密度の時間発展と線形応答を扱い、光吸収スペクトルなどへ接続する。TDDFTの精度は交換相関カーネルの近似に依存し、長距離相互作用や励起子を精密に扱うには工夫が必要になる。
強相関系では、局在した
のように書け、部分占有を罰することで局在化を促す。ここで二重計数補正(DFT側で既に含まれている相互作用を引く項)の選び方が結果に影響し、
DMFT(動的平均場理論)は、局所相互作用を動的自己エネルギーとして扱い、温度依存やスペクトル関数の形を改善する。DFT+DMFTは、DFTのバンド構造とDMFTの局所相関を結合し、金属絶縁体転移や有限温度磁性などをより物理的に扱える。ただし、計算負担が大きく、モデル部分空間の定義や二重計数の扱いが難しい。
これらの拡張は、何を追加し、何を近似しているかを明確にしたうえで選ぶ必要がある。基底状態の構造安定性だけなら標準DFTで十分な場合も多いが、光学応答、準粒子ギャップ、局在相関、有限温度スペクトルなどを議論するなら、拡張法が必要になる。
4.5 欠陥、界面、有限温度
材料物性の多くは、完全結晶ではなく欠陥や界面で決まる。DFTで欠陥を扱う基本はスーパーセル法であり、周期セル内に欠陥を導入して周期境界で繰り返す。欠陥形成エネルギーは、化学ポテンシャルと電荷状態を含めて
のように定義される。
界面・表面では、エネルギーを面積で割った量が重要になる。表面エネルギーはスラブの全エネルギー
のように定義できる(両面を持つ対称スラブの場合)。ここで
有限温度では、基底状態エネルギー
を比較する必要がある。
4.6 計算コードと国内外の基盤
DFT計算は多様なコードで実行でき、実装の選択は基底、擬ポテンシャル、機能(DFPT、SOC、非共線、ハイブリッド、GWなど)、並列性能に依存する。平面波擬ポテンシャルDFTは周期系に強く、材料科学で標準的に使われる機能が整備されていることが多い。局在基底DFTは大規模系や線形スケーリングに接続しやすく、欠陥・界面・ナノ構造の大規模計算で有利になる場合がある。どのコードを用いる場合でも、基底と擬ポテンシャルの整合、収束試験、相互検証が物性の信頼性を支える。
国内の計算基盤を利用する場合、利用可能なソフトウェア、コンパイラやMPI、推奨設定が環境ごとに異なる。大規模計算では、疎行列ソルバやFFTライブラリなどの性能が計算時間を左右するため、コードの実装特徴と計算機の構成(CPU、メモリ帯域、ネットワーク)との相性も重要になる。計算は再現可能であるべきなので、入力ファイル、擬ポテンシャルの版、汎関数、
また、物性へ接続するには、単一の計算結果だけでなく、構造の同定、相の比較、欠陥濃度の仮定、温度・化学ポテンシャルの設定など、熱力学的な枠組みの中で整理する必要がある。コードはその一部を担う道具であり、材料科学としての結論は、何を比較し、どの自由度を含めた自由エネルギーを議論したか、という定義の上に立つ。
小まとめ
第3章では、HK定理と変分原理からKohn–Sham方程式が導かれ、未知の多体相互作用が交換相関汎関数へ集約される構造を示した。第4章では、DFTの出力を物性へ変換するために、基底とスケーリング、BZ積分、構造最適化とフォノン、励起状態拡張、欠陥・界面・有限温度といった論点を式の形で整理した。今後は、これらのDFTで得たエネルギーや力、応答関数を、モンテカルロ法、分子動力学、フェーズフィールド、有限要素法などの多階層計算へ接続し、有限温度・有限濃度・非平衡の材料現象を記述する方向へ進むことになる。
第5章 モンテカルロ法の基礎
本章では、熱平衡の分布を数値的にサンプリングするという立場から、モンテカルロ法を統計力学と確率過程の言葉で組み立てる。材料科学で必要になるのは、単に乱数で平均を取ることではなく、所望の分布に従う標本列を設計し、その誤差を制御して物理量へ結び付けることである。
5.1 統計集団とボルツマン重み
温度
である。ここで分配関数
(離散状態の場合)または
(連続状態の場合)で定義される。モンテカルロ法の核心は、
(または積分)を
として推定する点にある。
材料科学では、温度以外の制御変数が重要になる。例えば、化学ポテンシャル
となる。圧力
さらに、材料の相転移や秩序化は、秩序変数の揺らぎとして観測される。秩序変数
を推定できれば、平均値だけでなく相共存や自由エネルギー障壁も議論できる。したがって、モンテカルロ法は「平均を出す道具」であると同時に「分布を推定し、相空間の形を読む道具」である。
5.2 Metropolis法と詳細釣り合い
正準分布
で与えられる。目標は、十分大きい
である。これを直接満たす
を課すと、不変分布の条件が保証される(両辺を
実装上、遷移は「提案」と「受理」に分けて書くと理解しやすい。まず、提案分布
であり、棄却された場合は
となる。詳細釣り合いは
であるから、受理確率の比は
となる。これを満たす具体形として、Metropolis–Hastingsの受理確率
が得られる。特に提案が対称
がMetropolis法の基本形である。エネルギーが下がる遷移(
もう一つの要件はエルゴード性である。これは、許される状態空間の任意の状態へ、有限回の遷移で到達できることを要求する。例えばスピン模型で「必ず1個だけ反転する」更新は多くの場合エルゴードであるが、制約条件(粒子数固定、電荷中性など)を入れると更新が到達不能な領域を作ることがある。したがって、詳細釣り合いだけでなく、更新規則が状態空間を十分に連結しているかも設計の一部である。
5.3 誤差評価と自己相関
モンテカルロで得られる標本列
の分散は、独立標本の
を定義し、正規化した自己相関
を考える。ここで
と書ける。ここで
を統合自己相関時間と呼ぶ。独立標本なら
に減ると解釈できる。材料の相転移近傍や巨大欠陥集団の更新では
実務上は、
さらに、初期条件の影響(バーンイン)も誤差源である。初期状態が平衡分布から遠いと、最初の区間は平衡サンプルとして扱えない。バーンイン長は、エネルギーや秩序変数の時系列が定常になったか、複数の初期条件(秩序相・無秩序相など)から始めても同じ分布へ収束するか、といった観点で評価する必要がある。誤差評価は、数値計算の末尾処理ではなく、サンプリング設計と同じ重みを持つ。
5.4 クラスタ更新
臨界点近傍では相関長が発散的に大きくなり、局所更新(1スピン反転など)では大域的な緩和が極端に遅くなる。この現象が臨界減速であり、自己相関時間が系サイズ
のように増大する(
クラスタ更新は、相関の強い自由度をまとめて反転することで、長波長の緩和を速める方法である。最も基本的な例は強磁性Ising模型
である。温度が低いと同符号のスピンがドメインを形成し、境界だけが揺らぐ。局所更新は境界を少しずつ動かすしかないため遅いが、クラスタ更新では同符号の連結成分を一括で反転できる。
Wolff法では、ある種の確率で結合(ボンド)を張り、クラスタを成長させて反転する。隣接スピンが同符号である辺に対して、ボンドを張る確率を
とし、1点から始めて同符号隣接へ確率
Swendsen–Wang法は、系全体に対して同様にボンドを張り、得られた全クラスタを独立に確率
ただし、クラスタ法は万能ではない。強い異方性、外場、フラストレーション、連続自由度(ベクトルスピン)などでは拡張が必要であり、場合によってはクラスタが巨大化しすぎて効率が落ちる。材料科学では「相関の構造を見て更新を選ぶ」ことが基本であり、クラスタ更新はその代表例である。
5.5 拡張アンサンブル
自由エネルギー地形が多峰性を持つと、単一温度の正準サンプリングは特定の谷に閉じ込められやすい。相共存、一次相転移、核生成、ガラス的遅い緩和などでは、この問題が顕在化する。拡張アンサンブルは、分布そのものを拡張・混合し、状態空間をより広く探索することで、障壁越えを促進する枠組みである。
レプリカ交換法(平行テンパリング)は、異なる温度
である。ここで温度交換(
となる。高温レプリカは障壁を越えやすく、低温レプリカは目的温度での分布を精密に記述する。交換により低温レプリカも間接的に障壁越えを行えるため、一次転移近傍などで有効である。温度列の設計は受理率と拡散速度を左右し、エネルギー揺らぎの大きさに応じて温度間隔を狭める必要がある。
多重カノニカル法は、エネルギー分布を平坦化するような重み
でサンプリングする。理想的には状態密度
傘サンプリング(umbrella sampling)は、秩序変数
でサンプリングする。目的は、普段はほとんど出現しない
拡張アンサンブルは「分布を変えて探索し、再重み付けで物理量へ戻す」という思想で統一できる。したがって、どの分布をどの領域で平坦化したいのか(温度、エネルギー、秩序変数、反応座標)を明確にすることが設計の中心となる。
5.6 自由エネルギー推定
自由エネルギーは、相安定性や相境界、欠陥濃度、界面・核生成障壁を議論する上で中心的である。正準集団の自由エネルギーは
で定義されるが、
自由エネルギー摂動(free energy perturbation)は、基準系
を与える。ここで
熱力学積分(thermodynamic integration)は、連続的な結合パラメータ
のように補間して
で評価する。補間経路を適切に選べば、各
再重み付け(reweighting)は、ある温度や重みで得たサンプルから近傍条件の物理量を推定する。正準分布の温度変更では
が成り立つ。これは便利であるが、
自由エネルギー推定法は、どれか一つを暗記するのではなく、(i) どの分布でサンプルしているか、(ii) 目的の自由エネルギー差がどの対数比として定義されるか、(iii) 分布重なりが十分か、という観点で選ぶべきである。以下に、材料科学で出番が多い推定法を比較する。
| 推定法 | 基本式 | 強み | 注意点 |
|---|---|---|---|
| 自由エネルギー摂動 | 形式が簡潔で実装しやすい | 分布重なりが悪いと破綻しやすい | |
| 熱力学積分 | 経路を作れば安定しやすい | ||
| 再重み付け | 比の形で | 近傍条件の補間に有効である | 離れた条件には使いにくい |
| 傘サンプル+統合 | 障壁近傍を直接サンプルできる | 反応座標選定と窓設計が効く |
本章の結論は、モンテカルロ法は「所望分布のサンプリング設計」と「相関を含む誤差評価」の両輪で成立しているという点である。次章では、この枠組みが材料科学の具体問題(合金秩序、相図、拡散・析出、粒成長など)へどう接続されるかを、モデル化と手法選択として整理する。
第6章 材料科学におけるモンテカルロ
本章では、材料の相転移・秩序化・拡散・成長といった現象を、確率変数と有効ハミルトニアンの選び方として定式化し、モンテカルロ法で何が計算できるかを具体化する。DFTやMDが局所エネルギーやミクロな運動を扱うのに対し、MCは平衡分布と統計揺らぎを正面から扱うため、対象スケールと目的物性に応じた連携が鍵になる。
6.1 格子模型と秩序変数
相転移を理解する第一歩は、自由度を最小限に抽出した格子模型である。代表はIsing模型であり、格子点
で定義される。
秩序変数は、相を区別する量として定義される。Ising模型では磁化
が秩序変数である。材料の秩序化(例えばB2秩序、L1_0秩序など)でも、サブラティス占有の差として同様の秩序パラメータが定義できる。秩序変数の平均値だけでなく揺らぎも重要であり、磁化の分散から磁化率が
で与えられる。これは、線形応答(
相関長を議論するには、二点相関関数
を用いる。無秩序相では
Potts模型は、
と書ける。これは多結晶の粒成長モデル(粒の方位を
6.2 クラスター展開と置換合金
置換合金では、原子が格子点を占有し、どの元素がどこに入るか(配置)が自由度になる。温度により短距離秩序が生じ、相分離や秩序化が起きるが、これを第一原理のエネルギーと統計力学でつなぐ代表がクラスター展開である。二元合金
で表す。すると、配置
と書ける。結晶対称性を用いると、同等なクラスターは同じ係数を共有するため、実際にはクラスターの種類
と表すことが多い。
ここで
クラスター展開が得られると、正準MCで温度依存の秩序変数や比熱を計算できる。例えば内部エネルギーの熱平均
で評価でき、相転移点ではピークとして現れる。組成が変動する大正準MCを使えば、化学ポテンシャル
この枠組みの落とし穴は、エネルギーだけでなく、格子歪みや磁気自由度が重要な材料では、有効ハミルトニアンが単純な置換自由度だけでは閉じない点である。必要に応じて、格子緩和を取り込んだ拡張クラスター展開や、磁気自由度を含めたスピン模型、連続自由度を含む有効モデルへ拡張することになる。重要なのは、どの自由度を「確率変数」に含め、どの自由度を「エネルギー関数」に吸収したのかを明確にすることである。
6.3 Kinetic Monte Carlo
拡散、析出、相変態の核生成、表面成長などの非平衡過程は、時間発展そのものが重要であり、平衡MCでは扱えない。Kinetic Monte Carlo(kMC)は、遷移率に基づく連続時間マルコフ過程として、イベント駆動で時間を進める方法である。状態
で確率分布が時間発展すると考える。kMCは、この確率過程の標本過程(軌道)を生成し、平均挙動や分布を推定する。
具体的には、状態
を定義すると、次に起こるイベントの種類は確率
に従う。したがって一様乱数
と生成すればよい。これがBortz–Kalos–Lebowitz(n-fold way)として知られる基本であり、イベント数が多い場合でも「実際に起こるイベントだけで時間を進める」ため、拡散のように長い時間スケールへ到達できる。
レート
が用いられ、
を満たすようにレートを設計する必要がある。これが満たされないと、長時間後に誤った定常状態へ落ちる。
材料科学のkMCでは、(i) どのイベント集合を採用するか、(ii) レートが局所環境でどう変わるか、(iii) 長距離弾性場や電場などの非局所効果をどう扱うか、が中心課題になる。kMCは「時間を持つMC」であるが、時間はレートモデルから生まれるため、モデルの妥当性がそのまま時間スケールの妥当性になる。
6.4 粒界、欠陥、成長過程
多結晶材料では、粒界の移動、粒成長、再結晶などが機械特性や輸送特性を大きく左右する。これらは幾何学(曲率、トポロジー)とエネルギー(界面エネルギー、欠陥相互作用)が強く結び付いた問題であり、MCは粗視化モデルとして有効である。代表はPotts模型を用いた粒成長モデルであり、格子点
とする。隣接点のラベルが異なるとエネルギーを持つため、粒界長(2次元)または粒界面積(3次元)に比例するエネルギーが得られる。ここでの
更新は局所的な再配向として実装される。ある格子点を選び、隣接点のラベルのいずれかへ変更を提案し、エネルギー差
で受理する。
欠陥集団(空孔クラスター、析出核、転位密度場の粗視化)も、確率変数の選び方次第でMCに落とし込める。例えば、格子上の占有変数
さらに、連続体場(例えば相場
でサンプリングすることに相当し、界面粗視化や熱揺らぎを含む相場の分布を求めるのに使われる。フェーズフィールドと併用する場合は、動力学(時間発展)と統計(分布)の役割分担を明確にする必要がある。
6.5 量子統計のモンテカルロ
軽元素、低温、あるいは強い量子揺らぎが支配する材料では、古典的なボルツマン重みだけでは不十分であり、量子統計を扱う必要がある。量子系の分配関数は
であり、状態和が波動関数ではなく演算子のトレースになっている点が古典と異なる。経路積分(path integral)に基づくモンテカルロは、この量子分配関数を「高次元の古典統計問題」へ写像してサンプリングする方法である。
基本はTrotter分解であり、
と近似する。
材料科学で出番があるのは、ゼロ点振動、同位体効果、水素拡散、軽元素の相安定などである。例えば、古典近似では振動の基底状態エネルギーがゼロと見なされるが、量子ではゼロ点エネルギー
一方で、量子モンテカルロには困難もある。フェルミオンの符号問題により、一般の電子系では統計誤差が指数的に増大しやすい。材料科学で電子相関を量子MCで扱う場合は、模型ハミルトニアン(Hubbard模型など)や、符号問題が軽減される特殊条件に限定されることが多い。したがって、どの自由度(核か電子か)を量子的に扱い、どの近似(ボルン–オッペンハイマー、調和近似など)を採用しているかを明確にする必要がある。
6.6 他手法との連携
MCは平衡分布と揺らぎを扱うのが得意であり、DFTやMD、フェーズフィールド、有限要素法と役割が異なる。連携の要点は、(i) ミクロから有効ハミルトニアンを作る、(ii) MCで有限温度平均を取る、(iii) 必要ならさらに粗視化して連続体へ渡す、という段階である。ここでは、代表的な結び方を「何を入出力として受け渡すか」という観点で整理する。
第一に、DFTとMCの連携である。DFTは
第二に、MDとMCの補完である。MDは時間発展(運動方程式)を持つため動的現象に強いが、長時間平衡化が難しい場合がある。MCは時間の解釈を持たない代わりに分布の探索に集中できるため、平衡構造探索や相共存の評価に向く。例えば、MDで局所的な緩和や振動特性を調べ、MCで配置自由度(置換、欠陥、方位)を熱平衡化する、といった役割分担が成立する。kMCはこの中間に位置し、レートモデルが用意できるなら長時間の拡散・析出を扱える。
第三に、連続体手法への受け渡しである。フェーズフィールドは自由エネルギー汎関数に基づく時間発展を扱い、有限要素法は複雑形状の力学場や拡散場を解く。これらへMCから渡す量は、自由エネルギー密度、界面エネルギー、移動度、相図(化学ポテンシャルと濃度の関係)などの粗視化パラメータである。例えば、MCで求めた平衡組成や化学ポテンシャル曲線を、フェーズフィールドの自由エネルギー関数形にフィットし、相分離ダイナミクスへつなぐことができる。重要なのは、粗視化で失われる自由度(短距離秩序など)が、連続体パラメータにどの程度埋め込まれているかを理解することである。
最後に、全体としての確認点は「同じ物理量を二通りで評価して整合を見る」ことである。例えば、相境界をMCで求めたなら、別の自由エネルギー推定(熱力学積分、再重み付け)でも一致するかを確かめる。連携は便利であるが、受け渡し量の定義が曖昧だと誤差が増幅されるため、何を平均し、何を固定し、どの集団の量として定義したかを明示する姿勢が必要である。
小まとめ
第5章では、統計集団の分布をサンプリングするという立場から、詳細釣り合いに基づくMetropolis–Hastings、自己相関を含む誤差評価、臨界減速を抑えるクラスタ更新、障壁越えのための拡張アンサンブル、自由エネルギー差の推定法を数式で組み立てた。第6章では、格子模型と秩序変数、クラスター展開による第一原理エネルギーの圧縮とMCによる有限温度相図、kMCによる非平衡時間発展、粒界・成長の粗視化モデル、量子統計への拡張、そして他手法との役割分担を整理した。今後は、MCで得られる平衡分布や自由エネルギーを、フェーズフィールドや有限要素法の連続体記述へ整合的に埋め込み、ミクロ起源を保ったままメゾ・マクロの材料現象を予測する方向へ進むことになる。
第7章 分子動力学の基礎と数値積分
本章では、原子の運動を古典力学として記述し、数値積分によって時間発展を追跡する分子動力学(MD)の骨格を式から組み立てる。MDは乱数で分布を直接サンプリングする手法ではなく、運動方程式を解いた軌道から統計量を評価する手法であるため、力学と統計の接続を意識して学ぶ必要がある。
7.1 ニュートン方程式とハミルトニアン
古典MDの出発点は、各原子
に従うという仮定である。ここで力
として与えられることが多く、MDは結局「
保存量を見通しよく扱うために、ハミルトニアン
を導入する。ハミルトン方程式は
であり、これを成分ごとに書けばニュートン方程式と同値である。実際、
となり、速度と力が自然に現れる。
ハミルトニアン形式の利点は、エネルギー保存の条件が明確になる点である。外場が時間に依存せず、数値積分が厳密解であれば
が成り立ち、全エネルギーは保存される。現実の計算では離散時間の近似誤差があるため、エネルギーは完全には一定にならないが、後述するVerlet系の積分法は「長時間での保存量の崩れ方」を抑えるように設計されている。
統計との接続としては、外部と熱交換も圧力交換もしない理想化では、系は微視的には一定エネルギー面上を運動し、微小正準集団(
7.2 Verlet系積分法
MDで支配的になる数値誤差は、時間の離散化によって生じる。Verlet積分は、テイラー展開から自然に導出でき、時間反転対称性と長時間安定性の点で基本となる。
まず、ある粒子の位置
これらを加えると奇数次の項が消え、
である。よって
を得る。
ただし位置Verletだけでは速度が直接得られないため、実務では速度Verlet(velocity Verlet)が広く用いられる。導出は同様に、速度のテイラー展開
を使い、
まずこの式で位置を更新して新しい力
で速度を更新する。力評価が1ステップに1回で済む点が計算効率上の利点である。
Verlet系が長時間安定になりやすい理由は、時間反転対称性と、離散写像が位相空間の幾何を保ちやすい(シンプレクティック性)ことに関係する。直感的には、エネルギーが単調に漂うのではなく、ある有効ハミルトニアンの周りで小さく揺らぐ形になりやすい。したがって、短時間の誤差が小さいだけでなく、長時間計算でのドリフトが抑えられる。
時間刻み
のような条件で支配される。固体では原子振動の高周波成分(軽元素の伸縮、硬い結合)により
7.3 温度制御と圧力制御
MDの軌道から温度や圧力を定義するには、まず運動エネルギーと自由度の関係を用いる。古典系で拘束がないとき、平衡での等分配則より
が成り立つ。ここで
で定義できるが、これは瞬間量であり揺らぐ。物性評価では時間平均やブロック平均で温度の統計を扱う必要がある。
温度一定の条件を模倣するには、熱浴を導入して正準集団(
を実現するのが目標である。熱浴には複数の考え方があり、どれを選ぶかで動力学の性質が変わる点を理解する必要がある。
Nosé–Hoover法は、拡張変数を導入して決定論的に熱浴効果を表現する。代表的には摩擦係数に相当する変数
のように運動方程式を拡張する(
確率的な熱浴としてはLangevin法がある。これは
の形で摩擦項とランダム力を加える。ランダム力は平均ゼロで
かつ揺らぎ散逸関係により
を満たすように与えられる。これにより長時間で正準分布に収束するが、摩擦係数
圧力一定(
(周期境界で原点の取り方に注意しつつ)と書ける。より一般には応力テンソル
の形を持つ(第一項が運動項、第二項がビリアル項である)。等方圧は
Parrinello–Rahman法はセル行列
7.4 ポテンシャルエネルギー面
古典MDの予測能力は、ポテンシャルエネルギー
で閉じているからである。したがって「どの結合をどの程度再現したいか」に応じて
経験力場は、結合伸縮・角度・ねじれ・非結合相互作用を足し合わせる形で
のように構成されることが多い。分子系ではこの分解が直感的であり、パラメータの意味が追いやすい。一方、金属では電子の局在結合というより集団的な結合が支配的であり、埋め込み原子法(EAM)のように
の形で「局所電子密度」
共有結合性が強い材料(Si、C、窒化物など)では、結合角や局所配位を反映するためにTersoff型やStillinger–Weber型のような角度依存項を含むモデルが使われる。これらは
のように2体・3体で構成され、結晶とアモルファスの違いを「角度分布」として表現できる利点がある。反応過程を含む場合は、結合の切断・生成を表現する反応力場(例:ReaxFFのような結合次数に基づく形式)が用いられるが、計算負荷が大きく、パラメータも多いため、適用条件と再現範囲を慎重に見極める必要がある。
第一原理MDは、
を用いる。Car–Parrinello MD(CPMD)では電子波動関数に仮想慣性を与えて拡張ラグランジアンを作り、電子と核を同時に時間発展させることで、毎ステップの厳密なSCFを避けるという考え方である。いずれも、電子の計算誤差が力へ直結するため、時間積分だけでなく電子収束の誤差も物理量へ反映されることを理解する必要がある。
ポテンシャルエネルギー面という言葉は、原子配置空間におけるエネルギー地形を意味する。平衡構造はその極小であり、遷移状態は鞍点である。MDはこの地形上を熱エネルギーで走る軌道であるため、観測したい現象(拡散、相転移、破壊)がどの程度の障壁を持つかを見積もると、必要な温度や時間、手法選択(加速法や別手法との接続)を論理的に決めやすい。
7.5 長距離相互作用
イオン性結晶や極性材料、荷電欠陥、界面の電気二重層などでは、クーロン相互作用
が本質的であり、単純な距離カットオフは系統誤差を生む。クーロン相互作用は
Ewald和は、長距離相互作用を実空間で急速に収束する項と逆空間(フーリエ空間)で収束する項に分けて計算する。概念は、点電荷をガウス分布で「ぼかした電荷」と「補正項」に分解することに対応する。エネルギーは形式的に
と分解され、
粒子メッシュ法(PME、PPPM)は、Ewaldの逆空間和を格子上のフーリエ変換(FFT)で加速する考え方である。電荷をメッシュへ割り当て(補間)して密度場を作り、ポアソン方程式に相当する処理をFFTで解いて逆空間寄与を得る。計算量は概ね
へ改善し、大規模系で有利である。補間次数、メッシュ間隔、実空間カットオフ、
周期境界と電荷の総和にも注意が必要である。総電荷がゼロでない系を周期境界で扱うと、数学的にエネルギーが発散するため、多くの実装では一様な中和背景を仮定して計算を成立させる。この仮定は、欠陥形成エネルギーや界面電位の解釈に影響しうるため、電荷状態を扱うときは「何を仮定して無限系を定義しているか」を理解する必要がある。同様に、双極子・四重極子の補正や、2次元周期(スラブ系)でのEwaldの派生形式など、境界条件に応じた取り扱いが存在する。
長距離相互作用はクーロンだけではない。磁気双極子相互作用、分散力(
7.6 物性評価と統計
MDの出力は軌道
が基本であり、液体やアモルファスの短距離秩序を定量化できる。散乱実験と対応する構造因子
として与えられ、実空間と逆空間の両方から構造を理解できる。
輸送係数の多くは、時間相関関数から得られる。自己拡散係数は、アインシュタイン関係により平均二乗変位(MSD)を用いて
で評価される。
と書けるため、両者が一致するかを確かめると数値の信頼性が増す。
粘性係数は応力テンソルの非対角成分(せん断応力)
と与えられる。熱伝導率は熱流
となる。ここで熱流の定義はポテンシャルの形式に依存し、多体ポテンシャルでは2体和の単純な形では書けないため、実装が与える熱流の定義を理解して使う必要がある。
有限サイズと有限時間は、MDの誤差構造を決める。周期境界下の拡散係数は、流体力学的な自己相互作用により系サイズ
のような補正が現れる(
統計誤差は、MDでも自己相関に支配される。時間相関関数の積分は、積分上限を有限に切ることによる系統誤差と、ノイズの積分による統計誤差が競合する。したがって、ブロック平均、複数独立初期条件、相関時間の見積もりを用い、誤差の減り方が想定通りかを確認しながら評価する必要がある。MDは軌道が連続であるため一見データ量が多いが、有効独立サンプル数は自己相関時間で決まる点を忘れてはならない。
第8章 分子動力学の材料応用と機械学習ポテンシャル
本章では、MDが材料のどの問題に強く、どの制約を持つかを、変形・輸送・反応という観点から整理し、近年の機械学習ポテンシャルによる拡張を理論と実装の接点で説明する。特に、第一原理の精度と大規模MDの計算量の間を埋める方法として、学習モデルと不確かさ評価が不可欠になりつつある。
8.1 変形、転位、破壊
塑性変形や破壊は、原子配列の不連続や欠陥の生成を伴うため、連続体の場だけでは直接見えにくい素過程が支配する。MDでは、外部ひずみを与えたときに転位がどこで核生成し、どの面で滑り、双晶がどう成長するかを原子配置として追跡できる。したがって、材料強度の起源を「欠陥の生成と運動」という機構へ分解して理解するのに向く。
計算では、ひずみテンソル
のように更新し、ひずみ速度
転位の同定には、バーガースベクトルや不整合の指標が必要である。結晶中の転位は、結晶学的には閉曲線に沿った変位の不一致として
で定義される。MDの離散原子配置では連続変位場
MDの制約は時間スケールにある。時間刻みが
加速MDの考え方は、遷移状態理論(TST)と結び付く。ある盆地からの脱出レートが
で与えられるとき、
破壊については、原子結合の切断や空孔成長、き裂先端での原子再配列が中心となる。経験力場では化学結合の切断が表現できない場合があるため、反応力場や第一原理MD、あるいは機械学習ポテンシャルの導入が重要になる。破壊は局所現象である一方で応力場は長距離に及ぶため、境界条件とセルサイズの影響が大きく、連続体との対応付け(応力拡大係数やJ積分の概念)を意識すると議論が整理される。
8.2 熱輸送と非平衡MD
熱伝導は、電子寄与と格子寄与に分けて考えることが多いが、古典MDが直接扱えるのは主に格子(原子振動)による熱輸送である。熱伝導率
で定義され、熱流
平衡MD(EMD)では、前章のグリーン・久保関係
により評価する。これは線形応答理論の結果であり、外部駆動を与えずに内部の揺らぎだけから輸送係数が得られる点が美しい。ただし、積分の収束は遅く、長時間尾部や有限サイズの長波長フォノンに敏感である。特に高熱伝導材料では自己相関の減衰が遅く、統計誤差の扱いが難しくなるため、ブロック平均と積分上限の選び方が重要である。
非平衡MD(NEMD)では、系の一端を高温、他端を低温に保つなどして定常温度勾配を作り、定常熱流から
である。温度は局所運動エネルギーから定義されるが、局所平衡が成立する程度の空間分割が必要である。熱浴領域の厚さが薄すぎると境界で非物理な散乱が起き、厚すぎると有効輸送領域が短くなるため、設計は物理と数値の両面に依存する。
Müller–Plathe法は、熱流を直接制御する逆NEMDの一種である。系の高温側と低温側で原子の速度を交換し、外から与えた運動エネルギー交換量により熱流を既知にする。交換によって生じた温度勾配を測り、
で求める。熱流
界面熱抵抗(Kapitza抵抗)は、異種材料界面での温度不連続として観測され、
で定義される。MDでは、界面を含む系に熱流を流して
熱輸送では、サイズ効果が物理そのものとして現れる。フォノンの平均自由行程が系サイズに近いと、散乱が境界で支配され、バルク熱伝導率とは異なる有効値が得られる。したがって、MDで得た
8.3 第一原理MDと反応過程
化学反応、電荷移動、結合の組み替えが本質である材料現象では、固定の古典ポテンシャルでは不十分になり、第一原理MDが必要になる。第一原理MDでは、原子配置
反応過程や拡散の議論では、自由エネルギーが中心となる。反応座標(集団変数)
と自由エネルギー曲線(平均力ポテンシャル)を定義できる。ここで
で与えられ、
傘サンプリングは、バイアス
でサンプリングし、後で再重み付けして
を用いれば、狙った
メタダイナミクスは、探索済み領域へ時間とともに反発バイアスを積み上げ、自由エネルギー障壁を乗り越えさせる方法である。集団変数
のようにガウス関数で加え、系が同じ領域に留まることを避ける。十分長時間後に、適切な条件では
第一原理MDは計算量が大きく、時間スケールの制約がさらに厳しい。したがって、反応・拡散の障壁を議論する場合は、(i) 最小エネルギー経路(NEB等)で鞍点を求める、(ii) 調和遷移状態理論でレートを推定する、(iii) 必要に応じて自由エネルギーを傘サンプリング等で評価する、というように、目的に応じて複数の方法を組み合わせることが合理的である。MDは万能ではないが、適切な自由エネルギー記述へ接続できれば、温度依存の反応性や相安定を一貫した枠組みで議論できる。
8.4 機械学習ポテンシャル
機械学習ポテンシャルは、第一原理に近い精度と古典MDに近い速度を両立することを目標に、原子配置からエネルギーと力を予測するモデルである。基本は教師あり学習であり、DFTなどで得たデータ集合
で定義できるため、エネルギーを学習しておけば自動微分により力が得られるという構造がある。これにより、エネルギー保存則と整合した力が得られ、NVE計算での挙動が安定しやすい。
多原子系のエネルギーを直接高次元関数として学習するのは難しいため、局所性を仮定して
のように原子環境
物理的に必須なのは対称性である。エネルギーは並進・回転・同種原子の置換に対して不変でなければならない。すなわち、任意の回転
が成り立つ必要がある。一方、力はベクトルであるため回転に対して共変であり、
を満たす必要がある。E(3)等変性(並進・回転・反転を含むユークリッド群の等変性)を組み込んだグラフニューラルネット型ポテンシャルは、この要請をモデル構造として実現する。メッセージパッシングにより原子間の相互作用を伝播させ、スカラー(エネルギー)とベクトル(力)の変換則を保つように設計されるため、少量データでも汎化しやすい傾向がある。
一方で、等変性モデルだけが選択肢ではない。対称関数(Behler–Parrinello型)、SOAP記述子に基づくガウス過程回帰(GAP)、多項式近似(MTP)、ビスペクトラム(SNAP)など、記述子と回帰器の組み合わせは多様である。重要なのは、(i) 局所環境の表現が連続で微分可能であること、(ii) 対称性が満たされること、(iii) 学習データの範囲で誤差が小さいだけでなく、MDで必要になる力の精度と滑らかさが確保されることである。
学習の損失関数は、エネルギーと力を同時に合わせる形が多い:
ここで重み
8.5 外挿と不確かさ
機械学習ポテンシャルの最大の弱点は、訓練データの外側での誤りである。学習は経験的リスク最小化であり、未知領域での誤差はデータが保証しない。MDは時間発展により状態空間を歩くため、初期状態が訓練領域にあっても、熱揺らぎや変形で外側へ出る可能性がある。このとき、わずかな力の誤りが軌道を大きく変え、最終的に非物理な構造へ進むことがある。
不確かさを定量化する代表はアンサンブルである。異なる初期値や異なるデータサブセットで学習した複数モデル
を不確かさの指標として用いる。力についても成分ごとに分散を計算でき、分散が大きい局所環境は訓練不足である可能性が高い。ガウス過程に基づくモデルでは予測分散が理論的に得られるため、同様の判断に使える。
不確かさの指標は、学習ループではなくMD運用の安全装置として機能する。例えば、NVEでの全エネルギーが短時間で崩れる、局所温度が異常に上昇する、原子間距離が非物理な値へ入るといった挙動は、外挿による破綻の兆候になりうる。これらは不確かさ推定と併用すると強く、分散が大きい状態に入った瞬間に計算を止め、追加の第一原理計算でデータを補強するという運用が可能になる。
外挿検知のもう一つの考え方は、局所環境の距離を定義し、訓練集合との類似度を測ることである。記述子空間で最近傍距離を測り、ある閾値を超える状態を未知領域とみなす。これはモデル非依存に近いが、記述子の選び方で結果が変わるため、物理的に意味のある表現(回転不変で連続な記述子)を用いる必要がある。
不確かさを含む設計は、材料応用に直結する。例えば破壊や化学反応は、平衡近傍とは異なる高エネルギー・低配位の状態が現れやすい。そこは訓練データが不足しがちであり、外挿が起きやすい領域である。したがって、目的が破壊であれば、き裂先端の低配位状態や引張下の遷移状態を意図的にデータへ含める必要がある。同様に高温液体やアモルファスを扱うなら、広い温度・密度のサンプルが必要になる。このように、データ設計と外挿検知は一体であり、学習精度の数字だけでは安全性は判断できない。
8.6 マルチスケール接続
MDで得られる量を連続体モデルへ渡すと、実験スケールのプロセスやデバイス条件に接続できる。重要なのは、渡す物理量が連続体方程式の中でどの係数として現れるかを明確にすることである。例えば拡散であれば、連続体では濃度場
に従うとモデル化され、MDから得た
界面現象では、界面エネルギー
(
弾性は有限要素法と接続しやすい。MDから弾性定数
を測るか、平衡揺らぎから弾性定数を推定する方法を用いる。結晶の異方性弾性を得られれば、FEM側の構成式へそのまま埋め込める。欠陥周りの局所弾性や界面の弾性不連続など、連続体では与えにくい情報をMDが補う形になる。
熱輸送では、MDで得た
の係数として組み込める。ナノ複合材料では、バルク
を与える形の接続が現実的である。ここでも、単位系(W/mK、m(^2)K/Wなど)と定義(どの温度差を界面温度差とみなすか)を明示しないと、連続体側で解釈がずれる。
機械学習ポテンシャルは、この接続を加速する。第一原理で得られる高精度情報を学習モデルで大規模へ拡張できれば、界面・欠陥・多相構造のような大系で係数を推定しやすくなる。一方で、マルチスケール接続では「平均化」が不可避であり、MDで観測される局所揺らぎや非平衡効果が、連続体の有効係数へどう写るかが難所である。したがって、係数の推定だけでなく、その係数がどのスケール・どの条件で有効かを併記する姿勢が必要である。
小まとめ
第7章では、ハミルトニアン形式からMDの運動方程式を導き、Verlet系の時間積分がテイラー展開から自然に得られること、温度・圧力制御が統計集団の実現として定式化できること、ポテンシャル選択と長距離相互作用の扱いが精度を支配すること、物性評価が時間相関と有限サイズ・有限時間誤差に敏感であることを式とともに整理した。第8章では、変形・破壊・熱輸送・反応の各問題でMDが見せる素過程と時間スケール制約を明確にし、第一原理MDによる自由エネルギー評価、機械学習ポテンシャルの対称性と学習式、不確かさ推定と外挿検知、連続体手法への係数受け渡しの論理を示した。今後は、等変性モデルと不確かさ推定を基盤として、必要な領域にデータを追加する学習戦略と、MDで得たミクロ情報を連続体の係数として矛盾なく埋め込む枠組みを発展させることで、材料の設計と予測の精度と適用範囲が拡張されていく方向にある。
第9章 フェーズフィールド法の方程式
本章では、フェーズフィールド(Phase-Field; PF)法が「界面を追いかける」のではなく「場(field)の方程式として界面を含めて解く」枠組みであることを、自由エネルギー汎関数と変分法から丁寧に導く。PFの強みは、(i) 界面のトポロジー変化(分岐・合体)が自然に扱えること、(ii) 濃度・温度・弾性・電磁場などの連成を同一のエネルギー原理で組み立てられること、(iii) 数値解法がPDEの枠組みに落ちることにある。一方で、界面幅や移動度などの「モデルパラメータ」が物理量にどう対応するかを理解しないと、結果の解釈が曖昧になる。ここでは、方程式の形がどこから来るのかを式で追い、どこに物理仮定が入っているかを明確にする。
9.1 秩序変数と自由エネルギー汎関数
PF法では、相を表す秩序変数(order parameter)
が相A(母相)、 が相B(析出相・固相など) のように、相を連続的に表す。界面は が 0 と 1 の間を滑らかに遷移する領域として表現され、界面位置を明示的に追跡しない。
秩序変数と組成場
と置く。ここで
:局所自由エネルギー密度(バルク項) : の勾配エネルギー(界面を有限幅にする) :濃度勾配に対するペナルティ(Cahn–Hilliard型で現れる) である。
なぜ勾配項が界面エネルギーを与えるのかを、1次元で確認する。
とする(
を最小化すると、オイラー=ラグランジュ方程式
を得る。両辺に
境界で
が得られる。これは「界面では勾配エネルギーと局所エネルギーが釣り合う」ことを意味する。
界面エネルギー(単位面積あたり)は、1次元の積分として
(上の釣り合い式を使った)と表せる。したがって
のように見積もれる(厳密な係数はポテンシャル形で変わる)。ここがPFの重要点で、
さらに、複数相や多相のときは秩序変数が複数になり、局所自由エネルギーは相ごとの自由エネルギーを補間する形で
のように設計することが多い。
9.2 Cahn–Hilliard方程式
組成場
は(外部から出入りがない限り)一定である。保存則の一般形は
であり、
と置く(
となる。ここまでは「保存則+線形応答」の議論であり、どんな自由エネルギーでも成り立つ枠組みである。
PFの要点は、
例えば
なら、変分計算により
を得る(勾配項の変分でラプラシアンが出る)。したがってCahn–Hilliard(CH)方程式は
となり、4階のPDEになる。4階になるのは
CH方程式が相分離(スピノーダル分解)を表現する理由を、線形安定性で見る。均一状態
とする。
定数項は勾配で消えるので、フーリエモード
が得られる。ここで
で与えられ、初期の相分離パターンの代表長さ
移動度
のように解釈できる(厳密には熱力学因子を含む)。したがって、
9.3 Allen–Cahn方程式
秩序変数
と定義する。これがAllen–Cahn(AC)方程式であり、
ACの重要な性質は「エネルギー減少性」である。
となり、
なら
である。従って
となる。界面(
(
注意点として、ACは非保存なので、保存量を表すべき変数(濃度など)にそのまま使うと物理が破綻する。濃度はCH、相指標はACという役割分担が基本である。
9.4 連成場と駆動力
PFは自由エネルギーに項を追加することで連成を統一的に扱える。例えば弾性場を加えると
のようになる。
で定義される。弾性エネルギー密度は線形弾性なら
と書ける。
弾性場は準静的とみなすことが多く、力のつり合い
を満たす
の右辺に弾性由来の駆動力が現れる。重要なのは、弾性場が「全域に瞬時に伝わる長距離相互作用」である点で、数値的には連成が硬くなり、単純な陽解法が不安定になりやすい。
電場・磁場も同様に、例えば電気エネルギー項として
のような形を入れたり、磁性なら
のように追加できる。PFの方程式は「追加したエネルギーがそのまま駆動力になる」ため、連成を設計する際はエネルギー項の次元と意味を先に整えることが重要である。
9.5 異方性と結晶方位
界面エネルギーの異方性が強い材料では、界面形状の選択(特定面方位が出やすい)や樹枝状成長の枝の向きが決まる。PFで異方性を入れる方法は大きく2系統ある。
(1) 勾配エネルギー係数を方位依存にする: 界面法線
とする。例えば結晶対称性に合わせて
のように置き、
(2) 界面移動度や界面エネルギーを、界面項で直接異方化する: いずれにせよ、異方性が入ると界面は面方位を選びやすくなり、数値計算では「メッシュが勝手に方位を作る」偽異方性(grid anisotropy)が問題になる。偽異方性を避けるには、
- メッシュを十分細かくする(界面幅に対して十分な点数)
- 等方性の良い離散化(高次要素、等方差分、スペクトル)
- 界面法線評価の精度を上げる(スムージング、正規化の安定化) などが必要になる。異方性の議論は物理だけでなく数値解析の要素が強いため、「得られた形が物理異方性か、メッシュ由来か」を必ず検証する(格子回転、格子幅変更、手法変更)姿勢が重要である。
9.6 熱ゆらぎと核生成
相変態の初期(核生成)は本質的に確率過程である。決定論的PF(CH/AC)だけでは、核生成が起きるきっかけがなく、初期ゆらぎを人工的に与えないと変態が始まらない場合がある。核生成を扱うためには、熱ゆらぎ項(ノイズ)を加える確率偏微分方程式(stochastic PF)を考える。
保存場(CH)では揺らぎも保存則を満たす必要があるため、
のように流束にノイズを入れる形になる。
のような相関を持つように定める(白色雑音の理想化)。これにより、平衡分布がボルツマン分布と整合する。
非保存場(AC)なら
のように直接ノイズを加える形が基本である。
ただし実務では、連続体ノイズをそのまま離散化するとメッシュ依存性が出やすく、核生成率の定量は難しい。PFで核生成を議論する場合は、(i) 古典核生成理論(CNT)との整合、(ii) 反応座標・自由エネルギー障壁の見積もり、(iii) MD/MCと併用して核生成率を別手法で検証、などを組み合わせるのが安全である。PFは「成長と形態選択」には強いが、「稀な核生成イベントの絶対レート」は別枠で扱う必要が出やすい、という役割分担を理解する。
第10章 フェーズフィールドによる組織形成
本章では、第9章で導いた方程式(CH/ACと連成)を用いて、材料組織がどのように形成されるかを、凝固・粒成長・析出・反応拡散などの典型問題に沿って整理する。PFの重要点は「組織を決める因子(界面エネルギー、拡散、弾性、潜熱、異方性)」が、自由エネルギーと移動度として方程式中に明示的に現れることである。つまり、パラメータを変えると何が変わるかを論理的に追える。
10.1 凝固と結晶成長
凝固は、相(固相/液相)の秩序変数
:AC型(非保存) :CH型または拡散方程式(保存) :熱拡散方程式(保存) が結び付く。
温度場は、潜熱を含めると
の形で、潜熱項が
ここで
溶質拡散も重要で、固相と液相で溶質の平衡濃度が異なる(分配係数)ため、界面で偏析が起きる。PFでは局所自由エネルギーを相ごとに用意し、
のように補間して、相平衡(化学ポテンシャル一致)を連続的に表現する。すると成長前線での溶質境界層、成長速度と拡散の競合が自動的に出る。
界面形態(セル状・樹枝状)が出る理由は、界面のわずかなゆらぎが、熱/溶質場の拡散と競合して成長しうるからである。簡単には、界面前方に溶質が溜まると局所的に凝固が抑制され、界面が不安定になる(組成過冷却に相当する)。PFはこれを場の方程式として解くため、形態選択を直接シミュレーションできる。
また、界面エネルギー異方性があると、枝の成長方向が結晶方位で固定される。異方性が弱いと枝は丸まり、強いと特定方位のファセットが現れる。ここで重要なのは、異方性の強さと界面幅・メッシュが絡むため、数値条件を変えて「枝の向きが物理的に安定か」を必ず検証することである。
10.2 粒成長と再結晶
粒成長は、粒界エネルギーによる曲率駆動で進む。最も基本的な理解は、粒界の法線速度が
に比例するという関係である(
PFでは、多結晶を表すために複数の秩序変数
のようにして、複数相(複数粒)が重なりにくいようにする。勾配項を入れれば粒界が有限幅で表現され、AC型の進化
により粒界が移動する。
再結晶では、単なる曲率駆動だけでなく、加工で蓄積したひずみエネルギー(転位密度)を減らす駆動力が中心となる。そこで内部状態変数として転位密度
のように組み込む。さらに
を併用することがある。再結晶PFはモデル自由度が増える分、パラメータ同定が難しくなるため、実験(EBSD、硬さ、XRD線幅)と対応する観測量を先に定義し、逆問題として組み立てるのが重要である。
10.3 析出、相分離、スピノーダル
析出と相分離はCH方程式が中心となるが、材料では弾性場の影響が非常に大きい。整合析出では、母相と析出相の格子定数差が整合ひずみとして蓄積し、析出物形状(球状・板状・針状)や配向を決める。これは「界面エネルギーを下げたい」という傾向と「弾性エネルギーを下げたい」という傾向の競合として理解できる。
弾性連成を入れた自由エネルギー
に対して
となり、CH方程式は
で進化する。弾性項
弾性場の解法としては、(i) FEMで
スピノーダル分解の初期波長が
10.4 反応拡散と多相系
化学反応や相変態が拡散と同時に進む場合、反応拡散系としてモデル化する。例えば濃度
のように反応項
のように置き、
多相系(3相以上)では秩序変数が増え、しばしば「総和が1」の制約が必要になる。例えば相分率
を保ちたい。これを実現する方法としては、
- ラグランジュ未定乗数で制約を入れる
- 罰則項(ペナルティ)を自由エネルギーに追加する
- 独立変数を
個に減らして最後を で定義する などがある。どの方法でも、数値的に制約が崩れると「存在しない相が負になる」などの非物理が出るため、離散化・時間積分の段階で制約を保つ工夫が必要になる。
さらに多相では三重点(トリプルジャンクション)の角度条件が重要になる。等方界面エネルギーならYoungの関係により、三重点で界面張力が釣り合う角度が決まる。PFでは界面エネルギーを正しく入れると、三重点の角度が自動的に再現されるが、界面幅が大きすぎると三重点の力学がぼやけるため、界面幅とメッシュ分解能の関係を確認する必要がある。
10.5 パラメータ同定とデータ同化
PFモデルは「方程式は美しいが、係数が不明」という逆問題を必ず伴う。代表的なパラメータは
- 界面エネルギー
( に対応) - 拡散係数
または移動度 - 界面移動度(ACの
に対応) - 弾性定数
、固有ひずみ - 潜熱、熱伝導率
などである。
まず「観測量」と「モデル量」の対応を明確にする。例えば界面エネルギーは、1次元平衡界面から
で与えられるので、界面幅
逆問題としては、観測データ
の形に落とし込むのが基本である。
- まず独立に決められるもの(弾性定数、拡散)を外部から固定
- 次に界面物性(
)を合わせる - 最後に動力学係数(
)を時系列データで合わせる という段階的同定が現実的である。
データ同化の観点では、PFの状態(
10.6 ソフトウェア基盤
PFはPDEであるため、実装では
- 空間離散(FEM/FVM/差分/スペクトル)
- 時間積分(陰解法、半陰解法、分割法)
- 疎行列ソルバ(Krylov、マルチグリッド、前処理) が不可欠になる。教育段階では差分やFFTのコードで全体像を掴み、研究段階ではFEM基盤や大規模ソルバを活用して連成・複雑形状へ拡張するのが自然な流れである。
MOOSEのPhase Fieldモジュールのような基盤は、弱形式・非線形ソルバ・自動微分・並列ソルバが統合されており、「方程式をどう書くか」が実装に直結している。重要なのは、ソフトの使い方そのものよりも、
- 連立方程式がどの弱形式に落ちているか
- どの項が陰的に、どの項が陽的に扱われているか
- 収束判定(非線形残差、線形残差)が物理的に十分か を読めるようになることである。PFは見た目の結果がそれらしく出やすい反面、ソルバ未収束や偽拡散で「もっともらしい誤り」になりやすい。したがって、実装と方程式の対応を常に確認する姿勢が、研究での信頼性を支える。
章末の要点(理解のチェック用)
- PFの核心は、自由エネルギー汎関数
を立て、 を駆動力として進化方程式を作る点にある。 - CHは保存則から導かれ、
、 により4階PDEとなる。 - ACは勾配降下であり、
は を保証する。 - 弾性・熱・電磁場などは
に項を追加することで統一的に連成できるが、数値的剛性が増す。 - 異方性は物理異方性とメッシュ偽異方性を区別し、格子回転やメッシュ変更で検証する必要がある。
- 核生成は確率過程であり、ノイズPFや他手法との併用で「どこまでをPFで扱うか」を整理する。
第11章 有限要素法の理論
本章では、有限要素法(Finite Element Method; FEM)が、偏微分方程式を弱形式へ移し、近似空間上で残差を最小化することで数値解を得る方法であることを、式変形から順に示す。材料計算でFEMが強い理由は、複雑形状・異方性・連成・境界条件を、同一の定式化で扱える点にある。
11.1 弱形式と重み付き残差法
有限要素法の出発点は、強形式(微分方程式そのもの)をそのまま点ごとに満たすのではなく、残差を重み付きで平均してゼロにするという発想である。具体例として、拡散方程式(ポアソン方程式の定常版)を考える。 領域
を満たし、境界
ここで
強形式を満たさないときの残差を
と定義し、任意の重み関数(試験関数)
を課す。これが重み付き残差法の基本である。FEMでは、後に
上式を弱形式に変形するため、部分積分(発散定理)を使う。以下の恒等式を用いる:
よって
となる。境界積分は
である。ここでディリクレ境界上では試験関数に
である。
この変形には大きな意味がある。第一に、強形式では
材料力学でも同様である。小ひずみ線形弾性では、釣合式
と、構成式
を用いる。試験関数
が得られる。これが仮想仕事の原理であり、FEMの材料力学がエネルギー原理と直結する入口である。
11.2 形状関数と要素
弱形式が得られたら、次は未知関数
である。
と置く。これにより連続問題が代数方程式へ落ちる。
一次要素(線形要素)では形状関数は一次多項式であり、計算が軽く頑健である。一方、応力集中や曲率の大きい場では高次要素が有利になる。高次要素は同じメッシュでも自由度が増え、滑らかな解に対して高精度になりやすいが、数値積分や要素歪みの影響が表に出やすい。
要素の形としては、2次元なら三角形・四角形、3次元なら四面体・六面体が基本である。一般に
- 四面体は複雑形状へ自動メッシュ生成がしやすい
- 六面体は同自由度あたりの精度が高くなりやすい という傾向があるが、最終的には目的の場(境界層、界面、き裂)に合わせて選ぶ必要がある。
代表的な要素の整理を示す。
| 要素 | 自由度配置 | 長所 | 注意点 |
|---|---|---|---|
| 三角形一次(P1) | 3節点で線形 | 自動メッシュに強い | 曲げや応力で硬くなりやすい場合がある |
| 四角形一次(Q1) | 4節点で双線形 | 構造格子で効率が良い | 形状が歪むと精度が落ちやすい |
| 四面体一次(P1) | 4節点で線形 | 3D複雑形状に強い | 低次では誤差が残りやすい |
| 六面体一次(Q1) | 8節点で三線形 | 効率と精度が良い | メッシュ生成が難しい場合がある |
| 二次要素(P2/Q2) | 中間節点を追加 | 滑らかな場に高精度 | 積分・要素品質の影響が大きい |
さらに、不連続ガラーキン(DG)では要素間で連続性を課さず、数値流束で結合する。DGは移流支配問題や不連続解に強いが、自由度が増えやすく、安定化項の設計が重要になる。材料計算でDGを使う場面は、衝撃波や破壊の不連続、あるいは高Péclet数の輸送などで現れる。
形状関数は実装上、参照要素(reference element)上で定義し、実要素へ写像する。参照座標
が現れる。勾配は
で変換され、積分は
となる。異方性材料や大変形では、この写像とヤコビアンが数値安定性の中心となる。
11.3 構成式と非線形性
材料の構成式は、非線形性の主要な供給源である。弾塑性、粘弾性、損傷、相変態などでは、応力がひずみだけでなく内部変数にも依存する。 一般形として
と書き、
非線形FEMでは、弱形式から得られる残差方程式
をニュートン法で解くのが基本である。反復
である。ここで
塑性の基本例として、降伏関数
と書く(
が課される。数値的には、各時間ステップで試行応力を計算し、降伏していれば戻し写像(return mapping)で内部変数を更新する。
ここで重要なのがアルゴリズム的一貫接線(consistent tangent)である。これは、内部変数更新を含めた離散化のもとで
となる有効接線
大ひずみ問題では、回転の取り扱い(客観性)や応力率の選択が重要になる。例えば変形勾配
(
から弱形式を構成する。材料モデルの式そのものが大きく変わるため、線形弾性の直感をそのまま持ち込まないことが重要である。
11.4 時間依存問題と安定化
拡散、粘性、波動、電磁場などは時間依存問題である。空間をFEMで離散化すると、一般に
のような常微分方程式(または微分代数方程式)になる。
拡散方程式の単純な例として
を弱形式化して離散化すると
となる。ここで
である。陽的オイラーを使うと
であり、安定条件が厳しい。陰的オイラーなら
となり、拡散支配では一般に無条件安定であるが、毎ステップ線形系を解く必要がある。材料計算では拡散係数の温度依存や強不連続があり、陰解法と反復解法の組み合わせが基本になる。
移流項が強い輸送(高Péclet数)では、標準ガラーキン法が数値振動を起こしやすい。移流拡散方程式
に対し、離散化で非物理振動が出るのは、移流が一方向情報伝播であるのに対して、標準ガラーキンが中心差分的に振る舞うためである。これを抑えるため、SUPG(Streamline Upwind/Petrov-Galerkin)などの安定化を導入する。SUPGでは試験関数を修正し
のように流線方向へバイアスを与える。
波動問題では質量行列が本質的である。例えば弾性波なら
を離散化して
となる。時間積分はNewmark法などが使われ、数値分散と減衰の制御が重要になる。材料の高周波応答や衝撃問題では、安定性だけでなく位相精度(波の速度の再現)を意識する必要がある。
11.5 メッシュと幾何
FEMの精度は、要素次数だけでなくメッシュ品質に強く依存する。要素がつぶれている、角度が極端に小さい、アスペクト比が大きいなどの場合、ヤコビアンが悪条件になり、勾配評価が不安定になる。特に異方性材料や大変形では、要素品質がそのまま収束性へ影響する。
複雑形状ではCAD形状とメッシュの整合が課題になる。境界が曲線・曲面である場合、低次要素では幾何誤差が支配的になりやすい。等価的に、解を近似する前に形状自体を粗く近似しているためである。曲面境界が重要な問題では、幾何の近似次数を上げる(高次要素、アイソパラメトリック要素)ことが有効である。
局所的に急峻な場が現れる領域では局所細分化が必要になる。き裂先端、界面端、転位に相当する特異場では、解が滑らかでないため高次要素よりも細分化が効く場合が多い。適応メッシュの考え方は、誤差推定子
- 解を計算する
- 各要素の誤差指標
を評価する が大きい要素を細分化する を繰り返すことである。材料計算では、応力やひずみの二次量、あるいはエネルギー密度の勾配が指標として使われることが多い。
誤差推定には大別して2系統ある。残差型推定は、離散解を強形式に代入した残差を使い、要素内部残差と要素境界でのフラックス不連続を合成して評価する。回復型推定は、要素内で不連続になりやすい勾配(応力など)を滑らかに再構成し、その差を誤差とみなす。どちらも、厳密な上界評価が得られる場合と、経験的に有効な場合があるため、目的量(例えば最大応力、変位、J積分)に対して妥当な推定になっているかを確認する必要がある。
11.6 誤差評価と収束
FEMの誤差は、メッシュ幅
のように減少する(ノルムの種類により指数が変わる)。勾配を含むエネルギーノルムでは
となることが多い。材料力学ではエネルギーに基づく誤差評価が自然であり、内部仕事の差として整理できる利点がある。
重要なのは、誤差が場そのものではなく、目的の物性値へどう伝播するかである。例えば応力集中の最大値は、局所的な特異性に敏感で、メッシュを細かくしても収束が遅い場合がある。逆に、全体平均のエネルギーや平均ひずみは比較的早く収束する。したがって、何を予測したいかを明確にし、その量に対してメッシュ収束を検証する姿勢が必要になる。
目的量
- 主問題で
を解く - 双対問題で
を解く - 残差と双対解から
を推定する という流れになる。材料同定や最適化では、目的量に対する精度が最重要であるため、この考え方が直接役に立つ。
第12章 有限要素法による連成材料計算
本章では、FEMが複数の物理場を同一の弱形式枠組みで結合できる点を活かし、熱・応力・拡散、電磁場、フェーズフィールドなどの連成をどのように組み立てるかを示す。連成では方程式が増えるだけでなく、時間スケールや支配パラメータが異なる場が結び付くため、解法設計と検証が材料科学的な結論の信頼性を左右する。
12.1 熱・応力・拡散の連成
熱と応力の連成の基本は熱膨張である。温度変化
とすると、線形弾性の構成式は
となる。ここで
であり、
拡散と応力の連成では、化学ポテンシャルに応力項が入る点が要である。組成
と分け、体積変化を伴う溶質(化学膨張)では
のように静水圧応力
となり、応力勾配が拡散を駆動しうる。これが応力誘起拡散や応力腐食割れの基本機構の一部になる。
弱形式では、変位
のようになる。ここで
連成解法にはモノリシック解法と分割解法がある。モノリシックでは未知ベクトルを
のようにまとめ、非線形方程式
| 解法 | 特徴 | 注意点 |
|---|---|---|
| モノリシック | 強連成でも頑健になりやすい | ヤコビアン構築と前処理が難しい |
| 分割 | 実装が分かりやすい | 強連成で不安定になりやすく反復が必要 |
材料プロセスでは温度と変形が強く結び付く場合が多く、分割解法では反復(固定点反復)を入れて収束を確保する設計が必要になる。どちらを選ぶかは、連成の強さ、時間刻み、計算規模、利用できるソルバ基盤で決めるのが現実的である。
12.2 電磁場と材料応答
電磁場連成ではマクスウェル方程式が出発点である。時間領域の基本形は
および
である。材料応答は構成式
で入る。導電率
渦電流問題(準静的近似)では、周波数が高すぎない範囲で変位電流を無視し、
とし、ベクトルポテンシャル
と置く。ゲージ条件を適切に選ぶと、未知数は
となる。ここで物性が複素化する意味が明確になる。例えば導電体では、電流が位相遅れを持ち、損失(ジュール熱)を生む。損失密度は
であり、調和場では時間平均が
で与えられる。FEMはこのような複素場の弱形式にも自然に対応する。
電磁場と材料変形の連成としては、ローレンツ力
が力学方程式の外力項に入る。これにより電磁成形や渦電流ブレーキが記述できる。磁歪などの磁気弾性連成では、自由エネルギーに磁気エネルギーと弾性エネルギーの結合項を入れ、磁化(または磁気ポテンシャル)と変位を同時に解く設計になる。周波数領域での複素透磁率は、損失を含む実効応答を表しており、実部が蓄積、虚部が散逸に対応する。材料計算では、どの物理を複素物性へ押し込めているのかを明確にする姿勢が重要である。
12.3 フェーズフィールドと有限要素法の統合
フェーズフィールド(PF)方程式は弱形式に落とせるため、FEMと自然に統合できる。例えばAllen–Cahn型
を、自由エネルギー
で定めると、化学ポテンシャルに対応する変分は
である。試験関数
右辺のラプラシアン項を部分積分して
となる。境界条件により境界積分を整理すれば、FEM離散化は通常の手順で行える。これによりPFを複雑形状へ適用しやすくなる。
Cahn–Hilliard型は4階PDEであるため、扱い方に工夫が要る。方法は2つある。 1つ目は
とし、未知数を
弾性場との連成は、PFの自由エネルギーに弾性項を加えて行う。例えば
である。すると
ここで数値誤差の扱いが重要になる。PFは界面幅
12.4 逆問題と最適化
材料定数同定や形状最適化は、観測量と設計変数を結ぶ最適化問題として定式化できる。設計変数を
のように定める。
問題は、
離散化後の支配方程式を
と書き、ラグランジュ関数
を導入する。
を満たす
で計算できる。ここで重要なのは、随伴方程式を1回解けば、設計変数の次元に依らず勾配が得られる点である。材料同定(多パラメータ)や形状最適化(多数自由度)でこの差は決定的である。
ただし、随伴法は目的関数の定義に強く依存する。例えば変位の誤差を最小化するのか、応力の誤差を最小化するのか、固有振動数を合わせるのかで随伴方程式が変わる。実験が提供する観測量と、モデルが出す量を同一の意味に合わせたうえで目的関数を設計する必要がある。
12.5 大規模計算基盤
連成材料計算では自由度が大きくなり、疎行列線形系・非線形反復・時間積分を効率良く回す基盤が不可欠である。離散化で得られる線形系は
の形であり、
のようになる。ここで重要なのは、単にKrylov法を使うだけでなく、ブロック構造を意識した前処理が収束を支配する点である。例えば力学ブロックにはAMG(代数マルチグリッド)、拡散ブロックにはマルチグリッド、強連成にはSchur補完に基づく前処理を用いるなど、物理に合わせた設計が必要になる。
時間積分も計算規模に直結する。陰解法は安定であるが、各ステップで大規模線形系を解く。適応時間刻みは、硬い過渡現象を小刻みに解き、緩やかな領域を大刻みに解くことで総計算量を下げる。ただし、連成では時間刻みが場ごとに異なる要求を持つため、単一の誤差指標で刻みを決めると非効率になりやすい。目的量や残差の支配項を見ながら、刻み制御の指標を設計する必要がある。
並列計算では、メッシュ分割と通信が鍵になる。有限要素の行列組立とソルバは、分散メモリ並列(MPI)を前提に設計されることが多い。連成問題では未知量が増えるため、単に要素数を増やすよりも、1節点あたりの未知数増加がメモリを圧迫しやすい。したがって、自由度見積もり、前処理メモリ、チェックポイント(途中保存)の設計を含めて、計算資源と解法を同時に考える必要がある。
12.6 相互検証とモデルの信頼性
連成材料計算は自由度が多く、同じ観測量に対して複数のモデルが合うことがある。したがって、単に計算が回ったことをもって結論とせず、数値と物理の両面から信頼性を積み上げる必要がある。
数値面では、メッシュ収束・時間刻み収束・ソルバ収束(残差と許容値)を分けて検証する。特に非線形問題では、見かけ上の収束(残差が下がった)が、物理量の収束(応力や相分率が安定)と一致しない場合がある。ゆえに、残差だけでなく保存則やエネルギー収支のチェックが重要になる。例えば熱連成なら、外部から入った熱と内部エネルギー変化、散逸がバランスしているかを確認する。
物理面では、複数の離散化・複数のコードで照合する姿勢が有効である。例えば、同じ問題を差分法とFEMで解き、主要な物理量が一致するかを確認する。あるいは、単純形状・単純境界で解析解がある問題(1D拡散、片持ち梁など)に落として一致を確認する。このような検証は、連成が複雑になるほど価値が増す。
さらに、スケールの切り出しも重要である。FEMは連続体であるため、原子スケール現象を有効物性に押し込めている。押し込めた物理が温度・ひずみ速度・組織で変化する場合、モデルの適用範囲が狭くなる。したがって、どの物性を定数として扱い、どれを場依存にするかを明示し、必要ならDFT・MD・実験から物性の依存性を持ち込む設計が求められる。
小まとめ
本章までで、FEMが弱形式に基づく統一的な枠組みであり、材料の非線形構成式や複数物理場の連成を、同じ数理構造として扱えることを示した。特に、連成問題の成否は、方程式の立て方だけでなく、ソルバ・前処理・収束判定・誤差評価・相互検証の設計に強く依存することが分かる。
今後の発展としては、PF-FEMのような場の連成に加え、逆問題・最適化・不確かさ評価がより重要になる。材料設計では、予測値そのものだけでなく、どの仮定が結果を支配しているか、どの不確かさが支配的かを定量化し、実験や下位スケール計算へ戻して更新する循環を構築することが鍵になる。FEMはその循環の中核となる計算器として、物理モデルと数値解法の両方を理解して使いこなすことが求められる。
第13章 Pymatgenの活用方法I
本章では、Pymatgenを材料データの共通言語として位置づけ、研究室内に散在する計算・実験データと、公開データベース由来のデータを、同じ表現に揃えて扱えるようになることを目的とする。材料データがビッグデータとして価値を持つのは、量が多いからではなく、同じ定義・同じ単位・同じ最小単位で比較できる状態に整えられているからである。Pymatgenは、構造・組成・エントリ(計算結果の要約)といった中核概念を、再利用可能なオブジェクトとして実装しており、データ選択・収集の入口で強い効果を発揮する。
13.1 Pymatgenが提供するデータモデル
Pymatgenの最重要な価値は、結晶構造をStructure、組成をComposition、元素をElementとして表現し、材料データを機械が扱える形で統一できる点にある。例えば、化学式は見かけ上同じでも、還元式・整数比・サイト占有・欠陥によって表現差が混入しやすいが、Compositionを介すことで、同一視すべきものと区別すべきものを規則として扱えるようになる。これは、後段の統計解析や学習で、表記揺れを誤って規則として学習する事態を抑える役割を持つ。
さらに、構造を単に原子座標の集合として扱うのではなく、格子(lattice)、周期境界、対称性、近接環境といった、材料物性と直結する情報を同一オブジェクトに保持できる点が重要である。結晶系では、局所環境の違いが相安定性や電子状態に影響するため、後で特徴量を設計する前段階として、まず構造表現が破綻しないことが必要である。Pymatgenのコアは、この壊れにくい表現に揃えることに照準が合っている。
13.2 ローカルデータを読む・書く
材料データの収集は、公開DBからの取得だけでは完結しない。研究室内には、VASPのPOSCAR/CONTCAR、CIF、XYZ、各種入出力ファイルなどが混在しており、まずそれらを同じStructureへ変換できることが出発点になる。Pymatgenは多様な形式の読み書きを支援し、from_fileとtoに相当する流れで、構造データを共通表現へ揃える設計になっている。
この段階で重要なのは、構造だけを集めて満足しないことである。後でどの条件で得た構造かを追跡できなければ、同じ化学式の別相や、同一相の格子定数違いを混同しやすくなる。したがって、ファイル由来(計算条件、実験条件、試料ID、作成日、解析バージョン)を、構造と結びついたメタ情報として保管する設計が必要である。Pymatgenは構造操作だけでなく解析モジュールも広く持つため、読み込み直後に一貫した基準で検査(格子の妥当性、サイト数、占有、元素種)を行う癖を付けると、後段の再解析が大幅に安定する。
13.3 Materials Projectからの取得
公開DBは多数あるが、材料科学で最も一般的な入口の一つがMaterials Projectである。現在、Materials Projectは公式Pythonクライアントとしてmp-apiパッケージのMPResterを推奨しており、プロパティ条件による検索や、ID指定での取得が行える。特に、summary(要約)データは広い物性空間の検索に適しており、まず候補を絞ってから詳細を取りに行くという構成が取りやすい。
ここで注意すべきなのは、APIの世代である。Materials Projectの旧APIは推奨されておらず、過去の互換性のために残っていたが、公式ドキュメントでは旧APIは2025年9月に停止予定であると明記されている。再現性目的で旧APIへ依存する場合は、取得日・DBバージョン・クライアントの版を明示し、同じ問い合わせが再現できる形で記録する必要がある。
また、ビッグデータから必要情報だけを集めるには、最初から何を要らないと判断するかを決めることが効く。例えば、最初の候補選別では、構造(Structure)、化学式、空間群、バンドギャップ、安定性指標(例:hullからの距離)、主要なエネルギー量など、後段の選別に必要な最小限に絞る設計が基本である。Materials Projectは検索用のフィルタ条件を用意しているため、最初はsummaryのフィールドで絞り込み、必要になった候補だけ追加のエンドポイントで深掘りするのが合理的である。
13.4 OPTIMADEを使った横断検索
材料データベースを横断して収集する際、各DB固有のAPIを個別に覚えると、収集規則がDBごとに変わって比較が壊れやすい。OPTIMADEは、材料DBを共通のREST形式で検索するための仕様であり、同じフィルタ文で複数のOPTIMADE対応DBを検索できることを狙っている。PymatgenにはOPTIMADEクライアント(OptimadeRester)や、OPTIMADE表現とStructureを相互変換する仕組みがあり、横断検索からPymatgen表現への収束が可能である。
OPTIMADEが役立つのは、収集対象が結晶構造と基本メタ情報に近いときである。DBごとの高度な計算物性(特殊な分光シミュレーション、特定ワークフロー出力)まで統一して取得することは難しいが、候補材料の構造集合を広く集め、後で同じ規則でフィルタしていく用途には向く。したがって、まずOPTIMADEで候補の母集団を作り、その後にMaterials ProjectやNOMADなどの個別DBで詳細プロパティを補う、という二段構えが取りやすい。
13.5 国内データベースとPymatgenの接続
国内にも材料データの基盤が存在し、代表例としてNIMSが整備するMaterials Data Repository(MDR)などがある。これらは、データ公開・共有・再利用を支える設計思想が明確であり、材料データを研究コミュニティの資産として扱う方向性と整合する。Pymatgen側で重要なのは、取得したデータを最終的にStructureやCompositionへ揃え、同じ検査規則で品質を確認できる状態に持ち込むことである。
国内DBは、データの種類が計算構造だけでなく実験条件・計測データ・プロセス情報を含みやすい。したがって、構造はCIF等で取り出してPymatgenへ揃えつつ、同時に条件情報(温度、雰囲気、装置、試料履歴)を別レコードとして紐づける設計が必要になる。ここで構造だけPymatgen、条件はバラバラのままにすると比較が破綻するため、構造オブジェクトとメタ情報を同一の識別子で結び続ける運用が重要である。
第14章 Pymatgenの活用方法II
本章では、収集した膨大な材料データから、研究目的に合致する情報だけを選び、同じ定義で比較できるデータ集合へ整える手順を記述する。ポイントは、検索条件と選別条件を、物理・化学の量として宣言し、第三者が同じ規則で再現できる状態にすることである。そのために、Pymatgenの解析機能と、データモデル(構造・組成・エントリ)を中心に据える。
14.1 選別条件を物理量として定義する
ビッグデータからの選別は、感覚的な絞り込みではなく、条件を数式と変数で表すところから始めるべきである。例えば安定な新材料候補を集めたいなら、形成エネルギーやhullからの距離(energy above hull)に基づく条件が自然である。形成エネルギーの基本形は、元素参照エネルギー
であり、比較は通常原子あたりで揃える。ここで
また、狙いが磁性・触媒・電池などで異なるなら、選別条件も変わる。例えば磁性材料なら、磁気秩序や磁化の大きさだけでなく、結晶対称性(異方性の候補)や、安定性と合成容易性の関係も同時に見ることが多い。したがって何を最大化し、何を下限として課すかを、最初に変数として整理し、後段の多目的選別へ接続できる形にしておくべきである。
14.2 構造と組成で重複を除き、比較単位を揃える
公開DB由来のビッグデータには、重複や近似的に同じ構造が混入しやすい。例えば、同一材料でも最適化条件やセル選びの違いで格子がわずかに変わり、別エントリとして残る場合がある。ここで重複を放置すると、同じ材料が何度も学習データに現れ、汎化性能の評価が壊れやすくなる。Pymatgenは構造表現を統一できるため、比較の前に同一視の規則を決め、同じ規則でデータ集合を正規化することが重要である。
比較単位も同様に揃える必要がある。例えば、バンドギャップはeV、体積は
14.3 DB検索を二段に分け、転送量を抑える
ビッグデータ収集で最初に詰まりやすいのは、必要以上に大きなデータを一度に取得してしまうことである。Materials Projectはsummaryを中心に検索し、ID集合を得てから必要な情報だけを追加取得する設計が取りやすい。検索条件はスキーマに基づいて組み立て、どのフィールドを使ったかを記録することが再現性の中心になる。
加えて、DBの版と取得日時を必ず記録すべきである。DBは更新され、同じ材料IDでもフィールドが追加・修正され得るため、いつのDBで選別した結果かを残しておかないと、後で同じ研究判断を再現できなくなる。したがって、学習用データセットを作るときは、数値列だけを残すのではなく、参照ID、取得エンドポイント、検索条件、DB版、取得日を同時に保存し、後から同じ条件で再生成できるようにするべきである。
14.4 取得データをPymatgen表現へ収束させる
データを選ぶことと集めることは別である。選別条件を満たした材料集合が得られても、構造表現が揃っていなければ、後の特徴量設計や追加解析で破綻する。ここでPymatgenが効くのは、OPTIMADE構造、ローカルCIF、VASP出力など異なる入り口から来たデータを、最終的にStructureへ揃えられる点にある。こうして構造表現が揃うことで、空間群や局所配位の比較、距離に基づく分割設計、類似構造の探索など、後段の処理が同じ道具で実行できる。
また、材料データの再利用性を上げるには、構造だけでなく由来を同時に保持する必要がある。データベースのフィールド仕様を確認し、どのフィールドが何を意味するかを理解した上で、必要なメタ情報を落とさずに保存することが重要である。これは後で、定義揺れや参照状態の違いを見抜く際の根拠になる。
14.5 収集したビッグデータを研究目的へ接続する
最後に重要なのは、収集と選別を研究の問いに結びつけることである。例えば磁性材料の探索であれば、選別で得た候補集合に対し、結晶対称性、局所配位、元素置換の傾向などを整理し、なぜその集合が現れたのかを説明できる形へ変換する必要がある。Pymatgenは構造操作と解析モジュールを広く持つため、候補集合を同じ基準で処理し、差分を物理量として比較できる状態に整えやすい。
ビッグデータ活用で陥りやすい誤りは、検索条件を増やしすぎて、条件の意味が研究者自身にも説明できなくなることである。選別条件は増やすほど強くなるが、同時に何が効いたのかが見えにくくなる。したがって、条件は物理・化学の妥当性、再現可能性、第三者説明可能性の三つを満たす範囲で増やし、増やした条件は必ずなぜ必要かを文章で残すべきである。これにより、収集は単なる作業ではなく、学術的に追跡可能な研究手続きとして成立する。
参考文献
- Pymatgen Documentation(modules/usage/coreなど)
- Materials Project API Documentation(Using the API, Getting Started, Querying Data, Swagger UI, Database Versions, Legacy API)
- MateriApps「pymatgenを使ったMaterials Projectのデータ収集」(日本語)
- OPTIMADEと関連ツール(仕様・Python tools・チュートリアル、総説論文)
- NOMAD / OQMD / AFLOW / JARVIS などの公開DB API(横断取得戦略の検討に有用)
- NIMS Materials Data Repository(国内のデータ共有基盤の例)
第17章 まとめと展望
本書(ここまでの章)では、材料計算を「何を未知量として置くか」「どのスケールを捨てて、どのスケールを残すか」という視点で整理し、DFT・MC・MD・PF・FEMを“別々の手法”としてではなく、“同じ問題を異なる解像度で見るための言語”として学べるように構成してきた。本章では、それらを一度ひとつの地図に戻し、手法間の情報変換、信頼性の作り方、データ駆動との統合、計算機環境、研究倫理、学習順序を、初学者が自力で研究へ接続できる形にまとめ直す。
17.1 手法間の情報変換
材料計算の手法は、扱う自由度(未知量)と平均化の仕方が異なる。DFTは電子密度
この違いは、「何を平均して捨てたか」に対応する。したがって手法間の接続は、単に数値を渡すのではなく、平均化の定義(どの自由度を消したか)を明示しないと物理がずれる。典型的な“情報変換”を、式の形で整理する。
(1) DFT → 連続体材料定数(弾性・磁気・拡散など) DFTから得られるのは、エネルギー
例:弾性定数(小ひずみ) ひずみテンソル
で弾性定数が定まる。実際の計算では、ひずみを微小に振って
例:拡散(活性化エネルギー) 拡散係数はしばしばアレニウス形
で整理され、
(2) DFT → クラスター展開 → MC(合金相図) 置換合金の配置エネルギーを格子上の有効ハミルトニアンで表し、
(
(3) MD → 連続体(輸送係数・界面物性) MDは時間相関から輸送係数を与える。例えば拡散係数は平均二乗変位(Einstein関係)で
粘性係数はグリーン久保関係で
のように与えられる。ここから連続体へ渡すときは、(i) 有限サイズ補正(特に拡散)、(ii) 平衡到達と統計誤差(自己相関)、(iii) 力場/MLポテンシャルの妥当性(外挿)を同時に評価する必要がある。
(4) PF ↔ FEM(組織進化と複雑形状・連成) PFは自由エネルギー汎関数
まとめると、手法間の接続は「量を渡す」ではなく「定義を渡す」である。どの平均化で得た量か、どの境界条件・拘束の下で定義された量かを、式の形で記録して初めて、マルチスケールが科学として成立する。
17.2 計算結果の信頼性の作り方
計算の信頼性は「再現できる」だけでは十分ではない。研究としては、(i) モデル仮定、(ii) 数値誤差、(iii) 統計誤差、(iv) 実装差(コード・擬ポテンシャル等)を分解し、観測量(結論)へどう影響するか説明できる必要がある。ここでは“信頼性の作り方”を、最低限のチェック項目ではなく、一本の手順(体系)として整理する。
(1) 収束試験:離散化・基底・サンプリング どの手法でも「解像度を上げたら答えが安定するか」を確認する。
- DFT:カットオフ
、 点密度、スメアリング幅、セルサイズ、擬ポテンシャル/PAWの選択 - MD:時間刻み
、カットオフ半径、系サイズ、熱浴・圧力制御の時定数、平衡化時間と観測時間 - MC:熱化ステップ数、サンプル数、自己相関時間、レプリカ数(拡張アンサンブル)
- PF/FEM:メッシュ幅
、時間刻み 、要素次数 、線形/非線形ソルバ許容誤差、前処理
重要なのは、誤差が単純にゼロへ行くわけではなく、目的量によって収束が違う点である。例えば最大応力のような局所量は収束が遅く、全エネルギーや平均量は収束が早いことが多い。したがって「何を結論として使う量か」を先に決め、その量の収束を優先して見る必要がある。
(2) 相互検証:別経路で同じ量を確かめる 同じ物理量を、異なる定義や異なる手法で評価して一致を確認する。
例:熱伝導
- 平衡MD(Green–Kubo)
- 非平衡MD(温度勾配を与える) の2系統で照合し、サイズ依存や境界条件依存を切り分ける。
例:PFの界面エネルギー PFの勾配係数
(3) 感度解析:仮定・パラメータが結論を支配していないか 結論が特定のパラメータに過度に依存していないかを見る。一次感度は
で定義できるが、実際には無次元化して
のように相対感度で比較すると、どの仮定が支配的かが見えやすい。連成問題では、ソルバ許容誤差や時間刻みが物性値を動かしてしまう場合もあるため、物理パラメータだけでなく数値パラメータも感度対象に含めるのが実務的である。
(4) 不確かさ(統計誤差・モデル誤差)を“数で”出す MC/MDでは統計誤差が本質的であり、自己相関を考慮した有効サンプル数
と見積もられ、誤差は概ね
モデル誤差(力場・汎関数・粗視化の誤差)は、単一指標で出しにくいが、少なくとも
- 参照データ(DFT/実験)との比較
- 外挿検知(MLポテンシャル)
- 物理制約(保存則・対称性)の破れの監視 をセットで行うと、説明可能性が上がる。
信頼性とは、単に「正しい」ではなく「どこまで正しいと言えるか」を定量化し、読者が判断できる材料を揃えることだと捉えるのがよい。
17.3 データ駆動との融合
データ駆動は材料計算を加速するが、万能ではない。むしろ初学者にとって重要なのは、「どこをデータで置き換え、どこを物理で固定するか」を設計できることにある。ここでは、データ駆動を“追加の便利機能”ではなく、“不確かさを伴う近似器”として取り扱う。
(1) MLポテンシャルは「近似した力学系」である MLポテンシャルは、原子配置
として得る。したがって、エネルギーが滑らかでない・学習誤差が大きいと、力の誤差は増幅されやすい。ここから、学習時に「エネルギーだけでなく力も教師信号に入れる」「保存則(エネルギー保存)を監視する」などの設計が必要になる理由が理解できる。
(2) 外挿を検出しないと科学にならない 学習モデルは訓練分布の外側で誤る。外挿検出の代表はアンサンブルで、複数モデルのばらつき
を不確かさ指標として用いる。これをMDの途中で評価し、指標が大きい構造をDFTで追加計算して学習データへ戻す(アクティブラーニング)と、外挿を系統的に潰せる。ここで重要なのは、誤差が小さく見える範囲でも、観測量(例えば拡散係数)が大きく変わる場合があるため、“目的量に対する不確かさ”を意識することである。
(3) 逆問題は「同定できないことがある」問題である 材料定数同定やモデル同定は、観測
- 観測量を増やす(異なる条件、異なる測定)
- 正則化を入れる(事前知識で解を絞る)
- 感度の高い実験条件を選ぶ(実験計画) が必要になる。データ駆動は、ここを“自動で解く”のではなく、“何が同定できていないか”を可視化する道具として使うのが安全で強い。
(4) 物理制約とデータの役割分担 データ駆動を入れるときの基本戦略は、
- 物理で守る:保存則、対称性、境界条件、単位系、極限挙動(高温・低温など)
- データで補う:未知の構成式、経験則、複雑な相互作用の近似 という分担である。これにより、外挿や偽相関の危険を抑えつつ、探索空間を拡張できる。
17.4 計算機環境の発展
材料計算は、アルゴリズムと計算機の両方の制約の上に成立する。今後はGPUや異種混在計算が増え、計算量(flops)よりもメモリ帯域・通信が支配的になる傾向が強い。初学者が押さえるべきポイントを、手法ごとに整理する。
(1) DFT:対角化と通信が重い 平面波DFTでは、Kohn–Sham方程式の解法がボトルネックであり、多くの時間が固有値問題(対角化)に費やされる。系サイズが増えると、計算量は概ね急増し、通信も増える。したがって、単純に“コア数を増やす”だけで速くならない領域がある。並列化は
点並列 - バンド並列
- 平面波(FFT)並列 などの分割を組み合わせて設計される。
(2) MD:粒子数にほぼ比例、ただし長距離相互作用が効く 短距離相互作用なら計算量は粒子数
(3) FEM/PF:疎行列ソルバと前処理が支配 大規模FEMは疎行列の反復解法が中心で、前処理の出来が速度を決める。GPUは疎行列演算が得意だが、前処理が難しいことも多い。したがって「どの前処理が使えるか」を前提に、離散化と未知量の並べ方(ブロック構造)を決めるのが実務的である。
(4) 再現性と高速化の緊張関係 並列計算では削減演算の順序が変わり、丸め誤差が非決定性を生む場合がある。これはバグではなく数値の性質であり、対策は
- 収束判定を残差だけでなく物理量でも行う
- 乱数種、並列数、ライブラリ版数を記録する
- 重要結果は条件を変えて再現する である。GPU時代ほど、この“再現性をどう作るか”が重要になる。
17.5 研究倫理と公開可能性
計算研究の倫理は、捏造や改ざんだけではなく、「他者が追試できる状態で提示しているか」に強く関わる。材料計算は入力パラメータが多く、ちょっとした違いで結果が動くため、公開可能性(reproducibility)を意識した記録が不可欠である。
最低限記録すべき情報は、手法により異なるが、共通して以下が重要である。
- 入力ファイル(モデル・境界条件・初期条件)
- ソフトウェア名とバージョン、コンパイル条件
- 擬ポテンシャル/PAWデータ、汎関数名
- 乱数種(MC、MDの初期速度など)
- 収束条件(残差、許容誤差、反復上限)
- 解析スクリプト(後処理も含む)
- 単位系と物理量の定義(特に“エネルギー密度”など)
ここで重要なのは、結果の図だけでは不十分だという点である。例えば「拡散係数」と言っても、Einstein関係を使ったのか、速度自己相関を使ったのか、有限サイズ補正を入れたのかで値が変わり得る。つまり“量の定義”が公開情報に含まれて初めて科学になる。
さらに、公開の粒度も設計が必要である。全てを出すのが難しい場合でも、
- 最小再現パッケージ(入力+解析手順+要点の出力)
- 代表ケースの完全再現(全条件のうち一部を完全公開) のどちらかを提供できると、コミュニティでの信頼が大きく上がる。教育の段階でこれを習慣化すると、研究室運営や共同研究でもメリットが大きい。
17.6 学習順序の提案
初学者が材料計算を自力で組み立てられるようになるには、「物理」と「数値」を別々に覚えるのではなく、同じ方程式が“どの形で解かれているか”をセットで理解するのが近道である。ここまでの内容を踏まえ、学習順序を具体的に提案する。
(1) まず共通基盤:モデル化と離散化(第1・2章) 保存則
と、自由エネルギーの変分
を“読める”ようになると、DFT/MD/PF/FEMの文章が急に繋がり始める。ここで無次元化・スケール見積もり・安定性条件まで含めて理解すると、計算条件が自分で設計できるようになる。
(2) 次にミクロ:DFT → MD/MC DFTは基礎に見えるが、数値面では巨大な線形代数・収束判定の集合体である。したがって、まずは
- どの物理量がエネルギー微分で定義されるか
- 収束と誤差がどこから来るか を押さえ、次にMDやMCへ進むと、統計誤差や時間相関の意味が理解しやすい。
(3) その後メゾ・マクロ:PF → FEM → 連成 PFは連続体場としての組織進化を扱い、FEMは複雑形状・境界条件・連成を扱う。PFとFEMを学ぶ順としては、
- PFで変分原理から方程式が出ることを理解する
- FEMで弱形式と離散化を理解する
- PF-FEM連成や熱・応力・拡散連成へ進む が自然である。ここまで来ると、材料プロセス(熱処理・凝固・拡散・変形)の“計算による説明”が可能になる。
(4) 最後に統合:マルチスケールとデータ駆動 統合のゴールは「全部の手法を使う」ではなく、「必要な情報だけを正しく変換して渡す」ことである。DFTで求めた定数をFEMへ渡す、MDで得た輸送係数をPFへ渡す、といった接続点を、量の定義と誤差の見積もり込みで説明できることが、研究としての完成度を決める。
おわりに
本書では、材料計算学における主要な数値手法――第一原理計算(DFT)、分子動力学(MD)、モンテカルロ法(MC)、フェーズフィールド法(PF)、有限要素法(FEM)――を、単なる計算法の羅列としてではなく、「どのスケールを、どの未知量で記述するのか」という視点から体系的に整理してきた。原子・電子スケールから連続体スケールへと記述が移り変わる中で、計算手法の違いは精度や計算コストの差ではなく、物理量の定義と平均化の仕方の違いとして理解できることが、本書を通じて明確になったはずである。
また、異なる手法を接続する際に本質的に重要なのは、数値そのものを無理に受け渡すことではなく、「何を平均し、何を捨て、どの量を有効変数として残すのか」という定義の選択であることを繰り返し強調してきた。この視点を持つことで、マルチスケール計算は特別な技術ではなく、記述の階層を意識的に行き来するための論理的操作として捉えられるようになる。
計算結果の信頼性についても、本書では単一の基準で語ることを避け、収束性、相互検証、感度解析、不確かさの評価という複数の要素に分解して考える立場を採った。信頼性とは結果に付随する性質ではなく、計算の設計・実行・解釈の過程を通じて構築されるものである。この分解的な理解は、新しい手法や未経験の計算に直面した際にも、自ら判断軸を立てるための基盤となる。
さらに、近年重要性を増しているデータ駆動型手法については、万能な代替手段としてではなく、物理制約、外挿領域の検出、逆問題に内在する非同定性といった限界を意識したうえで使うべき道具として位置づけた。計算科学とデータ科学は対立するものではなく、記述能力と推論能力を補完し合う関係にあることを理解することが、今後の材料研究において不可欠である。
最終的に本書が目指したのは、特定の手法を使えるようになることではなく、学習の順序そのものを自分で設計できる研究者・技術者を育てることである。いま自分はどのスケールを理解しており、どの未知量が曖昧なのか、次に学ぶべき理論や手法は何か――それらを自ら判断できるようになったとき、材料計算学は単なる計算技術ではなく、思考の道具として機能し始める。本書が、そのための出発点となることを期待したい。
参考ドキュメント
- 東京大学物性研究所 計算物質科学研究センター, 第一原理計算プログラムとOpenMXについて(PDF): https://ccms.issp.u-tokyo.ac.jp/wp-content/uploads/2019/12/openmx_text.pdf
- HPCI, PHASE/0(アプリケーション紹介): https://www.hpci-office.jp/for_users/appli_software/appli_phase0
- 小山敏幸, フェーズフィールド法による組織形成シミュレーション(日本語解説, J-STAGE): https://www.jstage.jst.go.jp/article/jinstmet/73/12/73_12_955/_article/-char/ja/
参考文献(その他)
- P. Hohenberg and W. Kohn, Inhomogeneous Electron Gas, Physical Review 136, B864 (1964).
- W. Kohn and L. J. Sham, Self-Consistent Equations Including Exchange and Correlation Effects, Physical Review 140, A1133 (1965).
- J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized Gradient Approximation Made Simple, Physical Review Letters 77, 3865 (1996).
- J. Sun, A. Ruzsinszky, and J. P. Perdew, Strongly Constrained and Appropriately Normed Semilocal Density Functional, Physical Review Letters 115, 036402 (2015).
- J. W. Furness et al., Accurate and Numerically Efficient r2SCAN Meta-GGA, Physical Review Materials 4, 083801 (2020).
- N. Metropolis et al., Equation of State Calculations by Fast Computing Machines, J. Chem. Phys. 21, 1087 (1953).
- R. H. Swendsen and J.-S. Wang, Nonuniversal critical dynamics in Monte Carlo simulations, Phys. Rev. Lett. 58, 86 (1987).
- U. Wolff, Collective Monte Carlo Updating for Spin Systems, Phys. Rev. Lett. 62, 361 (1989).
- A. B. Bortz, M. H. Kalos, and L. J. Lebowitz, A New Algorithm for Monte Carlo Simulation of Ising Spin Systems, J. Comput. Phys. 17, 10 (1975).
- L. Verlet, Computer Experiments on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules, Phys. Rev. 159, 98 (1967).
- S. Nosé, A unified formulation of the constant temperature molecular dynamics methods, J. Chem. Phys. 81, 511 (1984).
- W. G. Hoover, Canonical dynamics: Equilibrium phase-space distributions, Phys. Rev. A 31, 1695 (1985).
- M. Parrinello and A. Rahman, Polymorphic transitions in single crystals: A new molecular dynamics method, J. Appl. Phys. 52, 7182 (1981).
- J. W. Cahn and J. E. Hilliard, Free Energy of a Nonuniform System. I. Interfacial Free Energy, J. Chem. Phys. 28, 258 (1958).
- S. M. Allen and J. W. Cahn, A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening, Acta Metall. 27, 1085 (1979).
- L.-Q. Chen, Phase-Field Models for Microstructure Evolution, Annu. Rev. Mater. Res. 32, 113 (2002).
- A. Batzner et al., E(3)-equivariant Graph Neural Networks for Data-Efficient and Accurate Interatomic Potentials (NequIP), Nature Communications (2022).
- MOOSE Phase Field module documentation: https://mooseframework.inl.gov/modules/phase_field/index.html