Skip to content

材料計算学

本書は、材料物理・材料工学で用いられる計算手法を、数式に基づいて一貫して学ぶための教科書として設計するものである。連続体から原子スケールまでの記述を往復しながら、計算結果の信頼性を担保するための理論的根拠と手法選択の判断軸を身につける構成である。

第1章 数値モデル化と支配方程式

本章では、材料現象を計算可能な形にするために、どの量を変数として選び、どの方程式を立てるのかを数式の形で整理する。連続体としての記述と原子・粒子としての記述を往復しながら、モデル化の前提、支配方程式の導出、条件設定の意味を一つずつ確認する。

1.1 連続体モデルと粒子モデル

材料計算では、状態を場として表す連続体モデルと、原子や粒子の集まりとして表す粒子モデルの二つが基本である。連続体モデルでは、位置 r と時刻 t の関数として、変位 u(r,t)、温度 T(r,t)、濃度 c(r,t)、磁化 M(r,t) のような場を未知量にする。一方、粒子モデルでは、原子 i の座標 ri(t) と速度 vi(t) を未知量にし、運動方程式で時間発展を追う。

両者の違いは、何を「最小単位」とみなすかにある。連続体モデルは、材料内部の微細構造を直接には追わず、ある体積要素の平均量で状態を表す考え方である。これは、観測したい現象の空間スケールが原子間距離より十分大きく、平均化した量の変化で現象を説明できる場合に有効である。例えば、厚さ数mmの試験片の温度分布を議論するなら、原子の個数を数えるよりも、温度場 T(r,t) を扱う方が自然である。

この平均化の考え方を数式で表すと、連続体の場は粒子の情報から粗視化して得られる量として定義できる。例えば、密度場 ρ(r,t) は、各原子の質量 mi を滑らかなカーネル関数 W(rri) で平均化して

ρ(r,t)=imiW(rri(t))

のように書ける。ここで W は、W(r)dr=1 を満たす平滑化関数であり、粗視化の長さスケールを含む。粗視化長さを大きくするとノイズは減るが細部は失われ、逆に小さくすると原子性の揺らぎが強く残る。したがって、連続体と粒子の間には、観測スケールに応じた「どの程度まで平均するか」という自由度がある。

粒子モデルは、局所構造や熱揺らぎ、欠陥の原子レベル機構を直接追える点が強い。例えば拡散は、原子のランダムなジャンプの積み重ねであり、粒子モデルではその起点をそのまま記述できる。一方で粒子数 N が大きいほど計算量が増え、時間刻みも原子振動周期に合わせて小さくせざるを得ないため、長時間・大域スケールの現象には不利になりやすい。

連続体モデルは、長時間・大域の挙動を扱いやすい点が強い。例えば拡散方程式や弾性方程式は、空間に分布する量の滑らかな変化として現象を表すため、メッシュ幅を現象の代表長さに合わせて選べる。ただし、材料定数(弾性定数、拡散係数、界面エネルギーなど)が入力として必要になり、それらがどの条件で有効な値かを理解しないと結果の解釈を誤りやすい。連続体モデルは「何を捨てたか」を明示しながら使う姿勢が重要である。

材料計算の多くは、この二つを使い分けるだけではなく、接続する。粒子モデルで得た局所物性を連続体モデルに渡し、連続体モデルで得た場を粒子モデルの外場として与える。以降の章で扱うDFT、MD、MC、フェーズフィールド、有限要素法は、この「どの自由度で書くか」という選択の上に配置される。

1.2 保存則と生成消滅項

材料現象の支配方程式は、保存則から作るのが最も堅牢である。保存則は、ある領域に含まれる量の増減が、境界を通る流れと内部での生成消滅で決まるという主張である。これをまず積分形で書くと、領域 Ω に対して

ddtΩq(r,t)dV=ΩJ(r,t)ndS+Ωs(r,t)dV

となる。ここで q は体積あたりの保存量密度、J は流束、s は体積あたりの生成消滅率、n は外向き法線である。左辺は領域内総量の時間変化、右辺第1項は境界から外へ出ていく正味の流れ、右辺第2項は領域内部での生成消滅である。

次に、発散定理

ΩJndS=ΩJdV

を用いると、積分形は

Ω(qt+Js)dV=0

となる。これは任意の領域 Ω で成り立つ必要があるので、被積分関数がゼロでなければならない。したがって局所形として

qt+J=s

が得られる。この導出は、材料計算で繰り返し現れる基本操作であり、方程式の形がなぜこの形になるのかを理解する最短経路である。

この一般形の強みは、物理を qJs に分解して考えられる点にある。未知なのは多くの場合 J の構成式である。例えば拡散で保存量を濃度 c とし、経験則としてフィックの法則

J=Dc

を採用すると、

ct=(Dc)+s

が得られる。D が一定なら c/t=D2c+s となり、これが拡散方程式である。ここで 2= はラプラシアンであり、濃度の凹凸を平らにする方向へ時間発展が進むことを表す。

同じ枠組みで熱伝導も書ける。保存量を単位体積あたりの内部エネルギー q=ρcpT とし、フーリエの法則 J=kT を用いれば

ρcpTt=(kT)+s

が得られる。s は発熱(ジュール熱や反応熱)などである。形は拡散方程式と同型であり、同じ数値手法が流用できる理由がここにある。

運動量の保存則では、q=ρv を取り、流束は応力テンソル σ を含む形になる。外力密度 f を生成項として

(ρv)t+(ρvvσ)=f

のように書ける。弾性体では σ をひずみから決める構成式(フック則)で閉じ、流体では粘性応力で閉じる。ここでも、保存則は形を固定し、材料の性質は構成式として入ってくる。

生成消滅項 s は、反応、相変態、放射線損傷、析出などで現れる。例えば一次反応で濃度が減るなら s=kc のように置ける。このとき拡散方程式は反応拡散方程式となり、拡散による空間平滑化と反応による局所減衰が競合する。材料組織の形成やパターン形成は、まさにこの競合で生じる場合が多い。

以降の章で扱う多くの方程式は、見かけは異なっても、この保存則の枠組みか、次節の変分原理の枠組みに分類できる。どちらの出発点から来ているかを意識すると、式変形や境界条件の意味が理解しやすくなる。

1.3 変分原理と自由エネルギー

材料の平衡状態は、自由エネルギーが最小となる状態として決まる場合が多い。ここで言う自由エネルギーは、温度一定ならヘルムホルツ自由エネルギー、圧力一定ならギブズ自由エネルギーなどであり、条件に応じて適切な熱力学ポテンシャルを選ぶ必要がある。連続体の場 u(r) を用いると、自由エネルギーは汎関数(関数を入力として数を返す写像)として

F[u]=ΩL(u,u,r)dr

のように書ける。L は局所自由エネルギー密度であり、u やその勾配に依存する形を取る。

平衡条件は、u を微小に変えたときの F の一次変化がゼロになることで与えられる。uu+ϵη と変化させ、ϵ に関する一次の項だけを取り出すと

δF=ddϵF[u+ϵη]|ϵ=0=Ω(Luη+L(u)η)dr

となる。第二項に部分積分(多次元では発散定理に相当)を適用すると

ΩL(u)ηdr=Ω(L(u))ηdr+Ω(L(u)n)ηdS

が得られる。境界での変分の扱いが適切に定まっているとすると、領域内部での平衡条件は

Lu(L(u))=0

となる。これがオイラー・ラグランジュ方程式であり、変分原理から支配方程式が導かれる基本形である。

材料で特に重要なのは、勾配エネルギー項を含む自由エネルギーである。例えば

F[u]=Ω(f(u)+κ2|u|2)dr

を考える。ここで f(u) は局所的な自由エネルギー(多井戸形など)であり、κ>0 は勾配エネルギー係数である。このとき

Lu=f(u),L(u)=κu

なので、平衡条件は

f(u)κ2u=0

となる。すなわち、勾配エネルギー項はラプラシアン 2u を通じて現れ、場を滑らかにする方向へ働く。この構造はフェーズフィールドで界面幅が有限になる理由を与え、界面張力に相当する量が勾配項から生じることにつながる。

変分原理は平衡条件だけでなく、時間発展のモデル化にも使える。平衡へ向かって自由エネルギーが減少する緩和を表したいなら、非保存量 u の場合に

ut=LδFδu

と置くのが自然である。ここで L>0 は緩和係数であり、δF/δu は汎関数微分である。先ほどの F[u] なら

δFδu=f(u)κ2u

であり、結果として

ut=L(f(u)κ2u)

が得られる。これは非保存秩序変数の緩和方程式の基本形である。

一方、保存される量(濃度 c など)では、単に c/t=LδF/δc とすると総量が保存されない。保存を保ったまま自由エネルギーを減らすには、流束を

J=Mμ,μ=δFδc

と置き、保存則 c/t+J=0 と組み合わせて

ct=(MδFδc)

とする。ここで M>0 は移動度である。保存則と変分原理が一つの方程式の中で結びつく点が重要であり、材料の相分離や析出の記述の核になる。

この節で導いた「自由エネルギーから駆動力が出る」という構造は、DFTでは電子密度に対する変分、フェーズフィールドでは秩序変数に対する変分、有限要素法では弱形式(エネルギー原理)として再登場する。式の見た目が変わっても、根にある考え方は共通である。

1.4 初期条件と境界条件

支配方程式は、方程式だけでは解が決まらず、初期条件と境界条件が必要である。初期条件は t=0 の状態を与え、境界条件は領域境界 Ω で何が起きているかを指定する。材料計算では、境界が「現実の試験片の表面」「周期セルのつなぎ目」「計算領域の人工境界」のいずれであるかにより、適切な条件の形が変わる。

代表的な境界条件は次の三つである。ディリクレ条件は場の値を固定する条件であり、例えば温度固定なら

T(r,t)=T0(rΩ)

である。ノイマン条件は法線方向の勾配や流束を固定する条件であり、熱流束 qn を与えるなら

kTn=qn(rΩ)

である。ロビン条件は値と流束の線形結合を与え、対流境界

kTn=h(TT)

のように、外部との熱交換を表すのに使われる。

周期境界条件は、無限に繰り返す材料やバルクの一部を表したいときに用いる。例えば一次元なら

u(0,t)=u(L,t),ux(0,t)=ux(L,t)

のように、値と勾配の一致を課す。これは境界の影響を消すための道具であり、欠陥や界面がある場合は、周期像との相互作用が生じないようセルサイズを十分大きく取る必要がある。

吸収境界条件は、波動や弾性波、電磁波などで計算領域の端から反射が戻ってくるのを抑える目的で用いる。例えば波動方程式で、領域端で外向きに進む波を想定し

ut+cun=0

のような形を置くと、端での反射が減る場合がある。吸収の良さは周波数や入射角に依存するため、問題に応じて設計が必要である。

初期条件は、時間依存問題では必須である。拡散方程式なら c(r,0)=c0(r)、弾性の静的問題なら初期条件は不要だが境界条件が増える。ここで重要なのは、方程式の種類によって「何をどれだけ与えれば解が一意に定まるか」が異なる点である。例えば拡散方程式は一次の時間微分なので初期条件を1つ与えるが、波動方程式は二次の時間微分なので u(r,0)u/t(r,0) を両方与える必要がある。

境界条件は数値計算の実装とも深く結びつく。有限要素法の弱形式では、ディリクレ条件は本質境界条件として試験関数空間を制約する形で入る。一方、ノイマン条件は自然境界条件として境界積分項に現れ、弱形式に組み込まれる。これは、方程式の導出(部分積分で境界項が出る)と境界条件の種類が対応していることを意味する。

材料計算における条件設定は、物理の再現だけでなく、問題が数学的に解ける形になっているかにも関わる。例えば拡散方程式で全境界をノイマン条件(流束指定)にすると、総量保存のため定数分だけ解が不定になる場合がある。そのときは平均値を別途固定するなど、解の自由度と条件数の関係を意識して設計する必要がある。

1.5 無次元化とスケールの切り出し

無次元化は、式の中に隠れている重要な比を取り出し、何が支配的かを見える形にする操作である。さらに、数値計算では適切な格子幅 Δx と時間刻み Δt を選ぶための指針にもなる。やることは単純で、代表長さ L、代表時間 τ、代表量 u0 を決めて

r=Lr,t=τt,u=u0u

とおき、方程式を の付いた無次元変数で書き直す。

例として拡散方程式

ct=D2c

を無次元化する。=(1/L)/t=(1/τ)/t を用いると

u0τct=Du0L22c

であり、u0 は消えて

ct=DτL2Fo2c

となる。Fo=Dτ/L2 は無次元数であり、拡散が長さ L にわたって進むために必要な時間が τL2/D であることを示す。したがって、拡散の代表時間は

τdiffL2D

で見積もれる。

次に移流拡散方程式

ct+vc=D2c

を考える。速度の代表値を V として同様に無次元化すると

ct+VτLCovc=DτL2Fo2c

となる。ここで Co=Vτ/L は移流が長さ L を進む時間の尺度である。τ=L/V を選べば移流項の係数が1になり、拡散項の係数は

DVL=1Pe

となる。Pe=VL/D はペクレ数であり、移流と拡散の相対的重要度を表す。Pe1 なら移流が支配的、Pe1 なら拡散が支配的である。

このような無次元数は、数値条件にも直結する。例えば陽解法で拡散方程式を解く場合、空間を格子幅 Δx で離散化すると、安定性のために

Δt(Δx)22dD

のような上限が現れる(d は次元)。これは、拡散の代表時間が L2/D であることと整合しており、格子を細かくするほど Δt を二乗で小さくしないといけないことを意味する。移流が支配的なら、別の制約として

ΔtΔxV

が現れる。したがって、現象の支配項が何であるかを無次元化で把握すると、どの制約が支配的になるかも見通せる。

さらに、無次元化はモデルの単純化にも使える。ある無次元数が十分小さいなら、その項を無視した近似が成立しうる。例えば Pe1 なら移流を落として拡散だけのモデルに近づく。ただし、境界層や局所構造では小さいはずの項が重要になる場合があるため、無視の判断は空間領域全体で一律に行うのではなく、現象の局所性も考慮する必要がある。

材料計算では、同じ名前の方程式でも材料定数が桁違いになることがある。無次元化は、材料ごとに式の意味を揃える作法であり、異なる系の結果を比較するための共通のものさしになる。

1.6 連続体と離散体の接続

連続体モデルと粒子モデルは、どちらか一方を選べば終わる関係ではない。材料物性は階層的であり、原子スケールの相互作用がメソスケールの組織を作り、それがマクロな応答(強度、導電率、磁性など)を決める。したがって、計算でも階層をまたぐ情報の受け渡しが必要になる。

第一の方向は、離散体から連続体への粗視化である。代表例は、原子計算から材料定数を抽出して連続体モデルに渡す操作である。例えば弾性定数 Cijkl は、ひずみ ε を与えたときのエネルギー変化から求められる。体積 V を持つ結晶で、微小ひずみに対してエネルギーが

E(ε)=E0+V2i,j,k,lCijklεijεkl+O(ε3)

と展開できることを用い、いくつかの独立なひずみ状態を計算して二次係数をフィットする。ここで重要なのは、微小ひずみ領域で二次近似が妥当であることを確認し、数値ノイズや応力の収束誤差が係数に混ざらないようにする点である。

拡散係数 D も同様に、粒子モデル(分子動力学)から連続体の係数として求められる。Einsteinの関係により、d 次元で

D=limt12dt|r(t)r(0)|2

で定義できる。ここで は粒子や初期条件に関する平均である。これにより、拡散方程式の D が「原子のランダム運動の長時間極限の傾き」として具体化される。どの温度、どの欠陥濃度で得た D なのかを明記しないと、連続体モデルに入れたときの意味が曖昧になるため、条件の対応付けが必須である。

第二の方向は、連続体から離散体への情報供給である。例えば、連続体側で得られた応力場や温度場を、原子計算の境界条件や外場として与え、局所構造変化を調べることがある。温度勾配下の拡散、応力誘起の相変態、電場下のイオン移動などでは、マクロな場がミクロな遷移率を変えるため、この方向の接続が重要になる。

連続体と離散体の接続では、同じ名前の量でも定義が異なる場合がある点に注意が必要である。例えば応力は、連続体では面に作用する力の密度として定義されるが、原子系ではビリアル応力として

σ=1V(imivivi+i<jrijfij)

のように表される。ここで rij=rirjfij は相互作用力である。温度が高いと速度項の寄与も無視できず、どの体積 V を採るかも結果に影響する。したがって、接続では量の定義と平均化の手続きをセットで扱う必要がある。

さらに、界面や粒界のように連続体の仮定が破れやすい領域では、接続の設計が計算結果の解釈を左右する。例えばフェーズフィールドで界面エネルギー γ と界面幅 w をパラメータとして使うなら、これらを原子計算や実験から同定する必要がある。勾配エネルギー係数 κ と二重井戸形 f の形が与えられたとき、γw がどう決まるかは変分原理の枠組みで追跡できるため、単なる調整パラメータとして扱わず、どの物理を表しているかを説明できる形にしておくことが重要である。

連続体と離散体を接続する作業は、材料計算の中心である。支配方程式の形を理解し、保存則と自由エネルギーという二つの出発点を押さえたうえで、どの物理量をどの階層で定義し、どう平均し、どう渡すのかを一貫させることが、計算結果を材料科学の議論へつなげる基本動作である。

第2章 離散化と偏微分方程式の数値解法

本章では、偏微分方程式をコンピュータで解くために、連続な式を有限個の未知数へ落とし込む手順を、式変形から順に整理する。空間離散と時間離散の設計が、解の安定性・精度・計算量を同時に決めることを、最小の例から一般形へ広げて理解する。

2.1 空間離散の設計

偏微分方程式は、未知量が連続な位置 r の関数であるため、そのままでは有限時間で計算できない。そこで、空間を格子点、セル、要素などに分割し、未知量を有限個の自由度(数値)として近似する。この操作が空間離散であり、差分法・有限体積法・有限要素法・スペクトル法は、どの量を未知数にし、どの形で方程式を満たすかが異なる。

まず、一次元で直感を作る。例えば u(x) の導関数を格子点 xi=iΔx 上で近似する。テイラー展開を用いると

u(xi+1)=u(xi)+u(xi)Δx+12u(xi)Δx2+O(Δx3),u(xi1)=u(xi)u(xi)Δx+12u(xi)Δx2+O(Δx3)

である。これらを引き算すると

u(xi+1)u(xi1)=2u(xi)Δx+O(Δx3)

となり、中心差分

u(xi)ui+1ui12Δx

が得られる。誤差が O(Δx2) であることは、格子幅を半分にすると誤差が約1/4になることを意味する。二階微分も同様に足し算から

u(xi)ui+12ui+ui1Δx2

が得られ、拡散方程式の空間離散の基礎になる。

差分法は、このような導関数近似を多次元に拡張して方程式を点ごとに満たす。例えば二次元のラプラシアンは

2uui+1,j2ui,j+ui1,jΔx2+ui,j+12ui,j+ui,j1Δy2

となる。格子が直交で、境界形状が単純な場合に実装が簡潔で高速に動きやすい。一方で、材料の実形状が複雑である場合や、界面が曲線・曲面である場合には、格子に形状を合わせることが難しくなる。

有限体積法は、保存則の積分形から出発し、セル平均を未知数として扱う。保存則

qt+J=s

をセル Vk で積分すると

ddtVkqdV+VkJdV=VksdV

である。発散定理により

VkJdV=VkJndS

となるので、セル内総量の変化はセル境界を通る流束の和で決まる形になる。これを離散化すると、一次元では

ddt(qkΔx)+(Jk+1/2Jk1/2)=skΔx

のようになり、数値的に「出入りの差で増減する」構造がそのまま残る。材料の拡散・熱伝導・流体輸送などで、保存が重要な量を扱うときに強い。ただし、高次精度を狙う場合や界面近傍での流束評価(再構成)を丁寧に設計しないと、数値拡散や振動が出やすい。

有限要素法は、微分方程式を弱形式に変換してから近似する。例としてポアソン方程式

(ku)=f(rΩ)

を考える。試験関数 v を掛けて領域積分し、

Ωv(ku)dV=ΩvfdV

とする。左辺に部分積分(発散定理)を適用すると

ΩkuvdVΩvkundS=ΩvfdV

が得られる。境界項はノイマン境界条件として自然に入り、ディリクレ境界条件は u の関数空間を制約する形で入る。次に、近似として

uh(r)=a=1NUaϕa(r)

と展開する。ϕa は形状関数であり、Ua が未知係数である。これを弱形式に代入し、試験関数として v=ϕb を選ぶと

aUaΩkϕaϕbdV=ΩfϕbdV+Ωϕbt¯dS

となる。ここで t¯=kun は境界流束である。この式は行列形式 KU=F になり、複雑形状や異方性 k(r) を自然に扱える。材料計算で頻出する複雑な境界、局所的な物性変化、連成問題に対して強い一方、要素形状の品質、数値積分の精度、境界条件の実装が結果に直結する。

スペクトル法は、未知関数をフーリエ級数や直交多項式で展開する。周期境界の一次元で

u(x)=m=MMu^meikmx,km=2πmL

と書けるとき、微分は係数空間で

dudx=m(ikm)u^meikmx

となり、高い精度で微分が計算できる。解が十分滑らかである場合、誤差が非常に速く減る性質を持つが、界面の不連続や局所的な急峻勾配があるとギブス現象が出やすい。材料では界面や欠陥が本質の場合が多いため、スペクトル法は適用条件を見極めて使う必要がある。

材料現象では、界面近傍の急峻な勾配、結晶方位に依存する異方性、複雑な形状、複数場(温度・応力・濃度など)の連成が同時に現れる。そのため、どの離散化が何を得意とし、どの条件で不利になるのかを、式の構造と幾何の扱いから説明できることが重要である。

手法強み注意点
差分法構造が単純な領域で高速化しやすく、導関数近似の意味が追いやすい形状が複雑だと格子生成や境界近傍の精度が難しくなりやすい
有限体積法保存則が離散形でも保たれ、輸送問題で物理の整合性を保ちやすい流束再構成や界面条件の設計が結果を左右しやすい
有限要素法複雑形状・異方性・連成を弱形式で統一でき、境界条件も整理しやすい要素品質と数値積分、拘束条件の実装が精度と安定性に直結する
スペクトル法解が滑らかな場合に非常に高精度になりやすい不連続・界面を含む問題では振動が出やすく工夫が要る

2.2 時間積分と安定性

空間離散により、偏微分方程式は有限次元の常微分方程式系に落ちる。未知ベクトルを u(t) とすると、多くの場合

dudt=F(u,t)

の形になる。時間積分は、この系を離散時間 tn=nΔt 上で更新式に変換する操作である。材料計算では、拡散型(緩やかに平滑化)、波動型(伝播と反射)、反応型(局所増減)などが混ざるため、安定性の制約が支配的になることが多い。

最も基本の陽解法は前進オイラー法で、

un+1=un+ΔtF(un,tn)

である。実装が簡潔で並列化もしやすい一方、安定性が時間刻みを強く制約する。安定性を理解するために、線形の試験方程式

dydt=λy

を考える。前進オイラー法では

yn+1=(1+λΔt)yn

であり、増幅率 g=1+λΔt|g|1 を満たすときに安定である。拡散問題では λ は負の実数になりやすいので、例えば λ=αα>0)なら

|1αΔt|10αΔt2

となり、Δt2/α の上限が現れる。ここで α は空間離散により Δx2 に比例することが多く、結果として Δt(Δx)2 に比例して小さくなる。これが拡散型方程式で陽解法が厳しい理由である。

具体例として一次元拡散方程式

ut=D2ux2

を中心差分で離散化すると

duidt=Dui+12ui+ui1Δx2

となる。前進オイラー法で時間離散すると

uin+1=uin+r(ui+1n2uin+ui1n),r=DΔtΔx2

である。フォン・ノイマン解析により、安定には r1/2(一次元の場合)が必要であり、

ΔtΔx22D

が得られる。格子を細かくして界面を解像しようとすると、時間刻みが二乗で小さくなるため、計算回数が急増する。

波動や移流では別の制約が出る。例えば移流方程式

ut+cux=0

は情報が速度 c で伝播する。陽的スキームでは、時間刻みが「1ステップで情報が何格子進むか」により制約され、CFL条件として

ΔtCΔx|c|

が現れる(C はスキームに依存する定数)。これは、数値が物理的に伝わる速さより速く更新されると不安定になるという意味であり、拡散の Δx2 制約とは性質が異なる。

これらの制約を緩めるのが陰解法である。後退オイラー法は

un+1=un+ΔtF(un+1,tn+1)

であり、右辺が未知の un+1 を含むため、各ステップで方程式を解く必要がある。しかし試験方程式 y=λy に対しては

yn+1=11λΔtyn

となり、λ が負の実数であれば任意の Δt|g|1 となる。拡散型方程式では、陰解法により安定性制約が大幅に緩和され、実用的な時間刻みが選べる。代わりに、各ステップで疎行列方程式を解く計算が必要になり、次節の線形ソルバが重要になる。

精度の観点では、時間積分の次数も意識する必要がある。前進オイラー法と後退オイラー法は一次精度であり、時間刻みを半分にすると誤差はおよそ半分になる。二次精度の代表としてクランク・ニコルソン法があり、線形問題では

un+1=un+Δt2[F(un,tn)+F(un+1,tn+1)]

である。拡散型では安定である一方、強い非線形や不連続があると振動的な誤差が残る場合があるため、安定性と減衰性の両方を見て選ぶ必要がある。

材料計算では、拡散のように硬い項と、比較的緩やかな項が混ざることが多い。このとき、硬い項だけ陰に、その他を陽にするIMEX(Implicit-Explicit)法や、演算子分割が有効になる。例えば

dudt=Au+N(u)

A が硬い線形項なら

un+1=(IΔtA)1[un+ΔtN(un)]

のようにして、安定性と計算量を両立させる設計ができる。時間積分は、単に式を進める操作ではなく、安定性と計算コストの配分を決める設計問題である。

2.3 疎行列と線形ソルバ

空間離散と時間離散を行うと、多くの場合で線形方程式

Ax=b

を解くことになる。有限要素法の剛性行列、陰解法の時間発展行列、ポアソン方程式の離散化行列などが代表である。ここで重要なのは、A が疎行列であることである。疎行列とは、成分の大半がゼロで、非ゼロが局所近傍に限られる行列であり、格子やメッシュの局所結合を反映している。

直接法(ガウス消去、LU分解、Cholesky分解)は頑健であり、中規模問題では扱いやすい。しかし三次元で自由度 N が増えると、分解に伴うフィルイン(ゼロだった成分が非ゼロになる)が増え、メモリと計算量が急増する。特に大規模材料計算では、反復法が中心となる。

反復法の理解のために、線形方程式を残差で見る。近似解 xk に対し

rk=bAxk

を残差と呼び、rk を小さくすることが目標である。Krylov部分空間法は、A の作用を繰り返し使って

Km(A,r0)=span{r0,Ar0,A2r0,,Am1r0}

という部分空間上で最良近似を探す方法である。行列の全要素を扱うのではなく、疎行列の行列ベクトル積だけを中心に計算が進むため、大規模問題で有利になる。

A が対称正定値(SPD)である場合は共役勾配法(CG)が使える。ポアソン方程式や弾性の線形問題ではSPDになることが多い。CGはエネルギー最小化として解釈でき、誤差の減衰は条件数 κ(A) に依存する。大まかには、条件数が大きいほど収束が遅くなる。条件数はメッシュを細かくすると増えやすく、材料定数の強い不均一でも悪化する。

非対称の場合(移流を含む、非自己共役な離散化、あるいは連成によりブロック構造になる場合)はGMRESなどが用いられる。GMRESは残差ノルムを最小化するが、反復回数とともに記憶量が増えるため、再直交化やリスタートなどの工夫を用いる。材料計算では、物性の強い異方性や不連続があると、行列が悪条件になりやすく、ソルバの選択が計算時間を支配する。

ここで決定的に重要なのが前処理である。前処理は、元の問題 Ax=b を「解きやすい形」に変換する操作であり、左前処理なら

M1Ax=M1b

を解く。理想は MA で、M1A の条件数が小さくなれば収束が速くなる。単純な例として、Jacobi前処理は M=diag(A) を用い、スケーリングにより収束を改善する。より強い前処理として不完全LU分解(ILU)があり、フィルインを制限した近似分解である。ブロック構造の連成問題では、物理に合わせたブロック前処理(Schur補完に基づく分割)を設計することで、桁違いに収束が改善することがある。

線形ソルバの停止判定は、単に反復回数を固定すればよいわけではない。残差が十分小さいかを基準にするが、残差の大きさは右辺 b の大きさにも依存するため、相対残差

rkb

を用いることが多い。さらに、物理量として意味のある誤差(例えばエネルギー誤差や保存量誤差)を併せて監視しないと、数値的には収束していても物理的には不適切な解になる場合がある。線形ソルバは、PDEの解法の中心部であり、離散化と同じくらい「設計」の対象である。

2.4 非線形ソルバと変分問題

材料の相変態、塑性、接触、化学反応、非線形拡散などでは、離散化後に非線形方程式

R(u)=0

を解く必要がある。ここで R は残差ベクトルであり、例えば陰解法で時間発展を離散化した場合には

R(un+1)=un+1unΔtF(un+1)=0

の形になる。非線形性は、構成式が非線形であること、自由エネルギーが多井戸形であること、あるいは制約条件が入ることから生じる。

基本となるのはニュートン法である。現在の近似 uk の周りでテイラー展開し、

R(uk+δu)R(uk)+J(uk)δu

とする。ここで J=R/u はヤコビアンである。R(uk+1)=0 を目指すので、線形化により

J(uk)δu=R(uk),uk+1=uk+δu

を繰り返す。解の近くでは二次収束するため非常に速いが、初期値が悪いと発散することがある。また、各反復で線形方程式を解く必要があり、その線形ソルバの性能が全体を支配する。

ニュートン法を安定に使うために、線形探索や減速を入れる。例えば更新を

uk+1=uk+αkδu,0<αk1

とし、残差が確実に減るように αk を選ぶ。材料の非線形問題は、強い非線形性や多解性を持つことがあり、数値的には「どの解に収束するか」も問題になる。したがって、初期値の物理的意味、境界条件、負荷経路の与え方が重要である。

ニュートン法が重い場合、固定点反復(ピカール反復)を使うこともある。例えば

u=G(u)

の形に書けるとき

uk+1=G(uk)

で更新する。収束は遅いが、ヤコビアンを作らずに済み、実装が簡単になることがある。塑性の内部変数更新や、弱い非線形の拡散係数更新などで用いられる。ただし、収束条件は写像の縮小性に依存するため、緩和係数を入れて

uk+1=(1ω)uk+ωG(uk)

のようにする場合も多い。

変分問題として書ける場合は、非線形ソルバの見通しがよくなる。例えばエネルギー汎関数 Π(u) を最小化する問題

minuΠ(u)

は、最適性条件 Π(u)=0 が非線形方程式になっている。弾性の平衡、フェーズフィールドの平衡、あるいは時間離散後のエネルギー安定スキームなどがこの枠組みに入る。変分の利点は、エネルギーが単調に減少する、制約が自然に組み込める、といった物理的な性質を離散化の段階で保ちやすい点にある。

例えば、時間発展がエネルギーを減らすべき現象では、離散化後も

F(un+1)F(un)

が成立するように時間積分を設計すると、数値的不安定や非物理的振動を抑えやすい。この考え方は、フェーズフィールドの安定化スキームや、非線形拡散の陰解法設計で特に重要である。非線形ソルバは、単に収束させるだけでなく、物理として守るべき性質を壊さない形で解くことが目標になる。

2.5 マルチグリッドと前処理

大規模PDEの計算時間の多くは線形ソルバが占めるため、収束を速くする工夫が中心課題になる。前処理はその中核であり、マルチグリッドは前処理の最も強力な考え方の一つである。発想は、誤差には波長の短い成分と長い成分があり、反復法は短波長誤差を減らすのが得意だが、長波長誤差を減らすのが苦手である、という観察に基づく。

例えば一次元ポアソン方程式を考えると、Jacobi法やGauss-Seidel法は局所的な平均化なので、格子1つ2つで振動する誤差はすぐ減る。一方、領域全体でゆっくり変化する誤差は、局所平均化ではなかなか消えない。ここで、粗い格子へ移ると、その「ゆっくり変化する誤差」が粗格子では短波長に見えるため、少ない反復で減らせる。これがマルチグリッドの基本原理である。

マルチグリッドの構成は、滑らか化(smoothing)、制限(restriction)、粗格子補正(coarse-grid correction)、補間(prolongation)からなる。線形方程式 Ax=b に対し、現在の近似 x の残差 r=bAx を計算し、これを粗格子に制限して rH を作る。粗格子行列 AH 上で

AHeH=rH

を解き(あるいは近似的に解き)、誤差補正 eH を細格子へ補間して xx+PeH のように更新する。これにより、細格子では減りにくい長波長誤差が効率よく補正される。

幾何学的マルチグリッドは、格子やメッシュが階層構造(2倍粗くするなど)を持つ場合に適用しやすい。材料計算では複雑形状の有限要素メッシュが多く、規則的な粗格子が作りにくい場合がある。そのときは代数的マルチグリッド(AMG)が有効であり、行列 A の情報だけから粗格子(粗変数)を構成する。AMGはブラックボックス的に使える場合が多いが、物性不連続が極端な場合や連成ブロック系では、物理に合わせた設定が必要になることもある。

マルチグリッドは、それ自体をソルバとして使うだけでなく、Krylov法の前処理として使うことが多い。例えばCGやGMRESの前処理にAMGを入れると、反復回数が大きく減る場合がある。ここでの狙いは、M1A のスペクトルを改善し、条件数や固有値分布を良くして収束を速めることである。

材料計算では、メッシュの局所細分化、異方性、物性の桁違いの変化がしばしば現れ、これらがソルバを難しくする。マルチグリッドと前処理は、単なる計算高速化ではなく、計算自体を成立させるための中心技術である。離散化と同じく、どの誤差をどの階層で減らすのか、という視点で理解すると、手法の選択理由が明確になる。

2.6 並列計算と再現性

材料計算は自由度が大きくなりやすく、単一計算機では時間もメモリも足りないことが多い。そこで分散メモリ並列(MPI)による領域分割が標準となる。領域 Ω を部分領域 Ωp に分け、各プロセスが自分の領域の自由度を保持する。有限差分・有限体積では、隣接格子点や隣接セルとの結合があるため、境界近傍の値を交換する通信(ハロー交換)が必要になる。有限要素法でも、共有節点や要素境界に対応する自由度の整合のために通信が必要になる。

疎行列計算では、計算の中心は行列ベクトル積 Ax である。領域分割すると、各プロセスは自分の行(自分の自由度)に対応する部分を計算できるが、隣接領域の未知量が必要になる場合がある。したがって、通信量はメッシュの境界面積に比例し、計算量は体積に比例する。プロセス数を増やすと通信の割合が増え、加速が頭打ちになることがある。強スケーリング(問題サイズ固定でプロセス数増加)と弱スケーリング(プロセスあたりの問題サイズ固定)の見方を区別して、どこが律速かを整理する必要がある。

並列計算で重要なのは、再現性が揺れる可能性である。浮動小数点の加算は結合法則を満たさず、

(a+b)+ca+(b+c)

が一般に成り立つ。並列では、総和(内積、ノルム、残差評価など)を複数プロセスで分担し、最後に集約する。この集約の順序が実行ごとに変わると丸め誤差の入り方が変わり、最終結果がわずかに変化する場合がある。線形・非線形反復は、停止判定や分岐によりこの微小差を増幅することがあるため、数値の揺れを完全にゼロにするのは難しい。

この問題への対応は、目的に応じて複数ある。第一に、停止判定を物理量に基づいて頑健にし、丸め誤差程度の揺れで分岐しないようにする。例えば相対残差の閾値を過度に厳しくしすぎず、必要精度に合わせる。第二に、可能なら決定的な還元(同じ順序で総和を取る)や高精度和(補償和)を用いて揺れを減らす。第三に、検証では一度の実行結果だけでなく、メッシュ・時間刻み・ソルバ設定を変えたときに物理的結論が不変であることを示す。

疎行列ライブラリ(PETScなど)は、並列疎行列、Krylovソルバ、前処理、非線形ソルバを統一的に提供し、現代の材料計算の基盤となっている。重要なのは、ライブラリを使うこと自体ではなく、数値解法の構造(残差、安定性、前処理、スケーリング)を理解したうえで、設定が何を意味しているかを説明できることである。並列化は計算を速くする手段であると同時に、数値誤差の扱い方を再点検する契機でもある。

小まとめ

離散化は、連続な材料方程式を有限次元の代数方程式へ変換する操作であり、その選択が安定性・精度・計算量を決める。今後の章で扱うDFT、モンテカルロ、分子動力学、フェーズフィールド、有限要素法のいずれも、内部では本章の考え方に沿って数値的に解かれているため、方程式から解法までを一つの流れとして捉えることが材料計算学の基礎になる。

第3章 密度汎関数理論の理論基盤

本章では、量子多体系としての電子問題を、電子密度という場の変数へ落とし込む論理を、定理と変分計算として整理する。計算で現れるKohn–Sham方程式が、どのエネルギー汎関数の最小化問題として導かれるのかを、式のつながりが追える形で示す。

3.1 電子密度と基礎定理

電子の基底状態は、多体波動関数 Ψ(r1,,rN) によって記述されるが、これは 3N 次元の関数であり、材料科学の立場で直に扱うには重すぎる。DFTは、この困難を電子密度

n(r)=N|Ψ(r,r2,,rN)|2dr2drN

を主変数に据えることで回避する理論である。密度 n(r) は3次元の関数であり、実空間における電荷の分布として直感的にも解釈しやすい。

Hohenberg–Kohn(HK)の第一定理は、外場ポテンシャル Vext(r) と基底状態密度 n(r) が一対一に対応することを主張する。ここで外場とは、原子核が作るクーロンポテンシャル(および外部電場など)であり、電子系のハミルトニアン

H^=T^+V^ee+Vext(r)n^(r)dr

を定める。T^ は電子の運動エネルギー演算子、V^ee は電子間相互作用、n^(r) は密度演算子である。第一定理の重要点は、外場が決まれば基底状態の密度が決まり、逆に密度が分かれば外場(定数を除いて)が決まり、したがってハミルトニアンが一意に定まるという論理である。これは、密度だけからエネルギーを評価できる汎関数 E[n] が存在することを示唆する。

HKの第二定理は、密度に関する変分原理を与える。すなわち、ある外場 Vext のもとで許される密度の集合(N 電子の v-representable な密度)に対し、エネルギー汎関数

E[n]=F[n]+Vext(r)n(r)dr

を最小にする密度が基底状態密度 n0(r) である、という主張である。ここで

F[n]=minΨnΨ|T^+V^ee|Ψ

は普遍汎関数であり、外場に依存しない形で定義される。Ψn は、その波動関数が密度 n を与えるという制約である。外場依存は最後の項 Vextn にだけ入るため、材料が変わっても共通に使える骨格が F[n] にある、という見通しが得られる。

変分原理の計算としては、粒子数保存

n(r)dr=N

を拘束条件として、ラグランジュ未定乗数(化学ポテンシャル) μ を導入し、

δ(E[n]μn(r)dr)=0

を満たす密度が平衡密度になる。ここでの汎関数微分は、後にKohn–Sham方程式の形に直結する。HK定理は、密度という場を用いた変分問題として電子基底状態が定まることを保証し、計算に向けた数学的な入口を与える。

ただし、HK定理は存在定理であり、F[n] の具体形を与えない。DFTの実用化は、次節のKohn–Sham写像によって、未知の多体系を「非相互作用系+交換相関」という形に分解し、計算可能な形に落とすことにある。

3.2 Kohn–Sham方程式

Kohn–Sham(KS)形式の要点は、相互作用する電子系の基底状態密度を、ある非相互作用電子系の密度として再現できるようにすることである。非相互作用系では、N 電子の波動関数は一電子軌道のスレーター行列式で表され、未知量は一電子軌道 {ψi(r)} になる。ここで、電子密度は占有数 fi を用いて

n(r)=ifi|ψi(r)|2

と表される(スピンを明示する場合は後述する)。

KSのエネルギー汎関数は、普遍汎関数を

F[n]=Ts[n]+EH[n]+Exc[n]

と分解して定義する。Ts[n] は非相互作用系の運動エネルギー、EH[n] は古典的クーロン反発(Hartree項)、

EH[n]=12n(r)n(r)|rr|drdr

である。残りの差分が交換相関エネルギー

Exc[n]=(T[n]Ts[n])+(Vee[n]EH[n])

としてまとめられる。ここで T[n]Vee[n] は本来の相互作用系での運動エネルギーと相互作用エネルギーの密度汎関数表現であり、厳密には未知であるが、差分を Exc として吸収することで計算可能な枠組みにする。

次に、軌道 ψi に関する変分計算でKS方程式を導く。外場を含めた全エネルギーを

E[{ψi}]=Ts[{ψi}]+EH[n]+Exc[n]+Vext(r)n(r)dr

と書く。ここで

Ts[{ψi}]=ifiψi(r)(122)ψi(r)dr

は原子単位系での表現である。軌道は正規直交条件

ψi(r)ψj(r)dr=δij

を満たす必要があるため、ラグランジュ未定乗数行列 Λij を用いて

δ(E[{ψi}]i,jΛij[ψiψjdrδij])=0

をとる。ψi で変分すると、一般に

[122+Vext(r)+VH(r)+Vxc(r)]ψi(r)=jΛijψj(r)

が得られる。ここで

VH(r)=n(r)|rr|dr,Vxc(r)=δExc[n]δn(r)

である。Λij はエルミートであり、軌道空間のユニタリ変換で対角化できるので、固有値 ϵi を用いて

[122+Vext(r)+VH[n](r)+Vxc[n](r)]ψi(r)=ϵiψi(r)

と書ける。これがKohn–Sham方程式である。見た目は一電子シュレディンガー方程式に似ているが、有効ポテンシャルが密度 n(r) に依存するため、非線形な自己無撞着問題になっている点が本質である。

自己無撞着(SCF)の構造は次のように理解できる。密度 n を仮定すると VH[n]Vxc[n] が定まり、有効ハミルトニアン

H^KS[n]=122+Vext+VH[n]+Vxc[n]

を解いて軌道 {ψi} を得る。そこから新しい密度

nout(r)=ifi|ψi(r)|2

を作り、nout=nin になるまで反復する。収束判定は、密度差 noutnin、全エネルギー差、力の大きさなどで行う。金属ではフェルミ準位近傍に多くの準位が存在するため、占有数 fi をスメアリングで連続化し、有限温度の電子自由エネルギーを用いてSCFを安定化することが多い(第4章で整理する)。

KS方程式は、量子力学の基礎式と変分原理から導かれるため、式の意味を追えることが強みである。以降の節では、この枠組みの中で未知である交換相関汎関数と、計算上の工夫(擬ポテンシャル、相対論、精度管理)を体系的に整理する。

3.3 交換相関汎関数の近似

交換相関汎関数 Exc[n] は厳密形が未知であるが、実用DFTでは近似の階層が整備されている。近似の作り方は、局所的に一様電子気体を参照する、密度の勾配情報を加える、軌道由来の情報(運動エネルギー密度など)を加える、厳密交換を混合する、といった形で体系化される。重要なのは、近似は恣意的な当てはめではなく、既知の厳密条件(自己相互作用の構造、スケーリング、均一極限、和則など)や、参照系での高精度データに基づいて構成されるという点である。

LDA(局所密度近似)は、各点で電子密度が一様電子気体であるとみなし、

ExcLDA[n]=n(r)εxcUEG(n(r))dr

と置く。εxcUEG(n) は一様電子気体の交換相関エネルギー密度であり、交換は解析的に

εxUEG(n)=34(3π)1/3n1/3

で与えられる(原子単位系)。相関は量子モンテカルロなどの参照データをパラメータ化して与える。LDAは、密度がゆっくり変化する金属的結合に強い一方、分子や界面など密度変化が大きい領域では誤差が増えやすい。

GGA(一般化勾配近似)は密度勾配を導入し、

ExcGGA[n]=f(n(r),n(r))dr

の形を取る。交換部分では、無次元勾配

s(r)=|n(r)|2kF(r)n(r),kF(r)=(3π2n(r))1/3

を用いて、LDA交換に増強因子 Fx(s) を掛ける形がよく使われる。PBEはGGAの代表であり、既知の制約条件と均一極限の整合を重視して構成されている。GGAは結合長や格子定数、表面エネルギーなどでLDAより改善する場合が多いが、物質群によって得意不得意が残るため、物性ごとに収束と妥当性を評価する姿勢が必要である。

meta-GGAは、密度と勾配に加えて、運動エネルギー密度

τ(r)=12ifi|ψi(r)|2

などの軌道由来の量を使うことで、結合様式(共有結合、金属結合、弱結合など)の識別力を上げる。SCANは多くの厳密条件を満たすように設計され、meta-GGAの中で広く参照されてきた。数値安定性を改善する工夫を加えたr2SCANのような改良も提案されており、実装と収束性の観点から選択されることがある。meta-GGAは表面・界面、遷移金属化合物などで改善する例がある一方で、ポテンシャルの非線形性が増すためSCFが不安定になりやすい場合もあり、混合法や初期密度の工夫が効きやすい。

ハイブリッド汎関数は、厳密交換(Fock交換)を混合する。厳密交換エネルギーは占有軌道から

ExHF=12i,jfifjψi(r)ψj(r)ψj(r)ψi(r)|rr|drdr

で与えられ、局所・半局所近似より非局所的な交換を取り込む。一般形として

Exchyb=aExHF+(1a)Exsemilocal+Ecsemilocal

のように書ける。バンドギャップや局在電子の記述が改善する例がある一方、計算量が増え、金属では収束が難しくなる場合がある。スクリーン型(短距離だけ厳密交換を入れる)などの工夫は、この計算負担と収束性を緩和するために導入される。

弱い分散相互作用(ファンデルワールス力)を扱うには、半局所汎関数だけでは不十分な場合がある。経験的な分散補正(C6/R6 型)や、非局所相関を含む汎関数などが用いられるが、ここでも重要なのは、どの相互作用を追加したのかを明確にし、対象物性(吸着エネルギー、層間距離など)に対して妥当な選択をすることである。

近似の階層参照量取り込む情報計算上の性質
LDAn局所一様電子気体収束しやすいことが多い
GGAn,n密度の空間変化多くの材料で基準として使われる
meta-GGAn,n,τ局所環境の識別力向上非線形性が増し収束が難しくなる場合がある
ハイブリッド厳密交換混合非局所交換計算量増大、金属で難しい場合がある

交換相関汎関数は、物性予測の中心誤差源になりやすい。したがって、同じ構造でも汎関数を変えて結果の変化を見て、議論の頑健性を確かめる姿勢が、材料物性へ接続するうえで不可欠である(第4章で相互検証の観点を整理する)。

3.4 擬ポテンシャルとPAW

平面波基底で周期系を扱う場合、原子核近傍の内殻電子は波動関数が急峻に振動するため、これをそのまま展開すると非常に高い平面波カットオフが必要になる。材料物性の多くは価電子が担うため、内殻電子を凍結し、価電子だけを有効ポテンシャルで扱う考え方が擬ポテンシャルである。目的は、価電子の散乱特性を保ったまま、波動関数を滑らかにして基底数を減らすことである。

擬ポテンシャルの基本要求は、あるカットオフ半径 rc の外側で、擬波動関数 ψ~ が全電子波動関数 ψ と一致し、価電子のエネルギーや散乱位相が再現されることである。ノルム保存擬ポテンシャルでは、さらに

0rc|ψ~l(r)|2r2dr=0rc|ψl(r)|2r2dr

のように、角運動量 l ごとにノルム(内側の確率)を一致させる条件を課す。これにより転送可能性(異なる化学環境でも再現性が保たれる性質)が高まりやすいが、滑らかさの制約が強く、カットオフが高くなる場合がある。

ultrasoft擬ポテンシャルはノルム保存条件を緩めることで、より低いカットオフで計算できるようにする。その代わり、重なり演算子が単位ではなくなり、一般化固有値問題

H^ψi=ϵiS^ψi

として解く必要がある。ここで S^ は重なり演算子であり、電荷密度も補正電荷を含む形で定義される。仕組みは複雑になるが、遷移金属や重元素を含む系で計算効率が大きく向上する場合がある。

PAW(projector augmented-wave)法は、全電子波動関数と擬波動関数を線形変換で結ぶ枠組みである。概念的には

|Ψ=T^|Ψ~

という変換演算子 T^ を導入し、擬波動関数 |Ψ~ は平面波で滑らかに表現しつつ、原子球内では部分波と射影子により全電子的な性質を復元する。計算としては、滑らかな部分で効率を確保しながら、内側の情報も取り戻せるため、ultrasoftと全電子法の中間に位置づく有力な方法である。

擬ポテンシャルやPAWでは、どの電子を価電子に含めるか(半芯状態を含めるか)が物性に影響することがある。例えば磁性や結合に寄与する準位が核近傍にある場合、半芯を凍結すると誤差が増えることがある。逆に価電子を増やすと計算量が増えるため、対象物性(格子定数、磁気異方性、電場勾配など)に応じて選ぶ必要がある。

また、擬ポテンシャルは生成条件(参照配置、交換相関汎関数、相対論取り扱い)に依存するため、DFT計算では汎関数と擬ポテンシャルの整合が重要である。汎関数だけを変えて擬ポテンシャルを据え置くと、厳密には理論整合が崩れる場合がある。材料計算ではこの点をどの程度厳密に扱うかが現実的な判断になるが、少なくとも差の出やすい物性では慎重に相互検証を行うべきである。

3.5 スピン、相対論、磁性

材料機能においてスピン自由度は中心的である。スピン分極DFTでは、電子密度をスピンごとに分けて

n(r)=n(r)+n(r),m(r)=n(r)n(r)

とし、交換相関エネルギーを Exc[n,n] として扱う。Kohn–Sham方程式もスピンごとに

[122+Vext+VH+Vxc,σ]ψiσ=ϵiσψiσ,Vxc,σ=δExcδnσ

となる。ここで σ{,} である。交換は同スピン間で強く働くため、スピン分極は交換エネルギーの利得と運動エネルギー(バンド占有)の増加の競合で決まる。この競合を粗く捉えるモデルとしてStoner模型があり、バンド構造の密度状態と相互作用の積で強磁性が生じやすい条件が理解できる(第4章で金属の扱いと併せて述べる)。

非共線磁性では、スピン量子化軸が空間的に変化するため、スピン密度はベクトル場 m(r) で扱う。Kohn–Sham方程式はスピノル軌道 Ψi(r) を未知量とし、2成分(相対論を含むと4成分)を持つ方程式になる。交換相関ポテンシャルも、スカラーではなくスピン行列として現れ、局所磁化方向に沿った形で作用する。

スピン軌道相互作用(SOC)は、磁気異方性、Dzyaloshinskii–Moriya相互作用、異常ホール効果などに直接つながる。相対論的な出発点はDirac方程式であるが、固体計算では多くの場合、Pauli近似としてSOC項が有効ハミルトニアンに現れる。中心ポテンシャルを仮定すると、原子内SOCは

H^SOC=ξ(r)L^S^

の形で書け、ξ(r) はポテンシャル勾配に比例する係数である。結晶中では、非球対称成分や混成により有効SOCが変化し、磁気異方性エネルギーの差として現れる。磁気異方性はエネルギー差が非常に小さいため、k 点、カットオフ、収束条件が強く効きやすく、数値精度の管理が重要になる。

相対論効果はSOCだけではない。重元素ではスカラー相対論(質量増大や収縮効果)が結合長やバンド位置に影響する。擬ポテンシャルやPAWの生成時に相対論を含めたものが提供されるのは、この効果を価電子有効ポテンシャルに取り込むためである。SOCを後から摂動的に入れる方法(第二変分法)を用いる実装も多く、計算負担と精度のバランスを取る設計になっている。

磁性材料を扱うときは、スピン分極の初期条件(初期磁気モーメント、反強磁性配置など)が収束先を変えることがある。これはエネルギー地形が多谷性を持つためであり、どの磁気秩序が安定かを比較するには、複数の初期条件とセル設定で相互に確認する必要がある。SOCを含めるとさらに自由度が増え、SCFが難しくなる場合があるため、段階的に計算設定を強めていくのが合理的である。

3.6 精度管理と相互検証

DFTは原理計算であるが、数値計算である以上、近似と離散化の影響を管理しなければならない。誤差の源は大きく分けて、交換相関汎関数の近似、擬ポテンシャルやPAWの近似、基底と離散化(平面波カットオフ、実空間グリッド、数値積分)、ブリルアンゾーン積分(k 点)、占有数の扱い(スメアリング)、有限サイズ(スーパーセル)である。どの誤差が対象物性に効くかを理解し、収束試験で定量化する必要がある。

平面波DFTでは、波動関数の展開をエネルギーカットオフ Ecut で打ち切る。平面波数は概ね体積と Ecut3/2 に比例して増えるため、Ecut を上げるほど計算量が増える。全エネルギーは比較的早く収束する場合が多いが、応力や力、特に微小差(磁気異方性や欠陥形成エネルギーの差)はより厳しい収束が必要になることがある。したがって、目的物性ごとに収束基準を定める姿勢が重要である。

k 点は周期系の積分を数値的に近似するために必要であり、金属や低次元系では特に敏感である。k 点密度を上げるとフェルミ面近傍の寄与が精密になるが、SCFの不安定性が増す場合もある。スメアリング幅 σ は収束を助けるが、エネルギーを人工的に丸めるため、σ の依存性を見て物理値へ外挿する考え方も重要になる(第4章で式を示す)。

スーパーセル計算では、周期境界により欠陥像どうしが相互作用する。電荷欠陥では、長距離クーロン相互作用により収束が遅く、電荷補正やポテンシャル整列が必要になることがある。界面や表面でも、真空層厚さ、スラブ厚さ、双極子補正などの設計が結果を左右する。これらは数値条件の選択に見えるが、根は周期境界という数学的仮定にあるため、相互作用の長さスケールを見積もることが本質である。

相互検証は、誤差の切り分けに有効である。例えば、汎関数を変える、擬ポテンシャル(あるいはPAWデータセット)を変える、コードを変えるといった比較により、物理由来の差と数値由来の差を分離できる。全エネルギーの絶対値はコード間で直接比較できない場合が多いが、同一手順で評価したエネルギー差(形成エネルギー、相安定差など)は比較しやすい。比較は単に一致を確認する作業ではなく、どの仮定が結論に効いているかを理解する作業である。

本章の結論は、DFTが定理と変分原理で立ち上がる厳密な枠組みを持つ一方で、実用では交換相関近似と数値設定が物性値を左右するという点である。次章では、DFT計算を材料物性へ結びつけるために、スケーリング、ブリルアンゾーン積分、構造最適化、格子力学、欠陥・界面、有限温度といった具体的課題を、式と計算の観点から整理する。

第4章 DFT計算から材料物性へ

本章では、Kohn–Sham方程式を解いた結果を、材料科学で議論する物性(構造安定性、弾性、フォノン、電子状態、欠陥形成、界面エネルギー、有限温度自由エネルギーなど)へ変換する手順を、数学的定義と数値上の要点として整理する。DFTの出力をそのまま読むのではなく、何を計算し、どの差を物性として解釈するのかを、式の形で明確にする。

4.1 基底関数と計算スケーリング

Kohn–Sham方程式は連続な偏微分方程式であるため、何らかの基底に展開して有限次元の固有値問題として解く。平面波基底は周期境界と相性が良く、基底をエネルギーカットオフで系統的に増やせるため、収束の扱いが比較的分かりやすい。周期セル体積を V とし、最大平面波ベクトルの大きさ Gmax を取ると、平面波数は概ね

NPWVGmax3VEcut3/2

に比例して増える。固有値問題の計算は概ね行列ベクトル積と直交化が中心であり、バンド数が増えるほど計算負担が増大する。一般に、バンド法のコストは系サイズに対して O(N3) 的に増える部分を含み、大規模系では計算が難しくなる。

局在基底(原子軌道基底、数値原子軌道、ガウス基底など)は、原子の周りに局在した関数で波動関数を展開する。基底数は原子数に比例しやすく、疎行列構造が得られるため、大規模系で有利になりやすい。数値局在基底を用いると、行列要素の評価が実空間積分として整理でき、線形スケーリング法(密度行列の局在性を利用して O(N) に近づける手法)へ接続しやすい。ここでの核心は、絶縁体やギャップのある系では密度行列

ρ(r,r)=ifiψi(r)ψi(r)

|rr| に対して指数的に減衰しやすいという性質である。局在性を仮定できると、遠距離の結合を切り捨てる近似が成立し、計算量が大きく減る。

ただし、局在基底は基底の選び方が収束性に影響しやすい。平面波では Ecut を上げれば収束するという単調性があるが、局在基底では基底の半径、基底関数の数、分極関数の追加など、複数の選択がある。したがって、局在基底を使う場合は、対象物性に対して基底の系統的改良ができるように設計し、収束の見通しを持って使う必要がある。

計算スケーリングの理解は、材料物性の対象に応じて手法を選ぶ判断につながる。例えば、フォノンや欠陥でスーパーセルが必要な場合、平面波で NPW が増えすぎると現実的でなくなり、局在基底や混合法(平面波+局在補助基底など)を検討する余地が出る。逆に、表面や界面で高い精度のエネルギー差が必要な場合、平面波の系統的収束が強みになる場合がある。基底は単なる道具ではなく、誤差の形を決める要素である。

4.2 ブリルアンゾーン積分と金属の収束

周期系の電子状態はブロッホ波で記述され、波数 k を持つ状態として整理される。全エネルギーや電子密度は、ブリルアンゾーン(BZ)での積分として表される。例えば、密度は

n(r)=nBZfnk|ψnk(r)|2dkΩBZ

であり、エネルギーも同様に k 積分を含む。数値計算ではこの積分を離散化し、k 点サンプリングで近似する:

BZg(k)dkkwkg(k)

ここで wk は重みである。Monkhorst–Pack格子は対称性を活かしつつ一様サンプリングを行う基本的な方法であり、結晶対称性を用いて既約ブリルアンゾーンへ縮約することで計算量を減らせる。

金属の難しさは、フェルミ面近傍で占有数が急に変わり、積分の被積分関数が滑らかでなくなる点にある。絶縁体では占有がギャップで分離されるため、k 点収束は比較的穏やかであるが、金属では k 点密度を上げても揺れが残りやすい。そこで、占有数を有限温度のフェルミ分布

f(ϵ)=11+exp(ϵμkBT)

で滑らかにし、数値積分を安定化する方法が用いられる。これはスメアリングと呼ばれ、実際には T を物理温度として解釈するより、数値平滑化の幅として使うことが多い。スメアリングを導入すると、評価すべき熱力学量は内部エネルギー E ではなく、電子自由エネルギー

F=ETS

になる。ここでエントロピー項 S は占有数から

S=kBnkwk[fnklnfnk+(1fnk)ln(1fnk)]

で与えられる。T0 の極限で F は基底状態エネルギーへ戻るので、スメアリング幅を変えたときに F の変化が小さい領域を探すことで、金属の収束を評価できる。

テトラへドロン法は、BZを小さな四面体に分割し、エネルギー分散を線形補間して積分する方法である。絶縁体や半導体では高精度に積分でき、金属でも改良版(占有の扱いを工夫した方法)が用いられることがある。ただし、磁性金属や準二次元系など、フェルミ面が複雑で状態密度が鋭く変化する場合には、k 点とスメアリングの両方が強く効く。これらの系では、磁化や交換分裂の微小差が安定性を左右するため、収束をエネルギーだけでなく磁気モーメントや応力でも確認する必要がある。

低次元系(薄膜、表面、界面)では、周期方向が限られるため、k 点サンプリングは方向ごとに異なる密度が必要になる。例えば、スラブ計算では面内は密に、真空方向は1点でよい、という設計になる。BZ積分は単なる数値設定ではなく、結晶対称性と電子状態の滑らかさに根ざした問題であり、対象系の次元性と金属性を見て設計する必要がある。

4.3 構造最適化と格子力学

DFT計算で最初に行うことが多いのは、原子位置 {RI} と格子ベクトルを最適化し、力学的平衡構造を得ることである。基底状態エネルギー E({RI}) が定義できるとき、力は

FI=ERI

で定義され、平衡条件は FI=0 である。Kohn–Sham方程式が満たされているとき、力はHellmann–Feynmanの形で表されるが、基底が原子位置に依存する場合(局在基底など)にはPulay力が加わる。平面波基底では基底が原子位置に依存しないためPulay力が生じにくいが、有限カットオフでは応力にPulay成分が残ることがある。

構造最適化は、数値的には多変数最小化問題であり、共役勾配法、準ニュートン法(BFGSなど)、減衰付き分子動力学などが用いられる。これらは、力と(場合により)近似ヘッセ行列から更新量を決め、エネルギーを下げる方向へ進む。重要なのは、SCF収束が不十分だと力がノイズを含み、最適化が不安定になる点である。したがって、力収束閾値とSCF収束閾値は連動して設定する必要がある。

格子力学は、平衡構造の周りで原子を微小変位させたときの二次の応答からフォノンを求める枠組みである。原子 I の変位を uI とすると、エネルギーは

EE0+12Iα,JβΦIα,JβuIαuJβ

と展開できる。Φ は力定数であり、

ΦIα,Jβ=2EuIαuJβ=FIαuJβ

で定義される。質量 MI を用いて力定数を正規化すると、動力学行列

DIα,Jβ(q)=1MIMJRΦIα,Jβ(R)eiqR

が得られ、固有値問題

JβDIα,Jβ(q)eJβ(s)(q)=ωs2(q)eIα(s)(q)

を解くことでフォノン分散 ωs(q) が得られる。ω2<0(虚数周波数)が現れる場合、平衡構造がその波数の変位に対して不安定であり、相転移や歪み緩和の駆動を示す。

力定数は有限変位法(原子を少し動かして力差分から求める)でも、密度汎関数摂動論(DFPT)でも求められる。有限変位法では、変位量 δ を与えたとき

ΦIα,JβFIα(uJβ=δ)FIα(uJβ=δ)2δ

のように数値微分する。δ が小さすぎると数値ノイズが増え、大きすぎると非調和成分が混ざるため、適切な範囲を選ぶ必要がある。DFPTは線形応答として導くため、原理的に効率が良いが、実装や擬ポテンシャルとの整合が重要になる。

フォノンから熱力学へ接続するには、調和近似の振動自由エネルギー

Fvib(T)=kBTq,sln[2sinh(ωs(q)2kBT)]

を用いる。これにより、温度依存の自由エネルギー F(T)=E0+Fvib(T) を比較して相安定性を議論できる。さらに体積依存 ωs(q;V) を考慮する準調和近似により、熱膨張や温度依存弾性へ接続できる。

4.4 励起状態と強相関の拡張

標準DFTは基底状態理論であり、Kohn–Sham固有値 ϵnk は形式的には補助系の量である。バンド構造として可視化されるが、光学ギャップや準粒子ギャップを厳密に与える量ではない。励起状態を扱うには、基底状態DFTを拡張する必要がある。

GW法は、多体摂動論に基づき、自己エネルギー Σ

Σ(r,r;ω)=i2πG(r,r;ω+ω)W(r,r;ω)dω

で近似する(G はグリーン関数、W は遮蔽クーロン相互作用)。準粒子エネルギー EnkQP は、概念的には

EnkQP=ϵnkKS+ψnk|Σ(EnkQP)Vxc|ψnk

を解くことで得られる。これによりバンドギャップが改善する例が多いが、計算量が大きく、収束(空虚準位数、誘電関数のカットオフ、k 点など)が難所になる。

BSE(Bethe–Salpeter方程式)は、電子と正孔の相互作用を取り込み、励起子効果を含む光学応答を記述する。線形応答として、電子正孔対の有効ハミルトニアンを作り、その固有値として励起子準位を得る枠組みである。TDDFT(時間依存DFT)は、密度の時間発展と線形応答を扱い、光吸収スペクトルなどへ接続する。TDDFTの精度は交換相関カーネルの近似に依存し、長距離相互作用や励起子を精密に扱うには工夫が必要になる。

強相関系では、局在した d,f 電子の相互作用が大きく、半局所DFTが金属化や磁気秩序の誤りを生む場合がある。DFT+Uは、局在軌道の占有行列 nmmσ を導入し、Hubbard相互作用を平均場的に補正する。簡略化した表現では

EDFT+U=EDFT+Ueff2m,σ(nmσ(nmσ)2)

のように書け、部分占有を罰することで局在化を促す。ここで二重計数補正(DFT側で既に含まれている相互作用を引く項)の選び方が結果に影響し、U の値も物性に依存して調整が必要になる。

DMFT(動的平均場理論)は、局所相互作用を動的自己エネルギーとして扱い、温度依存やスペクトル関数の形を改善する。DFT+DMFTは、DFTのバンド構造とDMFTの局所相関を結合し、金属絶縁体転移や有限温度磁性などをより物理的に扱える。ただし、計算負担が大きく、モデル部分空間の定義や二重計数の扱いが難しい。

これらの拡張は、何を追加し、何を近似しているかを明確にしたうえで選ぶ必要がある。基底状態の構造安定性だけなら標準DFTで十分な場合も多いが、光学応答、準粒子ギャップ、局在相関、有限温度スペクトルなどを議論するなら、拡張法が必要になる。

4.5 欠陥、界面、有限温度

材料物性の多くは、完全結晶ではなく欠陥や界面で決まる。DFTで欠陥を扱う基本はスーパーセル法であり、周期セル内に欠陥を導入して周期境界で繰り返す。欠陥形成エネルギーは、化学ポテンシャルと電荷状態を含めて

Ef(Dq)=Etot(Dq)Etot(bulk)iniμi+q(EF+Ev)+Ecorr

のように定義される。Etot(Dq) は欠陥を含むセルの全エネルギー、Etot(bulk) は欠陥なしの全エネルギー、ni は出入りした原子数、μi は化学ポテンシャル、q は欠陥電荷、EF はフェルミ準位、Ev は基準(価電子帯上端など)、Ecorr は有限サイズ補正である。電荷欠陥では、周期像間のクーロン相互作用が長距離で残りやすく、ポテンシャル整列や補正が必要になる。補正の具体形はモデルに依存するが、どの相互作用が周期境界で人工的に入ったかを意識して設計することが本質である。

界面・表面では、エネルギーを面積で割った量が重要になる。表面エネルギーはスラブの全エネルギー Eslab を用いて

γ=EslabNEbulk2A

のように定義できる(両面を持つ対称スラブの場合)。ここで A は表面面積である。界面エネルギーも同様に、界面を含むセルのエネルギーからバルク成分を差し引いて面積で規格化する。界面では応力緩和や電荷移動、分極などが起き、セル厚さや真空層、双極子補正が結果に効きやすい。したがって、面積当たりエネルギーが厚さに対して収束していることを確認し、界面が二つ入る場合は配置が独立であることを確かめる必要がある。

有限温度では、基底状態エネルギー E0 だけでは相安定性が決まらない。自由エネルギー

F(T)=E0+Fvib(T)+Fel(T)+

を比較する必要がある。Fvib は前節のフォノンから評価でき、Fel は金属で重要になる電子自由エネルギーである。スメアリングで導入した F=ETS の枠組みは、有限温度の電子寄与として解釈できるが、数値平滑化として使っている場合は、物理温度との対応に注意が必要である。熱的無秩序や相転移のように調和近似が破れる場合、第一原理分子動力学(AIMD)でサンプリングし、平均エネルギーや構造相関を調べる方法が有効になる。自由エネルギーを精密に求めるには熱力学積分などの方法が用いられるが、基礎としては、どの自由度(格子振動、電子励起、磁気揺らぎなど)を含めているかを明確にすることが重要である。

4.6 計算コードと国内外の基盤

DFT計算は多様なコードで実行でき、実装の選択は基底、擬ポテンシャル、機能(DFPT、SOC、非共線、ハイブリッド、GWなど)、並列性能に依存する。平面波擬ポテンシャルDFTは周期系に強く、材料科学で標準的に使われる機能が整備されていることが多い。局在基底DFTは大規模系や線形スケーリングに接続しやすく、欠陥・界面・ナノ構造の大規模計算で有利になる場合がある。どのコードを用いる場合でも、基底と擬ポテンシャルの整合、収束試験、相互検証が物性の信頼性を支える。

国内の計算基盤を利用する場合、利用可能なソフトウェア、コンパイラやMPI、推奨設定が環境ごとに異なる。大規模計算では、疎行列ソルバやFFTライブラリなどの性能が計算時間を左右するため、コードの実装特徴と計算機の構成(CPU、メモリ帯域、ネットワーク)との相性も重要になる。計算は再現可能であるべきなので、入力ファイル、擬ポテンシャルの版、汎関数、k 点、カットオフ、収束条件を記録し、別環境で同等の結果が得られることを確認する姿勢が必要である。

また、物性へ接続するには、単一の計算結果だけでなく、構造の同定、相の比較、欠陥濃度の仮定、温度・化学ポテンシャルの設定など、熱力学的な枠組みの中で整理する必要がある。コードはその一部を担う道具であり、材料科学としての結論は、何を比較し、どの自由度を含めた自由エネルギーを議論したか、という定義の上に立つ。

小まとめ

第3章では、HK定理と変分原理からKohn–Sham方程式が導かれ、未知の多体相互作用が交換相関汎関数へ集約される構造を示した。第4章では、DFTの出力を物性へ変換するために、基底とスケーリング、BZ積分、構造最適化とフォノン、励起状態拡張、欠陥・界面・有限温度といった論点を式の形で整理した。今後は、これらのDFTで得たエネルギーや力、応答関数を、モンテカルロ法、分子動力学、フェーズフィールド、有限要素法などの多階層計算へ接続し、有限温度・有限濃度・非平衡の材料現象を記述する方向へ進むことになる。

第5章 モンテカルロ法の基礎

本章では、熱平衡の分布を数値的にサンプリングするという立場から、モンテカルロ法を統計力学と確率過程の言葉で組み立てる。材料科学で必要になるのは、単に乱数で平均を取ることではなく、所望の分布に従う標本列を設計し、その誤差を制御して物理量へ結び付けることである。

5.1 統計集団とボルツマン重み

温度 T の熱平衡で、系が状態 x(原子配置、スピン配置、格子占有、連続変数の集合など)にある確率は、統計集団により決まる。最も基本となる正準集団(カノニカル集団)では、エネルギー E(x) に対して

π(x)=1Zexp[βE(x)],β=1kBT

である。ここで分配関数 Z は規格化定数であり、

Z=xexp[βE(x)]

(離散状態の場合)または

Z=exp[βE(x)]dx

(連続状態の場合)で定義される。モンテカルロ法の核心は、Z を直接求めなくても π(x) に従う標本 x(1),x(2), を生成し、観測量 A(x) の熱平均

A=xA(x)π(x)

(または積分)を

A1Nt=1NA(x(t))

として推定する点にある。

材料科学では、温度以外の制御変数が重要になる。例えば、化学ポテンシャル μ を固定する大正準集団では、粒子数 N が変動し、

π(x,N)=1Ξexp[β(E(x,N)μN)]

となる。圧力 P を固定する等温等圧集団(NPT)では体積 V が変動し、エンタルピーに相当する量が重みに入る。これらは、相図や組成揺らぎ、欠陥濃度を議論するときに必要になるため、まず「どの集団の分布をサンプリングしたいのか」を定めることが出発点である。

さらに、材料の相転移や秩序化は、秩序変数の揺らぎとして観測される。秩序変数 m(x)(例:磁化、長距離秩序パラメータ、結晶方位ラベルの分布など)の分布

P(m)=xπ(x)δ(mm(x))

を推定できれば、平均値だけでなく相共存や自由エネルギー障壁も議論できる。したがって、モンテカルロ法は「平均を出す道具」であると同時に「分布を推定し、相空間の形を読む道具」である。

5.2 Metropolis法と詳細釣り合い

正準分布 π(x) を直接サンプルするのは一般に難しいため、マルコフ連鎖(Markov chain)を用いて標本列を生成する。状態 x から x への遷移確率を P(xx) とすると、分布 ρt(x) の時間発展は

ρt+1(x)=xρt(x)P(xx)

で与えられる。目標は、十分大きい tρt(x)π(x) となるように P を設計することである。このとき、π が不変分布(stationary distribution)である条件は

π(x)=xπ(x)P(xx)

である。これを直接満たす P を設計するのは難しいが、より強い条件として詳細釣り合い(detailed balance)

π(x)P(xx)=π(x)P(xx)

を課すと、不変分布の条件が保証される(両辺を x で和を取ればよい)。材料計算でMetropolis法が広く使われるのは、この詳細釣り合いを満たす形で受理確率が簡潔に構成できるからである。

実装上、遷移は「提案」と「受理」に分けて書くと理解しやすい。まず、提案分布 q(xx) に従って候補 x を生成し、次に受理確率 a(xx) で受理する。したがって

P(xx)=q(xx)a(xx)(xx)

であり、棄却された場合は x に留まるため

P(xx)=1xxq(xx)a(xx)

となる。詳細釣り合いは

π(x)q(xx)a(xx)=π(x)q(xx)a(xx)

であるから、受理確率の比は

a(xx)a(xx)=π(x)q(xx)π(x)q(xx)=exp[β(E(x)E(x))]q(xx)q(xx)

となる。これを満たす具体形として、Metropolis–Hastingsの受理確率

a(xx)=min{1,exp[β(E(x)E(x))]q(xx)q(xx)}

が得られる。特に提案が対称 q(xx)=q(xx) のとき、比が1となり、

a(xx)=min{1,exp[βΔE]},ΔE=E(x)E(x)

がMetropolis法の基本形である。エネルギーが下がる遷移(ΔE0)は必ず受理され、上がる遷移も確率的に受理されるため、局所最小に固定されずに熱揺らぎを表現できる。

もう一つの要件はエルゴード性である。これは、許される状態空間の任意の状態へ、有限回の遷移で到達できることを要求する。例えばスピン模型で「必ず1個だけ反転する」更新は多くの場合エルゴードであるが、制約条件(粒子数固定、電荷中性など)を入れると更新が到達不能な領域を作ることがある。したがって、詳細釣り合いだけでなく、更新規則が状態空間を十分に連結しているかも設計の一部である。

5.3 誤差評価と自己相関

モンテカルロで得られる標本列 x(t) はマルコフ連鎖であり、一般に独立同分布ではない。したがって、推定量

A¯=1Nt=1NA(x(t))

の分散は、独立標本の 1/N スケーリングからずれる。これを定量化するために、自己相関関数を導入する。平衡に達した後の時系列 At=A(x(t)) に対し、

C(k)=(AtA)(At+kA)

を定義し、正規化した自己相関

ρ(k)=C(k)C(0)

を考える。ここで C(0)=Var(A) である。すると、十分大きい N に対して推定量の分散は

Var(A¯)C(0)N[1+2k=1ρ(k)]=C(0)N2τint

と書ける。ここで

τint=12+k=1ρ(k)

を統合自己相関時間と呼ぶ。独立標本なら ρ(k)=0k1)であり τint=1/2 となるので、Var(A¯)C(0)/N が回復する。相関がある場合は 2τint>1 となり、実効的な独立サンプル数が

NeffN2τint

に減ると解釈できる。材料の相転移近傍や巨大欠陥集団の更新では τint が大きくなりやすく、単に試行回数を増やすだけでは誤差が下がりにくい。

実務上は、ρ(k) を無限和で評価する代わりにブロック平均がよく用いられる。長さ L のブロックに分け、各ブロック平均 A¯(b) を計算すると、L が自己相関時間より十分大きければブロック平均どうしはほぼ独立となる。ブロック数を B とすると、ブロック平均の分散から誤差を推定できる。L を増やして推定誤差が安定する領域を探すことで、自己相関の影響を抑えた評価が可能になる。

さらに、初期条件の影響(バーンイン)も誤差源である。初期状態が平衡分布から遠いと、最初の区間は平衡サンプルとして扱えない。バーンイン長は、エネルギーや秩序変数の時系列が定常になったか、複数の初期条件(秩序相・無秩序相など)から始めても同じ分布へ収束するか、といった観点で評価する必要がある。誤差評価は、数値計算の末尾処理ではなく、サンプリング設計と同じ重みを持つ。

5.4 クラスタ更新

臨界点近傍では相関長が発散的に大きくなり、局所更新(1スピン反転など)では大域的な緩和が極端に遅くなる。この現象が臨界減速であり、自己相関時間が系サイズ L に対して

τintLz

のように増大する(z は動的臨界指数)。材料の相転移や秩序化を有限サイズで調べる場合、臨界減速は統計誤差を支配しうるため、更新規則そのものを変える必要がある。

クラスタ更新は、相関の強い自由度をまとめて反転することで、長波長の緩和を速める方法である。最も基本的な例は強磁性Ising模型

E({si})=Ji,jsisj,si=±1,J>0

である。温度が低いと同符号のスピンがドメインを形成し、境界だけが揺らぐ。局所更新は境界を少しずつ動かすしかないため遅いが、クラスタ更新では同符号の連結成分を一括で反転できる。

Wolff法では、ある種の確率で結合(ボンド)を張り、クラスタを成長させて反転する。隣接スピンが同符号である辺に対して、ボンドを張る確率を

p=1exp(2βJ)

とし、1点から始めて同符号隣接へ確率 p でボンドを追加してクラスタを作る。得られたクラスタのスピンを全て反転する。重要なのは、この p の選び方が詳細釣り合いを満たすように設計されている点である。直感的には、クラスタ境界に現れるエネルギー変化を、ボンド生成の確率と反転の確率が打ち消すように調整している。

Swendsen–Wang法は、系全体に対して同様にボンドを張り、得られた全クラスタを独立に確率 1/2 で反転する。Wolff法は1クラスタ、Swendsen–Wang法は全クラスタという違いがあるが、いずれも相関長に沿った更新を行う点が本質である。臨界近傍での動的指数 z が局所更新より小さくなることが知られており、同じ計算時間で有効独立サンプル数を大きくできる。

ただし、クラスタ法は万能ではない。強い異方性、外場、フラストレーション、連続自由度(ベクトルスピン)などでは拡張が必要であり、場合によってはクラスタが巨大化しすぎて効率が落ちる。材料科学では「相関の構造を見て更新を選ぶ」ことが基本であり、クラスタ更新はその代表例である。

5.5 拡張アンサンブル

自由エネルギー地形が多峰性を持つと、単一温度の正準サンプリングは特定の谷に閉じ込められやすい。相共存、一次相転移、核生成、ガラス的遅い緩和などでは、この問題が顕在化する。拡張アンサンブルは、分布そのものを拡張・混合し、状態空間をより広く探索することで、障壁越えを促進する枠組みである。

レプリカ交換法(平行テンパリング)は、異なる温度 T1<T2<<TM の複数レプリカを同時に走らせ、隣接温度間で状態を交換する。レプリカ m の状態を xm とすると、全体系の目標分布は積分布

Π({xm})=m=1M1Z(βm)exp[βmE(xm)]

である。ここで温度交換(xmxm+1 を入れ替える)を提案したときの受理確率は、詳細釣り合いから

Pacc=min{1,exp[(βmβm+1)(E(xm+1)E(xm))]}

となる。高温レプリカは障壁を越えやすく、低温レプリカは目的温度での分布を精密に記述する。交換により低温レプリカも間接的に障壁越えを行えるため、一次転移近傍などで有効である。温度列の設計は受理率と拡散速度を左右し、エネルギー揺らぎの大きさに応じて温度間隔を狭める必要がある。

多重カノニカル法は、エネルギー分布を平坦化するような重み W(E) を導入し、

πmu(x)W(E(x))

でサンプリングする。理想的には状態密度 g(E) を用いて W(E)1/g(E) とすると、エネルギー分布 P(E)g(E)W(E) が一定となり、相共存のエネルギー障壁を横断しやすくなる。実際には g(E) は未知であるため、反復的に重みを学習する手法(密度状態推定)と組み合わせることが多い。材料の一次相転移や核生成障壁を自由エネルギーとして引き出したい場合、この考え方が直接効く。

傘サンプリング(umbrella sampling)は、秩序変数 m(x) に対してバイアス Ub(m) を加え、

πb(x)exp[β(E(x)+Ub(m(x)))]

でサンプリングする。目的は、普段はほとんど出現しない m の領域(遷移状態付近など)を意図的にサンプルし、後で再重み付けして元の分布へ戻すことである。バイアスをどう設計するかが効率を支配するため、複数窓に分けて重ね合わせる方法がよく使われる。

拡張アンサンブルは「分布を変えて探索し、再重み付けで物理量へ戻す」という思想で統一できる。したがって、どの分布をどの領域で平坦化したいのか(温度、エネルギー、秩序変数、反応座標)を明確にすることが設計の中心となる。

5.6 自由エネルギー推定

自由エネルギーは、相安定性や相境界、欠陥濃度、界面・核生成障壁を議論する上で中心的である。正準集団の自由エネルギーは

F(T)=kBTlnZ

で定義されるが、Z は巨大な状態和であり、直接評価できない。そこで、自由エネルギー差 ΔF を推定する方法が多数用意されている。重要なのは、自由エネルギー差は「平均エネルギー差」ではなく「対数の差」であるため、推定には分布の重なりやサンプリングの偏りが強く効く点である。

自由エネルギー摂動(free energy perturbation)は、基準系 0 と対象系 1 のエネルギー E0(x),E1(x) を用いて

ΔF=F1F0=kBTlnexp[β(E1(x)E0(x))]0

を与える。ここで 0 は基準系の分布での平均である。この式は厳密であるが、E1E0 の分布が広いと指数平均が少数の希な事象に支配され、分散が爆発しやすい。したがって、二つの系の分布が十分重なる範囲で使う必要がある。

熱力学積分(thermodynamic integration)は、連続的な結合パラメータ λ[0,1] を導入し、

E(x;λ)=(1λ)E0(x)+λE1(x)

のように補間して

ΔF=01E(x;λ)λλdλ

で評価する。補間経路を適切に選べば、各 λ で分布の重なりが保たれ、安定した推定が可能になる。材料科学では、結晶相間の自由エネルギー差、界面生成、欠陥生成などを「仮想経路」でつなぎ、積分として評価する考え方が有効である。

再重み付け(reweighting)は、ある温度や重みで得たサンプルから近傍条件の物理量を推定する。正準分布の温度変更では

Aβ=A(x)exp[(ββ)E(x)]βexp[(ββ)E(x)]β

が成り立つ。これは便利であるが、β が離れるほど指数重みが極端になり、実効サンプル数が減るため、適用範囲は狭い。多温度データを統合する方法としてWHAM(重み付きヒストグラム解析)などが用いられ、傘サンプリングの窓統合にも使われる。

自由エネルギー推定法は、どれか一つを暗記するのではなく、(i) どの分布でサンプルしているか、(ii) 目的の自由エネルギー差がどの対数比として定義されるか、(iii) 分布重なりが十分か、という観点で選ぶべきである。以下に、材料科学で出番が多い推定法を比較する。

推定法基本式強み注意点
自由エネルギー摂動ΔF=kTlneβΔE0形式が簡潔で実装しやすい分布重なりが悪いと破綻しやすい
熱力学積分ΔF=01E/λλdλ経路を作れば安定しやすいλ 点数と経路設計が必要である
再重み付け比の形で β へ変換近傍条件の補間に有効である離れた条件には使いにくい
傘サンプル+統合P(m) を窓で推定障壁近傍を直接サンプルできる反応座標選定と窓設計が効く

本章の結論は、モンテカルロ法は「所望分布のサンプリング設計」と「相関を含む誤差評価」の両輪で成立しているという点である。次章では、この枠組みが材料科学の具体問題(合金秩序、相図、拡散・析出、粒成長など)へどう接続されるかを、モデル化と手法選択として整理する。

第6章 材料科学におけるモンテカルロ

本章では、材料の相転移・秩序化・拡散・成長といった現象を、確率変数と有効ハミルトニアンの選び方として定式化し、モンテカルロ法で何が計算できるかを具体化する。DFTやMDが局所エネルギーやミクロな運動を扱うのに対し、MCは平衡分布と統計揺らぎを正面から扱うため、対象スケールと目的物性に応じた連携が鍵になる。

6.1 格子模型と秩序変数

相転移を理解する第一歩は、自由度を最小限に抽出した格子模型である。代表はIsing模型であり、格子点 i にスピン si=±1 を置き、

E({si})=Ji,jsisjhisi

で定義される。J>0 は強磁性的相互作用、h は外場である。温度が高いと熱揺らぎが支配して si0(無秩序相)になり、温度が低いと相互作用が支配して si が揃い si0(秩序相)になる。

秩序変数は、相を区別する量として定義される。Ising模型では磁化

m=1Nisi

が秩序変数である。材料の秩序化(例えばB2秩序、L1_0秩序など)でも、サブラティス占有の差として同様の秩序パラメータが定義できる。秩序変数の平均値だけでなく揺らぎも重要であり、磁化の分散から磁化率が

χ=βN(m2m2)

で与えられる。これは、線形応答(h に対する m の変化)と揺らぎが結び付く揺らぎ散逸関係の一例である。

相関長を議論するには、二点相関関数

G(r)=s0srs2

を用いる。無秩序相では G(r) は指数的に減衰し、秩序相では長距離相関が残る。臨界点近傍では相関長 ξ が大きくなり、有限サイズ効果が顕著になる。有限サイズ L の系では、真の発散の代わりに「ピーク位置が L に依存して動く」形で現れるため、温度走査とサイズ走査を組み合わせて臨界点を推定する必要がある。

Potts模型は、q 状態のスピン σi{1,,q} を用いて

E({σi})=Ji,jδσi,σj

と書ける。これは多結晶の粒成長モデル(粒の方位を q 状態で表す)や、多相系の粗視化にも自然に現れる。格子模型は材料を直接再現するのではなく、対称性、秩序変数、相関長という言葉を数式で扱う訓練として有効であり、後の合金秩序や粒界問題へ直結する。

6.2 クラスター展開と置換合金

置換合金では、原子が格子点を占有し、どの元素がどこに入るか(配置)が自由度になる。温度により短距離秩序が生じ、相分離や秩序化が起きるが、これを第一原理のエネルギーと統計力学でつなぐ代表がクラスター展開である。二元合金 AB を考え、格子点 i の占有を擬スピン

σi={+1(A)1(B)

で表す。すると、配置 {σi} に対するエネルギーを、クラスター(点、対、三体…)の基底で展開して

E({σ})=J0+iJiσi+i<jJijσiσj+i<j<kJijkσiσjσk+

と書ける。結晶対称性を用いると、同等なクラスターは同じ係数を共有するため、実際にはクラスターの種類 α とその相関関数 Φα(σ) を用いて

E(σ)=αJαΦα(σ)

と表すことが多い。Jα は有効クラスター相互作用(ECI)であり、エネルギーランドスケープの圧縮表現である。

Jα は、複数の代表配置 σ(p) に対してDFTで全エネルギー EpDFT を計算し、最小二乗などで同定する:

min{Jα}p(EpDFTαJαΦα(σ(p)))2+λR({Jα})

ここで R は正則化(過学習抑制)であり、λ は重みである。配置数が有限である以上、モデル選択(どのクラスターまで入れるか)と検証(学習に使わない配置で誤差を見る)が重要になる。クラスター展開が材料科学で価値を持つのは、(i) DFTの高精度エネルギーを用いながら、(ii) 任意の大規模配置のエネルギーを高速に評価でき、(iii) その上でMCで有限温度の熱平均を取れる点にある。

クラスター展開が得られると、正準MCで温度依存の秩序変数や比熱を計算できる。例えば内部エネルギーの熱平均 E と比熱は

CV=ET=kBβ2(E2E2)

で評価でき、相転移点ではピークとして現れる。組成が変動する大正準MCを使えば、化学ポテンシャル μ を制御して吸着等温線や相境界を求めることもできる。短距離秩序は対相関 σiσj の距離依存として測れ、実験の散漫散乱や局所構造と比較できる形になる。

この枠組みの落とし穴は、エネルギーだけでなく、格子歪みや磁気自由度が重要な材料では、有効ハミルトニアンが単純な置換自由度だけでは閉じない点である。必要に応じて、格子緩和を取り込んだ拡張クラスター展開や、磁気自由度を含めたスピン模型、連続自由度を含む有効モデルへ拡張することになる。重要なのは、どの自由度を「確率変数」に含め、どの自由度を「エネルギー関数」に吸収したのかを明確にすることである。

6.3 Kinetic Monte Carlo

拡散、析出、相変態の核生成、表面成長などの非平衡過程は、時間発展そのものが重要であり、平衡MCでは扱えない。Kinetic Monte Carlo(kMC)は、遷移率に基づく連続時間マルコフ過程として、イベント駆動で時間を進める方法である。状態 x から x への遷移率を k(xx) とし、マスター方程式

dP(x,t)dt=xx[P(x,t)k(xx)P(x,t)k(xx)]

で確率分布が時間発展すると考える。kMCは、この確率過程の標本過程(軌道)を生成し、平均挙動や分布を推定する。

具体的には、状態 x で可能なイベント i(原子ジャンプ、吸着・脱離、反応、空孔交換など)を列挙し、各イベントのレート ki を与える。総レート

R=iki

を定義すると、次に起こるイベントの種類は確率 ki/R で選ばれ、待ち時間 Δt は指数分布

P(Δt)=Rexp(RΔt)

に従う。したがって一様乱数 u(0,1) を用いて

Δt=lnuR

と生成すればよい。これがBortz–Kalos–Lebowitz(n-fold way)として知られる基本であり、イベント数が多い場合でも「実際に起こるイベントだけで時間を進める」ため、拡散のように長い時間スケールへ到達できる。

レート ki の物理モデルがkMCの本体である。しばしばArrhenius形

ki=νiexp(EikBT)

が用いられ、νi は試行頻度、Ei は活性化障壁である。障壁はNEB法などでDFTから求めることもあるが、大規模に列挙するのは難しいため、有効モデル(局所環境に応じた障壁推定)を組むことが多い。詳細釣り合いに相当する条件として、平衡へ戻るべき系では前後遷移の比がボルツマン因子に一致する

k(xx)k(xx)=exp[β(E(x)E(x))]

を満たすようにレートを設計する必要がある。これが満たされないと、長時間後に誤った定常状態へ落ちる。

材料科学のkMCでは、(i) どのイベント集合を採用するか、(ii) レートが局所環境でどう変わるか、(iii) 長距離弾性場や電場などの非局所効果をどう扱うか、が中心課題になる。kMCは「時間を持つMC」であるが、時間はレートモデルから生まれるため、モデルの妥当性がそのまま時間スケールの妥当性になる。

6.4 粒界、欠陥、成長過程

多結晶材料では、粒界の移動、粒成長、再結晶などが機械特性や輸送特性を大きく左右する。これらは幾何学(曲率、トポロジー)とエネルギー(界面エネルギー、欠陥相互作用)が強く結び付いた問題であり、MCは粗視化モデルとして有効である。代表はPotts模型を用いた粒成長モデルであり、格子点 i に方位ラベル σi{1,,q} を割り当て、

E({σ})=i,jJ(1δσi,σj)

とする。隣接点のラベルが異なるとエネルギーを持つため、粒界長(2次元)または粒界面積(3次元)に比例するエネルギーが得られる。ここでの q は方位状態数であり、現実の連続方位を離散化したものと解釈できる。

更新は局所的な再配向として実装される。ある格子点を選び、隣接点のラベルのいずれかへ変更を提案し、エネルギー差 ΔE に対してMetropolis受理確率

Pacc=min{1,exp(βΔE)}

で受理する。β は「熱ゆらぎ」そのものというより、界面移動の粗視化におけるノイズ強度と解釈されることが多い。温度を下げると界面が滑らかになり、温度を上げると粗い揺らぎが増える。ここでの重要点は、モデルが「曲率駆動の界面縮退」を再現するように設計されていることであり、粒界エネルギーが等方的なら平均粒径が時間の平方根に比例して成長する法則などが得られる。

欠陥集団(空孔クラスター、析出核、転位密度場の粗視化)も、確率変数の選び方次第でMCに落とし込める。例えば、格子上の占有変数 ni{0,1}(欠陥の有無)と相互作用エネルギーを定義すれば、欠陥の凝集や相分離が平衡MCで議論できる。非平衡過程(拡散・凝集の時間発展)を扱うならkMCへ移る。粒界や欠陥の問題では、エネルギーに加えて「幾何学的制約」が効くため、更新規則がその制約を尊重しているかが重要になる。

さらに、連続体場(例えば相場 ϕ(r))を離散化して、その場を確率変数としてMCを行う方法もある。これは、場の自由エネルギー F[ϕ] に対して

π(ϕ)exp[βF[ϕ]]

でサンプリングすることに相当し、界面粗視化や熱揺らぎを含む相場の分布を求めるのに使われる。フェーズフィールドと併用する場合は、動力学(時間発展)と統計(分布)の役割分担を明確にする必要がある。

6.5 量子統計のモンテカルロ

軽元素、低温、あるいは強い量子揺らぎが支配する材料では、古典的なボルツマン重みだけでは不十分であり、量子統計を扱う必要がある。量子系の分配関数は

Z=Tr[eβH^]

であり、状態和が波動関数ではなく演算子のトレースになっている点が古典と異なる。経路積分(path integral)に基づくモンテカルロは、この量子分配関数を「高次元の古典統計問題」へ写像してサンプリングする方法である。

基本はTrotter分解であり、H^=T^+V^ と分けて

eβ(T^+V^)=(eΔτT^eΔτV^)M+O(Δτ2),Δτ=βM

と近似する。M を大きくすると誤差が減り、M で厳密に近づく。これにより、粒子の量子揺らぎが「虚時間方向に連結した多数のビーズ(リングポリマー)」として表現され、古典的な多体問題としてMCサンプリングできる。例えば、1粒子の位置 rM 個のスライス r(1),,r(M) に拡張され、隣接スライス間にばね項が現れるという構造になる。

材料科学で出番があるのは、ゼロ点振動、同位体効果、水素拡散、軽元素の相安定などである。例えば、古典近似では振動の基底状態エネルギーがゼロと見なされるが、量子ではゼロ点エネルギー 12ω が残るため、相の自由エネルギー差が変わりうる。この差が小さい材料では、量子核効果を見積もれることが相安定性の議論に効く。

一方で、量子モンテカルロには困難もある。フェルミオンの符号問題により、一般の電子系では統計誤差が指数的に増大しやすい。材料科学で電子相関を量子MCで扱う場合は、模型ハミルトニアン(Hubbard模型など)や、符号問題が軽減される特殊条件に限定されることが多い。したがって、どの自由度(核か電子か)を量子的に扱い、どの近似(ボルン–オッペンハイマー、調和近似など)を採用しているかを明確にする必要がある。

6.6 他手法との連携

MCは平衡分布と揺らぎを扱うのが得意であり、DFTやMD、フェーズフィールド、有限要素法と役割が異なる。連携の要点は、(i) ミクロから有効ハミルトニアンを作る、(ii) MCで有限温度平均を取る、(iii) 必要ならさらに粗視化して連続体へ渡す、という段階である。ここでは、代表的な結び方を「何を入出力として受け渡すか」という観点で整理する。

第一に、DFTとMCの連携である。DFTは T=0 のエネルギー面や局所相互作用を与えるが、温度で平均する機能は持たない。クラスター展開により E(σ) を構成し、MCで ECV、秩序変数分布、短距離秩序などを評価すれば、有限温度相図や秩序化転移温度へ接続できる。ここでの要点は、DFTの系サイズ制約を「有効モデル」によって回避し、統計力学の計算へ移すことである。

第二に、MDとMCの補完である。MDは時間発展(運動方程式)を持つため動的現象に強いが、長時間平衡化が難しい場合がある。MCは時間の解釈を持たない代わりに分布の探索に集中できるため、平衡構造探索や相共存の評価に向く。例えば、MDで局所的な緩和や振動特性を調べ、MCで配置自由度(置換、欠陥、方位)を熱平衡化する、といった役割分担が成立する。kMCはこの中間に位置し、レートモデルが用意できるなら長時間の拡散・析出を扱える。

第三に、連続体手法への受け渡しである。フェーズフィールドは自由エネルギー汎関数に基づく時間発展を扱い、有限要素法は複雑形状の力学場や拡散場を解く。これらへMCから渡す量は、自由エネルギー密度、界面エネルギー、移動度、相図(化学ポテンシャルと濃度の関係)などの粗視化パラメータである。例えば、MCで求めた平衡組成や化学ポテンシャル曲線を、フェーズフィールドの自由エネルギー関数形にフィットし、相分離ダイナミクスへつなぐことができる。重要なのは、粗視化で失われる自由度(短距離秩序など)が、連続体パラメータにどの程度埋め込まれているかを理解することである。

最後に、全体としての確認点は「同じ物理量を二通りで評価して整合を見る」ことである。例えば、相境界をMCで求めたなら、別の自由エネルギー推定(熱力学積分、再重み付け)でも一致するかを確かめる。連携は便利であるが、受け渡し量の定義が曖昧だと誤差が増幅されるため、何を平均し、何を固定し、どの集団の量として定義したかを明示する姿勢が必要である。

小まとめ

第5章では、統計集団の分布をサンプリングするという立場から、詳細釣り合いに基づくMetropolis–Hastings、自己相関を含む誤差評価、臨界減速を抑えるクラスタ更新、障壁越えのための拡張アンサンブル、自由エネルギー差の推定法を数式で組み立てた。第6章では、格子模型と秩序変数、クラスター展開による第一原理エネルギーの圧縮とMCによる有限温度相図、kMCによる非平衡時間発展、粒界・成長の粗視化モデル、量子統計への拡張、そして他手法との役割分担を整理した。今後は、MCで得られる平衡分布や自由エネルギーを、フェーズフィールドや有限要素法の連続体記述へ整合的に埋め込み、ミクロ起源を保ったままメゾ・マクロの材料現象を予測する方向へ進むことになる。

第7章 分子動力学の基礎と数値積分

本章では、原子の運動を古典力学として記述し、数値積分によって時間発展を追跡する分子動力学(MD)の骨格を式から組み立てる。MDは乱数で分布を直接サンプリングする手法ではなく、運動方程式を解いた軌道から統計量を評価する手法であるため、力学と統計の接続を意識して学ぶ必要がある。

7.1 ニュートン方程式とハミルトニアン

古典MDの出発点は、各原子 i の座標 ri(t) と運動量 pi(t)=mivi(t) がニュートン方程式

mid2ridt2=Fi

に従うという仮定である。ここで力 Fi はポテンシャルエネルギー U({r}) から

Fi=riU({r})

として与えられることが多く、MDは結局「U を与えれば運動が決まる」という形式で閉じる。

保存量を見通しよく扱うために、ハミルトニアン

H({r},{p})=ipi22mi+U({r})

を導入する。ハミルトン方程式は

dridt=Hpi,dpidt=Hri

であり、これを成分ごとに書けばニュートン方程式と同値である。実際、

dridt=pimi=vi,dpidt=riU=Fi

となり、速度と力が自然に現れる。

ハミルトニアン形式の利点は、エネルギー保存の条件が明確になる点である。外場が時間に依存せず、数値積分が厳密解であれば

dHdt=i(Hridridt+Hpidpidt)=0

が成り立ち、全エネルギーは保存される。現実の計算では離散時間の近似誤差があるため、エネルギーは完全には一定にならないが、後述するVerlet系の積分法は「長時間での保存量の崩れ方」を抑えるように設計されている。

統計との接続としては、外部と熱交換も圧力交換もしない理想化では、系は微視的には一定エネルギー面上を運動し、微小正準集団(NVE)の記述に対応する。すなわち、軌道の時間平均が位相空間上の平均と一致する(エルゴード性が十分である)なら、時間平均で物性が評価できるという立場になる。したがってMDでは、運動方程式の解法だけでなく「どの熱力学条件を模倣しているか」を常に意識する必要がある。

7.2 Verlet系積分法

MDで支配的になる数値誤差は、時間の離散化によって生じる。Verlet積分は、テイラー展開から自然に導出でき、時間反転対称性と長時間安定性の点で基本となる。

まず、ある粒子の位置 r(t) を時刻 t の前後でテイラー展開する:

r(t+Δt)=r(t)+r˙(t)Δt+12r¨(t)Δt2+16r(3)(t)Δt3+O(Δt4),r(tΔt)=r(t)r˙(t)Δt+12r¨(t)Δt216r(3)(t)Δt3+O(Δt4).

これらを加えると奇数次の項が消え、

r(t+Δt)+r(tΔt)=2r(t)+r¨(t)Δt2+O(Δt4)

である。よって

r(t+Δt)=2r(t)r(tΔt)+r¨(t)Δt2+O(Δt4)

を得る。r¨=F/m を用いれば、提示された位置Verlet更新式に一致し、局所誤差が O(Δt4)、全体誤差が O(Δt2) の二次精度であることも見通せる。

ただし位置Verletだけでは速度が直接得られないため、実務では速度Verlet(velocity Verlet)が広く用いられる。導出は同様に、速度のテイラー展開

v(t+Δt)=v(t)+v˙(t)Δt+12v¨(t)Δt2+O(Δt3)

を使い、v˙(t)=a(t)=F(t)/m を代入して構成する。具体的には

r(t+Δt)=r(t)+v(t)Δt+12a(t)Δt2,

まずこの式で位置を更新して新しい力 F(t+Δt) を計算し、

v(t+Δt)=v(t)+12[a(t)+a(t+Δt)]Δt

で速度を更新する。力評価が1ステップに1回で済む点が計算効率上の利点である。

Verlet系が長時間安定になりやすい理由は、時間反転対称性と、離散写像が位相空間の幾何を保ちやすい(シンプレクティック性)ことに関係する。直感的には、エネルギーが単調に漂うのではなく、ある有効ハミルトニアンの周りで小さく揺らぐ形になりやすい。したがって、短時間の誤差が小さいだけでなく、長時間計算でのドリフトが抑えられる。

時間刻み Δt の選択は、最も速い運動モードで決まる。最も単純な調和振動子 mx¨=kx(角振動数 ω=k/m)に対して、数値的安定性はおおむね

Δt<2ω

のような条件で支配される。固体では原子振動の高周波成分(軽元素の伸縮、硬い結合)により ω が大きくなり、Δt を小さくする必要がある。したがって、温度や応力でゆっくり変わる現象を見たくても、時間刻みは高速振動に縛られるという構造がある。

7.3 温度制御と圧力制御

MDの軌道から温度や圧力を定義するには、まず運動エネルギーと自由度の関係を用いる。古典系で拘束がないとき、平衡での等分配則より

K=ipi22mi=f2kBT

が成り立つ。ここで f は自由度であり、3次元の N 粒子系で重心運動を除くなら f=3N3 がよく使われる。したがって瞬間温度は

T(t)=2K(t)fkB

で定義できるが、これは瞬間量であり揺らぐ。物性評価では時間平均やブロック平均で温度の統計を扱う必要がある。

温度一定の条件を模倣するには、熱浴を導入して正準集団(NVT)の分布

π({r},{p})exp[βH({r},{p})]

を実現するのが目標である。熱浴には複数の考え方があり、どれを選ぶかで動力学の性質が変わる点を理解する必要がある。

Nosé–Hoover法は、拡張変数を導入して決定論的に熱浴効果を表現する。代表的には摩擦係数に相当する変数 ξ を導入し、

r˙i=pimi,p˙i=Fiξpi,ξ˙=1Q(ipi2mifkBT)

のように運動方程式を拡張する(Q は熱浴の慣性に相当するパラメータである)。ipi2/mi が目標値 fkBT より大きいと ξ が増えて減速し、小さいと ξ が減って加速するため、運動エネルギーが調整される。Q の選び方が強すぎると温度が激しく振動し、弱すぎると温度制御が効きにくいという性質がある。

確率的な熱浴としてはLangevin法がある。これは

mir¨i=Fiγmir˙i+Ri(t)

の形で摩擦項とランダム力を加える。ランダム力は平均ゼロで

Ri(t)=0

かつ揺らぎ散逸関係により

Riα(t)Rjβ(t)=2γmikBTδijδαβδ(tt)

を満たすように与えられる。これにより長時間で正準分布に収束するが、摩擦係数 γ が大きいと動力学が過度に減衰し、拡散や振動スペクトルが変わりうる。したがって、構造探索なのか、輸送係数評価なのかで熱浴の選択が変わる。

圧力一定(NPT)のためには体積やセル形状を動的自由度として追加する。圧力の瞬間値は、運動量と相互作用から定まるビリアル式として

P=NkBTV+13ViriFi

(周期境界で原点の取り方に注意しつつ)と書ける。より一般には応力テンソル σ を導入し、成分は

σαβ=1V(imiviαviβ+i<jrij,αFij,β)

の形を持つ(第一項が運動項、第二項がビリアル項である)。等方圧は P=13Tr(σ) として得られる。

Parrinello–Rahman法はセル行列 h(格子ベクトルを列に持つ行列)を自由度として導入し、セル変形を含む運動方程式を与える。これにより、相転移や弾性変形を自然に表現できるが、セル自由度にも時間刻みと安定性が関わるため、熱浴・圧力浴・数値積分を一体として設計する必要がある。固体の弾性定数や相変態を扱う場合、等方体積だけを変えるバロスタットでは不十分になりうるため、目的に応じて等方制御と異方制御を使い分けることが重要である。

7.4 ポテンシャルエネルギー面

古典MDの予測能力は、ポテンシャルエネルギー U({r}) の妥当性にほぼ依存する。これは、運動方程式が

mir¨i=riU({r})

で閉じているからである。したがって「どの結合をどの程度再現したいか」に応じて U の形式を選び、どの条件で校正されているかを理解する必要がある。

経験力場は、結合伸縮・角度・ねじれ・非結合相互作用を足し合わせる形で

U=Ubond+Uangle+Utorsion+UvdW+UCoulomb

のように構成されることが多い。分子系ではこの分解が直感的であり、パラメータの意味が追いやすい。一方、金属では電子の局在結合というより集団的な結合が支配的であり、埋め込み原子法(EAM)のように

U=iF(ρi)+12ijϕ(rij),ρi=jif(rij)

の形で「局所電子密度」 ρi を媒介に多体性を表現するモデルが有効である。ここで F は埋め込みエネルギー、ϕ は対相互作用、f は密度寄与である。多体項を明示的に入れているため、金属の凝集エネルギーや弾性、欠陥形成などを比較的よく再現できる。

共有結合性が強い材料(Si、C、窒化物など)では、結合角や局所配位を反映するためにTersoff型やStillinger–Weber型のような角度依存項を含むモデルが使われる。これらは

U=i<jV2(rij)+i<j<kV3(rij,rik,θjik)

のように2体・3体で構成され、結晶とアモルファスの違いを「角度分布」として表現できる利点がある。反応過程を含む場合は、結合の切断・生成を表現する反応力場(例:ReaxFFのような結合次数に基づく形式)が用いられるが、計算負荷が大きく、パラメータも多いため、適用条件と再現範囲を慎重に見極める必要がある。

第一原理MDは、U を固定の解析式で与える代わりに、電子状態計算(多くはDFT)で各時刻の力を求める。Born–Oppenheimer MD(BOMD)では、各時間ステップで電子系を基底状態へ自己無撞着に収束させ、得られた全エネルギー EDFT({r}) をポテンシャルとして

Fi=riEDFT({r})

を用いる。Car–Parrinello MD(CPMD)では電子波動関数に仮想慣性を与えて拡張ラグランジアンを作り、電子と核を同時に時間発展させることで、毎ステップの厳密なSCFを避けるという考え方である。いずれも、電子の計算誤差が力へ直結するため、時間積分だけでなく電子収束の誤差も物理量へ反映されることを理解する必要がある。

ポテンシャルエネルギー面という言葉は、原子配置空間におけるエネルギー地形を意味する。平衡構造はその極小であり、遷移状態は鞍点である。MDはこの地形上を熱エネルギーで走る軌道であるため、観測したい現象(拡散、相転移、破壊)がどの程度の障壁を持つかを見積もると、必要な温度や時間、手法選択(加速法や別手法との接続)を論理的に決めやすい。

7.5 長距離相互作用

イオン性結晶や極性材料、荷電欠陥、界面の電気二重層などでは、クーロン相互作用

UCoul=14πε0i<jqiqjrij

が本質的であり、単純な距離カットオフは系統誤差を生む。クーロン相互作用は 1/r で減衰が遅く、周期境界条件では無限個の画像セルとの相互作用を含むため、取り扱いが数学的に繊細になる。

Ewald和は、長距離相互作用を実空間で急速に収束する項と逆空間(フーリエ空間)で収束する項に分けて計算する。概念は、点電荷をガウス分布で「ぼかした電荷」と「補正項」に分解することに対応する。エネルギーは形式的に

U=Ureal(α)+Urecip(α)+Uself(α)

と分解され、α は分割のパラメータである。Ureal は実空間での補正付き 1/r をカットオフ内で足し上げ、Urecip は逆空間格子ベクトル k の和として評価する。Uself は自己相互作用を打ち消す項である。α とカットオフと逆空間の打ち切りを連動させることで、所望の精度に対する計算量を制御できる。

粒子メッシュ法(PME、PPPM)は、Ewaldの逆空間和を格子上のフーリエ変換(FFT)で加速する考え方である。電荷をメッシュへ割り当て(補間)して密度場を作り、ポアソン方程式に相当する処理をFFTで解いて逆空間寄与を得る。計算量は概ね

O(NlogN)

へ改善し、大規模系で有利である。補間次数、メッシュ間隔、実空間カットオフ、α の選択が誤差と計算量を左右するため、電気的性質(誘電率や界面電場)を議論する場合は設定の影響が大きい。

周期境界と電荷の総和にも注意が必要である。総電荷がゼロでない系を周期境界で扱うと、数学的にエネルギーが発散するため、多くの実装では一様な中和背景を仮定して計算を成立させる。この仮定は、欠陥形成エネルギーや界面電位の解釈に影響しうるため、電荷状態を扱うときは「何を仮定して無限系を定義しているか」を理解する必要がある。同様に、双極子・四重極子の補正や、2次元周期(スラブ系)でのEwaldの派生形式など、境界条件に応じた取り扱いが存在する。

長距離相互作用はクーロンだけではない。磁気双極子相互作用、分散力(C6/r6)の長距離補正、弾性場の長距離性などもあり、何が支配的かでモデル化が変わる。MDでは計算を回せる形に落とすために近似を入れがちであるが、その近似が物性に与える影響を式のスケールから見積もる姿勢が重要である。

7.6 物性評価と統計

MDの出力は軌道 {ri(t),vi(t)} であり、物性はそこから統計的に評価する。構造の評価としては、動径分布関数

g(r)=14πr2ρNijδ(rrij)

が基本であり、液体やアモルファスの短距離秩序を定量化できる。散乱実験と対応する構造因子 S(q)

S(q)=1N|j=1Neiqrj|2

として与えられ、実空間と逆空間の両方から構造を理解できる。

輸送係数の多くは、時間相関関数から得られる。自己拡散係数は、アインシュタイン関係により平均二乗変位(MSD)を用いて

D=limt16t|ri(t)ri(0)|2

で評価される。t の極限が重要であり、短時間では弾道運動(t2)から拡散(t)へ移るため、十分長い時間窓で線形領域を取る必要がある。同じ D はグリーン・久保関係として速度自己相関関数からも

D=130vi(0)vi(t)dt

と書けるため、両者が一致するかを確かめると数値の信頼性が増す。

粘性係数は応力テンソルの非対角成分(せん断応力)Pαβαβ)の時間相関で

η=VkBT0Pαβ(0)Pαβ(t)dt

と与えられる。熱伝導率は熱流 JQ の自己相関として

κ=VkBT20JQ(0)JQ(t)dt

となる。ここで熱流の定義はポテンシャルの形式に依存し、多体ポテンシャルでは2体和の単純な形では書けないため、実装が与える熱流の定義を理解して使う必要がある。

有限サイズと有限時間は、MDの誤差構造を決める。周期境界下の拡散係数は、流体力学的な自己相互作用により系サイズ L に依存してずれ、代表的に

D(L)=D()ckBT6πηL

のような補正が現れる(c は境界条件に依存する定数である)。粘性や熱伝導も同様に長距離のモードに敏感であり、セルサイズを変えた比較や、複数手法(平衡相関と非平衡駆動)での照合が有効である。

統計誤差は、MDでも自己相関に支配される。時間相関関数の積分は、積分上限を有限に切ることによる系統誤差と、ノイズの積分による統計誤差が競合する。したがって、ブロック平均、複数独立初期条件、相関時間の見積もりを用い、誤差の減り方が想定通りかを確認しながら評価する必要がある。MDは軌道が連続であるため一見データ量が多いが、有効独立サンプル数は自己相関時間で決まる点を忘れてはならない。

第8章 分子動力学の材料応用と機械学習ポテンシャル

本章では、MDが材料のどの問題に強く、どの制約を持つかを、変形・輸送・反応という観点から整理し、近年の機械学習ポテンシャルによる拡張を理論と実装の接点で説明する。特に、第一原理の精度と大規模MDの計算量の間を埋める方法として、学習モデルと不確かさ評価が不可欠になりつつある。

8.1 変形、転位、破壊

塑性変形や破壊は、原子配列の不連続や欠陥の生成を伴うため、連続体の場だけでは直接見えにくい素過程が支配する。MDでは、外部ひずみを与えたときに転位がどこで核生成し、どの面で滑り、双晶がどう成長するかを原子配置として追跡できる。したがって、材料強度の起源を「欠陥の生成と運動」という機構へ分解して理解するのに向く。

計算では、ひずみテンソル ε(t) を時間的に変化させることで変形を与えることが多い。例えば単軸伸長ならセル長 Lz

Lz(t)=Lz(0)[1+ε˙t]

のように更新し、ひずみ速度 ε˙ を設定する。周期境界下でせん断を与える場合はLees–Edwards条件などの移動周期境界を用いる。応力は前章のビリアル応力から求まり、応力–ひずみ曲線を作れば弾性域、降伏、加工硬化・軟化、破断に至る過程が定量化できる。

転位の同定には、バーガースベクトルや不整合の指標が必要である。結晶中の転位は、結晶学的には閉曲線に沿った変位の不一致として

b=du

で定義される。MDの離散原子配置では連続変位場 u を直接は持たないため、近傍の局所格子を参照し、原子の局所配位や結晶構造識別量(例:CNAや局所対称性指標)から転位芯や積層欠陥を抽出する。ここで重要なのは、転位は「点欠陥の集合」ではなく「連続体の特異線」に相当するため、解析も幾何学的な見方を取り入れる必要がある点である。

MDの制約は時間スケールにある。時間刻みが Δt1 fs 程度であれば、107 ステップでも 10 ns 程度であり、室温で起こる拡散や時効、遅い破壊過程の多くは直接には届かない。したがって、観測したい素過程が「高速で起きるか」「障壁越えが支配して遅いか」を区別し、遅い場合は加速法や別手法との接続が必要になる。

加速MDの考え方は、遷移状態理論(TST)と結び付く。ある盆地からの脱出レートが

kνexp(ΔEkBT)

で与えられるとき、ΔE が数十 kBT に達すると待ち時間が急増する。Hyperdynamicsは盆地内部にバイアスを加えて滞在時間を短縮し、温度加速法は高温で遷移を起こしてレートの温度依存から低温へ外挿する。並列レプリカ法は同一状態から多数のレプリカを走らせ、最初の脱出イベントで時間を加速する。これらは「運動方程式をそのまま速く解く」のではなく、「待ち時間統計を保存しつつ稀なイベントを効率化する」方法であり、kMCとの接続と同じ思想である。

破壊については、原子結合の切断や空孔成長、き裂先端での原子再配列が中心となる。経験力場では化学結合の切断が表現できない場合があるため、反応力場や第一原理MD、あるいは機械学習ポテンシャルの導入が重要になる。破壊は局所現象である一方で応力場は長距離に及ぶため、境界条件とセルサイズの影響が大きく、連続体との対応付け(応力拡大係数やJ積分の概念)を意識すると議論が整理される。

8.2 熱輸送と非平衡MD

熱伝導は、電子寄与と格子寄与に分けて考えることが多いが、古典MDが直接扱えるのは主に格子(原子振動)による熱輸送である。熱伝導率 κ はフーリエの法則

q=κT

で定義され、熱流 q と温度勾配 T の比例係数として測られる。MDでは、平衡揺らぎから求める方法と、実際に温度勾配や熱流を与える方法の二通りがあり、両者の一致が信頼性の根拠になりうる。

平衡MD(EMD)では、前章のグリーン・久保関係

κ=VkBT20JQ(0)JQ(t)dt

により評価する。これは線形応答理論の結果であり、外部駆動を与えずに内部の揺らぎだけから輸送係数が得られる点が美しい。ただし、積分の収束は遅く、長時間尾部や有限サイズの長波長フォノンに敏感である。特に高熱伝導材料では自己相関の減衰が遅く、統計誤差の扱いが難しくなるため、ブロック平均と積分上限の選び方が重要である。

非平衡MD(NEMD)では、系の一端を高温、他端を低温に保つなどして定常温度勾配を作り、定常熱流から κ を求める。例えば長さ L の系で z 方向に温度差 ΔT を与え、定常熱流 qz が得られれば

κqzdT/dzqzLΔT

である。温度は局所運動エネルギーから定義されるが、局所平衡が成立する程度の空間分割が必要である。熱浴領域の厚さが薄すぎると境界で非物理な散乱が起き、厚すぎると有効輸送領域が短くなるため、設計は物理と数値の両面に依存する。

Müller–Plathe法は、熱流を直接制御する逆NEMDの一種である。系の高温側と低温側で原子の速度を交換し、外から与えた運動エネルギー交換量により熱流を既知にする。交換によって生じた温度勾配を測り、

κ=JT

で求める。熱流 J が制御できるため統計的に安定しやすいが、交換が系の動力学へ与える影響を小さくするには交換頻度を適切に選ぶ必要がある。

界面熱抵抗(Kapitza抵抗)は、異種材料界面での温度不連続として観測され、

RK=ΔTinterfaceq

で定義される。MDでは、界面を含む系に熱流を流して ΔT を測定し、界面の散乱機構(音響ミスマッチ、拡散ミスマッチ、界面粗さ、欠陥)を原子スケールで解析できる。ナノ構造では界面が支配的になるため、バルクの κ だけでは記述できず、界面抵抗と有効媒体の考え方が必要になる。

熱輸送では、サイズ効果が物理そのものとして現れる。フォノンの平均自由行程が系サイズに近いと、散乱が境界で支配され、バルク熱伝導率とは異なる有効値が得られる。したがって、MDで得た κ(L) をそのままバルクとみなすのではなく、フォノンスケールとセルスケールの関係を見積もり、必要ならサイズ外挿や別理論(ボルツマン輸送方程式)との接続を考える必要がある。

8.3 第一原理MDと反応過程

化学反応、電荷移動、結合の組み替えが本質である材料現象では、固定の古典ポテンシャルでは不十分になり、第一原理MDが必要になる。第一原理MDでは、原子配置 {r} に対してDFT等で全エネルギー E({r}) と力 Fi を求め、核の運動を時間発展させる。BOMDでは各ステップで電子を基底状態へ収束させるため、核の運動は断熱ポテンシャル面上の運動として解釈できる。

反応過程や拡散の議論では、自由エネルギーが中心となる。反応座標(集団変数)s({r}) を導入し、s の分布 P(s) から

F(s)=kBTlnP(s)+const.

と自由エネルギー曲線(平均力ポテンシャル)を定義できる。ここで P(s) は正準分布に対して

P(s)=1Zd{r}d{p} δ(ss({r}))exp[βH]

で与えられ、s の希な領域(遷移状態付近)のサンプルが少ないことが問題になる。

傘サンプリングは、バイアス Ub(s) を加えて

πb({r},{p})exp[β(H+Ub(s))]

でサンプリングし、後で再重み付けして P(s) を復元する。例えば調和バイアス

Ub(s)=k2(ss0)2

を用いれば、狙った s0 近傍を集中的にサンプルできる。複数の窓を s 軸に沿って配置し、各窓の分布を統合して滑らかな F(s) を得る。

メタダイナミクスは、探索済み領域へ時間とともに反発バイアスを積み上げ、自由エネルギー障壁を乗り越えさせる方法である。集団変数 s(t) の履歴に対して、バイアスを

V(s,t)=t<twexp[(ss(t))22σ2]

のようにガウス関数で加え、系が同じ領域に留まることを避ける。十分長時間後に、適切な条件では V(s,t)F(s) を近似するという考え方により、自由エネルギーが得られる。実際には、集団変数の選び方が支配的であり、遷移の本質自由度を外すと、バイアスがかかっても障壁を越えられないか、誤った経路を強調してしまう。

第一原理MDは計算量が大きく、時間スケールの制約がさらに厳しい。したがって、反応・拡散の障壁を議論する場合は、(i) 最小エネルギー経路(NEB等)で鞍点を求める、(ii) 調和遷移状態理論でレートを推定する、(iii) 必要に応じて自由エネルギーを傘サンプリング等で評価する、というように、目的に応じて複数の方法を組み合わせることが合理的である。MDは万能ではないが、適切な自由エネルギー記述へ接続できれば、温度依存の反応性や相安定を一貫した枠組みで議論できる。

8.4 機械学習ポテンシャル

機械学習ポテンシャルは、第一原理に近い精度と古典MDに近い速度を両立することを目標に、原子配置からエネルギーと力を予測するモデルである。基本は教師あり学習であり、DFTなどで得たデータ集合 {({r}(n),E(n),{F}(n))} を用いて、モデル Eθ({r}) を学習する。力はエネルギーの勾配として

Fi,θ=riEθ({r})

で定義できるため、エネルギーを学習しておけば自動微分により力が得られるという構造がある。これにより、エネルギー保存則と整合した力が得られ、NVE計算での挙動が安定しやすい。

多原子系のエネルギーを直接高次元関数として学習するのは難しいため、局所性を仮定して

Eθ({r})=iεθ(Ei)

のように原子環境 Ei の和として表すことが多い。Ei はカットオフ半径内の近傍配置であり、これにより計算量が O(N) に近い形でスケールする。どの情報を Ei に入れ、どの不変性・共変性を満たすかがモデルの性能を左右する。

物理的に必須なのは対称性である。エネルギーは並進・回転・同種原子の置換に対して不変でなければならない。すなわち、任意の回転 RSO(3) に対して

E({Rri})=E({ri})

が成り立つ必要がある。一方、力はベクトルであるため回転に対して共変であり、

Fi({Rr})=RFi({r})

を満たす必要がある。E(3)等変性(並進・回転・反転を含むユークリッド群の等変性)を組み込んだグラフニューラルネット型ポテンシャルは、この要請をモデル構造として実現する。メッセージパッシングにより原子間の相互作用を伝播させ、スカラー(エネルギー)とベクトル(力)の変換則を保つように設計されるため、少量データでも汎化しやすい傾向がある。

一方で、等変性モデルだけが選択肢ではない。対称関数(Behler–Parrinello型)、SOAP記述子に基づくガウス過程回帰(GAP)、多項式近似(MTP)、ビスペクトラム(SNAP)など、記述子と回帰器の組み合わせは多様である。重要なのは、(i) 局所環境の表現が連続で微分可能であること、(ii) 対称性が満たされること、(iii) 学習データの範囲で誤差が小さいだけでなく、MDで必要になる力の精度と滑らかさが確保されることである。

学習の損失関数は、エネルギーと力を同時に合わせる形が多い:

L(θ)=wEn|Eθ(n)E(n)|2+wFni|Fi,θ(n)Fi(n)|2+

ここで重み wE,wF は、エネルギーと力の次元・重要度のバランスを取るために導入される。材料の相安定や欠陥形成エネルギーを重視するならエネルギー誤差が重要になり、MDでの動力学や輸送を重視するなら力誤差とその滑らかさが重要になる。目的物性から逆算して学習目標を設計することが必要である。

8.5 外挿と不確かさ

機械学習ポテンシャルの最大の弱点は、訓練データの外側での誤りである。学習は経験的リスク最小化であり、未知領域での誤差はデータが保証しない。MDは時間発展により状態空間を歩くため、初期状態が訓練領域にあっても、熱揺らぎや変形で外側へ出る可能性がある。このとき、わずかな力の誤りが軌道を大きく変え、最終的に非物理な構造へ進むことがある。

不確かさを定量化する代表はアンサンブルである。異なる初期値や異なるデータサブセットで学習した複数モデル θ(m) を用意し、予測の分散

Var[E({r})]1Mm(Eθ(m)E¯)2,E¯=1MmEθ(m)

を不確かさの指標として用いる。力についても成分ごとに分散を計算でき、分散が大きい局所環境は訓練不足である可能性が高い。ガウス過程に基づくモデルでは予測分散が理論的に得られるため、同様の判断に使える。

不確かさの指標は、学習ループではなくMD運用の安全装置として機能する。例えば、NVEでの全エネルギーが短時間で崩れる、局所温度が異常に上昇する、原子間距離が非物理な値へ入るといった挙動は、外挿による破綻の兆候になりうる。これらは不確かさ推定と併用すると強く、分散が大きい状態に入った瞬間に計算を止め、追加の第一原理計算でデータを補強するという運用が可能になる。

外挿検知のもう一つの考え方は、局所環境の距離を定義し、訓練集合との類似度を測ることである。記述子空間で最近傍距離を測り、ある閾値を超える状態を未知領域とみなす。これはモデル非依存に近いが、記述子の選び方で結果が変わるため、物理的に意味のある表現(回転不変で連続な記述子)を用いる必要がある。

不確かさを含む設計は、材料応用に直結する。例えば破壊や化学反応は、平衡近傍とは異なる高エネルギー・低配位の状態が現れやすい。そこは訓練データが不足しがちであり、外挿が起きやすい領域である。したがって、目的が破壊であれば、き裂先端の低配位状態や引張下の遷移状態を意図的にデータへ含める必要がある。同様に高温液体やアモルファスを扱うなら、広い温度・密度のサンプルが必要になる。このように、データ設計と外挿検知は一体であり、学習精度の数字だけでは安全性は判断できない。

8.6 マルチスケール接続

MDで得られる量を連続体モデルへ渡すと、実験スケールのプロセスやデバイス条件に接続できる。重要なのは、渡す物理量が連続体方程式の中でどの係数として現れるかを明確にすることである。例えば拡散であれば、連続体では濃度場 c(r,t)

ct=(Dc)

に従うとモデル化され、MDから得た D(T,c) を与えれば時間発展の速度が定まる。ここで D は温度や組成に依存しうるため、MD側でその依存性を系統的に評価することが接続の質を決める。

界面現象では、界面エネルギー γ や界面移動度 M がフェーズフィールドや曲率流の係数として現れる。例えば界面の曲率駆動移動が

vn=Mγκ

vn は法線速度、κ は曲率)で書けるとき、γM をMDや第一原理から推定できれば、粒成長や相変態の速度論をより根拠ある形で扱える。M は原子ジャンプの素過程と結び付くため、kMCやTSTとの組み合わせで推定する方が自然な場合もある。

弾性は有限要素法と接続しやすい。MDから弾性定数 Cijkl を得るには、微小ひずみ εkl を与えたときの応力応答

σij=Cijklεkl

を測るか、平衡揺らぎから弾性定数を推定する方法を用いる。結晶の異方性弾性を得られれば、FEM側の構成式へそのまま埋め込める。欠陥周りの局所弾性や界面の弾性不連続など、連続体では与えにくい情報をMDが補う形になる。

熱輸送では、MDで得た κ や界面熱抵抗 RK を用いて、熱拡散方程式

ρcpTt=(κT)+Q

の係数として組み込める。ナノ複合材料では、バルク κ だけでは不十分であり、界面抵抗が支配的になるため、界面条件として

q=ΔTRK

を与える形の接続が現実的である。ここでも、単位系(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)ϕ(r,t) を導入し、例えば

  • ϕ=0 が相A(母相)、
  • ϕ=1 が相B(析出相・固相など) のように、相を連続的に表す。界面は ϕ が 0 と 1 の間を滑らかに遷移する領域として表現され、界面位置を明示的に追跡しない。

秩序変数と組成場 c(r,t) を同時に扱う典型として、自由エネルギー汎関数を

F[ϕ,c]=Ω[f(ϕ,c)+κϕ2|ϕ|2+κc2|c|2]dr

と置く。ここで

  • f(ϕ,c):局所自由エネルギー密度(バルク項)
  • κϕ2|ϕ|2ϕ の勾配エネルギー(界面を有限幅にする)
  • κc2|c|2:濃度勾配に対するペナルティ(Cahn–Hilliard型で現れる) である。

なぜ勾配項が界面エネルギーを与えるのかを、1次元で確認する。c は一旦固定し、ϕ のみを考える。f(ϕ) を二重井戸(double-well)型

f(ϕ)=Wϕ2(1ϕ)2

とする(W>0 は障壁の高さ)。x 方向の界面(ϕ()=0, ϕ(+)=1)で

F[ϕ]=[f(ϕ)+κϕ2(dϕdx)2]dx

を最小化すると、オイラー=ラグランジュ方程式

δFδϕ=fϕκϕd2ϕdx2=0

を得る。両辺に dϕ/dx を掛けて積分すると

f(ϕ)κϕ2(dϕdx)2=const.

境界で dϕ/dx0 かつ f(ϕ)0 とすれば定数は0になり、

κϕ2(dϕdx)2=f(ϕ)

が得られる。これは「界面では勾配エネルギーと局所エネルギーが釣り合う」ことを意味する。

界面エネルギー(単位面積あたり)は、1次元の積分として

γ=[f(ϕ)+κϕ2(dϕdx)2]dx=κϕ(dϕdx)2dx

(上の釣り合い式を使った)と表せる。したがって κϕW により界面エネルギーと界面幅が決まる。例えば二重井戸の場合、界面厚さの代表スケールは

κϕW

のように見積もれる(厳密な係数はポテンシャル形で変わる)。ここがPFの重要点で、κϕW は単なる数値パラメータではなく、界面物性(γ と幅)へ対応する。

さらに、複数相や多相のときは秩序変数が複数になり、局所自由エネルギーは相ごとの自由エネルギーを補間する形で

f(ϕ,c)=h(ϕ)fB(c)+[1h(ϕ)]fA(c)+Wϕ2(1ϕ)2

のように設計することが多い。h(ϕ)h(0)=0, h(1)=1 を満たす補間関数であり、界面での熱力学一貫性(化学ポテンシャルの連続性など)に影響するため、形を選ぶ理由を理解して使う必要がある。

9.2 Cahn–Hilliard方程式

組成場 c(r,t) は「保存量」である。つまり領域 Ω 全体での総量

Ωc(r,t)dr

は(外部から出入りがない限り)一定である。保存則の一般形は

ct+J=0

であり、J は組成の流束である。拡散の基本仮定は「流束は化学ポテンシャルの勾配で駆動される」という形で

J=Mμ

と置く(M は移動度、μ は化学ポテンシャル)。これを保存則へ代入すると

ct=(Mμ)

となる。ここまでは「保存則+線形応答」の議論であり、どんな自由エネルギーでも成り立つ枠組みである。

PFの要点は、μ を自由エネルギー汎関数の汎関数微分として定義する点にある:

μ(r)=δFδc(r).

例えば

F[c]=(f(c)+κc2|c|2)dr

なら、変分計算により

δFδc=fcκc2c

を得る(勾配項の変分でラプラシアンが出る)。したがってCahn–Hilliard(CH)方程式は

ct=[M(fcκc2c)]

となり、4階のPDEになる。4階になるのは 2c をさらに で微分するためで、これが数値的な難しさ(剛性、境界条件)につながる。

CH方程式が相分離(スピノーダル分解)を表現する理由を、線形安定性で見る。均一状態 c=c0 に微小揺らぎ δc を入れて

c(r,t)=c0+δc(r,t)

とする。M を定数、fc0 周りで2次まで展開すると

fcfc|c0+2fc2|c0δc.

定数項は勾配で消えるので、フーリエモード δcexp(ikr+ωt) を代入すると

ω(k)=Mk2[f(c0)+κck2]

が得られる。ここで f(c0)<0(自由エネルギーが凹、スピノーダル領域)なら、小さな k に対して ω>0 となり揺らぎが成長する。さらに最も成長する波数は

kmax=f(c0)2κc

で与えられ、初期の相分離パターンの代表長さ 2π/kmax を決める。PFは「パターン長さが自由エネルギーの曲率と勾配係数で決まる」ことを式で示せる点が強い。

移動度 M は、拡散係数 D と関係して

DM2fc2

のように解釈できる(厳密には熱力学因子を含む)。したがって、M を定数にするか、濃度依存にするか、相依存にするかは「拡散の物理をどこまで入れるか」に対応する。

9.3 Allen–Cahn方程式

秩序変数 ϕ は一般に「保存されない」。相が入れ替われば局所の ϕ は変わるが、ϕdr が一定である必要はない。非保存場の最も基本的な緩和は「自由エネルギーの勾配降下(gradient descent)」として

ϕt=LδFδϕ

と定義する。これがAllen–Cahn(AC)方程式であり、L>0 は緩和係数(動力学係数)である。

ACの重要な性質は「エネルギー減少性」である。F[ϕ] の時間微分は

dFdt=δFδϕϕtdr=δFδϕ(LδFδϕ)dr=L(δFδϕ)2dr0

となり、F は単調減少する。つまり、ACは「物理的に自然な緩和過程」を保証する形になっている。数値解法でもこの性質が保たれていると、発散しにくく、結果が解釈しやすい。

δF/δϕ の具体形は、例えば

F[ϕ]=(Wϕ2(1ϕ)2+κϕ2|ϕ|2)dr

なら

δFδϕ=ϕ(Wϕ2(1ϕ)2)κϕ2ϕ=2Wϕ(1ϕ)(12ϕ)κϕ2ϕ

である。従って

ϕt=L[2Wϕ(1ϕ)(12ϕ)κϕ2ϕ]

となる。界面(ϕ0)ではラプラシアン項が効き、界面の曲率駆動(平均曲率に比例した移動)へつながる。よく知られた結果として、等方界面で外場がない場合、界面の法線速度 vn

vnMgbγκ

κ は曲率)に比例する形で現れる。ここで Mgb は粒界移動度で、PFの L と界面幅・勾配係数の組み合わせから有効的に決まる。つまり「ACの係数 L は界面の運動学(速度尺度)を規定する」という対応がある。

注意点として、ACは非保存なので、保存量を表すべき変数(濃度など)にそのまま使うと物理が破綻する。濃度はCH、相指標はACという役割分担が基本である。

9.4 連成場と駆動力

PFは自由エネルギーに項を追加することで連成を統一的に扱える。例えば弾性場を加えると

F[ϕ,c,u]=(f(ϕ,c)+κϕ2|ϕ|2+κc2|c|2+fel(ϕ,c,ε))dr

のようになる。u は変位場、ひずみは

ε=12(u+(u)T)

で定義される。弾性エネルギー密度は線形弾性なら

fel=12(εε0(ϕ,c)):C:(εε0(ϕ,c))

と書ける。ε0 は固有ひずみ(整合ひずみ、変態ひずみ)であり、相や濃度に依存させることで「析出に伴う格子不整合」や「マルテンサイト変態」などの駆動力を表現できる。

弾性場は準静的とみなすことが多く、力のつり合い

σ=0,σ=felε

を満たす u を各時刻で解く。すると ϕc の進化方程式には fel の変分が加わり、

ϕt=LδFδϕ,ct=(MδFδc)

の右辺に弾性由来の駆動力が現れる。重要なのは、弾性場が「全域に瞬時に伝わる長距離相互作用」である点で、数値的には連成が硬くなり、単純な陽解法が不安定になりやすい。

電場・磁場も同様に、例えば電気エネルギー項として

felec=P(ϕ,c)E+ϵ2|E|2

のような形を入れたり、磁性なら

fmag=μ0M(ϕ)H+fanis(ϕ,m)+fex(m)

のように追加できる。PFの方程式は「追加したエネルギーがそのまま駆動力になる」ため、連成を設計する際はエネルギー項の次元と意味を先に整えることが重要である。

9.5 異方性と結晶方位

界面エネルギーの異方性が強い材料では、界面形状の選択(特定面方位が出やすい)や樹枝状成長の枝の向きが決まる。PFで異方性を入れる方法は大きく2系統ある。

(1) 勾配エネルギー係数を方位依存にする: 界面法線 n=ϕ/|ϕ| を用いて

Fgrad=12κϕ(n)|ϕ|2dr

とする。例えば結晶対称性に合わせて

κϕ(n)=κ0[1+ϵag(n)]

のように置き、g(n) を立方晶なら nx4+ny4+nz4 のような形で作る。ϵa が異方性の強さを表す。

(2) 界面移動度や界面エネルギーを、界面項で直接異方化する: いずれにせよ、異方性が入ると界面は面方位を選びやすくなり、数値計算では「メッシュが勝手に方位を作る」偽異方性(grid anisotropy)が問題になる。偽異方性を避けるには、

  • メッシュを十分細かくする(界面幅に対して十分な点数)
  • 等方性の良い離散化(高次要素、等方差分、スペクトル)
  • 界面法線評価の精度を上げる(スムージング、正規化の安定化) などが必要になる。異方性の議論は物理だけでなく数値解析の要素が強いため、「得られた形が物理異方性か、メッシュ由来か」を必ず検証する(格子回転、格子幅変更、手法変更)姿勢が重要である。

9.6 熱ゆらぎと核生成

相変態の初期(核生成)は本質的に確率過程である。決定論的PF(CH/AC)だけでは、核生成が起きるきっかけがなく、初期ゆらぎを人工的に与えないと変態が始まらない場合がある。核生成を扱うためには、熱ゆらぎ項(ノイズ)を加える確率偏微分方程式(stochastic PF)を考える。

保存場(CH)では揺らぎも保存則を満たす必要があるため、

ct=(Mμ+η(r,t))

のように流束にノイズを入れる形になる。η は平均ゼロで、揺らぎ散逸関係により

ηα(r,t)ηβ(r,t)=2kBTMδαβδ(rr)δ(tt)

のような相関を持つように定める(白色雑音の理想化)。これにより、平衡分布がボルツマン分布と整合する。

非保存場(AC)なら

ϕt=LδFδϕ+ξ(r,t),ξ(r,t)ξ(r,t)=2kBTLδ(rr)δ(tt)

のように直接ノイズを加える形が基本である。

ただし実務では、連続体ノイズをそのまま離散化するとメッシュ依存性が出やすく、核生成率の定量は難しい。PFで核生成を議論する場合は、(i) 古典核生成理論(CNT)との整合、(ii) 反応座標・自由エネルギー障壁の見積もり、(iii) MD/MCと併用して核生成率を別手法で検証、などを組み合わせるのが安全である。PFは「成長と形態選択」には強いが、「稀な核生成イベントの絶対レート」は別枠で扱う必要が出やすい、という役割分担を理解する。

第10章 フェーズフィールドによる組織形成

本章では、第9章で導いた方程式(CH/ACと連成)を用いて、材料組織がどのように形成されるかを、凝固・粒成長・析出・反応拡散などの典型問題に沿って整理する。PFの重要点は「組織を決める因子(界面エネルギー、拡散、弾性、潜熱、異方性)」が、自由エネルギーと移動度として方程式中に明示的に現れることである。つまり、パラメータを変えると何が変わるかを論理的に追える。

10.1 凝固と結晶成長

凝固は、相(固相/液相)の秩序変数 ϕ、溶質濃度 c、温度 T が連成する典型問題である。最小モデルでも

  • ϕ:AC型(非保存)
  • c:CH型または拡散方程式(保存)
  • T:熱拡散方程式(保存) が結び付く。

温度場は、潜熱を含めると

ρcpTt=(kT)+Qlatent

の形で、潜熱項が ϕ の変化に比例して現れることが多い:

Qlatent=Lhϕt.

ここで Lh は潜熱密度に対応する係数である。固化が進むと潜熱が放出され、温度が上がり、固化駆動力が弱まるため、熱輸送と界面運動がフィードバックする。

溶質拡散も重要で、固相と液相で溶質の平衡濃度が異なる(分配係数)ため、界面で偏析が起きる。PFでは局所自由エネルギーを相ごとに用意し、

f(ϕ,c,T)=h(ϕ)fs(c,T)+[1h(ϕ)]fl(c,T)+Wϕ2(1ϕ)2

のように補間して、相平衡(化学ポテンシャル一致)を連続的に表現する。すると成長前線での溶質境界層、成長速度と拡散の競合が自動的に出る。

界面形態(セル状・樹枝状)が出る理由は、界面のわずかなゆらぎが、熱/溶質場の拡散と競合して成長しうるからである。簡単には、界面前方に溶質が溜まると局所的に凝固が抑制され、界面が不安定になる(組成過冷却に相当する)。PFはこれを場の方程式として解くため、形態選択を直接シミュレーションできる。

また、界面エネルギー異方性があると、枝の成長方向が結晶方位で固定される。異方性が弱いと枝は丸まり、強いと特定方位のファセットが現れる。ここで重要なのは、異方性の強さと界面幅・メッシュが絡むため、数値条件を変えて「枝の向きが物理的に安定か」を必ず検証することである。

10.2 粒成長と再結晶

粒成長は、粒界エネルギーによる曲率駆動で進む。最も基本的な理解は、粒界の法線速度が

vn=Mgbγgbκ

に比例するという関係である(κ は曲率)。曲率が大きい(小さな粒)ほど縮みやすく、結果として平均粒径が増大する。

PFでは、多結晶を表すために複数の秩序変数 {ηα} を使い、各粒 αηα1、他が0となるように設計することが多い。局所エネルギーとして

f({η})=αWηα2(1ηα)2+α<βWαβηα2ηβ2

のようにして、複数相(複数粒)が重なりにくいようにする。勾配項を入れれば粒界が有限幅で表現され、AC型の進化

ηαt=LαδFδηα

により粒界が移動する。

再結晶では、単なる曲率駆動だけでなく、加工で蓄積したひずみエネルギー(転位密度)を減らす駆動力が中心となる。そこで内部状態変数として転位密度 ρd(r,t) や蓄積エネルギー Estore(ρd) を導入し、

FF+Estore(ρd)dr

のように組み込む。さらに ρd 自体も回復や再配列で変化するため、反応拡散型の進化式(生成・消滅を含む)

ρdt=krecρd+

を併用することがある。再結晶PFはモデル自由度が増える分、パラメータ同定が難しくなるため、実験(EBSD、硬さ、XRD線幅)と対応する観測量を先に定義し、逆問題として組み立てるのが重要である。

10.3 析出、相分離、スピノーダル

析出と相分離はCH方程式が中心となるが、材料では弾性場の影響が非常に大きい。整合析出では、母相と析出相の格子定数差が整合ひずみとして蓄積し、析出物形状(球状・板状・針状)や配向を決める。これは「界面エネルギーを下げたい」という傾向と「弾性エネルギーを下げたい」という傾向の競合として理解できる。

弾性連成を入れた自由エネルギー

F[c,u]=(f(c)+κc2|c|2+fel(c,ε))dr

に対して

μ=δFδc=fcκc2c+felc

となり、CH方程式は

ct=(Mμ)

で進化する。弾性項 fel/c が「濃度が変わると弾性がどう変わるか」を通じて駆動力となり、弾性の長距離性により析出同士の相互作用も生じる。

弾性場の解法としては、(i) FEMで σ=0 を直接解く、(ii) 周期境界でGreen関数法(Khachaturyanのマイクロエラスティシティ)を用いてフーリエ空間で高速に解く、の2系統が有名である。周期系で均一な弾性定数を仮定する場合、Green関数法はFFTと相性が良い。一方、異方弾性や不均一弾性、複雑形状ではFEM連成が自然になる。

スピノーダル分解の初期波長が κcf で決まることは第9章で示したが、弾性が入ると有効的な曲率が変わり、成長するモードが選択される。つまり、弾性は「相分離のパターンの方向性とスケール」を変えうる。この点がPFの面白いところで、形態をエネルギーの項の競合として説明できる。

10.4 反応拡散と多相系

化学反応や相変態が拡散と同時に進む場合、反応拡散系としてモデル化する。例えば濃度 c が拡散しつつ反応で変化するなら

ct=(Dc)+R(c,ϕ,T)

のように反応項 R が入る。PFではさらに自由エネルギー駆動の形を保つために、反応を化学ポテンシャル差で駆動するように設計することがある。例えば二種の種 c1,c2 が反応で入れ替わるなら

R(μ1μ2)

のように置き、μi=δF/δci として熱力学整合性(詳細釣り合い)を保つ方向に設計する。

多相系(3相以上)では秩序変数が増え、しばしば「総和が1」の制約が必要になる。例えば相分率 ϕα を使い

α=1Nϕα=1,0ϕα1

を保ちたい。これを実現する方法としては、

  • ラグランジュ未定乗数で制約を入れる
  • 罰則項(ペナルティ)を自由エネルギーに追加する
  • 独立変数を N1 個に減らして最後を 1 で定義する などがある。どの方法でも、数値的に制約が崩れると「存在しない相が負になる」などの非物理が出るため、離散化・時間積分の段階で制約を保つ工夫が必要になる。

さらに多相では三重点(トリプルジャンクション)の角度条件が重要になる。等方界面エネルギーならYoungの関係により、三重点で界面張力が釣り合う角度が決まる。PFでは界面エネルギーを正しく入れると、三重点の角度が自動的に再現されるが、界面幅が大きすぎると三重点の力学がぼやけるため、界面幅とメッシュ分解能の関係を確認する必要がある。

10.5 パラメータ同定とデータ同化

PFモデルは「方程式は美しいが、係数が不明」という逆問題を必ず伴う。代表的なパラメータは

  • 界面エネルギー γκϕ,W に対応)
  • 拡散係数 D または移動度 M
  • 界面移動度(ACの L に対応)
  • 弾性定数 C、固有ひずみ ε0
  • 潜熱、熱伝導率 k などである。

まず「観測量」と「モデル量」の対応を明確にする。例えば界面エネルギーは、1次元平衡界面から

γ=κϕ(dϕdx)2dx

で与えられるので、界面幅 と合わせて κϕ,W を決められる。拡散係数はMDや実験から得た D(T) を用いて、熱力学因子を含めた関係で M を推定する。界面移動度は、単純な曲率駆動のテスト問題(円形粒が縮む速度など)で、PFの L を合わせることができる。

逆問題としては、観測データ yobs とモデル出力 y(θ)θ はパラメータ)のずれを最小化する

minθ J(θ)=12y(θ)yobs2+λR(θ)

の形に落とし込むのが基本である。R(θ) は正則化で、未同定(情報不足)を抑える。PFは計算が重いので、全パラメータを同時にフィットするより、

  • まず独立に決められるもの(弾性定数、拡散)を外部から固定
  • 次に界面物性(γ,)を合わせる
  • 最後に動力学係数(L,M)を時系列データで合わせる という段階的同定が現実的である。

データ同化の観点では、PFの状態(ϕ,c)自体を観測から更新する(カルマンフィルタ的)ことも可能であり、時間発展の予測と観測の融合ができる。ただし、観測ノイズとモデル誤差を分離しないと、パラメータを誤って吸収してしまうため、何を推定し、何を固定するかの設計が重要である。

10.6 ソフトウェア基盤

PFはPDEであるため、実装では

  • 空間離散(FEM/FVM/差分/スペクトル)
  • 時間積分(陰解法、半陰解法、分割法)
  • 疎行列ソルバ(Krylov、マルチグリッド、前処理) が不可欠になる。教育段階では差分やFFTのコードで全体像を掴み、研究段階ではFEM基盤や大規模ソルバを活用して連成・複雑形状へ拡張するのが自然な流れである。

MOOSEのPhase Fieldモジュールのような基盤は、弱形式・非線形ソルバ・自動微分・並列ソルバが統合されており、「方程式をどう書くか」が実装に直結している。重要なのは、ソフトの使い方そのものよりも、

  • 連立方程式がどの弱形式に落ちているか
  • どの項が陰的に、どの項が陽的に扱われているか
  • 収束判定(非線形残差、線形残差)が物理的に十分か を読めるようになることである。PFは見た目の結果がそれらしく出やすい反面、ソルバ未収束や偽拡散で「もっともらしい誤り」になりやすい。したがって、実装と方程式の対応を常に確認する姿勢が、研究での信頼性を支える。

章末の要点(理解のチェック用)

  • PFの核心は、自由エネルギー汎関数 F を立て、δF/δ() を駆動力として進化方程式を作る点にある。
  • CHは保存則から導かれ、tc=(Mμ)μ=δF/δc により4階PDEとなる。
  • ACは勾配降下であり、tϕ=LδF/δϕdF/dt0 を保証する。
  • 弾性・熱・電磁場などは F に項を追加することで統一的に連成できるが、数値的剛性が増す。
  • 異方性は物理異方性とメッシュ偽異方性を区別し、格子回転やメッシュ変更で検証する必要がある。
  • 核生成は確率過程であり、ノイズPFや他手法との併用で「どこまでをPFで扱うか」を整理する。

第11章 有限要素法の理論

本章では、有限要素法(Finite Element Method; FEM)が、偏微分方程式を弱形式へ移し、近似空間上で残差を最小化することで数値解を得る方法であることを、式変形から順に示す。材料計算でFEMが強い理由は、複雑形状・異方性・連成・境界条件を、同一の定式化で扱える点にある。

11.1 弱形式と重み付き残差法

有限要素法の出発点は、強形式(微分方程式そのもの)をそのまま点ごとに満たすのではなく、残差を重み付きで平均してゼロにするという発想である。具体例として、拡散方程式(ポアソン方程式の定常版)を考える。 領域 Ω で未知関数 u(x)

(ku)=fin Ω

を満たし、境界 Ω はディリクレ境界 ΓD とノイマン境界 ΓN に分割されるとする。

u=u¯on ΓD,(ku)n=t¯on ΓN

ここで n は外向き法線である。

強形式を満たさないときの残差を

R(u)=(ku)f

と定義し、任意の重み関数(試験関数)w(x) に対して

ΩwR(u)dΩ=0

を課す。これが重み付き残差法の基本である。FEMでは、後に w を近似空間と同じ空間(ガラーキン法)から取る点が重要である。

上式を弱形式に変形するため、部分積分(発散定理)を使う。以下の恒等式を用いる:

Ωw(ku)dΩ=Ωw(ku)ndΓΩkwudΩ

よって

ΩkwudΩ=ΩwfdΩ+Ωw(ku)ndΓ

となる。境界積分は ΓDΓN に分けて

Ωw(ku)ndΓ=ΓDw(ku)ndΓ+ΓNwt¯dΓ

である。ここでディリクレ境界上では試験関数に w=0 を課す(許容関数の定義)ことで、ΓD の境界積分が消える。したがって弱形式は

ΩkwudΩ=ΩwfdΩ+ΓNwt¯dΓ

である。

この変形には大きな意味がある。第一に、強形式では u に2階微分が必要であったが、弱形式では uw に1階微分があればよくなる。第二に、ノイマン境界条件が自然に右辺へ現れ、境界条件の組み込みが統一される。第三に、材料力学の釣合式や電磁場方程式でも、同様の部分積分により弱形式が得られ、同じ枠組みで離散化できる。

材料力学でも同様である。小ひずみ線形弾性では、釣合式

σ+b=0

と、構成式

σ=C:ε,ε=12(u+(u)T)

を用いる。試験関数 w を掛けて積分し、部分積分すると、内部仕事と外力仕事の等価性

Ωε(w):σ(u)dΩ=ΩwbdΩ+ΓNwt¯dΓ

が得られる。これが仮想仕事の原理であり、FEMの材料力学がエネルギー原理と直結する入口である。

11.2 形状関数と要素

弱形式が得られたら、次は未知関数 u を有限次元の近似空間に制限する。領域 Ω を要素(element)に分割し、要素内部で多項式近似を行う。 スカラー場の基本形は

uh(x)=a=1NdofNa(x)Ua

である。Na は形状関数(shape function)、Ua は節点値(自由度)である。ガラーキン法では試験関数も同じ形で

wh(x)=aNa(x)Wa

と置く。これにより連続問題が代数方程式へ落ちる。

一次要素(線形要素)では形状関数は一次多項式であり、計算が軽く頑健である。一方、応力集中や曲率の大きい場では高次要素が有利になる。高次要素は同じメッシュでも自由度が増え、滑らかな解に対して高精度になりやすいが、数値積分や要素歪みの影響が表に出やすい。

要素の形としては、2次元なら三角形・四角形、3次元なら四面体・六面体が基本である。一般に

  • 四面体は複雑形状へ自動メッシュ生成がしやすい
  • 六面体は同自由度あたりの精度が高くなりやすい という傾向があるが、最終的には目的の場(境界層、界面、き裂)に合わせて選ぶ必要がある。

代表的な要素の整理を示す。

要素自由度配置長所注意点
三角形一次(P1)3節点で線形自動メッシュに強い曲げや応力で硬くなりやすい場合がある
四角形一次(Q1)4節点で双線形構造格子で効率が良い形状が歪むと精度が落ちやすい
四面体一次(P1)4節点で線形3D複雑形状に強い低次では誤差が残りやすい
六面体一次(Q1)8節点で三線形効率と精度が良いメッシュ生成が難しい場合がある
二次要素(P2/Q2)中間節点を追加滑らかな場に高精度積分・要素品質の影響が大きい

さらに、不連続ガラーキン(DG)では要素間で連続性を課さず、数値流束で結合する。DGは移流支配問題や不連続解に強いが、自由度が増えやすく、安定化項の設計が重要になる。材料計算でDGを使う場面は、衝撃波や破壊の不連続、あるいは高Péclet数の輸送などで現れる。

形状関数は実装上、参照要素(reference element)上で定義し、実要素へ写像する。参照座標 ξ と実座標 x を結ぶ写像を x=x(ξ) とすると、ヤコビアン

J=xξ

が現れる。勾配は

xN=JTξN

で変換され、積分は

Ωeg(x)dΩ=Ω^g(x(ξ))|detJ|dΩ^

となる。異方性材料や大変形では、この写像とヤコビアンが数値安定性の中心となる。

11.3 構成式と非線形性

材料の構成式は、非線形性の主要な供給源である。弾塑性、粘弾性、損傷、相変態などでは、応力がひずみだけでなく内部変数にも依存する。 一般形として

σ=σ(ε,q)

と書き、q は塑性ひずみ、硬化変数、損傷変数などである。これらは各積分点で更新され、全体の平衡反復と結合する。

非線形FEMでは、弱形式から得られる残差方程式

R(U)=0

をニュートン法で解くのが基本である。反復 k から k+1 への更新は

Ktan(U(k))ΔU=R(U(k)),U(k+1)=U(k)+ΔU

である。ここで Ktan は接線剛性(ヤコビアン)であり、これが正しくないと収束が著しく悪化する。

塑性の基本例として、降伏関数 f(σ,κ)0 を導入する。f=0 が降伏面であり、塑性流動則を

ε˙p=λ˙gσ

と書く(g は塑性ポテンシャル、λ˙ は塑性乗数)。整合条件として

λ˙0,f0,λ˙f=0

が課される。数値的には、各時間ステップで試行応力を計算し、降伏していれば戻し写像(return mapping)で内部変数を更新する。

ここで重要なのがアルゴリズム的一貫接線(consistent tangent)である。これは、内部変数更新を含めた離散化のもとで

δσ=Calg:δε

となる有効接線 Calg を導き、これを Ktan に使うことを意味する。Calg を正しく入れると、ニュートン法が理論通り2次収束し、連成問題でも計算が成立しやすくなる。逆に、弾性接線のまま回すと、収束しないか、過度に小さい時間刻みに頼ることになりやすい。

大ひずみ問題では、回転の取り扱い(客観性)や応力率の選択が重要になる。例えば変形勾配 F を用いて

F=xX

X が参照配置)とし、エネルギー密度を W(F) で定める超弾性では、第一Piola応力

P=WF

から弱形式を構成する。材料モデルの式そのものが大きく変わるため、線形弾性の直感をそのまま持ち込まないことが重要である。

11.4 時間依存問題と安定化

拡散、粘性、波動、電磁場などは時間依存問題である。空間をFEMで離散化すると、一般に

MU˙+KU=F

のような常微分方程式(または微分代数方程式)になる。M は質量行列や容量行列、K は剛性行列に相当する。ここから時間積分法の選択が、安定性と計算量を決める。

拡散方程式の単純な例として

ut=(Du)+s

を弱形式化して離散化すると

MU˙+KU=F

となる。ここで

Mab=ΩNaNbdΩ,Kab=ΩDNaNbdΩ

である。陽的オイラーを使うと

MUn+1UnΔt+KUn=Fn

であり、安定条件が厳しい。陰的オイラーなら

MUn+1UnΔt+KUn+1=Fn+1

となり、拡散支配では一般に無条件安定であるが、毎ステップ線形系を解く必要がある。材料計算では拡散係数の温度依存や強不連続があり、陰解法と反復解法の組み合わせが基本になる。

移流項が強い輸送(高Péclet数)では、標準ガラーキン法が数値振動を起こしやすい。移流拡散方程式

ut+vu=(Du)

に対し、離散化で非物理振動が出るのは、移流が一方向情報伝播であるのに対して、標準ガラーキンが中心差分的に振る舞うためである。これを抑えるため、SUPG(Streamline Upwind/Petrov-Galerkin)などの安定化を導入する。SUPGでは試験関数を修正し

ww+τvw

のように流線方向へバイアスを与える。τ は要素サイズと速度、拡散係数から決まるパラメータであり、過大にすると過剰な数値拡散となるため、物理拡散との区別が必要である。

波動問題では質量行列が本質的である。例えば弾性波なら

ρu¨=σ+b

を離散化して

MU¨+KU=F

となる。時間積分はNewmark法などが使われ、数値分散と減衰の制御が重要になる。材料の高周波応答や衝撃問題では、安定性だけでなく位相精度(波の速度の再現)を意識する必要がある。

11.5 メッシュと幾何

FEMの精度は、要素次数だけでなくメッシュ品質に強く依存する。要素がつぶれている、角度が極端に小さい、アスペクト比が大きいなどの場合、ヤコビアンが悪条件になり、勾配評価が不安定になる。特に異方性材料や大変形では、要素品質がそのまま収束性へ影響する。

複雑形状ではCAD形状とメッシュの整合が課題になる。境界が曲線・曲面である場合、低次要素では幾何誤差が支配的になりやすい。等価的に、解を近似する前に形状自体を粗く近似しているためである。曲面境界が重要な問題では、幾何の近似次数を上げる(高次要素、アイソパラメトリック要素)ことが有効である。

局所的に急峻な場が現れる領域では局所細分化が必要になる。き裂先端、界面端、転位に相当する特異場では、解が滑らかでないため高次要素よりも細分化が効く場合が多い。適応メッシュの考え方は、誤差推定子 η を用いて

  1. 解を計算する
  2. 各要素の誤差指標 ηe を評価する
  3. ηe が大きい要素を細分化する を繰り返すことである。材料計算では、応力やひずみの二次量、あるいはエネルギー密度の勾配が指標として使われることが多い。

誤差推定には大別して2系統ある。残差型推定は、離散解を強形式に代入した残差を使い、要素内部残差と要素境界でのフラックス不連続を合成して評価する。回復型推定は、要素内で不連続になりやすい勾配(応力など)を滑らかに再構成し、その差を誤差とみなす。どちらも、厳密な上界評価が得られる場合と、経験的に有効な場合があるため、目的量(例えば最大応力、変位、J積分)に対して妥当な推定になっているかを確認する必要がある。

11.6 誤差評価と収束

FEMの誤差は、メッシュ幅 h と要素次数 p に対して体系的に議論できる。滑らかな解 u に対して、p 次要素での近似誤差は概ね

uuhChp+1

のように減少する(ノルムの種類により指数が変わる)。勾配を含むエネルギーノルムでは

uuhEChp

となることが多い。材料力学ではエネルギーに基づく誤差評価が自然であり、内部仕事の差として整理できる利点がある。

重要なのは、誤差が場そのものではなく、目的の物性値へどう伝播するかである。例えば応力集中の最大値は、局所的な特異性に敏感で、メッシュを細かくしても収束が遅い場合がある。逆に、全体平均のエネルギーや平均ひずみは比較的早く収束する。したがって、何を予測したいかを明確にし、その量に対してメッシュ収束を検証する姿勢が必要になる。

目的量 J(u) の誤差評価では、双対問題(随伴問題)を用いたゴール指向推定が有効である。概念としては、誤差が観測される方向を双対解で重み付けして評価する。実装では

  • 主問題で uh を解く
  • 双対問題で λh を解く
  • 残差と双対解から J(u)J(uh) を推定する という流れになる。材料同定や最適化では、目的量に対する精度が最重要であるため、この考え方が直接役に立つ。

第12章 有限要素法による連成材料計算

本章では、FEMが複数の物理場を同一の弱形式枠組みで結合できる点を活かし、熱・応力・拡散、電磁場、フェーズフィールドなどの連成をどのように組み立てるかを示す。連成では方程式が増えるだけでなく、時間スケールや支配パラメータが異なる場が結び付くため、解法設計と検証が材料科学的な結論の信頼性を左右する。

12.1 熱・応力・拡散の連成

熱と応力の連成の基本は熱膨張である。温度変化 ΔT による熱ひずみを

εth=αΔTI

とすると、線形弾性の構成式は

σ=C:(εεth)

となる。ここで α は線膨張係数、I は単位テンソルである。熱伝導は

ρcpTt=(kT)+Q

であり、Q が塑性発熱や相変態潜熱を含む場合、温度場は変形と双方向に結合する。

拡散と応力の連成では、化学ポテンシャルに応力項が入る点が要である。組成 c の化学ポテンシャルを

μ=μ0(c,T)+μmech

と分け、体積変化を伴う溶質(化学膨張)では

μmech=ΩσH

のように静水圧応力 σH=13tr(σ) が入る(Ω は部分モル体積に対応する係数)形で表されることがある。すると拡散流束は

J=Mμ

となり、応力勾配が拡散を駆動しうる。これが応力誘起拡散や応力腐食割れの基本機構の一部になる。

弱形式では、変位 u、温度 T、濃度 c を未知関数として同時に扱い、各方程式の試験関数 wsq を用意する。例えば拡散の弱形式は

ΩqctdΩ+ΩqMμdΩ=ΓNcqj¯dΓ

のようになる。ここで ΓNc は濃度のフラックス境界、j¯ は外部からの供給である。μu,T,c に依存すれば、これだけで連成が成立する。

連成解法にはモノリシック解法と分割解法がある。モノリシックでは未知ベクトルを

X=[UTC]

のようにまとめ、非線形方程式 R(X)=0 を一括で解く。分割解法では、例えば T を解いてから u を解くなど、場を順に更新する。

解法特徴注意点
モノリシック強連成でも頑健になりやすいヤコビアン構築と前処理が難しい
分割実装が分かりやすい強連成で不安定になりやすく反復が必要

材料プロセスでは温度と変形が強く結び付く場合が多く、分割解法では反復(固定点反復)を入れて収束を確保する設計が必要になる。どちらを選ぶかは、連成の強さ、時間刻み、計算規模、利用できるソルバ基盤で決めるのが現実的である。

12.2 電磁場と材料応答

電磁場連成ではマクスウェル方程式が出発点である。時間領域の基本形は

×E=Bt,×H=J+Dt

および

B=0,D=ρe

である。材料応答は構成式

B=μH,D=ϵE,J=σE

で入る。導電率 σ、透磁率 μ、誘電率 ϵ は材料・温度・周波数で変化し得るため、連成が強くなる。

渦電流問題(準静的近似)では、周波数が高すぎない範囲で変位電流を無視し、

×H=J,J=σE

とし、ベクトルポテンシャル A を導入して

B=×A,E=Atϕ

と置く。ゲージ条件を適切に選ぶと、未知数は Aϕ に整理できる。周波数領域(調和定常)では時間微分が iω に置き換わり、例えば

E=iωAϕ

となる。ここで物性が複素化する意味が明確になる。例えば導電体では、電流が位相遅れを持ち、損失(ジュール熱)を生む。損失密度は

p=JE

であり、調和場では時間平均が

p=12Re(JE)

で与えられる。FEMはこのような複素場の弱形式にも自然に対応する。

電磁場と材料変形の連成としては、ローレンツ力

f=J×B

が力学方程式の外力項に入る。これにより電磁成形や渦電流ブレーキが記述できる。磁歪などの磁気弾性連成では、自由エネルギーに磁気エネルギーと弾性エネルギーの結合項を入れ、磁化(または磁気ポテンシャル)と変位を同時に解く設計になる。周波数領域での複素透磁率は、損失を含む実効応答を表しており、実部が蓄積、虚部が散逸に対応する。材料計算では、どの物理を複素物性へ押し込めているのかを明確にする姿勢が重要である。

12.3 フェーズフィールドと有限要素法の統合

フェーズフィールド(PF)方程式は弱形式に落とせるため、FEMと自然に統合できる。例えばAllen–Cahn型

ϕt=LδFδϕ

を、自由エネルギー

F[ϕ]=(f(ϕ)+κ2|ϕ|2)dΩ

で定めると、化学ポテンシャルに対応する変分は

δFδϕ=fϕκ2ϕ

である。試験関数 w を掛けて積分すると

ΩwϕtdΩ=ΩwL(fϕκ2ϕ)dΩ

右辺のラプラシアン項を部分積分して

ΩwϕtdΩ=ΩwLfϕdΩΩκLwϕdΩ+ΓκLwϕndΓ

となる。境界条件により境界積分を整理すれば、FEM離散化は通常の手順で行える。これによりPFを複雑形状へ適用しやすくなる。

Cahn–Hilliard型は4階PDEであるため、扱い方に工夫が要る。方法は2つある。 1つ目は C1 連続な要素(高次の特殊要素)を使って直接離散化する方法である。2つ目は補助変数を導入して2つの2階方程式に分解する方法である。例えば

μ=fcκ2c,ct=(Mμ)

とし、未知数を (c,μ) に増やすと、各式は2階である。FEMではこの分割が実装上有利で、連成問題でも標準的な C0 要素で扱える。

弾性場との連成は、PFの自由エネルギーに弾性項を加えて行う。例えば

F[ϕ,u]=(f(ϕ)+κ2|ϕ|2+12(εε0(ϕ)):C:(εε0(ϕ)))dΩ

である。すると ϕ 方程式の駆動力に fel/ϕ が加わり、逆に弾性方程式には固有ひずみが入る。PF-FEMは析出形態、応力誘起相変態、相界面のピン止めなどの議論へ直接つながる。

ここで数値誤差の扱いが重要になる。PFは界面幅 をモデル内で持つため、メッシュ幅 h に対して十分小さい必要がある。一般に界面を N 点で解像するとき h/N であり、N が小さいと界面エネルギーや移動度が数値的に変化してしまう。したがって、PFパラメータと離散化誤差を独立に議論するため、界面幅の収束試験とメッシュ収束試験を分けて行う必要がある。

12.4 逆問題と最適化

材料定数同定や形状最適化は、観測量と設計変数を結ぶ最適化問題として定式化できる。設計変数を θ、FEMで得られる状態変数を U(θ)、観測量を y(θ) とし、目的関数を

J(θ)=12y(θ)yobs2+λR(θ)

のように定める。R は正則化であり、測定ノイズや未同定性に対して解を安定化する役割を持つ。

問題は、y(θ) がFEM解を通じて定まるため、勾配 θJ をどう計算するかである。設計変数の次元が大きいと、有限差分で勾配を取るのは計算量的に不可能になりやすい。ここで随伴法(adjoint method)が有効である。

離散化後の支配方程式を

R(U,θ)=0

と書き、ラグランジュ関数

L(U,θ,λ)=J(U,θ)+λTR(U,θ)

を導入する。λ は随伴変数である。LU 微分がゼロになる条件

JU+(RU)Tλ=0

を満たす λ を解けば、勾配は

dJdθ=Jθ+λTRθ

で計算できる。ここで重要なのは、随伴方程式を1回解けば、設計変数の次元に依らず勾配が得られる点である。材料同定(多パラメータ)や形状最適化(多数自由度)でこの差は決定的である。

ただし、随伴法は目的関数の定義に強く依存する。例えば変位の誤差を最小化するのか、応力の誤差を最小化するのか、固有振動数を合わせるのかで随伴方程式が変わる。実験が提供する観測量と、モデルが出す量を同一の意味に合わせたうえで目的関数を設計する必要がある。

12.5 大規模計算基盤

連成材料計算では自由度が大きくなり、疎行列線形系・非線形反復・時間積分を効率良く回す基盤が不可欠である。離散化で得られる線形系は

Ax=b

の形であり、A は疎行列である。連成ではブロック構造を持つことが多く、例えば

[AuuAuTAucATuATTATcAcuAcTAcc][ΔUΔTΔC]=[RuRTRc]

のようになる。ここで重要なのは、単に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_filetoに相当する流れで、構造データを共通表現へ揃える設計になっている。

この段階で重要なのは、構造だけを集めて満足しないことである。後でどの条件で得た構造かを追跡できなければ、同じ化学式の別相や、同一相の格子定数違いを混同しやすくなる。したがって、ファイル由来(計算条件、実験条件、試料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側で重要なのは、取得したデータを最終的にStructureCompositionへ揃え、同じ検査規則で品質を確認できる状態に持ち込むことである。

国内DBは、データの種類が計算構造だけでなく実験条件・計測データ・プロセス情報を含みやすい。したがって、構造はCIF等で取り出してPymatgenへ揃えつつ、同時に条件情報(温度、雰囲気、装置、試料履歴)を別レコードとして紐づける設計が必要になる。ここで構造だけPymatgen、条件はバラバラのままにすると比較が破綻するため、構造オブジェクトとメタ情報を同一の識別子で結び続ける運用が重要である。

第14章 Pymatgenの活用方法II

本章では、収集した膨大な材料データから、研究目的に合致する情報だけを選び、同じ定義で比較できるデータ集合へ整える手順を記述する。ポイントは、検索条件と選別条件を、物理・化学の量として宣言し、第三者が同じ規則で再現できる状態にすることである。そのために、Pymatgenの解析機能と、データモデル(構造・組成・エントリ)を中心に据える。

14.1 選別条件を物理量として定義する

ビッグデータからの選別は、感覚的な絞り込みではなく、条件を数式と変数で表すところから始めるべきである。例えば安定な新材料候補を集めたいなら、形成エネルギーやhullからの距離(energy above hull)に基づく条件が自然である。形成エネルギーの基本形は、元素参照エネルギーμiを用いて

Eform=Etotiniμiini

であり、比較は通常原子あたりで揃える。ここでμiの定義(基準相、温度、計算設定)を混在させると、選別が物理ではなく基準差を拾ってしまうため、同じDB内の整備された値を使うか、参照状態を固定して再計算する必要がある。

また、狙いが磁性・触媒・電池などで異なるなら、選別条件も変わる。例えば磁性材料なら、磁気秩序や磁化の大きさだけでなく、結晶対称性(異方性の候補)や、安定性と合成容易性の関係も同時に見ることが多い。したがって何を最大化し、何を下限として課すかを、最初に変数として整理し、後段の多目的選別へ接続できる形にしておくべきである。

14.2 構造と組成で重複を除き、比較単位を揃える

公開DB由来のビッグデータには、重複や近似的に同じ構造が混入しやすい。例えば、同一材料でも最適化条件やセル選びの違いで格子がわずかに変わり、別エントリとして残る場合がある。ここで重複を放置すると、同じ材料が何度も学習データに現れ、汎化性能の評価が壊れやすくなる。Pymatgenは構造表現を統一できるため、比較の前に同一視の規則を決め、同じ規則でデータ集合を正規化することが重要である。

比較単位も同様に揃える必要がある。例えば、バンドギャップはeV、体積は\AA3、磁化はμB/f.u.A/mなど表現が分かれるため、同じ量を比較しているかを常に確認しなければならない。Pymatgenは構造を通じて原子数、化学式単位、密度などの換算に必要な情報を保持するため、単位換算を設計に組み込み、後から場当たり的に直す状況を避けるべきである。

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は電子密度 n(r) を主変数とし、MDは原子座標 {ri(t)} を時間発展させ、MCは状態 x を確率分布 π(x) でサンプリングし、PFは場(秩序変数・濃度)ϕ(r,t),c(r,t) を連続体として扱い、FEMは複雑形状上で多物理場(変位・温度・電磁場など)を弱形式で統合する。

この違いは、「何を平均して捨てたか」に対応する。したがって手法間の接続は、単に数値を渡すのではなく、平均化の定義(どの自由度を消したか)を明示しないと物理がずれる。典型的な“情報変換”を、式の形で整理する。

(1) DFT → 連続体材料定数(弾性・磁気・拡散など) DFTから得られるのは、エネルギー E(基底状態エネルギー)を中心とする量である。連続体材料定数は、エネルギーの微分として定義されることが多い。

例:弾性定数(小ひずみ) ひずみテンソル ε に対して、単位体積あたりのエネルギー密度 W(ε) を考えると、

σij=Wεij,Cijkl=2Wεijεkl

で弾性定数が定まる。実際の計算では、ひずみを微小に振って E(ε) を複数点計算し、2次多項式フィットで2階微分を推定する。このとき重要なのは「ひずみの定義」「体積一定かどうか」「内部原子緩和を許すか(clamped/relaxed)」により、得られる定数が変わる点である。

例:拡散(活性化エネルギー) 拡散係数はしばしばアレニウス形

D=D0exp(EakBT)

で整理され、Ea は遷移状態(NEBなど)から得る活性化障壁である。このとき D0 は振動数(試行頻度)や跳躍距離を含むため、DFTだけでなくフォノン近似やMD/kMCの補助が必要になる。

(2) DFT → クラスター展開 → MC(合金相図) 置換合金の配置エネルギーを格子上の有効ハミルトニアンで表し、

E(σ)=J0+αJαΦα(σ)

σ は占有変数、Φα はクラスター相関、Jα は有効相互作用)と書ければ、MCで有限温度の相図や短距離秩序を計算できる。ここでは「どの格子を採用するか」「クラスターの打ち切り」「DFT計算点の選び方」が、精度と外挿性を左右する。単にフィットが良いだけでは不十分で、物理対称性や外挿の健全性をチェックする必要がある。

(3) MD → 連続体(輸送係数・界面物性) MDは時間相関から輸送係数を与える。例えば拡散係数は平均二乗変位(Einstein関係)で

D=limt16t|r(t)r(0)|2

粘性係数はグリーン久保関係で

η=VkBT0Pxy(0)Pxy(t)dt

のように与えられる。ここから連続体へ渡すときは、(i) 有限サイズ補正(特に拡散)、(ii) 平衡到達と統計誤差(自己相関)、(iii) 力場/MLポテンシャルの妥当性(外挿)を同時に評価する必要がある。

(4) PF ↔ FEM(組織進化と複雑形状・連成) PFは自由エネルギー汎関数 F の変分から駆動力を作る手法であり、FEMは弱形式で複雑形状・連成を扱う基盤である。両者は自然に統合できるが、接続点は「界面幅」「メッシュ幅」「物理パラメータ」の関係である。 例えば界面幅 をPFで導入したなら、数値的には を十分解像する必要があり、h(少なくとも数点以上)を満たさないと、界面エネルギーや移動度が数値依存になる。つまりPFの“モデルパラメータ”とFEMの“離散化誤差”を混ぜない設計が必須である。

まとめると、手法間の接続は「量を渡す」ではなく「定義を渡す」である。どの平均化で得た量か、どの境界条件・拘束の下で定義された量かを、式の形で記録して初めて、マルチスケールが科学として成立する。

17.2 計算結果の信頼性の作り方

計算の信頼性は「再現できる」だけでは十分ではない。研究としては、(i) モデル仮定、(ii) 数値誤差、(iii) 統計誤差、(iv) 実装差(コード・擬ポテンシャル等)を分解し、観測量(結論)へどう影響するか説明できる必要がある。ここでは“信頼性の作り方”を、最低限のチェック項目ではなく、一本の手順(体系)として整理する。

(1) 収束試験:離散化・基底・サンプリング どの手法でも「解像度を上げたら答えが安定するか」を確認する。

  • DFT:カットオフ Ecutk 点密度、スメアリング幅、セルサイズ、擬ポテンシャル/PAWの選択
  • MD:時間刻み Δt、カットオフ半径、系サイズ、熱浴・圧力制御の時定数、平衡化時間と観測時間
  • MC:熱化ステップ数、サンプル数、自己相関時間、レプリカ数(拡張アンサンブル)
  • PF/FEM:メッシュ幅 h、時間刻み Δt、要素次数 p、線形/非線形ソルバ許容誤差、前処理

重要なのは、誤差が単純にゼロへ行くわけではなく、目的量によって収束が違う点である。例えば最大応力のような局所量は収束が遅く、全エネルギーや平均量は収束が早いことが多い。したがって「何を結論として使う量か」を先に決め、その量の収束を優先して見る必要がある。

(2) 相互検証:別経路で同じ量を確かめる 同じ物理量を、異なる定義や異なる手法で評価して一致を確認する。

例:熱伝導

  • 平衡MD(Green–Kubo)
  • 非平衡MD(温度勾配を与える) の2系統で照合し、サイズ依存や境界条件依存を切り分ける。

例:PFの界面エネルギー PFの勾配係数 κ と二重井戸ポテンシャルのパラメータから界面エネルギー γ が解析的に得られる場合、その値と数値的に測った界面張力が一致するかを確認する。

(3) 感度解析:仮定・パラメータが結論を支配していないか 結論が特定のパラメータに過度に依存していないかを見る。一次感度は

Sθ=Qθ

で定義できるが、実際には無次元化して

S~θ=θQQθ

のように相対感度で比較すると、どの仮定が支配的かが見えやすい。連成問題では、ソルバ許容誤差や時間刻みが物性値を動かしてしまう場合もあるため、物理パラメータだけでなく数値パラメータも感度対象に含めるのが実務的である。

(4) 不確かさ(統計誤差・モデル誤差)を“数で”出す MC/MDでは統計誤差が本質的であり、自己相関を考慮した有効サンプル数 Neff を意識する。統合自己相関時間 τint を用いると

NeffN2τint

と見積もられ、誤差は概ね 1/Neff で減る。ここを無視すると「計算は長く回したが精度が上がっていない」状態になる。

モデル誤差(力場・汎関数・粗視化の誤差)は、単一指標で出しにくいが、少なくとも

  • 参照データ(DFT/実験)との比較
  • 外挿検知(MLポテンシャル)
  • 物理制約(保存則・対称性)の破れの監視 をセットで行うと、説明可能性が上がる。

信頼性とは、単に「正しい」ではなく「どこまで正しいと言えるか」を定量化し、読者が判断できる材料を揃えることだと捉えるのがよい。

17.3 データ駆動との融合

データ駆動は材料計算を加速するが、万能ではない。むしろ初学者にとって重要なのは、「どこをデータで置き換え、どこを物理で固定するか」を設計できることにある。ここでは、データ駆動を“追加の便利機能”ではなく、“不確かさを伴う近似器”として取り扱う。

(1) MLポテンシャルは「近似した力学系」である MLポテンシャルは、原子配置 {ri} からエネルギー E を予測し、力を

Fi=Eri

として得る。したがって、エネルギーが滑らかでない・学習誤差が大きいと、力の誤差は増幅されやすい。ここから、学習時に「エネルギーだけでなく力も教師信号に入れる」「保存則(エネルギー保存)を監視する」などの設計が必要になる理由が理解できる。

(2) 外挿を検出しないと科学にならない 学習モデルは訓練分布の外側で誤る。外挿検出の代表はアンサンブルで、複数モデルのばらつき

Var[E(r)]

を不確かさ指標として用いる。これをMDの途中で評価し、指標が大きい構造をDFTで追加計算して学習データへ戻す(アクティブラーニング)と、外挿を系統的に潰せる。ここで重要なのは、誤差が小さく見える範囲でも、観測量(例えば拡散係数)が大きく変わる場合があるため、“目的量に対する不確かさ”を意識することである。

(3) 逆問題は「同定できないことがある」問題である 材料定数同定やモデル同定は、観測 y に対してパラメータ θ を推定する問題であるが、写像 y(θ) が一対一でない場合、複数の θ が同じ y を説明できる(非同定)。これを避けるには

  • 観測量を増やす(異なる条件、異なる測定)
  • 正則化を入れる(事前知識で解を絞る)
  • 感度の高い実験条件を選ぶ(実験計画) が必要になる。データ駆動は、ここを“自動で解く”のではなく、“何が同定できていないか”を可視化する道具として使うのが安全で強い。

(4) 物理制約とデータの役割分担 データ駆動を入れるときの基本戦略は、

  • 物理で守る:保存則、対称性、境界条件、単位系、極限挙動(高温・低温など)
  • データで補う:未知の構成式、経験則、複雑な相互作用の近似 という分担である。これにより、外挿や偽相関の危険を抑えつつ、探索空間を拡張できる。

17.4 計算機環境の発展

材料計算は、アルゴリズムと計算機の両方の制約の上に成立する。今後はGPUや異種混在計算が増え、計算量(flops)よりもメモリ帯域・通信が支配的になる傾向が強い。初学者が押さえるべきポイントを、手法ごとに整理する。

(1) DFT:対角化と通信が重い 平面波DFTでは、Kohn–Sham方程式の解法がボトルネックであり、多くの時間が固有値問題(対角化)に費やされる。系サイズが増えると、計算量は概ね急増し、通信も増える。したがって、単純に“コア数を増やす”だけで速くならない領域がある。並列化は

  • k 点並列
  • バンド並列
  • 平面波(FFT)並列 などの分割を組み合わせて設計される。

(2) MD:粒子数にほぼ比例、ただし長距離相互作用が効く 短距離相互作用なら計算量は粒子数 N にほぼ比例する。一方、クーロン相互作用はEwald/PMEで NlogN 程度になり、FFTが支配する。GPUでは粒子計算が非常に速いが、FFTや通信が律速になるため、アルゴリズムのメモリ局所性が重要になる。

(3) FEM/PF:疎行列ソルバと前処理が支配 大規模FEMは疎行列の反復解法が中心で、前処理の出来が速度を決める。GPUは疎行列演算が得意だが、前処理が難しいことも多い。したがって「どの前処理が使えるか」を前提に、離散化と未知量の並べ方(ブロック構造)を決めるのが実務的である。

(4) 再現性と高速化の緊張関係 並列計算では削減演算の順序が変わり、丸め誤差が非決定性を生む場合がある。これはバグではなく数値の性質であり、対策は

  • 収束判定を残差だけでなく物理量でも行う
  • 乱数種、並列数、ライブラリ版数を記録する
  • 重要結果は条件を変えて再現する である。GPU時代ほど、この“再現性をどう作るか”が重要になる。

17.5 研究倫理と公開可能性

計算研究の倫理は、捏造や改ざんだけではなく、「他者が追試できる状態で提示しているか」に強く関わる。材料計算は入力パラメータが多く、ちょっとした違いで結果が動くため、公開可能性(reproducibility)を意識した記録が不可欠である。

最低限記録すべき情報は、手法により異なるが、共通して以下が重要である。

  • 入力ファイル(モデル・境界条件・初期条件)
  • ソフトウェア名とバージョン、コンパイル条件
  • 擬ポテンシャル/PAWデータ、汎関数名
  • 乱数種(MC、MDの初期速度など)
  • 収束条件(残差、許容誤差、反復上限)
  • 解析スクリプト(後処理も含む)
  • 単位系と物理量の定義(特に“エネルギー密度”など)

ここで重要なのは、結果の図だけでは不十分だという点である。例えば「拡散係数」と言っても、Einstein関係を使ったのか、速度自己相関を使ったのか、有限サイズ補正を入れたのかで値が変わり得る。つまり“量の定義”が公開情報に含まれて初めて科学になる。

さらに、公開の粒度も設計が必要である。全てを出すのが難しい場合でも、

  • 最小再現パッケージ(入力+解析手順+要点の出力)
  • 代表ケースの完全再現(全条件のうち一部を完全公開) のどちらかを提供できると、コミュニティでの信頼が大きく上がる。教育の段階でこれを習慣化すると、研究室運営や共同研究でもメリットが大きい。

17.6 学習順序の提案

初学者が材料計算を自力で組み立てられるようになるには、「物理」と「数値」を別々に覚えるのではなく、同じ方程式が“どの形で解かれているか”をセットで理解するのが近道である。ここまでの内容を踏まえ、学習順序を具体的に提案する。

(1) まず共通基盤:モデル化と離散化(第1・2章) 保存則

qt+J=s

と、自由エネルギーの変分

μ=δFδc,ϕt=LδFδϕ

を“読める”ようになると、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)――を、単なる計算法の羅列としてではなく、「どのスケールを、どの未知量で記述するのか」という視点から体系的に整理してきた。原子・電子スケールから連続体スケールへと記述が移り変わる中で、計算手法の違いは精度や計算コストの差ではなく、物理量の定義と平均化の仕方の違いとして理解できることが、本書を通じて明確になったはずである。

また、異なる手法を接続する際に本質的に重要なのは、数値そのものを無理に受け渡すことではなく、「何を平均し、何を捨て、どの量を有効変数として残すのか」という定義の選択であることを繰り返し強調してきた。この視点を持つことで、マルチスケール計算は特別な技術ではなく、記述の階層を意識的に行き来するための論理的操作として捉えられるようになる。

計算結果の信頼性についても、本書では単一の基準で語ることを避け、収束性、相互検証、感度解析、不確かさの評価という複数の要素に分解して考える立場を採った。信頼性とは結果に付随する性質ではなく、計算の設計・実行・解釈の過程を通じて構築されるものである。この分解的な理解は、新しい手法や未経験の計算に直面した際にも、自ら判断軸を立てるための基盤となる。

さらに、近年重要性を増しているデータ駆動型手法については、万能な代替手段としてではなく、物理制約、外挿領域の検出、逆問題に内在する非同定性といった限界を意識したうえで使うべき道具として位置づけた。計算科学とデータ科学は対立するものではなく、記述能力と推論能力を補完し合う関係にあることを理解することが、今後の材料研究において不可欠である。

最終的に本書が目指したのは、特定の手法を使えるようになることではなく、学習の順序そのものを自分で設計できる研究者・技術者を育てることである。いま自分はどのスケールを理解しており、どの未知量が曖昧なのか、次に学ぶべき理論や手法は何か――それらを自ら判断できるようになったとき、材料計算学は単なる計算技術ではなく、思考の道具として機能し始める。本書が、そのための出発点となることを期待したい。

参考ドキュメント

  1. 東京大学物性研究所 計算物質科学研究センター, 第一原理計算プログラムとOpenMXについて(PDF): https://ccms.issp.u-tokyo.ac.jp/wp-content/uploads/2019/12/openmx_text.pdf
  2. HPCI, PHASE/0(アプリケーション紹介): https://www.hpci-office.jp/for_users/appli_software/appli_phase0
  3. 小山敏幸, フェーズフィールド法による組織形成シミュレーション(日本語解説, 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