第1章マテリアルに応用される計算科学と求まる物理量
本書は「原子がどのように並んでいるか」という情報だけを出発点にして，その物質がどんな性質を示すのかを紙と計算機の上で言い当てる，という営みを扱う．第1章はその全体像を眺める章である．計算科学と一口に言っても，原子1個の大きさ（$10^{-10}\ \mathrm{m}$）を扱うものから橋やビル（$10^{2}\ \mathrm{m}$）を扱うものまで，対象とするスケールごとにまったく別の手法が使われている．本章ではまずスケールの地図を描き，そのうえで原子スケールの中心的な手法である第一原理計算・分子動力学計算・モンテカルロ法を紹介する．続いて，結晶構造を入力すると何が出てくるのかを「型1：全エネルギー」「型2：電子バンド」「型3：フォノンバンド」の三つに整理し，実際の材料研究の例で確かめる．数式が本格的に登場するのは第2章からだが，本章でも単振動の運動方程式だけは完全に解いてみせる．そこから，振動数が虚数になるという一見不可解な状況が，じつは「その構造は不安定だ」という強力な予言になっていることを理解する．
- 唯物論（materialism）と要素還元主義（reductionism）という立場，および「Input＝結晶構造，Output＝物性」という計算科学の基本構図
- 材料現象の長さスケールの階層（原子・メゾ・マクロ）と，各スケールで用いられる計算手法の対応
- 第一原理計算・分子動力学計算・モンテカルロ法・フェーズフィールド法・有限要素法が，それぞれ何を計算する道具なのか
- 型1：全エネルギー $E_{\mathrm{tot}}$ から，結晶多型の安定性・新物質探索・表面エネルギー・分子の解離が議論できること
- 型2：電子バンドから，金属か絶縁体か，結合性軌道か反結合性軌道か（COHP・COOP・COBI の符号）が読み取れること
- 型3：フォノンバンドから振動数 $\omega=\sqrt{k/m}$ が得られ，$k\lt 0$ のとき $\omega=i\sqrt{\abs{k}/m}$ という虚数振動数になること，その物理的意味
- 単振動の運動方程式 $m\ddot{x}=-kx$ と $m\ddot{x}=+kx$ を初期条件つきで完全に解き，$\sin$ 解と $\sinh$ 解を得ること
- 本書全14章がどのようにつながって「分子の量子化学計算」から「固体の電子状態」へ到達するのか
1.1 なぜ計算科学を学ぶのか — 原子の並び方から物性を演繹する
1.1.1 唯物論と要素還元主義
本書を読み進めるあいだ，次の思想を頭の片隅に置いておいてほしい．
言い換えれば「どのような元素を，どのように並べたら，どんな性質が出てくるのか」という問いである．この立場を支えているのが二つのキーワード，唯物論（materialism）と要素還元主義（reductionism）である．唯物論とは，物質世界で起きる現象は物質そのものの構成と運動によって説明されるという立場であり，要素還元主義とは，複雑な系のふるまいを，それを構成する要素とその相互作用にまで分解して理解しようとする立場である．
材料科学に翻訳すればこうなる．ある物質が硬いか柔らかいか，電気を通すか通さないか，赤く光るか光らないか——これらはすべて，その物質を構成する原子核と電子，およびそれらのあいだに働く相互作用にまで還元して説明できるはずだ，ということである．そして原子核と電子の相互作用を支配する法則は，すでに20世紀前半に量子力学として書き下されている．つまり原理的には，結晶構造さえ与えれば，あとは方程式を解くだけで物性が出てくる．
なぜ？：どうして「原子の並び方」だけを入力にしてよいのか
物質を構成する部品は原子核と電子である．原子核の種類（＝元素）と位置が決まれば，電子が感じる静電ポテンシャルが完全に決まる．すると電子の従うSchrödinger方程式が確定し，原理的には電子の状態がすべて決まってしまう．すなわち，「どの元素の原子核がどこにあるか」という情報だけが，物質を指定するために必要な唯一の入力である．これが計算科学の出発点であり，なぜ入力が「結晶構造（分子なら分子構造）」ひとつで済むのかの理由である．第8章で学ぶBorn–Oppenheimer近似は，この「原子核の位置を先に固定してしまう」という手続きを正当化する近似である．
1.1.2 Input は結晶構造，Output は物性
この思想を計算の言葉に直すと，次の単純な図式になる．
ここで重要なのは，入力に実験値を一切使っていないという点である．使うのは元素の種類（原子番号）と原子の位置，そして物理定数（電子質量 $m$，電気素量 $e_0$，真空の誘電率 $\varepsilon_0$，Planck定数 $h$）だけである．実験で測った量をパラメータとして与えないので，この種の計算を第一原理計算（first-principles calculations，あるいは ab initio calculations）と呼ぶ．「第一原理」とは「基本法則だけから出発する」という意味である．
1.1.3 導入例：ダイヤモンドとグラファイト
この図式の威力を実感するのに，これ以上ないほど適した例がある．ダイヤモンドとグラファイトである．両者はまったく同じ炭素だけからできている．含まれる元素は $\mathrm{C}$ ただ一種類で，組成式はどちらも $\mathrm{C}$ である．違うのは並び方だけだ．
- ダイヤモンド：各炭素原子が4個の隣接炭素と正四面体状に $\mathrm{sp}^3$ 混成軌道で結合し，三次元的な共有結合網をつくる．結合長は $1.54\times10^{-10}\ \mathrm{m}$．
- グラファイト：各炭素原子が3個の隣接炭素と平面六角形状に $\mathrm{sp}^2$ 混成軌道で結合し，シート（グラフェン層）をつくる．層内の結合長は $1.42\times10^{-10}\ \mathrm{m}$ と短く強いが，層と層のあいだは $3.35\times10^{-10}\ \mathrm{m}$ 離れており，弱いvan der Waals力（第13章）だけで積み重なっている．
この「並び方の違い」だけで，物性は劇的に変わる．ダイヤモンドは知られている物質のなかで最も硬く，無色透明で，電気をほとんど流さない絶縁体である．一方グラファイトは鉛筆の芯として紙の上を滑るほど柔らかく，黒色不透明で，層に平行な方向にはよく電気を流す．同素体（allotrope）という言葉ではこの落差を説明したことにならない．なぜそうなるのかを，原子の並び方から演繹して説明したい．それが計算科学の仕事である．
注意：グラファイトは「金属」か「半金属」か
講義では手短に「グラファイトは金属」と言うことがあるが，より正確には半金属（semimetal）である．半金属とは，バンドギャップが開いていない（あるいは負の値になっている）ものの，Fermi準位における状態密度が普通の金属に比べて桁違いに小さい物質を指す．グラファイトの価電子帯と伝導帯の重なりはたかだか $0.04\ \mathrm{eV}$ 程度で，キャリア密度は銅のような典型的金属より4〜5桁小さい．また，層に平行な方向の電気伝導度は層に垂直な方向の $10^{3}$〜$10^{4}$ 倍にもなり，強い異方性をもつ．本書では「ギャップが開いていない＝電気を流しうる」という点を強調するときに金属的と呼び，厳密な分類が必要な場面では半金属と書く．
1.1.4 理論・実験・計算の三本柱
ここで，計算科学が科学の営みのなかでどこに位置するのかを確認しておきたい．長らく自然科学は理論と実験の二本柱で進んできた．理論は紙と鉛筆で解ける問題を扱い，実験は自然に直接問いかける．しかし，量子力学の方程式は水素原子のようなごく限られた場合を除いて解析的に解けない（第2章）．一方で実験には，試料をつくるコスト，危険物や希少元素の扱い，そして「なぜそうなるのか」が直接は見えないという限界がある．
計算は，この二つの隙間を埋める第三の柱である．解析的に解けない方程式を数値的に解き，実験では取り出せない量（たとえば結晶中のある一つの化学結合がどれだけエネルギーを下げているか）を取り出す．さらに，まだ誰も合成していない物質について「つくったらどうなるか」を先に答えることさえできる．1.5節で紹介する新物質探索は，まさにその実例である．
1.2 材料の現象はスケールで考える
1.2.1 長さのスケールと学問分野
材料に関わる現象を考えるとき，まずスケールを考えるという習慣を身につけてほしい．同じ「材料」という言葉を使っていても，原子1個の中を問題にしているのか，結晶粒の集合体を問題にしているのか，橋やタービンブレードという構造物を問題にしているのかで，必要な理論も手法もまったく異なるからである．
長さのスケールを $10$ の冪でたどると，次のような階層が見えてくる．そして階層ごとに，それを主戦場とする学問分野が対応している．
数学ノート：スケールの数字の読み方
講義で用いた模式図では原子の位置に $10^{-12}\ \mathrm{m}$ と書かれているが，これは「原子の内部構造まで含めた領域」を大まかに指した目安である．数値をきちんと押さえておこう．原子核の直径は $10^{-15}$〜$10^{-14}\ \mathrm{m}$（フェムトメートル），原子の直径は $10^{-10}\ \mathrm{m}$ 程度である．この $10^{-10}\ \mathrm{m}$ という長さにはオングストローム（$1\ \text{Å}=10^{-10}\ \mathrm{m}$）という専用の単位があり，材料科学ではきわめて頻繁に使われる．水素原子の Bohr半径（第2章）$a_0 = \dfrac{\varepsilon_0 h^2}{\pi m e_0^2} = 0.529\times10^{-10}\ \mathrm{m} = 0.529\ \text{Å}$ が，原子スケールの自然な物差しになる．原子核と原子で5桁も違うのに，化学と材料科学が原子核の中身をいっさい気にしないで済むのは，この隔たりのおかげである．
1.2.2 スケールと計算手法の対応
スケールが変われば，そもそも「何を粒子とみなすか」が変わる．原子スケールでは電子と原子核を扱う．メゾスケールでは個々の原子を追うのをあきらめ，「この領域はどの相か」「分極はどちらを向いているか」といった秩序変数（order parameter）の場を扱う．マクロスケールでは物質を連続体とみなし，応力とひずみを扱う．この読み替えに応じて，使う手法も次のように対応する．
| スケール | 手法 | 扱う対象 | 典型的な規模 | 主な出力 |
|---|---|---|---|---|
| 原子 $10^{-10}$〜$10^{-8}\ \mathrm{m}$ | 第一原理計算（DFT） | 電子密度・原子核位置 | $10^{2}$〜$10^{3}$ 原子，0 K | 全エネルギー，電子バンド，フォノン |
| 第一原理分子動力学 | 電子＋原子核の運動 | $10^{2}$ 原子，$10\ \mathrm{ps}$ | 有限温度の構造，拡散，相転移 | |
| 古典分子動力学 | 原子（点粒子） | $10^{6}$〜$10^{9}$ 原子，$\mathrm{ns}$〜$\mu\mathrm{s}$ | 熱膨張，粘性，力学特性 | |
| モンテカルロ法 | 配置・スピン・占有 | $10^{4}$〜$10^{6}$ サイト | 相図，規則–不規則転移，固溶体 | |
| メゾ $10^{-8}$〜$10^{-4}\ \mathrm{m}$ | フェーズフィールド法 | 秩序変数の場 | $\mu\mathrm{m}$ 領域，$\mathrm{s}$ | 組織形成，ドメイン運動，析出 |
| 転位動力学 | 転位線 | $\mu\mathrm{m}$ 領域 | 加工硬化，降伏 | |
| マクロ $10^{-4}$〜$10^{0}\ \mathrm{m}$ | 有限要素法 | 連続体（要素の集合） | $\mathrm{mm}$〜$\mathrm{m}$，$10^{6}$ 要素 | 応力分布，変形，振動，破壊 |
例題1.1 原子スケールの計算が「小さい」ことを実感する
ケイ素（Si）の密度は $2.33\ \mathrm{g/cm^3}$，モル質量は $28.09\ \mathrm{g/mol}$ である．
(1) $1\ \mathrm{cm^3}$ の Si に含まれる原子数を求めよ．
(2) 第一原理計算で扱えるのがせいぜい $10^{3}$ 原子だとすると，それは一辺何 nm の立方体に相当するか．
解答 (1) 単位体積あたりの物質量は $2.33/28.09 = 8.29\times10^{-2}\ \mathrm{mol/cm^3}$ である．Avogadro定数 $N_{\mathrm{A}}=6.022\times10^{23}\ \mathrm{mol^{-1}}$ を掛けて
$$ n = 8.29\times10^{-2}\times 6.022\times10^{23} = 4.99\times10^{22}\ \mathrm{cm^{-3}} \simeq 5.0\times10^{22}\ \text{原子}/\mathrm{cm^3}. $$(2) $10^{3}$ 原子が占める体積は $V = 10^{3}/(4.99\times10^{22}) = 2.00\times10^{-20}\ \mathrm{cm^3}$．一辺は
$$ L = V^{1/3} = (2.00\times10^{-20})^{1/3}\ \mathrm{cm} = 2.71\times10^{-7}\ \mathrm{cm} = 2.71\ \mathrm{nm}. $$つまり第一原理計算が直接扱えるのは，一辺 $3\ \mathrm{nm}$ に満たない極小の立方体にすぎない．$1\ \mathrm{cm^3}$ の試料を丸ごと計算しようとすれば原子数は $10^{19}$ 倍になり，計算量が原子数の3乗で増えることを考えると $10^{57}$ 倍の手間になる．どれだけ計算機が速くなっても届く数字ではない．だからこそスケールを分け，各スケールに適した手法を使い分ける必要がある．表1.1と図1.3の意味はここにある．
数学ノート：ちょっとした余談 — 雑誌 Physical Review
本章にはこれから多くの論文が引用されるが，そのかなりの割合が Physical Review という雑誌に載っている．物理学で最も権威のある雑誌群のひとつで，本書に関係する記念碑的な論文がここに集まっている．たとえば——
- E. Schrödinger, "An Undulatory Theory of the Mechanics of Atoms and Molecules", Phys. Rev. 28, 1049 (1926)：本書第2章以降の主役であるSchrödinger方程式．1933年ノーベル物理学賞．
- J. Bardeen, L. N. Cooper, J. R. Schrieffer, "Theory of Superconductivity", Phys. Rev. 108, 1175 (1957)：超伝導のBCS理論．1972年ノーベル物理学賞．
- W. Kohn, L. J. Sham, "Self-Consistent Equations Including Exchange and Correlation Effects", Phys. Rev. 140, A1133 (1965)：本章1.3節で紹介する第一原理計算の土台．1998年ノーベル化学賞．
- B. P. Abbott et al. (LIGO), "Observation of Gravitational Waves from a Binary Black Hole Merger", Phys. Rev. Lett. 116, 061102 (2016)：重力波の初観測．2017年ノーベル物理学賞．
いずれも「発表されたときはただの1本の論文だった」ものである．教科書に載っている式にも，それが世に出た最初の日がある，ということを心に留めておくとよい．
1.3 原子スケールの計算手法
ここからは，本書が主戦場とする原子スケールの手法を，代表的な三つ——第一原理計算，分子動力学計算，モンテカルロ法——に絞って紹介する．「計算科学」と一口に言っても本当に色々あるのだ，ということをまず知ってほしい．
1.3.1 密度汎関数理論に基づく第一原理計算
定義：第一原理計算（first-principles calculations）
実験値をパラメータとして用いず，量子力学の基本法則と物理定数のみから物質の性質を計算する手法の総称．現在の材料科学で第一原理計算といえば，ほとんどの場合密度汎関数理論（density functional theory, DFT）に基づく計算を指す．密度汎関数理論とは，運動エネルギーやポテンシャルエネルギーといった量を，多電子波動関数の代わりに電子密度 $\rho(\bm{r})$ の関数（正確には汎関数）として表現する理論である．Input は結晶構造，Output は種々の物性であり，主として絶対零度における静的な物性を与える．
「汎関数」という言葉に身構える必要はない．関数が数を数に対応させるものだとすれば，汎関数は関数を数に対応させる規則のことである．電子密度 $\rho(\bm{r})$ という関数を入れると全エネルギーという1個の数が出てくる，という対応が $E[\rho]$ という汎関数である．この理論の基礎は前掲のKohn–Shamの論文（1965年）に置かれ，W. Kohn は1998年にノーベル化学賞を受賞した．理論の中身は本書の範囲を超えるが，その入口にあたる変分原理（第3章），Slater行列式と交換エネルギー（第7・9章），Born–Oppenheimer近似（第8章）は，本書で一つずつ丁寧に扱う．
第一原理計算から具体的にどんな物理量が得られるのか，代表例を挙げておく．
| 物理量 | 何がわかるか | 関連 |
|---|---|---|
| 全エネルギー $E_{\mathrm{tot}}$ | 結晶多型のうちどれが安定か | 1.5節，第9章 |
| 生成エネルギー・化学ポテンシャル依存の相図 | どの組成・雰囲気でどの相が現れるか | 1.5節 |
| 欠陥形成エネルギー | 空孔・不純物がどれだけ入りやすいか，電荷状態はどうか | 1.5節 |
| 光吸収係数 $\alpha(h\nu)$ | どの波長の光を吸収・発光するか | 1.5節 |
| 電子バンド構造 | 金属か絶縁体か，バンドギャップの大きさ | 1.6節，第14章 |
| COHP・COOP・COBI | 個々の化学結合が結合性か反結合性か | 1.6節，付録1.A |
| フォノンバンド | 構造が動的に安定か，振動の振動数はいくらか | 1.7節 |
| 熱膨張係数・体積弾性率 | 温めるとどれだけ伸びるか，圧すとどれだけ縮むか | 1.7節 |
表1.2の下2つ——熱膨張係数や体積弾性率——は，一見すると「有限温度の量」「力学的な量」であって，0 K の静的計算からは出てこないように思える．ところが，全エネルギーを体積の関数として計算しておけば体積弾性率 $B = V\,\partial^2 E/\partial V^2$ が求まり，フォノンの振動数の体積依存性まで計算すれば熱膨張係数が求まる．0 K の全エネルギーと振動状態がわかれば，有限温度の熱力学量にまで手が届くのである．実例としては，MgO や CaO の体積熱膨張係数と体積弾性率を第一原理計算から求めた研究がある〔J. Phys. Chem. C 128, 525 (2024)〕．また欠陥形成エネルギーの計算については ZnO を対象とした古典的な研究〔Phys. Rev. B 77, 245202 (2008)〕が知られており，酸素が豊富な雰囲気（O-rich）と乏しい雰囲気（O-poor）で，どの欠陥がどの電荷状態で安定になるかが定量的に議論されている．
1.3.2 分子動力学計算
定義：分子動力学計算（molecular dynamics calculations, MD）
各原子について Newton の運動方程式 $M_I \ddot{\bm{R}}_I = \bm{F}_I$ を数値的に時間積分し，原子の軌跡を追跡することで，主として動的な（有限温度における）物性を計算する手法．原子に働く力 $\bm{F}_I$ をどうやって得るかによって，次の二種類に分かれる．
- 古典分子動力学計算（classical MD）：あらかじめ関数形を決めた原子間ポテンシャル模型（Lennard-Jones型，埋め込み原子法，Buckingham型など）から力を計算する．1ステップが軽いので，$10^{6}$ 個以上の原子を $\mathrm{ns}$〜$\mu\mathrm{s}$ 追跡できる．
- 第一原理分子動力学計算（ab initio MD, AIMD）：各時刻ごとに電子状態を第一原理計算で解き直し，そこから力を求める．模型を仮定しないので化学反応や結合の組み替えも扱えるが，計算量は桁違いに大きく，$10^{2}$ 原子・$10\ \mathrm{ps}$ 程度が実際的な上限である．
古典MDは速いが，原子間ポテンシャル模型の出来にすべてが左右される．第一原理MDは信頼できるが遅い．この二律背反を大きく変えつつあるのが機械学習ポテンシャル（machine-learning interatomic potential, MLIP）である．あらかじめ第一原理計算で多数の構造についてエネルギーと力を計算しておき，その対応関係をニューラルネットワークなどで学習させる．学習が済めば，第一原理計算に近い精度の力を古典MD並みの速度で評価できる．近年のMD計算がすさまじい勢いで進歩しているのは，主にこの技術による．
MDの成果としては，構造相転移の理論予測〔Phys. Rev. Lett. 122, 225701 (2019)〕や，格子定数の温度依存性の予測がある．たとえば $\mathrm{KNbSi_2O_7}$ では，まず第一原理計算と機械学習ポテンシャルによって，$a$ 軸と $b$ 軸が温度を上げるほど縮むこと——すなわち二軸性の負熱膨張（negative thermal expansion, NTE）——が理論的に予測され，そのあとで合成と回折実験によって確かめられた．実験で測られた格子定数の温度変化は，計算があらかじめ示していた傾き・曲率をよく再現しており，計算の精度がきわめて高いことが裏づけられた〔J. Phys. Chem. C 129, 9509 (2025)〕．
この順序——計算が先，実験があと——には，材料探索の道具としての大きな意味がある．ほとんどの物質は温めれば膨らむのだから，負熱膨張はまれな例外であり，新しい負熱膨張材料はこれまで「作って測ってみたら縮んでいた」という形でしか見つからなかった．どの物質がどの方向に縮むかを机上で言い当てることは，長らく不可能とさえ思われていたのである．ところが，フォノンの体積依存性を第一原理計算で押さえ，機械学習ポテンシャルによって有限温度のMDにまで手を伸ばせば，合成に取りかかる前に「この物質はこの軸方向に縮む」と言い切れてしまう．新規の負熱膨張材料を理論から予測できるようになったことは，計算科学が実験の後追いではなく実験を導く道具になりつつあることを示す，わかりやすい実例である．「温度を上げると格子がどう伸び縮みするか」は実験でしか測れない量に見えるが，原子の運動をそのまま追えば計算でも——しかも実験より先に——出てくるのである．
例題1.2 MD の時間刻みと到達できる時間
典型的な光学フォノンの振動数を $\nu = 12\ \mathrm{THz}$ とする．
(1) 振動の周期を求めよ．
(2) MD の時間刻み（タイムステップ）を $\Delta t = 1\ \mathrm{fs}$ とすると，1周期は何ステップで表現されるか．
(3) $10^{6}$ ステップ計算したとき，追跡できる実時間はどれだけか．
解答 (1) 周期は振動数の逆数だから
$$ T = \frac{1}{\nu} = \frac{1}{12\times10^{12}\ \mathrm{s^{-1}}} = 8.33\times10^{-14}\ \mathrm{s} = 83.3\ \mathrm{fs}. $$(2) $T/\Delta t = 83.3\ \mathrm{fs} / 1\ \mathrm{fs} = 83$ ステップ．1周期を80点あまりでサンプルすることになり，数値積分の精度としては十分である．逆に $\Delta t$ を $10\ \mathrm{fs}$ などに大きくすると1周期が8点しかとれず，エネルギーが保存しなくなる．タイムステップは系に存在する最も速い振動の周期の $1/20$ 以下にとる，というのが実務上の目安である．
(3) $10^{6}\times 1\ \mathrm{fs} = 10^{6}\times10^{-15}\ \mathrm{s} = 10^{-9}\ \mathrm{s} = 1\ \mathrm{ns}$．
100万ステップという大変な計算をしても，追えるのはたった1ナノ秒である．原子スケールの計算がいかに「短い時間しか見られない」かがわかる．秒の単位で起きるゆっくりした組織変化を追うには，1.4節のフェーズフィールド法のように，原子を追うことを潔くあきらめた手法が必要になる．
1.3.3 モンテカルロ法
定義：モンテカルロ法（Monte Carlo method, MC）
熱平衡状態における物理量の平均値を，確率的なサンプリングによって計算する手法．ある状態 $s$ の出現確率が Boltzmann 分布 $P(s)\propto \exp[-E(s)/k_{\mathrm{B}}T]$ に従うように乱数を用いて状態を次々と生成し，生成された状態の集合について物理量を平均する．原子間ポテンシャル模型を用いる古典モンテカルロ法と，多体波動関数そのものを確率的にサンプルする量子モンテカルロ法がある．
MD と MC は，どちらも「熱平衡状態での平均値」を求める手法でありながら，平均のとり方が根本的に違う．ここで，この違いを一発で覚えられる比喩を紹介しよう．
MD は，原子が実際に動いていく様子を1コマ1コマ，時間順につないだ動画である．コマとコマは因果でつながっており，前のコマの位置と速度から次のコマが決まる．だから「どれくらいの速さで拡散するか」「どのくらいの時間で相転移するか」といった時間に関する情報が得られる．
一方 MC は，Boltzmann 分布に従うようにばらばらに描かれた絵を集めた漫画である．1ページ目と2ページ目のあいだに「実際の時間」は流れていない．各コマは独立に，正しい確率で選ばれた状態のスナップショットである．したがって時間の情報は得られないが，そのぶん状態空間を効率よく渡り歩ける．とくに，原子がどのサイトを占めるかという配置の問題——たとえば固溶体 $\mathrm{A}_{1-x}\mathrm{B}_x$ で A と B がどう混ざるか，酸素空孔がどこにどれだけ入るか——では，原子がサイトからサイトへ跳び移る事象は $\mathrm{ns}$ の MD ではまず起こらないため，MD では永遠に平衡に達しない．MC ならば「A と B を入れ替える」という試行を確率的に受理・棄却するだけでよく，一気に平衡配置に到達できる．
物理的意味：なぜ二つの平均が同じ答えを与えるのか
「時間平均」と「確率平均（集団平均）」が一致するという主張は，統計力学のエルゴード仮説（ergodic hypothesis）と呼ばれる．十分長い時間ほうっておけば，系はエネルギー的に許されるあらゆる状態を，Boltzmann 分布の重みどおりに訪れる，という仮説である．この仮説が成り立つ限り，MD の答えと MC の答えは一致しなければならない．逆に言えば，両者が食い違うときは，MD の計算時間が短すぎて系が状態空間を渡りきっていないことを疑うべきである．固溶体の配置問題で MD が使えないのは，まさにこの理由による．
MC の実例としては，$\mathrm{ZrH_2}$ の格子定数の温度依存性〔Phys. Rev. B 88, 214111 (2013)〕や，$\mathrm{ZrO}_x$ における酸素空孔濃度別の相図〔Phys. Rev. B 88, 094108 (2013)〕がある．後者は，酸素の組成比を横軸，温度を縦軸にとった相図の上に $\mathrm{ZrO_{1/6}}$，$\mathrm{ZrO_{1/3}}$，$\mathrm{ZrO_{1/2}}$ といった秩序相の領域が描かれるもので，「どの酸素サイトが埋まるか」という純粋に配置の問題を扱っている．MD ではまず到達できない結果である．
1.4 メゾ・マクロスケールの計算手法
1.4.1 フェーズフィールド法
定義：フェーズフィールド法（phase-field method）
自由エネルギーのモデルから，物質・材料の相転移や材料組織の変化を計算する手法．個々の原子を追う代わりに，位置 $\bm{r}$ と時刻 $t$ の関数である秩序変数 $\eta(\bm{r},t)$（相の種類，電気分極，組成など）を導入し，全自由エネルギー $F[\eta]$ が減少する向きに $\eta$ を時間発展させる．典型的には保存量でない秩序変数に対して Allen–Cahn 方程式 $\partial \eta/\partial t = -L\,\delta F/\delta \eta$，組成のような保存量に対して Cahn–Hilliard 方程式を用いる．
「界面がどこにあるか」を陽に追跡しなくてよいのがこの方法の強みである．$\eta$ が $0$ から $1$ へなめらかに変化する薄い領域が自動的に界面になる．おかげで，複雑に入り組んだ組織の生成・消滅・合体をそのまま計算できる．
材料研究での代表例を二つ挙げる．ひとつは強誘電体の電気分極反転である．強誘電体薄膜に探針で電圧をかけると，分極が反転した領域（分極ドメイン）が核形成し，ドメイン壁が動いて広がっていく．この一連の過程をフェーズフィールド法で追跡し，透過電子顕微鏡で観察された分極マップと直接比較した研究がある〔Acta Materialia 112, 285 (2016)〕．もうひとつは合金の組織変化で，Ni–Al 系の相図上の各点（Ni 22% から 92% まで）について，時間発展の末にどんな組織——析出物が球状に散らばるのか，板状に成長するのか，迷路のような二相共存組織になるのか——が現れるかを網羅的に計算した研究が知られている〔Nature Communications 10, 3451 (2019)〕．
なぜ？：原子を追うのをやめてよい理由
分極ドメインの大きさは $10^{-7}$〜$10^{-6}\ \mathrm{m}$，動く時間は $10^{-6}$〜$10^{0}\ \mathrm{s}$ である．例題1.1・1.2で見たとおり，原子スケールの計算は $10^{-9}\ \mathrm{m}$・$10^{-9}\ \mathrm{s}$ の世界にとどまる．桁が3〜9も違うのだから，同じ土俵で戦うのは不可能である．そこで「この領域の分極はどちらを向いているか」という粗視化した情報だけを残し，原子の細かい振動は自由エネルギーの中に押し込めてしまう．捨てた情報の代償は，自由エネルギー模型のパラメータ（分極の何乗の項の係数，界面エネルギーなど）を外から与えなければならない点にある．そのパラメータを第一原理計算から供給するのが，図1.3で示したマルチスケール計算の考え方である．
1.4.2 有限要素法
定義：有限要素法（finite-element method, FEM）
物体や構造物を小さな要素（三角形・四面体・六面体など）に分割し，要素と要素をばねで繋いだ連続体とみなして，その性質を数値化して計算することで全体の挙動を解析する手法．各要素の内部では変位が簡単な関数（多くは1次または2次の多項式）で表されると仮定し，全体の弾性エネルギーが最小になる条件から連立方程式を立てて解く．
有限要素法は，構造物の垂直応力・せん断応力の分布，変形，固有振動，破壊の起点の予測などに用いられる．橋梁の設計は代表的な応用先で，たとえば連続剛性橋の桁に生じる応力分布と，実際に生じた損傷（下フランジ接合部の損傷，負曲げ領域の舗装ひび割れ）との対応を解析した研究がある〔arXiv:2501.03093〕．
物理的意味：タコマ橋の崩落が教えたこと
1940年，米国ワシントン州のタコマナローズ橋は，開通からわずか4か月後に，風速 $19\ \mathrm{m/s}$ という決して強くない風のもとで大きくねじれ振動し，崩落した．橋は静的な荷重に対しては十分な強度をもっていた．崩落の原因は，風と橋のねじれ変形が結びついて振動が自励的に成長する現象（空力弾性フラッター）であり，当時の設計はこれを考慮していなかった．
ここで注目したいのは，この事故の本質が「構造物の振動モードと，その振動が成長するか減衰するか」にあった点である．これは1.7節で扱う，結晶の振動モードとその安定性の問題と，数学的にはまったく同じ構造をしている．振動が時間とともに減衰するか成長するかは，運動方程式に現れる係数の符号で決まる．スケールが $10^{-10}\ \mathrm{m}$ の結晶であれ $10^{3}\ \mathrm{m}$ の吊り橋であれ，「平衡点のまわりで揺らしたときに戻ってくるか，離れていくか」という問いの形は変わらないのである．
ここから1.7節までは，原子スケールの第一原理計算に話を絞り，結晶構造を入力すると何が出てくるのかを三つの「型」に整理して，実際の研究例とともに見ていく．型1は全エネルギー，型2は電子バンド，型3はフォノンバンドである．この三つを押さえれば，第一原理計算の論文を読むときの見取り図ができあがる．
1.5 型1：結晶構造から全エネルギーを得る
1.5.1 全エネルギーとは何か
第一原理計算がまず出力するのは，与えられた結晶構造の全エネルギー（total energy）$E_{\mathrm{tot}}$ である．これは，原子核を指定の位置に固定したうえで（Born–Oppenheimer近似，第8章），その静電場のなかで電子が基底状態に落ち着いたときの，系全体のエネルギーである．原子核の運動エネルギーは含まれないので，これは絶対零度における内部エネルギーにあたる（厳密には，後述する零点振動エネルギーを除いたもの）．
全エネルギーそのものの絶対値には，あまり意味がない．意味があるのは差である．
式 \eqref{eq:1-energy-difference} が負ならば構造Aのほうが安定である．この単純な比較が，計算科学のもっとも基本的で，もっとも強力な道具になる．比較の対象は何でもよい．同じ組成で結晶構造だけ違うもの（結晶多型，polymorph）どうしを比べれば，どの多型が安定かがわかる．原子を1個抜いた構造と抜かない構造を比べれば，空孔の形成エネルギーがわかる．分子が表面に吸着した構造と，離れている構造を比べれば，吸着エネルギーがわかる．
なぜ？：全エネルギーの差だけで「安定性」を語ってよいのか
厳密には，有限温度で自発的に起こる変化を決めるのは自由エネルギー $F = U - TS$（定積の場合）である．全エネルギー $E_{\mathrm{tot}}$ は $T=0$ での $U$ にすぎず，エントロピー項 $-TS$ も，原子核の零点振動エネルギーも含まれていない．それでも $\Delta E_{\mathrm{tot}}$ の比較が意味をもつのは，結晶多型のエネルギー差が典型的には $10$〜$1000\ \mathrm{meV}$／化学式単位という大きさをもつのに対し，室温での $k_{\mathrm{B}}T$ が $25.9\ \mathrm{meV}$ と小さく，多くの場合エントロピー項が順位をひっくり返さないからである．逆に，エネルギー差が $k_{\mathrm{B}}T$ 程度しかないときには，フォノン計算からエントロピーを求めて自由エネルギーで比べ直す必要がある．「差が大きいときは $E_{\mathrm{tot}}$ で十分，差が小さいときは自由エネルギーまで踏み込む」という使い分けを覚えておこう．
1.5.2 新物質探索：希少金属を使わない半導体をつくれるか
全エネルギーの比較がもっとも劇的な形で使われるのが，まだ存在しない物質の探索である．計算機の中では，実験室ではつくれない組成・構造でも自由に組み立ててエネルギーを評価できる．この「まず計算で当たりをつけ，有望なものだけ実験する」という進め方をマテリアルズ・インフォマティクス（materials informatics）と呼ぶ．
典型例が，赤色発光ダイオードなどに使われる化合物半導体である．従来の高効率な赤色発光体には In や Ga といった希少金属が必須であった．そこで，豊富に存在する元素（Ca, Zn, N など）だけで同等の性能をもつ半導体をつくれないかを，計算で探索する研究が行われた．三元系 Ca–Zn–N の状態図上で安定に存在しうる組成と構造を計算で洗い出し，$\mathrm{CaZn_2N_2}$ という未報告の化合物が有望であることを予測，実際に高圧合成で実現してみせた研究がある〔Nature Communications 7, 11962 (2016)〕．
同じ発想を，$\mathrm{A_3BN}$ という組成の逆ペロブスカイト（inverse perovskite，あるいは反ペロブスカイト）に対して網羅的に適用した研究が，担当教員らによる〔Phys. Rev. Materials 4, 044601 (2020)〕である．通常のペロブスカイト $\mathrm{ABO_3}$ では陰イオン（酸素）が八面体の頂点を占めるが，逆ペロブスカイトでは陰イオンと陽イオンの役割が入れ替わり，N を中心に $\mathrm{A}$ 陽イオンが八面体を組む．ここで $\mathrm{A} = \mathrm{Mg, Ca, Sr, Ba}$，$\mathrm{B} = \mathrm{P, As, Sb, Bi}$ の全16通りの組成について，考えうる7種類の結晶型すべての全エネルギーを計算し，それぞれの組成で最も低いエネルギーをもつ構造を決定した．
| 予測された最安定構造 | 空間群 | 該当する組成 |
|---|---|---|
| 立方晶ペロブスカイト型 | $Pm\bar{3}m$ | $\mathrm{Mg_3SbN}$, $\mathrm{Ca_3SbN}$, $\mathrm{Sr_3SbN}$, $\mathrm{Mg_3BiN}$, $\mathrm{Ca_3BiN}$, $\mathrm{Sr_3BiN}$ |
| $\mathrm{CaTiO_3}$ 型（斜方晶） | $Pbnm$ | $\mathrm{Mg_3PN}$, $\mathrm{Ca_3PN}$, $\mathrm{Sr_3PN}$, $\mathrm{Ba_3PN}$, $\mathrm{Mg_3AsN}$, $\mathrm{Ca_3AsN}$, $\mathrm{Sr_3AsN}$ |
| $\mathrm{BaNiO_3}$ 型（六方晶） | $P6_3/mmc$ | $\mathrm{Ba_3AsN}$, $\mathrm{Ba_3SbN}$, $\mathrm{Ba_3BiN}$ |
この研究では，あわせて光吸収係数 $\alpha(h\nu)$ も計算されている．$\mathrm{Mg_3PN}$ と $\mathrm{Sr_3PN}$ は，実用材料である GaAs や CdTe に匹敵する立ち上がりの急峻な吸収端をもち，しかも希少金属も毒性元素も含まない．まだ誰も合成していない物質の光学特性を，先に言い当てているのである．これが型1の計算の到達点のひとつである．
例題1.3 網羅探索の計算コストを見積もる
上記の逆ペロブスカイト探索では，$\mathrm{A}$ 4種類 × $\mathrm{B}$ 4種類 = 16組成について，7種類の結晶型を検討した．
(1) 全エネルギー計算は何回必要か．
(2) 1回の計算に平均3時間かかり，128 並列で計算機を1台占有できるとすると，全体で何日かかるか．ただし各計算は同時に実行できないものとする．
(3) 実験でこの16組成すべてを合成し，構造を決定するとしたら，どれくらいの規模の仕事になると想像されるか．
解答 (1) $16\times 7 = 112$ 回である．ただし実際には，各構造について原子位置と格子定数を緩和させる構造最適化（第8章）を行うため，1「回」の中で何十回もエネルギーと力を評価している．
(2) $112\times 3 = 336$ 時間 $= 14.0$ 日．1台の計算機を2週間まわせば終わる規模である．
(3) 16組成の新物質を合成し，そのうち何割かは高圧合成や特殊な雰囲気制御が必要で，さらに単結晶X線回折で構造を決めるとなれば，研究室ひとつが数年がかりで取り組む仕事になる．しかも，4つの組成はそもそも誰も合成したことがない．計算が「実験に先立って当たりをつける」道具として決定的に有用であることが，この見積もりからよくわかる．
1.5.3 表面エネルギーと，分子の結合が切れるかどうか
全エネルギーの差で議論できるのは，バルクの結晶多型だけではない．結晶を切って表面をつくったときのエネルギー増加分が表面エネルギー（surface energy）$\gamma$ である．厚さ有限の薄板（スラブ）を計算セルに入れ，その全エネルギー $E_{\mathrm{slab}}$ と，同じ原子数のバルクのエネルギー $N E_{\mathrm{bulk}}$ を比べればよい．
式 \eqref{eq:1-surface-energy} の分母にある $2$ は，スラブに表と裏の二つの表面ができることに対応し，$A$ は表面の面積である．$\gamma$ が小さい面ほど「つくるのに損をしない」面であり，実際の結晶で現れやすい．単位としては $\mathrm{meV/\text{Å}^2}$ や $\mathrm{J/m^2}$ が用いられる．
この計算が思わぬ分野で役に立った例を紹介しよう．銅系の材料は抗ウイルス活性をもつことが知られているが，$\mathrm{Cu_2O}$ は大気中で $\mathrm{CuO}$ へと変化してしまい，活性が低下するという問題があった．ところが $\mathrm{CuO}$ に $\mathrm{La_2O_3}$ を組み合わせると，高い安定性と高い抗ウイルス活性の両方が維持されることが実験的に見出された．なぜ効くのかを第一原理計算で調べた研究が，担当教員らによる〔ACS Applied Materials & Interfaces 17, 61872 (2025)〕である．
手順はこうである．まず $\mathrm{La_2CuO_4}$ の (110)，(010)，(001)，(021) といったさまざまな表面について式 \eqref{eq:1-surface-energy} で表面エネルギーを計算し，酸素が乏しい条件・豊富な条件のそれぞれで最も安定な面を特定した．その結果，(110) 面が最安定であることがわかった．次に，その (110) 面の上に，タンパク質のジスルフィド結合 $\mathrm{S{-}S}$ の模型分子であるジメチルジスルフィド $(\mathrm{CH_3})_2\mathrm{S_2}$ を置き，$\mathrm{S{-}S}$ 結合が切れた状態と切れていない状態の全エネルギーを比較した．
式 \eqref{eq:1-dissociation} が負であれば，表面の上では $\mathrm{S{-}S}$ 結合が切れたほうが安定，すなわちこの表面はジスルフィド結合を切るということになる．ウイルスの表面タンパク質は，ジスルフィド結合によって立体構造を保っている．それが切れれば構造は壊れ，感染力を失う．こうして，抗ウイルス活性の起源がジスルフィド結合の切断にあることが，第一原理計算によって初めて実証された．
物理的意味：「なぜ効くのか」に答えるということ
実験は「$\mathrm{CuO}+\mathrm{La_2O_3}$ は効く」という事実を教えてくれる．しかし，なぜ効くのかは教えてくれない．全エネルギーの差 \eqref{eq:1-dissociation} という，たった1本の引き算が「表面上で $\mathrm{S{-}S}$ 結合が切れるから」という機構を与える．機構がわかれば，次にどんな物質を試すべきかが指針をもって決まる．$\mathrm{S{-}S}$ 結合を切る能力が高い表面をもつ酸化物を，計算で先に探せばよいからである．1.1.4項で述べた「計算という第三の柱」の役割が，ここに端的に現れている．
1.6 型2：結晶構造から電子バンドを得る
1.6.1 電子バンドとは何か
孤立した原子では，電子のとりうるエネルギーは $1s$, $2s$, $2p$, … といったとびとびの準位である（第2章）．ところが原子を $10^{22}$ 個も並べて結晶をつくると，隣り合う原子の軌道どうしが重なり合い，準位は無数に分裂して連続的な帯（バンド）になる．この帯を，電子の波数 $\bm{k}$ の関数として描いたものが電子バンド構造（electronic band structure）$\varepsilon_n(\bm{k})$ である．図1.1の下段に模式的に描いたものがそれにあたる．
バンド構造を見れば，その物質が金属か絶縁体かがただちにわかる．電子が詰まっている最も高いエネルギー（Fermi準位 $\varepsilon_{\mathrm{F}}$）の位置に電子の状態があれば金属，なければ絶縁体・半導体である．ダイヤモンドでは電子の詰まった価電子帯と空の伝導帯のあいだに約 $5.5\ \mathrm{eV}$ の禁制帯（バンドギャップ）が開いており，可視光（$1.6$〜$3.3\ \mathrm{eV}$）では電子を励起できない．だから無色透明で絶縁体なのである．グラファイトでは価電子帯と伝導帯がわずかに重なるためギャップがなく，電気が流れる．
1.6.2 化学結合は波動関数の重なり方で説明できる
バンド構造は「どのエネルギーに状態があるか」を教えてくれるが，材料設計の観点からはもう一歩踏み込みたい．その状態は，どの原子とどの原子の，どの軌道からできているのか．そしてその重なり方は結合を強めているのか，弱めているのか．第一原理計算は波動関数そのものを出力するので，これに答えることができる．
原子AとBの原子軌道 $\phi_A$，$\phi_B$ が近づいたとき，二つの重ね合わせ方がある．同位相で足し合わせた $\psi_{\mathrm{S}} \propto \phi_A + \phi_B$ と，逆位相で引いた $\psi_{\mathrm{A}} \propto \phi_A - \phi_B$ である．
- 結合性軌道（bonding state）$\psi_{\mathrm{S}}$：二つの原子核のあいだで波動関数が強め合う．電子密度 $\abs{\psi_{\mathrm{S}}}^2$ が原子核のあいだに厚く分布し，正電荷をもつ二つの原子核を電子が糊のようにつなぎとめる．エネルギーはもとの原子軌道より低い．
- 反結合性軌道（anti-bonding state）$\psi_{\mathrm{A}}$：二つの原子核のちょうど中間で波動関数が打ち消し合い，節（node）ができる．電子密度は核間で $0$ になり，糊の役目を果たさない．エネルギーはもとの原子軌道より高い．
なぜ？：いま数式なしで図だけ見せる理由
図1.5の内容は，本書のあとの章で完全に数式化される．第5章では二つの $1s$ 軌道がどれだけ重なっているかを表す重なり積分 $S$ を実際に積分して求め，第11章では永年方程式を解いて $\varepsilon_{\pm} = \dfrac{\alpha \pm \beta}{1 \pm S}$ という結合性・反結合性軌道のエネルギーを導出する．いまここで図だけを見せているのは，これから積む長い数式の階段が，最後にどんな景色に到達するのかを先に見ておいてほしいからである．ゴールの絵を頭に入れておけば，第3章の変分原理も第5章の積分計算も，何のためにやっているのかがぶれない．
例題1.4 2準位模型で結合の安定化を見積もる
簡単のため重なり積分を無視し（$S \simeq 0$），原子軌道のエネルギーを $\alpha$，二つの軌道を結びつける共鳴積分を $\beta\ (\lt 0)$ とすると，結合性・反結合性軌道のエネルギーは
$$ \begin{equation} \varepsilon_{\mathrm{bonding}} = \alpha + \beta, \qquad \varepsilon_{\mathrm{anti}} = \alpha - \beta \label{eq:1-two-level} \end{equation} $$となる（第11章で導出する）．$\beta = -2.5\ \mathrm{eV}$ とする．
(1) 電子が2個あり，両方が結合性軌道に入るとき，孤立原子の状態（各原子に1個ずつ）と比べてエネルギーはどれだけ下がるか．
(2) 電子が4個あり，結合性軌道と反結合性軌道の両方が満たされるとき，エネルギーはどうなるか．
解答 (1) 孤立原子2個に電子が1個ずつあるときの全エネルギーは $2\alpha$ である．結合してできた分子で2個の電子が結合性軌道を占めると，全エネルギーは $2(\alpha+\beta) = 2\alpha + 2\beta$．よって変化は
$$ \Delta E = (2\alpha + 2\beta) - 2\alpha = 2\beta = 2\times(-2.5\ \mathrm{eV}) = -5.0\ \mathrm{eV}. $$$5.0\ \mathrm{eV}$ だけ安定化する．これが共有結合が生じる理由である．
(2) 電子4個では結合性軌道に2個，反結合性軌道に2個入る．全エネルギーは
$$ 2(\alpha+\beta) + 2(\alpha-\beta) = 4\alpha, $$孤立原子4電子分 $4\alpha$ とまったく同じで，$\Delta E = 0$ である．つまり結合による安定化が反結合による不安定化で完全に相殺され，結合が生じない．He 原子どうしが $\mathrm{He_2}$ 分子をつくらないのは，この理屈による（実際には重なり積分 $S\gt 0$ の効果で反結合性軌道の押し上げのほうが大きくなり，むしろ反発する）．結合性軌道だけが占有されているか，反結合性軌道まで占有されているかが，結合の有無を決める——これが次項の COHP の出発点である．
1.6.3 COHP・COOP・COBI — 結合性の度合いを数値化する
例題1.4の考え方を，原子が $10^{22}$ 個並んだ結晶に対して，しかもエネルギーの関数として実行できるようにしたのが COHP である．
定義：COHP・COOP・COBI
- COHP（crystal orbital Hamilton populations，結晶軌道ハミルトン占有数）：端的に言えば共鳴積分の密度である．エネルギー $\varepsilon$ における，ある原子対のあいだの結合が電子系のエネルギーをどれだけ変化させているかを表す．無次元量である．符号の約束は $$\mathrm{COHP}\lt 0 \;\Longrightarrow\; \text{結合性},\qquad \mathrm{COHP}\gt 0 \;\Longrightarrow\; \text{反結合性}.$$
- COOP（crystal orbital overlap populations，結晶軌道重なり占有数）：隣り合う原子間の重なり積分の密度．単位は $\mathrm{eV^{-1}}$．符号の約束は COHP と逆で $$\mathrm{COOP}\gt 0 \;\Longrightarrow\; \text{結合性},\qquad \mathrm{COOP}\lt 0 \;\Longrightarrow\; \text{反結合性}.$$
- COBI（crystal orbital bond index，結晶軌道結合指数）：COHP を電子数（占有数）まで考慮したもので，化学でいう結合次数に対応する量．エネルギー積分値 ICOBI が $1$ なら単結合，$2$ なら二重結合に近い意味をもつ．
注意：COHP と COOP で符号の向きが逆であること
COHP が負なら結合性，COOP が正なら結合性である．まぎらわしいが，定義を思い出せば理由は明快である．COHP はハミルトニアン行列要素 $H_{ij}$（＝共鳴積分 $\beta$，これは負）を含むので，結合性の状態では負になる．COOP は重なり積分 $S_{ij}$（同位相で重なれば正）を含むので，結合性の状態では正になる．符号を丸暗記せず「エネルギーを下げるほう（負）が結合性なのが COHP，重なりが正のほうが結合性なのが COOP」と，由来から思い出せるようにしておこう．なお論文の図では，直感に合わせるため COHP の符号を反転した $-\mathrm{COHP}$ を横軸にとることが多い．$-\mathrm{COHP}$ のグラフでは右側（正）が結合性である．軸のラベルを必ず確認してほしい．
この解析の威力を示す例が，担当教員らによる $\mathrm{N_2}$ と $\mathrm{O_2}$ の結合解析である〔J. Phys. Chem. C 125, 7959 (2021)〕．高校化学では $\mathrm{N_2}$ は三重結合，$\mathrm{O_2}$ は二重結合と習う．分子軌道図を描けば，$\mathrm{N_2}$ では $\sigma$ と $\pi$ の結合性軌道がすべて満たされ反結合性軌道は空，$\mathrm{O_2}$ では $\pi^{*}$ 軌道に電子が2個入る，という違いが見える．COHP・COOP・COBI を計算すると，この定性的な描像が数値として再現される．すなわち $\sigma$ 結合と $\pi$ 結合それぞれの寄与を分離して，どちらがどれだけ結合を支えているかまで定量的に言えるのである．
1.6.4 エピタキシャル歪みによる金属絶縁体転移
電子バンドが読めると，「構造をこう変えれば電気を通さなくなる」という予言ができる．担当教員らの研究〔Phys. Rev. Materials 2, 125001 (2018)〕はその好例である．
$\mathrm{La_3Ni_2O_7}$ という層状ニッケル酸化物を，格子定数のわずかに違う基板の上に薄膜として成長させると，膜は基板の格子に合わせて面内方向に引き伸ばされたり縮められたりする．これをエピタキシャル歪み（epitaxial strain）という．歪み $\varepsilon_{\mathrm{s}} = -2\%$ を加えた場合の電子バンドを計算したところ，歪みのないバルクではバンドが Fermi 準位を横切る金属であったものが，歪みを加えるとギャップが開いて絶縁体になることが予測された．結晶を数パーセント歪ませるだけで，電気を通す物質が通さない物質に変わるのである．
なぜそんなことが起こるのか．鍵はFermi面ネスティング（Fermi surface nesting）とPeierls不安定性（Peierls instability）にある．Fermi面とは，波数空間で Fermi 準位に対応する等エネルギー面のことである．この Fermi 面の一部が，ある一定のベクトル $\bm{q}$ だけ平行移動すると別の部分にぴったり重なる——つまり平らな面どうしが向かい合っている——とき，Fermi面がネスティングしているという．$\mathrm{La_3Ni_2O_7}$ の場合，Fermi面が第一Brillouinゾーンの半分でうまく表現できる形になっており，強いネスティングが生じていた．
物理的意味：Peierls不安定性はなぜ金属を絶縁体に変えるのか
最も簡単な1次元の例で考えよう．原子間隔 $a$ で等間隔に並んだ1次元鎖に，電子が各原子あたり1個ある（半分だけ詰まった＝ハーフフィリング）とする．強束縛近似（第14章）でバンドは
$$ \begin{equation} \varepsilon(k) = -2t\cos(ka) \label{eq:1-tb-band} \end{equation} $$となる（$t\gt 0$ は隣接原子への飛び移りの積分）．ハーフフィリングでは Fermi 波数が $k_{\mathrm{F}} = \pi/2a$ となり，バンドは $k=\pm\pi/2a$ で Fermi 準位を横切る．これは金属である．
ところがここで，原子が2個ずつ組をつくるように交互に近づいたり離れたりする変位（二量体化，dimerization）が起きたとしよう．周期が $2a$ に倍化するので，Brillouinゾーンの境界は $\pi/2a$ に移り，ちょうど $k_{\mathrm{F}}$ の位置に来る．変位の振幅に比例する大きさ $\Delta$ の周期ポテンシャルが加わると，そこでバンドが分裂して
$$ \begin{equation} \varepsilon_{\pm}(k) = \pm\sqrt{4t^2\cos^2(ka) + \Delta^2} \label{eq:1-peierls-gap} \end{equation} $$となり，$2\Delta$ の大きさのギャップが開く．ここが肝心である．ギャップが開くとき，Fermi 準位のすぐ下にあった占有状態はエネルギーが下がり，すぐ上にあった空の状態はエネルギーが上がる．上がった状態には電子が入っていないので，損はしない．得だけが残る．したがって電子系は必ずエネルギーを下げる．この「格子を歪ませてでもギャップを開けたい」という駆動力が Peierls不安定性である．格子を歪ませる弾性エネルギーの損より電子系の得が勝てば，金属は自発的に絶縁体へ転移する．Fermi面ネスティングが強いほど，この機構で得られるエネルギーが大きい．
1.6.5 遺伝的アルゴリズムによる表面構造の探索
ここまでは「与えられた結晶構造」を計算する話であったが，そもそも構造がわからない場合はどうするか．固体の表面がまさにそれである．バルクをどこかで切っただけの理想的な表面は，実際にはほとんど存在しない．表面の原子は組成を変え，位置をずらし，しばしばバルクとは異なる周期構造（表面再構成）をつくる．候補構造は組み合わせ爆発的に多く，人間の勘で全部試すことはできない．
そこで，遺伝的アルゴリズム（genetic algorithm）と第一原理計算を組み合わせる．表面の原子配置を「遺伝子」とみなし，エネルギーの低い個体を選択して交叉・突然変異させる操作を繰り返すことで，安定な表面構造を自動的に探し出す方法である．担当教員らはこの手法をペロブスカイト $\mathrm{ABO_3}$（$\mathrm{SrTiO_3}$，$\mathrm{KTaO_3}$，$\mathrm{BaTiO_3}$ など8種類）の表面に適用し，Cation exchange，Checker，Stripe，Zigzag と名づけられる4種類の再構成パターンを系統的に見出した〔Chem. Mater. 35, 2047 (2023)〕．
これが何の役に立つのか．同じ物質でも，表面構造が違えば真空準位から測った価電子帯上端・伝導帯下端の位置が数 eV も変わることが計算で示された．すなわち，固体から電子を取り出すのに必要なエネルギーが表面構造によってまったく違うということであり，光触媒や電極材料としての性能を左右する．表面で起こる化学反応の選択性も変わる．バルクの性質だけを見ていても材料は設計できないという教訓である．
なぜ？：機械にやらせるべき仕事とは
候補構造を1つずつ人間が思いついて計算する，という進め方には限界がある．人間の勘は，これまで見てきた構造の範囲を超えないからである．遺伝的アルゴリズムや，1.5節で見た網羅探索は，人間が思いつかない構造にたどり着くための道具である．計算科学を学ぶ意義のひとつは，「これは人間がやるべき仕事か，機械にやらせるべき仕事か」を見分けられるようになることにある．物理を考えるのは人間の仕事，可能性をしらみつぶしに試すのは機械の仕事である．
1.7 型3：結晶構造からフォノンバンドを得る
1.7.1 フォノンバンドとは何か
結晶中の原子は，絶対零度であってもその平衡位置に静止しているわけではなく，たがいにばねでつながれたおもりのように振動している．この格子振動を量子化したものがフォノン（phonon）である．電子バンドが電子のエネルギー $\varepsilon$ を波数 $\bm{k}$ の関数として描いたものであったのと同じように，フォノンバンド（フォノン分散関係）は格子振動の振動数 $\omega$ を波数 $\bm{q}$ の関数として描いたものである．
ダイヤモンドとグラファイトのフォノンバンドを比べると，どちらも最高振動数が $40$〜$47\ \mathrm{THz}$ に達する．炭素は軽く，しかも共有結合が強いためである．しかしグラファイトには，層間の弱い van der Waals 結合に対応する低振動数の枝があり，層に垂直な振動は非常に柔らかい．フォノンバンドを見れば，その物質のどの方向がどれだけ「硬い」かが一目でわかる．
1.7.2 振動エネルギーとは，要するに振動数のことである
「振動エネルギー準位」といういかめしい名前がついているが，量子力学における調和振動子のエネルギー準位は $E_n = \left(n + \tfrac{1}{2}\right)\hbar\omega$（$n = 0, 1, 2, \dots$）であり，決めるべきものは振動数 $\omega$ だけである．そしてこの $\omega$ は，古典力学の単振動から求まる．まずはそこを復習しよう．
ある原子の質量を $m$，平衡位置からのずれを $x$ とする．この原子に，ずれと逆向きの復元力 $-kx$（$k \gt 0$ をばね定数，あるいは力定数という）が働くとすると，Newton の運動方程式は
$$ \begin{equation} m\ddot{x} = -kx \label{eq:1-newton-restoring} \end{equation} $$である．ここで $\ddot{x} = \dd^2 x/\dd t^2$ である．一方，単振動の運動方程式は
$$ \begin{equation} \ddot{x} = -\omega^2 x \label{eq:1-sho} \end{equation} $$と書けたことを思い出そう．式 \eqref{eq:1-newton-restoring} の両辺を $m$ で割ると $\ddot{x} = -(k/m)x$ となるから，式 \eqref{eq:1-sho} と見比べて $\omega^2 = k/m$，すなわち
を得る．力定数 $k$ がわかれば振動数がわかり，逆に振動数がわかれば力定数がわかる．第一原理計算では，原子を平衡位置からわずかにずらして力を計算し，力定数 $k$（結晶では力定数行列）を求めてから式 \eqref{eq:1-omega-real} を使う．これがフォノン計算の骨格である．
1.7.3 例題：単振動の運動方程式を完全に解く
式 \eqref{eq:1-omega-real} は「解の形が $\sin$ や $\cos$ であること」を前提にした議論だった．ここでは前提を置かずに，微分方程式そのものを初期条件つきで解いてみよう．次の節で登場する「復元力の符号が反転した場合」を扱うために，この作業が決定的に重要になる．
例題1.5 単振動と，その符号を反転した場合の完全解
質量 $m$ の質点が，平衡位置からのずれ $x$ に対して力を受けて運動する．初期条件を
$$ x(0) = 0, \qquad \left.\frac{\dd x}{\dd t}\right|_{t=0} = v_0 \;(\gt 0) $$とする．すなわち，時刻 $t=0$ に平衡位置を初速度 $v_0$ で通過する．線形微分方程式
$$ a_2\frac{\dd^2 x}{\dd t^2} + a_1\frac{\dd x}{\dd t} + a_0 x = 0 $$の解が $x(t) = A e^{\lambda_1 t} + B e^{\lambda_2 t}$（ただし $\lambda_1, \lambda_2$ は特性方程式 $a_2\lambda^2 + a_1\lambda + a_0 = 0$ の解）と書けることを用いて，次の二つを解け．
(1) $m\ddot{x} = -kx$（復元力が働く通常の場合，$k\gt 0$）
(2) $m\ddot{x} = +kx$（力の向きが反転した場合，$k\gt 0$）
解答 (1) 方程式を移項して $m\ddot{x} + kx = 0$ とすると，$a_2 = m$, $a_1 = 0$, $a_0 = k$ である．特性方程式は
$$ \begin{equation} m\lambda^2 + k = 0 \quad\Longrightarrow\quad \lambda^2 = -\frac{k}{m} \quad\Longrightarrow\quad \lambda = \pm i\sqrt{\frac{k}{m}} = \pm i\omega, \qquad \omega \equiv \sqrt{\frac{k}{m}} \;(\gt 0) \label{eq:1-char-stable} \end{equation} $$となる．$k\gt 0$，$m\gt 0$ なので $-k/m \lt 0$ であり，平方根の中身が負になるため解は純虚数である．よって一般解は
$$ x(t) = A e^{i\omega t} + B e^{-i\omega t} . $$初期条件を代入する．まず $x(0)=0$ より
$$ A + B = 0 \quad\Longrightarrow\quad B = -A . $$次に，両辺を $t$ で微分して
$$ \frac{\dd x}{\dd t} = i\omega A e^{i\omega t} - i\omega B e^{-i\omega t} $$であるから，$t=0$ とおいて
$$ \left.\frac{\dd x}{\dd t}\right|_{t=0} = i\omega A - i\omega B = i\omega (A - B) = v_0 . $$ここに $B = -A$ を代入すると $i\omega\,(A - (-A)) = 2i\omega A = v_0$，すなわち
$$ A = \frac{v_0}{2i\omega}, \qquad B = -\frac{v_0}{2i\omega} . $$したがって
$$ x(t) = \frac{v_0}{2i\omega}e^{i\omega t} - \frac{v_0}{2i\omega}e^{-i\omega t} = \frac{v_0}{2i\omega}\left(e^{i\omega t} - e^{-i\omega t}\right) . $$ここで Euler の公式から導かれる関係 $\sin\theta = \dfrac{e^{i\theta} - e^{-i\theta}}{2i}$，すなわち $e^{i\theta} - e^{-i\theta} = 2i\sin\theta$ を用いると
$$ x(t) = \frac{v_0}{2i\omega}\cdot 2i\sin(\omega t) $$となり，$2i$ が約分されて
を得る．振幅 $v_0/\omega$，角振動数 $\omega$ の単振動である．$x(0)=0$，$\dot{x}(0) = v_0\cos 0 = v_0$ となっており，確かに初期条件を満たしている．
解答 (2) こんどは $m\ddot{x} - kx = 0$ であるから，$a_2 = m$, $a_1 = 0$, $a_0 = -k$ である．特性方程式は
$$ \begin{equation} m\lambda^2 - k = 0 \quad\Longrightarrow\quad \lambda^2 = +\frac{k}{m} \quad\Longrightarrow\quad \lambda = \pm\sqrt{\frac{k}{m}} = \pm\Omega, \qquad \Omega \equiv \sqrt{\frac{k}{m}} \;(\gt 0) \label{eq:1-char-unstable} \end{equation} $$となる．今度は平方根の中身が正なので，解は純虚数ではなく実数である．ここが (1) との決定的な違いである．一般解は
$$ x(t) = A e^{\Omega t} + B e^{-\Omega t} . $$$x(0)=0$ より $A + B = 0$，すなわち $B = -A$．微分すると
$$ \frac{\dd x}{\dd t} = \Omega A e^{\Omega t} - \Omega B e^{-\Omega t} $$だから，$t=0$ で
$$ \Omega(A - B) = \Omega\,(A + A) = 2\Omega A = v_0 \quad\Longrightarrow\quad A = \frac{v_0}{2\Omega},\quad B = -\frac{v_0}{2\Omega}. $$したがって
$$ x(t) = \frac{v_0}{2\Omega}\left(e^{\Omega t} - e^{-\Omega t}\right) . $$双曲線正弦関数の定義 $\sinh\theta = \dfrac{e^{\theta} - e^{-\theta}}{2}$ を使えば $e^{\Omega t} - e^{-\Omega t} = 2\sinh(\Omega t)$ であるから
を得る．$\sinh$ は $t$ が大きいところで $\sinh(\Omega t) \simeq \tfrac{1}{2}e^{\Omega t}$ とふるまうから，$x(t)$ は時間とともに指数関数的に発散する．振動して戻ってくることは二度とない．
数学ノート：$\sin$ と $\sinh$ は同じ式の顔違いである
解 \eqref{eq:1-sol-stable} と \eqref{eq:1-sol-unstable} が別物に見えるかもしれないが，じつは同じ式である．式 \eqref{eq:1-sol-stable} で形式的に $\omega \to i\Omega$（つまり $k \to -k$）と置き換えてみよう．$\sin(i\theta) = i\sinh\theta$ という関係を使うと
$$ \frac{v_0}{\omega}\sin(\omega t) \;\xrightarrow{\;\omega \to i\Omega\;}\; \frac{v_0}{i\Omega}\sin(i\Omega t) = \frac{v_0}{i\Omega}\cdot i\sinh(\Omega t) = \frac{v_0}{\Omega}\sinh(\Omega t) $$となり，式 \eqref{eq:1-sol-unstable} にぴったり一致する．「振動数が虚数になる」ということと「解が指数関数的に発散する」ということは，同じ事実の言い換えにすぎない．次項の議論の核心はここにある．
1.7.4 復元力の符号が反転したら — 虚数振動数の登場
ここで，次のような問いを立ててみよう．もしも，平衡位置から $x$ だけずらしたときに，働く力のベクトルが逆方向を向いてしまったらどうなるのか．つまり
$$ \begin{equation} m\ddot{x} = +kx \label{eq:1-newton-inverted} \end{equation} $$という状況である．式 \eqref{eq:1-sho} と同じ形に書けば $\ddot{x} = +\omega^2 x$ ということになる．この式を単振動の公式 \eqref{eq:1-omega-real} に無理やり当てはめると，力定数が $-k$ に変わったのだから
となってしまう．振動数が虚数になったのである．これはもはや振動ではない．例題1.5(2)で確かめたとおり，この場合の解は $\sinh$ であり，時間とともに発散する．すなわち，平衡位置だと思っていた点からずれたほうが安定であることを示唆している．この状態を，虚数振動数を有する振動状態（imaginary phonon，あるいはソフトモード）という．
なぜ？：虚数の振動数などという奇妙なものを，なぜ大事にするのか
虚数振動数は「計算が失敗した印」ではなく，もっとも有益な出力のひとつである．理由はこうだ．第一原理計算では，まず何らかの結晶構造を仮定して原子に働く力がゼロになるまで構造最適化する．力がゼロになれば，その構造はエネルギー曲面上の停留点である．しかし，停留点は極小点とは限らない．峠（鞍点）や極大点かもしれない．力がゼロというだけでは区別がつかないのである．
そこで，その点のまわりで原子を微小変位させ，力定数（＝エネルギーの2階微分）の符号を調べる．すべての方向で $k\gt 0$ なら本物の極小点であり，その構造は動的に安定である．どこか1つの方向で $k\lt 0$ なら，その方向にずらすとエネルギーが下がる——つまりより安定な構造が，その方向の先に存在する．虚数振動数は，それを探すべき方向を名指しで教えてくれる矢印なのである．1.7.6項で見るように，この矢印をたどることで新しい強誘電体が発見された．
1.7.5 非調和ポテンシャルと二重井戸
いま述べたことを，力ではなくポテンシャルエネルギーの形で見ておこう．原子変位 $x$ に対する原子間ポテンシャルエネルギー $U(x)$ を，$x=0$ のまわりで Taylor 展開する．$x=0$ が停留点なので1次の項は消え，対称性から奇数次の項も落ちるとすれば
と書ける．$x^2$ の項を調和項，$x^4$ 以上の項を非調和項（anharmonic term）と呼ぶ．原子間力は $F(x) = -\dd U/\dd x$ だから
$$ \begin{equation} F(x) = -\frac{\dd U}{\dd x} = -kx - 4Cx^3 \label{eq:1-force-stable} \end{equation} $$である．$k\gt 0$ かつ $x$ が小さい領域では $F \simeq -kx$ となって復元力が働き，式 \eqref{eq:1-omega-real} により実の振動数 $\omega = \sqrt{k/m}$ が得られる．これが動的に安定な場合である．
一方，調和項の係数の符号が反転した場合，すなわち
のときは
$$ \begin{equation} F(x) = -\frac{\dd U}{\dd x} = +kx - 4Cx^3 \label{eq:1-force-unstable} \end{equation} $$となり，$x$ が小さいところでは $F \simeq +kx$，つまり変位を拡大する向きに力が働く．これが式 \eqref{eq:1-newton-inverted} の状況であり，振動数は式 \eqref{eq:1-omega-imag} により虚数になる．これが虚数振動数を有するフォノンである．
ただし $x$ が大きくなると $x^3$ の項が効いてきて力の向きが再び戻り，変位は無限には成長しない．$U(x)$ が極小になる位置は $\dd U/\dd x = 0$ から
$$ -kx + 4Cx^3 = 0 \quad\Longrightarrow\quad x\left(4Cx^2 - k\right) = 0 $$より $x=0$（これは極大）と
$$ \begin{equation} x_{\min} = \pm\sqrt{\frac{k}{4C}} \label{eq:1-double-well-min} \end{equation} $$である．$x=0$ を挟んで左右対称に二つの極小をもつこの形を二重井戸ポテンシャル（double-well potential）という．極小におけるエネルギーは，式 \eqref{eq:1-double-well-min} を \eqref{eq:1-potential-unstable} に代入して
$$ \begin{align} U(x_{\min}) &= -\frac{1}{2}k\cdot\frac{k}{4C} + C\left(\frac{k}{4C}\right)^2 \notag\\ &= -\frac{k^2}{8C} + C\cdot\frac{k^2}{16C^2} \notag\\ &= -\frac{k^2}{8C} + \frac{k^2}{16C} = -\frac{k^2}{16C} \label{eq:1-double-well-depth} \end{align} $$となる．式 \eqref{eq:1-double-well-depth} は負であり，もとの $x=0$ の構造より，変位した構造のほうがエネルギーが低いことを定量的に示している．
例題1.6 力定数から振動数を計算する
質量 $m = 16\ \mathrm{u}$（酸素原子，$1\ \mathrm{u} = 1.6605\times10^{-27}\ \mathrm{kg}$）の原子が，力定数 $k = 10\ \mathrm{eV/\text{Å}^2}$ のポテンシャルの中で振動している．
(1) $k$ を SI 単位（$\mathrm{N/m}$）に直せ．
(2) 角振動数 $\omega$ と振動数 $\nu = \omega/2\pi$ を求めよ．
(3) 力定数の符号が反転して $k = -10\ \mathrm{eV/\text{Å}^2}$ になったとき，振動数はどうなるか．
解答 (1) $1\ \mathrm{eV} = 1.602\times10^{-19}\ \mathrm{J}$，$1\ \text{Å} = 10^{-10}\ \mathrm{m}$ であるから
$$ 1\ \mathrm{eV/\text{Å}^2} = \frac{1.602\times10^{-19}\ \mathrm{J}}{(10^{-10}\ \mathrm{m})^2} = \frac{1.602\times10^{-19}}{10^{-20}}\ \mathrm{J/m^2} = 16.02\ \mathrm{N/m}. $$よって $k = 10\times 16.02 = 160.2\ \mathrm{N/m}$ である．
(2) 質量は $m = 16\times 1.6605\times10^{-27} = 2.657\times10^{-26}\ \mathrm{kg}$．式 \eqref{eq:1-omega-real} より
$$ \omega = \sqrt{\frac{k}{m}} = \sqrt{\frac{160.2}{2.657\times10^{-26}}} = \sqrt{6.030\times10^{27}} = 7.766\times10^{13}\ \mathrm{rad/s}, $$ $$ \nu = \frac{\omega}{2\pi} = \frac{7.766\times10^{13}}{6.2832} = 1.236\times10^{13}\ \mathrm{Hz} = 12.4\ \mathrm{THz}. $$参考までに分光学でよく使う波数に直すと $\tilde{\nu} = \nu/c = 1.236\times10^{13}/(2.998\times10^{10}\ \mathrm{cm/s}) = 412\ \mathrm{cm^{-1}}$ となり，酸化物の光学フォノンとして妥当な値である．
(3) 式 \eqref{eq:1-omega-imag} より
$$ \omega = \sqrt{\frac{-160.2}{2.657\times10^{-26}}} = i\times 7.766\times10^{13}\ \mathrm{rad/s}, \qquad \nu = 12.4\,i\ \mathrm{THz}. $$絶対値は同じで，虚数になる．フォノンバンドの図では，慣習としてこれを $-12.4\ \mathrm{THz}$ の位置，すなわち横軸より下にプロットする．論文の図で「振動数が負の枝」を見たら，それは負の振動数ではなく虚数振動数の意味である．
1.7.6 フォノンバンドにおける虚数振動数と，新規強誘電体の発見
結晶では，原子1個ではなく単位胞内のすべての原子が波数 $\bm{q}$ をもつ波として集団的に振動する．したがって力定数 $k$ は行列（力定数行列）になり，振動数はその固有値の平方根として，波数ごとに複数本の枝（ブランチ）として求まる．ある波数 $\bm{q}$ である枝の固有値が負になれば，その $\bm{q}$ とその変位パターンに沿って構造を変形させたほうがエネルギーが下がる，ということになる．
この考え方から新しい材料が実際に生まれた例を紹介する．担当教員らによる層状ペロブスカイト $\mathrm{Li_2SrNb_2O_7}$ の研究〔Chem. Mater. 33, 1257 (2021)〕である．
まず，対称性の高い $Amam$ 構造を仮定してフォノン分散を計算したところ，Brillouinゾーン境界の Y 点において虚数振動数をもつ枝が見つかった（$\mathrm{Y_2^-}$ モードと名づけられる変位パターン）．次に，その変位パターンに沿って原子を実際にずらしながら全エネルギーを計算すると，図1.7の青い曲線とまったく同じ二重井戸型のグラフが得られた．極小に落ち込んだ先の構造は，対称性の低い $Pna2_1$ 構造であり，これは自発分極をもつ，すなわち強誘電体である．
この機構には名前がついている．$\mathrm{A_3B_2O_7}$ 型の層状ペロブスカイトでは，それ自体は分極を生まない二つのゾーン境界フォノンモード（酸素八面体の傾き $\mathrm{X_3^-}$ と回転 $\mathrm{X_2^+}$）が同時に凍結することで，二次的な効果として分極が発生する．これをハイブリッド不正則強誘電性（hybrid improper ferroelectricity）という〔N. A. Benedek and C. J. Fennie, Phys. Rev. Lett. 106, 107204 (2011)〕．通常の強誘電体のように分極そのものが不安定になるのではなく，分極とは無関係に見える二つの回転モードの積として分極が現れるので「不正則（improper）」と呼ばれる．
物理的意味：勘と経験に頼らない構造探索
従来，新しい強誘電体を探すには「この元素を置換したらどうか」という化学的な勘と，膨大な試行錯誤が必要だった．ところがフォノン計算を使えば，手順は次のように機械化される．
- 対称性の高い構造を仮定してフォノン分散を計算する．
- 虚数振動数をもつモードを探す．なければ，その構造は安定であり，話はここで終わる．
- 虚数振動数のモードが指し示す変位方向に沿って構造を歪ませ，エネルギーが下がることを確認する．
- 落ち込んだ先の構造の対称性を調べる．中心対称性が破れていれば，その物質は強誘電体の候補である．
虚数振動数という「一見すると計算のエラーに見える量」が，実は安定な結晶構造を探索するための最も信頼できる指針になっているのである．第1章の冒頭で掲げた「原子の並び方から物性を演繹する」という思想が，ここでは「物性が欲しい並び方を，計算に探させる」という逆向きの使い方にまで発展している．
1.8 本書の地図
ここまで見てきたのは，第一原理計算という完成した道具で何ができるかという話であった．本書の残りの13章は，その道具をどうやって作るかを，一段ずつ積み上げていく．到達点は，分子の量子化学計算を自分で理解して実行でき，固体の電子状態を定性的に説明できるようになることである．
全体の流れは次のとおりである．まず第2章で水素原子のSchrödinger方程式に戻る．水素原子は，多電子系に進む前の唯一の「厳密に解ける」よりどころである．ところが電子が2個以上になると方程式は解けなくなるので，第II部（第3〜6章）で近似解を求める二つの道具，変分原理と摂動論を学び，実際に $1s$ 軌道の重なり積分を計算して基底関数 STO-3G を導出し，その近似解を厳密解と比べる．第III部（第7〜10章）では電子を増やしていく．Pauliの排他律とSlater行列式から交換エネルギーが現れ，Born–Oppenheimer近似によって原子核と電子の運動が分離され，全エネルギーの内訳が明らかになる．ここまでくれば実際に分子の量子化学計算が実行できる．第IV部（第11〜12章）では永年方程式を解いて共有結合とイオン結合を数式から理解し，群論によって軌道の対称性から結合を定性的に判断する術を得る．最後の第V部（第13〜14章）で分子から固体へ渡り，Madelungエネルギーとイオン結晶の電子状態にたどりつく．そこで見るバンド構造は，本章の図1.1で眺めたものと同じものである．
なぜ？：なぜ第2章で「水素原子」に戻るのか
本章で紹介した研究は，どれも原子が数十から数百個も含まれる複雑な系である．それなのに次の章で電子1個の水素原子に戻るのは，遠回りに見えるかもしれない．しかし理由は明快である．水素原子はSchrödinger方程式が厳密に解ける唯一の実在原子であり，そこで得られる $1s$, $2s$, $2p$ といった軌道の形と，エネルギー $E_n = \epsilon_{1s}/n^2$（$\epsilon_{1s} = -13.606\ \mathrm{eV}$）という結果が，これ以降のすべての近似の出発点であり検算の基準になるからである．第5章では水素の $1s$ 軌道どうしの重なり積分を実際に計算し，第6章ではその近似解を厳密解と比べて誤差を測る．基準になる厳密解を知らなければ，近似が良いのか悪いのかを判断することすらできない．
付録1.A COHP の定義
1.6.3項で紹介した COHP について，定義式を挙げておく．本章の他の部分を読むために必要な内容ではないので，興味のある読者だけ目を通せばよい．式の意味が完全にわかるのは，第11章で永年方程式を，第14章でバンド構造を学んだ後になる．
数学ノート：COHP の定義式
原子サイト $R$ に置かれた，方位量子数と磁気量子数の組 $L$ で指定される基底関数を $\chi_{RL}$ と書く．ハミルトニアン行列の要素は
$$ \begin{equation} H_{RL,R'L'} = \bra{\chi_{RL}}\hat{H}\ket{\chi_{R'L'}} \label{eq:1-cohp-hamiltonian} \end{equation} $$である．これは第11章で導入する共鳴積分を，結晶の基底関数について書いたものにほかならない．一方，電子系のバンドエネルギー（一電子準位の和）は，$j$ 番目の一電子準位を $\epsilon_j$，その占有数を $f_j$（$0$〜$2$）として
$$ \begin{equation} E^{\mathrm{band}} = \int^{\epsilon_{\mathrm{F}}} \dd\epsilon \sum_j f_j\,\epsilon_j\,\delta(\epsilon_j - \epsilon) \label{eq:1-cohp-band} \end{equation} $$と書ける．$\delta$ は Dirac のデルタ関数で，「エネルギー $\epsilon$ のところに準位 $\epsilon_j$ があるときだけ数え上げる」という役目をもつ．ここで，被積分関数を基底関数の対 $(RL, R'L')$ に分解すると
$$ \begin{equation} \sum_j f_j\,\epsilon_j\,\delta(\epsilon_j - \epsilon) = \sum_{RL}\sum_{R'L'} H_{RL,R'L'}\cdot N_{RL,R'L'}(\epsilon) \equiv \sum_{RL}\sum_{R'L'} \mathrm{COHP}_{RL,R'L'}(\epsilon) \label{eq:1-cohp-def} \end{equation} $$となる．$N_{RL,R'L'}(\epsilon)$ は状態密度行列であり，エネルギー $\epsilon$ における状態が基底 $\chi_{RL}$ と $\chi_{R'L'}$ にどれだけ由来するかを表す．式 \eqref{eq:1-cohp-def} の各項が COHP である．
意味を読み取ろう．$H_{RL,R'L'}$ はハミルトニアン行列要素（＝共鳴積分）で，$N_{RL,R'L'}(\epsilon)$ は状態密度であるから，COHP は共鳴積分の「密度」である．式 \eqref{eq:1-cohp-band} と \eqref{eq:1-cohp-def} を合わせれば，COHP をエネルギーについて Fermi 準位まで積分した量（ICOHP）は，その原子対がバンドエネルギーに与えている寄与そのものになる．結合性の状態では $H_{RL,R'L'}\lt 0$ であるから COHP は負となり，エネルギーを下げている．これが「$\mathrm{COHP}\lt 0$ なら結合性」という符号の約束の由来である．
なお，式 \eqref{eq:1-cohp-band} の左辺 $E^{\mathrm{band}}$ はエネルギーの次元をもち，右辺の積分測度 $\dd\epsilon$ もエネルギーの次元をもつから，COHP 自身は無次元量である．一方 COOP は重なり積分（無次元）と状態密度（$\mathrm{eV^{-1}}$）の積なので $\mathrm{eV^{-1}}$ の次元をもつ．定義：R. Dronskowski and P. E. Blöchl, J. Phys. Chem. 97, 8617 (1993).
1.9 まとめと演習
1.9.1 この章のまとめ
- 思想：唯物論と要素還元主義に立てば，原子の並び方から物性を演繹的に説明できるはずである．計算科学の基本構図は「Input＝結晶構造／分子構造，Output＝物性」であり，入力に実験値を使わない計算を第一原理計算と呼ぶ．同じ炭素からできたダイヤモンド（絶縁体，$E_g\simeq 5.5\ \mathrm{eV}$）とグラファイト（半金属）の対比が，この構図の威力を端的に示す．
- スケール：材料の現象は，まずスケールを考える．原子核 $10^{-15}\ \mathrm{m}$，原子 $10^{-10}\ \mathrm{m}$，分子 $10^{-9}\ \mathrm{m}$，高分子・組織 $10^{-6}\ \mathrm{m}$，細胞 $10^{-3}\ \mathrm{m}$，巨視的物体 $10^{0}\ \mathrm{m}$ という階層があり，それぞれ素粒子物理学・物性物理学／物理化学・高分子化学・分子生物学が対応する．
- 原子スケールの手法：第一原理計算（DFT，0 K の静的物性），分子動力学計算（古典MD／第一原理MD／機械学習ポテンシャル，有限温度の動的物性），モンテカルロ法（配置の熱平衡）．MD は時間平均（動画）を取り，MC は確率平均（漫画）を取る．
- メゾ・マクロの手法：フェーズフィールド法（組織変化・分極ドメイン），有限要素法（応力分布・構造振動）．下位スケールの結果を上位のパラメータとして渡すのがマルチスケール計算である．
- 型1（全エネルギー）：$\Delta E = E_{\mathrm{tot}}(\mathrm{A})-E_{\mathrm{tot}}(\mathrm{B})$ の比較だけで，結晶多型の安定性，新物質探索（逆ペロブスカイト $\mathrm{A_3BN}$），表面エネルギー \eqref{eq:1-surface-energy}，分子の解離 \eqref{eq:1-dissociation}（抗ウイルス活性の機構）が議論できる．
- 型2（電子バンド）：金属か絶縁体かがわかる．波動関数の重なり方から結合性軌道と反結合性軌道が区別でき，COHP（$\lt 0$ なら結合性）・COOP（$\gt 0$ なら結合性）・COBI（結合次数）で定量化される．Fermi面ネスティングとPeierls不安定性は，格子を歪ませてギャップを開くことで電子系のエネルギーを下げる機構である．
- 型3（フォノンバンド）：振動数は $\omega=\sqrt{k/m}$ で決まる．$m\ddot{x}=-kx$ の解は $x(t)=(v_0/\omega)\sin(\omega t)$ で振動するが，復元力の符号が反転した $m\ddot{x}=+kx$ の解は $x(t)=(v_0/\Omega)\sinh(\Omega t)$ となって発散する．このとき $\omega = \sqrt{-k/m}=i\sqrt{k/m}$ が虚数振動数であり，平衡点だと思っていた構造が不安定であることを意味する．ポテンシャルは二重井戸 $U(x)=-\tfrac12 kx^2+Cx^4$ となり，極小は $x=\pm\sqrt{k/4C}$，深さは $-k^2/16C$ である．$\mathrm{Li_2SrNb_2O_7}$ の $\mathrm{Y_2^-}$ ソフトモードは，この指針に従って新規強誘電体（ハイブリッド不正則強誘電性）を見出した実例である．
1.9.2 演習問題
演習1.1 初期条件を変えた単振動
例題1.5と同じ二つの運動方程式について，初期条件を
$$ x(0) = x_0\;(\gt 0), \qquad \left.\frac{\dd x}{\dd t}\right|_{t=0} = 0 $$（平衡位置から $x_0$ だけずらして，静かに手を放す）に変えて解け．すなわち
(1) $m\ddot{x}=-kx$ の解 $x(t)$ を求めよ．
(2) $m\ddot{x}=+kx$ の解 $x(t)$ を求めよ．
(3) (2) の解について，$t\to\infty$ での挙動を述べ，$x_0$ を2倍にしたとき「変位がある値 $X$ に達するまでの時間」がどう変わるかを論じよ．
ヒント：いずれも一般解 $x(t)=Ae^{\lambda_1 t}+Be^{\lambda_2 t}$ から出発する．(1) では $\lambda=\pm i\omega$，(2) では $\lambda=\pm\Omega$ である（$\omega=\Omega=\sqrt{k/m}$）．$x(0)=x_0$ から $A+B=x_0$，$\dot{x}(0)=0$ から $A=B$ が出るので $A=B=x_0/2$．あとは $\cos\theta = \dfrac{e^{i\theta}+e^{-i\theta}}{2}$，$\cosh\theta = \dfrac{e^{\theta}+e^{-\theta}}{2}$ を使えばよい．(3) では $\cosh(\Omega t)\simeq \tfrac12 e^{\Omega t}$ と近似し，$X = x_0 e^{\Omega t}/2$ を $t$ について解いて $x_0$ 依存性を見る．$\ln$ が現れるので，初期変位を2倍にしても到達時間はごくわずかしか短くならない．「不安定性は初期条件の大きさにあまりよらず，$\Omega$ の大きさで決まる」ことを確認せよ．
演習1.2 二重井戸ポテンシャルの数値計算
式 \eqref{eq:1-potential-unstable} の二重井戸ポテンシャルについて，$k = 111\ \mathrm{eV/\text{Å}^2}$，$C = 309\ \mathrm{eV/\text{Å}^4}$ とする．
(1) 極小の位置 $x_{\min}$ と，そのときのエネルギー $U(x_{\min})$ を数値で求め，図1.7と一致することを確かめよ．
(2) $x=0$（対称性の高い構造）を基準にしたとき，変位した構造がどれだけ安定化するかを $\mathrm{meV}$ 単位で答えよ．この値を室温の熱エネルギー $k_{\mathrm{B}}T = 25.9\ \mathrm{meV}$ と比べ，この構造変化が室温で保たれるかどうか論じよ．
(3) 質量 $m = 16\ \mathrm{u}$ の原子について，$x=0$ における虚数振動数の大きさ $\abs{\omega}$ を求め，$\mathrm{THz}$ 単位で答えよ．
ヒント：(1) は式 \eqref{eq:1-double-well-min} と \eqref{eq:1-double-well-depth} にそのまま代入する．(2) では $1\ \mathrm{eV}=1000\ \mathrm{meV}$．得られる値は $k_{\mathrm{B}}T$ の約100倍であることに注目せよ．(3) は例題1.6と同じ手順で，$k$ を $\mathrm{N/m}$ に直してから $\abs{\omega}=\sqrt{\abs{k}/m}$ を計算する．$\nu = \omega/2\pi$ に直すのを忘れないこと．
演習1.3 手法の選択
次の(a)〜(e)のそれぞれについて，表1.1の中から最も適した計算手法を1つ選び，理由を2〜3行で述べよ．
(a) まだ合成されていない組成 $\mathrm{Ba_3SbN}$ が，どの結晶構造で最も安定になるかを知りたい．
(b) $\mathrm{Cu_{1-x}Zn_x}$ 合金で，Cu と Zn が規則配列するか無秩序に混ざるかを，組成と温度の関数として知りたい．
(c) 高分子の融体が細い管を流れるときの粘性を，$100\ \mathrm{ns}$ にわたって調べたい．
(d) 強誘電体薄膜に電圧をかけたとき，分極ドメイン壁が $\mu\mathrm{m}$ の距離を動く様子を知りたい．
(e) 自動車のサスペンション部品に繰り返し荷重をかけたときの，応力集中箇所を知りたい．
ヒント：判断の軸は二つある．ひとつは長さのスケール（何 m の現象か），もうひとつは時間のスケール（何 s の現象か，そもそも時間が要るのか）である．(b) では「原子がサイトを入れ替わる」事象が問題であり，これは MD の到達時間では起こらないことを思い出そう（1.3.3項）．(a) では時間の情報がまったく要らないことに注意．
演習1.4 COHP・COOP の符号を読む
ある化合物の第一原理計算から，隣り合う金属原子 M と酸素原子 O の対について，次のような結果が得られたとする．
- $-8\ \mathrm{eV}$ 付近のエネルギー領域で $\mathrm{COHP}\lt 0$
- $-1\ \mathrm{eV}$ 付近のエネルギー領域（Fermi 準位のすぐ下）で $\mathrm{COHP}\gt 0$
(1) それぞれの領域の状態は結合性か反結合性か．
(2) 同じ計算で COOP を求めた場合，それぞれの領域で COOP の符号はどうなるか．
(3) この化合物から電子を少し抜く（正孔ドープする）と，M–O 結合は強くなるか弱くなるか．理由とともに述べよ．
(4) 逆に電子を少し加える（電子ドープする）と，M–O 結合はどうなるか．
ヒント：(1)(2) は定義（1.6.3項）をそのまま適用する．COHP と COOP で符号の向きが逆であることに注意．(3)(4) は例題1.4を思い出そう．Fermi 準位直下の状態は反結合性である．正孔ドープとは，その反結合性の状態から電子を抜くことに相当する．反結合性軌道の電子が減ればどうなるかを考えればよい．この「Fermi 準位付近に反結合性状態がある物質は，正孔ドープで結合が強くなる」という論法は，実際の材料設計で頻繁に用いられる．
演習1.5 網羅探索の設計
ペロブスカイト $\mathrm{ABO_3}$ について，$\mathrm{A}$ に4種類，$\mathrm{B}$ に6種類の元素を割り当て，それぞれ5種類の結晶型を検討する網羅探索を計画する．
(1) 必要な全エネルギー計算は何回か．
(2) 1回の構造最適化に平均5時間かかるとして，20件を同時並列で流せる計算機環境があるとき，全体の所要時間は何日か．
(3) さらに，得られた最安定構造すべてについてフォノン計算（1件あたり構造最適化の20倍のコスト）を行うと，追加でどれだけの時間がかかるか．
(4) (3) の結果を踏まえ，フォノン計算をすべての構造について行うべきか，それとも絞るべきかを論じよ．絞るとすればどんな基準が考えられるか．
ヒント：(1) は $4\times6\times5$．(2) は総時間を並列数で割る．(3) では「最安定構造」は組成ごとに1つなので $4\times6=24$ 件であることに注意．(4) では，フォノン計算の目的が「その構造が動的に安定かどうかの確認」と「虚数振動数からより安定な構造を探すこと」の二つであることを思い出そう．全エネルギーが明らかに高い構造まで調べる意味はあるだろうか．
演習1.6 同素体を Input／Output の言葉で説明する
ダイヤモンドとグラファイトは，どちらも炭素のみからなる同素体である．
(1) 「同素体である」という説明と，「原子の並び方が違うから物性が違う」という説明は，どこがどう違うか．要素還元主義の立場から述べよ．
(2) 第一原理計算でダイヤモンドとグラファイトの全エネルギーを比較したとき，常温常圧ではグラファイトのほうがわずかに安定であることが知られている．それにもかかわらずダイヤモンドが室温で存在し続けられるのはなぜか．
(3) 炭素の第3の形態としてグラフェン（グラファイトの1層）がある．これは半金属だが，グラファイトとは異なる特異なバンド構造をもつ．本章で学んだ枠組みで，グラフェンの性質を予測するには何を Input とし，何を Output として計算すればよいか述べよ．
ヒント：(1) 「同素体だから」は分類の名前を与えただけで，因果を説明していない．1.1.1項の議論に戻ろう．(2) は熱力学的安定性と速度論的（動力学的）安定性の区別の問題である．エネルギーの高い状態から低い状態へ移るには，途中の峠（活性化障壁）を越えなければならない．1.7.5項の二重井戸の図で，井戸の深さと山の高さが別の量であることを思い出すとよい．(3) Input は「グラフェン1層分の周期構造（層間の相互作用がない）」であり，Output としては電子バンドを計算することになる．$\mathrm{K}$ 点でのバンドの交わり方に注目せよ．
1.9.3 参考文献
- 原田 義也『量子化学（上巻）』裳華房（2007）. — 本書第2章以降の内容を，より詳しく数式で追いたいときの標準的な日本語教科書．
- C. Kittel（宇野良清ほか訳）『固体物理学入門（上・下）』第8版, 丸善（2005）. — 電子バンド，フォノン，Peierls不安定性など，本章1.6・1.7節の背景となる固体物理の基礎．
- P. A. Cox（魚崎浩平ほか訳）『固体の電子構造と化学』技報堂出版（1989）. — 化学の言葉で固体のバンド構造を語る良書．結合性・反結合性の議論を固体に拡張する視点が得られる．
- R. Dronskowski and P. E. Blöchl, "Crystal orbital Hamilton populations (COHP): energy-resolved visualization of chemical bonding in solids based on density-functional calculations", J. Phys. Chem. 97, 8617 (1993). — 付録1.A の原論文．
- Y. Hinuma, F. Oba et al., "Discovery of earth-abundant nitride semiconductors by computational screening and high-pressure synthesis", Nature Communications 7, 11962 (2016). — 1.5.2項，$\mathrm{CaZn_2N_2}$ の予測と合成．
- Y. Mochizuki, F. Oba et al., "Theoretical exploration of mixed-anion antiperovskite semiconductors $\mathrm{M_3XN}$", Phys. Rev. Materials 4, 044601 (2020). — 1.5.2項，逆ペロブスカイト $\mathrm{A_3BN}$ の網羅探索と光吸収係数．
- Y. Mochizuki, F. Oba et al., "Strain-engineered metal-to-insulator transition and orbital polarization in nickelate superlattices", Phys. Rev. Materials 2, 125001 (2018). — 1.6.4項，エピタキシャル歪みによる金属絶縁体転移とFermi面ネスティング．
- Y. Mochizuki et al., "Cation-size effect on ferroelectricity in layered perovskites", Chem. Mater. 33, 1257 (2021). — 1.7.6項，$\mathrm{Li_2SrNb_2O_7}$ の $\mathrm{Y_2^-}$ ソフトモードと新規強誘電体の予測．
- Y. Mochizuki et al., "Chemical bonding analysis based on crystal orbital indices", J. Phys. Chem. C 125, 7959 (2021). — 1.6.3項，COHP・COOP・COBI による結合解析．
- Y. Mochizuki et al., "Surface reconstruction and band alignment of perovskite oxides explored by genetic algorithm", Chem. Mater. 35, 2047 (2023). — 1.6.5項，遺伝的アルゴリズムによる $\mathrm{SrTiO_3}$ 等の表面構造探索．
- R. Kiribayashi, Y. Mochizuki, A. Nakajima et al., ACS Applied Materials & Interfaces 17, 61872 (2025). — 1.5.3項，$\mathrm{La_2CuO_4}$ 表面によるジスルフィド結合切断と抗ウイルス活性．
- N. A. Benedek and C. J. Fennie, "Hybrid improper ferroelectricity: a mechanism for controllable polarization-magnetization coupling", Phys. Rev. Lett. 106, 107204 (2011). — 1.7.6項，ハイブリッド不正則強誘電性の提唱．