6  波をどこまで細かく刻むか

電子の波を、計算機は有限個の波の重ね合わせで表す。
では、どこまで細かい波まで数えれば「足りた」と言えるのか。

この章のゴール:ecutwfc を「波の解像度」として理解し、収束を近似誤差の制御として説明でき、収束表を自分の手で読めるようになる。

第4章で、入力ファイルの各行を物理の供述書として読み終えた。ただ一か所、値を空けたまま先送りにした約束がある。ecutwfcK_POINTS――「どこまで細かく刻むか」を決める二つである。この章で、その前半 ecutwfc にけりをつける。

止め時をどう決めるか、結論の形だけ先に言っておこう。判断のよりどころは、エネルギーの絶対値ではない。刻みを上げたときの変化(差)である。第3章で、全エネルギーは差にだけ意味があると確かめた。ここでもまた、差が主役になる。

1 電子を波の重ね合わせで表す

計算機の中で電子の状態を扱うには、まず波動関数 \(\psi(\mathbf{r})\) を有限のデータとして持たせなければならない。素直なやり方のひとつが、平面波の重ね合わせで表すことである。音をいくつもの正弦波(基音と倍音)の足し合わせで作れるのと同じように、結晶の中の電子の波も、波数のそろった単純な波をたくさん足し合わせて組み立てる。

ここで言う単純な波――平面波――は \(e^{i\mathbf{G}\cdot\mathbf{r}}\) という形をしていて、波数ベクトル \(\mathbf{G}\) がその波の細かさを決める。波長は \(\lambda = 2\pi/|\mathbf{G}|\) だから、\(|\mathbf{G}|\) が大きいほど波長は短く、空間で細かく振動する。少ない波だけを足せば、ゆるやかにうねる形しか作れない。細かい波を足すほど、鋭い凹凸まで再現できる。ちょうど、解像度の低い絵に細部を描き込んでいく作業に似ている。

この「細かさ」に、エネルギーという値段がついているのが要点である。波数 \(\mathbf{G}\) の平面波がもつ運動エネルギーは

\[ E_{\mathrm{kin}} = \frac{\hbar^2 |\mathbf{G}|^2}{2m} \]

で与えられる。\(|\mathbf{G}|\) が大きい波ほど、すなわち細かく振動する波ほど、運動エネルギーが高い。これは難しい話ではなく、波長の短い波ほど勢いが要る、という日常の感覚そのままである。

そして Quantum ESPRESSO が使う Ry という単位で運動エネルギーを測ると、この関係はとても見やすくなる。\(|\mathbf{G}|\) を bohr⁻¹ で測れば、平面波のエネルギーはちょうど \(|\mathbf{G}|^2\)(単位 Ry)になる。「波の細かさ」と「エネルギー」が、二乗ひとつで直結しているのである。

2 ecutwfc は波の細かさの上限

波動関数を厳密に表そうとすれば、いくらでも細かい波――無限個の \(\mathbf{G}\)――が要る。計算機はそれを扱えない。そこで、ある高さで線を引く。運動エネルギーがある上限を超える波は、もう数えないと決めるのである。この上限が ecutwfc(波動関数のエネルギーカットオフ、wavefunction の cutoff)で、Ry で与える。条件を式で書けば

\[ \frac{\hbar^2 |\mathbf{G}|^2}{2m} \le E_{\mathrm{cut}} \quad\Longleftrightarrow\quad |\mathbf{G}|^2 \le \texttt{ecutwfc} \]

であり、右側は先ほどの Ry 単位の言い換えである。図形で言えば、波数空間に半径 \(\sqrt{\texttt{ecutwfc}}\) の球を描き、その球の中に入る平面波だけを使う、ということだ。ecutwfc を上げれば球は大きくなり、数える波が増え、表せる形は細かくなる。ecutwfc とは、まさに「どこまで細かい(高エネルギーの)波まで数えるか」の上限なのである。

この線引きが、いくらでも上げてよいものではないことも、同じ式から見える。球の中の波の本数は球の体積に比例し、体積は半径の3乗、すなわち \(\texttt{ecutwfc}^{3/2}\) に比例する。カットオフを2倍にすれば波の本数はおよそ \(2^{3/2}\approx 2.8\) 倍に膨らむ。解像度を上げる代金は、思いのほか急に高くなる。

では、必要な刻みの細かさは何が決めるのか。波動関数がもつ最も鋭い空間変化である。原子と原子のあいだの結合領域では、電子の波はゆるやかにうねっていて、低い \(\mathbf{G}\) の波だけで足りる。鋭い変化が起きるのは、原子核のごく近くだ。第4章で「便利だが約束つき」と位置づけた擬ポテンシャルは、まさにこの芯の近くを意図的になめらかにならし、本来なら高い \(\mathbf{G}\) まで動員しなければならなかった部分を肩代わりしてくれる。だからシリコンは、数百 Ry ではなく数十 Ry という穏やかな ecutwfc で間に合う。必要なカットオフは、物質と擬ポテンシャルが二人で決める、と言ってもよい。

では、ecutwfc を 30 Ry と置くと、空間ではどれくらい細かい構造まで描けるのか。まず見積もってみよう。上限の波数は \(|\mathbf{G}|_{\max} = \sqrt{30} \approx 5.48\ \mathrm{bohr^{-1}}\)\(1\ \mathrm{bohr} = 0.529\ \mathrm{\text{Å}}\) を使って Å⁻¹ に直すと \(|\mathbf{G}|_{\max} \approx 5.48 / 0.529 \approx 10.4\ \mathrm{\text{Å}^{-1}}\)。これに対応する最短波長は

\[ \lambda_{\min} = \frac{2\pi}{|\mathbf{G}|_{\max}} \approx \frac{6.28}{10.4} \approx 0.60\ \mathrm{\text{Å}}. \]

つまり 30 Ry の基底は、約 0.6 Å の細かさの凹凸まで表せる。シリコンの隣り合う原子の間隔(結合長)は \(a\sqrt{3}/4 = 5.431\times1.732/4 \approx 2.35\ \mathrm{\text{Å}}\) だから、原子間隔のおよそ4分の1まで刻めることになる。原子と原子のあいだに広がる結合の電子雲を描くには十分に細かく、しかも余裕がある。30 Ry という値が乱暴な設定ではないことが、計算機を回す前に、長さの感覚として確かめられた。

3 収束 ― どこで止めてよいか

ecutwfc を上げると全エネルギーはどう動くか。下がっていって、やがて頭打ちになる。これは経験則ではなく、理由のある振る舞いである。

鍵は、計算がエネルギーを「下げられるだけ下げた値」を返すことにある。基底(数える波の集合)を広げれば、波動関数をこれまでと同じか、それ以上にうまく表せる。表現力が上がった分、エネルギーを最小化する余地も広がる。だから ecutwfc を上げて損をすることはなく、全エネルギーは下がるか、変わらないかのどちらかにしかならない。上限の波数を増やしていくと、エネルギーは真の値(選んだ近似の枠内での真の値)へ、つねに上から近づいていく。これが「下がって頭打ち」の正体である。

ありがたいのは、近似の誤差が片側に出ると分かっていることだ。いまの値は真の値より必ず高い側にあり、刻みを細かくすればその差は縮む。誤差がどちらを向いているか分かっていれば、それは制御できる誤差である。

そこで判断はこうなる。頭打ちの極限を執拗に追う必要はない。第3章で見たとおり、全エネルギーの絶対値そのものには意味がなく、しかもその絶対値は擬ポテンシャルの取り方しだいで変わってしまう。私たちが本当に欲しいのは、格子定数を変えたときのエネルギー差や、ばらばらの原子と結晶の差――要するにである。だから止め時も、関心のある差が動かなくなったところでよい。

さらに追い風がある。同じ ecutwfc どうしで差をとると、基底が足りないことによる誤差の大部分が引き算で打ち消し合う。その結果、差は絶対値よりも早く収束することが多い。第3章で「比べるときは土俵をそろえよ」と言ったことの、目に見える御利益である。絶対値を極限まで詰めなくても、差の精度は先に確保できる。

4 収束表を読む:エネルギー差で判断する

言葉だけでは止め時は決まらない。表で読む。下に、シリコンの全エネルギーを ecutwfc を変えながら実際に求めた収束表を示す。計算条件は Quantum ESPRESSO(PWSCF v7.5)、k点は \(8\times8\times8\) に固定、シリコンの2原子セル、擬ポテンシャルは PAW 型の Si.pbe-n-kjpaw_psl.1.0.0.UPF である。電荷密度のカットオフ ecutrho は、後で述べる目安どおり ecutwfc の4倍にそろえてある。

ecutwfc (Ry) ecutrho (Ry) 全エネルギー (Ry) 隣接差 (meV/atom) 80 Ry基準の差 (meV/atom)
40 160 −93.45288506 0.728
45 180 −93.45291006 0.170 0.558
50 200 −93.45294415 0.232 0.326
60 240 −93.45296728 0.157 0.169
70 280 −93.45298074 0.092 0.077
80 320 −93.45299208 0.077 0.000

第3列の全エネルギーは2原子セルの値で、差は原子あたりに直してある(セルの差を2で割る)。まずその全エネルギーに目をやると、ほとんど驚くほど動かない。40 Ry でも 80 Ry でも、小数第3位までは同じ −93.452 Ry で、動くのは第4位以下のごく小さな桁だけである。全域での変化はおよそ 0.0001 Ry――絶対値のおよそ100万分の1にすぎない。第3章で確かめたとおり、全エネルギーの絶対値そのものは擬ポテンシャル次第の数字で、単独では何も語らない。動くのは小数以下の差だけだという、この章の冒頭の予告が、そのまま目の前にある。

読むべきは差のほうだ。第4列は刻みを一段上げたときの隣接差、第5列は最も細かい 80 Ry を仮の到達点として、そこからの隔たりをとったものである。隣接差は 0.170, 0.232, 0.157, 0.092, 0.077 meV/atom と、どれも 1 meV/atom をはるかに下回り、おおむね縮んでいく。ecutwfc がいかに速く収束するかが、ひと目で見てとれる。40 Ry の時点で、最も細かい 80 Ry の値からの隔たりは、わずか 0.728 meV/atom しかない。

ここに止め時のものさしを当てる。古典的には「全エネルギーで 1 mRy/atom」という緩い目安があり、\(1\ \mathrm{mRy} = 13.6\ \mathrm{meV}\) だから原子あたり 14 meV 程度にあたる。だがこの擬ポテンシャルでは、その線は 40 Ry の時点でとうに下回っている。研究で使うなら、もう一桁締めて meV/atom の単位で考えるとよい。許容を 1 meV/atom とすれば、第5列が示すとおり 40 Ry ですでに 0.728 meV/atom――最小の刻みで、もう収束値から 1 meV/atom 以内にいる。これは偶然ではない。使った PAW 擬ポテンシャルの推奨 ecutwfc は約 44 Ry であり、40〜50 Ry でほぼ足りるのは筋が通っている。

安全側に倒すなら、許容をさらに締めればよい。0.5 meV/atom を基準にすると、第5列で 50 Ry が 0.326 で下回り、45 Ry は 0.558 でまだ上にある。よって 50 Ry。0.2 meV/atom まで締めれば 60 Ry(0.169)まで上げることになる。落とし所は 40〜60 Ry――どこで止めるかは、求める精度が決める。70 Ry や 80 Ry まで積んでも、買えるのは 0.1 meV/atom 前後の上積みであり、\(\texttt{ecutwfc}^{3/2}\) で増える波の本数という代金には見合わない。

収束は一発で決まる設定ではなく、手続きである。低めの刻みから始め、一段上げて差を測り、差が許容を下回ったら手前で止める。この流れを図にしておく。

flowchart TB
    A["低めの ecutwfc から始める"] --> B["一段上げて再計算し<br>差 ΔE を測る"]
    B --> C{"ΔE は許容を<br>下回ったか<br>(例 1 meV/atom)"}
    C -- "いいえ:まだ動く" --> B
    C -- "はい:止まった" --> D["その手前の刻みで<br>収束とみなして採用"]
図 1: 収束は手続きである。刻みを上げ、差を測り、差が許容を下回ったところで止める。

この読み方は、自分の計算でも同じように使える。手元のシリコンで ecutwfc を 40, 45, 50, … と振り、出てきた全エネルギーの差をこの目で確かめてみてほしい。絶対値が手元で何 Ry になるかは擬ポテンシャル次第だが、「差が縮んで頭打ちになる」という形は、必ず再現するはずだ。

小課題 5.1

別のシリコン計算で、ecutwfc を変えながら全エネルギー(原子あたり、eV 単位)を求めたところ、代表値として次が得られたとする。

ecutwfc (Ry) 全エネルギー /atom (eV)
40 −215.0000
45 −215.0014
50 −215.0021
55 −215.0025
60 −215.0027
65 −215.0028
  1. 前段からの隣接差を meV/atom で求めよ(\(1\ \mathrm{eV} = 1000\ \mathrm{meV}\))。

  2. 全エネルギーを 0.5 meV/atom の精度で収束させたい。採用すべき ecutwfc はいくらか。

  3. 結合長の比較など、もっと粗い 1 meV/atom で足りる用途なら、ecutwfc はどこまで下げられるか。(2) と何が変わったかを一言で述べよ。

解答例

  1. 前段との差をとる(eV の差を1000倍して meV にする):
区間 隣接差 (meV/atom)
40→45 −1.4
45→50 −0.7
50→55 −0.4
55→60 −0.2
60→65 −0.1
  1. 「それ以降の差がすべて許容を下回る、最も小さい刻み」を探す。0.5 meV/atom を基準にすると、50→55 は 0.4(<0.5)だが、45→50 は 0.7(>0.5)でまだ動く。よって採用は 50 Ry

  2. 1 meV/atom なら、45→50 が 0.7(<1)で、それ以降の差もすべて 1 を下回る。一方 40→45 は 1.4(>1)。よって 45 Ry まで下げられる。求める精度をゆるめると、必要な刻みは粗くなり、計算は軽くなる――止め時は精度の要求と一対である。

よくある誤解

ecutwfc は大きいほど安心だから、迷ったら大きく取ればよい」と考えたくなる。確かに大きくして精度が悪くなることはない。だが代金を思い出してほしい。波の本数は \(\texttt{ecutwfc}^{3/2}\) で増え、計算時間はそれ以上の勢いで膨らむ。収束した後の上積みは、自分が決めた許容より小さい変化でしかない。必要な精度で止める判断こそが腕の見せどころであって、やみくもに上げることではない。

もうひとつ、忘れてはならない土俵の話がある。二つの計算を比べてエネルギー差を読むなら、両者の ecutwfc をそろえること――これは大きさそのものより大切だ。差で見れば誤差は打ち消し合うので、「同じ刻みで比べる」かぎり、そこそこの ecutwfc でも信用できる差が得られる。

なお入力にはもう一つ、電荷密度を表す波のカットオフ ecutrho がある。密度は波動関数の二乗だから、より細かい波まで含み、目安として ecutwfc の数倍――PAW や超ソフト型(US)では 4 倍――に取る。先の収束表でも 40→160、60→240 と、つねに4倍にそろえてあった。これは ecutwfc と一緒に決まる相棒であって、単独で吊り上げるつまみではない。

5 この章のまとめ

ecutwfc の正体と、その決め方が手に入った。

  • 計算機は電子の波を平面波の重ね合わせで表す。細かく振動する波ほど運動エネルギーが高く、ecutwfc は「どこまで細かい波まで数えるか」の上限(波数空間の球の半径の二乗)である。
  • ecutwfc を上げると全エネルギーは上から単調に下がって頭打ちになる。これは基底を広げれば表現力が上がるという変分の理屈による、片側の制御できる誤差である。
  • 止め時は絶対値ではなくで読む。「それより上げても差が許容(研究なら 1 meV/atom、安全側で 0.5 meV/atom 程度)を超えなくなる、最も小さい刻み」を採る。同じ刻みで差をとれば誤差は打ち消し合い、差は絶対値より早く収束する(第3章の御利益)。
  • シリコンの実測では 40〜60 Ry が落とし所になり、使った PAW 擬ポテンシャルの推奨値(約 44 Ry)とも整合する。なお 30 Ry でも実空間で約 0.6 Å、結合長 2.35 Å のおよそ4分の1を刻めるから、幾何学的には数十 Ry で十分に細かい。ただしこの数字は暗記する量ではなく、収束表で各自が確かめる量である。

これで保留の前半は片付いた。残るは K_POINTS である。実空間で「波をどこまで細かく刻むか」は決着したが、結晶にはもう一つの連続が潜んでいる。周期的な結晶の電子は、連続に広がる波数をもつ。その連続を、いくつの点で代表させれば、答えはぐらつかないのか。次章で、もう一つの「刻み」に向き合う。