4 全エネルギーは比較の量

計算が返す「全エネルギー」は、たとえば −93.45 Ry のような大きな負の数だ。
この数字ひとつを見て、私たちは何か言えるのだろうか。

この章のゴール：全エネルギーの意味は単独の絶対値ではなく「差」に宿ることを理解し、二つの計算を正しく比べるために何をそろえるべきかを述べられるようになる。

1 大きな負の数、それ自体は何を語るか

前章で、シリコンの結晶を単位胞・格子ベクトル・原子位置という数件のデータとして計算機に置けるようになった。ibrav=2 を選び、celldm(1) \(\approx 10.26\) bohr と書いた、あの数行である。では、その配置を第一原理という変換器に通したとき、最初に返ってくる数字は何か。全エネルギーである。

画面には、たとえば

\[ E_\text{tot} \approx -93.453\ \mathrm{Ry} \]

のような、大きな負の数がひとつ現れる（具体的な値は擬ポテンシャルの選び方で変わる。ここでは代表的な桁として見ておけばよい）。\(1\ \mathrm{Ry} = 13.606\ \mathrm{eV}\) をかけて直せば、およそ \(-1270\ \mathrm{eV}\)。ずいぶん深い谷の底にいるように見える。負号は、結晶が「計算のたまたま選んだ原点」より下に落ち着いていることを表しているが、その原点がどこなのかは、この数字だけでは何も言っていない。

この一個の数字を見つめて、私たちはシリコンについて何か言えるだろうか。正直に答えれば、ほとんど何も言えない。この数字は「どこを高さの基準（ゼロ）に選んだか」に丸ごと依存しているからである。

高校で習った位置エネルギーを思い出してほしい。\(mgh\) の \(h\) をどこから測るかは、こちらの勝手だった。机の上をゼロにしてもよいし、床をゼロにしてもよい。同じ物体の位置エネルギーが、基準しだいで \(+5\ \mathrm{J}\) にも \(-3\ \mathrm{J}\) にもなる。それでも誰も困らなかった。物理に効くのは高さの差――落ちるときに解放される仕事――のほうだったからである。

第一原理計算が返す全エネルギーも、これと同じ立場にある。ある約束ごとの原点から測った「深さ」なのだ。そして、その原点を決めているのが、前章でファイル名だけ顔を出した擬ポテンシャルである。擬ポテンシャルは、原子核と内側の電子をひとまとめにした「芯」を、扱いやすい有効ポテンシャルで置き換える道具だが、その芯の組み立て方には流儀がいくつもある。流儀が変われば、エネルギーを測りはじめる原点も一緒に動く。だから、まったく同じ物理のシリコンでも、ある擬ポテンシャルでは \(-93.453\ \mathrm{Ry}\)、別の擬ポテンシャルでは \(-300\ \mathrm{Ry}\) というふうに、絶対値はいくらでもずれてよい。富士山の標高が、海面から測れば \(3776\ \mathrm{m}\)、地球の中心から測れば約 \(6378\ \mathrm{km}\) になるのと同じことである。山は動いていない。ものさしの根もとが違うだけだ。

では、この数字は役立たずなのか。そうではない。原点さえ固定してしまえば、これは立派な座標になる。同じ原点から測った二つの深さを引き算したとき、全エネルギーははじめて物理を語りだす。

2 意味は「差」に宿る

差をとる、とは具体的にどういうことか。やり方は大きく二つある。

ひとつは、格子定数を少しだけ変えること。前章の celldm(1) をほんの少し縮めたり伸ばしたりして、そのたびに全エネルギーを計算する。たとえば \(a = 5.40,\ 5.43,\ 5.47\ \mathrm{\text{Å}}\) と振ってみる。返ってくる絶対値はどれも \(-93.453\ \mathrm{Ry}\) 前後の大きな負の数で、桁を眺めても区別はつかない。ところが差をとると、\(a\) を変えるとエネルギーがどちらへ動くかという谷の形が現れる。その谷の底にあたる \(a\) こそ、計算機が「いちばん安定」と考える格子定数である。この谷を実際に描いて格子定数を読み取る作業は、第8章で正面から行う。

もうひとつは、原子を動かすこと。結晶の中の一個の原子を平衡位置から少しずらし、ずらす前との全エネルギーの差をとる。この差が、原子を元へ引き戻そうとする力の出どころになる。配置を入れると力が返る、と第1章で言ったが、その力もまた、突きつめれば全エネルギーの差から来ている。

差をとる視点に慣れると、読み取れるものが一段増える。格子定数を振って描いた谷は、底のまわりで放物線に近い形をしている。その放物線の開き具合――底の鋭さ――は、結晶をぐっと縮めるのにどれだけのエネルギーが要るか、つまり結晶の硬さを映している。底の位置（格子定数）だけでなく、底の曲がり方（硬さ）まで、すべては全エネルギーの差から読み取れるのだ。一個の絶対値が黙っていたぶん、差は何重にも雄弁である。

なぜ差をとると物理が浮かび上がるのか。仕組みは一行で書ける。同じ約束ごとで計算した二つの全エネルギーを、それぞれ

\[ E_1 = (\text{物理の本体})_1 + C, \qquad E_2 = (\text{物理の本体})_2 + C \]

と書いてみよう。\(C\) は、原点のとり方と擬ポテンシャルの流儀で決まる「共通のゲタ」である。同じ設定で計算しているかぎり、\(C\) は両方にそっくり同じだけ乗っている。だから差をとれば

\[ E_2 - E_1 = (\text{物理の本体})_2 - (\text{物理の本体})_1 \]

となって、\(C\) はきれいに消える。さっきまで絶対値を支配していた、あの大きくて気まぐれなゲタが、引き算ひとつで姿を消すのである。

消えるのはゲタだけではない。後の章で出てくる計算の粗さ――波をどこまで細かく表すか、周期をどれだけ細かく標本化するか――から来る誤差も、二つの計算で条件をそろえておけば、やはり同じだけ両方に乗る。すると差の中で大半が打ち消し合う。これが効いてくるのが第5章・第6章で、絶対値はまだ落ち着いていないのに、差のほうはとっくに安定している、という場面に何度も出会う。絶対値より差のほうが、はるかに信用できる。これは第一原理計算を読むときの、いちばん大事な勘所のひとつである。

flowchart LR
    E1["計算 1<br>E₁ = 物理₁ + ゲタ C + 誤差 δ"]
    E2["計算 2<br>E₂ = 物理₂ + ゲタ C + 誤差 δ"]
    SUB["差をとる<br>E₂ - E₁"]
    OUT["物理₂ - 物理₁<br>（C と δ は消えた）"]
    E1 --> SUB
    E2 --> SUB
    SUB --> OUT
    style OUT fill:#d4edda,stroke:#28a745

図 1: 同じ土俵で計算した二つの全エネルギーから差をとると、約束ごとで決まる共通のゲタ C と、刻みの粗さからくる共通の誤差 δ がそろって消え、物理の差だけが残る。

3 同じ土俵で比べる ― 何をそろえるべきか

ゲタ \(C\) と共通の誤差が差の中で消えるのは、「同じ約束ごとで計算したかぎり」という条件があってのことだった。約束ごとが食い違った二つの計算を引き算すれば、\(C\) は消え残り、差は物理を語ってくれない。だから、何をそろえれば比べてよいのかをはっきりさせておく必要がある。

そろえるべきものは、つきつめれば次の四つである。

交換相関の近似（汎関数）。本書が標準で使う近似――後の章で名前を挙げる――を、両方の計算で同じにする。
擬ポテンシャル。原点を決める当人だから、同じ元素ファイルを使う。これが違えばゲタ \(C\) が食い違う。
波の解像度（ecutwfc、必要なら ecutrho）。電子の波をどこまで細かく表すかの刻みである。
周期の標本化（K_POINTS で指定する k 点メッシュ、たとえば \(6\times6\times6\)）。無限結晶をどれだけ細かく見るかの刻みである。

この四つが選ばれているのには理由がある。前節の図で差の中から消えてくれたゲタ \(C\) と共通の誤差 \(\delta\) は、まさにこの四つが決めているからだ。汎関数と擬ポテンシャルが原点（ゲタ \(C\)）を握り、ecutwfc と k 点が刻みの粗さ（誤差 \(\delta\)）を握っている。四つを両方の計算でそろえておけば、\(C\) も \(\delta\) も両辺に同じだけ乗り、引き算で消えてくれる。消えてほしい量をつくっている当人たちだからこそ、ここを動かさずにおくのである。

そして、これらをそろえたうえで、変えてよいのはただ一つ、物理として知りたい変数だけにする。格子定数を知りたいなら格子定数だけを動かし、力を知りたいなら原子位置だけを動かす。残りはすべて固定する。これはまさに、実験で言う「一つだけ条件を変えて比べる」対照実験の作法そのものである。

ここで一つ、もっともらしく見える落とし穴を見ておこう。二つの格子定数 \(a_1,\ a_2\) の安定性を比べたいのに、うっかり \(a_1\) を ecutwfc \(=30\ \mathrm{Ry}\) で、\(a_2\) を \(40\ \mathrm{Ry}\) で計算してしまったとする。すると差 \(E(a_2)-E(a_1)\) の中には、格子定数の違いという知りたい物理と、刻みを変えたことによる余計な効果とが混ざり込み、より分けられなくなる。刻みをそろえてさえいれば差の中で消えてくれたはずの共通の誤差が、刻みを変えたせいで消え残るのである。だから比較の前に、知りたい変数以外はすべて固定する。これを数字で確かめるのが、章末の小課題 3.2 である。なお、刻みを「どこまで細かくすれば十分か」という問いそのものは、第5章で手を動かして見届ける。

そろえるべき四つは、どれも次章で開く一枚の入力ファイルの中に、それぞれの行として書き込まれている。比較の土俵は、入力ファイルを書いた時点ですでに宣言されているのだ。

4 凝集エネルギー：ばらばらの原子と結晶を比べる

差が物理を語るという話の、いちばん分かりやすい実例が凝集エネルギーである。問いはこうだ。ばらばらに離れていたシリコン原子が集まって結晶になると、エネルギーはどれだけ下がるのか。下がった分が、結晶を結びつけている糊の強さにあたる。

これを差として定義する。結晶の中の原子一個あたりの全エネルギーを \(E_\text{結晶}/N\)、遠くに孤立した原子一個の全エネルギーを \(E_\text{原子}\) とすると、凝集エネルギーは

\[ E_\text{coh} = E_\text{原子} - \frac{E_\text{結晶}}{N} \]

で与えられる。結晶のほうが安定（深い谷）だから \(E_\text{結晶}/N\) のほうが低く、差 \(E_\text{coh}\) は正になる。この正の値が大きいほど、結晶は強く結ばれている。

ここで「同じ土俵」の作法が効いてくる。\(E_\text{原子}\) を出すには、大きな空の箱にシリコン原子を一個だけ置いて計算する。結晶と孤立原子はそもそも別の系だから、単位胞の大きさも k 点のとり方も違ってよい――孤立原子の箱は十分に大きくとり、k 点は一点で足りる。そろえるべきは、一原子あたりのエネルギーの原点を決めている擬ポテンシャルと汎関数、そして ecutwfc の三つである。ここさえ同じなら、二つの全エネルギーは同じ原点から測った深さになり、引き算がゲタを消してくれる。

ある擬ポテンシャルと汎関数で、シリコンのダイヤモンド型 2 原子セルの全エネルギーが \(E_\text{結晶} \approx -93.453\ \mathrm{Ry}\) と出たとする。シリコンの凝集エネルギーは、実験では約 \(4.6\ \mathrm{eV/atom}\) である。

結晶の原子一個あたりの全エネルギー \(E_\text{結晶}/N\) を求めよ。
実験値の凝集エネルギー \(4.6\ \mathrm{eV/atom}\) を、\(1\ \mathrm{Ry}=13.606\ \mathrm{eV}\) を使って \(\mathrm{Ry/atom}\) に直し、(1) で求めた一原子あたりの全エネルギーの大きさと比べて、凝集エネルギーが全エネルギーの絶対値の何割にあたるかを見積もれ。

ダイヤモンド型 2 原子セルには原子が 2 個入るから \[ \frac{E_\text{結晶}}{N} = \frac{-93.453}{2} \approx -46.73\ \mathrm{Ry/atom}. \]
実験値を Ry に直すと \[ 4.6\ \mathrm{eV/atom} \div 13.606 \approx 0.338\ \mathrm{Ry/atom}. \] これを (1) の一原子あたりの全エネルギーの大きさと比べると \[ \frac{0.338}{46.73} \approx 0.007, \] すなわち凝集エネルギーは全エネルギーの絶対値のわずか約 \(0.7\ \%\) にあたる。

ここで注目してほしいのは桁である。一原子あたりの全エネルギーは \(-46.73\ \mathrm{Ry}\)、eV に直せば \(-636\ \mathrm{eV}\) ほどの深い谷だった。孤立した原子の全エネルギーもやはり同じくらい深い谷にあるはずで、凝集エネルギー \(0.338\ \mathrm{Ry}\) はその二つの大きな数のわずかな差、全体の約 \(0.7\ \%\) にすぎない。化学のすべて――結晶が結晶でいられる理由――が、その小さな差の中に畳み込まれている。だからこそ、絶対値を数 \(\%\) の精度でしか出せなくても、孤立原子まで同じ土俵で計算して差をとれば、実験の \(4.6\ \mathrm{eV}\) にちゃんと届く。これが「意味は差に宿る」の、いちばん具体的な姿である。

なお、この凝集エネルギーを計算だけで出すには、結晶とは別に、孤立原子の全エネルギー \(E_\text{原子}\)（さきに触れた、大きな空の箱の計算）をもう一度求める必要がある。それは別の一手なので、本章では \(4.6\ \mathrm{eV/atom}\) を実験値として受け取っておく。

シリコン結晶を二通り計算した。格子定数 \(a_1\) では ecutwfc \(=30\ \mathrm{Ry}\) を使って \(E\approx -93.45\ \mathrm{Ry}\)、別の格子定数 \(a_2\) では ecutwfc \(=40\ \mathrm{Ry}\) を使って \(E\approx -93.75\ \mathrm{Ry}\) を得た。

「\(a_2\) のほうがエネルギーが低いから、\(a_2\) のほうが安定だ」と結論してよいか。理由とともに答えよ。
\(a_1\) と \(a_2\) のどちらが安定かを正しく比べるには、計算をどう設定し直せばよいか。

結論できない。二つの計算は、格子定数も ecutwfc も同時に変えてしまっている。差 \(-0.3\ \mathrm{Ry}\) の中には、格子定数の違いによる物理の差と、刻みを細かくしたことによる見かけのエネルギー低下とが混ざっていて、より分けられない。変えてよいのは一度に一つ、という対照実験の作法を破っているのである。
ecutwfc を両方そろえて（たとえばどちらも \(40\ \mathrm{Ry}\) にして）、\(a_1\) と \(a_2\) を計算し直す。刻みがそろえば共通の誤差は差の中で打ち消し合い、残った差は格子定数の違いだけを映す。これが「同じ土俵で比べる」ことの実際である。

「全エネルギーが \(-93.453\ \mathrm{Ry}\) も \(-300\ \mathrm{Ry}\) も大きな負の数なのだから、桁が大きいほど強く結合していて安定なのだろう」と読みたくなる。だが、その桁の大きさを決めているのは結合の強さではなく、擬ポテンシャルが選んだ原点の位置である。同じシリコンが、擬ポテンシャルを替えただけで \(-93.453\) にも \(-300\) にもなる。だから、設定の違う二つの絶対値を並べて安定性を比べることはできない。

ではこの数字に意味はないのかというと、そうではない。原点を固定したまま――同じ擬ポテンシャル・同じ汎関数・同じ刻みのまま――格子定数や原子位置を動かして差をとれば、その差は結合の強さも、格子の安定な間隔も、はっきり語ってくれる。桁そのものではなく、同じものさしで測った差を読む。それがこの章で身につける目である。

5 この章のまとめ

全エネルギーという最初の出力を、正しく読む構えができた。

全エネルギーの絶対値（たとえば \(-93.453\ \mathrm{Ry}\)）は、原点のとり方と擬ポテンシャルの流儀で決まる「ゲタ」を含んでいて、単独ではほとんど物理を語らない。位置エネルギーの基準が勝手なのと同じである。
意味は差に宿る。同じ設定で計算した二つの全エネルギーを引き算すると、共通のゲタも、刻みの粗さからくる共通の誤差も大半が消え、物理の差だけが残る。だから絶対値より差のほうがはるかに信用できる。
比べてよいのは同じ土俵の二つだけである。汎関数・擬ポテンシャル・ecutwfc・k 点をそろえ、知りたい物理変数だけを動かす。
凝集エネルギーは差の代表例で、\(E_\text{coh}=E_\text{原子}-E_\text{結晶}/N\) という、大きな二つの数のわずかな差にあたる。実験ではおよそ \(4.6\ \mathrm{eV/atom}\)、一原子あたりの全エネルギー（約 \(-46.73\ \mathrm{Ry}\)）のわずか \(0.7\ \%\) にすぎない。

そろえるべき土俵――汎関数・擬ポテンシャル・刻み・k 点――は、どれも次章で開く一枚の入力ファイルの各行として書き込まれている。その一枚は、計算機への命令の羅列であると同時に、「この物理を、この近似で、この細かさで解く」という物理の宣言文でもある。次章では、その各行がどんな仮定を置いているのかを、一行ずつ読み解いていこう。