第一原理分子動力学法の基礎（AIMD）

第一原理分子動力学法（ab initio molecular dynamics, AIMD）は、原子核の運動を分子動力学で追いながら、各時刻の原子に働く力を第一原理電子状態計算からその都度求める方法である。経験的ポテンシャルの関数形を仮定せず、電子の量子力学を通じて化学結合の生成・切断や金属結合の集団的応答を同一の枠組みで扱える点に強みがある。

目次

• 位置づけ
• 古典分子動力学の出発点
• 電子と原子核の分離
• エネルギー地形と運動方程式
• 力の導出
• 応力とセル自由度
• Born–Oppenheimer分子動力学
• Car–Parrinello分子動力学
• 有限温度DFTと占有数
• 時間積分と時間刻み
• 温度制御
• 圧力制御
• 周期境界条件と電荷・表面
• 電子状態計算の設定がダイナミクスへ与える影響
• 初期条件と平衡化
• 観測量の計算
• 統計誤差と相関時間
• 計算規模と到達時間
• 機械学習ポテンシャルとの併用
• 手法の比較
• まとめと展望
• 参考文献

参考ドキュメント

1. D. Marx and J. Hutter, Ab initio molecular dynamics: Theory and Implementation (NIC Series, 2000, PDF) https://www.theochem.rub.de/images/theochem/research/marx/marx.pdf
2. R. Car and M. Parrinello, Unified Approach for Molecular Dynamics and Density-Functional Theory, Phys. Rev. Lett. 55, 2471 (1985) https://link.aps.org/doi/10.1103/PhysRevLett.55.2471
3. 繁田秀哉, 第一原理ダイナミクスによる材料の物性・反応解析（HPCI講習資料, PDF） https://www.hpci-office.jp/invite2/documents2/ws_material_191018_shigeta.pdf

位置づけ

分子動力学（MD）は、原子核（原子）を質点とみなし、力学方程式を数値積分して原子配置の時間発展を得る方法である。古典MDでは力を経験的ポテンシャルから与えるが、AIMDでは力を電子の量子力学（多くは密度汎関数理論）で決めるため、力の根拠がより直接的に物理へ結び付く。

AIMDは、固体・液体・表面・界面などの凝縮系で、温度や圧力の下で原子が実際に動く様子を、電子構造を保ったまま追うために用いられる。静的な構造最適化が一つの安定構造を探すのに対し、AIMDは有限温度でのゆらぎや拡散、相転移の核形成など、時間方向の情報が本質となる現象へ踏み込む役割を持つ。

古典分子動力学の出発点

古典MDの核は、各原子$I$についてのニュートン方程式

$$M_I\frac{d^2{\mathbf{R}}_I}{d{t}^2}={\mathbf{F}}_I$$

である。ここで$M_I$は質量、${\mathbf{R}}_I$は位置、${\mathbf{F}}_I$は力である。

多くの場合、力はポテンシャル$U(\mathbf{R})$の勾配として

$${\mathbf{F}}_I=-{\mathrm{\nabla }}_{{\mathbf{R}}_I}U(\mathbf{R})$$

と定義される。$U$が時間に依らないなら、外力や温度制御を入れない限り、全エネルギー（運動＋ポテンシャル）が一定になるという保存則が見通しよく成立する。

古典MDで本質となるのは、力がどれだけ正確に$U$の勾配になっているかである。AIMDでも同じであり、電子状態計算が与える力があるエネルギー汎関数の勾配として一貫しているかが、長時間積分の安定性と観測量の信頼性を左右する。

電子と原子核の分離

電子と原子核の運動を同時に量子力学で解くことは一般に困難であるため、原子核を古典的に、電子を量子的に扱う近似が導入される。原子核は電子よりはるかに重く、時間スケールが遅いので、原子核の配置$\mathbf{R}$が変わる間に電子はほぼ準静的に応答する、という見方が成り立ちやすい。

この考え方は断熱近似（Born–Oppenheimerの考え方）として知られる。AIMDでは、原子核配置$\mathbf{R}$を固定した瞬間ごとに電子の基底状態（または所定の電子状態）を求め、そのエネルギーと力を用いて原子核を一歩進める、という構造が基本になる。

エネルギー地形と運動方程式

原子核配置$\mathbf{R}$に対して、電子状態計算が全エネルギー$E(\mathbf{R})$を与えるとき、原子核はそのエネルギー地形の上を動くとみなせる。すなわち

$$M_I{\ddot{\mathbf{R}}}_I=-{\mathrm{\nabla }}_{{\mathbf{R}}_I}E(\mathbf{R})$$

である。

この式は見た目が単純であるが、含意は重い。第一に、原子核の運動は多体かつ量子由来の相互作用を、単一のスカラー関数$E(\mathbf{R})$へ畳み込んだものとして表現される。第二に、$E(\mathbf{R})$がその場で計算されるため、AIMDは計算の各ステップで電子状態計算を反復する構造を持つ。

力の導出

AIMDで用いられる力は、電子状態計算の全エネルギーを原子位置で微分したものである。概念的には

$${\mathbf{F}}_I=-\frac{\partial E}{\partial {\mathbf{R}}_I}$$

であるが、実際には電子密度やKohn–Sham軌道、擬ポテンシャル、基底の取り方に依存して具体形が分かれる。

Hellmann–Feynmanの見方

ハミルトニアン$\hat{H}(\mathbf{R})$が原子位置に依存し、固有状態$|\mathrm{\Psi }(\mathbf{R})⟩$が得られているとき、期待値$E(\mathbf{R})=⟨\mathrm{\Psi }|\hat{H}|\mathrm{\Psi }⟩$の微分は、条件が整えば

$$\frac{\partial E}{\partial {\mathbf{R}}_I}=⟨\mathrm{\Psi }\left|\frac{\partial \hat{H}}{\partial {\mathbf{R}}_I}\right|\mathrm{\Psi }⟩$$

の形（Hellmann–Feynman型）で書ける。ここでは波動関数の微分が見えなくなるため、力の計算が直感的になる。

ただし、この簡単な形がそのまま数値実装に落ちるとは限らない。基底が原子位置に依存する場合（局在基底など）には追加項（Pulay項）が現れ、平面波基底でも有限カットオフや数値近似により力の精度が揺らぐことがある。AIMDでは一回の誤差が次のステップへ伝搬するため、力の一貫性が静的計算以上に重要になる。

密度汎関数理論の立場での力

Kohn–Sham DFTでは、エネルギーは電子密度$\rho (\mathbf{r})$（およびKohn–Sham軌道）から計算され、原子位置依存は外部ポテンシャル（原子核による）を通じて入る。したがって、力は核-電子相互作用、核-核反発、交換相関汎関数の寄与などを含み、数値的にはそれらを整合的に評価する必要がある。

AIMDで力が不安定なとき、原因はしばしば電子状態計算が十分収束していないこと、基底や離散化が粗いこと、占有数（スメアリング）の扱いが一貫していないことにある。力はエネルギーの微分であるため、エネルギー自体が見かけ上よく見えても、力の誤差が顕在化することがある。

応力とセル自由度

固体や高密度液体では、セル体積や形状が変化する状況（熱膨張、圧力下、ひずみ下）を扱いたくなる。そこで、原子位置だけでなくセルの自由度も変数として導入し、応力テンソル$\mathit{\sigma }$を用いてセルを更新する枠組みが必要になる。

応力は、エネルギーのひずみ微分として

$${\sigma }_{\alpha \beta }=\frac{1}{\mathrm{\Omega }}\frac{\partial E}{\partial {ϵ}_{\alpha \beta }}$$

（$\mathrm{\Omega }$は体積、${ϵ}_{\alpha \beta }$はひずみ）で定義される。AIMDでは応力も電子状態計算から得るが、応力は力以上に数値誤差（カットオフ、k点、占有数）に敏感であることが多い。したがって、圧力制御やセル変形を伴う計算では、応力の収束性を意識した設定が必要になる。

Born–Oppenheimer分子動力学

Born–Oppenheimer分子動力学（BOMD）は、各時刻の原子配置$\mathbf{R}(t)$に対して電子を基底状態へ十分に収束させ、そのエネルギー$E(\mathbf{R})$と力${\mathbf{F}}_I=-{\mathrm{\nabla }}_{{\mathbf{R}}_I}E$を用いて原子核を進める方法である。原子核は常に同一のエネルギー地形$E(\mathbf{R})$の上を動くとみなせるため、物理像が明確である。

BOMDの反復構造

時刻$t_n$で${\mathbf{R}}^{(n)}$が与えられると、まず電子状態計算（自己無撞着計算）を行い、電子密度と全エネルギーを求める。次に力を計算し、時間積分により${\mathbf{R}}^{(n+1)}$へ進む。これを繰り返すことで軌道が得られる。

このとき、電子状態計算は前ステップの波動関数を初期値にすることで速く収束することが多い。AIMDが進むほど構造変化が小さくなり、電子の初期推定が良くなるため、計算が安定することがある一方、急激な相転移や反応が起こると電子状態が大きく変わり、収束が難しくなることがある。

エネルギー保存とドリフト

BOMDをNVE（外部から温度制御や圧力制御を入れない）で行えば、理想的には全エネルギーが保存される。しかし数値計算では、時間積分誤差に加えて電子の未収束が力の非保存性として現れ、全エネルギーがゆっくり変化することがある。

この現象は、力が厳密にあるスカラー関数の勾配になっていない、という意味である。したがって、時間刻み$\mathrm{\Delta }t$を小さくすると改善するか、電子収束条件を厳しくすると改善するか、という観察は原因の切り分けに有効である。AIMDでは、この種の整合性が軌道の品質に直結するため、エネルギーと力の両方を同じ基準で収束させる必要がある。

Car–Parrinello分子動力学

Car–Parrinello分子動力学（CPMD）は、各ステップで電子を完全に収束させる代わりに、電子自由度にも仮想的な運動を与え、原子核と同時に連成して時間発展させる方法である。出発点は拡張ラグランジアンであり、代表的には

$$L=\sum _I\frac{1}{2}{M}_I{\dot{\mathbf{R}}}_I^2+\sum _i\frac{1}{2}\mu ⟨{\dot{\psi }}_i|{\dot{\psi }}_i⟩-E[\{{\psi }_i\},\mathbf{R}]+\sum _{ij}{\mathrm{\Lambda }}_{ij}(⟨{\psi }_i|{\psi }_j⟩-{\delta }_{ij})$$

と書ける。$\mu$が仮想質量、${\mathrm{\Lambda }}_{ij}$が直交制約である。

CPMDの物理像

CPMDでは、電子軌道が基底状態の近傍で小さく振動しながら原子核に追随するように、周波数分離を作ることが意図される。すなわち、電子自由度の固有振動が原子核の振動より十分高く、かつ混合（エネルギー移動）が抑えられるように$\mu$と$\mathrm{\Delta }t$を選ぶ、という設計である。

この設計が崩れると、仮想電子自由度が原子核自由度へ熱を渡し、見かけの温度が変化するなどの問題が生じる。したがってCPMDは、BOMDよりも設定依存性が強く現れやすい一方で、SCF反復を減らしやすいという計算上の利点を持つ。

BOMDとCPMDの関係

BOMDは各ステップで電子を基底状態へ戻すのに対し、CPMDは電子を基底状態近傍に束縛して同時に動かす。両者は同じ断熱的地形の上の運動を狙うが、数値誤差の構造が異なるため、同じ物理量でも得やすさが変わることがある。

例えば、熱力学平均（構造分布や$g(r)$）が主目的であれば、時間相関の厳密さより分布の再現が重要になる。逆に、振動スペクトルや粘性のような時間相関が本質となる量では、温度制御や仮想自由度が相関へ与える影響を意識して条件を選ぶ必要がある。

有限温度DFTと占有数

金属系や高温系では、電子占有がフェルミ準位近傍で急峻であり、自己無撞着計算が不安定になりやすい。そのため、占有数を滑らかにする目的で電子温度（スメアリング）を導入することが多い。ここで重要なのは、これは原子核温度とは異なる道具であり、数値安定化のための自由エネルギー最小化として位置付く点である。

Merminの有限温度DFT

有限温度DFTでは、内部エネルギーではなく自由エネルギー

$$\mathcal{F}[\rho ]=E[\rho ]-TS[\{f_i\}]$$

を最小化する枠組みが導入される。$S$は占有数$\{f_i\}$に対応する電子エントロピーであり、フェルミ・ディラック分布に近い形で占有が決まる。したがって、計算が返す量が内部エネルギーか自由エネルギーか、あるいはエントロピー補正を含むかは、力や応力の一貫性にも関係する。

AIMDでは、原子核の運動方程式へ入る力が自由エネルギーの勾配になっているか、という点が本質である。占有数を変えると力の滑らかさが変わり、エネルギードリフトや相転移の見え方に影響することがあるため、温度条件に応じて占有の設定を一貫させる必要がある。

時間積分と時間刻み

原子核運動の時間積分には、エネルギー保存や位相体積保存の性質をなるべく壊しにくい方法が好まれる。速度Verlet法はその代表であり、加速度$\mathbf{a}(t)=\mathbf{F}(t)/M$を用いて

$$\mathbf{R}(t+\mathrm{\Delta }t)=\mathbf{R}(t)+\mathbf{v}(t)\mathrm{\Delta }t+\frac{1}{2}\mathbf{a}(t)\mathrm{\Delta }{t}^2$$$$\mathbf{v}(t+\mathrm{\Delta }t)=\mathbf{v}(t)+\frac{1}{2}[\mathbf{a}(t)+\mathbf{a}(t+\mathrm{\Delta }t)]\mathrm{\Delta }t$$

と更新する。

時間刻みを決める要因

時間刻み$\mathrm{\Delta }t$は、系の最も速い振動モード（高周波モード）をどこまで解像するかで上限が決まる。水素を含む結合伸縮、硬い共有結合、軽元素の局所振動は周波数が高く、$\mathrm{\Delta }t$を大きくすると数値的不安定や温度の異常上昇として現れやすい。

一方、$\mathrm{\Delta }t$を小さくすると、同じ物理時間を進めるのに必要なステップ数が増え、電子状態計算回数が増える。AIMDでは一ステップが高価であるため、$\mathrm{\Delta }t$の選択は計算可能時間と精度の釣り合いとして現れる。

質量スケーリングと拘束

水素の運動が時間刻みを制約する場合、質量スケーリング（重水素化に相当）や結合長拘束を用いる発想がある。ただし第一原理MDでは、電子状態が同じでも核運動の慣性が変わるため、得たい物理量（拡散、反応速度、振動）によっては解釈が変わる。何を保存し、何を近似するかを意識して導入する必要がある。

温度制御

温度一定条件を実現するには、系を熱浴と結合させてカノニカル分布へ導く方法を用いる。温度制御は単に温度を一定値に固定する操作ではなく、運動エネルギーのゆらぎを許した上で、長時間平均として所望の分布を再現する仕組みである。

Nosé–Hoover法

Noséの考え方は、拡張自由度を導入してカノニカル分布を生成する点にある。Hooverの形式では、運動量${\mathbf{p}}_I$に摩擦変数$\xi$が結合し、

$${\dot{\mathbf{p}}}_I={\mathbf{F}}_I-\xi {\mathbf{p}}_I$$$$\dot{\xi }=\frac{1}{Q}(\sum _I\frac{{\mathbf{p}}_I^2}{M_I}-g{k}_{B}T)$$

の形で表されることが多い。$Q$は熱浴の慣性に相当し、結合の強さと応答速度を支配する。

Nosé–Hooverは決定論的であり、適切に設定すればカノニカル分布を与えるが、系によってはエルゴード性が弱く、温度が十分混合しないことがある。そのため、鎖型（Nosé–Hoover chain）を用いて熱浴自由度を増やし、混合を改善する考え方が広く用いられる。

Langevin法と確率的熱浴

Langevin法は、摩擦とランダム力を加えることで温度を制御する。運動方程式は

$$M_I{\ddot{\mathbf{R}}}_I={\mathbf{F}}_I-\gamma M_I{\dot{\mathbf{R}}}_I+{\mathit{\eta }}_I(t)$$

であり、${\mathit{\eta }}_I(t)$は平均ゼロの確率過程で、ゆらぎ散逸関係により

$$⟨{\eta }_{I,\alpha }(t){\eta }_{J,\beta }(t^{\prime })⟩=2\gamma M_I{k}_{B}T{\delta }_{IJ}{\delta }_{\alpha \beta }\delta (t-t^{\prime })$$

を満たす。確率的熱浴はエルゴード性を改善しやすい一方、時間相関を直接議論する量（粘性、振動線幅など）では、摩擦項の影響を意識した解釈が必要になる。

温度制御と観測量の関係

構造分布（$g(r)$など）を得る目的では、分布が正しく生成されることが主眼であり、熱浴の種類は比較的選びやすい。これに対し、拡散係数や振動スペクトルのように時間相関が重要な量では、熱浴の影響を避けるために、平衡化は温度制御で行い、本計測はNVEで行うという考え方が用いられることがある。

ただしAIMDでは到達時間が短く、平衡化と計測を分けると統計が不足しやすい。目的量が分布か時間相関かを先に定め、その目的に対して温度制御の入れ方を選ぶことが、短い計算時間で意味ある結論を得るために重要になる。

圧力制御

圧力一定条件（NPT）を実現するには、セル体積やセル形状を動的に変化させ、目標圧力へ緩和させる必要がある。固体では等方体積変化だけでなく、せん断を含むセル形状変化が物性へ直結することがあるため、可変セルMDが重要になる。

Parrinello–Rahmanの考え方

Parrinello–Rahman法は、セル行列$\mathbf{h}$を自由度として導入し、${\mathbf{R}}_I=\mathbf{h}{\mathbf{s}}_I$（${\mathbf{s}}_I$は分率座標）と表して、$\mathbf{h}$の運動方程式を与える。セルが変形できることで、結晶構造の自発的な対称性変化や相転移が自然に表現される。

圧力制御では応力テンソルが直接関与するため、応力の収束が不十分だと、セルが不自然に振動したり、平均圧力がずれたりすることがある。AIMDで圧力制御を行う際は、力の収束だけでなく応力の収束を同時に満たす設定が必要になる。

周期境界条件と電荷・表面

凝縮系のAIMDでは周期境界条件が基本であり、有限サイズの超格子が無限系の近似になる。周期境界条件は計算効率を高める一方、系のサイズが小さいと周期像が相互作用し、拡散や欠陥の弾性場、界面の電場が人工的に強く見えることがある。

帯電系と補正

周期境界条件下で系が帯電すると、クーロン相互作用が発散するため、背景電荷で中和するなどの手続きが入ることが多い。これは物理的な境界条件を変更する操作でもあるため、帯電欠陥の形成エネルギーや電位整列などを議論する場合には、有限サイズ補正や電位基準の取り方が重要になる。

表面・スラブ系

表面を扱うスラブ模型では、真空層を入れ、法線方向に非対称な電荷分布が生じる場合がある。電気双極子補正を導入して人工的な電場を抑える手法が用いられることがあるが、AIMDでは温度ゆらぎで双極子が揺らぐため、平均電位や仕事関数を議論する際に時間平均の取り方が重要になる。

電子状態計算の設定がダイナミクスへ与える影響

AIMDの核は力であり、力の品質は電子状態計算の設定に支配される。ここでは、静的計算で許容できる粗さが、AIMDでは問題化する理由を、いくつかの軸で整理する。

基底（カットオフ）と力の滑らかさ

平面波基底ではカットオフエネルギーを上げるほど基底が増え、エネルギーと力が収束する。AIMDでは力が時系列として使われるため、基底が粗いと力のノイズが増え、温度ゆらぎやエネルギードリフトとして現れることがある。エネルギーの収束よりも力の収束が遅い場合がある点が、AIMDでは顕在化しやすい。

k点と金属系

金属系ではk点密度と占有数の扱いが重要である。k点が粗いとフェルミ面近傍の寄与が不正確になり、力や応力が振動したり、セル自由度の揺らぎが大きくなったりすることがある。スメアリングは数値安定化に有効だが、自由エネルギーを最小化する枠組みであることを踏まえ、温度条件や物性評価と整合する値に固定する必要がある。

スピンと磁性

磁性材料ではスピン分極を含めることが多いが、有限温度では局所磁気モーメントが時間的に揺らぐことがある。強磁性体であっても、短時間・小サイズのAIMDではスピン状態の初期条件や収束の道筋が系の平均磁化へ影響する場合がある。非共線スピンやスピン軌道相互作用まで入れると、計算コストだけでなく収束性も変化し、力のノイズ構造が変わるため、目的（構造・拡散・磁気ゆらぎのどれを主に見たいか）を先に定める必要がある。

DFT+Uやハイブリッド汎関数

強相関や局在電子を扱うためにDFT+Uやハイブリッド汎関数を用いる場合、エネルギー地形が変わり、力の計算コストも増える。AIMDでは各ステップでその重い計算を行うことになるため、到達時間が短くなる。したがって、どの相互作用が観測量に本質的かを見極め、必要な範囲で高精度手法を導入する設計が重要になる。

初期条件と平衡化

AIMDは初期条件に敏感であり、短い計算時間では初期構造の記憶が残りやすい。したがって、初期構造の準備と平衡化の設計が、結果の解釈と不可分になる。

初期構造の作り方

結晶では、まず0 Kの構造最適化で力が十分小さい状態を作り、その後に温度を導入することが多い。構造最適化が不十分なまま温度を上げると、単に残留応力が緩和しているだけの過程がダイナミクスとして見えてしまい、温度効果と区別しにくくなる。

液体・アモルファスでは、加熱して溶融させ、所望の温度へ冷却し、密度を整えるという手順が用いられる。ここでの加熱・冷却速度は実験と比べると非常に速く、得られる構造が準安定であることが多い。したがって、どの短距離秩序を再現したいか、どの温度域の性質を見たいかに応じて、複数の初期条件からの再現性を確認することが望ましい。

初速度と運動量の処理

初速度はMaxwell–Boltzmann分布から与えられることが多い。分布からのサンプリングは確率的であるため、短い計算では初速度のばらつきが観測量へ影響することがある。系全体の重心運動は内部自由度と独立であるため、全運動量をゼロにする処理が加えられることが多い。

平衡化時間と観測時間

温度制御を入れていても、分布がすぐに平衡へ達するとは限らない。拡散が遅い固体、ガラス状構造、界面の再配列などでは、原子が本質的に遅い自由度を持つため、数psから数十psでは十分に混合しないことがある。AIMDでは到達時間が限られるため、平衡化不足の可能性を前提に、観測量が時間窓に依存していないかを複数窓で確認する姿勢が重要になる。

観測量の計算

AIMDの出力は原子位置と速度の時系列であり、そこから構造・輸送・振動の情報を抽出する。以下では、AIMDの時間スケールでも比較的扱いやすい量から順に整理する。

動径分布関数と配位数

動径分布関数$g(r)$は、距離$r$にある原子対の確率密度を無相関系で規格化したものである。定義は

$$g(r)=\frac{1}{4\pi r^2{\rho }_0}⟨\sum _{i\ne j}\delta (r-|{\mathbf{r}}_i-{\mathbf{r}}_j|)⟩$$

であり、${\rho }_0$は平均数密度である。

$g(r)$の第一ピークは最近接殻の距離と分布幅を与える。配位数$N_c$は第一極小$r_{min}$まで積分して

$$N_c=4\pi {\rho }_0{\int }_0^{r_{min}}{r}^{2}g(r)dr$$

で定義され、液体・アモルファスの局所構造の比較に直接用いられる。

角度分布と局所多面体

共有結合性や局所多面体を議論するには、結合角分布$P(\theta )$が有効である。例えば3原子$i$-$j$-$k$の角度$\theta =\mathrm{\angle }ijk$を多数サンプルし、分布のピーク位置と幅を見ることで、四面体的秩序や平面三角形的秩序の混合を評価できる。

AIMDでは構造ゆらぎが自然に含まれるため、0 K最適化での単一角度ではなく、分布として角度を議論できる点が大きい。特にアモルファスでは、平均値より分布の非対称性や多峰性が構造の本質を示すことがある。

構造因子

散乱実験と対応させるには構造因子$S(\mathbf{q})$を計算する。数密度のフーリエ成分$\rho (\mathbf{q})=\sum _j{e}^{-i\mathbf{q}\cdot {\mathbf{r}}_j}$を用いると

$$S(\mathbf{q})=\frac{1}{N}⟨\rho (\mathbf{q})\rho (-\mathbf{q})⟩$$

で定義される。等方系では$q=|\mathbf{q}|$で平均し、$S(q)$として実験の散乱プロファイルと比較する。

拡散係数

輸送の基本量として拡散係数$D$がある。平均二乗変位（MSD）

$$\mathrm{MSD}(t)=⟨|{\mathbf{r}}_i(t)-{\mathbf{r}}_i(0){|}^2⟩$$

を計算し、十分長時間の拡散領域で

$$D=\underset{t\to \mathrm{\infty }}{lim}\frac{\mathrm{MSD}(t)}{6t}$$

で定義される。短時間では振動やケージングでMSDが直線にならないため、どの時間窓が拡散領域かを見極める必要がある。

AIMDでは到達時間が短いことが多く、$t\to \mathrm{\infty }$の極限へ達していない可能性が常にある。そこで、複数の原点時間（time origins）を平均して統計を稼ぎ、MSDの直線性が時間窓で安定しているかを見ることが重要になる。

振動スペクトル

振動状態密度の手掛かりは速度自己相関関数

$$C_v(t)=\frac{⟨\mathbf{v}(t)\cdot \mathbf{v}(0)⟩}{⟨\mathbf{v}(0)\cdot \mathbf{v}(0)⟩}$$

から得られる。フーリエ変換

$$I(\omega )={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{C}_v(t)e^{-i\omega t}dt$$

により周波数成分が得られ、固体ではフォノン的ピーク、液体では広がったバンドとして現れる。

温度制御が強いと速度相関が人工的に減衰することがあるため、振動線幅や寿命を議論する場合には、温度制御の入れ方と観測時間窓を慎重に選ぶ必要がある。ピーク位置（周波数）だけを議論するなら影響は比較的小さいが、線幅（散逸）まで含めると影響が顕在化しやすい。

粘性や熱伝導率

より進んだ輸送係数はGreen–Kubo関係で表される。例えばせん断粘性$\eta$は応力テンソルの自己相関から

$$\eta =\frac{V}{k_{B}T}{\int }_0^{\mathrm{\infty }}⟨{\sigma }_{\alpha \beta }(t){\sigma }_{\alpha \beta }(0)⟩dt(\alpha \ne \beta )$$

で与えられる。AIMDでは応力計算が高価で統計も不足しがちであるため、粘性のような長時間相関を要する量は難度が高いが、短時間でも相対比較や温度依存の傾向が得られることがある。

統計誤差と相関時間

AIMDで得られる時系列データは互いに独立ではなく、相関時間を持つ。したがって、単にサンプル数（ステップ数）が多いことと、有効な独立サンプル数が多いことは一致しない。

ある観測量$A(t)$の自己相関関数$C_A(t)$がゆっくり減衰するなら、有効サンプル数は小さくなる。統計誤差を見積もるには、ブロック平均によりブロック長を伸ばし、分散が収束する領域を探す方法がよく用いられる。AIMDでは到達時間が短いので、誤差が小さいという主張よりも、誤差を含めた上で結論が変わらない形の議論（定性的傾向、複数条件の差の頑健性）が重要になりやすい。

計算規模と到達時間

AIMDの計算コストは、電子状態計算の支配的処理（行列対角化など）により、原子数$N$に対して急速に増加することが多い。結果として、原子数を少し増やしただけで、同じ物理時間を進める計算時間が大きく伸びる。

時間方向は、$\mathrm{\Delta }t$がfs程度で、数万ステップが現実的な上限として現れやすい。したがって、psから数十psの現象は扱いやすいが、遅い拡散や相分離、欠陥の希少事象はAIMD単独では捉えにくいことが多い。この制約をどう補うかが、AIMDの使い方の設計そのものになる。

機械学習ポテンシャルとの併用

AIMDの結果（エネルギーと力、応力）を教師データとして機械学習ポテンシャル（MLIP）を構築し、より大規模・長時間のMDへ拡張する流れが強い。AIMDはその際の基準器として働き、学習モデルが再現すべきエネルギー地形を与える。

さらに、MDの進行中に予測不確かさを評価し、不確かさが大きい配置だけ第一原理計算で追加学習するオンザフライ学習も広く用いられる。この方向では、AIMDの短時間性が弱点でなく、学習データの質を保証する長所として機能する。

手法の比較

手法	力の起源	主に得やすい情報	到達しやすい系サイズ・時間	設定が効きやすい点
古典MD	経験的ポテンシャル	大規模統計、長時間輸送	大・長	ポテンシャルの適用範囲
BOMD（AIMD）	DFTなどの基底状態エネルギー勾配	反応、相転移、有限温度構造	小〜中・短〜中	電子収束、基底、占有数
CPMD（AIMD）	拡張ラグランジアンの連成	構造ゆらぎ、液体構造	小〜中・短〜中	仮想質量、時間刻み、熱浴
非断熱第一原理MD	電子の時間発展を含む	励起・電荷移動・散逸	小・短	電子状態の混合、遷移の扱い
MLIP-MD	学習した近似ポテンシャル	大規模・長時間 + 反応も場合により	大・長	学習データの網羅性

BOMDとCPMDは、同じ第一原理のエネルギー地形を目指しつつ、電子自由度の扱いが異なる。どちらを選ぶかは、時間相関をどれだけ厳密に扱いたいか、電子収束コストと設定依存性のどちらを許容するか、という観点で整理できる。

まとめと展望

第一原理分子動力学法は、原子核の運動方程式に必要な力を第一原理電子状態計算から逐次与え、経験的ポテンシャルに頼らずに有限温度ダイナミクスを扱う方法である。BOMDは電子を逐次収束させて地形上の運動を実現し、CPMDは拡張自由度で電子を基底状態近傍へ束縛して連成運動として進める。

今後は、有限温度DFTや占有数の扱いを含む力・応力の一貫性を高めつつ、機械学習ポテンシャルとの併用で系サイズと時間の制約を乗り越える方向がさらに進むと考えられる。加えて、非断熱効果や散逸をより直接に取り込む第一原理ダイナミクスが、励起下物性や極限条件材料へ接続され、AIMDの適用範囲は広がり続ける見込みである。

参考文献

• S. Nosé, A unified formulation of the constant temperature molecular dynamics methods, J. Chem. Phys. 81, 511 (1984) https://pubs.aip.org/aip/jcp/article/81/1/511/607222/A-unified-formulation-of-the-constant-temperature
• W. G. Hoover, Canonical dynamics: Equilibrium phase-space distributions, Phys. Rev. A 31, 1695 (1985) https://link.aps.org/doi/10.1103/PhysRevA.31.1695
• M. Parrinello and A. Rahman, Crystal Structure and Pair Potentials: A Molecular-Dynamics Study, Phys. Rev. Lett. 45, 1196 (1980) https://link.aps.org/doi/10.1103/PhysRevLett.45.1196
• N. D. Mermin, Thermal Properties of the Inhomogeneous Electron Gas, Phys. Rev. 137, A1441 (1965) https://link.aps.org/doi/10.1103/PhysRev.137.A1441
• VASP Wiki, Molecular-dynamics calculations（温度制御・圧力制御の実装整理） https://vasp.at/wiki/Molecular-dynamics_calculations
• VASP Wiki, Nosé-Hoover thermostat https://www.vasp.at/wiki/Nosé-Hoover_thermostat
• VASP Wiki, Nosé-Hoover chain thermostat https://www.vasp.at/wiki/Nosé-Hoover_chain_thermostat