Si で学ぶ第一原理計算

QE がやっているのは「結晶の中の電子についてシュレーディンガー方程式を解く」ことです。ただし、そのままでは解けないので、DFT という考え方で問題を小さくし、SCF という繰り返し計算で答えを出します。

1. そのままでは解けない理由：電子が多すぎる

講義で習ったシュレーディンガー方程式 Hψ = Eψ は、電子が1個なら3次元（x, y, z）の問題です。ところが電子どうしは反発し合うので、N 個の電子をまとめて表す波動関数は 3N 個の変数の関数になります。各変数を10点ずつに区切って数値で持つだけでも、必要な数は次のように爆発します。

2. DFT の発想：主役を「電子密度」に替える

DFT（密度汎関数理論）は、波動関数の代わりに電子密度 n(r)（場所ごとの電子の濃さ）を主役にします。電子密度は電子が何個あっても x, y, z の3変数だけの関数です。「一番エネルギーの低い状態（基底状態）は、電子密度だけで決まる」ことが証明されている（ホーエンベルグ・コーンの定理）ので、これで足ります。

実際の計算では、コーン・シャム方程式を解きます。電子どうしの反発を「ほかの電子がつくる平均の場」に置き換えて、電子1個ずつの方程式にしたものです。

V_xc が近似であることが、PBE でバンドギャップが 0.58 eV（実験値 1.12 eV）と小さく出る原因です。スライド6とスライド9の工夫④の「理論の限界」はここにつながります。詳しくは次の「PBE とは何か」で扱います。

3. SCF：答えを使って答えを求める（自己無撞着場）

上の式には「鶏と卵」の関係があります。V_eff を作るには電子密度 n が要り、n を求めるには方程式を解いて ψ を出す必要があり、ψ を出すには V_eff が要ります。そこで、仮の密度から始めて、入れた密度と出てきた密度が一致するまで繰り返します。一致した状態を自己無撞着（self-consistent）と呼び、この計算を SCF（self-consistent field）計算と呼びます。第1回の成果物「計算手順図」は、この流れの図です。

下の図は、手順①〜⑦をブラウザの中で実際に計算しながら1つずつ進めます。Si の3次元の計算は重いので、同じ仕組みの1次元のモデル（長さ10の箱に原子2個・電子4個、原子単位）を解いています。手順ごとに、図の箱・式・コードの行・数値・グラフが一緒に切り替わります。

実際の収束の様子（この PC で計算した Si の SCF）

⑥で新しい密度をどれだけ混ぜるかが mixing_beta です。値を変えて同じ Si の SCF を4回実行し、繰り返しごとの「誤差の目安」（出力の estimated scf accuracy）を並べました。灰色の点線が収束の基準 conv_thr = 1.0d-8 です。

4. 平面波と擬ポテンシャル：計算を現実的な大きさにする工夫

• 平面波基底：波動関数を、いろいろな細かさの波（平面波）の足し合わせで表します。どこまで細かい波まで使うかの上限が ecutwfc（単位 Ry）です。大きいほど正確で、計算は重くなります。どこで十分かを第3回で受講者自身が調べ、30 Ry を選びます。
• 擬ポテンシャル：Si の電子14個のうち、内側の10個（1s, 2s, 2p）は結合にほとんど関わりません。これを原子核とひとまとめにして、外側の価電子4個だけを計算します。2原子で 4×2 ＝ 8個なので、出力には number of electrons = 8.00 と出ます。この講義の擬ポテンシャルはウルトラソフト型（USPP）で、そのため電子密度用のカットオフ ecutrho を ecutwfc の8倍（240 Ry）にしています。

5. Quantum ESPRESSO：計算プログラムの集まり

Quantum ESPRESSO（QE）は、平面波と擬ポテンシャルで DFT 計算をする無償のオープンソースのプログラム群です。1つのソフトではなく、役割ごとのプログラムを順に使います。中心は pw.x で、ほかは pw.x の結果から必要な量を取り出す「後処理」です。

計算の順番：まず pw.x の SCF で電子密度を作る → それを土台に、バンド図（pw.x の bands → bands.x）、DOS・PDOS（pw.x の nscf → dos.x / projwfc.x）、電子密度（pp.x）へ枝分かれします。SCF が全部の出発点なので、第2・3回でその条件をしっかり決めます。

6. 入力ファイルの3つの「段落」と3つの「カード」

入力ファイルは、&名前 で始まり / で終わる段落（ネームリスト）と、大文字のカードでできています。講義の Si.scf.in を分けると次の通りです。結晶構造の行（ibrav など）の意味は図1で詳しく扱っています。

&CONTROL                         ! 何をするか
    calculation = 'scf'          ! SCF 計算（ほかに 'bands', 'nscf' など）
    prefix = 'silicon'           ! 結果ファイルの名前の頭
    outdir = './tmp/'            ! 結果の置き場所（後処理がここを読む）
    pseudo_dir = '../pseudo'     ! 擬ポテンシャルの置き場所
/
&SYSTEM                          ! 何を計算するか（結晶と精度）
    ibrav = 2,  celldm(1) = 10.26,  nat = 2,  ntyp = 1
    ecutwfc = 30.0               ! 平面波の上限（第3回で選ぶ）
    ecutrho = 240.0              ! 電子密度の上限（USPP なので 8 倍）
/
&ELECTRONS                       ! どう解くか（SCF の進め方）
    conv_thr = 1.0d-8            ! 手順⑤の基準（Ry）
    mixing_beta = 0.7            ! 手順⑥の混ぜる割合
/
ATOMIC_SPECIES                   ! 元素・質量・擬ポテンシャル
    Si  28.0855  Si.pbe-n-rrkjus_psl.1.0.0.UPF
ATOMIC_POSITIONS (alat)          ! 原子の位置（図1の Si₁, Si₂）
K_POINTS (automatic)             ! k 点の網目（図2の下の説明）
    6 6 6  1 1 1

7. 出力ファイルで最初に見る行

単位の目安：エネルギー 1 Ry（リュードベリ）＝ 13.606 eV、長さ 1 bohr ＝ 0.5292 Å。QE の入力と出力は、この「原子単位」が基本です。

PBE は、コーン・シャム方程式の中で正確な式が分かっていない「交換相関」の項 V_xc を近似する方法の1つです。1996年に Perdew・Burke・Ernzerhof の3人が発表したので、頭文字をとって PBE と呼びます。電子密度だけでなく電子密度の傾きも使う「GGA（一般化勾配近似）」という種類で、固体の計算で最もよく使われています。

1. そもそも「交換相関」とは

コーン・シャム方程式では、電子どうしの反発を「電荷の雲としての平均」（ハートリー項）で表しました。しかし実際の電子は、平均の雲のように振る舞うだけでなく、お互いを避け合っています。この避け合いの効果が交換相関です。

• 交換：同じ向きのスピンの電子は、パウリの排他原理により同じ場所に来られません。そのため自然に距離をとり、反発が平均より小さくなります（エネルギーが下がる）。
• 相関：スピンの向きに関係なく、電子はクーロン反発によって互いに近づきにくくなります。平均の雲では表せない、この「よけ合い」の残りの部分です。

どちらも、ある電子のまわりに「ほかの電子が少ない領域（交換相関の穴）」ができる効果と考えられます。この穴の正確な形を求めるのは、もとの多電子問題を解くのと同じくらい難しいので、近似が必要になります。

2. 近似の段階：LDA から PBE へ

LDA（局所密度近似）は、「その場所の電子密度 n(r) と同じ濃さの、一様な電子ガス」の交換相関エネルギーをそのまま使う方法です。一様な電子ガスの値は精密に分かっているので、場所ごとの密度だけで計算できます。

ところが、実際の結晶では電子密度は場所によって大きく変わります。原子の近くは濃く、原子の間は薄い。GGA は、密度の傾き ∇n も使って LDA を補正します。PBE の交換エネルギーは次の形です。

F_x(s) は「LDA を何倍に強めるか」を表す増強因子です。密度がほとんど変わらない所（s ＝ 0）では 1 で LDA と同じになり、密度が急に変わる所ほど大きくなります。κ と μ は、合わせたい実験値ではなく、量子力学として満たすべき条件（一様な電子ガスの極限や、交換エネルギーの上限など）から決められています。これが PBE が「特定の物質に合わせ込んでいない」と言われる理由です。相関エネルギーにも同じ考え方の傾き補正が入っています。

3. 近似の「はしご」：PBE はどの段にいるか

交換相関の近似は、使う情報を増やすほど上の段に上がる「はしご」として整理されます（Perdew の「ヤコブのはしご」）。上の段ほど一般に正確ですが、計算は重くなります。

右側のギャップは、PBE のみこの PC で計算した値です。ほかは文献でよく報告される代表的な値の目安で、擬ポテンシャルや計算条件によって 0.1 eV 程度は変わります。

4. PBE でバンドギャップが小さく出る理由

この講義の Si のギャップは 0.57 eV で、実験値 1.12 eV の約半分です。これは計算の失敗ではなく、近似の性質です。理由は大きく2つあります。

• 自己相互作用の誤差：ハートリー項には、電子が自分自身の電荷の雲と反発する分が含まれています。正確な交換ならこれをぴったり打ち消しますが、LDA や PBE では打ち消しきれません。その結果、電子は実際より広がりすぎ、占有された状態（価電子帯）のエネルギーが持ち上がって、ギャップが狭まります。
• コーン・シャムのギャップは、もともと本当のギャップではない：コーン・シャム方程式は「基底状態の電子密度」を正しく出すための仕組みで、その固有値の差は、電子を1個出し入れするときの本当のエネルギー差と同じではありません。正確な交換相関なら、電子数が整数をまたぐ所で V_xc がとびを持ち、その分がギャップに加わります。LDA や PBE はなめらかな関数なので、このとびが入りません。

この過小評価は Si だけでなく、半導体・絶縁体に共通の傾向です。たとえばゲルマニウムでは PBE のギャップがほぼ 0 になり、ダイヤモンドでも実験値 5.5 eV に対して 4 eV 台にとどまります。

5. PBE で格子定数を求めるとどうなるか（この PC で計算）

講義の入力は、格子定数に実験値（celldm(1) = 10.26 bohr ＝ 5.43 Å）を使っています。PBE で一番エネルギーが低くなる格子定数を調べるため、格子定数を6通り変えて SCF 計算をしました（30 Ry、6×6×6）。

6. PBE が得意なこと・苦手なこと

7. 入力と出力のどこに PBE が出てくるか

• 入力：入力ファイルには PBE と書く行はありません。QE は擬ポテンシャルのファイルに書かれた汎関数を使います。Si.pbe-n-rrkjus_psl.1.0.0.UPF の名前の pbe がその印で、ファイルの中にも functional="PBE" と書かれています。擬ポテンシャルは作るときに汎関数を決めているので、汎関数を変えたいときは、その汎関数で作った擬ポテンシャルに替えるのが原則です（input_dft で上書きもできますが、組み合わせがずれるので勧められません）。
• 出力：Exchange-correlation= PBE の次の行の ( 1 4 3 4 0 0 0) は、内訳の番号です。順に、交換の LDA 部分（スレーター）、相関の LDA 部分（Perdew–Wang）、交換の傾き補正（PBE）、相関の傾き補正（PBE）を表します。

格子定数の計算の入力と結果は 配布物の qe/04_lattice/ に保存しました。

基本的にはテキストエディタで書くファイルです。QE 本体には入力を作る画面がありません。ただ、白紙から書く人はほとんどいません。既存の例をコピーして、必要な行を書き換えるのが普通です。

よく使われる作り方

• 例をコピーして書き換える：QE に付属している例題、チュートリアル、論文の付録などにある入力を元にします。Si は定番の例なので、講義の入力とほぼ同じものが多く出回っています。
• 結晶構造ファイルから変換する：別の材料を扱うときは、CIF などの構造ファイルから ASE や pymatgen（Python ライブラリ）で入力を書き出します。ASE は pip install ase で入れられます。
• Web やGUIのツールを使う：Materials Cloud の入力生成ツールや、BURAI のような GUI があります。k 点経路だけを自動で決めるツール（SeeK-path など）もあります。

この講義での扱い

講義では、各回のフォルダ（02_scf_crystal、04_band_calc など）に入力ファイルが配布されています。受講者は中身を読み、ecutwfc や k 点の数などの一部を書き換えて実行します。第3回の収束性の調査がその例です。白紙から書くことは前提にしていません。

課題①の対策「結晶構造の図と入力ファイルを1対1で対応させる」は、配布された入力の各行が、構造のどこに対応しているかを図で確かめる、という意味です。下の図1と図2がちょうどその対応になります。

入力ファイルの行との対応

ibrav = 2           ! 面心立方（fcc）の基本格子を使う ＝ 図の琥珀色の箱
celldm(1) = 10.26   ! 灰色の立方体の一辺 a（bohr）。10.26 × 0.5292 ≒ 5.43 Å
nat = 2             ! 箱の中の原子 ＝ 大きく描いた 2 個
ATOMIC_POSITIONS (alat)
    Si  0.00  0.00  0.00   ! Si₁：立方体の角
    Si  0.25  0.25  0.25   ! Si₂：対角線を 1/4 進んだ位置（Si₁ と結合）

補足：図の基本格子は、立方体の中に収まって見やすい3本のベクトル (0,½,½)a、(½,0,½)a、(½,½,0)a で描いています。QE の ibrav = 2 が内部で使うベクトルは (−½,0,½)a などで、向きの取り方が少し違います。ただし表している格子はまったく同じです。

入力の数値と図の位置の対応

K_POINTS (crystal_b) の数値は、基本格子に合わせた斜めの軸で測った座標です。教科書は立方体の軸（直交座標）で書くので、同じ点でも数値が違って見えます。

変換は「入力の3つの数 × QE の逆格子ベクトル b₁=(−1,−1,1)、b₂=(1,1,1)、b₃=(−1,1,−1) の和」です。例：X は 0×b₁ + 0.5×b₂ + 0.5×b₃ = (0, 1, 0)。図の点にカーソルを合わせると、両方の座標が表示されます。

k 点はもう1種類ある

SCF 計算の K_POINTS (automatic) 6 6 6 は、この立体の中を6×6×6の網目で均等に調べる指定です。全エネルギーや電子密度を求めるためのもので、経路ではありません。バンド図用の crystal_b は、上の経路の線上だけを調べます。

図1と図2の「結晶方位」ボタンで、方位の矢印を重ねて表示できます。濃い矢印が代表の1本、淡い線が同じ仲間（等価な方位）です。2つの図で同じボタンが連動します。

記号の読み方

• [100]：角かっこは特定の1方向。立方体の x 軸の向き。
• <100>：山かっこは等価な方向の仲間全部。[100]、[010]、[001] と、その逆向きを含めて6方向。<110> は12方向、<111> は8方向。
• (100)：丸かっこは面（[100] に垂直な面）。「Si(100) ウェハ」の (100) はこちら。

立方晶では、実空間の方向 [hkl] と k 空間の方向 (h, k, l) が同じ向きになります。このため、バンド図の区間を結晶方位の言葉で呼べます。

図2の代表矢印は、講義の経路で使っている X (0,1,0)、K (0,¾,¾)、L (−½,½,½) に向けています。図1の代表矢印は立方体の角 Si₁ から [100]、[110]、[111] に向けています。図1で <111> を表示すると、Si₁ から Si₂ への結合が [111] の矢印の上に乗っていることがわかります。

方位を正面に向けて眺める

「正面に向ける」の [100]・[110]・[111] を押すと、その方位がまっすぐ手前を向くように図が回ります。図1と図2は同時に同じ向きへ回ります。正面に向けている間は、奥行きで形がゆがまない投影（平行投影）に切り替え、座標軸を隠します。同じ結晶でも、見る向きによって原子の並びがまったく違って見えます。実空間は「2×2×2 に広げる」をオンにすると、並びの模様がわかりやすくなります。

実空間と k 空間で、同じ向きから見ると同じ回転対称性（4回・2回・3回）が現れます。図1と図2を見比べると確かめられます。

計算には「細かさ」のつまみが2つあります。波動関数をどこまで細かい波で表すか（ecutwfc）と、k 空間をどれだけ細かい網目で調べるか（K_POINTS）です。上げるほど正確になり、計算は重くなります。そこで、上げても全エネルギーがほとんど変わらなくなる（±2 meV/atom 以内に収まる）一番小さい値を選びます。これが収束性の確認です。

1. カットオフ ecutwfc：何個の波で表すか

QE は波動関数を平面波（いろいろな細かさの波）の足し合わせで表します。細かい波ほど運動エネルギーが大きいので、「運動エネルギーが ecutwfc 以下の波だけを使う」という形で、使う波の数を決めます。下のデモは、1次元の山形の関数を、細かさの上限を決めた波の足し合わせで表すものです（実際に計算しています）。

2. k 点の網目：ブリルアンゾーンを何点で調べるか

全エネルギーや電子密度は、ブリルアンゾーン全体（図2の立体）の足し合わせ（積分）で決まります。実際には網目の点で代表させて足し合わせます。K_POINTS (automatic) 6 6 6 1 1 1 は、各方向を6等分した網目（6×6×6 ＝ 216 点）を半目盛りずらして置く指定です。ただし対称性で同じになる点はまとめられるので、実際に計算するのは 28 点だけです（出力の number of k points= 28）。Si は絶縁体でエネルギーが k に対してなめらかに変わるので、粗い網目でもすぐに正確になります。金属では、何倍も細かい網目が必要になります。

3. この PC で実行した収束性の計算

講義と同じ条件（PBE・USPP、ecutrho ＝ 8 × ecutwfc）で、カットオフを9通り（k 点 6×6×6 固定）、k 点を9通り（30 Ry 固定）変えて SCF 計算をしました。縦軸は、全エネルギーの差を1原子あたりの meV に直したものです（Si は2原子なので2で割る。1 Ry ＝ 13.606 eV）。

「平面波の数の目安」は出力の Smooth grid: … G-vectors の値で、ecutwfc の 3/2 乗に比例します（30 Ry を基準にすると、60 Ry で 2^3/2 ＝ 2.83 倍。表の値もちょうど 2.83 倍）。カットオフを2倍にすると波の数は約2.8倍になり、計算はそれ以上に重くなります。

4. 読み取れること

• カットオフ：上げるほど全エネルギーは下がり続け、ある値で変化が止まります（平面波が増えるほど答えが真の値に近づく、という性質があるため、一方向に近づく）。±2 meV/atom の基準では 30 Ry（差 1.77 meV/atom）が最小で、25 Ry（5.08 meV/atom）では足りません。スライド5の結論と一致します。
• k 点：網目を細かくすると差が急速に小さくなります（2→3→4 で 81 → 9.1 → 1.3 meV/atom）。この計算では 4×4×4 の時点で ±2 meV/atom に入っています。講義が 6×6×6 を採用しているのは、差が 0.04 meV/atom と十分に余裕があり、計算する k 点も28点で負担が小さいからです。後の DOS 計算では、さらに細かい網目（12×12×12）を使います。
• 差の取り方：「一番細かい条件との差」と「1つ前との差」で、選ばれる値が変わることがあります。ボタンで切り替えて確かめられます。

計算の入力と結果は 配布物の qe/03_convergence/ に保存しました。講義の 03_convergence の値（50 Ry までの5点）とも、同じ条件の点では全エネルギーが8桁一致します（例：30 Ry・6×6×6 で −22.83957491 Ry）。

できます。ただし k が3次元なので、エネルギーまで含めると4次元になり、そのままでは描けません。k 空間の平面を1枚選び、その上のエネルギーを高さにした曲面として描くのが現実的です。

下の図は、Γ・X・K を含む平面（k_z = 0）の上で、価電子帯の一番上のバンドと伝導帯の一番下のバンドを描いたものです。

立体を重ねて表示する意味はあるか

• 結晶（原子の立体）を重ねる意味はありません。この図の床は k（電子の波の向きと細かさ）で、高さはエネルギーです。原子の位置（実空間）とは別の空間で、単位も違います。重ねると「原子の真上でエネルギーが高い」のような誤った読み方を招きます。
• ブリルアンゾーンの立体を丸ごと重ねるのも適しません。高さ方向は k_z ではなくエネルギーだからです。ゾーンとの関係は、床の八角形（ゾーンの断面）として描いています。
• 意味があるのは「並べて連動させる」ことです。ブリルアンゾーンのどの断面を見ているのか、2D のバンド図のどこに当たるのかを、同じ k 点で同時に示します。下の図で、点が経路をたどる様子を確かめられます。

検討結果

計算の方法：k_z = 0 平面の 1/4（k_x, k_y ≥ 0）を K_POINTS (tpiba) で 961 点並べて calculation = 'bands' を実行し、立方対称性（E(±k_x, ±k_y) が等しい）で全体に広げました。

バンド図は「どの k で、どのエネルギーか」を見る図でした。DOS（状態密度）は k をまとめてしまい、「あるエネルギーに状態がいくつあるか」だけを数えた図です。PDOS（部分状態密度）は、その状態を s 軌道・p 軌道などの成分に分けたものです。

1. DOS は「ブリルアンゾーン全体のエネルギーを数えたヒストグラム」

ブリルアンゾーンの中に網目状に k 点を置き、各 k 点の8本のバンドのエネルギーを全部集めて、エネルギーごとに数えます。数えるときは、各エネルギーに幅を持たせた山（ガウス関数）を置いて足し合わせます。下の図は、QE の nscf 計算で出したエネルギーを、ページの中で実際に数えたものです。

2. バンド図と DOS・PDOS を並べる

同じエネルギーの高さで、左にバンド図、右に DOS・PDOS を並べました。縦の位置をそろえてあるので、カーソルを合わせると横に点線が伸びます。

3. 読み取れること

• DOS がゼロの範囲 ＝ バンド図のギャップ：0 eV（VBM）から約 0.6 eV まで、状態が1つもありません。バンド図で Γ の VBM と Γ–X の CBM の間に線がない範囲と同じです。
• 深い価電子帯（−12〜−8 eV）は 3s が主：この範囲の状態の約7割が s 成分です。バンド図では、Γ 点の一番下のバンド（−12.0 eV）が純粋な s（100%）です。
• VBM 付近（−4〜0 eV）は 3p が主：この範囲の約95%が p 成分です。Γ 点の VBM は3本が重なっていて（縮退）、3本とも純粋な p です。これは p_x・p_y・p_z の3つに対応します。
• 伝導帯も p が主だが、s も混ざる：0.5〜4 eV の範囲は p が約7割です。CBM（Γ–X の約85%の位置）の状態は s : p ＝ 38 : 62 で、s もかなり混ざっています。
• sp³ 混成の数値的な証拠：価電子帯全体で足し合わせると、s 成分が 2.32 個、p 成分が 5.60 個です。1原子あたりでは s 1.2 個・p 2.8 個で、ほぼ s¹p³ の形になっています。孤立した Si 原子は 3s²3p²（s 2個・p 2個）なので、結晶になると s の電子が約1個 p に移り、1つの s と3つの p が混ざった4本の結合（sp³ 混成軌道）をつくっていることが分かります。中ほどの価電子帯（−8〜−4 eV）で s と p が同じくらい混ざって現れるのも、この混成の表れです。

4. QE での手順

• pw.x（scf）：6×6×6 で電子密度を作る。
• pw.x（nscf）：calculation = 'nscf'、occupations = 'tetrahedra'、K_POINTS 12 12 12 1 1 1。電子密度は固定したまま、細かい網目でエネルギーだけ計算し直す。
• dos.x：DOS を出す（silicon.dos.dat、列は エネルギー・DOS・積分した数）。
• projwfc.x：各状態を原子の軌道（3s, 3p）に射影して PDOS を出す（silicon.pdos.pdos_atm#1(Si)_wfc#1(s) など、原子ごと・軌道ごとのファイル）。

projwfc.x は、結晶の状態を「原子の 3s と 3p の組み合わせ」で表そうとします。表しきれなかった割合が出力の Spilling Parameter で、この計算では 0.0091（約1%） です。小さいほど、s・p だけで状態がよく表せていることを意味します。左の図のファットバンドは、バンド計算の後に projwfc.x を実行して求めました（lsym = .false.）。

電子密度 ρ(r) は、SCF 計算の答えそのもの（第1回の手順⑦で出てくる密度）です。そこから「孤立した原子の密度を重ねただけのもの」を引くと、結晶になったことで電子がどこに移動したかだけが残ります。これが差電子密度 Δρ です。

1. 断面で見る：結合を含む面・水平の面・結合に垂直な面

Si–Si 結合は <111> 方向を向いているので、結合をまっすぐ含む (1−10) 面で切ります。横が [110] 方向、縦が [001] 方向で、Si 原子がジグザグにつながった鎖が見えます（図1で [110] から見た並びの断面です）。「断面」のボタンで、水平の (001) 面と、結合に垂直な (111) 面にも切り替えられます。右のグラフ（結合に沿った密度）はどの断面でも同じです。

2. 3次元で見る（スライド8 と同じ等値面）

Δρ が同じ値になる面（等値面）を、立方体（慣用格子）の中に描きました。スライド8 と同じく、値を 0.003 → 0.007 e/bohr³ と上げていくと、弱い成分が消えても結合の中央の電子のかたまりが残ります。

3. QE での手順

• pw.x（scf）：電子密度を作る（講義と同じ 6×6×6）。
• pp.x：plot_num = 0 で ρ、plot_num = 9 で Δρ（原子の密度の重ね合わせを引いたもの）を取り出す。iflag = 3（3次元）、output_format = 5（XSF）または 6（cube）で、VESTA に読み込める形にする。
• VESTA：等値面の値を変えて表示し、(110) 面などの断面も作れる。

このページの図は、講義の 07_charge_density と同じ入力で pp.x を実行した cube ファイル（基本格子の 36×36×36 の網目）を、立方体と (1−10) 面の網目に補間して描いています。電子密度を全体で足し合わせると 8.000 個（価電子 4 個 × 2 原子）、差電子密度を足し合わせるとほぼ 0 になり（電子が移動しただけで、増えも減りもしない）、どちらも定義どおりです。擬ポテンシャルを使っているので、ここでの電子密度は価電子だけのもので、内殻の 1s・2s・2p 電子は含みません。