9  出した数字を、実在と突き合わせる

計算機に渡す前に、私たちは実験に近い格子定数を入力へ置いた。手で見積もった値も、そこに重なる。
だが――その心強い一致は、計算が選び取ったのか、それとも私たちが初めから置いたのか。

この章のゴール：入力に置いた格子定数と、計算が返した数字とを見分けられるようになる。第1章の手計算・実験値・入力に置いた値の三つを一本の線で結び、「計算自身に格子を選ばせる」には何をすればよいか（次の一手）まで言えるようになる。

前章で、SCF が電子と場のつじつま合わせだと分かった。反復が閉じ、画面には収束した全エネルギーがひとつ載っている。準備は整った。いよいよこの数字を、本物のシリコンに突き合わせる番である。最初の標的は、第1章で鉛筆だけで当てにいった、あの格子定数だ。ただし突き合わせの前に、ひとつ取り違えやすい段差がある。画面の数字は、私たちが入力に置いたものなのか、計算が返してきたものなのか。この区別を、格子定数はいちばん最初に問うてくる。

1 第1章の見積もりを、覚えているか

第1章で、私たちは計算機を開く前にひとつ実験をした。シリコンの密度と原子量だけを手がかりに、原子一個の体積を出し、立方単位胞に 8 個という事実を重ねて、格子定数を \(a \approx 5.4\ \mathrm{\text{Å}}\) と見積もった。あのとき「後で計算機と突き合わせる」と約束した。その「後で」が、いま来た。

ところが、いざ突き合わせようとすると、ひとつ手前で引っかかる。第1章で確かめたとおり、第一原理計算という変換器が直接返すのは、ある配置での力と全エネルギーである。「シリコンの格子定数は○○ Å です」という答えを、変換器はそのままの形では手渡してくれない。こちらが格子定数 \(a\) をひとつ決めて入力すれば、その \(a\) に対する全エネルギー \(E(a)\) と、結晶を縮めようとする（あるいは広げようとする）力が返ってくる。それだけである。

ならば、問いの立て方を変えればよい。「格子定数はいくらか」を、変換器が答えられる形へ翻訳するのだ。結晶がいちばん落ち着く間隔とは、それ以上縮めても広げてもエネルギーが上がってしまう間隔のことである。式で言えば、

\[ \frac{dE}{da} = 0 \quad(\text{かつ下に凸}) \]

を満たす \(a\) ――すなわち \(E(a)\) 曲線の底である。そこでは、結晶を縮めたい力と広げたい力がちょうど打ち消し合い、各原子にはたらく正味の力がゼロになる。平衡の格子定数とは、エネルギーが最小になり、力が消える間隔だ。これなら変換器に尋ねられる。

だから手順はこうなる。全エネルギーを格子定数の関数として何点か計算し、最小を与える \(a\) を読む。この \(E(a)\) 曲線には名前がついていて、状態方程式（あるいは E–V 曲線）と呼ぶ。ありがたいことに、底のごく近くではこの曲線はほぼ放物線になる。だから数点あれば、頂点として平衡 \(a\) を読み取れる。

なぜ放物線でよいのか。底では定義により \(dE/da = 0\) だから、平衡 \(a_0\) のまわりで \(E(a)\) をテイラー展開すると、一次の項が消えて、生き残る先頭は二次の項になる。

\[ E(a) \approx E_0 + \tfrac{1}{2}\,k\,(a-a_0)^2 . \]

底の近くなら、より高次の項は無視できるほど小さい。つまり数点を放物線で結ぶのは、近似の手抜きではなく、最小の近くでは曲線が本当に放物線に見える、という物理の帰結である。しかもこの \(k\)――放物線の開き具合――は捨て物ではない。\(E(a)\) 曲線がどれだけ急に立ち上がるかは、結晶を押し縮めるのにどれだけ力がいるか、すなわち体積弾性率（かたさ）そのものだ。同じ一本の曲線が、格子定数だけでなく、もうひとつ実験と突き合わせられる量まで抱えている。本章では格子定数に集中するが、頭の隅に置いておきたい。

ただし、ここまではすべて「もし格子を振ったら」という方法の話だ。ならば私たちの実行は、この振る作業を実際にやったのだろうか。次節で、入力に置いたものと計算が返したものを、正直に棚卸ししてみよう。

2 私たちが置いた格子と、出てきた数字

答えを先に言おう。私たちの実行は、いま言った手順を踏んでいない。第4章で読んだ入力ファイルを思い出そう。celldm(1) = 10.26――Å に戻せば 5.43 Å、実験値 5.431 Å にほとんど重なる格子を、私たちは初めから入力に置いた。そして calculation = 'scf' は、その固定した格子のまま電子状態と全エネルギーを一度だけ解く指示だった。今回走らせたのは、格子を振って谷を探す計算ではなく、実験に近い一点での SCF である。

だから、いま胸を張って言えるのは、格子定数そのものではない。言えるのは、その一点で計算がよく収束した、ということだ。第5章と第6章で確かめたとおり、ecutwfc を上げても k 点を細かくしても、全エネルギーが動くのは meV/atom の幅にとどまり、値は −93.45 Ry 付近に落ち着いた。電子状態の側も走り切り、最高被占準位（highest occupied level）として 6.21 eV が印字されている。これらは「走り切って、刻みのせいでは揺らがないところまで収束した」という、確かな出力である。

ここで、格子定数について手元にそろった数を、素性まで添えて並べてみよう。

道筋 格子定数 \(a\) 素性
手の見積もり（第1章・密度から） 約 5.4 Å 紙と鉛筆で出した
実験値（X 線回折） 5.431 Å 測定で出した
入力に置いた値（celldm(1)=10.26） 5.43 Å 私たちが入力に置いた

三つが 1% の幅に集まっているのは、たしかに心強い。鉛筆だけから当てた 5.4 Å も、実験室の 5.431 Å も、入力に据えた 5.43 Å も、一辺の長さでぴたりと寄り添う。ただし、いちばん下の行だけは、上の二つと同じ目で読まないでおきたい。手の見積もりと実験値は、それぞれの道筋が出してきた答えだ。いっぽう 5.43 Å は、計算が選んだ値ではなく、私たちが初めから置いた前提である。第4章で「celldm(1) はまだ計算の結果ではなく、前提の側にある」と確かめた、あの一行がここで効いてくる。

ここに、本書の背骨がある。入れたもの（入力に置いた前提）と、出てきたもの（計算が返した結果）を、取り違えない。 一点の SCF が返した −93.45 Ry や 6.21 eV は、まぎれもなく計算の出力だ。いっぽう 5.43 Å は、計算の出力ではなく、私たちが据えた前提の側にある。もし「計算がシリコンの格子定数を 5.43 Å と当てた」と書けば、前提を一つ、結果の欄にこっそり書き込むことになる。三つの一致は喜んでよい。そのうえで「三つ目は計算の手柄ではない」と一言添える――この惜しみなさが、検算の正直さである。

3 計算に格子を選ばせるには ― 次の一手

では、計算「自身」にシリコンの格子定数を選ばせるには、どうすればよいか。答えは、もう前々節で出ている。格子定数 \(a\) を何点か振り、全エネルギーが最小になる谷――\(E(a)\) 曲線の底――を探すのだ。入力の言葉でいえば、一点に固定していた celldm(1) を、こんどはつまみとして 5.43 Å の近傍で何通りも回し、それぞれ SCF を収束させて \(E(a)\) を並べ、頂点を読む。これは本書がまだ走らせていない計算であり、ペンを自分の手に持ち替える第11章で、最初の設計課題として正面から取り上げる。

道具立ては、もうできている。底のごく近くで \(E(a)\) がほぼ放物線になることは前節で見た。だから数点を放物線で結べば、頂点として平衡 \(a\) が読める。実際にどう読むのか、仮の数字で一度なぞっておこう。

小課題 8.1

下の表は、方法を示すための仮の代表値である。本書の実行が返した出力ではなく、谷を放物線で読む手つきを練習するための数字だと思ってほしい。ある PBE 計算で \(a\) を振ったところ、最小値からの差がおおよそ次のようになったとしよう。

\(a\) (Å) 5.40 5.45 5.50 5.55
\(\Delta E\) (meV/atom) 9.8 0.8 1.8 12.8

  1. 表のうち、もっともエネルギーが低い \(a\) はどれか。

  2. 底をはさむ 3 点 \(a=5.40,\,5.45,\,5.50\ \mathrm{\text{Å}}\)（間隔 \(h=0.05\ \mathrm{\text{Å}}\)）を放物線で結び、その頂点 \(a_0\) を求めよ。等間隔 3 点 \((E_1,E_2,E_3)\) の放物線の頂点は \[ a_0 = a_2 + \frac{h}{2}\cdot\frac{E_1-E_3}{E_1-2E_2+E_3} \] で与えられる（\(a_2\) は中央の点）。

  3. 得た \(a_0\) を実験値 5.431 Å と並べ、ずれを百分率で言え。

解答例

  1. 表の中では \(a=5.45\ \mathrm{\text{Å}}\) がもっとも低い（\(\Delta E = 0.8\) meV/atom）。だが両隣を見ると、5.45 と 5.50 のあいだはわずか 1.0 meV しか上がらないのに、5.45 から 5.40 へは 9.0 meV も跳ね上がる。底は 5.45 のすぐ右、5.45 と 5.50 のあいだに隠れていそうだ。

  2. 公式に \(a_2=5.45\), \(h=0.05\), \((E_1,E_2,E_3)=(9.8,\,0.8,\,1.8)\) を入れる。分子は \(E_1-E_3 = 9.8-1.8 = 8.0\)、分母は \(E_1-2E_2+E_3 = 9.8-1.6+1.8 = 10.0\)。よって \[ a_0 = 5.45 + \frac{0.05}{2}\cdot\frac{8.0}{10.0} = 5.45 + 0.025\times0.8 = 5.47\ \mathrm{\text{Å}}. \] 表の最低点（5.45）より、底はわずかに右の 5.47 Å にある。

  3. 実験値 5.431 Å と並べると、\(5.47/5.431 = 1.0072\)、すなわち約 0.7% の過大評価。くり返すが、これは仮の表から読んだ谷であって、本書の実行値ではない。それでも向きと大きさ――実験のわずかに上、1% 弱――は、次に見る PBE のくせとよく合う。実際に谷がどこへ落ちるかは、自分で celldm(1) を振って確かめる番だ。

仮の数字とはいえ、谷が実験値のわずかに右、\(5.47\ \mathrm{\text{Å}}\) あたりに来たのは、でたらめではない。私たちが使った PBE という汎関数（GGA の一種）には、よく知られたくせがある。多くの固体で結合をほんの少しやわらかく見積もり、平衡の格子をわずかに広げ気味に返すのだ。文献を見渡すと、PBE はシリコンの格子定数を実験より 1% 前後だけ上に置くことが多い。ひと世代前の LDA は逆向きで、結合を締めすぎて格子を小さめに出す。だから、実際に上の \(E(a)\) スキャンを走らせれば、谷はおそらく \(5.47\ \mathrm{\text{Å}}\) 前後、実験を 0.7% ほど上回るあたりに落ちるだろう――と予告できる。ただしこれは文献既知の傾向であって、本書の実行が返した値ではない。確かめる道はひとつ、自分で celldm(1) を振ることだ。

そして、いざ谷を読んだとき、実験値との残りのずれ（たとえば +0.7%）は、二つの出どころに腑分けできる。ひとつは近似――選んだ汎関数のくせだ。PBE が格子をやや広げるのは、向きも大きさも PBE の指紋どおりで、たまたまの数字ではない。もうひとつは刻み――ecutwfc（第5章）や k 点（第6章）、そして \(a\) を何点で振ったかの粗さである。これらが足りなければ \(E(a)\) 曲線ごとずれ、谷の位置も動く。順番が肝心で、まず刻みを潰し、そのうえで残ったずれを汎関数に帰属する。先に収束を詰めておくからこそ、生き残ったずれを「PBE のくせ」と胸を張って言える。この腑分けを一組の計算で実際に決着させる手立ては、第11章の設計課題に取ってある。

flowchart TB
    IN["私たちが置いた（前提・入力）<br>celldm(1)=10.26 → a = 5.43 Å<br>実験 5.431 Å に近い格子"]
    RUN["固定した格子で<br>SCF を一度回す"]
    OUT["計算が返した（出力・収束済み）<br>全エネルギー ≈ −93.45 Ry<br>最高被占準位 6.21 eV"]
    NEXT["次の一手：a を数点振り<br>E(a) 曲線の谷を探す<br>＝計算自身に格子を選ばせる"]
    IN --> RUN --> OUT
    RUN -. "まだ走らせていない" .-> NEXT
    classDef placed fill:#d4edda,stroke:#28a745,color:#000
    classDef todo fill:#fff3cd,stroke:#ffc107,color:#000
    class IN,OUT placed
    class NEXT todo
図 1: 一回の SCF で、私たちが置いたものと、計算が返したものを分けて読む。格子定数 5.43 Å は入力に置いた前提（左）。収束した全エネルギーと最高被占準位は、計算の出力（右）。計算自身がこの格子を選ぶかは、a を振って E(a) の谷を探して初めて分かる（黄、次の一手）。

4 検算の儀式 ― 理論値・実験値と必ず比べる

ここまでやってきたことを、ひとつの作法として言葉に固定しておきたい。本書ではこれを検算の儀式と呼ぶ。計算機が数字をひとつ返すたび、それを信じる前に、次の四手を踏む。

  1. これは何の量か、単位は何か、そして入力に置いた前提か、計算が返した出力かを言う。
  2. すぐ隣に既知の値を置く（実験値・理論値・第1章の手の見積もり）。
  3. 比と百分率、そしてオーダーで一致を測る。
  4. 一致やずれの出どころを帰属する（前提を映しただけか、近似のくせか、刻みの残差か）。

第3章で確かめたように、全エネルギーは単独では何も語らない。差を取ってはじめて意味が生まれる量だった。格子定数もまた、たった一つの数字としては宙に浮いている。隣にもう一つの数字を置いて、はじめて「どこまで本物か」を語りだす。だから儀式の心臓は、第2手――必ず既知の値の隣に並べる――にある。

この四手を、本章の二つの数字にかけてみよう。まず格子定数。量は長さ、単位は Å、素性は入力に置いた前提だ。隣に実験値 5.431 Å を置くと \(5.43/5.431 \approx 0.9998\) でほぼ完全に重なるが、第1手で素性を「前提」と言い切ったから、この一致を「計算の的中」とは読まない。近い値を置いたのだから、重なって当然である。手の見積もり 5.4 Å（\(-0.6\%\)）まで並べれば、三つが 1% に集まる気持ちよさは確かに味わえる。ただしその一致のいくらかは、私たちが初めから据えた前提の映りだと、正直に割り引いておく。

もう一つは、計算が返した数字――最高被占準位 6.21 eV である。こちらは出力だ。ところが第2手で隣に既知の値を置こうとすると、手が止まる。6.21 eV は、それ単独では比べる相手をもたない。第3章の全エネルギーと同じで、ゼロ点が計算の内側で決まる絶対値だからだ。この数を物理の量――たとえばバンドギャップ――へ向けるには、すぐ上の空準位まで尋ね、その差を取らねばならない。検算の儀式は、ここで「比べる相手がまだいない」という空白そのものを指し示す。その空白の先に、次章の山場が待っている。

私たちはいま、やさしい題材で儀式を一通り回した。同じ四手を、次章ではバンドギャップに対して回す。そのとき返ってくる数字は、これほど素直ではない。やさしい場面で手順を体になじませておくのは、難しい場面が来たときに、自分の期待ではなく手順のほうを信じられるようにするためだ。

よくある誤解

計算値が実験値から数 % ずれていると、「計算が失敗した」「どこかにバグがある」と受け取りたくなる。だが、収束まで詰めたうえでの数 % のずれは、計算が壊れた証拠ではなく、選んだ汎関数の系統的なくせの表れである。PBE が格子をやや過大評価するのは向きも大きさも織り込み済みで、別の汎関数を選べば改善しうる種類のものだ。だから数 % のずれは「失敗」ではなく、その近似の指紋として読める。ずれを失敗と切り捨てず、近似と刻みに帰属できること――それを読み取れる目こそ、ここで養いたいものだ。

5 この章のまとめ

格子定数という最初の標的に、検算の儀式を一周かけ終えた。

  • 計算「自身」に格子定数を選ばせるとは、\(E(a)\) 曲線の底（力・応力がゼロになる間隔）を読むことである。変換器は格子定数を直接返さないので、こちらが最小化の問いに翻訳する。
  • 私たちの実行は、実験に近い 5.43 Å を入力に置いた一点の SCF だった。よく収束し（全エネルギー −93.45 Ry 付近、刻みによる差は meV/atom）、最高被占準位 6.21 eV を返した。これらは出力である。いっぽう 5.43 Å は、計算の出力ではなく、私たちが据えた前提だ。
  • 手 5.4 Å／実験 5.431 Å／置いた 5.43 Å の三つが 1% 以内に集まるのは心強い。ただし三つ目は計算の手柄ではなく、私たちが置いた値である。入れたものと出てきたものを取り違えない――これが本章の核である。
  • 計算自身に格子を選ばせる次の一手は、\(a\) を振って \(E(a)\) の谷を探すこと（第11章）。PBE なら谷は \(5.47\ \mathrm{\text{Å}}\) 前後、実験を 0.7% ほど上回るだろうと予告できる（文献既知）。走らせて初めて、その残差を近似（PBE のくせ）と刻み（収束の残差）に腑分けできる。先に刻みを潰すから、残差を汎関数に帰属できる。
  • 検算の儀式：数字を出したら、まず前提か出力かを言い、既知の値の隣に置き、% とオーダーで一致を測り、出どころを帰属する。本書のすべての量に、この四手をかけていく。

格子定数をめぐる三つの数は、心強い近さで寄り添った。だが本章がいちばん伝えたかったのは、その一致そのものより、一致の読み方だ。出した数字を既知の値の隣に置き、それが前提か出力かを見分ける――この一手さえ手放さなければ、心強さに足をすくわれずに済む。

さて、同じ一回の計算は、もう一つの出力をまだ手元に残している。最高被占準位 6.21 eV――それ単独では比べる相手をもたない、あの数だ。ここですぐ上の空準位まで尋ね、その隔たり、すなわちバンドギャップを読みにいったら、計算は何と答えるだろうか。教科書は約 1.1 eV と言う。格子では指一本ぶんに収まったずれが、この問いでは桁ごと跳ね上がる。

次章の入口で、私たちはその答えの前に立ち尽くすことになる。計算は最後まで走り切り、エネルギーも収束した。格子の側は――入力に置いた 5.43 Å も、振って読めば落ちるはずの 5.47 Å も――どう転んでも実験から 1% と離れない。ところが同じ計算にバンドギャップを尋ねると、読み取れるのは教科書のおよそ半分。壊れた計算が、格子をここまで言い当てるだろうか。何が起きているのか。