第15章擬ポテンシャル

前章では，コーン・シャム方程式を計算機上で解くための枠組み—基底関数の選択，一般化固有値問題，自己無撞着ループ—を学んだ．しかしそこでは，原子核が作る $-Z/r$ という深いクーロンポテンシャルと，その底に沈んでいる内殻電子をどう扱うかという問題を保留にしていた．本章はこの問題に正面から取り組む．内殻電子は化学結合にほとんど関与しないにもかかわらず，その存在は価電子の波動関数に細かい節構造を強制し，平面波展開の収束を絶望的に悪化させる．そこで「内殻を消去し，核+内殻の効果を価電子だけが感じる有効ポテンシャル—擬ポテンシャル—に置き換える」という戦略が生まれた．本章では，この置き換えがなぜ許されるのか(散乱理論とノルム保存)，どう構成するのか(Troullier–Martins法，Kleinman–Bylander分離形，ウルトラソフト化とそのノルム保存化)，そしてどう検証するのか(対数微分・バンド構造の全電子計算との比較)を，途中の計算を省かずに追う．次章では，こうして固まった計算の枠組みの上で，原子に働く力と応力を計算し，構造最適化と第一原理分子動力学へ進む．

この章で学ぶこと

内殻電子の凍結が正当化される理由と，全電子計算で1s軌道を平面波展開するのに必要なカットオフの定量的見積り
OPW法からPhillips–Kleinman変換への道筋:芯状態への射影がつくる反発ポテンシャルとキャンセレーション定理
散乱理論から見た擬ポテンシャルの合格条件:対数微分の一致と移植性(transferability)
本章の核心:対数微分のエネルギー微分と波動関数のノルムを結ぶ恒等式のWronskianによる完全導出，およびノルム保存条件の根拠
Troullier–Martins構成法とスクリーニング解除，Kleinman–Bylander分離形とゴースト状態，ウルトラソフト擬ポテンシャルとそのノルム保存化(MBK法)
相対論的擬ポテンシャル(スカラー相対論とスピン軌道)と，作った擬ポテンシャルを検証する実際の手続き

本章でもHartree原子単位系($\hbar = m_e = e = 4\pi\varepsilon_0 = 1$)を用いる．エネルギーの単位は Hartree(1 Ha $=$ 27.211 eV)，長さの単位は Bohr(1 Bohr $=$ 0.5292 Å)である．この単位系では光速は $c = 1/\alpha \simeq 137.036$ となることを15.8節で使う．

15.1 なぜ擬ポテンシャルか

15.1.1 内殻電子は化学的に不活性である

原子中の電子は，エネルギー的にも空間的にも，はっきり二つの階層に分かれている．シリコン原子($Z=14$，電子配置 $1s^2 2s^2 2p^6 3s^2 3p^2$)を例にとろう．コーン・シャム固有値で見ると，1s 準位はおよそ $-65$ Ha($\approx -1.8$ keV)，2s と 2p はおよそ $-5$ Ha と $-3.5$ Ha($-100$ eV 前後)という深さにあるのに対し，価電子である 3s と 3p は $-0.4$〜$-0.1$ Ha($-10$〜$-3$ eV)程度の浅さにある．内殻と価電子の間には2桁を超えるエネルギーギャップがある．

空間的にも同じことが言える．1s 軌道は核から $1/Z \approx 0.07$ Bohr 程度の距離に局在し，2s・2p も 0.3 Bohr 程度までにほぼ収まる．一方，化学結合を担う 3s・3p は 1〜3 Bohr のスケールに広がる．結晶中で隣の原子と重なり合うのは価電子の裾だけであり，内殻軌道同士の重なりは指数関数的に小さい．したがって固体や分子を作っても内殻電子の状態は孤立原子のときからほとんど変化しない．これを積極的に近似として採用し，内殻軌道を原子のまま凍結するのが凍結内殻近似(frozen-core approximation)である．多くの元素で，凍結内殻近似が結合エネルギーや格子定数に与える誤差は化学精度(1 kcal/mol $\approx$ 1.6 mHa)以下であることが全電子計算との比較で確かめられている．

図15.1 シリコン原子の内殻・価電子のエネルギースケール(左)と動径波動関数 $u(r)=rR(r)$ の空間スケール(右)の模式図．内殻 1s は $1/Z$ 程度に鋭く局在し，価電子 3s は広がっているが，内殻との直交性のため核近傍に2つの節を持つ．

15.1.2 それでも計算は軽くならない:平面波収束の悪夢

内殻を凍結すれば対角化すべき軌道の数は減る．しかしそれだけでは計算は軽くならない．問題は二つ残る．

核のポテンシャル $-Z/r$ が残ること．ポテンシャルが $r\to0$ で発散するため，そこに束縛される波動関数は核近傍で急峻に変化する．
価電子の波動関数が核近傍で節を持つこと．3s 軌道は 1s・2s と直交しなければならないため($\braket{\psi_{1s}|\psi_{3s}}=0$ など)，内殻が局在する 0.3 Bohr 程度の領域内で符号を2回変える(図15.1右)．滑らかな関数同士の内積は消えないから，直交性は激しい振動によってしか実現できないのである．

この「細かい構造」が平面波展開にとってどれほど致命的か，具体的に見積もってみよう．まず，体積 $\Omega$ の単位胞で運動エネルギーカットオフ $E_{\rm cut}$($|\kk+\bm{G}|^2/2 \le E_{\rm cut}$)を課したときの平面波の本数を数えておく．逆格子点は $\bm{G}$ 空間で1点あたり $(2\pi)^3/\Omega$ の体積を占める(第12章)．したがって半径 $q_{\max}=\sqrt{2E_{\rm cut}}$ の球に入る平面波数は，球の体積を1点あたりの体積で割って

$$ \begin{equation} N_{\rm PW} \simeq \frac{\dfrac{4\pi}{3}q_{\max}^3}{\dfrac{(2\pi)^3}{\Omega}} = \frac{\Omega\, q_{\max}^3}{6\pi^2} \label{eq:15-npw} \end{equation} $$

となる．1行目から2行目へは分母・分子を整理しただけである($4\pi/3 \div 8\pi^3 = 1/(6\pi^2)$)．重要なのは $N_{\rm PW} \propto q_{\max}^3 \propto E_{\rm cut}^{3/2}$ というスケーリングである．

例:シリコン 1s 軌道に必要なカットオフ

水素様原子の 1s 軌道は $R_{1s}(r) \propto e^{-Zr}$ で，特徴的スケールは $1/Z$ である．シリコンでは $1/Z = 1/14 \approx 0.071$ Bohr．この関数のフーリエ変換は $\tilde{\phi}_{1s}(q) \propto (q^2+Z^2)^{-2}$ という裾の重い形をしており(章末演習15.1で導出する)，ノルムの 99.9% を回収するには $q_{\max} \approx 5Z$ 程度まで波数を含める必要がある．シリコンなら $$q_{\max} \approx 5 \times 14 = 70~\text{Bohr}^{-1},\qquad E_{\rm cut} = \frac{q_{\max}^2}{2} \approx 2450~\text{Ha} \approx 67{,}000~\text{eV}$$ である．一方，後述する擬ポテンシャルを使ったシリコンの計算では $E_{\rm cut} \approx 15$ Ha(約 400 eV)で十分収束する．シリコンのダイヤモンド構造の基本単位胞($\Omega \approx 270$ Bohr³，2原子)に式 \eqref{eq:15-npw} を適用すると， $$N_{\rm PW}(\text{全電子}) \approx \frac{270 \times 70^3}{6\pi^2} \approx 1.6\times10^{6}, \qquad N_{\rm PW}(\text{擬ポテンシャル}) \approx \frac{270 \times (\sqrt{30})^3}{6\pi^2} \approx 7.5\times10^{2}.$$ 基底の数で約2000倍，対角化のコストは $N^3$ に比例するから約 $10^{10}$ 倍の差になる．1s をあからさまに平面波で展開する計算が現実的でないことがわかる．しかも価電子 3s も核近傍の節のせいで同程度の細かさを要求するから，内殻を凍結するだけでは不十分で，価電子の波動関数そのものを核近傍で「滑らかに作り替える」必要がある．

15.1.3 戦略:芯を消して価電子だけを扱う

そこで次の戦略をとる．

定義:擬ポテンシャル法の基本戦略

(1) 内殻電子は原子のまま凍結し，計算の自由度から消去する．
(2) 「核 $+$ 内殻電子」の系が価電子に及ぼす影響を，あらかじめ原子計算から構成した有効ポテンシャル(擬ポテンシャル)で置き換える．
(3) その際，価電子の波動関数を，核近傍(カットオフ半径 $r_c$ の内側)では節のない滑らかな擬波動関数に置き換える．$r_c$ の外側では真の波動関数と厳密に一致させる．

こうすれば，解くべきコーン・シャム方程式には浅く滑らかなポテンシャルと滑らかな波動関数しか現れず，平面波でも少数の基底で収束する．もちろん，この置き換えが物理を壊さないことを保証する条件が必要である．それを与えるのが本章の主題であり，答えを先に言えば「$r_c$ の外側で真の原子と同じ散乱をすること」(15.3節)，そしてそれを1つのエネルギー点の周りで有限のエネルギー幅にわたって保証する処方が「ノルム保存」(15.4節)である．

なお，内殻を消去せずに扱う全電子法(LAPW法など．第14章で数値解法の分類として整理した)も広く使われている．全電子法は内殻由来の物性(内殻準位シフト，メスバウアー分光など)を直接扱える一方，計算コストは高い．擬ポテンシャル法の精度は最終的に全電子法との比較で検証される(15.9節)．両者は競合というより，精度の基準系とその軽量な代理という補完関係にある．

15.2 OPWからPhillips–Kleinmanへ

擬ポテンシャルという概念は，天下りに導入されたものではなく，平面波基底の収束改善という実務的な工夫(OPW法)から論理必然的に「発見」されたものである．この節ではその道筋をたどり，擬ポテンシャルの原型であるPhillips–Kleinman(PK)変換を完全に導出する．ここでの目標は，(i) 芯状態への射影が反発的なポテンシャルを生むこと，(ii) それが核の引力とほぼ相殺すること(キャンセレーション定理)，(iii) しかしPK形式のままでは実用にならない3つの欠点があること，の3点を式の上で確認することである．

第12章との関係(記号の対応)

OPW 法と PK 変換は，第12章12.8節でも「単純金属でなぜ自由電子模型が働くのか」という問いへの答えとして導出済みである．本節は擬ポテンシャル法の出発点として同じ変換を，本章の記号でもう一度たどる(本章だけを読んでも閉じるようにしてある)．導出の骨格と結論——反発ポテンシャル $\hat V_R$，キャンセレーション定理，PK 形式の3つの欠点——は両章で完全に一致しており，記号だけが次のように対応する．

芯状態:第12章 $\left|\phi_c\right\rangle$，固有値 $\varepsilon_c$ ⇔ 本章 $\left|\psi_c\right\rangle$，固有値 $E_c$
真の価電子状態:第12章 $\left|\phi_v\right\rangle$，固有値 $\varepsilon_v$ ⇔ 本章 $\left|\Psi\right\rangle$，固有値 $E$
滑らかな擬波動関数:第12章 $\left|\tilde\phi\right\rangle$ ⇔ 本章 $\left|\phi\right\rangle$

したがって $\hat V_R=\sum_c(\varepsilon_v-\varepsilon_c)\left|\phi_c\right\rangle\left\langle\phi_c\right|$(第12章)と \eqref{eq:15-vr} は同一の演算子である．第12章を読んだ読者は 15.2.3 項(3つの欠点)から読み始めてもよい．

15.2.1 直交化平面波(OPW)法

第12章で学んだように，周期系の波動関数はBloch関数であり，平面波 $|{\rm PW},\kk\rangle \propto e^{i\kk\cdot\rr}$ の重ね合わせで書ける．前節で見たとおり，素朴な平面波展開は内殻との直交性が強制する振動のために収束しない．Herringは1940年，基底の側にあらかじめ直交性を組み込むという解決を提案した(OPW法)．ハミルトニアン $\hat{H}$ の固有状態を内殻($c$)と価電子($v$)に分け，内殻軌道 $|\psi_c\rangle$($\hat{H}|\psi_c\rangle = E_c|\psi_c\rangle$)は原子計算から既知だとする．平面波から内殻成分を射影して除いた

$$ \begin{equation} |{\rm OPW},\kk\rangle = |{\rm PW},\kk\rangle - \sum_c |\psi_c\rangle\langle\psi_c|{\rm PW},\kk\rangle \label{eq:15-opw-def} \end{equation} $$

を新しい基底(直交化平面波)とする．これが実際にすべての内殻軌道と直交することは，次のように1行ずつ確かめられる．任意の内殻軌道 $|\psi_{c'}\rangle$ との内積をとると

$$ \begin{align} \langle\psi_{c'}|{\rm OPW},\kk\rangle &= \langle\psi_{c'}|{\rm PW},\kk\rangle - \sum_c \langle\psi_{c'}|\psi_c\rangle\langle\psi_c|{\rm PW},\kk\rangle \label{eq:15-opw-orth1}\\ &= \langle\psi_{c'}|{\rm PW},\kk\rangle - \sum_c \delta_{c'c}\,\langle\psi_c|{\rm PW},\kk\rangle \label{eq:15-opw-orth2}\\ &= \langle\psi_{c'}|{\rm PW},\kk\rangle - \langle\psi_{c'}|{\rm PW},\kk\rangle = 0. \label{eq:15-opw-orth3} \end{align} $$

1行目は定義 \eqref{eq:15-opw-def} を代入しただけである．2行目では内殻軌道同士の正規直交性 $\langle\psi_{c'}|\psi_c\rangle=\delta_{c'c}$ を使った(同一ハミルトニアンの束縛固有状態だから直交にとれる)．3行目ではクロネッカーデルタで和を潰した．OPWは「遠方では平面波，核近傍では内殻に直交するための振動を最初から持つ関数」であり，これを基底に使うと少数の基底で価電子状態を表現できる．実際，1950年代までの半導体のバンド計算はOPW法で行われた．

15.2.2 Phillips–Kleinman変換:反発ポテンシャルの出現

Phillips と Kleinman は1959年，OPW法の構造を演算子の言葉で読み替えると「擬ポテンシャル」という概念が現れることを見抜いた．これから，その変換を導出する．出発点は，価電子の真の波動関数 $|\Psi\rangle$ を「滑らかな成分 $|\phi\rangle$」と「内殻に直交させるための補正」に分けて書くことである．滑らかな成分を平面波の重ね合わせ $|\phi\rangle = \sum_{\bm G} C_{\bm G}|{\rm PW},\bm{G}\rangle$ とし，OPWの構造 \eqref{eq:15-opw-def} をそのまま線形結合に持ち上げて

$$ \begin{equation} |\Psi\rangle = |\phi\rangle - \sum_c |\psi_c\rangle\langle\psi_c|\phi\rangle \label{eq:15-pk-ansatz} \end{equation} $$

とおく．$|\Psi\rangle$ は構成により内殻と直交している(前小節と同じ計算)．これをシュレーディンガー方程式 $\hat{H}|\Psi\rangle = E|\Psi\rangle$ に代入し，滑らかな成分 $|\phi\rangle$ が満たす方程式を導く．まず左辺を計算する．

$$ \begin{align} \hat{H}|\Psi\rangle &= \hat{H}|\phi\rangle - \sum_c \hat{H}|\psi_c\rangle\langle\psi_c|\phi\rangle \label{eq:15-pk-lhs1}\\ &= \hat{H}|\phi\rangle - \sum_c E_c|\psi_c\rangle\langle\psi_c|\phi\rangle. \label{eq:15-pk-lhs2} \end{align} $$

1行目はansatz \eqref{eq:15-pk-ansatz} を代入して $\hat{H}$ の線形性を使った．$\langle\psi_c|\phi\rangle$ はただの数なので $\hat{H}$ の外に出せることに注意．2行目では内殻軌道が固有状態であること $\hat{H}|\psi_c\rangle = E_c|\psi_c\rangle$ を使った．次に右辺は

$$ \begin{equation} E|\Psi\rangle = E|\phi\rangle - \sum_c E|\psi_c\rangle\langle\psi_c|\phi\rangle \label{eq:15-pk-rhs} \end{equation} $$

である(ansatzに数 $E$ を掛けただけ)．式 \eqref{eq:15-pk-lhs2} と \eqref{eq:15-pk-rhs} を等しいと置き，両辺から共通項を整理する．

$$ \begin{align} \hat{H}|\phi\rangle - \sum_c E_c|\psi_c\rangle\langle\psi_c|\phi\rangle &= E|\phi\rangle - \sum_c E|\psi_c\rangle\langle\psi_c|\phi\rangle \label{eq:15-pk-eq1}\\ \hat{H}|\phi\rangle + \sum_c (E - E_c)|\psi_c\rangle\langle\psi_c|\phi\rangle &= E|\phi\rangle. \label{eq:15-pk-eq2} \end{align} $$

2行目へは，右辺の $-\sum_c E|\psi_c\rangle\langle\psi_c|\phi\rangle$ を左辺に移項し，左辺の $-\sum_c E_c|\psi_c\rangle\langle\psi_c|\phi\rangle$ と同類項としてまとめた($-E_c + E = E - E_c$)．ここで

$$ \begin{equation} \hat{V}_R \equiv \sum_c (E - E_c)\,|\psi_c\rangle\langle\psi_c| \label{eq:15-vr} \end{equation} $$

と定義すれば，式 \eqref{eq:15-pk-eq2} は次のようにまとまる．これがPhillips–Kleinman方程式である．

$$ \begin{equation} \left(\hat{H} + \hat{V}_R\right)|\phi\rangle = \left(\hat{T} + \hat{v}_{\rm ext} + \hat{V}_R\right)|\phi\rangle = E|\phi\rangle, \qquad \hat{V}_R = \sum_c (E-E_c)|\psi_c\rangle\langle\psi_c| \label{eq:15-pk} \end{equation} $$

何が得られたか．滑らかな成分 $|\phi\rangle$ は，それ自身が固有値方程式を満たし，しかもその固有値は真の価電子エネルギー $E$ と同一である．ただし $|\phi\rangle$ が感じるポテンシャルは，真の外場 $\hat{v}_{\rm ext}$ に射影演算子の項 $\hat{V}_R$ を加えた

$$ \begin{equation} \hat{v}_{\rm PK} = \hat{v}_{\rm ext} + \sum_c (E-E_c)|\psi_c\rangle\langle\psi_c| \label{eq:15-pk-veff} \end{equation} $$

である．価電子のエネルギー $E$ は内殻のエネルギー $E_c$ よりずっと高いから $E - E_c > 0$ であり，$\hat{V}_R$ の期待値 $\langle\phi|\hat{V}_R|\phi\rangle = \sum_c (E-E_c)|\langle\psi_c|\phi\rangle|^2 \ge 0$ は必ず非負，すなわち $\hat{V}_R$ は反発的である．しかも $|\psi_c\rangle$ が核近傍に局在しているから，$\hat{V}_R$ が効くのも核近傍だけである．つまり「核近傍でだけ働く強い反発」が，深いクーロン引力の上に自動的に乗る．

物理的意味:キャンセレーション定理

核近傍で価電子が感じる正味のポテンシャル $\hat{v}_{\rm PK}$ は，深い引力 $\hat{v}_{\rm ext} \sim -Z/r$ と，内殻直交性に由来する強い反発 $\hat{V}_R$ の和であり，両者は大きく打ち消し合って浅く弱い正味ポテンシャルだけが残る(図15.2)．反発の起源はパウリ原理である:価電子は内殻と直交しなければならず，内殻が占める領域に「入り込めない」．この排除効果がエネルギー的には反発ポテンシャルとして現れる．この相殺(キャンセレーション定理)は，「なぜナトリウムの3s電子がほとんど自由電子のように振る舞うのか」「なぜ単純金属で自由電子モデル(第3章)が驚くほどうまくいくのか」という古い疑問への答えでもある．核の裸の電荷 $+11$ は内殻電子と直交性反発によってほぼ完全に遮蔽され，価電子が見るのは弱い残差ポテンシャルだけなのである．

図15.2 キャンセレーション定理の模式図．核のクーロン引力(赤)と内殻直交性が生む反発 $\hat{V}_R$(緑)は核近傍で大きく打ち消し合い，価電子が感じる正味のポテンシャル(青)は浅く滑らかになる．

15.2.3 PK形式の3つの欠点

PK方程式 \eqref{eq:15-pk} は擬ポテンシャルの概念的出発点だが，そのまま実用に使うには3つの深刻な欠点がある．順に式で確認する．

(1) エネルギー依存性．$\hat{V}_R = \sum_c(E-E_c)|\psi_c\rangle\langle\psi_c|$ には，これから求めるべき固有値 $E$ 自身が入っている．つまり状態ごとにポテンシャルが違い，方程式は非線形になる．また異なる $E$ に属する解 $|\phi\rangle$ 同士は(異なる演算子の固有状態なので)直交しない．

(2) 非局所性．$\hat{V}_R$ は $|\psi_c\rangle\langle\psi_c|$ という射影演算子であり，位置表示では $\langle\rr|\hat{V}_R|\phi\rangle = \sum_c (E-E_c)\psi_c(\rr)\int \dd^3 r'\, \psi_c^*(\rr')\phi(\rr')$ という積分核として働く．一点 $\rr$ での値が波動関数の全空間の値に依存する．これ自体は本質的な障害ではなく，後で見るように現代の擬ポテンシャルも非局所性を積極的に利用する(15.6節)．ただし扱いには工夫が要る．

(3) ノルムの非保存．これが最も重要な欠点である．射影演算子 $\hat{P} \equiv \sum_c|\psi_c\rangle\langle\psi_c|$ を使うとansatz \eqref{eq:15-pk-ansatz} は $|\Psi\rangle = (1-\hat{P})|\phi\rangle$ と書ける．$\hat{P}$ は射影演算子だから $\hat{P}^2 = \hat{P}$，$\hat{P}^\dagger = \hat{P}$ である(内殻軌道の正規直交性 $\langle\psi_c|\psi_{c'}\rangle=\delta_{cc'}$ から $\hat{P}^2 = \sum_{cc'}|\psi_c\rangle\delta_{cc'}\langle\psi_{c'}| = \hat{P}$ と直接確かめられる)．すると $|\Psi\rangle$ のノルムは

$$ \begin{align} \langle\Psi|\Psi\rangle &= \langle\phi|(1-\hat{P})^\dagger(1-\hat{P})|\phi\rangle \label{eq:15-norm1}\\ &= \langle\phi|(1-\hat{P})(1-\hat{P})|\phi\rangle \label{eq:15-norm2}\\ &= \langle\phi|1 - 2\hat{P} + \hat{P}^2|\phi\rangle \label{eq:15-norm3}\\ &= \langle\phi|1 - \hat{P}|\phi\rangle = \langle\phi|\phi\rangle - \sum_c \left|\langle\psi_c|\phi\rangle\right|^2. \label{eq:15-norm4} \end{align} $$

1行目はノルムの定義に $|\Psi\rangle=(1-\hat{P})|\phi\rangle$ を代入した．2行目はエルミート性 $(1-\hat{P})^\dagger = 1-\hat{P}$，3行目は展開，4行目は冪等性 $\hat{P}^2=\hat{P}$ により $-2\hat{P}+\hat{P}=-\hat{P}$ とまとめ，最後に $\hat{P}$ の定義を代入した．真の波動関数を $\langle\Psi|\Psi\rangle=1$ と規格化すると，滑らかな成分は $\langle\phi|\phi\rangle = 1 + \sum_c|\langle\psi_c|\phi\rangle|^2 > 1$ となる．擬波動関数 $\phi$ は核近傍に余分な電荷を持ち，$|\phi|^2$ は正しい電子密度を与えない．密度が主役のDFTにとってこれは致命的である:密度が違えばHartreeポテンシャルも交換相関ポテンシャルも狂い，自己無撞着計算全体が汚染される．

まとめると，PK変換は「芯を消して浅いポテンシャルに置き換えられる」ことを原理的に示したが，(1) エネルギー依存，(2) 非局所，(3) ノルム非保存という3つの問題を残した．現代のノルム保存擬ポテンシャルは，(3)を設計原理として明示的に解決し，その副産物として(1)も実用上十分な精度で解決する(15.4節)．(2)は解決すべき欠点というより，計算コストを下げる道具として洗練される(15.6節)．次節ではまず，「良い擬ポテンシャル」の合格基準そのものを散乱理論の言葉で定式化する．

15.3 散乱理論の見方:対数微分

擬ポテンシャルは真のポテンシャルの「偽物」である．では，偽物が本物の代わりを務めるための合格基準は何か．この節では，その基準が散乱理論の言葉で「カットオフ半径の外側で同じ散乱を引き起こすこと」，より具体的には「波動関数の対数微分が一致すること」として定式化できることを示す．これが次節のノルム保存条件の土台になる．

15.3.1 発想:結晶の中の原子は散乱体である

結晶や分子の中の1個の原子に注目しよう．その原子のポテンシャルが強く働くのは半径 $r_c$ 程度の球の内側だけであり，球の外側では価電子は他の原子や電子の作る環境の中を伝播している．外側の環境から見ると，この原子球は「入ってきた電子波を散乱して送り返す装置」である．もし偽物の原子球(擬ポテンシャル)が，興味あるエネルギー範囲のすべての電子波に対して本物と同じ散乱波を送り返すなら，外側の世界は両者を区別できない．したがって結合エネルギーもバンド構造も正しく再現される．これが擬ポテンシャルの合格基準である．

15.3.2 球対称ポテンシャルの散乱と位相のずれ

基準を数式にするため，球対称ポテンシャル $V(r)$($r > r_0$ で十分速くゼロになるとする)による散乱を考える．エネルギー $\varepsilon = q^2/2$ の電子の波動関数は，球対称性により角運動量 $l$ ごとに分離できる(部分波展開)．第14章14.6節と同様に，動径波動関数を $R_l(r)$，$u_l(r) \equiv rR_l(r)$ とおくと，$u_l$ は動径シュレーディンガー方程式

$$ \begin{equation} \left[-\frac{1}{2}\frac{\dd^2}{\dd r^2} + \frac{l(l+1)}{2r^2} + V(r)\right]u_l(r) = \varepsilon\, u_l(r) \label{eq:15-radial} \end{equation} $$

を満たす($u_l(0)=0$，原点近傍で $u_l \propto r^{l+1}$)．遠心力項 $l(l+1)/2r^2$ は3次元ラプラシアンの角度部分から来るもので，導出は第14章で行った．

数学ノート:球ベッセル関数と位相のずれ

ポテンシャルがない領域($V=0$)では，式 \eqref{eq:15-radial} は $R_l$ について $$R_l'' + \frac{2}{r}R_l' + \left[q^2 - \frac{l(l+1)}{r^2}\right]R_l = 0$$ となる．これは変数 $x = qr$ の球ベッセル方程式であり，2つの独立解は球ベッセル関数 $j_l(x)$(原点で正則)と球ノイマン関数 $n_l(x)$(原点で発散)である．最低次は $$j_0(x) = \frac{\sin x}{x},\qquad n_0(x) = -\frac{\cos x}{x}$$ で，一般の $l$ でも遠方($x \gg l$)では $$j_l(x) \to \frac{\sin(x - l\pi/2)}{x},\qquad n_l(x) \to -\frac{\cos(x - l\pi/2)}{x}$$ と振る舞う(この漸近形は $j_0, n_0$ の表式と，隣接する $l$ を結ぶ漸化式から帰納的に得られる)．ポテンシャルの外側の解はこの2つの線形結合で書くしかない．散乱理論の中心的事実は，エネルギー $\varepsilon$ における散乱の観測可能な情報(微分断面積など)が，各部分波の混合比，すなわち下で定義する位相のずれ $\eta_l(\varepsilon)$ にすべて集約されることである．

ポテンシャルの外側 $r > r_0$ では，一般解を係数 $A_l$ と混合角 $\eta_l$ を使って

$$ \begin{equation} R_l(r) = A_l\left[\cos\eta_l\, j_l(qr) - \sin\eta_l\, n_l(qr)\right] \qquad (r > r_0) \label{eq:15-outside} \end{equation} $$

と書く(2つの独立解の任意の線形結合はこの形に規格化できる)．この $\eta_l$ が「位相のずれ」と呼ばれる理由は，遠方の漸近形を見ればわかる．数学ノートの漸近形を \eqref{eq:15-outside} に代入すると

$$ \begin{align} R_l(r) &\to \frac{A_l}{qr}\left[\cos\eta_l \sin\!\left(qr - \frac{l\pi}{2}\right) + \sin\eta_l \cos\!\left(qr - \frac{l\pi}{2}\right)\right] \label{eq:15-asymp1}\\ &= \frac{A_l}{qr}\,\sin\!\left(qr - \frac{l\pi}{2} + \eta_l\right). \label{eq:15-asymp2} \end{align} $$

1行目は漸近形の代入($n_l$ の符号に注意)，2行目は正弦の加法定理 $\sin\alpha\cos\beta + \cos\alpha\sin\beta = \sin(\alpha+\beta)$ を逆向きに使った．つまり，ポテンシャルがあることの唯一の痕跡は，遠方の正弦波の位相が自由な場合($\eta_l = 0$)から $\eta_l$ だけずれることである(図15.3)．引力ポテンシャルは波を内側に「引き込む」ので $\eta_l > 0$ になる．

図15.3 位相のずれの意味．ポテンシャル領域($r < r_0$，薄青)の外側では，波動関数は自由な正弦波と同じ形をもち，唯一の違いは位相が $\eta_l$ だけずれていることである．引力は波を引き込み($\eta_l > 0$)，外の世界にはこの位相ずれだけが伝わる．

15.3.3 位相のずれは対数微分で決まる

位相のずれ $\eta_l$ は，ポテンシャル内部の情報をどのように「受け取る」のだろうか．答えは，境界 $r_0$ における波動関数の対数微分

$$ \begin{equation} D_l(\varepsilon, r) \equiv \frac{u_l'(r)}{u_l(r)} = \frac{\dd}{\dd r}\ln u_l(r) \label{eq:15-Ddef} \end{equation} $$

を通してである($'$ は $r$ 微分)．対数微分は波動関数の規格化によらない量であることに注意しよう($u_l \to \lambda u_l$ としても分子分母で $\lambda$ が消える)．文献では無次元化した $r\,\dd\ln u_l/\dd r$ や，$R_l$ に基づく $R_l'/R_l$ もよく使われるが，$u_l = rR_l$ の対数をとって微分すれば $\dd \ln u_l/\dd r = 1/r + \dd \ln R_l/\dd r$ となるから，両者は既知の項 $1/r$ の差しかなく，どちらを使っても以下の議論は変わらない．

導出:位相のずれと対数微分の関係

内側($r \le r_0$)でシュレーディンガー方程式を解いて得た解の，境界での対数微分を $\gamma \equiv R_l'(r_0)/R_l(r_0)$ とする．外側の解 \eqref{eq:15-outside} と，値および1階微分が $r_0$ で接続する条件(シュレーディンガー方程式が2階だから，これで解は一意に繋がる)を書く．値の比を取る形で，対数微分の連続性として書くと

$$ \begin{equation} \gamma = \frac{R_l'(r_0)}{R_l(r_0)} = \frac{q\left[\cos\eta_l\, j_l'(qr_0) - \sin\eta_l\, n_l'(qr_0)\right]}{\cos\eta_l\, j_l(qr_0) - \sin\eta_l\, n_l(qr_0)}. \label{eq:15-match1} \end{equation} $$

ここで $j_l', n_l'$ は引数 $x=qr$ についての微分であり，$\dd j_l(qr)/\dd r = q\,j_l'(qr)$ という連鎖律で $q$ が前に出た．分母を払う(両辺に分母を掛ける):

$$ \begin{equation} \gamma\cos\eta_l\, j_l(qr_0) - \gamma\sin\eta_l\, n_l(qr_0) = q\cos\eta_l\, j_l'(qr_0) - q\sin\eta_l\, n_l'(qr_0). \label{eq:15-match2} \end{equation} $$

$\cos\eta_l$ を含む項を左辺に，$\sin\eta_l$ を含む項を右辺に集める:

$$ \begin{equation} \cos\eta_l\left[\gamma\, j_l(qr_0) - q\, j_l'(qr_0)\right] = \sin\eta_l\left[\gamma\, n_l(qr_0) - q\, n_l'(qr_0)\right]. \label{eq:15-match3} \end{equation} $$

両辺を $\cos\eta_l\left[\gamma n_l - q n_l'\right]$ で割れば

$$ \begin{equation} \tan\eta_l(\varepsilon) = \frac{q\, j_l'(qr_0) - \gamma\, j_l(qr_0)}{q\, n_l'(qr_0) - \gamma\, n_l(qr_0)} \label{eq:15-tan-eta} \end{equation} $$

を得る(最後に分子分母の符号を同時に反転して見やすく並べ替えた)．

何が得られたか．式 \eqref{eq:15-tan-eta} の右辺に現れる内部の情報は対数微分 $\gamma = D_l(\varepsilon, r_0)$ ただ一つである($j_l, n_l$ は普遍的な既知関数)．つまり:

物理的意味:対数微分がすべてを語る

半径 $r_0$ の球の外の世界にとって，球の中身の情報は境界での対数微分 $D_l(\varepsilon, r_0)$ に完全に集約される．中のポテンシャルがどんな形をしていても，すべての $l$ と，問題となるエネルギー範囲のすべての $\varepsilon$ で $D_l(\varepsilon, r_0)$ が本物と一致していれば，外の世界(=他の原子，化学環境)は本物と偽物を区別できない．

定義:移植性(transferability)

擬ポテンシャルは1つの参照原子配置(通常は基底状態の中性原子)で構成される．それを異なる化学環境(分子，固体，表面，イオン結晶…)に持ち込んでも全電子計算と同じ結果を与える性質を移植性という．上の議論から，移植性の定量的な中身は「価電子が関与するエネルギー窓(占有状態から数eV上の非占有状態まで)にわたって，各 $l$ の対数微分 $D_l(\varepsilon, r_c)$ が全電子原子のものと一致すること」である．環境が変わることは，原子球に入射する波のエネルギーと混合比が変わることに相当するからである．

ここで難題が現れる．擬ポテンシャルの構成(次節以降)では，参照エネルギー $\varepsilon_l$(通常は原子の価電子固有値)において $r \ge r_c$ で擬波動関数を全電子波動関数に一致させる．すると $D_l(\varepsilon_l, r_c)$ はその1点では自動的に一致する．しかし移植性が要求するのは有限のエネルギー幅での一致である．1点での一致から幅のある一致へ—この飛躍を可能にするのが，次節のノルム保存条件である．

15.4 ノルム保存条件

この節が本章の核心である．これから証明するのは，次の一見不思議な恒等式である:対数微分のエネルギー微分は，その半径より内側の波動関数のノルム(確率の総量)で決まる．すなわち

$$-\frac{\partial D_l(\varepsilon, r_c)}{\partial \varepsilon} \propto \int_0^{r_c} |u_l(r)|^2\, \dd r.$$

この恒等式が成り立つなら，話は一気に片付く．擬波動関数が $r_c$ の内側で全電子波動関数と同じノルムを持つように作れば(ノルム保存)，$D_l$ の値だけでなくエネルギー微分 $\partial D_l/\partial\varepsilon$ まで参照点で一致し，対数微分は参照エネルギーの近傍で1次のオーダーまで一致する．つまり1点での工作が有限幅での散乱の一致を保証する．この処方を設計原理として確立したのが Hamann, Schlüter, Chiang (1979) のノルム保存擬ポテンシャルである．

15.4.1 準備:動径方程式とWronskian

数学ノート:動径方程式とWronskian

2階線形微分方程式 $u'' = f(r)\,u$(1階微分の項がない形)の2つの解 $u_1, u_2$ に対し， $$W[u_1, u_2](r) \equiv u_1(r)u_2'(r) - u_1'(r)u_2(r)$$ をWronskian(ロンスキアン)という．その微分は $$W' = (u_1u_2')' - (u_1'u_2)' = u_1'u_2' + u_1u_2'' - u_1''u_2 - u_1'u_2' = u_1u_2'' - u_1''u_2$$ となる(積の微分則を各項に適用すると $u_1'u_2'$ が相殺する)．同じ方程式の2解なら $u_1u_2'' - u_1''u_2 = u_1(fu_2) - (fu_1)u_2 = 0$ で $W$ は定数になる(これが「2解の独立性の判定」にWronskianが使われる理由である)．ここで使うのは少しひねった状況で，$f$ がわずかに異なる2つの方程式($f_1 \ne f_2$)の解を組ませる．すると $W' = u_1u_2'' - u_1''u_2 = (f_2 - f_1)u_1u_2$ となり，$W$ は定数でなくなる．その「ずれ」の積分がちょうど波動関数のノルムを拾い上げる—これが以下の導出のからくりである．

動径方程式 \eqref{eq:15-radial} は，両辺に $-2$ を掛けて移項すると $$u_l'' = \left[\frac{l(l+1)}{r^2} + 2V(r) - 2\varepsilon\right]u_l$$ という $u'' = f u$ の形になる($f$ にエネルギーが入っていることに注目)．以下ではこの形を使う．

15.4.2 中心的恒等式の導出

これから，同じポテンシャル・同じ $l$ でエネルギーだけが異なる2つの動径解を組み合わせ，Wronskianの微分を積分することで恒等式を導く．長い導出だが，一歩ずつ進めば各行は単純である．

導出:対数微分のエネルギー微分とノルムの関係

ステップ1:2つのエネルギーで動径方程式を書く．エネルギー $\varepsilon_1$，$\varepsilon_2$ に対する動径方程式の解を $u_1(r) \equiv u_l(\varepsilon_1, r)$，$u_2(r) \equiv u_l(\varepsilon_2, r)$ とする．どちらも原点で正則な解($u \propto r^{l+1}$，$r\to0$)を選ぶ．固有値である必要はなく，各エネルギーで原点正則解は(規格化を除き)一意に定まる．ポテンシャルと $l$ は共通である．数学ノートの形で書くと

$$ \begin{align} u_1'' &= \left[\frac{l(l+1)}{r^2} + 2V(r) - 2\varepsilon_1\right]u_1, \label{eq:15-u1}\\ u_2'' &= \left[\frac{l(l+1)}{r^2} + 2V(r) - 2\varepsilon_2\right]u_2. \label{eq:15-u2} \end{align} $$

ステップ2:交差して引く．式 \eqref{eq:15-u2} に $u_1$ を掛け，式 \eqref{eq:15-u1} に $u_2$ を掛けて，前者から後者を引く:

$$ \begin{align} u_1 u_2'' - u_2 u_1'' &= u_1\left[\frac{l(l+1)}{r^2} + 2V - 2\varepsilon_2\right]u_2 - u_2\left[\frac{l(l+1)}{r^2} + 2V - 2\varepsilon_1\right]u_1 \label{eq:15-cross1}\\ &= \left(-2\varepsilon_2 + 2\varepsilon_1\right)u_1 u_2 = 2(\varepsilon_1 - \varepsilon_2)\,u_1 u_2. \label{eq:15-cross2} \end{align} $$

2行目では，遠心力項とポテンシャル項が $u_1u_2$ の係数として完全に共通なので相殺し，エネルギー項だけが生き残った．ポテンシャルの詳細がここで消えるのがこの計算の要点である．

ステップ3:左辺をWronskianの微分と同定する．数学ノートで確認したとおり $$\frac{\dd}{\dd r}\left[u_1 u_2' - u_1' u_2\right] = u_1 u_2'' - u_1'' u_2$$ である(積の微分で $u_1'u_2'$ が相殺)．よって式 \eqref{eq:15-cross2} は

$$ \begin{equation} \frac{\dd}{\dd r} W[u_1,u_2](r) = 2(\varepsilon_1 - \varepsilon_2)\, u_1(r) u_2(r), \qquad W[u_1,u_2] = u_1u_2' - u_1'u_2 \label{eq:15-wronsk-der} \end{equation} $$

と書ける．

ステップ4:$0$ から $r_c$ まで積分する．微積分学の基本定理より

$$ \begin{equation} W[u_1,u_2](r_c) - W[u_1,u_2](0) = 2(\varepsilon_1 - \varepsilon_2)\int_0^{r_c} u_1(r)\, u_2(r)\, \dd r. \label{eq:15-wronsk-int} \end{equation} $$

ここで原点の境界項が消えることを確認する．原点近傍で $u_1, u_2 \propto r^{l+1}$，したがって $u_1', u_2' \propto r^l$ であり， $$W(r) = u_1u_2' - u_1'u_2 \propto r^{l+1}\cdot r^l = r^{2l+1} \xrightarrow{\ r\to0\ } 0$$ となる($l \ge 0$ だから $2l+1 \ge 1$)．よって $W(0) = 0$ であり，

$$ \begin{equation} u_1(r_c)u_2'(r_c) - u_1'(r_c)u_2(r_c) = 2(\varepsilon_1 - \varepsilon_2)\int_0^{r_c} u_1 u_2\, \dd r. \label{eq:15-wronsk-rc} \end{equation} $$

この式は後で(15.7節のMBK法の導出で)このままの形でもう一度使うので，覚えておいてほしい．

ステップ5:対数微分で書き直す．対数微分 $D_i \equiv u_i'(r_c)/u_i(r_c)$($i=1,2$)を使うと，左辺は $$u_1(r_c)u_2'(r_c) - u_1'(r_c)u_2(r_c) = u_1(r_c)u_2(r_c)\left[\frac{u_2'(r_c)}{u_2(r_c)} - \frac{u_1'(r_c)}{u_1(r_c)}\right] = u_1(r_c)u_2(r_c)\left[D_2 - D_1\right]$$ と因数分解できる(第1式の各項を $u_1(r_c)u_2(r_c)$ でくくっただけである)．よって

$$ \begin{equation} u_1(r_c)\,u_2(r_c)\left[D_l(\varepsilon_2, r_c) - D_l(\varepsilon_1, r_c)\right] = 2(\varepsilon_1 - \varepsilon_2)\int_0^{r_c} u_1 u_2\, \dd r. \label{eq:15-wronsk-D} \end{equation} $$

ステップ6:エネルギー差を無限小にする．$\varepsilon_1 = \varepsilon$，$\varepsilon_2 = \varepsilon + \delta\varepsilon$ とおき，$\delta\varepsilon \to 0$ の極限をとる．解はエネルギーに滑らかに依存するから $$u_2(r) = u_1(r) + O(\delta\varepsilon), \qquad D_l(\varepsilon_2, r_c) - D_l(\varepsilon_1, r_c) = \frac{\partial D_l(\varepsilon, r_c)}{\partial\varepsilon}\,\delta\varepsilon + O(\delta\varepsilon^2)$$ と展開できる(1つ目はエネルギーに関する連続性，2つ目はTaylor展開の1次)．これらを式 \eqref{eq:15-wronsk-D} に代入すると $$ u_l(\varepsilon, r_c)^2\,\frac{\partial D_l}{\partial\varepsilon}\,\delta\varepsilon + O(\delta\varepsilon^2) = -2\,\delta\varepsilon\int_0^{r_c} u_l(\varepsilon, r)^2\, \dd r + O(\delta\varepsilon^2) $$ となる(右辺では $\varepsilon_1 - \varepsilon_2 = -\delta\varepsilon$ を使い，$u_1u_2 \to u_l^2$ と置き換えた;置き換えの誤差は $\delta\varepsilon$ の高次)．両辺を $\delta\varepsilon$ で割り，$\delta\varepsilon \to 0$ とすれば，求める恒等式が得られる．

結果を章の最重要式として掲げる．

$$ \begin{equation} \frac{\partial D_l(\varepsilon, r_c)}{\partial\varepsilon} = -\frac{2}{u_l(\varepsilon, r_c)^2}\int_0^{r_c} \left|u_l(\varepsilon, r)\right|^2 \dd r \label{eq:15-normid} \end{equation} $$

何が得られたか．対数微分のエネルギー依存性(左辺)が，半径 $r_c$ より内側の確率の総量(右辺)だけで決まるという恒等式である．いくつか確認しておく．(i) 右辺は必ず負である:エネルギーを上げると対数微分は必ず減る(束縛状態がエネルギー順に並ぶことの背後にもこの単調性がある)．(ii) 両辺とも波動関数の規格化によらない:$u_l \to \lambda u_l$ とすると右辺は分母・分子とも $\lambda^2$ 倍で不変，左辺はもともと不変である．(iii) $R_l = u_l/r$ で書き直せば $\partial_\varepsilon D_l = -2\left[r_c R_l(r_c)\right]^{-2}\int_0^{r_c}r^2|R_l|^2\dd r$ となり，分子は3次元の動径確率密度そのものである(演習15.2)．無次元対数微分 $rD_l$ を使う流儀では両辺に $r_c$ が掛かるだけで内容は同じである．

15.4.3 ノルム保存擬ポテンシャルの根拠

恒等式 \eqref{eq:15-normid} を，全電子原子と擬原子の比較に適用しよう．参照エネルギー $\varepsilon_l$(原子の価電子固有値)で，擬波動関数 $u_l^{\rm PS}$ が次を満たすように構成したとする．

定義:ノルム保存擬ポテンシャルの条件(Hamann–Schlüter–Chiang)

擬原子の価電子固有値が全電子原子の固有値と一致する:$\varepsilon_l^{\rm PS} = \varepsilon_l^{\rm AE}$．
擬波動関数は節を持たない(節がなければ滑らかにでき，また最低固有状態として得られる)．
$r \ge r_c$ で擬波動関数は全電子波動関数と一致する:$u_l^{\rm PS}(r) = u_l^{\rm AE}(r)$．
ノルム保存: $$\int_0^{r_c}\left|u_l^{\rm PS}(r)\right|^2\dd r = \int_0^{r_c}\left|u_l^{\rm AE}(r)\right|^2\dd r.$$

このとき対数微分について何が言えるか．条件3により，$r_c$ での値と微分が一致するから $$D_l^{\rm PS}(\varepsilon_l, r_c) = D_l^{\rm AE}(\varepsilon_l, r_c)$$ がまず成り立つ(参照点での一致)．さらに条件3と4を恒等式 \eqref{eq:15-normid} に代入すると，右辺の分母 $u_l(r_c)^2$ は条件3で，積分は条件4で一致するから

$$ \begin{equation} \left.\frac{\partial D_l^{\rm PS}}{\partial\varepsilon}\right|_{\varepsilon_l} = \left.\frac{\partial D_l^{\rm AE}}{\partial\varepsilon}\right|_{\varepsilon_l} \label{eq:15-dderiv-match} \end{equation} $$

も成り立つ．したがってTaylor展開すれば

$$ \begin{equation} D_l^{\rm PS}(\varepsilon, r_c) - D_l^{\rm AE}(\varepsilon, r_c) = O\!\left((\varepsilon - \varepsilon_l)^2\right). \label{eq:15-taylor} \end{equation} $$

対数微分—したがって位相のずれ，したがって散乱—は，参照エネルギーの近傍で1次のオーダーまで自動的に一致する．これがノルム保存擬ポテンシャルの理論的根拠であり，「原子1配置で作ったポテンシャルがなぜ分子でも固体でも通用するのか」という移植性の問いへの定量的な答えである．価電子のバンド幅は数eV($\sim 0.1$–$0.3$ Ha)程度だから，参照点を価電子準位に取れば，1次までの一致で実用上十分な精度が出ることが多い．

物理的意味:ノルム保存のもう一つのご利益 — 静電気学

ノルム保存には散乱とは独立のboonがある．$r_c$ 内の電荷量が正しいので，Gaussの定理により球の外に作る静電ポテンシャルが正しくなる．半径 $r_c$ の球内の全電荷を $Q$ とすると，球対称電荷が外に作るポテンシャルは $Q/r$ だけで決まり，分布の詳細によらない．つまりノルム保存擬原子は，散乱(量子力学)と静電場(電磁気学)の両方で本物になりすますことができる．自己無撞着計算ではHartreeポテンシャルが密度から作られるから，この性質は精度に直結する．

図15.4 ノルム保存擬波動関数の構成．$r_c$ の外側では全電子波動関数(青)と厳密に一致し，内側では節のない滑らかな関数(緑・破線)に置き換える．内側のノルム(確率の総量)を保存するように作ると，恒等式 \eqref{eq:15-normid} により対数微分のエネルギー微分まで一致する．

残る問題は「では条件1〜4を満たす擬波動関数と擬ポテンシャルを具体的にどう作るか」である．次節で，現在も広く使われている Troullier–Martins の構成法を見る．

15.5 Troullier–Martins構成法

前節で，擬ポテンシャルが満たすべき条件は出そろった．しかし条件を並べただけでは物は作れない．実際に擬波動関数 $u_l^{\rm PS}(r)$ と擬ポテンシャル $v_l^{\rm PS}(r)$ を数値的に生成する処方が要る．ここで採用するのは逆問題の発想である．ふつうは「ポテンシャルを与えて波動関数を解く」が，擬ポテンシャルの構成では順序を逆にし，まず望ましい擬波動関数を関数形として書き下し，それを厳密に固有関数として持つポテンシャルを動径方程式から逆算する．この方針はKerker (1980) が提案し，TroullierとMartins (1991) が「平面波展開にとって最も柔らかい(必要カットオフの小さい)」形に磨き上げた．現在も標準的に使われるTroullier–Martins(TM)法を，係数決定の連立条件まで省略せずに追う．

15.5.1 擬波動関数に課す条件

まず，作るべき動径擬波動関数 $u_l^{\rm PS}(r) = r R_l^{\rm PS}(r)$ に何を要求するかを，理由とともに列挙する．1〜4は15.4節のノルム保存条件そのものだが，5と，4を「4階微分まで」と具体化する部分がTM法の特徴である．

定義:TM擬波動関数が満たすべき5条件

$r_c$ の外で全電子波動関数と一致する: $r \ge r_c$ で $u_l^{\rm PS}(r) = u_l^{\rm AE}(r)$．これにより，$r_c$ の外では真の原子と同一の物理が保たれ，参照エネルギー $\varepsilon_l$ での対数微分が厳密に一致する．
節を持たない: $0 < r < r_c$ で $u_l^{\rm PS}(r) \ne 0$．節があればそこで波動関数は急峻に符号を変え，滑らかさという目的そのものが失われる．また節がなければ，その擬波動関数はその $l$ チャネルの最低固有状態になり(節の数と固有値の順序の対応)，自己無撞着計算で安定に得られる．
ノルム保存: $\int_0^{r_c}|u_l^{\rm PS}|^2\dd r = \int_0^{r_c}|u_l^{\rm AE}|^2\dd r$．恒等式 \eqref{eq:15-normid} により，対数微分のエネルギー1次微分まで一致する(15.4節)．
$r_c$ で4階微分まで連続: $\left.\dfrac{\dd^k u_l^{\rm PS}}{\dd r^k}\right|_{r_c} = \left.\dfrac{\dd^k u_l^{\rm AE}}{\dd r^k}\right|_{r_c}$($k = 0,1,2,3,4$)．
原点で正則: $r\to0$ で $u_l^{\rm PS}(r) \propto r^{l+1}$，かつ逆算した擬ポテンシャルが原点で発散しないこと．TM法ではさらに $\left.\dfrac{\dd^2 v_l^{\rm scr}}{\dd r^2}\right|_{r=0} = 0$ という平滑条件を課す．

条件4の「4階微分まで」という数字の出どころを説明しておく．あとで見るように，擬ポテンシャルは $v \sim \varepsilon_l - \frac{l(l+1)}{2r^2} + \frac{u''}{2u}$ という形で波動関数の2階微分から逆算される．したがって

$v_l^{\rm PS}$ が $r_c$ で連続 $\Leftarrow$ $u$ の2階微分が連続，
$v_l^{\rm PS}{}'$ が連続 $\Leftarrow$ $u$ の3階微分が連続，
$v_l^{\rm PS}{}''$ が連続 $\Leftarrow$ $u$ の4階微分が連続

という対応になる．ポテンシャルの $n$ 階微分が不連続だと，そのフーリエ変換は $\tilde v(q)\sim q^{-(n+1)}$ という遅い減衰しかせず，平面波カットオフを上げても収束が遅い．4階微分まで揃えれば $v_l^{\rm PS}$ は2階微分まで連続になり，$\tilde v(q)$ は $q^{-4}$ より速く落ちる．条件4は「滑らかさ」を定量的に要求する条件なのである．

15.5.2 多項式アンザッツ

TM法は，上の5条件を同時に満たすように，次の関数形(アンザッツ)を置く．

$$ \begin{equation} u_l^{\rm PS}(r) = \begin{cases} u_l^{\rm AE}(r), & r \ge r_c,\\[4pt] r^{\,l+1}\exp\left[p(r)\right], & r \le r_c, \end{cases} \qquad p(r) = c_0 + c_2 r^2 + c_4 r^4 + c_6 r^6 + c_8 r^8 + c_{10}r^{10} + c_{12}r^{12} \label{eq:15-tm-ansatz} \end{equation} $$

この形が選ばれた理由を一つずつ確認する．

因子 $r^{l+1}$: 動径方程式 \eqref{eq:15-radial} の原点正則解は $u \propto r^{l+1}$ である(遠心力項 $l(l+1)/2r^2$ と2階微分項が釣り合う指数)．この因子を最初から括り出しておけば条件5の前半が自動的に満たされる．
指数関数 $\exp[p(r)]$: 指数関数は決してゼロにならない．したがって $0 < r < r_c$ で $u_l^{\rm PS} \ne 0$，すなわち条件2(節を持たない)が構成上自動的に保証される．もし多項式そのもので $u$ を表すと，係数を決めた後で節ができていないかを毎回検査しなければならない．
偶数冪だけ: $r$ の奇数冪($c_1 r$，$c_3 r^3$ など)を含めない．理由は次の数学ノートで述べる．
$r^{12}$ まで(係数7個): 後で数えるように，課す条件はちょうど7個である．未知数と方程式の数を合わせた結果が $c_0, c_2, \ldots, c_{12}$ の7係数なのである．Kerkerの原型は $r^4$ までの4係数で，条件4を2階微分までしか課していなかった．

数学ノート:なぜ $r$ の偶数冪だけなのか

理由は2つあり，どちらも「原点での滑らかさ」に関わる．

(i) 3次元関数としての滑らかさ．球対称な関数 $f(r)$ を3次元の関数 $f(\abs{\rr})$ として見ると，$r = \sqrt{x^2+y^2+z^2}$ は原点で微分不可能である(円錐の頂点のような尖りがある)．一方 $r^2 = x^2+y^2+z^2$ は $x,y,z$ の多項式であり，原点でも解析的である．したがって$f$ が原点で解析的であるためには，$f$ が $r^2$ の冪級数でなければならない．$r$ の奇数冪が混じると擬ポテンシャルは核の位置に尖り(cusp)を持ち，そのフーリエ変換は冪的にしか減衰しなくなる．平面波計算にとってこれは致命的である．

(ii) 逆算式に現れる $p'/r$ の正則性．次項で見るように，擬ポテンシャルには $\dfrac{(l+1)p'(r)}{r}$ という項が現れる．$p$ が偶数冪のみなら $p'(r) = 2c_2 r + 4c_4 r^3 + \cdots$ は $r$ の1次以上で消えるので，$p'/r = 2c_2 + 4c_4 r^2 + \cdots$ は原点で有限である．もし $c_1 r$ という項があれば $p' \to c_1 \ne 0$ となり，$p'/r \sim c_1/r$ が発散してしまう．

15.5.3 ポテンシャルの逆算

アンザッツ \eqref{eq:15-tm-ansatz} を固有関数として持つポテンシャルを求める．やることは単純で，動径方程式をポテンシャルについて解き直すだけである．

導出:遮蔽擬ポテンシャルの表式

ステップ1:動径方程式を $V$ について解く．式 \eqref{eq:15-radial} で $\varepsilon = \varepsilon_l$，$u = u_l^{\rm PS}$ とし，$V$ を含む項だけを左辺に残す:

$$ \begin{align} -\frac{1}{2}u'' + \frac{l(l+1)}{2r^2}u + V u &= \varepsilon_l u \label{eq:15-tm-inv1}\\ V(r) &= \varepsilon_l - \frac{l(l+1)}{2r^2} + \frac{u''(r)}{2u(r)}. \label{eq:15-tm-inv2} \end{align} $$

2行目へは，両辺を $u$ で割り($u \ne 0$ が条件2で保証されているからこの割り算が許される—ここで節を持たないという条件が本質的に効いている)，$V$ 以外を右辺に移しただけである．節があればその点で右辺が発散し，逆算が破綻する．

ステップ2:アンザッツの1階微分．$u = r^{l+1}e^{p(r)}$ に積の微分則を適用する:

$$ \begin{align} u'(r) &= (l+1)r^{l}e^{p} + r^{l+1}p'(r)e^{p} \label{eq:15-tm-d1}\\ &= r^{l+1}e^{p}\left[\frac{l+1}{r} + p'\right] = u(r)\left[\frac{l+1}{r} + p'(r)\right]. \label{eq:15-tm-d1b} \end{align} $$

1行目は $r^{l+1}$ と $e^p$ それぞれを微分した2項(後者には連鎖律で $p'$ が出る)．2行目では共通因子 $r^{l+1}e^p = u$ を括り出した．

ステップ3:2階微分．式 \eqref{eq:15-tm-d1b} をもう一度微分する．$u' = u\,g$，$g \equiv \frac{l+1}{r}+p'$ と書けば $u'' = u'g + ug'$ であり，$u' = ug$ を代入すると $u'' = u g^2 + u g'$，すなわち

$$ \begin{align} \frac{u''(r)}{u(r)} &= \left[\frac{l+1}{r} + p'\right]^2 + \frac{\dd}{\dd r}\left[\frac{l+1}{r} + p'\right] \label{eq:15-tm-d2a}\\ &= \frac{(l+1)^2}{r^2} + \frac{2(l+1)p'}{r} + (p')^2 - \frac{l+1}{r^2} + p'' \label{eq:15-tm-d2b}\\ &= \frac{(l+1)\left[(l+1)-1\right]}{r^2} + \frac{2(l+1)p'}{r} + (p')^2 + p'' \label{eq:15-tm-d2c}\\ &= \frac{l(l+1)}{r^2} + \frac{2(l+1)p'}{r} + (p')^2 + p''. \label{eq:15-tm-d2d} \end{align} $$

2行目は括弧の2乗を展開($(a+b)^2 = a^2+2ab+b^2$)し，$\frac{\dd}{\dd r}\frac{l+1}{r} = -\frac{l+1}{r^2}$ を使った．3行目は $1/r^2$ の係数をまとめ，4行目で $(l+1)l = l(l+1)$ とした．

ステップ4:代入して遠心力項を消す．式 \eqref{eq:15-tm-d2d} を \eqref{eq:15-tm-inv2} に入れる:

$$ \begin{align} v_l^{\rm scr}(r) &= \varepsilon_l - \frac{l(l+1)}{2r^2} + \frac{1}{2}\left[\frac{l(l+1)}{r^2} + \frac{2(l+1)p'}{r} + (p')^2 + p''\right] \label{eq:15-tm-v1}\\ &= \varepsilon_l + \frac{(l+1)\,p'(r)}{r} + \frac{1}{2}\left[p''(r) + \left(p'(r)\right)^2\right]. \label{eq:15-tm-v2} \end{align} $$

2行目では，$-\frac{l(l+1)}{2r^2}$ と $+\frac{1}{2}\cdot\frac{l(l+1)}{r^2}$ が完全に相殺した．これは偶然ではなく，$u \propto r^{l+1}$ という因子を最初に括り出しておいたおかげで遠心力項がちょうど吸収されたのである．残った表式には $1/r^2$ の発散が一切ない．

結果を掲げる．$v_l^{\rm scr}$ を「遮蔽(screened)擬ポテンシャル」と呼ぶ理由は15.5.5項で説明する．

$$ \begin{equation} v_l^{\rm scr}(r) = \varepsilon_l + \frac{(l+1)\,p'(r)}{r} + \frac{p''(r) + \left[p'(r)\right]^2}{2} \qquad (r \le r_c) \label{eq:15-tm-vscr} \end{equation} $$

この式から，原点での値と平滑条件が直ちに読み取れる．$p' = 2c_2 r + 4c_4 r^3 + 6c_6 r^5 + \cdots$，$p'' = 2c_2 + 12c_4 r^2 + 30c_6r^4 + \cdots$，$(p')^2 = 4c_2^2r^2 + \cdots$ を \eqref{eq:15-tm-vscr} に代入して $r$ の冪で整理すると

$$ \begin{align} v_l^{\rm scr}(r) &= \varepsilon_l + (l+1)\left[2c_2 + 4c_4r^2 + \cdots\right] + \frac{1}{2}\left[2c_2 + 12c_4r^2+\cdots\right] + \frac{1}{2}\left[4c_2^2r^2+\cdots\right] \label{eq:15-tm-exp1}\\ &= \underbrace{\varepsilon_l + (2l+3)c_2}_{r^0} + \underbrace{2\left[(2l+5)c_4 + c_2^2\right]}_{r^2\text{ の係数}}\,r^2 + O(r^4). \label{eq:15-tm-exp2} \end{align} $$

2行目の定数項は $2(l+1)c_2 + c_2 = (2l+3)c_2$，$r^2$ の係数は $4(l+1)c_4 + 6c_4 + 2c_2^2 = (4l+10)c_4 + 2c_2^2 = 2\left[(2l+5)c_4+c_2^2\right]$ とまとめた．第一の帰結は

$$ \begin{equation} v_l^{\rm scr}(0) = \varepsilon_l + (2l+3)\,c_2 \label{eq:15-tm-v0} \end{equation} $$

で，擬ポテンシャルは原点で有限である($-Z/r$ の発散が消えた!)．第二の帰結は，$r^2$ の係数を消せば $v_l^{\rm scr}$ の原点での2階微分がゼロになるということである．$v(r) = v(0) + \frac{1}{2}v''(0)r^2 + \cdots$ と比べれば $v''(0) = 4\left[(2l+5)c_4+c_2^2\right]$ だから，TMの平滑条件は

$$ \begin{equation} c_2^2 + (2l+5)\,c_4 = 0 \label{eq:15-tm-smooth} \end{equation} $$

と書ける．これは「原点近傍で擬ポテンシャルをできるだけ平らにする」ための条件であり，TM擬ポテンシャルが同種の他の構成法より小さい $E_{\rm cut}$ で収束する主因である．$r^4$ の項までは消せない(係数の自由度をすべて使い切っているため)ことに注意しよう．

15.5.4 係数を決める連立条件

7個の係数 $c_0, c_2, c_4, c_6, c_8, c_{10}, c_{12}$ を決める7個の条件を書き下す．

導出:7つの条件の具体形

(A) $r_c$ での0〜4階微分の一致(5条件)．$u^{\rm PS}$ の対数をとると $\ln u^{\rm PS} = (l+1)\ln r + p(r)$ だから，$u$ の $k$ 階微分を揃えることは $p$ の $k$ 階微分を揃えることと同値である(逆に解ける関係で結ばれている)．順に見る．

$k=0$:値の一致 $u^{\rm PS}(r_c) = u^{\rm AE}(r_c)$ から，両辺の対数をとって

$$ \begin{equation} p(r_c) = \sum_{i=0}^{6}c_{2i}\,r_c^{2i} = \ln\left[\frac{u_l^{\rm AE}(r_c)}{r_c^{\,l+1}}\right]. \label{eq:15-tm-c0} \end{equation} $$

$k=1$:式 \eqref{eq:15-tm-d1b} を $u'/u$ の形にすると $\frac{u'}{u} = \frac{l+1}{r} + p'$ だから，$r_c$ での対数微分 $D_l^{\rm AE}(\varepsilon_l, r_c) = u^{\rm AE}{}'(r_c)/u^{\rm AE}(r_c)$ を使って

$$ \begin{equation} p'(r_c) = D_l^{\rm AE}(\varepsilon_l, r_c) - \frac{l+1}{r_c}. \label{eq:15-tm-c1} \end{equation} $$

$k=2,3,4$:これらは，逆算式 \eqref{eq:15-tm-vscr} を $p''$ について解いた式と，その微分から得られる．まず \eqref{eq:15-tm-vscr} を $p''$ について解くと

$$ \begin{equation} p'' = 2\left(v - \varepsilon_l\right) - \frac{2(l+1)p'}{r} - (p')^2. \label{eq:15-tm-c2} \end{equation} $$

これを $r$ で1回，2回微分すると(商の微分則 $\left(\frac{p'}{r}\right)' = \frac{p''}{r} - \frac{p'}{r^2}$ を繰り返し使う)

$$ \begin{align} p''' &= 2v' - \frac{2(l+1)p''}{r} + \frac{2(l+1)p'}{r^2} - 2p'p'', \label{eq:15-tm-c3}\\ p'''' &= 2v'' - \frac{2(l+1)p'''}{r} + \frac{4(l+1)p''}{r^2} - \frac{4(l+1)p'}{r^3} - 2\left(p''\right)^2 - 2p'p'''. \label{eq:15-tm-c4} \end{align} $$

ここで $v, v', v''$ には $r_c$ における全電子の遮蔽ポテンシャルの値とその微分を代入する($r \ge r_c$ で両者は一致するから)．右辺は既知量と，すでに決まった低階の $p^{(k)}(r_c)$ だけで書けているので，$p''(r_c)$，$p'''(r_c)$，$p''''(r_c)$ が順に確定する．

そして左辺の $p^{(k)}(r_c)$ は多項式の微分だから，係数について線形である:

$$ p^{(k)}(r_c) = \sum_{i=1}^{6}\frac{(2i)!}{(2i-k)!}\,c_{2i}\,r_c^{\,2i-k}\qquad (k=1,2,3,4). $$

(B) ノルム保存(1条件)．アンザッツを代入すると

$$ \begin{align} \int_0^{r_c}\left|u_l^{\rm PS}\right|^2\dd r &= \int_0^{r_c} r^{2l+2}e^{2p(r)}\dd r \label{eq:15-tm-n1}\\ &= e^{2c_0}\int_0^{r_c} r^{2l+2}\exp\left[2\left(p(r)-c_0\right)\right]\dd r \label{eq:15-tm-n2} \end{align} $$

となる．2行目では $e^{2p} = e^{2c_0}e^{2(p-c_0)}$ と定数因子を括り出した．$p - c_0$ は $c_0$ を含まないので，右辺の積分は $c_0$ に依存しない．両辺の対数をとって全電子のノルムと等置すれば

$$ \begin{equation} \ln\left[\int_0^{r_c}\left|u_l^{\rm AE}(r)\right|^2\dd r\right] = 2c_0 + \ln\left[\int_0^{r_c} r^{2l+2}\exp\left(2\left[p(r)-c_0\right]\right)\dd r\right] \label{eq:15-tm-norm} \end{equation} $$

を得る．対数の形にしておく実務的な利点は2つある．(i) $e^{2p}$ は $p \sim -10$ 程度になることもあり，そのまま指数を計算すると桁あふれ・桁落ちを起こすが，$p - c_0$ なら $r_c$ 近傍で $O(1)$ に収まる．(ii) $c_0$ が右辺に陽に分離されるので，反復更新の式として直接使える．

(C) 原点での平滑条件(1条件)．式 \eqref{eq:15-tm-smooth}:$c_2^2 + (2l+5)c_4 = 0$．

条件は (A) 5個 $+$ (B) 1個 $+$ (C) 1個 $=$ 7個，未知数も7個で数が合う．ただし (B) は $c_0$ について超越的，(C) は $c_2$ について2次であり，系全体は非線形である．TMが与えた実際の解法は次の入れ子反復である．

定義:TM係数決定のアルゴリズム

$c_2$ に試行値を与える(外側ループ)．
$c_0$ に試行値を与える(内側ループ)．
条件 \eqref{eq:15-tm-c0} と $k=1,2,3,4$ の4条件を合わせた5元1次連立方程式を，未知数 $c_4, c_6, c_8, c_{10}, c_{12}$ について解く($c_0, c_2$ は既知として右辺に移す)．
得られた $p(r)$ で \eqref{eq:15-tm-norm} の右辺の積分を数値的に評価し，$c_0$ を更新する．3〜4を収束するまで繰り返す．
収束した $c_4$ を使って残差 $F(c_2) \equiv c_2^2 + (2l+5)c_4$ を評価する．$F(c_2) = 0$ となるまで $c_2$ を二分法(またはNewton法)で更新し，1に戻る．

収束したら，式 \eqref{eq:15-tm-vscr} で $r \le r_c$ の $v_l^{\rm scr}$ を，$r \ge r_c$ では全電子の遮蔽ポテンシャルをそのまま使って，$l$ ごとの遮蔽擬ポテンシャルが完成する．ここまでの作業は $l = 0, 1, 2, \ldots$ について独立に行うので，得られるのは$l$ ごとに異なるポテンシャルの組 $\{v_0^{\rm scr}, v_1^{\rm scr}, v_2^{\rm scr}, \ldots\}$ である．この「$l$ 依存性」の扱いが15.6節の主題になる．

15.5.5 スクリーニング解除(unscreening)

いま得た $v_l^{\rm scr}$ には，まだ使えない理由がある．逆算の出発点にした全電子原子の動径方程式のポテンシャルは，コーン・シャム有効ポテンシャル(第10章)

$$ \begin{equation} v^{\rm AE}(r) = -\frac{Z}{r} + v_{\rm H}\!\left[\rho^{\rm AE}\right]\!(r) + v_{\rm xc}\!\left[\rho^{\rm AE}\right]\!(r) \label{eq:15-ae-veff} \end{equation} $$

であって，核の引力だけでなくその原子の電子自身が作るHartree項と交換相関項を含んでいる．つまり $v_l^{\rm scr}$ は「参照原子の価電子分布によって遮蔽された(screened)」ポテンシャルである．これを別の環境(分子・固体)に持ち込むと，その環境での価電子密度から作られるHartree項・交換相関項が自己無撞着計算の中で改めて足されるので，参照原子の価電子の寄与が二重に数えられてしまう．

そこで，参照原子の価電子密度 $\rho_v$ が作る寄与をあらかじめ差し引いておく．これがスクリーニング解除である．

$$ \begin{equation} v_l^{\rm ion}(r) = v_l^{\rm scr}(r) - v_{\rm H}\!\left[\rho_v\right]\!(r) - v_{\rm xc}\!\left[\rho_v\right]\!(r) \label{eq:15-unscreen} \end{equation} $$

ここで $\rho_v(\rr) = \sum_{l} q_l \left|R_l^{\rm PS}(r)\right|^2/(4\pi)$ は擬価電子密度($q_l$ は占有数)であり，全電子の価電子密度ではないことに注意する．ノルム保存のおかげで $r \ge r_c$ では両者は一致し，$r < r_c$ でも球内の総電荷は一致する．得られた $v_l^{\rm ion}$ をイオン(裸の)擬ポテンシャルと呼び，これが元素ごとに1度だけ作られてファイルに保存され，あらゆる計算に使い回される量である．

物理的意味:$v^{\rm ion}$ の遠方形と価電荷

$r$ が大きいところで $v_l^{\rm ion}(r) \to -Z_v/r$ となる($Z_v$ は価電子数)．理由は次のとおりである．$r \gg r_c$ では $v^{\rm scr}$ は中性原子のコーン・シャム有効ポテンシャルであり，核 $+Z$ と全電子 $-Z$ が打ち消して急速にゼロに近づく．一方 $v_{\rm H}[\rho_v] \to Z_v/r$(球外から見た価電荷 $Z_v$ の作るポテンシャル)，$v_{\rm xc}[\rho_v]\to 0$(密度が指数的にゼロになるので)である．したがって \eqref{eq:15-unscreen} より $v_l^{\rm ion} \to 0 - Z_v/r - 0 = -Z_v/r$．擬原子は，遠くから見ると電荷 $+Z_v$ の点電荷に見える—これは物理的に当然の要請であり，実装のチェックに使える(炭素なら $-4/r$，シリコンなら $-4/r$，酸素なら $-6/r$)．

図15.5 Troullier–Martins擬ポテンシャルの模式図．全電子の有効ポテンシャル(赤・破線)は原点で $-Z/r$ として発散するのに対し，遮蔽擬ポテンシャル(青)は原点で有限値 $\varepsilon_l + (2l+3)c_2$ をとり，$r \ge r_c$ では全電子と厳密に一致する．スクリーニング解除して得たイオン擬ポテンシャル(緑)は遠方で $-Z_v/r$ に漸近する．

15.5.6 非線形コア補正の必要性

式 \eqref{eq:15-unscreen} には隠れた欠陥がある．Hartree項は密度の線形汎関数だから

$$ \begin{equation} v_{\rm H}\!\left[\rho_v + \rho_c\right] = v_{\rm H}\!\left[\rho_v\right] + v_{\rm H}\!\left[\rho_c\right] \label{eq:15-vh-linear} \end{equation} $$

が厳密に成り立つ($v_{\rm H}[\rho](\rr) = \int\dd^3r'\,\rho(\rr')/\abs{\rr-\rr'}$ は $\rho$ について線形)．したがって内殻電子のHartree寄与 $v_{\rm H}[\rho_c]$ は $v_l^{\rm ion}$ の中にきれいに畳み込まれ，環境が変わっても正しく働く．ところが交換相関ポテンシャルは密度の非線形汎関数であり

$$ \begin{equation} v_{\rm xc}\!\left[\rho_v + \rho_c\right] \ne v_{\rm xc}\!\left[\rho_v\right] + v_{\rm xc}\!\left[\rho_c\right] \label{eq:15-vxc-nonlinear} \end{equation} $$

である(たとえばLDA交換では $v_{\rm x} \propto \rho^{1/3}$ で，$(\rho_v+\rho_c)^{1/3} \ne \rho_v^{1/3}+\rho_c^{1/3}$)．誤差の大きさを見積もるために $\rho_c$ が小さいとしてTaylor展開すると

$$ \begin{align} v_{\rm xc}\!\left[\rho_v+\rho_c\right](\rr) &= v_{\rm xc}\!\left[\rho_v\right](\rr) + \int\dd^3r'\,\frac{\delta v_{\rm xc}(\rr)}{\delta\rho(\rr')}\,\rho_c(\rr') + O(\rho_c^2) \label{eq:15-nlcc-taylor} \end{align} $$

となり，線形項ですら $\rho_c$ に比例する有限の寄与を与える．この誤差が効くのは$\rho_c$ と $\rho_v$ が同程度の大きさを持つ領域である．核のごく近傍では $\rho_c \gg \rho_v$ で価電子は存在しないから問題にならず，遠方では $\rho_c \simeq 0$ だから問題にならない．危険なのはその中間，内殻の裾と価電子の内側が重なる $r \sim 0.3$–$1$ Bohr の殻である．アルカリ金属，3d遷移金属，そしてスピン分極計算(交換相関のスピン依存性がさらに非線形性を増幅する)で誤差が顕在化する．

対策がLouie–Froyen–Cohen (1982) の非線形コア補正(nonlinear core correction, NLCC;部分内殻補正 PCC ともいう)である．処方は単純で，内殻密度 $\rho_c$ を，核近傍だけ滑らかな関数で置き換えた部分内殻密度 $\tilde\rho_c$ を用意し

$$ \begin{equation} v_l^{\rm ion}(r) = v_l^{\rm scr}(r) - v_{\rm H}\!\left[\rho_v\right]\!(r) - v_{\rm xc}\!\left[\rho_v + \tilde\rho_c\right]\!(r) \label{eq:15-unscreen-nlcc} \end{equation} $$

としてスクリーニング解除し，その後の分子・固体計算でも交換相関エネルギーを $E_{\rm xc}[\rho + \tilde\rho_c]$ として評価する．こうすれば，非線形性の効く領域で内殻密度が常に「見えている」ことになり，参照原子と実際の系で同じ扱いになる．$\tilde\rho_c$ を真の $\rho_c$ でなく平滑化したもので置き換えてよいのは，核の直近では $\rho_v \simeq 0$ なので非線形性がそもそも効かないからである(平滑化しないと $\rho_c$ 自体が急峻すぎて平面波で表せない—本末転倒になる)．

例:炭素原子のTM擬ポテンシャル

炭素($1s^2 2s^2 2p^2$，$Z=6$，$Z_v=4$)を例にとる．内殻は $1s$ のみ，価電子は $2s$ と $2p$ である．$2s$ チャネルでは $r_c \approx 1.3$–$1.5$ Bohr にとる(全電子 $2s$ の最外節が $r\approx0.4$ Bohr，最大値が $r\approx1.2$ Bohr にあるので，その少し外側)．得られる $v_0^{\rm scr}$ の原点での値は式 \eqref{eq:15-tm-v0} より $\varepsilon_{2s} + 3c_2$ で，数値的に $-4$ Ha 程度の有限値になる．$2p$ チャネルでは全電子波動関数がそもそも節を持たない($2p$ は $l=1$ の最低状態)ので，擬波動関数は全電子波動関数を「少し外側に広げただけ」の形になり，$v_1^{\rm scr}$ は $v_0^{\rm scr}$ より浅い．イオン擬ポテンシャルはどちらも遠方で $-4/r$ に漸近する．部分内殻密度 $\tilde\rho_c$ は $r\approx0.35$ Bohr にピークを持ち，$r\approx1$ Bohr でほぼ消える．炭素は $1s$ が非常にコンパクトなので，NLCCなしでも実用上問題ないことが多いが，スピン分極を扱う場合には入れておくのが安全である．

これで，$l$ ごとの滑らかなイオン擬ポテンシャル $v_l^{\rm ion}(r)$ が手に入った．次の問題は，この「$l$ ごとに違うポテンシャル」を実際の計算にどう乗せるかである．

15.6 セミローカル形とKleinman–Bylander分離形

前節の構成法は角運動量 $l$ ごとに独立に走らせるので，成果物は1つのポテンシャルではなくポテンシャルの組 $\{v_0^{\rm ion}, v_1^{\rm ion}, v_2^{\rm ion},\ldots\}$ である．$l$ ごとに違うポテンシャルを1つの演算子にまとめるには射影演算子が要る．この節では，素直にまとめた「セミローカル形」がなぜ計算量的に高くつくのか，Kleinman–Bylander (1982) の変換がそれをどう解消するのか，そしてその代償として現れる「ゴースト状態」という病理を扱う．

15.6.1 セミローカル形

角運動量 $l$，磁気量子数 $m$ の球面調和関数 $Y_{lm}$ への射影演算子を

$$ \begin{equation} \hat P_{lm} = \left|Y_{lm}\right\rangle\left\langle Y_{lm}\right|, \qquad \left\langle \rr \middle| \hat P_{lm} \middle| \psi \right\rangle = Y_{lm}(\hat\rr)\int \dd\hat\rr'\, Y_{lm}^*(\hat\rr')\,\psi(r,\hat\rr') \label{eq:15-proj} \end{equation} $$

と書く($\hat\rr$ は方向単位ベクトル，$\int\dd\hat\rr$ は立体角積分)．角度部分にだけ作用し，動径座標 $r$ には触れないことに注意．これを使うと，$l$ ごとのポテンシャルは

$$ \begin{equation} \hat v^{\rm SL} = \sum_{l=0}^{\infty}\sum_{m=-l}^{l}\left|Y_{lm}\right\rangle v_l^{\rm ion}(r)\left\langle Y_{lm}\right| \label{eq:15-semilocal} \end{equation} $$

という1つの演算子にまとめられる．これをセミローカル(semilocal)ポテンシャルという．「セミ」なのは，動径方向には局所的(掛け算)だが角度方向には非局所的(積分核)だからである．位置表示では

$$ \begin{equation} \left\langle\rr\right|\hat v^{\rm SL}\left|\rr'\right\rangle = \sum_{lm} Y_{lm}(\hat\rr)\,\frac{\delta(r-r')}{r^2}\,v_l^{\rm ion}(r)\,Y_{lm}^*(\hat\rr') \label{eq:15-semilocal-kernel} \end{equation} $$

となる($\delta(r-r')/r^2$ は3次元デルタ関数の動径部分)．

実際には無限個の $l$ を扱う必要はない．$l$ が大きいと遠心力障壁 $l(l+1)/2r^2$ が電子を核から遠ざけるので，その電子は内殻をまったく感じない．したがって $l \gtrsim l_{\max}+1$ ではすべての $v_l^{\rm ion}$ が共通の形に収束する．そこで，ある $l$($l = l_{\rm loc}$，または適当に平滑化した関数)を局所ポテンシャル $v_{\rm loc}(r)$ として選び出し，差だけを非局所項として残す:

$$ \begin{equation} \hat v^{\rm SL} = v_{\rm loc}(r) + \sum_{l=0}^{l_{\max}}\sum_{m=-l}^{l}\left|Y_{lm}\right\rangle \delta v_l(r)\left\langle Y_{lm}\right|, \qquad \delta v_l(r) \equiv v_l^{\rm ion}(r) - v_{\rm loc}(r). \label{eq:15-sl-split} \end{equation} $$

$\sum_{lm}|Y_{lm}\rangle\langle Y_{lm}| = 1$(完全性)を使えば \eqref{eq:15-semilocal} と \eqref{eq:15-sl-split} が同じ演算子であることが確かめられる．ここで重要なのは，$\delta v_l(r)$ が $r \ge r_c$ でゼロになることである($r_c$ の外では全電子ポテンシャルと一致するので $l$ に依らない)．非局所項は半径 $r_c$ の球の中だけで働く短距離演算子である．

15.6.2 平面波基底でのコスト

セミローカル形の何が問題か．平面波基底での行列要素を計算してみればわかる．

数学ノート:Rayleighの平面波展開と球面調和関数の加法定理

以下で2つの公式を使う．第一はRayleighの展開 $$e^{i\bm q\cdot\rr} = 4\pi\sum_{l=0}^{\infty}\sum_{m=-l}^{l} i^{\,l}\, j_l(qr)\, Y_{lm}^*(\hat{\bm q})\,Y_{lm}(\hat\rr)$$ である．平面波を，原点まわりの球面波(球ベッセル関数 $\times$ 球面調和関数)で展開したものにほかならない．角度部分について両辺に $Y_{lm}^*(\hat\rr)$ を掛けて立体角積分すると，$Y$ の正規直交性から $$\int\dd\hat\rr\; e^{i\bm q\cdot\rr}\,Y_{lm}^*(\hat\rr) = 4\pi i^{\,l} j_l(qr)Y_{lm}^*(\hat{\bm q}), \qquad \int\dd\hat\rr\; e^{-i\bm q\cdot\rr}\,Y_{lm}(\hat\rr) = 4\pi (-i)^{\,l} j_l(qr)Y_{lm}(\hat{\bm q})$$ が得られる(2番目は1番目の複素共役をとって $Y_{lm}^* \to Y_{lm}$ と読み替えた)．

第二は加法定理 $$\sum_{m=-l}^{l} Y_{lm}(\hat{\bm q})\,Y_{lm}^*(\hat{\bm q}') = \frac{2l+1}{4\pi}\,P_l(\hat{\bm q}\cdot\hat{\bm q}')$$ である($P_l$ はLegendre多項式)．左辺は2つの方向の間の回転不変な量なので，両者のなす角の関数でなければならず，$l$ 次であることからLegendre多項式に比例する—係数は $\hat{\bm q}=\hat{\bm q}'$ とおいて $\sum_m|Y_{lm}|^2 = (2l+1)/4\pi$，$P_l(1)=1$ と比べれば決まる．

導出:セミローカル形の平面波行列要素

体積 $\Omega$ の単位胞で規格化した平面波 $\left\langle\rr|\bm q\right\rangle = \Omega^{-1/2}e^{i\bm q\cdot\rr}$ をとる．非局所項の行列要素は，核 \eqref{eq:15-semilocal-kernel} を挟んで

$$ \begin{align} \left\langle\bm q\right|\delta\hat v\left|\bm q'\right\rangle &= \frac{1}{\Omega}\sum_{lm}\int_0^{\infty}\!\! r^2\dd r\;\delta v_l(r) \left[\int\dd\hat\rr\, e^{-i\bm q\cdot\rr}Y_{lm}(\hat\rr)\right] \left[\int\dd\hat\rr'\, e^{i\bm q'\cdot\rr'}Y_{lm}^*(\hat\rr')\right] \label{eq:15-sl-me1}\\ &= \frac{1}{\Omega}\sum_{lm}\int_0^{\infty}\!\! r^2\dd r\;\delta v_l(r) \left[4\pi(-i)^l j_l(qr)Y_{lm}(\hat{\bm q})\right]\left[4\pi i^{\,l} j_l(q'r)Y_{lm}^*(\hat{\bm q}')\right] \label{eq:15-sl-me2}\\ &= \frac{16\pi^2}{\Omega}\sum_{l}\left[\sum_m Y_{lm}(\hat{\bm q})Y_{lm}^*(\hat{\bm q}')\right] \int_0^{r_c}\!\! r^2\dd r\; j_l(qr)\,\delta v_l(r)\, j_l(q'r) \label{eq:15-sl-me3}\\ &= \frac{4\pi}{\Omega}\sum_{l=0}^{l_{\max}}(2l+1)\,P_l(\hat{\bm q}\cdot\hat{\bm q}') \int_0^{r_c}\!\! r^2\dd r\; j_l(qr)\,\delta v_l(r)\, j_l(q'r). \label{eq:15-sl-me4} \end{align} $$

1行目は核を代入し，動径のデルタ関数 $\delta(r-r')/r^2$ で $r'$ 積分を実行した(残る $r^2\dd r$ は3次元体積要素の動径部分)．2行目は数学ノートの2つの角度積分を代入した．3行目では $(-i)^l i^{\,l} = 1$ を使い，$\delta v_l$ が $r\ge r_c$ で消えることから積分区間を $[0,r_c]$ に縮めた．4行目は加法定理で $m$ 和を潰し，$16\pi^2 \times \frac{2l+1}{4\pi} = 4\pi(2l+1)$ とまとめた．

ここで注目すべきは，式 \eqref{eq:15-sl-me4} の動径積分が $j_l(qr)$ と $j_l(q'r)$ を同じ積分の中に抱えていることである．すなわち行列要素は

$$ \left\langle\bm q\right|\delta\hat v\left|\bm q'\right\rangle \ne \sum_{\alpha}\left(\text{$\bm q$ だけの関数}\right)_\alpha\times\left(\text{$\bm q'$ だけの関数}\right)_\alpha $$

という形に因子化できない．したがって $N_{\rm PW}$ 個の平面波に対して $N_{\rm PW}^2$ 個の独立な行列要素をすべて計算・保持しなければならず，ハミルトニアンをベクトルに作用させる操作 $\hat H\psi$ は1バンドあたり $O(N_{\rm PW}^2)$ の演算を要する．15.1節の見積りのように $N_{\rm PW}\sim 10^4$–$10^5$ ともなれば，この項だけで計算全体を支配してしまう．

15.6.3 Kleinman–Bylander変換

KleinmanとBylander (1982) の解決は鮮やかである．非局所項を，参照擬波動関数から作った1本のベクトルの射影(ランク1演算子)で置き換える．各 $(l,m)$ チャネルについて，擬原子の参照擬波動関数を $\phi_{lm}(\rr) = R_l^{\rm PS}(r)Y_{lm}(\hat\rr)$ として

$$ \begin{equation} \hat v^{\rm KB} = v_{\rm loc}(r) + \sum_{l=0}^{l_{\max}}\sum_{m=-l}^{l} \frac{\left|\delta v_l\,\phi_{lm}\right\rangle\left\langle \delta v_l\,\phi_{lm}\right|} {\left\langle \phi_{lm}\right|\delta v_l\left|\phi_{lm}\right\rangle} \label{eq:15-kb} \end{equation} $$

と定義する．ここで $\left|\delta v_l\phi_{lm}\right\rangle$ は関数 $\delta v_l(r)\phi_{lm}(\rr)$ が表すケットである．分母 $E_l^{\rm KB} \equiv \left\langle\phi_{lm}\right|\delta v_l\left|\phi_{lm}\right\rangle = \int\dd^3r\,\left|\phi_{lm}\right|^2\delta v_l(r)$ はKBエネルギーと呼ばれ，$m$ に依らない実数である($\delta v_l$ が実の動径関数，$\int|Y_{lm}|^2\dd\hat\rr=1$ だから)．

定理:KB形は参照状態に対してセミローカル形と同一に作用する

任意の $(l',m')$ について $$\hat v^{\rm KB}\left|\phi_{l'm'}\right\rangle = \hat v^{\rm SL}\left|\phi_{l'm'}\right\rangle.$$ したがって擬原子の参照固有値 $\varepsilon_{l'}$ と参照擬波動関数 $\phi_{l'm'}$ は，KB形でも厳密に再現される．

導出:定理の証明

ステップ1:内積 $\left\langle\delta v_l\phi_{lm}\middle|\phi_{l'm'}\right\rangle$ を計算する．定義通りに書き下すと

$$ \begin{align} \left\langle\delta v_l\phi_{lm}\middle|\phi_{l'm'}\right\rangle &= \int\dd^3 r\;\left[\delta v_l(r)R_l^{\rm PS}(r)Y_{lm}(\hat\rr)\right]^{*}R_{l'}^{\rm PS}(r)Y_{l'm'}(\hat\rr) \label{eq:15-kb-p1}\\ &= \left[\int_0^{\infty}\!\! r^2\dd r\;\delta v_l(r)R_l^{\rm PS}(r)R_{l'}^{\rm PS}(r)\right] \left[\int\dd\hat\rr\; Y_{lm}^*(\hat\rr)Y_{l'm'}(\hat\rr)\right] \label{eq:15-kb-p2}\\ &= \delta_{ll'}\delta_{mm'}\int_0^{r_c}\!\! r^2\dd r\;\delta v_l(r)\left|R_l^{\rm PS}(r)\right|^2 = \delta_{ll'}\delta_{mm'}\,E_l^{\rm KB}. \label{eq:15-kb-p3} \end{align} $$

1行目は定義．2行目では被積分関数が動径部分と角度部分の積に分かれること($\delta v_l$ は $r$ のみの実関数，$R^{\rm PS}$ も実)を使って積分を分離した．3行目では球面調和関数の正規直交性 $\int\dd\hat\rr\,Y_{lm}^*Y_{l'm'} = \delta_{ll'}\delta_{mm'}$ を使い，クロネッカーデルタのおかげで $l'=l$ に固定されたので $R_{l'}^{\rm PS}=R_l^{\rm PS}$ とし，$\delta v_l$ の台が $[0,r_c]$ であることを使った．最後の等号は $E_l^{\rm KB}$ の定義そのものである．

ステップ2:KB演算子を作用させる．式 \eqref{eq:15-kb} の非局所部分を $\left|\phi_{l'm'}\right\rangle$ に当てる:

$$ \begin{align} \sum_{lm}\frac{\left|\delta v_l\phi_{lm}\right\rangle\left\langle\delta v_l\phi_{lm}\middle|\phi_{l'm'}\right\rangle}{E_l^{\rm KB}} &= \sum_{lm}\frac{\left|\delta v_l\phi_{lm}\right\rangle\,\delta_{ll'}\delta_{mm'}E_l^{\rm KB}}{E_l^{\rm KB}} \label{eq:15-kb-a1}\\ &= \left|\delta v_{l'}\phi_{l'm'}\right\rangle = \delta v_{l'}(r)\left|\phi_{l'm'}\right\rangle. \label{eq:15-kb-a2} \end{align} $$

1行目はステップ1の結果を代入しただけ．2行目ではクロネッカーデルタで和を潰し，$E_l^{\rm KB}$ が約分された．

ステップ3:セミローカル形と比べる．式 \eqref{eq:15-sl-split} を同じ状態に作用させると，射影 $\left|Y_{lm}\right\rangle\left\langle Y_{lm}\middle|\phi_{l'm'}\right\rangle$ が $\delta_{ll'}\delta_{mm'}$ を出すので

$$ \sum_{lm}\left|Y_{lm}\right\rangle\delta v_l(r)\left\langle Y_{lm}\middle|\phi_{l'm'}\right\rangle = \delta v_{l'}(r)\left|\phi_{l'm'}\right\rangle $$

となり，\eqref{eq:15-kb-a2} と一致する．局所部分 $v_{\rm loc}$ は両者共通なので，定理が示された．$\blacksquare$

何が得られたか．KB形は，参照状態という「1点」ではセミローカル形と完全に同じでありながら，演算子としてはランク1の分離形になっている．その利益は平面波行列要素に直ちに現れる:

$$ \begin{equation} \left\langle\bm q\right|\hat v^{\rm KB}_{\rm nl}\left|\bm q'\right\rangle = \sum_{lm}\frac{F_{lm}(\bm q)\,F_{lm}^*(\bm q')}{E_l^{\rm KB}}, \qquad F_{lm}(\bm q) \equiv \left\langle\bm q\middle|\delta v_l\,\phi_{lm}\right\rangle = \frac{4\pi(-i)^l}{\sqrt{\Omega}}Y_{lm}(\hat{\bm q})\!\int_0^{r_c}\!\!r^2\dd r\, j_l(qr)\delta v_l(r)R_l^{\rm PS}(r) \label{eq:15-kb-me} \end{equation} $$

($F_{lm}$ の具体形は，数学ノートの角度積分をそのまま使っただけである．)行列要素は $\bm q$ の関数と $\bm q'$ の関数の積に完全に因子化した．したがって，$N_{\rm PW}\times N_{\rm proj}$ 個の数 $F_{lm}(\bm q)$ だけを保持すればよく($N_{\rm proj} = \sum_{l\le l_{\max}}(2l+1)$ は普通10個以下)，$\hat H\psi$ の計算は

$$ \left(\hat v^{\rm KB}_{\rm nl}\psi\right)(\bm q) = \sum_{lm}\frac{F_{lm}(\bm q)}{E_l^{\rm KB}}\underbrace{\sum_{\bm q'}F_{lm}^*(\bm q')\psi(\bm q')}_{\text{先にこのスカラーを作る}} $$

という2段階で行える．内側の和が $O(N_{\rm PW})$，外側が $O(N_{\rm PW})$ で，合わせて $O(N_{\rm PW}N_{\rm proj})$—セミローカル形の $O(N_{\rm PW}^2)$ から $N_{\rm PW}/N_{\rm proj}$ 倍，実際には3〜4桁の高速化である．しかも記憶量も $N_{\rm PW}^2$ から $N_{\rm PW}N_{\rm proj}$ に落ちる．この一手によって，平面波DFTは大規模系に届く手法になった．

図15.6 セミローカル形とKleinman–Bylander分離形の計算コストの違い．セミローカル形の行列要素 \eqref{eq:15-sl-me4} は $q$ と $q'$ が動径積分の中で絡み合って因子化できないが，KB形 \eqref{eq:15-kb-me} は $F_{lm}(\bm q)F_{lm}^*(\bm q')$ という外積(ランク1)になり，行列を組む必要すらなくなる．

15.6.4 ゴースト状態

KB変換はただ乗りではない．定理が保証するのは参照状態に対してだけの一致であって，他の状態に対しては $\hat v^{\rm KB}$ と $\hat v^{\rm SL}$ はまったく別の演算子である．この違いが病的な形で現れることがある．

機構を見よう．参照状態 $\phi_{lm}$ と直交する任意の状態 $\psi$ を考えると，ステップ1の計算をなぞれば $\left\langle\delta v_l\phi_{lm}\middle|\psi\right\rangle$ は一般にゼロではないが，$\psi$ の作り方によってはほとんどゼロにできる．そのとき $\hat v^{\rm KB}$ の非局所項は $\psi$ にほとんど作用せず，$\psi$ は $v_{\rm loc}$ だけを感じる．もし $v_{\rm loc}$ がその $l$ チャネルにとって深すぎる引力なら，$\psi$ は参照エネルギーより低いところに束縛されてしまう．セミローカル形なら $\delta v_l$ の反発がこれを防いでいたのに，分離形ではその防壁が「参照状態方向」にしか立っていないのである．こうして現れる偽の束縛状態をゴースト状態(ghost state)と呼ぶ．

物理的意味:ゴースト状態が起きる条件

危険信号はKBエネルギーの符号である．$E_l^{\rm KB} = \left\langle\phi_l\right|\delta v_l\left|\phi_l\right\rangle$ が負(すなわち $\delta v_l$ が平均として引力)のとき，KB項 $\left|\delta v_l\phi_l\right\rangle\left\langle\delta v_l\phi_l\right|/E_l^{\rm KB}$ は負の係数を持つランク1演算子となり，変分的にエネルギーを下げる方向に働く．$\left|E_l^{\rm KB}\right|$ が小さいときは分母が小さいので係数の絶対値が大きくなり，いっそう危険である．$\delta v_l = v_l - v_{\rm loc}$ だから，これは局所ポテンシャル $v_{\rm loc}$ の選び方の問題である．$v_{\rm loc}$ として深いチャネル(たとえば $s$)を選ぶと他のチャネルの $\delta v_l$ が反発的になりゴーストは出にくく，浅いチャネルを選ぶと逆になる．Gonze らはこの条件を，セミローカル・ハミルトニアンの固有値と $\delta v_l$ の符号から系統的に判定する規準として定式化した(Gonze–Stumpf–Scheffler 1991)．

検査法は2つある．

分離形のまま原子問題を解き直す．$\hat T + v_{\rm loc} + \hat v^{\rm KB}_{\rm nl}$ を球対称原子に対して対角化し，各 $l$ の最低固有値が参照値 $\varepsilon_l$ と一致するか，それより低い状態がないかを見る．ゴーストがあれば $\varepsilon_l$ より低い節を持つ状態として現れる．
対数微分曲線を比較する．エネルギー $\varepsilon$ を掃引して $D_l(\varepsilon, r_d)$($r_d\gtrsim r_c$)を，全電子・セミローカル・分離形の3者について描く．$D_l$ は $\varepsilon$ の増加に対して単調に減少し($\partial D_l/\partial\varepsilon < 0$，式 \eqref{eq:15-normid})，$u_l(r_d)=0$ となるたびに $-\infty$ から $+\infty$ へ飛ぶ．この極の本数と位置が3者で一致しなければならない．ゴースト状態は「余分な極」として一目で見える(図15.7)．

治療法も同じ場所にある．$v_{\rm loc}$ の選択を変える(別の $l$ にする，あるいは全 $l$ から独立に平滑な局所ポテンシャルを作る)，$r_c$ を変える，そして1つの $l$ に対して射影子を複数本用意する．最後の処方は，ランク1では足りないなら階数を上げればよいという素直な発想であり，そのまま次節のVanderbilt型・MBK型の擬ポテンシャルにつながる．

図15.7 対数微分 $D_l(\varepsilon, r_d)$ によるゴースト状態の検出．(a) 正常な分離形では全電子(青・実線)と分離形(緑・破線)の曲線が，極の位置も含めて広いエネルギー範囲で重なる．(b) ゴースト状態があると，全電子には存在しない余分な極が現れ，その低エネルギー側で対数微分が全電子と食い違う．

15.7 ウルトラソフト擬ポテンシャルとノルム保存化

前節の末尾で「1つの $l$ に複数の射影子を持たせる」という方向が見えた．この方向を推し進めると，Vanderbilt (1990) のウルトラソフト擬ポテンシャル(ultrasoft pseudopotential, USPP)にたどり着く．ウルトラソフト法は，ノルム保存という枷を外すことで圧倒的な柔らかさを得るが，代わりに一般化固有値問題という構造変化を招く．本節ではその論理を追い，最後に「Vanderbiltの多エネルギー参照の利点だけを受け取り，ノルム保存に戻すことで実装を単純化する」というMorrison–Bylander–Kleinman (MBK, 1993) の考え方まで進む．

15.7.1 ノルム保存が課す下限

ノルム保存擬ポテンシャルにも限界がある．酸素の $2p$，窒素の $2p$，遷移金属の $3d$，希土類の $4f$ のように，その $l$ の最低状態でありもともと節を持たない価電子軌道を考えよう．この場合，全電子波動関数はすでに「これ以上滑らかにしようがない」形をしている．それでも波動関数が核近傍に鋭く局在しているのは，節のせいではなく単に軌道が小さいからである．

ここでノルム保存が効いてくる．条件3(ノルム保存)と条件1($r\ge r_c$ で一致)を同時に課すと，$r_c$ の内側に収める電荷量と，$r_c$ での値・傾きが両方固定される．すると擬波動関数は「決められた量の電荷を，決められた端点条件で，決められた球の中に詰める」ことになり，広げる自由度がない．実際，$r_c$ を大きくして無理に広げようとすると，$r_c$ の内側で波動関数がへこんだり盛り上がったりして，かえって高い波数成分が増えることさえある．酸素で $E_{\rm cut}\approx 70$–$80$ Ha，$3d$ 遷移金属で $E_{\rm cut}\approx 60$ Ha 程度が必要になるのはこのためである．

15.7.2 Vanderbiltの発想:ノルムを捨て，電荷を補償する

Vanderbiltの提案は思い切っている．ノルム保存条件を放棄する．そうすれば擬波動関数は $r_c$ を大きくとって好きなだけ滑らかにでき，$E_{\rm cut}$ は $E_{\rm cut}\approx 15$–$20$ Ha にまで落ちる．もちろん代償がある．

電子密度が足りなくなる．$\left|\phi^{\rm PS}\right|^2$ を積分しても $r_c$ 内の真の電荷にならない．
恒等式 \eqref{eq:15-normid} の保証が失われる．対数微分のエネルギー1次微分が自動的には一致しない．

2つ目は，参照エネルギーを1つの $l$ につき複数個($\varepsilon_{l1}, \varepsilon_{l2},\ldots$)とり，それぞれで対数微分の値を合わせることで直接解決する．1点で値と傾きを合わせる代わりに，複数点で値を合わせるのである(これは実は移植性をさらに広いエネルギー窓に広げる)．1つ目は，失った電荷を明示的に帳簿に書き戻すことで解決する．

定義:補償電荷(augmentation charge)

参照状態 $i=(n,l,m,\varepsilon_i)$ に対する全電子波動関数を $\psi_i$，擬波動関数を $\phi_i$ とするとき， $$Q_{ij}(\rr) \equiv \psi_i^*(\rr)\psi_j(\rr) - \phi_i^*(\rr)\phi_j(\rr)$$ を補償電荷密度という．条件1により $r \ge r_c$ で $\psi_i = \phi_i$ だから，$Q_{ij}(\rr)$ は半径 $r_c$ の球の内側にのみ台を持つ．その積分値を $$q_{ij} \equiv \int\dd^3 r\; Q_{ij}(\rr)$$ と書く．ノルム保存擬ポテンシャルとは，まさに $q_{ii}=0$(さらにMBKでは $q_{ij}=0$)が成り立つ場合にほかならない．

そして電子密度を，滑らかな部分と補償部分の和として定義し直す:

$$ \begin{equation} n(\rr) = \sum_{n}^{\rm occ}\left[\left|\phi_n(\rr)\right|^2 + \sum_{ij} Q_{ij}(\rr)\left\langle\phi_n\middle|\beta_i\right\rangle\left\langle\beta_j\middle|\phi_n\right\rangle\right] \label{eq:15-us-dens} \end{equation} $$

ここで $\left|\beta_i\right\rangle$ は次項で定義する射影子である．$\left\langle\beta_j|\phi_n\right\rangle$ は「結晶中のバンド $n$ の波動関数が，その原子サイトの参照状態 $j$ をどれだけ含んでいるか」を測る係数であり，その重みで補償電荷が原子サイトに置き戻される．

導出:重なり演算子と一般化固有値問題

ステップ1:電子数を数える．式 \eqref{eq:15-us-dens} を全空間で積分する:

$$ \begin{align} \int\dd^3 r\; n(\rr) &= \sum_n\left[\int\dd^3r\left|\phi_n\right|^2 + \sum_{ij}\left(\int\dd^3r\,Q_{ij}(\rr)\right)\left\langle\phi_n\middle|\beta_i\right\rangle\left\langle\beta_j\middle|\phi_n\right\rangle\right] \label{eq:15-us-n1}\\ &= \sum_n\left[\left\langle\phi_n\middle|\phi_n\right\rangle + \sum_{ij}q_{ij}\left\langle\phi_n\middle|\beta_i\right\rangle\left\langle\beta_j\middle|\phi_n\right\rangle\right] \label{eq:15-us-n2}\\ &= \sum_n\left\langle\phi_n\right|\underbrace{\left(1 + \sum_{ij}q_{ij}\left|\beta_i\right\rangle\left\langle\beta_j\right|\right)}_{\displaystyle \equiv\ \hat S}\left|\phi_n\right\rangle. \label{eq:15-us-n3} \end{align} $$

1行目は積分と和の順序を入れ替え，$\rr$ に依らない係数 $\left\langle\phi_n|\beta_i\right\rangle\left\langle\beta_j|\phi_n\right\rangle$ を積分の外に出した．2行目は $q_{ij}$ の定義．3行目では，括弧の中身を1つの演算子 $\hat S$ としてくくり出した(演算子の外積表示 $\left|\beta_i\right\rangle\left\langle\beta_j\right|$ を $\left\langle\phi_n\right|$ と $\left|\phi_n\right\rangle$ で挟むと確かに係数が再現される)．

ステップ2:規格化条件を読み取る．電子数が正しく $N_v$ になるためには，軌道の規格化条件は $\left\langle\phi_n|\phi_n\right\rangle=1$ ではなく

$$ \begin{equation} \left\langle\phi_n\right|\hat S\left|\phi_m\right\rangle = \delta_{nm}, \qquad \hat S = 1 + \sum_{ij}q_{ij}\left|\beta_i\right\rangle\left\langle\beta_j\right| \label{eq:15-us-S} \end{equation} $$

でなければならない．$\hat S$ を重なり演算子と呼ぶ．$q_{ij}=0$(ノルム保存)なら $\hat S=1$ に戻ることを確認しておこう．

ステップ3:変分原理を課す．全エネルギー $E[\{\phi_n\}]$ を，拘束条件 \eqref{eq:15-us-S} のもとで最小化する．第10章と同じくラグランジュ未定乗数 $\varepsilon_{nm}$ を導入し， $$\mathcal{L}[\{\phi\}] = E[\{\phi\}] - \sum_{nm}\varepsilon_{nm}\left(\left\langle\phi_n\right|\hat S\left|\phi_m\right\rangle - \delta_{nm}\right)$$ を $\left\langle\phi_n\right|$ について汎関数微分してゼロと置く(汎関数の記号に $\mathcal{L}$ を使うのは，本章で $\Omega$ を単位胞の体積 \eqref{eq:15-npw} に予約しているためである)．$\delta E/\delta\left\langle\phi_n\right| = \hat H\left|\phi_n\right\rangle$，$\delta\left(\left\langle\phi_n\right|\hat S\left|\phi_m\right\rangle\right)/\delta\left\langle\phi_n\right| = \hat S\left|\phi_m\right\rangle$ だから

$$ \hat H\left|\phi_n\right\rangle = \sum_m \varepsilon_{nm}\hat S\left|\phi_m\right\rangle. $$

$\varepsilon_{nm}$ はエルミート行列なので，それを対角化するユニタリ変換を軌道に施せば(第10章のコーン・シャム方程式の導出と同じ論法)，対角形

$$ \begin{equation} \hat H\left|\phi_n\right\rangle = \varepsilon_n\,\hat S\left|\phi_n\right\rangle \label{eq:15-us-gen} \end{equation} $$

が得られる．これは通常の固有値問題ではなく一般化固有値問題である．ここでの $\hat S$ は，第14章の非直交局在基底で現れた重なり行列とは起源が違う(基底の非直交性ではなく，擬変換そのものが持つ計量の変形である)ことに注意する．

ウルトラソフト法の実務上の代償をもう一つ挙げておく．$\hat H$ の非局所項の係数は

$$ \begin{equation} D_{ij} = B_{ij} + \int\dd^3r\; v_{\rm eff}(\rr)\,Q_{ij}(\rr) \label{eq:15-us-D} \end{equation} $$

という形になり，自己無撞着ポテンシャル $v_{\rm eff}$ を含む．したがって$D_{ij}$ は自己無撞着ループの中で毎回更新しなければならない．密度の式 \eqref{eq:15-us-dens} も，力の表式も，応力の表式も，すべて補償項の分だけ複雑になる．ウルトラソフト法の「速さ」は，この実装の複雑さと引き換えに得られている．

図15.8 ウルトラソフト擬ポテンシャルの考え方．左:ノルム保存を捨てると擬波動関数(赤・破線)を大きな $r_c$ まで極端に滑らかにでき，必要カットオフが下がる．右:失われた電荷を補償電荷 $Q_{ij}(\rr)$ として密度に足し戻す．その帰結として重なり演算子 $\hat S$ が現れ，方程式は一般化固有値問題 \eqref{eq:15-us-gen} になる．

数学ノート:PAW法との関係

Blöchl (1994) のprojector augmented-wave (PAW) 法は，滑らかな擬波動関数 $\left|\tilde\psi\right\rangle$ から真の全電子波動関数 $\left|\psi\right\rangle$ への線形変換 $$\left|\psi\right\rangle = \hat{\mathcal T}\left|\tilde\psi\right\rangle,\qquad \hat{\mathcal T} = 1 + \sum_i\left(\left|\psi_i\right\rangle - \left|\phi_i\right\rangle\right)\left\langle\beta_i\right|$$ を出発点にする．任意の演算子 $\hat A$ の期待値は $\left\langle\psi\right|\hat A\left|\psi\right\rangle = \left\langle\tilde\psi\right|\hat{\mathcal T}^\dagger\hat A\hat{\mathcal T}\left|\tilde\psi\right\rangle$ で計算され，$\hat A = 1$ とおけば $\hat{\mathcal T}^\dagger\hat{\mathcal T} = \hat S$ がまさにウルトラソフトの重なり演算子になる．すなわちウルトラソフト擬ポテンシャルはPAW変換から，全電子波動関数の情報を密度と全エネルギーの範囲に限って取り出した近似であると位置づけられる(Kresse–Joubert 1999)．PAWは節を持つ全電子波動関数を $\hat{\mathcal T}$ を通じて陽に保持するので，内殻近傍の量(超微細相互作用，電場勾配など)も扱える．

15.7.3 分離形の一般構成:$\chi$，$B$，$\beta$

ここで，Vanderbilt型の非局所演算子を一般的に構成し，そのエルミート性がノルム保存とどう結びつくかを見る．この構成はKB形($1$ チャネル1参照)をそのまま多参照に拡張したものであり，MBK法の土台でもある．

局所ポテンシャルとして，遮蔽された局所部分 $v^{\rm SL}(r) = v_{\rm loc}(r) + v_{\rm H}[\rho_v](r) + v_{\rm xc}[\rho_v+\tilde\rho_c](r)$ をとる．参照状態 $i$(エネルギー $\varepsilon_i$，擬波動関数 $\phi_i$)ごとに

$$ \begin{equation} \left|\chi_i\right\rangle \equiv \left(\varepsilon_i - \hat T - v^{\rm SL}\right)\left|\phi_i\right\rangle, \qquad B_{ij} \equiv \left\langle\phi_i\middle|\chi_j\right\rangle, \qquad \left|\beta_i\right\rangle \equiv \sum_j\left(B^{-1}\right)_{ji}\left|\chi_j\right\rangle \label{eq:15-chi-B-beta} \end{equation} $$

と定義し，非局所演算子を

$$ \begin{equation} \hat v^{\rm NL} = \sum_{ij} B_{ij}\left|\beta_i\right\rangle\left\langle\beta_j\right| \label{eq:15-vnl-def} \end{equation} $$

とする．$\left|\chi_i\right\rangle$ の意味を確認しよう．$\phi_i$ が真に $\left(\hat T + v^{\rm SL} + \hat v^{\rm NL}\right)\phi_i = \varepsilon_i\phi_i$ を満たしてほしいのだから，$\hat v^{\rm NL}\phi_i$ は $\left(\varepsilon_i - \hat T - v^{\rm SL}\right)\phi_i = \chi_i$ に等しくなければならない．$\chi_i$ は「非局所項が $\phi_i$ に対して出すべき答え」そのものである．また $r \ge r_c$ では $\phi_i = \psi_i$ かつ $v^{\rm SL}$ が真のポテンシャルに一致するので，$\chi_i(\rr) = 0$($r\ge r_c$)—$\chi_i$ も $\beta_i$ も球の内側だけの短距離関数である．

導出:$\hat v^{\rm NL}\left|\phi_k\right\rangle = \left|\chi_k\right\rangle$ の証明

$B$ が正則(逆行列を持つ)でエルミートだとする．まず射影子と参照状態の内積を計算する:

$$ \begin{align} \left\langle\beta_j\middle|\phi_k\right\rangle &= \sum_m \left(B^{-1}\right)_{mj}^{*}\left\langle\chi_m\middle|\phi_k\right\rangle \label{eq:15-bp1}\\ &= \sum_m \left(B^{-1}\right)_{jm}\,B_{km}^{*} \label{eq:15-bp2}\\ &= \sum_m B_{mk}\left(B^{-1}\right)_{jm} = \left(B^{-1}B\right)_{jk} = \delta_{jk}. \label{eq:15-bp3} \end{align} $$

1行目は $\left|\beta_j\right\rangle$ の定義のブラをとった(係数が複素共役になる)．2行目では $B^{-1}$ のエルミート性 $\left(B^{-1}\right)^*_{mj} = \left(B^{-1}\right)_{jm}$($B$ がエルミートなら $B^{-1}$ もエルミート)と，$\left\langle\chi_m|\phi_k\right\rangle = \left\langle\phi_k|\chi_m\right\rangle^* = B_{km}^*$ を使った．3行目では $B$ のエルミート性 $B_{km}^* = B_{mk}$ を使い，行列積の定義に戻した．

これを使って

$$ \begin{align} \hat v^{\rm NL}\left|\phi_k\right\rangle &= \sum_{ij}B_{ij}\left|\beta_i\right\rangle\left\langle\beta_j\middle|\phi_k\right\rangle = \sum_{ij}B_{ij}\,\delta_{jk}\left|\beta_i\right\rangle = \sum_i B_{ik}\left|\beta_i\right\rangle \label{eq:15-vnl-1}\\ &= \sum_i B_{ik}\sum_j\left(B^{-1}\right)_{ji}\left|\chi_j\right\rangle = \sum_j\left[\sum_i\left(B^{-1}\right)_{ji}B_{ik}\right]\left|\chi_j\right\rangle \label{eq:15-vnl-2}\\ &= \sum_j\delta_{jk}\left|\chi_j\right\rangle = \left|\chi_k\right\rangle. \label{eq:15-vnl-3} \end{align} $$

1行目は \eqref{eq:15-bp3} で $j$ 和を潰した．2行目は $\left|\beta_i\right\rangle$ の定義を代入し，$i$ 和と $j$ 和を入れ替えた．3行目は $\sum_i(B^{-1})_{ji}B_{ik} = (B^{-1}B)_{jk} = \delta_{jk}$．$\blacksquare$

したがって，参照状態に対しては

$$ \left(\hat T + v^{\rm SL} + \hat v^{\rm NL}\right)\left|\phi_k\right\rangle = \left(\hat T + v^{\rm SL}\right)\left|\phi_k\right\rangle + \left(\varepsilon_k - \hat T - v^{\rm SL}\right)\left|\phi_k\right\rangle = \varepsilon_k\left|\phi_k\right\rangle $$

が成り立つ．1つの $l$ に何個参照状態を置いても，そのすべてが厳密に固有値・固有関数として再現される．KB形 \eqref{eq:15-kb} は，$l$ ごとに参照が1個の場合($B$ が $1\times1$ で $B_{11} = \left\langle\phi\right|\delta v_l\left|\phi\right\rangle = E^{\rm KB}_l$，$\chi = \delta v_l\phi$)にほかならない．

15.7.4 $B$ のエルミート性とノルム保存の同値性

上の証明は $B$ のエルミート性を仮定していた．それはいつ成り立つのか．答えは驚くほどきれいで，15.4節で導いたWronskianの式 \eqref{eq:15-wronsk-rc} がそのまま効く．

定理:$B$ の非エルミート性は一般化ノルム差に等しい

$$ B_{ij} - B_{ji}^{*} = \left(\varepsilon_i - \varepsilon_j\right)Q_{ij}, \qquad Q_{ij} \equiv \int_0^{r_c}\!\left[u_i^{\rm AE}u_j^{\rm AE} - u_i^{\rm PS}u_j^{\rm PS}\right]\dd r $$ すなわち，一般化ノルム保存条件 $Q_{ij}=0$ と，$B$ のエルミート性は同値である(参照エネルギーが縮退していない限り)．

記号の注意:ここに現れた $Q_{ij}$(添字だけで引数を持たない量)は，15.7.2項で定義した補償電荷密度 $Q_{ij}(\rr)$ ではなく，その積分値 $q_{ij}=\int Q_{ij}(\rr)\dd^3r$ そのものである．実際，同じ $(l,m)$ に属する参照状態では角度積分が $1$ を与えるので $\int Q_{ij}(\rr)\dd^3r=\int_0^{r_c}\bigl[u_i^{\rm AE}u_j^{\rm AE}-u_i^{\rm PS}u_j^{\rm PS}\bigr]\dd r$ となり，両者は一致する．以下 $Q_{ij}=0$ と書いたら $q_{ij}=0$(重なり演算子が $\hat S=1$ に戻ること)と同じ意味である．

導出:定理の証明

同じ $l$ に属する2つの参照状態を考え，動径関数を実にとる($u_i \equiv u_i^{\rm PS}$)．積分は $\chi$ の台が $[0,r_c]$ であることから半径 $r_c$ の球に限られる．

ステップ1:$B_{ij}$ と $B_{ji}^*$ の差を書く．定義より

$$ \begin{align} B_{ij} &= \left\langle\phi_i\right|\left(\varepsilon_j - \hat T - v^{\rm SL}\right)\left|\phi_j\right\rangle, \qquad B_{ji}^{*} = \left\langle\phi_j\right|\left(\varepsilon_i - \hat T - v^{\rm SL}\right)\left|\phi_i\right\rangle^{*} \label{eq:15-B1}\\ B_{ij} - B_{ji}^{*} &= \left(\varepsilon_j - \varepsilon_i\right)\left\langle\phi_i\middle|\phi_j\right\rangle_{r_c} - \left[\left\langle\phi_i\right|\hat T\left|\phi_j\right\rangle_{r_c} - \left\langle\phi_j\right|\hat T\left|\phi_i\right\rangle_{r_c}\right]. \label{eq:15-B2} \end{align} $$

2行目では，$v^{\rm SL}$ が実の掛け算演算子なので $\left\langle\phi_i\right|v^{\rm SL}\left|\phi_j\right\rangle = \left\langle\phi_j\right|v^{\rm SL}\left|\phi_i\right\rangle^*$ となって完全に相殺したことを使った．$\left\langle\cdot\right\rangle_{r_c}$ は半径 $r_c$ の球内の積分を表す．

ステップ2:運動エネルギー項をWronskianに直す．動径表示では $\hat T \to -\frac{1}{2}\frac{\dd^2}{\dd r^2} + \frac{l(l+1)}{2r^2}$ で，遠心力項は共通なので差をとると消える．残るのは

$$ \begin{align} \left\langle\phi_i\right|\hat T\left|\phi_j\right\rangle_{r_c} - \left\langle\phi_j\right|\hat T\left|\phi_i\right\rangle_{r_c} &= -\frac{1}{2}\int_0^{r_c}\left(u_i u_j'' - u_j u_i''\right)\dd r \label{eq:15-B3}\\ &= -\frac{1}{2}\Big[W[u_i,u_j](r_c) - W[u_i,u_j](0)\Big] = -\frac{1}{2}W^{\rm PS}(r_c). \label{eq:15-B4} \end{align} $$

2行目では，15.4節の数学ノートで確認した $\frac{\dd}{\dd r}W = u_iu_j'' - u_i''u_j$ と微積分学の基本定理を使い，原点で $W(0)=0$(15.4節ステップ4と同じ理由:$u\propto r^{l+1}$ だから $W\propto r^{2l+1}\to0$)を使った．

ステップ3:境界のWronskianを全電子量で書き換える．擬波動関数は条件1・4により $r_c$ で全電子波動関数と値も1階微分も一致する．Wronskianは値と1階微分だけで作られるから

$$ W^{\rm PS}(r_c) = u_i^{\rm PS}(r_c)u_j^{\rm PS\prime}(r_c) - u_i^{\rm PS\prime}(r_c)u_j^{\rm PS}(r_c) = W^{\rm AE}(r_c). $$

そして全電子動径関数は同じ全電子ポテンシャルの下でエネルギーだけが違う2解だから，式 \eqref{eq:15-wronsk-rc} がそのまま使えて

$$ \begin{equation} W^{\rm AE}(r_c) = 2\left(\varepsilon_i - \varepsilon_j\right)\int_0^{r_c}u_i^{\rm AE}u_j^{\rm AE}\,\dd r. \label{eq:15-B5} \end{equation} $$

(擬波動関数の方はこの関係を満たさないことに注意せよ．擬波動関数が満たすのは非局所項を含む方程式であり，単純な2階常微分方程式ではない．だからこそ左辺と右辺に差が生じ，それが $Q_{ij}$ として現れる．)

ステップ4:まとめる．\eqref{eq:15-B4} と \eqref{eq:15-B5} を \eqref{eq:15-B2} に代入する:

$$ \begin{align} B_{ij}-B_{ji}^{*} &= \left(\varepsilon_j-\varepsilon_i\right)\int_0^{r_c}u_i^{\rm PS}u_j^{\rm PS}\dd r + \frac{1}{2}\cdot 2\left(\varepsilon_i-\varepsilon_j\right)\int_0^{r_c}u_i^{\rm AE}u_j^{\rm AE}\dd r \label{eq:15-B6}\\ &= \left(\varepsilon_i-\varepsilon_j\right)\left[\int_0^{r_c}u_i^{\rm AE}u_j^{\rm AE}\dd r - \int_0^{r_c}u_i^{\rm PS}u_j^{\rm PS}\dd r\right] \label{eq:15-B7}\\ &= \left(\varepsilon_i-\varepsilon_j\right)Q_{ij}. \label{eq:15-B8} \end{align} $$

1行目は代入(第2項の符号は $-\left(-\frac12 W\right) = +\frac12 W$ から来る)．2行目は $\left(\varepsilon_i-\varepsilon_j\right)$ でくくった．$\blacksquare$

物理的意味:$Q_{ij}$ が3つの顔を持つ

同じ量 $Q_{ij}$ が3か所に現れたことに注意しよう．(i) $i=j$ のときの $Q_{ii}$ はノルム保存条件そのもの(15.4節)である．(ii) $q_{ij}=\int Q_{ij}(\rr)\dd^3 r$ は重なり演算子 $\hat S$ を1からずらす量(15.7.2項)である．(iii) そしていま，それは$B$ 行列の非エルミート性でもあった．ウルトラソフト法は $Q_{ij}\ne0$ を積極的に許して滑らかさを買い，その代償として $\hat S \ne 1$ という一般化固有値問題と，エルミート性を回復するための追加の対称化を引き受ける．逆に $Q_{ij}=0$ を課せば，$\hat S=1$，$B$ はエルミート，非局所項は対角化可能—すべてが同時に単純になる．

15.7.5 MBK法:多参照のままノルム保存に戻す

Morrison, Bylander, Kleinman (1993) の着眼はここにある．Vanderbiltの構成のうち「1つの $l$ に複数の参照エネルギーを持てる」という利点だけを取り，$Q_{ij}=0$ を課してノルム保存に戻す．すると:

$\hat S = 1$，方程式は通常の固有値問題 $\hat H\phi = \varepsilon\phi$ に戻る．
$B$ はエルミートだから，ユニタリ行列 $U$ で対角化できる．$B = U\Lambda U^\dagger$ と書き，$\left|\tilde\beta_i\right\rangle = \sum_j U_{ji}^*\left|\beta_j\right\rangle$ と変換すれば $$\hat v^{\rm NL} = \sum_{ij}B_{ij}\left|\beta_i\right\rangle\left\langle\beta_j\right| = \sum_{i}\lambda_i\left|\tilde\beta_i\right\rangle\left\langle\tilde\beta_i\right|$$ という対角(多射影子KB)形になる．これはKB形 \eqref{eq:15-kb} の自然な多本化であり，平面波でも局在基底でも実装は既存のKBルーチンのわずかな拡張で済む．
密度は $n=\sum_n|\phi_n|^2$ のまま，力・応力の表式もノルム保存型と同一．$D_{ij}$ の自己無撞着更新も不要．

残る仕事は「$Q_{ij}=0$ をどうやって満たすか」である．ここで自由度が要る．1つの参照だけならTM法がすでに $Q_{ii}=0$ を実現しているが，参照を2個3個と増やすと $Q_{ij}$($i\ne j$)は一般にゼロにならない．そこでMBK法は，TM擬波動関数に補正関数を足す:

$$ \begin{equation} u_i^{\rm PS}(r) = u_i^{\rm TM}(r) + f_i(r), \qquad f_i(r) = \sum_{k=1}^{N_b}c_{ik}\;r\,j_l\!\left(q_{lk}r\right), \qquad q_{lk} = \frac{u_{lk}}{r_c} \label{eq:15-mbk-f} \end{equation} $$

ここで $u_{lk}$ は球ベッセル関数 $j_l$ の第 $k$ 零点である($j_l(u_{lk})=0$)．この基底の選び方は絶妙で，次の性質が自動的についてくる．

導出:球ベッセル補正関数の境界条件が自動的に満たされること

$f_k(r) \equiv r\,j_l(q_kr)$ とおく($q_k \equiv q_{lk}$)．球ベッセル関数の定義から，$f_k$ はポテンシャルのない動径方程式(エネルギー $q_k^2/2$)の原点正則解である:

$$ \begin{equation} f_k''(r) = g_k(r)\,f_k(r), \qquad g_k(r) \equiv \frac{l(l+1)}{r^2} - q_k^2. \label{eq:15-mbk-free} \end{equation} $$

(これは式 \eqref{eq:15-radial} で $V=0$，$\varepsilon = q_k^2/2$ とした形にほかならない．)まず $q_k = u_{lk}/r_c$ という選び方から

$$ f_k(r_c) = r_c\,j_l\!\left(\frac{u_{lk}}{r_c}r_c\right) = r_c\,j_l(u_{lk}) = 0 $$

がすべての $k$ について成り立つ．すなわち補正関数は $r_c$ で必ずゼロになり，$r\ge r_c$ で全電子と一致するという条件1を壊さない．次に \eqref{eq:15-mbk-free} を $r_c$ で評価すると

$$ f_k''(r_c) = g_k(r_c)f_k(r_c) = 0 $$

となり，2階微分の一致(条件4の $k=2$)も自動的に保たれる．さらに \eqref{eq:15-mbk-free} を1回，2回微分すると

$$ \begin{align} f_k''' &= g_k' f_k + g_k f_k' \;\;\longrightarrow\;\; f_k'''(r_c) = g_k(r_c)\,f_k'(r_c), \label{eq:15-mbk-d3}\\ f_k'''' &= g_k'' f_k + 2g_k' f_k' + g_k f_k'' \;\;\longrightarrow\;\; f_k''''(r_c) = 2g_k'(r_c)\,f_k'(r_c) \label{eq:15-mbk-d4} \end{align} $$

となる(いずれも $f_k(r_c)=f_k''(r_c)=0$ を使った)．ここで $g_k'(r) = -2l(l+1)/r^3$ は $k$ に依らないことに注意．したがって和 $f_i = \sum_k c_{ik}f_k$ については

$$ f_i''''(r_c) = 2g'(r_c)\sum_k c_{ik}f_k'(r_c) = 2g'(r_c)\,f_i'(r_c) $$

となり，$f_i'(r_c)=0$ さえ課せば4階微分の一致も自動的に従う．一方 $f_i'''(r_c) = \sum_k c_{ik}\,g_k(r_c)f_k'(r_c)$ は，重み $g_k(r_c) = l(l+1)/r_c^2 - q_{lk}^2$ が $k$ ごとに異なるため，$f_i'(r_c)=0$ からは従わない．

結局，独立に課すべき条件は $$f_i'(r_c) = 0,\qquad f_i'''(r_c) = 0,\qquad Q_{ij}=0\;(j=1,\ldots,N_{\rm ref})$$ の $N_{\rm ref}+2$ 個であり，係数 $c_{ik}$ の個数は $N_b = N_{\rm ref}+2$ でちょうど足りる．

$Q_{ij}=0$ の条件は $c_{ik}$ について2次(積 $u_i^{\rm PS}u_j^{\rm PS}$ が入るため)なので，実際には $c$ を小さいところから反復して解く．係数が決まれば $\chi_i$ も閉じた形で書ける:

$$ \begin{equation} \chi_i = \left(\varepsilon_i - \hat T - v^{\rm SL}\right)\left(u_i^{\rm TM}+f_i\right) = \underbrace{v_i^{\rm TM}u_i^{\rm TM} - v^{\rm SL}u_i^{\rm PS} + \varepsilon_i f_i}_{\text{既知量}} \;-\;\frac{1}{2}\sum_k c_{ik}\,q_{lk}^2\;r\,j_l(q_{lk}r) \label{eq:15-mbk-chi} \end{equation} $$

(第1項は $u_i^{\rm TM}$ がTMポテンシャル $v_i^{\rm TM}$ の固有関数であることを使い，最後の項は $f_k$ が自由粒子解であること \eqref{eq:15-mbk-free} から $\left(-\frac12\frac{\dd^2}{\dd r^2}+\frac{l(l+1)}{2r^2}\right)f_k = \frac{q_k^2}{2}f_k$ を使った．)数値微分を一切使わずに $\chi_i$ が得られるので，精度も安定する．

例:柔らかさと精度のトレードオフ

3種の擬ポテンシャルで酸素を扱ったときの典型的な必要カットオフを並べると，ノルム保存TM型で $E_{\rm cut}\approx 70$ Ha，ウルトラソフト型で $E_{\rm cut}\approx 15$ Ha，PAWで $E_{\rm cut}\approx 30$ Ha 程度である．式 \eqref{eq:15-npw} により基底数は $E_{\rm cut}^{3/2}$ に比例するから，TM型からウルトラソフト型へ移ると基底数は $(70/15)^{3/2}\approx 10$ 倍減る．対角化コストが基底数の3乗なら $10^3$ 倍—これがウルトラソフト法が広く使われる理由である．一方，精度の観点では，多参照ノルム保存型(MBK)は全電子法(FLAPW)との差(15.9節で定義する $\Delta$ ゲージ)で $2$ meV/atom 程度を達成し，全電子コード同士のばらつき($0.3$–$1$ meV/atom)に迫る．柔らかさと単純さと精度のどれを重視するかで選択が変わる，というのが現状である．

15.8 相対論的効果

ここまで，原子の全電子計算は非相対論的なシュレーディンガー方程式で行うものとしてきた．軽い元素ではそれでよい．しかし周期表を下へ降りていくと，この仮定は静かに，しかし決定的に破綻する．本節では，破綻がいつ起きるかを見積もり，Dirac方程式からスカラー相対論近似とスピン軌道項がどう分離するかを導き，$j$ 依存擬ポテンシャルの構成に結びつける．

15.8.1 いつ相対論が必要になるか

まず速度を見積もる．水素様原子の $1s$ 電子について，ビリアル定理(第2章)より運動エネルギーの期待値は全エネルギーの符号を変えたものに等しく，$\left\langle T\right\rangle = Z^2/2$(原子単位)である．非相対論的に $\left\langle T\right\rangle = \left\langle v^2\right\rangle/2$ とすれば

$$ \begin{equation} \sqrt{\left\langle v^2\right\rangle} = Z \quad\text{(原子単位)}, \qquad \frac{v}{c} = \frac{Z}{c} = Z\alpha \simeq \frac{Z}{137} \label{eq:15-vc} \end{equation} $$

となる．ここで $\alpha = 1/c \simeq 1/137.036$ は微細構造定数で，原子単位系では光速そのものの逆数である(15章冒頭)．具体的な数値を入れると

炭素($Z=6$):$v/c \approx 0.044$，相対論的補正は $O\!\left((v/c)^2\right)\approx 0.2\%$．無視できる．
銅($Z=29$):$v/c\approx0.21$，$(v/c)^2 \approx 4.5\%$．内殻準位に効くが価電子には間接的．
水銀($Z=80$):$v/c\approx0.58$．Lorentz因子は $\gamma = 1/\sqrt{1-0.58^2}\approx 1.23$ で，電子の実効質量が23%増える．
金($Z=79$)，鉛($Z=82$)，ビスマス($Z=83$)も同程度．

物理的意味:なぜ内殻の相対論効果が化学に効くのか

速度が大きいのは内殻電子であって価電子ではない．それなのに相対論効果が化学的性質を変えるのは，次の連鎖による．(1) 質量増加により $s$ 軌道(と $p_{1/2}$ 軌道)の実効Bohr半径 $a = 1/(m_{\rm eff}Z)$ が縮む．(2) 内殻 $s$ が縮むと，それと直交しなければならない価電子の $s$ 軌道も縮む．(3) 縮んだ $s$ 殻は核電荷をよく遮蔽するので，$d$ と $f$ 軌道は逆に広がって不安定化する．この「$s$ 収縮・$d$ 膨張」が，金の $5d\to6s$ 遷移エネルギーを可視域(約2.4 eV，青)まで下げて金色を生み，水銀の $6s$ を閉殻的に安定化して常温で液体にし，鉛蓄電池の起電力の大部分を担う．相対論は「重元素の化学そのもの」なのである．

実務上大切なのは，この効果を擬ポテンシャルの生成段階だけで取り込めばよいという点である．相対論的な原子計算(Dirac方程式)で全電子波動関数を作り，そこから擬ポテンシャルを構成すれば，固体の計算は非相対論的なシュレーディンガー方程式のままでよい．相対論の情報はすべて $v_l^{\rm ion}$ の形の中に凍結されている．

15.8.2 Dirac方程式から動径方程式へ

球対称ポテンシャル $V(r)$ 中の電子に対するDirac方程式は，角度部分を球スピノル(spin-spherical harmonics)$\mathcal{Y}_{\kappa m}$ で分離できる．分離後に残るのは，大成分 $G(r)$ と小成分 $F(r)$ に対する連立1階方程式である(静止エネルギーを除いたエネルギーを $\varepsilon$ とする):

$$ \begin{align} \frac{\dd G}{\dd r} &= -\frac{\kappa}{r}G + 2c\,M(r)\,F, \label{eq:15-dirac1}\\ \frac{\dd F}{\dd r} &= \frac{\kappa}{r}F - \frac{\varepsilon - V(r)}{c}\,G, \qquad M(r) \equiv 1 + \frac{\varepsilon - V(r)}{2c^2}. \label{eq:15-dirac2} \end{align} $$

$M(r)$ は相対論的な実効質量である($c\to\infty$ で $M\to1$，すなわち電子の静止質量)．深いポテンシャル($V\ll0$)の中では $M$ が1より大きくなる—これが「質量増加」の正体である．$\kappa$ は角運動量とスピンの結合を表す整数で

$$ \begin{equation} \kappa = \begin{cases} -(l+1), & j = l+\tfrac{1}{2},\\ +\,l, & j = l-\tfrac{1}{2}, \end{cases} \qquad \text{縮重度}\; 2j+1 = \begin{cases} 2l+2, & j = l+\tfrac12,\\ 2l, & j = l-\tfrac12. \end{cases} \label{eq:15-kappa} \end{equation} $$

と定義される．以下では小成分 $F$ を消去して，$G$ だけの2階方程式に直す．

導出:Dirac方程式から $G$ の2階方程式へ

ステップ1:$F$ を $G$ で表す．式 \eqref{eq:15-dirac1} を $F$ について解く:

$$ \begin{equation} F = \frac{1}{2cM}\left(G' + \frac{\kappa}{r}G\right) \equiv \frac{Y}{2cM}, \qquad Y \equiv G' + \frac{\kappa}{r}G. \label{eq:15-dr1} \end{equation} $$

ステップ2:これを \eqref{eq:15-dirac2} に代入する．

$$ \begin{align} \frac{\dd}{\dd r}\left(\frac{Y}{2cM}\right) &= \frac{\kappa}{r}\cdot\frac{Y}{2cM} - \frac{\varepsilon-V}{c}G \label{eq:15-dr2}\\ \frac{Y'}{2cM} - \frac{Y M'}{2cM^2} &= \frac{\kappa Y}{2cMr} - \frac{\varepsilon-V}{c}G. \label{eq:15-dr3} \end{align} $$

2行目は左辺に商の微分則を適用しただけである．両辺に $2cM$ を掛けると

$$ \begin{equation} Y' - \frac{M'}{M}Y = \frac{\kappa}{r}Y - 2M\left(\varepsilon - V\right)G. \label{eq:15-dr4} \end{equation} $$

ステップ3:$Y$ を $G$ で書き戻す．$Y = G' + \frac{\kappa}{r}G$ だから $Y' = G'' + \frac{\kappa}{r}G' - \frac{\kappa}{r^2}G$ である(第2項に商の微分則)．代入して

$$ \begin{align} G'' + \frac{\kappa}{r}G' - \frac{\kappa}{r^2}G - \frac{M'}{M}\left(G'+\frac{\kappa}{r}G\right) &= \frac{\kappa}{r}\left(G'+\frac{\kappa}{r}G\right) - 2M(\varepsilon-V)G \label{eq:15-dr5}\\ G'' - \frac{\kappa}{r^2}G - \frac{M'}{M}\left(G'+\frac{\kappa}{r}G\right) &= \frac{\kappa^2}{r^2}G - 2M(\varepsilon-V)G. \label{eq:15-dr6} \end{align} $$

2行目では両辺にある $\frac{\kappa}{r}G'$ を消去した．整理して

$$ \begin{equation} G'' = \frac{\kappa(\kappa+1)}{r^2}G + \frac{M'}{M}G' + \frac{M'}{M}\frac{\kappa}{r}G - 2M(\varepsilon-V)G. \label{eq:15-dr7} \end{equation} $$

ステップ4:$\kappa(\kappa+1)$ を評価する．ここが要点である．式 \eqref{eq:15-kappa} の2つの場合を代入すると

$$ \begin{align} j = l+\tfrac12\;(\kappa = -(l+1)):&\quad \kappa(\kappa+1) = -(l+1)\left[-(l+1)+1\right] = -(l+1)(-l) = l(l+1), \label{eq:15-dr8}\\ j = l-\tfrac12\;(\kappa = l):&\quad \kappa(\kappa+1) = l(l+1). \label{eq:15-dr9} \end{align} $$

どちらの場合も同じ $l(l+1)$ になる．つまり遠心力項には $j$ 依存性がまったく入らない．$\kappa$ 依存性，すなわちスピン軌道の効果は，\eqref{eq:15-dr7} の第3項 $\frac{M'}{M}\frac{\kappa}{r}G$ ただ一つに集約される．

ステップ5:標準形に整える．\eqref{eq:15-dr7} の両辺に $-\frac{1}{2M}$ を掛け，$\varepsilon G$ を右辺に集めると

$$ \begin{equation} \left[-\frac{1}{2M}\frac{\dd^2}{\dd r^2} + \frac{l(l+1)}{2Mr^2} + V(r) - \frac{V'(r)}{4M^2c^2}\frac{\dd}{\dd r} \;\underbrace{-\;\frac{V'(r)}{4M^2c^2}\frac{\kappa}{r}}_{\text{スピン軌道}}\right]G = \varepsilon\,G \label{eq:15-koelling} \end{equation} $$

を得る．ここで $M' = -V'/(2c^2)$(定義 \eqref{eq:15-dirac2} を $r$ で微分)を使い，$\frac{M'}{2M^2} = -\frac{V'}{4M^2c^2}$ と書き換えた．式 \eqref{eq:15-koelling} はKoelling–Harmon (1977) の相対論的動径方程式である．

15.8.3 スカラー相対論近似とスピン軌道分離

式 \eqref{eq:15-koelling} で $\kappa$ を含む最後の項だけがスピンに依存する．そこで，この項を $j = l\pm\frac12$ について縮重度の重みで平均してしまえば，スピンを一切含まない方程式が得られる．これがスカラー相対論近似である．平均すべき $\kappa$ の値は

$$ \begin{align} \kappa_{\rm av} &= \frac{(2l)\cdot l + (2l+2)\cdot\left[-(l+1)\right]}{2l + (2l+2)} \label{eq:15-kav1}\\ &= \frac{2l^2 - 2(l+1)^2}{4l+2} = \frac{2l^2 - 2l^2 - 4l - 2}{4l+2} = \frac{-(4l+2)}{4l+2} = -1 \label{eq:15-kav2} \end{align} $$

である．1行目は $j=l-\frac12$ の縮重度 $2l$ と $\kappa=l$，$j=l+\frac12$ の縮重度 $2l+2$ と $\kappa=-(l+1)$ を重み付き平均したもの．2行目は分子を展開して整理しただけである．$l$ に依らず $\kappa_{\rm av}=-1$ になるという事実は美しい．

$$ \begin{equation} \left[-\frac{1}{2M}\frac{\dd^2}{\dd r^2} + \frac{l(l+1)}{2Mr^2} + V(r) - \frac{V'(r)}{4M^2c^2}\left(\frac{\dd}{\dd r} - \frac{1}{r}\right)\right]G = \varepsilon\,G \label{eq:15-scalar-rel} \end{equation} $$

これがスカラー相対論方程式である($\kappa\to\kappa_{\rm av}=-1$ を代入すると $-\frac{V'}{4M^2c^2}\frac{\kappa}{r} \to +\frac{V'}{4M^2c^2}\frac{1}{r}$ となり，上のように括弧にまとめられる)．質量増加($M$)と，それに伴う軌道の収縮・膨張はすべて含まれ，スピン軌道分裂だけが落とされている．多くの物質の格子定数・凝集エネルギー・バンド幅は，これで十分な精度で得られる．

数学ノート:非相対論的極限—質量速度項とDarwin項

式 \eqref{eq:15-koelling} を $1/c^2$ の1次まで展開すると，教科書でおなじみの3項が現れる．3次元の形で書けば

$$ \hat H^{\rm rel} = \hat T + V \underbrace{-\frac{\hat p^4}{8c^2}}_{\text{質量速度項}} +\underbrace{\frac{1}{8c^2}\nabla^2 V}_{\text{Darwin項}} +\underbrace{\frac{1}{2c^2}\frac{1}{r}\frac{\dd V}{\dd r}\,\bm{L}\cdot\bm{S}}_{\text{スピン軌道項}} $$

質量速度項は相対論的運動エネルギーの展開から出る．静止エネルギーを引いた運動エネルギーは $$E_{\rm kin} = \sqrt{p^2c^2 + c^4} - c^2 = c^2\sqrt{1+\frac{p^2}{c^2}} - c^2$$ であり，$\sqrt{1+x} = 1 + \frac{x}{2} - \frac{x^2}{8} + O(x^3)$ を $x = p^2/c^2$ に適用すると $$E_{\rm kin} = c^2\left[1 + \frac{p^2}{2c^2} - \frac{p^4}{8c^4} + \cdots\right] - c^2 = \frac{p^2}{2} - \frac{p^4}{8c^2} + \cdots$$ となる．第2項が質量速度項で，常に負(運動エネルギーを下げる=束縛を強める)．$p$ が大きい核近傍で効くので，$s$ 軌道に最も強く働く．

Darwin項は，電子が正負エネルギー状態の干渉によりCompton波長 $\lambda_C = 1/c$ 程度の範囲にぼやける効果(Zitterbewegung)から生じる．ポテンシャルをその範囲で平均すると $\left\langle V(\rr+\bm\delta)\right\rangle \simeq V(\rr) + \frac{1}{6}\left\langle\delta^2\right\rangle\nabla^2V$ となり，$\left\langle\delta^2\right\rangle\sim 1/c^2$ から $\nabla^2V/(8c^2)$ が出る．裸のクーロンポテンシャル $V=-Z/r$ では $\nabla^2V = 4\pi Z\delta(\rr)$ なので，Darwin項は原点で有限の振幅を持つ $s$ 軌道($l=0$)にしか効かない．

スピン軌道項が式 \eqref{eq:15-koelling} の $\kappa$ 項に対応する．実際，$\bm L\cdot\bm S$ の固有値は次項で見るように $\kappa$ で書けて，$-\frac{V'}{4M^2c^2}\frac{\kappa}{r}$ が $\frac{V'}{2c^2 r}\bm L\cdot\bm S$ に対応する($M\to1$ の極限で)．

15.8.4 $j$ 依存擬ポテンシャルとその分解

相対論的擬ポテンシャルの作り方は原理的には簡単で，$l$ ごとではなく $(l,j)$ ごとに15.5節の手続きを実行するだけである．すなわちDirac方程式(または \eqref{eq:15-koelling})を $\kappa = -(l+1)$ と $\kappa = l$ の両方について解き，それぞれの大成分 $G_{lj}$ を全電子波動関数としてTM法を適用し，$v_{l,l+1/2}$ と $v_{l,l-1/2}$ という2本のポテンシャルを得る($l=0$ では $j=1/2$ しかないので1本)．

実際の固体計算では，この2本をそのまま使うのではなく，スピンに依らない部分とスピン軌道部分に分けるのが便利である．$\bm J = \bm L + \bm S$ から

$$ \begin{align} J^2 &= (\bm L + \bm S)^2 = L^2 + S^2 + 2\bm L\cdot\bm S \label{eq:15-ls1}\\ \bm L\cdot\bm S &= \frac{1}{2}\left[J^2 - L^2 - S^2\right] \;\longrightarrow\; \frac{1}{2}\left[j(j+1) - l(l+1) - \frac{3}{4}\right] \label{eq:15-ls2} \end{align} $$

(2行目で $s=1/2$ より $S^2 \to s(s+1) = 3/4$)．これに $j = l\pm\frac12$ を代入する．$j=l+\frac12$ では $j(j+1) = \left(l+\frac12\right)\left(l+\frac32\right) = l^2+2l+\frac34$ だから

$$ \bm L\cdot\bm S \to \frac{1}{2}\left[l^2+2l+\frac34 - l^2 - l - \frac34\right] = \frac{l}{2}, $$

$j=l-\frac12$ では $j(j+1) = \left(l-\frac12\right)\left(l+\frac12\right) = l^2-\frac14$ だから

$$ \bm L\cdot\bm S \to \frac{1}{2}\left[l^2-\frac14 - l^2 - l - \frac34\right] = -\frac{l+1}{2}. $$

そこで，スカラー相対論部分 $v_l^{\rm SR}$ とスピン軌道部分 $v_l^{\rm SO}$ を

$$ \begin{equation} v_l^{\rm SR}(r) = \frac{l\,v_{l,l-1/2}(r) + (l+1)\,v_{l,l+1/2}(r)}{2l+1}, \qquad v_l^{\rm SO}(r) = \frac{2\left[v_{l,l+1/2}(r) - v_{l,l-1/2}(r)\right]}{2l+1} \label{eq:15-sr-so} \end{equation} $$

と定義すると，$v_l^{\rm SR} + v_l^{\rm SO}\,\bm L\cdot\bm S$ が元の $j$ 依存ポテンシャルを正しく再現する．確かめてみよう．$j = l+\frac12$ のとき $\bm L\cdot\bm S \to l/2$ だから

$$ \begin{align} v_l^{\rm SR} + \frac{l}{2}v_l^{\rm SO} &= \frac{l\,v_{-} + (l+1)v_{+}}{2l+1} + \frac{l}{2}\cdot\frac{2\left(v_{+}-v_{-}\right)}{2l+1} \label{eq:15-check1}\\ &= \frac{l\,v_{-} + (l+1)v_{+} + l\,v_{+} - l\,v_{-}}{2l+1} = \frac{(2l+1)v_{+}}{2l+1} = v_{+} = v_{l,l+1/2}. \label{eq:15-check2} \end{align} $$

($v_\pm \equiv v_{l,l\pm1/2}$ と略記した．1行目は定義の代入，2行目は分子を展開して $l\,v_-$ が相殺したところ．)同様に $j=l-\frac12$ では $\bm L\cdot\bm S\to -(l+1)/2$ で

$$ v_l^{\rm SR} - \frac{l+1}{2}v_l^{\rm SO} = \frac{l\,v_{-}+(l+1)v_{+} - (l+1)v_{+} + (l+1)v_{-}}{2l+1} = \frac{(2l+1)v_{-}}{2l+1} = v_{l,l-1/2} $$

となり，両方とも正しく再現される．$v_l^{\rm SR}$ の重み $l:(l+1)$ は縮重度 $2l:(2l+2)$ の比にほかならず，$\kappa_{\rm av}=-1$ を導いたのと同じ平均である．したがって最終的な擬ポテンシャル演算子は

$$ \begin{equation} \hat v^{\rm ps} = v_{\rm loc}(r) + \sum_{l}\left[\delta v_l^{\rm SR}(r) + \delta v_l^{\rm SO}(r)\,\bm L\cdot\bm S\right]\hat P_l, \qquad \hat P_l = \sum_{m}\left|Y_{lm}\right\rangle\left\langle Y_{lm}\right| \label{eq:15-vps-rel} \end{equation} $$

という形になる(実装ではKB/MBK分離形にしてから使う)．スピン軌道項を落とせばスカラー相対論計算，残せばスピン軌道相互作用込みの計算になる．同じ擬ポテンシャルファイルで両方が扱えるのは，この分解の実務的な利点である．

例:半導体の価電子帯頂上のスピン軌道分裂

ダイヤモンド構造・閃亜鉛鉱構造の半導体では，価電子帯の頂上($\Gamma$ 点)は主に陰イオン $p$ 軌道からなる3重縮退状態である．スピン軌道相互作用を入れると，これが $j=3/2$ の4重項($\Gamma_8$)と $j=1/2$ の2重項($\Gamma_7$)に分裂する．式 \eqref{eq:15-ls2} から分裂幅は $$\Delta_{\rm SO} = \left(\frac{l}{2} - \left(-\frac{l+1}{2}\right)\right)\left\langle v_l^{\rm SO}\right\rangle = \frac{2l+1}{2}\left\langle v_l^{\rm SO}\right\rangle \quad(l=1\;\text{なら}\;\tfrac32\left\langle v_1^{\rm SO}\right\rangle)$$ である．実測値を並べると Si:$0.044$ eV，Ge:$0.29$ eV，$\alpha$-Sn:$0.80$ eV，GaAs:$0.34$ eV，InSb:$0.81$ eV となり，構成元素が重くなるにつれて急激に増大する(水素様の見積りでは $\Delta_{\rm SO}\propto Z^4/n^3$)．シリコンなら室温の熱エネルギー($0.026$ eV)と同程度でしばしば無視されるが，ゲルマニウム以降ではホール有効質量やバンドギャップの計算に決定的である．ビスマス系のトポロジカル絶縁体に至っては，スピン軌道相互作用そのものがバンド反転を引き起こして物質の性格を決めている．

図15.9 相対論的擬ポテンシャルとスピン軌道分裂．(a) $j=l\pm\frac12$ ごとに構成した擬ポテンシャルの差がスピン軌道項 $v_l^{\rm SO}$，縮重度平均がスカラー相対論項 $v_l^{\rm SR}$ を与える(式 \eqref{eq:15-sr-so})．(b) その結果，半導体の価電子帯頂上の3重縮退が $j=3/2$ と $j=1/2$ に分裂する．

15.9 検証の実際

擬ポテンシャルは近似である．しかもパラメータ($r_c$，局所チャネルの選択，参照配置，価電子に含める殻)を人間が選ぶ近似である．したがって作ったら必ず検証する．この節では，実務で行われる検証の階層を，原子レベル(対数微分)から固体レベル(バンド構造・状態方程式)まで順に述べ，最後に初学者が陥りやすい失敗を挙げる．

15.9.1 第一の関門:対数微分曲線

最初に，そして最も安価に行える検証が対数微分曲線の比較である．手続きは次のとおり．

診断半径 $r_d$ を選ぶ($r_d \ge \max_l r_{cl}$，通常は最大の $r_c$ より少し外)．
エネルギー $\varepsilon$ を，価電子帯の底から非占有状態の上まで(たとえば $-3$ Ha から $+3$ Ha まで)細かく掃引する．
各 $\varepsilon$ で，(i) 全電子(相対論的)動径方程式，(ii) セミローカル擬ポテンシャル，(iii) 実際に使う分離形(KB/MBK)，の3つについて原点正則解を数値積分し，$D_l(\varepsilon, r_d) = u_l'(r_d)/u_l(r_d)$ を求める．
$l=0,1,2,3$(相対論的なら $j$ ごと)について3本の曲線を重ねて描く．

物理的意味:対数微分曲線の読み方

15.3節で示したとおり，$r_d$ の外の世界にとって原子の中身の情報は $D_l(\varepsilon, r_d)$ に尽きる．したがってこの1枚のグラフは，移植性の直接の測定である．読むべき点は3つ．

(1) 一致の幅．3本が重なるエネルギー範囲が，その擬ポテンシャルの「有効な窓」である．ノルム保存型は参照エネルギー $\varepsilon_l$ のまわりで1次まで一致する(式 \eqref{eq:15-taylor})ので，$\varepsilon_l$ の近傍で必ず重なる．窓の広さは $r_c$ に強く依存し，$r_c$ を大きくすると窓が狭くなる．実際の結晶で価電子バンドが広がる幅(数 eV $\sim 0.1$–$0.3$ Ha)以上の窓があれば合格である．

(2) 極の本数と位置．$D_l$ は $\varepsilon$ の増加に対して単調減少し(式 \eqref{eq:15-normid} の右辺が負)，$u_l(r_d)=0$ になるたびに $-\infty$ から $+\infty$ へ跳ぶ．跳びの本数はその窓の中の束縛/共鳴状態の数に対応する．分離形だけに余分な極があればゴースト状態である(15.6.4項，図15.7)．

(3) セミローカルと分離形の差．(ii)と(iii)がずれていれば，それはKB/MBK変換そのものの問題(射影子の本数不足，局所ポテンシャルの選択ミス)である．(i)と(ii)がずれていれば擬ポテンシャルの生成そのものの問題である．原因の切り分けができる．

経験則として，対数微分が合っていればバンド構造も合う．逆に対数微分がずれていれば，ほぼ確実にバンド構造がずれる．数十秒で終わる原子計算だけで固体計算の成否をかなりの精度で予言できるので，擬ポテンシャルを作ったらまずこれを見る．

15.9.2 第二の関門:固体の物性を全電子法と比べる

対数微分は必要条件であって十分条件ではない．実際の固体で何が起きるかは，結局のところ固体で計算してみるしかない．標準的な検証項目は次の3つである．

(a) バンド構造．同じ交換相関汎関数(通常はPBE)を使い，全電子法—FLAPW(WIEN2k, FLEUR, exciting)や，LMTO系，あるいは数値基底全電子コード(FHI-aims, Elk)—と，擬ポテンシャル計算のバンド分散を重ねて描く．占有帯だけでなく，フェルミ準位から数 eV 上の非占有帯まで一致するかを見る．非占有帯は擬ポテンシャルの高エネルギー側の移植性を鋭敏に反映する．半芯状態を価電子に入れていない場合，その準位の位置がずれることが典型的な失敗である．

(b) 格子定数と体積弾性率．体積を平衡値のまわりで振ってエネルギー $E(V)$ を計算し，Birch–Murnaghanなどの状態方程式でフィットして平衡体積 $V_0$，体積弾性率 $B_0 = -V\left(\partial P/\partial V\right)_{V_0} = V_0\left(\partial^2 E/\partial V^2\right)_{V_0}$，その圧力微分 $B_0'$ を求める．全電子法との差が格子定数で $0.1\%$ 以内，$B_0$ で数 $\%$ 以内であれば良好である．これは対数微分より厳しい検査である:体積を変えるということは原子間距離を変えることで，擬ポテンシャルが作られたときとは異なる環境に置くことにほかならないからである．

(c) $\Delta$ ゲージ．(b)を1つの数値に凝縮した指標が $\Delta$ ゲージである．2つの計算手法 $a$，$b$ が与える $E(V)$ 曲線を，それぞれの最小値を基準にそろえた上で，平衡体積の $\pm6\%$ の区間で二乗平均差をとる:

$$ \begin{equation} \Delta_{ab} = \sqrt{\frac{1}{0.12\,V_0}\int_{0.94V_0}^{1.06V_0}\left[E_a(V) - E_b(V)\right]^2\dd V} \label{eq:15-delta-gauge} \end{equation} $$

単位は meV/atom である．前因子 $1/(0.12V_0)$ は積分区間の幅 $0.12V_0$ で割る規格化で，$\Delta$ が「区間内での典型的なエネルギー差」を表すようにしてある．71元素の単体結晶についてこの指標で15の主要コードを比較した大規模検証(Lejaeghere ら, 2016)によれば，全電子コード同士の $\Delta$ は $0.3$–$1$ meV/atom，よく作られたノルム保存・ウルトラソフト・PAWの擬ポテンシャルは全電子との差が $1$–$3$ meV/atom に収まる．参考までに，PBE汎関数そのものが実験と食い違う大きさは平均 $23.5$ meV/atom で1桁大きい．すなわち現代の擬ポテンシャルは，汎関数の誤差に比べれば無視できる水準まで来ている—ただしそれは「よく作られた」ものに限る，というのが本節の主題である．

15.9.3 半芯状態を価電子に入れるか

凍結内殻近似の前提は「内殻は環境によって変わらない」であった(15.1節)．この前提が怪しい殻を半芯状態(semicore state)という．カリウムの $3s3p$，チタンの $3s3p$，ガリウムの $3d$，インジウムの $4d$，バリウムの $5s5p$，希土類の $5s5p$ などが代表例である．価電子に含めれば精度は上がるが，これらの軌道は空間的にコンパクトなので必要カットオフが跳ね上がる(たとえば $3d$ を価電子に入れると $E_{\rm cut}$ が2〜3倍になる)．判断の基準を挙げる．

定義:半芯状態を価電子に含める判断基準

空間的重なり．その内殻軌道の動径確率密度 $r^2\left|R_{nl}\right|^2$ のピーク位置や裾が，価電子軌道のピーク位置と同程度，あるいは $r_c$ を越えて外に出ていないか．目安として，内殻軌道の平均半径 $\left\langle r\right\rangle$ が価電子の $r_c$ の $0.6$ 倍を超えたら要注意である．物理的には「隣の原子の価電子と重なるかどうか」であり，最近接原子間距離の半分と比べるのが本質的である．
エネルギー的近接．その準位が価電子準位から $2$–$3$ Ha 以内にあるか．ギャップが小さいほど化学環境による混成・シフトを受けやすい．
化学シフト．参照配置を変えて($\text{中性原子}\to\text{イオン}\to\text{励起配置}$)全電子計算を繰り返し，その準位が動く量を見る．数十 mHa 以上動くなら，その殻はもはや「凍っていない」．
非線形コア補正で足りるか．問題が交換相関の非線形性だけなら NLCC(15.5.6項)で救えることが多い．しかし混成や電荷移動が起きているなら，価電子に昇格させるしかない．
実地検証．結局は両方作って(b)(c)の検査に掛けるのが確実である．高圧下や酸化物のようにイオン間距離が短い系ほど半芯が必要になりやすい．

15.9.4 $r_c$ を大きくしすぎる誘惑

擬ポテンシャルの作成でいちばん多い失敗は，$r_c$ を大きくとりすぎることである．誘惑は強い．$r_c$ を大きくすれば擬波動関数はいくらでも滑らかになり，$E_{\rm cut}$ は下がり，計算は速くなる．しかも参照原子での固有値は定義により厳密に再現されるから，原子の段階では何の警告も出ない．危険は環境を変えたときに初めて現れる．

物理的意味:$r_c$ が支配するもの

$r_c$ は「本物を偽物に置き換える領域の大きさ」である．大きくするほど，

移植性が落ちる．ノルム保存が保証するのは参照点まわりの1次までの一致だけであり(式 \eqref{eq:15-taylor})，2次以降の誤差係数は $r_c$ が大きいほど大きくなる．対数微分の一致窓が目に見えて狭くなる．
電荷分布の詳細が失われる．球内の総電荷は保存するが，その分布が違う．$r_c$ が結合長の半分に近づくと，隣の原子の電子が「偽物の領域」に侵入し，静電気学的なごまかし(15.4節の物理的意味ボックス)が通用しなくなる．
複数チャネルの $r_c$ の不均衡も危ない．$s$ だけ大きく，$d$ だけ小さいというような組み合わせは，$\delta v_l = v_l - v_{\rm loc}$ の形を歪め，ゴースト状態を誘発しやすい．

実務上の目安は次のとおりである．(i) $r_c$ は，その $l$ の全電子価電子軌道の最外節より外，最大値の近くにとる(節を内側に取り込まないと擬波動関数を節なしにできない)．(ii) $r_c$ は想定する系の最近接原子間距離の半分を超えない．(iii) 迷ったら小さめから始め，対数微分の窓が十分広いことを確認しつつ，$E_{\rm cut}$ が許容範囲に入る最小限まで広げる．

例:移植性テストの実際

擬ポテンシャルの移植性を原子レベルで定量化する簡単なテストがある．参照配置(たとえば中性基底状態 $3s^23p^2$)で作った擬ポテンシャルを固定したまま，占有数だけを変えた複数の配置($3s^23p^1$，$3s^13p^3$，$3s^23p^{1}3d^{1}$，イオン化状態など)について，擬原子と全電子原子の両方で自己無撞着計算を回し，固有値と全エネルギー差を比べる．

合格の目安は，固有値のずれが $1$–$5$ mHa 以内，配置間の全エネルギー差(励起エネルギー)のずれが数 mHa 以内である．ずれが $10$ mHa を超えるようなら $r_c$ が大きすぎるか，半芯状態を入れるべきである．参照配置そのものでは誤差がゼロになるのが構成上当然なので，必ず参照とは違う配置で試すことがこのテストの要点である．同じ理由で，固体の検証も平衡格子定数だけでなく圧縮・膨張させた体積で行う．

15.9.5 検証のチェックリスト

段階	検査	合格の目安	不合格のときに疑うもの
原子	参照固有値の再現	厳密に一致(構成上)	実装のバグ
原子	擬波動関数が節を持たない	目視で確認	$r_c$ が節より内側
原子	$v^{\rm ion}\to -Z_v/r$	遠方で一致	スクリーニング解除の誤り
原子	対数微分(全電子 vs セミローカル)	$\pm 0.5$ Ha の窓で一致	$r_c$ が大きすぎる
原子	対数微分(セミローカル vs 分離形)	極の本数・位置が同一	ゴースト状態，$v_{\rm loc}$ の選択
原子	配置を変えた移植性テスト	固有値差 $\lesssim 5$ mHa	$r_c$，半芯の扱い，NLCC
固体	$E_{\rm cut}$ 収束	全エネルギー差 $<1$ mHa/atom	柔らかさ不足
固体	バンド構造 vs 全電子(FLAPW)	占有帯・数 eV 上まで一致	半芯状態，$l_{\max}$ 不足
固体	格子定数・体積弾性率	$a$ で $0.1\%$，$B_0$ で数 $\%$	移植性全般
固体	$\Delta$ ゲージ	$\lesssim 3$ meV/atom	同上

最後に強調しておく．擬ポテンシャルは「元素ごとに1度作れば終わり」の量ではあるが，その1度の作成が計算全体の精度の上限を決めてしまう．公開されている検証済みのライブラリ(SG15，PseudoDojo，GBRV，各コード付属のデータベースなど)を使うのが実務では第一選択であり，自作するなら本節の検証を全部通すのが最低限の作法である．

15.10 まとめ

動機．内殻電子は化学的に不活性だが，価電子は内殻との直交性のために核近傍で細かい節を持つ．この構造のせいで全電子波動関数を平面波展開すると $E_{\rm cut}\sim 10^3$ Ha が必要になり($N_{\rm PW}\propto E_{\rm cut}^{3/2}$，式 \eqref{eq:15-npw})，計算は非現実的になる．内殻を凍結するだけでは足りず，価電子の波動関数そのものを核近傍で滑らかに作り替える必要がある．
起源．OPW法で基底に直交性を組み込む工夫を演算子の言葉に翻訳すると，Phillips–Kleinman方程式 \eqref{eq:15-pk} が現れる．内殻への射影は反発ポテンシャル $\hat V_R = \sum_c(E-E_c)\left|\psi_c\right\rangle\left\langle\psi_c\right|$ を生み，核の引力と大きく相殺する(キャンセレーション定理)．ただしPK形はエネルギー依存・非局所・ノルム非保存という3つの欠点を残した．
合格基準．擬ポテンシャルが本物の代役を務めるには，境界 $r_c$ での対数微分 $D_l(\varepsilon, r_c)$ が，価電子が関与するエネルギー窓の全体で全電子のそれと一致すればよい(式 \eqref{eq:15-tan-eta})．これが移植性の定量的な中身である．
ノルム保存．本章の核心はWronskianから導いた恒等式 \eqref{eq:15-normid} $\;\partial D_l/\partial\varepsilon = -2u_l(r_c)^{-2}\int_0^{r_c}\left|u_l\right|^2\dd r\;$ である．$r_c$ 内のノルムを保存させれば，対数微分は値だけでなくエネルギー1次微分まで一致し，散乱は参照点まわりで $O\!\left((\varepsilon-\varepsilon_l)^2\right)$ の精度で一致する．副産物として球外の静電ポテンシャルも正しくなる．
構成法．Troullier–Martins法は $u^{\rm PS} = r^{l+1}e^{p(r)}$($p$ は $r^{12}$ までの偶数冪多項式)を置き，動径方程式を逆に解いて $v^{\rm scr} = \varepsilon_l + (l+1)p'/r + \left[p''+(p')^2\right]/2$ を得る(式 \eqref{eq:15-tm-vscr})．7個の係数は，$r_c$ での0〜4階微分の一致，ノルム保存，原点平滑条件 $c_2^2+(2l+5)c_4=0$ の7条件で決まる．最後に価電子のHartree項と交換相関項を差し引いて(スクリーニング解除)イオン擬ポテンシャルにするが，交換相関は非線形なので部分内殻密度による補正(NLCC)が必要になる．
分離形．$l$ 依存のセミローカル形は平面波行列要素が因子化せず $O(N_{\rm PW}^2)$ を要する．Kleinman–Bylander変換 \eqref{eq:15-kb} は非局所項をランク1の射影に置き換え，参照状態に対する作用を保ったまま行列要素を因子化して $O(N_{\rm PW}N_{\rm proj})$ に落とす．代償はゴースト状態で，局所ポテンシャルの選択に依存し，対数微分の余分な極として検出できる．
ウルトラソフトとMBK．ノルム保存を捨てれば擬波動関数はいくらでも滑らかにできるが，補償電荷 $Q_{ij}$，重なり演算子 $\hat S = 1 + \sum q_{ij}\left|\beta_i\right\rangle\left\langle\beta_j\right|$，一般化固有値問題 $\hat H\phi = \varepsilon\hat S\phi$ が現れる．同じ $Q_{ij}$ が $B$ 行列の非エルミート性でもあること($B_{ij}-B_{ji}^* = (\varepsilon_i-\varepsilon_j)Q_{ij}$)を使えば，多参照のままノルム保存に戻す道(MBK法)が開け，$\hat S=1$ の通常の固有値問題と多射影子KB形という単純な実装が得られる．
相対論．$1s$ 電子の速度は $v/c \simeq Z\alpha$ で，$Z=80$ では $0.58$ に達する．Dirac方程式から $G$ の2階方程式を導くと，遠心力項は $j$ に依らず $l(l+1)$ で，$\kappa$ 依存性はただ1項に集まる(式 \eqref{eq:15-koelling})．縮重度平均 $\kappa_{\rm av}=-1$ を取ればスカラー相対論，$j$ 依存ポテンシャルの差を取ればスピン軌道項になる(式 \eqref{eq:15-sr-so})．相対論の情報は擬ポテンシャル生成の段階だけで取り込め，固体計算は非相対論のままでよい．
検証．対数微分曲線(一致の窓・極の本数)，配置を変えた移植性テスト，全電子法とのバンド構造・格子定数・体積弾性率の比較，$\Delta$ ゲージ(式 \eqref{eq:15-delta-gauge})の順に厳しくなる．$r_c$ を大きくとりすぎる誘惑に注意し，半芯状態を価電子に入れるかは空間的重なりと化学シフトで判断する．

次章への橋渡し．本章で，コーン・シャム方程式を計算機で解くための道具立てはひととおり揃った．基底(第14章)，内殻の扱い(第15章)，そして自己無撞着ループ．しかし電子構造が求まっても，それだけでは物質の姿は決まらない．原子核はどこに落ち着くのか，結晶はどんな形に歪むのか，有限温度で原子はどう動くのか—これらに答えるには，全エネルギーを原子座標と格子で微分した力と応力が要る．次章ではHellmann–Feynmanの定理から出発し，基底が原子に付随して動く場合の補正(Pulay力)を導き，構造最適化と第一原理分子動力学へ進む．そこでは，本章で作った擬ポテンシャルの非局所項が力の表式にどう寄与するかも具体的に現れる．

15.11 演習問題

演習15.1 1s軌道のフーリエ変換と必要カットオフ

水素様原子の $1s$ 軌道 $\phi_{1s}(\rr) = \sqrt{Z^3/\pi}\;e^{-Zr}$ のフーリエ変換 $$\tilde\phi_{1s}(\bm q) = \int\dd^3 r\; e^{-i\bm q\cdot\rr}\phi_{1s}(\rr)$$ を計算し，$\tilde\phi_{1s}(q) \propto \left(q^2+Z^2\right)^{-2}$ となることを示せ．次に，運動エネルギーの期待値を $q$ 空間で $\left\langle T\right\rangle = \frac{1}{2}\int\frac{\dd^3q}{(2\pi)^3}q^2\left|\tilde\phi_{1s}\right|^2$ と書いたとき，これを $q_{\max}$ までで打ち切ったときの相対誤差が $1\%$ 以下になる $q_{\max}$ を $Z$ の関数として求め，シリコン($Z=14$)での $E_{\rm cut}=q_{\max}^2/2$ を評価せよ．

ヒント:球対称関数のフーリエ変換は，$\hat{\bm q}$ を極軸にとって $\int\dd\hat\rr\,e^{-i\bm q\cdot\rr} = 4\pi\frac{\sin qr}{qr}$(Rayleigh展開の $l=0$ 項)を使うと $\tilde\phi(q) = \frac{4\pi}{q}\int_0^\infty r\sin(qr)\phi(r)\dd r$ という1次元積分に帰着する．$\int_0^\infty r e^{-Zr}\sin(qr)\dd r$ は $\sin(qr) = \mathrm{Im}\,e^{iqr}$ と書いて $\int_0^\infty re^{-(Z-iq)r}\dd r = (Z-iq)^{-2}$ を使えばよい．運動エネルギーの積分は $\int_{q_{\max}}^{\infty}\frac{q^4\dd q}{(q^2+Z^2)^4}$ の形になり，$q\gg Z$ では被積分関数が $q^{-4}$ なので裾が重いことに注意．

演習15.2 恒等式の $R_l$ 表示と規格化不変性

中心的恒等式 \eqref{eq:15-normid} を，$u_l = rR_l$ を使って $R_l$ の言葉に書き直すと $$\frac{\partial D_l(\varepsilon, r_c)}{\partial\varepsilon} = -\frac{2}{\left[r_cR_l(r_c)\right]^2}\int_0^{r_c} r^2\left|R_l(r)\right|^2\dd r$$ となることを示せ．右辺の分子が3次元の動径確率密度の積分($\int_0^{r_c}\left|R_l\right|^2r^2\dd r$)であることを確認せよ．また，無次元対数微分 $\bar D_l \equiv r\,\dfrac{\dd}{\dd r}\ln R_l(r)$ を使うと恒等式がどう書き換わるかを求めよ．

ヒント:$\ln u_l = \ln r + \ln R_l$ を $r$ で微分すると $D_l = 1/r + \dd\ln R_l/\dd r$．$1/r$ は $\varepsilon$ に依らない定数なので $\partial D_l/\partial\varepsilon = \partial(\dd\ln R_l/\dd r)/\partial\varepsilon$ である．無次元版では両辺に $r_c$ が掛かる．

演習15.3 TM擬ポテンシャルの原点近傍

(a) 式 \eqref{eq:15-tm-vscr} から出発し，$p(r) = c_0 + c_2r^2 + c_4r^4 + c_6r^6 + \cdots$ を代入して $v_l^{\rm scr}(r)$ を $r^4$ の項まで展開せよ．$r^0$ と $r^2$ の係数が本文の \eqref{eq:15-tm-exp2} と一致することを確かめ，$r^4$ の係数を $c_2, c_4, c_6$ と $l$ で表せ．

(b) 平滑条件 \eqref{eq:15-tm-smooth} $c_2^2 + (2l+5)c_4 = 0$ から，$c_4$ の符号は必ず負であることを示せ．これは擬ポテンシャルの形についてどんな意味を持つか，\eqref{eq:15-tm-v0} と合わせて論じよ．

ヒント:(b) $c_2^2 \ge 0$ と $2l+5 > 0$ から直ちに従う．$c_4 < 0$ は $v''(0)=0$ を満たしつつ $r^4$ 以降で下に曲がることを意味する．(c) $p' \to c_1$，$p''\to 2c_2$ なので $(l+1)p'/r \to (l+1)c_1/r$ が発散する．$Z=0$ の「見かけ上のクーロン発散」が残ってしまうことになる．

演習15.4 Kleinman–Bylander形の限界とゴースト状態

(a) 参照擬波動関数 $\phi_l$ に対して $\left\langle\phi_l\right|\hat v^{\rm KB}_{\rm nl}\left|\phi_l\right\rangle = \left\langle\phi_l\right|\delta v_l\left|\phi_l\right\rangle$ が成り立つことを，式 \eqref{eq:15-kb} から直接示せ．

(b) いま，同じ $l$ チャネルの別の状態 $\psi$ が $\left\langle\delta v_l\phi_l\middle|\psi\right\rangle = 0$ を満たすとする．このとき $\left\langle\psi\right|\hat v^{\rm KB}_{\rm nl}\left|\psi\right\rangle$ と $\left\langle\psi\right|\delta v_l\left|\psi\right\rangle$ を比較し，KB形とセミローカル形の食い違いを述べよ．$\delta v_l$ が全域で反発的($\delta v_l > 0$)なら $\left\langle\psi\right|\delta v_l\left|\psi\right\rangle > 0$ であるのに対し，KB形は $0$ を返す．この差が偽の束縛状態(ゴースト)を許す機構であることを説明せよ．

(c) 局所ポテンシャル $v_{\rm loc}$ の選び方を $v_{\rm loc} = v_{l'}$ と変えると $\delta v_l = v_l - v_{l'}$ の符号が変わりうる．KBエネルギー $E_l^{\rm KB}$ の符号とゴースト発生のしやすさの関係を，(b)の議論をもとに論じよ．

ヒント:(a) 分子と分母に同じ $\left\langle\phi_l\right|\delta v_l\left|\phi_l\right\rangle$ が現れて約分される．(b) KB形はランク1なので，$\left|\delta v_l\phi_l\right\rangle$ と直交する部分空間には一切作用しない．(c) $E_l^{\rm KB}<0$ のときKB項は負の重みを持つランク1演算子になり，変分的にエネルギーを下げる．

参考文献

C. Herring, A new method for calculating wave functions in crystals, Phys. Rev. 57, 1169 (1940) — OPW法．
J. C. Phillips and L. Kleinman, New method for calculating wave functions in crystals and molecules, Phys. Rev. 116, 287 (1959) — Phillips–Kleinman変換．
D. R. Hamann, M. Schlüter, and C. Chiang, Norm-conserving pseudopotentials, Phys. Rev. Lett. 43, 1494 (1979) — ノルム保存条件の確立．
G. P. Kerker, Non-singular atomic pseudopotentials for solid state applications, J. Phys. C 13, L189 (1980) — 逆問題型構成法の原型．
G. B. Bachelet, D. R. Hamann, and M. Schlüter, Pseudopotentials that work: From H to Pu, Phys. Rev. B 26, 4199 (1982) — BHS，周期表全域の相対論的擬ポテンシャル表．
N. Troullier and J. L. Martins, Efficient pseudopotentials for plane-wave calculations, Phys. Rev. B 43, 1993 (1991) — 本章15.5節の構成法．
L. Kleinman and D. M. Bylander, Efficacious form for model pseudopotentials, Phys. Rev. Lett. 48, 1425 (1982) — KB分離形．
X. Gonze, R. Stumpf, and M. Scheffler, Analysis of separable potentials, Phys. Rev. B 44, 8503 (1991) — ゴースト状態の系統的解析．
S. G. Louie, S. Froyen, and M. L. Cohen, Nonlinear ionic pseudopotentials in spin-density-functional calculations, Phys. Rev. B 26, 1738 (1982) — 非線形コア補正．
D. Vanderbilt, Soft self-consistent pseudopotentials in a generalized eigenvalue formalism, Phys. Rev. B 41, 7892 (1990) — ウルトラソフト擬ポテンシャル．
I. Morrison, D. M. Bylander, and L. Kleinman, Nonlocal Hermitian norm-conserving Vanderbilt pseudopotential, Phys. Rev. B 47, 6728 (1993) — MBK法．
P. E. Blöchl, Projector augmented-wave method, Phys. Rev. B 50, 17953 (1994);G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999) — PAW法とUSPPとの関係．
D. R. Hamann, Optimized norm-conserving Vanderbilt pseudopotentials, Phys. Rev. B 88, 085117 (2013) — ONCV，現代的な多射影子ノルム保存型．
D. D. Koelling and B. N. Harmon, A technique for relativistic spin-polarised calculations, J. Phys. C 10, 3107 (1977) — スカラー相対論方程式．
L. Kleinman, Relativistic norm-conserving pseudopotential, Phys. Rev. B 21, 2630 (1980) — $j$ 依存擬ポテンシャルの $\rm SR/SO$ 分解．
K. Lejaeghere et al., Reproducibility in density functional theory calculations of solids, Science 351, aad3000 (2016) — $\Delta$ ゲージによる大規模検証．
W. E. Pickett, Pseudopotential methods in condensed matter applications, Comput. Phys. Rep. 9, 115 (1989) — 総説．
R. M. Martin, Electronic Structure: Basic Theory and Practical Methods (Cambridge University Press, 2004), 第11章 — 教科書としての標準的まとめ．