密度汎関数理論 (DFT) における擬ポテンシャルの基礎と定式化：規格保存型、超ソフト型、PAW 法、および解析的 GTH ポテンシャルの数理

導入と擬ポテンシャルの物理的起源

密度汎関数理論 (DFT) を用いた電子状態計算、特に結晶バルクや界面などの固体物理・材料科学の計算においては、空間的な周期性を活かして波動関数を平面波（Plane Wave）基底で展開する手法が広く用いられる。しかし、すべての電子（全電子：All-Electron）を平面波で直接展開しようとすると、計算コストが事実上無限大に発散する深刻なボトルネックに直面する。この課題を克服するために導入されたのが擬ポテンシャル (Pseudopotential) の概念である。

本稿では、擬ポテンシャルの物理的起源である内核・価電子分離および直交化平面波 (OPW) 法による定式化から出発し、規格保存型 (Norm-Conserving)、超ソフト型 (Ultrasoft: USPP)、そして現代の標準である Projector Augmented Wave (PAW) 法への数学的発展、さらには数学的に閉じている解析的 Goedecker-Teter-Hutter (GTH) 汎関数の定式化について、数式展開を一切省略せずに解説する。

---

1. 内核・価電子分離と直交化平面波 (OPW) の定式化

1.1 内核・価電子分離と凍結コア近似の物理的妥当性

固体や分子などの電子状態を量子力学的に記述する際、系を構成する電子は物理的役割と局在性の観点から内核電子 (Core electrons) と価電子 (Valence electrons) の2つに大別される。

• 内核電子: 原子核の極めて近くに局在し、非常に深い束縛エネルギー（エネルギー準位 $\epsilon_c$）を持つ。これらは化学結合の形成や環境の変化（固体の凝集状態、分子構造の歪みなど）に対してほぼ不活性（Chemically inert）であり、孤立原子の中にある時と実質的に同じ状態を保つ。
• 価電子: 空間的に広がっており、隣接する原子の電子と混成して化学結合（共有結合、金属結合など）を形成する。物性物理や材料化学における電気伝導度、光学特性、磁性、化学的安定性などのほぼすべての性質は、この価電子の挙動によって決定される。

この物理的な非対称性に基づき、内核電子の状態をあらかじめ孤立原子における状態に固定し、価電子のみを動的な変数として扱う近似を凍結コア近似 (Frozen core approximation) と呼ぶ。この近似は、系の自由度を劇的に減少させ、電子状態計算のコストを大幅に削減する上での物理的妥当な出発点となる。

1.2 平面波展開における高カットオフエネルギーのボトルネック

価電子のみを解く場合であっても、周期的な境界条件を持つ結晶系などでは、波動関数を単純な平面波 (Plane Wave, PW) 基底で展開するのが自然である：

\[\psi_v(\mathbf{r}) = \sum_{\mathbf{G}} c_{\mathbf{G}} e^{i(\mathbf{k} + \mathbf{G})\cdot\mathbf{r}}\]

ここで $\mathbf{k}$ は第一ブリルアンゾーン内の波数ベクトル、$\mathbf{G}$ は逆格子ベクトルである。しかし、この単純な平面波展開を実際の電子状態計算にそのまま適用しようとすると、壊滅的な計算コスト（ボトルネック）に直面する。その原因は、内核領域（原子核の近傍）における波動関数の急峻な振動挙動にある。

原子核の極めて近傍では、クーロン引力 $-Z/r$ が発散的に強くなる。また、量子力学におけるパウリの排他原理（あるいはハミルトニアンの自己共役性）により、価電子の波動関数 $\psi_v$ は、より低いエネルギーにあるすべての内核電子の波動関数 $\psi_c$ と直交しなければならない（$\langle \psi_c \vert \psi_v \rangle = 0$）。 この直交性条件を満たすため、価電子波動関数は内核領域において激しく振動する節（Node）構造を持つことになる。 平面波基底は空間的に平滑な関数であるため、このような内核領域での急峻な振動を正確に再現するためには、極めて波長の短い（高波数の）平面波成分が大量に必要となる。これは、必要とする平面波の最大キネティックエネルギー（カットオフエネルギー $E_{\text{cut}} = \frac{1}{2} \vert \mathbf{k} + \mathbf{G}_{\text{max}}\vert ^2$）が数千から数万 Rydberg に達することを意味し、実質的に計算を不可能にする。

1.3 直交化平面波 (Orthogonalized Plane Wave: OPW) 法の定式化

この「内核領域における激しい振動」と「空間的平滑性」のジレンマを解決するために、Herring (1940) は直交化平面波 (Orthogonalized Plane Wave: OPW) 法を提案した。 OPW法の基本的なアイデアは、空間的に滑らかな波動関数（擬波動関数、Pseudo-wavefunction）$\vert \phi\rangle$ に対し、あらかじめ内核軌道 $\vert \psi_c\rangle$ と直交するように補正を加えることで、真の価電子波動関数 $\vert \psi_v\rangle$ を構成することである。

具体的には、価電子の真の波動関数 $\vert \psi_v\rangle$ を以下のように定式化する：

\[\vert \psi_v\rangle = \vert \phi\rangle - \sum_c \vert \psi_c\rangle \langle \psi_c \vert \phi \rangle\]

ここで、和 $\sum_c$ はすべての内核状態 $c$ にわたって取られる。内核状態はあらかじめ解かれており、互いに正規直交化されていると仮定する（$\langle \psi_c \vert \psi_{c’} \rangle = \delta_{cc’}$）。

この構成により、任意の内核状態 $\langle \psi_{c’}\vert $ との重なり（内積）を計算すると、

\[\begin{aligned} \langle \psi_{c'} \vert \psi_v \rangle &= \langle \psi_{c'} \vert \phi \rangle - \sum_c \langle \psi_{c'} \vert \psi_c \rangle \langle \psi_c \vert \phi \rangle \\ &= \langle \psi_{c'} \vert \phi \rangle - \sum_c \delta_{c' c} \langle \psi_c \vert \phi \rangle \\ &= \langle \psi_{c'} \vert \phi \rangle - \langle \psi_{c'} \vert \phi \rangle \\ &= 0 \end{aligned}\]

となり、擬波動関数 $\vert \phi\rangle$ がどのような滑らかな関数であっても、構成された $\vert \psi_v\rangle$ はすべての内核状態と厳密に直交することが保証される。物理的には、右辺第二項の $\sum_c \vert \psi_c\rangle \langle \psi_c \vert \phi \rangle$ が内核領域において急峻に変化し、滑らかな $\vert \phi\rangle$ から内核軌道成分を差し引くことで、価電子波動関数の持つべき正しい節構造（振動挙動）を自動的に再現する役割を担う。

1.4 ハミルトニアンの作用と有効固有値方程式の導出

価電子の真の波動関数 $\vert \psi_v\rangle$ は、系の真のハミルトニアン $\hat{H}$ に対する固有値方程式を満たす：

\[\hat{H} \vert \psi_v\rangle = \epsilon_v \vert \psi_v\rangle\]

ここで $\epsilon_v$ は価電子の真のエネルギー固有値である。この方程式に、OPW法による $\vert \psi_v\rangle$ の表現式を代入する：

\[\hat{H} \left( \vert \phi\rangle - \sum_c \vert \psi_c\rangle \langle \psi_c \vert \phi \rangle \right) = \epsilon_v \left( \vert \phi\rangle - \sum_c \vert \psi_c\rangle \langle \psi_c \vert \phi \rangle \right)\]

左辺のハミルトニアンを分配法則に従って作用させる：

\[\hat{H} \vert \phi\rangle - \sum_c \left( \hat{H} \vert \psi_c\rangle \right) \langle \psi_c \vert \phi \rangle = \epsilon_v \vert \phi\rangle - \sum_c \epsilon_v \vert \psi_c\rangle \langle \psi_c \vert \phi \rangle\]

ここで、内核状態 $\vert \psi_c\rangle$ は真のハミルトニアン $\hat{H}$ の固有状態であり、対応する固有値は $\epsilon_c$ である（$\hat{H} \vert \psi_c\rangle = \epsilon_c \vert \psi_c\rangle$）。これを左辺第二項の括弧内に代入すると、以下を得る：

\[\hat{H} \vert \phi\rangle - \sum_c \epsilon_c \vert \psi_c\rangle \langle \psi_c \vert \phi \rangle = \epsilon_v \vert \phi\rangle - \sum_c \epsilon_v \vert \psi_c\rangle \langle \psi_c \vert \phi \rangle\]

この式を整理し、擬波動関数 $\vert \phi\rangle$ に対する方程式の形に書き換える。左辺第二項を右辺に移項して $\vert \phi\rangle$ で括ると、

\[\hat{H} \vert \phi\rangle + \sum_c (\epsilon_v - \epsilon_c) \vert \psi_c\rangle \langle \psi_c \vert \phi \rangle = \epsilon_v \vert \phi\rangle\]

これを演算子の形式で書き直すことで、擬波動関数 $\vert \phi\rangle$ に課される以下の有効固有値方程式（シュレーディンガー方程式）が導出される：

\[\left( \hat{H} + \hat{V}_{NL} \right) \vert \phi\rangle = \epsilon_v \vert \phi\rangle\]

ここで、$\hat{V}_{NL}$ は以下のように定義される非局所ポテンシャル演算子 (Non-local potential operator) である：

\[\hat{V}_{NL} = \sum_c (\epsilon_v - \epsilon_c) \vert \psi_c\rangle \langle \psi_c\vert\]

1.5 物理的相殺効果：核引力と斥力ポテンシャル $\hat{V}_{NL}$

導出された有効シュレーディンガー方程式において、有効ポテンシャルは真のポテンシャル（ハミルトニアン $\hat{H} = \hat{T} + \hat{V}$ に含まれるもの）と、追加された非局所演算子 $\hat{V}{NL}$ の和 $\hat{V}{\text{eff}} = \hat{V} + \hat{V}_{NL}$ として振る舞う。

この有効ポテンシャルにおいて働く物理的な本質は、ポテンシャルの相殺効果 (Cancellation effect) である。

1. 内核状態の束縛エネルギーの深さ: 内核状態 $\vert \psi_c\rangle$ は原子核の非常に近くに束縛されているため、そのエネルギー固有値 $\epsilon_c$ は価電子のエネルギー固有値 $\epsilon_v$ に比べて極めて低い（$\epsilon_c \ll \epsilon_v$）。
2. $\hat{V}{NL}$ の符号**: $\epsilon_c < \epsilon_v$ であることから、エネルギー差 $\epsilon_v - \epsilon_c$ は常に**正の値**をとる。したがって、演算子 $\hat{V}{NL}$ は内核軌道 $\vert \psi_c\rangle$ が存在する領域において強力な斥力（反発力）**として働く。
3. クーロン引力との相殺: 原子核近傍では、ハミルトニアンの中の真のポテンシャル $\hat{V}$ は強い負のクーロン引力ポテンシャル（$-Z/r$ の発散）を持つ。しかし、同じ内核領域で正の斥力として働く $\hat{V}_{NL}$ が加わることによって、引力と斥力が互いに打ち消し合う。

この相殺効果により、全有効ポテンシャル $\hat{V} + \hat{V}_{NL}$ は内核領域において非常に弱く、かつ滑らかなものとなる。この弱い有効ポテンシャルの存在こそが、急激な振動を持たない平滑な擬波動関数 $\vert \phi\rangle$ を用いて価電子状態を記述できる（すなわち、少ない平面波基底で展開できる）理論的保証（擬ポテンシャル理論の基礎）である。

---

2. 規格保存型擬ポテンシャルと Kleinman-Bylander 展開

2.1 規格保存型擬ポテンシャルの設計指針と境界条件

直交化平面波 (OPW) 法は擬ポテンシャルの概念的基礎を与えたが、実際の定量的な計算に耐えうる擬ポテンシャルを生成するための近代的な設計指針を確立したのが、Hamann, Schlüter, Chiang (HSC, 1979) らによる規格保存型擬ポテンシャル (Norm-Conserving Pseudopotential) の提唱である。

規格保存型擬ポテンシャルの設計において、原子の全電子（All-Electron, AE）波動関数 $\Psi^{\text{AE}}{nl}(r)$ と、対応する擬波動関数（Pseudo-wavefunction, PS） $\Psi^{\text{PS}}{nl}(r)$ （ここで $n$ は主量子数、$l$ は方位量子数）は、選択されたカットオフ半径 $r_c$ に対して以下の境界条件および要請を満たさなければならない：

1. 実空間での一致 (Spatial Matching): カットオフ半径 $r_c$ 以上の領域（外部領域）において、擬波動関数と全電子波動関数は完全に一致する：

\[\Psi^{\text{PS}}_{nl}(r) = \Psi^{\text{AE}}_{nl}(r) \quad (\text{for } r \ge r_c)\]

1. 滑らかな接続と $C^2$ 連続性 (Smoothness and Continuity): 内核領域（$r < r_c$）において擬波動関数は節（Node）を持たず、滑らかに原点へと接続する。また、境界 $r = r_c$ において値、1階微分、および2階微分が連続である（$C^2$ 連続）：

\[\left. \Psi^{\text{PS}}_{nl}(r) \right\vert _{r_c} = \left. \Psi^{\text{AE}}_{nl}(r) \right\vert _{r_c}, \quad \left. \frac{d}{dr}\Psi^{\text{PS}}_{nl}(r) \right\vert _{r_c} = \left. \frac{d}{dr}\Psi^{\text{AE}}_{nl}(r) \right\vert _{r_c}, \quad \left. \frac{d^2}{dr^2}\Psi^{\text{PS}}_{nl}(r) \right\vert _{r_c} = \left. \frac{d^2}{dr^2}\Psi^{\text{AE}}_{nl}(r) \right\vert _{r_c}\]

1. 固有値の一致 (Eigenvalue Matching): ある参照原子配置（通常は孤立原子の基底状態）において、価電子の擬波動関数の固有値 $\epsilon^{\text{PS}}{nl}$ は、全電子の固有値 $\epsilon^{\text{AE}}{nl}$ と厳密に一致する：

\[\epsilon^{\text{PS}}_{nl} = \epsilon^{\text{AE}}_{nl}\]

1. 規格保存条件 (Norm-Conservation Condition): カットオフ半径 $r_c$ の内部における電荷量が、擬波動関数と全電子波動関数で厳密に等しい：

\[\int_0^{r_c} \vert \Psi^{\text{AE}}_{nl}(r)\vert ^2 r^2 dr = \int_0^{r_c} \vert \Psi^{\text{PS}}_{nl}(r)\vert ^2 r^2 dr\]

2.2 規格保存条件と対数微分エネルギー恒等式の数理的証明

規格保存条件は、単に内核領域の電荷密度を正しく記述するだけでなく、異なる化学的環境（例えば分子や結晶中での電荷再配分）へのトランスファラビリティ（遷移性・転移性、Transferability）を保証する上で極めて重要な数学的役割を果たす。これを示すのが、対数微分のエネルギー微分に関する恒等式である。

球対称ポテンシャル $V(r)$ 中の動径シュレーディンガー方程式（原子単位系 $e = \hbar = m_e = 1$、ただしここでは $\hbar^2/2m = 1/2$）を考える：

\[-\frac{1}{2} \frac{d^2}{dr^2} (r\Psi(r, \epsilon)) + \left[ \frac{l(l+1)}{2r^2} + V(r) \right] (r\Psi(r, \epsilon)) = \epsilon (r\Psi(r, \epsilon))\]

ここで、簡単のために $u(r, \epsilon) = r\Psi(r, \epsilon)$ と置くと、上記方程式は以下のように簡略化される：

\[-\frac{1}{2} u''(r, \epsilon) + \left[ \frac{l(l+1)}{2r^2} + V(r) \right] u(r, \epsilon) = \epsilon u(r, \epsilon) \quad \cdots \text{(A)}\]

式 (A) をエネルギー $\epsilon$ で偏微分する。偏微分を表す記号として $\dot{u}(r, \epsilon) = \frac{\partial}{\partial \epsilon} u(r, \epsilon)$ を導入すると：

\[-\frac{1}{2} \dot{u}''(r, \epsilon) + \left[ \frac{l(l+1)}{2r^2} + V(r) \right] \dot{u}(r, \epsilon) = u(r, \epsilon) + \epsilon \dot{u}(r, \epsilon) \quad \cdots \text{(B)}\]

次に、式 (A) に $\dot{u}$ を掛けたものから、式 (B) に $u$ を掛けたものを差し引く：

\[\begin{aligned} & \left( -\frac{1}{2} u'' \dot{u} + \left[ \frac{l(l+1)}{2r^2} + V(r) \right] u \dot{u} \right) - \left( -\frac{1}{2} \dot{u}'' u + \left[ \frac{l(l+1)}{2r^2} + V(r) \right] \dot{u} u \right) = \epsilon u \dot{u} - (u^2 + \epsilon \dot{u} u) \\ \Rightarrow \quad & -\frac{1}{2} \left( u'' \dot{u} - u \dot{u}'' \right) = -u^2 \\ \Rightarrow \quad & \frac{1}{2} \frac{d}{dr} \left( u' \dot{u} - u \dot{u}' \right) = u^2 \end{aligned}\]

この両辺を $r = 0$ から $r = r_c$ まで積分する：

\[\int_0^{r_c} u(r, \epsilon)^2 dr = \frac{1}{2} \left[ u'(r, \epsilon) \dot{u}(r, \epsilon) - u(r, \epsilon) \dot{u}'(r, \epsilon) \right]_0^{r_c}\]

原点において物理的に妥当な正則解は $u(0, \epsilon) = 0$ および $\dot{u}(0, \epsilon) = 0$ を満たすため、下限からの寄与は消失する。したがって、

\[\int_0^{r_c} u(r, \epsilon)^2 dr = \frac{1}{2} \left[ u'(r_c, \epsilon) \frac{\partial u(r_c, \epsilon)}{\partial \epsilon} - u(r_c, \epsilon) \frac{\partial u'(r_c, \epsilon)}{\partial \epsilon} \right] \quad \cdots \text{(C)}\]

ここで、波動関数の対数微分（Logarithmic derivative） $x_l(\epsilon)$ を $r = r_c$ において以下のように定義する：

\[x_l(\epsilon) = \left. \frac{d}{dr} \ln u(r, \epsilon) \right\vert _{r_c} = \frac{u'(r_c, \epsilon)}{u(r_c, \epsilon)}\]

対数微分 $x_l(\epsilon)$ をエネルギー $\epsilon$ で偏微分すると：

\[\frac{\partial x_l(\epsilon)}{\partial \epsilon} = \frac{\partial}{\partial \epsilon} \left( \frac{u'(r_c, \epsilon)}{u(r_c, \epsilon)} \right) = \frac{\frac{\partial u'(r_c, \epsilon)}{\partial \epsilon} u(r_c, \epsilon) - u'(r_c, \epsilon) \frac{\partial u(r_c, \epsilon)}{\partial \epsilon}}{u(r_c, \epsilon)^2}\]

この分子の符号を反転させると、式 (C) の右辺の括弧内と厳密に一致することがわかる。すなわち、

\[u(r_c, \epsilon)^2 \frac{\partial x_l(\epsilon)}{\partial \epsilon} = - \left[ u'(r_c, \epsilon) \frac{\partial u(r_c, \epsilon)}{\partial \epsilon} - u(r_c, \epsilon) \frac{\partial u'(r_c, \epsilon)}{\partial \epsilon} \right]\]

これより、式 (C) は以下のように書き換えられる：

\[\int_0^{r_c} u(r, \epsilon)^2 dr = -\frac{1}{2} u(r_c, \epsilon)^2 \frac{\partial x_l(\epsilon)}{\partial \epsilon}\]

ここで $u(r, \epsilon) = r\Psi(r, \epsilon)$ であるから、

• 左辺の積分は $\int_0^{r_c} \vert \Psi(r, \epsilon)\vert ^2 r^2 dr$ となる。
• 右辺の係数は $u(r_c, \epsilon)^2 = r_c^2 \vert \Psi(r_c, \epsilon)\vert ^2$ となる。

また、対数微分 $x_l(\epsilon)$ と $\Psi(r, \epsilon)$ の対数微分の関係を調べる。

\[\frac{u'(r_c, \epsilon)}{u(r_c, \epsilon)} = \frac{\Psi(r_c, \epsilon) + r_c \Psi'(r_c, \epsilon)}{r_c \Psi(r_c, \epsilon)} = \frac{1}{r_c} + \frac{\Psi'(r_c, \epsilon)}{\Psi(r_c, \epsilon)}\]

これをエネルギー $\epsilon$ で微分すると、第1項の $1/r_c$ は定数であるため消失する：

\[\frac{\partial}{\partial \epsilon} \left( \frac{u'(r_c, \epsilon)}{u(r_c, \epsilon)} \right) = \frac{\partial}{\partial \epsilon} \left( \frac{\Psi'(r_c, \epsilon)}{\Psi(r_c, \epsilon)} \right)\]

これらを代入して整理することにより、求める恒等式が得られる：

\[- \frac{r_c^2}{2} \left[ \frac{\partial}{\partial \epsilon} \left( \frac{\Psi'(r_c, \epsilon)}{\Psi(r_c, \epsilon)} \right) \right] = \frac{1}{\vert \Psi(r_c, \epsilon)\vert ^2} \int_0^{r_c} \vert \Psi(r, \epsilon)\vert ^2 r^2 dr\]

2.3 散乱位相シフトとトランスファラビリティの物理的意味

この恒等式は、擬ポテンシャルの品質（トランスファラビリティ）を保証する上で極めて重大な物理的帰結をもたらす。

量子散乱理論において、ある球対称ポテンシャルによる外向き波動関数の漸近挙動は、境界における波動関数の対数微分のみによって決定される。対数微分が一致することは、散乱の位相のズレを示す散乱位相シフト (Scattering phase shift) $\eta_l(\epsilon)$ が一致することと同値である。

規格保存条件：

\[\int_0^{r_c} \vert \Psi^{\text{AE}}_{nl}(r)\vert ^2 r^2 dr = \int_0^{r_c} \vert \Psi^{\text{PS}}_{nl}(r)\vert ^2 r^2 dr\]

および境界での値の一致 $\Psi^{\text{AE}}(r_c, \epsilon_0) = \Psi^{\text{PS}}(r_c, \epsilon_0)$ は、上記の恒等式を通じて、以下の関係を強いる：

\[\left. \frac{\partial}{\partial \epsilon} \left( \frac{\Psi^{\text{AE}\prime}(r_c, \epsilon)}{\Psi^{\text{AE}}(r_c, \epsilon)} \right) \right\vert _{\epsilon_0} = \left. \frac{\partial}{\partial \epsilon} \left( \frac{\Psi^{\text{PS}\prime}(r_c, \epsilon)}{\Psi^{\text{PS}}(r_c, \epsilon)} \right) \right\vert _{\epsilon_0}\]

すなわち、規格保存型擬ポテンシャルは、参照エネルギー $\epsilon_0$ において全電子系と散乱位相シフトが一致するだけでなく、エネルギーの1次変化に対しても位相シフトが一致することを保証する。 これにより、原子が分子や結晶などの化学結合環境に置かれ、周囲の電子密度変化によって価電子の有効エネルギーが $\epsilon_0$ から多少ずれた場合であっても、擬ポテンシャルは全電子系とほぼ同等の散乱特性を維持する。これが規格保存型擬ポテンシャルの高いトランスファラビリティの源泉である。

2.4 Kleinman-Bylander (KB) 変換：半局所から非局所形式へ

HSC型をはじめとする初期の規格保存型擬ポテンシャルは、角運動量 $l$ ごとに異なるポテンシャルが作用する半局所的 (Semi-local) 形式で記述されていた。

\[\hat{V}_{\text{SL}} = V_{\text{loc}}(r) + \sum_l \Delta V_l(r) \hat{P}_l\]

ここで、$\Delta V_l(r) = V_l^{\text{PS}}(r) - V_{\text{loc}}(r)$ は $l$ チャネルに対するポテンシャルと局所ポテンシャル $V_{\text{loc}}(r)$ との差分であり、$\hat{P}l = \sum{m=-l}^l \vert Y_{lm}\rangle \langle Y_{lm}\vert $ は球調和関数による角運動量投影演算子である。

この半局所形式を平面波基底を用いたシュレーディンガー方程式の求解に適用する場合、ハミルトニアンの行列要素 $\langle \mathbf{k}+\mathbf{G} \vert \hat{V}_{\text{SL}} \vert \mathbf{k}+\mathbf{G}’ \rangle$ の計算において、逆格子ベクトル $\mathbf{G}$ と $\mathbf{G}’$ のすべての組み合わせについて 2次元の二重積分を実行する必要がある。平面波の数を $N_p$ とすると、この計算コストは $O(N_p^2)$ でスケールし、大規模な計算において極めて重いボトルネックとなる。

この問題を解決するため、Kleinman と Bylander (1982) は、半局所ポテンシャルを完全に分離可能 (Separable) な完全非局所 (Fully non-local) 形式に変換する Kleinman-Bylander (KB) 変換を提案した。

KB変換後の非局所ポテンシャル演算子 $\hat{V}_{NL}^{\text{KB}}$ は、以下のように定式化される：

\[\hat{V}_{NL}^{\text{KB}} = \sum_{lm} \frac{\vert \chi_{lm}\rangle \langle \chi_{lm}\vert }{\langle \phi_{lm}^{\text{PS}} \vert \Delta V_l \vert \phi_{lm}^{\text{PS}} \rangle} \quad \text{with} \quad \vert \chi_{lm}\rangle = \Delta V_l \vert \phi_{lm}^{\text{PS}}\rangle\]

ここで、$\vert \phi_{lm}^{\text{PS}}\rangle$ は孤立原子の参照状態における擬波動関数（固有状態）である。

作用の保存の証明

この構成が、参照エネルギーにおいて元の半局所ポテンシャルと同じ作用を波動関数に与えることを証明する。 定義された $\hat{V}{NL}^{\text{KB}}$ を参照擬波動関数 $\vert \phi{lm}^{\text{PS}}\rangle$ に作用させる：

\[\hat{V}_{NL}^{\text{KB}} \vert \phi_{lm}^{\text{PS}}\rangle = \frac{\vert \chi_{lm}\rangle \langle \chi_{lm} \vert \phi_{lm}^{\text{PS}}\rangle}{\langle \phi_{lm}^{\text{PS}} \vert \Delta V_l \vert \phi_{lm}^{\text{PS}} \rangle}\]

ここで、分子のブラベクトル $\langle \chi_{lm}\vert $ に $\vert \chi_{lm}\rangle = \Delta V_l \vert \phi_{lm}^{\text{PS}}\rangle$ のエルミート共役 $\langle \chi_{lm}\vert = \langle \phi_{lm}^{\text{PS}} \vert \Delta V_l$ を代入すると、

\[\langle \chi_{lm} \vert \phi_{lm}^{\text{PS}}\rangle = \langle \phi_{lm}^{\text{PS}} \vert \Delta V_l \vert \phi_{lm}^{\text{PS}}\rangle\]

となる。これを代入すると：

\[\hat{V}_{NL}^{\text{KB}} \vert \phi_{lm}^{\text{PS}}\rangle = \vert \chi_{lm}\rangle \frac{\langle \phi_{lm}^{\text{PS}} \vert \Delta V_l \vert \phi_{lm}^{\text{PS}}\rangle}{\langle \phi_{lm}^{\text{PS}} \vert \Delta V_l \vert \phi_{lm}^{\text{PS}} \rangle} = \vert \chi_{lm}\rangle = \Delta V_l \vert \phi_{lm}^{\text{PS}}\rangle\]

これにより、参照状態に対しては、KB形式の非局所ポテンシャルは元の半局所ポテンシャルと全く同じ作用をもたらすことが証明された。

計算量の削減効果

KB形式では、ハミルトニアンの行列要素は次のように積の形に分離される：

\[\langle \mathbf{k}+\mathbf{G} \vert \hat{V}_{NL}^{\text{KB}} \vert \mathbf{k}+\mathbf{G}' \rangle = \sum_{lm} \frac{\langle \mathbf{k}+\mathbf{G} \vert \chi_{lm}\rangle \langle \chi_{lm} \vert \mathbf{k}+\mathbf{G}' \rangle}{\langle \phi_{lm}^{\text{PS}} \vert \Delta V_l \vert \phi_{lm}^{\text{PS}} \rangle}\]

したがって、各平面波状態について射影（プロジェクション）ベクトル $\langle \mathbf{k}+\mathbf{G} \vert \dots \rangle$ をあらかじめ計算しておけば、行列要素の評価は 1次元の配列の積（外積）だけで完了する。これにより、計算コストは各電子状態あたり $O(N_p)$ に低減され、平面波基底を用いたDFT計算の速度が劇的に向上した。

2.5 ゴースト状態（Ghost States）の数理と物理的起源

Kleinman-Bylander展開は計算効率を圧倒的に高める一方で、ゴースト状態 (Ghost States) と呼ばれる深刻な数値的不安定性を誘起するリスクを孕んでいる。

物理的・数学的起源

ゴースト状態とは、KB型の非局所ハミルトニアンにおいて、参照状態のエネルギー（価電子エネルギー）よりも遥かに低いエネルギー領域（あるいは価電子領域の直下）に現れる、物理的に存在しないはずの偽の束縛状態である。

数学的には、KB変換の分母にある期待値 $E_{\text{KB}, l} = \langle \phi_{lm}^{\text{PS}} \vert \Delta V_l \vert \phi_{lm}^{\text{PS}} \rangle$ の符号と、$\Delta V_l(r)$ 自体の実空間での符号が不整合を起こすことが原因である。 通常、ある角運動量チャネル $l$ において、$\Delta V_l(r)$ が全領域または一部の領域で引力（負のポテンシャル）として働くとする。しかし、もし分母の $E_{\text{KB}, l}$ が正（斥力的）になってしまうと、KBポテンシャル演算子：

\[\hat{V}_{NL}^{\text{KB}} = \frac{\vert \chi_{lm}\rangle \langle \chi_{lm}\vert }{E_{\text{KB}, l}}\]

は、本来引力であるはずの領域において強力な斥力として働いてしまう。逆に、$\Delta V_l(r)$ が斥力的である領域で $E_{\text{KB}, l} < 0$ になると、ポテンシャルは強力な引力として働き、本来束縛状態を持たないチャネル（例えば $d$ 軌道や $f$ 軌道など）において unphysical な束縛状態を発生させる。

Wronskian や Sturm-Liouville の定理（ノード定理）に従う通常のエルミートポテンシャルでは、固有状態はそのエネルギー順に節（Node）の数が増加する。しかし、KBポテンシャルは射影演算子（1次元のランク1演算子）であるため、通常のシュレーディンガー方程式の境界値問題に対するノード定理が適用されない。このため、ハミルトニアンのスペクトルにおいて、節を持たないか、あるいは想定外の節を持つ異常な固有状態が、基底状態または非常に低いエネルギーの束縛状態として滑り込んでしまう。これがゴースト状態の数学的実体である。

ゴースト状態の検出：対数微分の極（Pole）

ゴースト状態が存在するかどうかは、擬ハミルトニアンから計算される動径波動関数の対数微分 $x_l(\epsilon)$ をエネルギー $\epsilon$ の関数としてプロットすることで厳密に判定できる。

1. 正常な状態: 対数微分 $x_l(\epsilon)$ はエネルギーの上昇とともに単調に減少する。全電子系の対数微分に極（無限大への発散）が現れるのは、そのエネルギーにおいて束縛状態または共鳴状態が存在するときだけである。
2. ゴースト状態が存在する場合: 全電子系では滑らかな領域であるにもかかわらず、擬ポテンシャル系の対数微分 $x_l(\epsilon)$ に突如として極（Pole）が出現し、発散挙動を示す。この極は、そのエネルギーにおいて擬ポテンシャルが自己無撞着に偽の束縛状態を作っている動かぬ証拠である。

回避策と対策

ゴースト状態を排除するためには、以下の手法が取られる：

1. 局所ポテンシャル $V_{\text{loc}}(r)$ の選定変更: KB変換の基準となる局所ポテンシャル $V_{\text{loc}}(r)$ にどの角運動量 $l$ を選ぶか（通常、最も外側の遮蔽が強い $l$ や、最も関与の薄い $l$ を局所成分とする）によって $\Delta V_l(r)$ の形状が変化し、分母の符号や強度が変わるため、ゴースト状態を消滅させることができる。
2. カットオフ半径 $r_c$ の調整: 各 $l$ チャネルの $r_c$ をわずかに変更することで、参照擬波動関数 $\vert \phi_{lm}^{\text{PS}}\rangle$ の広がりと $\Delta V_l(r)$ の重なり積分が変化し、異常な引力相互作用を抑えることができる。
3. 多参照点（Multi-projector）化: 後述する Vanderbilt 型の超ソフト擬ポテンシャルや PAW 法のように、プロジェクタを複数導入することで、単一プロジェクタ起因の特異性を解消する。

---

3. Goedecker-Teter-Hutter (GTH) 解析的擬ポテンシャルの数理

3.1 解析的擬ポテンシャルの設計思想と導入の動機

前節までに述べた HSC 型などの規格保存型擬ポテンシャルや、後述する Vanderbilt 型の超ソフト擬ポテンシャルの多くは、動径波動関数や有効ポテンシャルが動径格子（Radial grid）上に離散化された数値データ（Tabulated data）として提供される。しかし、これらの数値的擬ポテンシャルを実空間グリッドや平面波基底で扱う際には、以下の実用上の課題が生じる：

1. グリッド補間誤差 (Grid interpolation errors): 動径格子上の数値を実空間の直交格子や球座性格子へ射影・補間する際、線形補間やスプライン補間などの内挿アルゴリズムに起因する数値誤差が発生する。これは系のエネルギーや力の計算に微小なノイズをもたらし、特に構造最適化や分子動力学（MD）計算におけるエネルギー保存則の厳密性を損なう原因となる。
2. 逆空間（平面波）表現の複雑さ: 数値的に定義されたプロジェクターやポテンシャルを逆空間で評価するためには、球ベッセル関数を用いた数値的なハンケル変換（Fourier-Bessel 変換）を実行する必要があり、数値積分の精度管理やカットオフ依存性の問題が付きまとう。

これらの課題を根本的に解決するために、Goedecker, Teter, および Hutter (1996) によって提案されたのが Goedecker-Teter-Hutter (GTH) 擬ポテンシャル である。GTH 擬ポテンシャルは、局所ポテンシャルおよび非局所ポテンシャルのプロジェクターがすべて解析的な関数（ガウス関数、誤差関数、多項式の積）として定義されている点に最大の特徴がある。

GTH 擬ポテンシャルの導入には、以下のような極めて強力な動機が存在する：

• CP2K などのハイブリッド基底法との圧倒的な親和性: CP2K で採用されている GPW (Gaussian and Plane Waves) 法のように、電子密度を記述するために平面波基底を、波動関数を記述するために原子軌道類似のガウス基底（Gaussian-type orbitals, GTO）をそれぞれ併用する手法において、ハミルトニアンの行列要素（特にポテンシャルエネルギー項）を解析的に評価することが極めて容易になる。ポテンシャル自体がガウス関数で構成されているため、ガウス基底との重なり積分がすべて閉じた解析的形式で評価でき、数値積分グリッドに起因する誤差が完全に排除される。
• フーリエ変換の単純性と精度: すべての関数要素が実空間においてガウス関数ベースで記述されているため、そのフーリエ変換（逆空間表現）もまたガウス関数と多項式の積という非常に単純な解析的閉形式で求まる。これにより、平面波基底を用いたエネルギーおよび力・応力テンソルの計算が極めて高速かつ数値的に安定に行われる。
• グリッド補間誤差の完全な排除: ポテンシャルの値を補間なしに任意の実空間座標で厳密に評価できるため、離散化格子に起因する数値的ノイズが存在しない。これは高精度の構造最適化や長時間の分子動力学シミュレーションにおいて絶大な威力を発揮する。
• 高いトランスファラビリティの確保: 解析的関数のパラメータ（ガウス幅 $r_{\text{loc}}$ や多項式係数 $C_i$）は、全電子計算による孤立原子の固有値、波動関数の対数微分、電荷密度に適合するように入念に最適化されており、単純な関数形でありながら多様な化学結合環境において高い遷移性を維持する。

---

3.2 局所ポテンシャル $V_{\text{loc}}(r)$ の定式化と物理的解釈

GTH 擬ポテンシャルにおける局所ポテンシャル $V_{\text{loc}}(r)$ は、長距離クーロン相互作用を記述する誤差関数（Error function, $\text{erf}$）の項と、短距離の原子核近傍における斥力・ orthogonalization 効果を補正するガウス変調多項式の項の和として、以下のように定義される：

\[V_{\text{loc}}(r) = -\frac{Z_{\text{ion}}}{r} \text{erf}\left( \frac{r}{\sqrt{2}r_{\text{loc}}} \right) + \exp\left( -\frac{r^2}{2r_{\text{loc}}^2} \right) \left[ C_1 + C_2 \left(\frac{r}{r_{\text{loc}}}\right)^2 + C_3 \left(\frac{r}{r_{\text{loc}}}\right)^4 + C_4 \left(\frac{r}{r_{\text{loc}}}\right)^6 \right]\]

ここで、各記号および物理的構成要素の意味は以下の通りである：

• $Z_{\text{ion}}$: 擬ポテンシャルによって遮蔽されたイオン芯の有効電荷（価電子数）。
• $r_{\text{loc}}$: 局所ポテンシャルの有効範囲を決定するパラメータ（ガウス電荷分布の広がり）。
• $C_1, C_2, C_3, C_4$: 各元素および価電子配置に応じてフィッティングされる展開係数。

長距離項の物理的起源：電荷密度のポアソン方程式

第1項の $-\frac{Z_{\text{ion}}}{r} \text{erf}\left( \frac{r}{\sqrt{2}r_{\text{loc}}} \right)$ は、点電荷としての原子核ではなく、空間的にガウス分布したマクロなイオン芯電荷密度 $\rho_{\text{core}}(r)$ が作る静電ポテンシャルに対応している。

実際、以下のガウス型電荷密度分布：

\[\rho_{\text{core}}(r) = -\frac{Z_{\text{ion}}}{(2\pi r_{\text{loc}}^2)^{3/2}} \exp\left( -\frac{r^2}{2r_{\text{loc}}^2} \right)\]

に対するポアソン方程式：

\[\nabla^2 V(r) = -4\pi \rho_{\text{core}}(r)\]

を解くことで、この項が厳密に導出される。

動径対称性を仮定すると、ポアソン方程式の動径成分は以下のようになる：

\[\frac{1}{r^2} \frac{d}{dr} \left( r^2 \frac{dV}{dr} \right) = -4\pi \rho_{\text{core}}(r)\]

両辺を原点から $r$ まで 2回積分し、境界条件 $V(r \to \infty) \to -Z_{\text{ion}}/r$ を適用すると、まさに静電ポテンシャルとして以下が得られる：

\[V(r) = -\frac{Z_{\text{ion}}}{r} \text{erf}\left( \frac{r}{\sqrt{2}r_{\text{loc}}} \right)\]

この項の漸近挙動および特異性の解消は以下のように数学的に理解できる：

1. 長距離挙動 ($r \gg r_{\text{loc}}$): 誤差関数の性質上、極限 $x \to \infty$ において $\text{erf}(x) \to 1$ となるため、ポテンシャルは通常の点電荷クーロン引力：

\[V_{\text{loc}}(r) \to -\frac{Z_{\text{ion}}}{r}\]

に漸近し、長距離の静電相互作用を正しく再現する。

1. 原点近傍 ($r \to 0$): 誤差関数のマクローリン展開：

\[\text{erf}(x) = \frac{2}{\sqrt{\pi}} \left( x - \frac{x^3}{3} + \mathcal{O}(x^5) \right)\]

を用いると、原点近傍におけるポテンシャルの極限値は：

\[\lim_{r \to 0} -\frac{Z_{\text{ion}}}{r} \text{erf}\left( \frac{r}{\sqrt{2}r_{\text{loc}}} \right) = \lim_{r \to 0} -\frac{Z_{\text{ion}}}{r} \cdot \frac{2}{\sqrt{\pi}} \left( \frac{r}{\sqrt{2}r_{\text{loc}}} \right) = -Z_{\text{ion}} \frac{\sqrt{2}}{\sqrt{\pi} r_{\text{loc}}}\]

となり、クーロン力固有の原点における $1/r$ の発散（特異性）が有限の値へと完全にかつ滑らかに正則化（Regularize）されていることがわかる。

短距離補正項の物理的解釈

第2項のガウス変調多項式は、内核電子との重なりによって生じるパウリの排他原理的な斥力効果（非局所ポテンシャルで表現しきれない局所的な寄与）や、擬波動関数を滑らかにするための補正を担う。 多項式の展開において偶数乗（$r^2, r^4, r^6$）のみが採用されているのは、原点 $r = 0$ におけるポテンシャルの傾き（1階微分）を厳密にゼロにし、原子核位置での物理的な反転対称性とポテンシャルの滑らかさ（$C^\infty$ 連続性）を担保するためである。

---

3.3 非局所ポテンシャル $V_{\text{nl}}(\mathbf{r}, \mathbf{r}’)$ と解析的ガウスプロジェクター

GTH 擬ポテンシャルにおける方位量子数 $l$ チャネルごとの非局所成分 $V_{\text{nl}}(\mathbf{r}, \mathbf{r}’)$ は、規格保存型の要請に基づき、以下のような完全分離型の形式で定義される：

\[V_{\text{nl}}(\mathbf{r}, \mathbf{r}') = \sum_{l=0}^{l_{\text{max}}} \sum_{m=-l}^l \sum_{i,j} Y_{lm}(\hat{\mathbf{r}}) \beta_i^l(r) h_{ij}^l \beta_j^l(r') Y_{lm}^*(\hat{\mathbf{r}}')\]

ここで、各構成要素は以下のように定義される：

• $Y_{lm}(\hat{\mathbf{r}})$: 球調和関数。角度方向の対称性を決定する。
• $h_{ij}^l$: 方位量子数 $l$ チャネルにおける射影関数（プロジェクター）の相互作用の強さを表す結合行列（Coupling matrix）。実対称行列（$h_{ij}^l = h_{ji}^l$）である。
• $\beta_i^l(r)$: 動径方向の解析的ガウス射影関数（プロジェクター）。以下のように明示的に定式化される。

\[\beta_i^l(r) = \frac{\sqrt{2} r^{l + 2i - 2} \exp\left(-\frac{r^2}{2r_l^2}\right)}{r_l^{l + 2i - 0.5} \sqrt{\Gamma\left(l + 2i - 0.5\right)}}\]

ここで、$r_l$ は角運動量 $l$ に対するプロジェクターの広がり（有効半径）を表すパラメータであり、$\Gamma(z)$ はガンマ関数である。動径方向の射影インデックス $i, j$ は通常、価電子のエネルギー依存性を高精度にフィッティングするために、1つの角運動量チャネルにつき最大3つ程度（$i, j \in {1, 2, 3}$）まで導入される。

動径プロジェクターの規格化条件の数学的証明

この動径プロジェクター $\beta_i^l(r)$ は、動径積分において自乗積分が 1 に正規化されている。このことを厳密に証明する。

証明すべき関係式は以下の動径積分である：

\[I = \int_0^\infty \left[ \beta_i^l(r) \right]^2 r^2 dr = 1\]

定義式を代入して $I$ を書き下す：

\[I = \int_0^\infty \left( \frac{\sqrt{2} r^{l + 2i - 2} \exp\left(-\frac{r^2}{2r_l^2}\right)}{r_l^{l + 2i - 0.5} \sqrt{\Gamma\left(l + 2i - 0.5\right)}} \right)^2 r^2 dr\]

分子の自乗および分母の根号を展開して定数項を積分の外へ括り出すと：

\[I = \frac{2}{r_l^{2l + 4i - 1} \Gamma(l + 2i - 0.5)} \int_0^\infty r^{2l + 4i - 4} \exp\left(-\frac{r^2}{r_l^2}\right) r^2 dr\] \[I = \frac{2}{r_l^{2l + 4i - 1} \Gamma(l + 2i - 0.5)} \int_0^\infty r^{2l + 4i - 2} \exp\left(-\frac{r^2}{r_l^2}\right) dr\]

この定積分を評価するため、変数変換を行う。

\[x = \frac{r^2}{r_l^2} \quad \Rightarrow \quad r = r_l x^{1/2}\]

両辺を微分すると：

\[dr = r_l \cdot \frac{1}{2} x^{-1/2} dx\]

積分範囲は $r: 0 \to \infty$ から $x: 0 \to \infty$ へと変わらない。これらを積分項に代入すると：

\[\int_0^\infty r^{2l + 4i - 2} \exp\left(-\frac{r^2}{r_l^2}\right) dr = \int_0^\infty (r_l x^{1/2})^{2l + 4i - 2} e^{-x} \left( \frac{r_l}{2} x^{-1/2} \right) dx\] \[= \frac{r_l^{2l + 4i - 1}}{2} \int_0^\infty x^{l + 2i - 1} x^{-1/2} e^{-x} dx\] \[= \frac{r_l^{2l + 4i - 1}}{2} \int_0^\infty x^{l + 2i - 1.5} e^{-x} dx\]

ここで、オイラーのガンマ関数の定義：

\[\Gamma(z) = \int_0^\infty x^{z-1} e^{-x} dx\]

と比較すると、積分の被積分関数は $z - 1 = l + 2i - 1.5$ すなわち $z = l + 2i - 0.5$ の場合に対応していることがわかる。したがって：

\[\int_0^\infty x^{l + 2i - 1.5} e^{-x} dx = \Gamma(l + 2i - 0.5)\]

これらを元の式 $I$ に代入すると：

\[I = \frac{2}{r_l^{2l + 4i - 1} \Gamma(l + 2i - 0.5)} \cdot \left[ \frac{r_l^{2l + 4i - 1}}{2} \Gamma(l + 2i - 0.5) \right]\]

分子と分母で $r_l^{2l + 4i - 1}$ および $\Gamma(l + 2i - 0.5)$ が互いに相殺し、また係数 $2 \times \frac{1}{2} = 1$ となるため、最終的に以下が得られる：

\[I = 1\]

これにより、任意の角運動量 $l$、射影インデックス $i$、および広がりパラメータ $r_l$ に対し、動径プロジェクターが厳密に規格化されていることが数学的に証明された。

複数プロジェクター（多参照）と結合行列 $h_{ij}^l$ の役割

GTH 形式では、各 $l$ に対し複数のプロジェクター $\beta_i^l(r)$ (通常は $i = 1, 2$) を導入することで、従来の単一プロジェクター型の Kleinman-Bylander ポテンシャルが抱えていた「特定のエネルギーでのみ正確で、そこから離れるとトランスファラビリティが低下する」という弱点を解決している。結合行列 $h_{ij}^l$ は、これらのプロジェクター空間内での有効ハミルトニアンの行列要素に相当し、非局所斥力の強さを表現する。この対称行列の各成分は、局所項パラメータとともに全電子状態の散乱特性に最適に適合するように決定される。

このように、GTH 擬ポテンシャルはその極めて美しい解析的形式の中に、規格保存型の高い物理的要請と、多参照プロジェクターによる広範なエネルギー領域でのトランスファラビリティを両立させた、現代の第一原理電子状態計算（特に Gaussian and Plane Wave 法を基礎とするコード）において必要不可欠な数理的成果物である。

---

4. 超ソフト擬ポテンシャル (Ultrasoft Pseudopotential: USPP) の定式化

4.1 規格保存条件の限界と緩和の物理的動機

規格保存型擬ポテンシャル (Norm-Conserving Pseudopotential: NCPP) は、トランスファラビリティを保証する上で極めて優れた数学的基盤を持っていた。しかし、周期表の第1周期元素（窒素 $\text{N}$、酸素 $\text{O}$、フッ素 $\text{F}$ など）や、遷移金属（特に $3d$ 軌道を持つ元​​素）のように、価電子波動関数が原子核近傍で極めて鋭いピークを持つ（強く局在している）場合、規格保存条件を課したまま滑らかな擬波動関数を構成することは困難になる。 規格保存条件：

\[\int_0^{r_c} \vert \Psi^{\text{AE}}(r)\vert ^2 r^2 dr = \int_0^{r_c} \vert \Psi^{\text{PS}}(r)\vert ^2 r^2 dr\]

を満たしつつ、かつ波動関数の節を取り除くと、擬波動関数 $\Psi^{\text{PS}}(r)$ は必然的に $r_c$ 内部で急峻なピークを持たざるを得ない。この鋭いピークを平面波で展開するためには、非常に高いカットオフエネルギー（通常 $70\sim100\text{ Ry}$ 以上）が必要となり、計算コストが急増する。

この課題を克服するため、Vanderbilt (1990) は規格保存条件の緩和を提案した。すなわち、カットオフ半径 $r_c$ 内部における電荷量が等しいという制約を意図的に破り、擬波動関数の自乗積分を全電子波動関数のそれよりも小さくすることを許容する。これにより、擬波動関数を極めて滑らかに（平面波カットオフを $30\sim40\text{ Ry}$ 程度まで劇的に低減）することが可能となる。しかし、この緩和によって失われた「コア領域における電荷量」を数学的に補正し、かつ正規直交性を回復するための新たな枠組みが必要となる。

4.2 一般化された重なり演算子 $\hat{S}$ と一般化正規直交関係

規格保存条件を緩和した結果、得られる擬波動関数（以下、チルダを付して $\vert \tilde{\psi}i\rangle$ と表記）は、通常の意味での正規直交関係（すなわち、単純な重なり積分 $\langle \tilde{\psi}_i \vert \tilde{\psi}_j \rangle = \delta{ij}$）を満たさなくなる。この問題を解決するため、以下のような一般化された重なり演算子 (Generalized overlap operator) $\hat{S}$ を導入する：

\[\hat{S} = 1 + \sum_{I,\mu\nu} q_{\mu\nu}^I \vert \beta_\mu^I\rangle \langle \beta_\nu^I\vert\]

ここで、

• $I$ は各原子サイトのインデックスである。
• $\beta_\mu^I$ は、原子 $I$ に局在する Vanderbilt 型のプロジェクター関数であり、カットオフ半径 $r_c$ の外側ではゼロとなる。
• $q_{\mu\nu}^I$ は、規格保存条件の破れを補正するための電荷増強強度（あるいは欠損電荷、Charge augmentation strength）であり、全電子部分波 $\vert \phi_\mu^I\rangle$ と擬部分波 $\vert \tilde{\phi}_\mu^I\rangle$ の内積の差分として以下のように定義される：

\[q_{\mu\nu}^I = \langle \phi_\mu^I \vert \phi_\nu^I \rangle - \langle \tilde{\phi}_\mu^I \vert \tilde{\phi}_\nu^I \rangle = \int_{\Omega_I} \left[ \phi_\mu^{I*}(\mathbf{r}) \phi_\nu^I(\mathbf{r}) - \tilde{\phi}_\mu^{I*}(\mathbf{r}) \tilde{\phi}_\nu^I(\mathbf{r}) \right] d\mathbf{r}\]

ここで、$\Omega_I$ は原子 $I$ を中心とするコア領域（増強領域）である。 この重なり演算子 $\hat{S}$ を用いることで、超ソフト擬波動関数は以下の一般化正規直交関係を満たすように定式化される：

\[\langle \tilde{\psi}_i \vert \hat{S} \vert \tilde{\psi}_j \rangle = \delta_{ij}\]

4.3 一般化固有値方程式と系のハミルトニアン

チルダ付きの超ソフト擬波動関数 $\vert \tilde{\psi}_i\rangle$ が満たすべき有効シュレーディンガー方程式（Kohn-Sham 方程式）は、重なり演算子 $\hat{S}$ の導入に伴い、通常の固有値問題から以下の一般化固有値方程式 (Generalized eigenvalue equation) へと変更される：

\[\hat{H} \vert \tilde{\psi}_i\rangle = \epsilon_i \hat{S} \vert \tilde{\psi}_i\rangle\]

ここで、$\hat{H}$ は系の自己無撞着な有効ハミルトニアンである。この固有値方程式は、波動関数の直交化条件が重なり演算子 $\hat{S}$ を通じて記述されているため、通常の対角化アルゴリズムではなく、共役勾配法（CG法）や Davidson 法などを一般化固有値問題用に拡張した手法を用いて解かれる。

4.4 電荷密度の再構成と補償電荷密度 $Q_{\mu\nu}^I(\mathbf{r})$

規格保存条件が課されていない超ソフト擬波動関数 $\vert \tilde{\psi}_i\rangle$ から、単に $\vert \tilde{\psi}_i(\mathbf{r})\vert ^2$ を計算しただけでは、コア領域内の総電荷が不足するため、正しいポアソン方程式を解くことができず、またハートリー・交換相関エネルギーも正しく評価できない。したがって、物理的に正しい全価電子電荷密度 $\rho(\mathbf{r})$ を再構成する必要がある。

USPP において、電荷密度 $\rho(\mathbf{r})$ は滑らかな擬電荷密度と、コア領域に局在する補償電荷密度（あるいは電荷増強密度、Augmentation charge density）の和として以下のように再構成される：

\[\rho(\mathbf{r}) = \sum_i \vert \tilde{\psi}_i(\mathbf{r})\vert ^2 + \sum_{I,\mu\nu} \rho_{\mu\nu}^I Q_{\mu\nu}^I(\mathbf{r})\]

ここで、各変数および関数の定義は以下の通りである：

1. $\rho_{\mu\nu}^I$ は、チルダ付きの価電子状態が原子 $I$ のプロジェクターチャネル $\mu, \nu$ を占有する割合を示す密度行列（占有数）であり、以下のように定義される：

\[\rho_{\mu\nu}^I = \sum_i \langle \tilde{\psi}_i \vert \beta_\nu^I \rangle \langle \beta_\mu^I \vert \tilde{\psi}_i \rangle\]

1. $Q_{\mu\nu}^I(\mathbf{r})$ は、全電子部分波ペアと擬部分波ペアの積の差分として定義される増強電荷関数である：

\[Q_{\mu\nu}^I(\mathbf{r}) = \phi_\mu^{I*}(\mathbf{r})\phi_\nu^I(\mathbf{r}) - \tilde{\phi}_\mu^{I*}(\mathbf{r})\tilde{\phi}_\nu^I(\mathbf{r})\]

この増強電荷関数 $Q_{\mu\nu}^I(\mathbf{r})$ は原子コア近傍でのみ非ゼロとなる急峻な関数である。実際の数値計算（特に平面波基底）では、この急峻な関数をそのままグリッド上で扱うとカットオフエネルギーを下げられないため、その多重極モーメントを保存したまま滑らかな補助関数に置き換える「電荷増強の擬化（pseudization）」などの数学的処理が施される。

---

5. Projector Augmented Wave (PAW) 法の定式化と USPP との数理的比較

5.1 線形変換演算子 $\hat{\mathcal{T}}$ の導入と全電子波動関数へのマッピング

Projector Augmented Wave (PAW) 法は、Blöchl (1994) によって提唱された、全電子 (All-Electron: AE) 計算の厳密性と擬ポテンシャル法の計算効率の高さを統合した現代の決定的な電子状態計算手法である。 PAW法では、物理的に本物である激しく振動する全電子波動関数 $\vert \Psi\rangle$ と、空間的に極めて滑らかな（平面波展開に適した）擬波動関数 $\vert \tilde{\Psi}\rangle$ の間をつなぐ線形変換演算子 $\hat{\mathcal{T}}$ を導入する：

\[\vert \Psi\rangle = \hat{\mathcal{T}} \vert \tilde{\Psi}\rangle\]

この変換演算子 $\hat{\mathcal{T}}$ は、結晶中の結合領域（間隙領域、Interstitial region）においては恒等写像（すなわち $\vert \Psi\rangle = \vert \tilde{\Psi}\rangle$）として振る舞い、各原子 $a$ を中心とする球状のコア領域（増強領域 $\Omega_a$）の内部においてのみ、波動関数を激しく振動する AE 波動関数へと書き換える役割を持つ。したがって、$\hat{\mathcal{T}}$ は以下のように局所的な寄与の和として記述される：

\[\hat{\mathcal{T}} = 1 + \sum_a \hat{\mathcal{T}}^a\]

ここで、$\hat{\mathcal{T}}^a$ は原子 $a$ の増強球 $\Omega_a$ の外側では完全にゼロとなる局所変換演算子である。

5.2 部分波展開と双直交プロジェクター $\langle p_i^a \vert $ の数理

増強球 $\Omega_a$ の内部における物理を記述するため、原子 $a$ において全電子部分波 $\vert \phi_i^a\rangle$ と、対応する滑らかな擬部分波 $\vert \tilde{\phi}_i^a\rangle$ の正規直交ではない基底ペアを導入する。局所変換演算子 $\hat{\mathcal{T}}^a$ は、擬部分波を全電子部分波へ写像するように定義される：

\[\left( 1 + \hat{\mathcal{T}}^a \right) \vert \tilde{\phi}_i^a\rangle = \vert \phi_i^a\rangle \quad \Rightarrow \quad \hat{\mathcal{T}}^a \vert \tilde{\phi}_i^a\rangle = \vert \phi_i^a\rangle - \vert \tilde{\phi}_i^a\rangle\]

この変換を、任意の滑らかな擬波動関数 $\vert \tilde{\Psi}\rangle$ に適用するために、増強球の内部で定義され、擬部分波と双直交関係（Biorthogonality）を満たすプロジェクター関数 $\langle p_i^a \vert $ を導入する：

\[\langle p_i^a \vert \tilde{\phi}_j^a \rangle = \delta_{ij}\]

このプロジェクターの完全性（あるいは局所的な展開可能性）を利用することで、局所変換演算子は以下のように記述できる：

\[\hat{\mathcal{T}}^a = \sum_i \left( \vert \phi_i^a\rangle - \vert \tilde{\phi}_i^a\rangle \right) \langle p_i^a \vert\]

これを全体の線形変換式に代入することで、PAW法における最も基本的な関係式である部分波展開（波動関数の再構成式）が得られる：

\[\vert \Psi\rangle = \vert \tilde{\Psi}\rangle + \sum_a \sum_i \left( \vert \phi_i^a\rangle - \vert \tilde{\phi}_i^a\rangle \right) \langle p_i^a \vert \tilde{\Psi} \rangle\]

この式は、全電子波動関数 $\vert \Psi\rangle$ が、

1. 空間全体で定義された滑らかな擬波動関数 $\vert \tilde{\Psi}\rangle$
2. 各増強球内での全電子部分波による局所展開 $\sum_i \vert \phi_i^a\rangle \langle p_i^a \vert \tilde{\Psi} \rangle$
3. および対応する擬部分波による局所展開 $\sum_i \vert \tilde{\phi}_i^a\rangle \langle p_i^a \vert \tilde{\Psi} \rangle$ の差し引き の3つの成分に完全に分解できることを示している。

5.3 物理量演算子の変換と期待値の評価

PAW変換の最大の特徴は、任意の物理量演算子 $\hat{A}$ の全電子状態に対する期待値 $\langle \Psi \vert \hat{A} \vert \Psi \rangle$ を、滑らかな擬波動関数 $\vert \tilde{\Psi}\rangle$ と、変換された有効演算子 $\hat{\tilde{A}}$ を用いて厳密に計算できる点にある：

\[\langle \Psi \vert \hat{A} \vert \Psi \rangle = \langle \tilde{\Psi} \vert \hat{\mathcal{T}}^\dagger \hat{A} \hat{\mathcal{T}} \vert \tilde{\Psi} \rangle = \langle \tilde{\Psi} \vert \hat{\tilde{A}} \vert \tilde{\Psi} \rangle\]

ここで、有効演算子 $\hat{\tilde{A}}$ を展開すると、異なる原子サイト間の増強球が重ならない（$\Omega_a \cap \Omega_b = \emptyset$ for $a \ne b$）という仮定のもとで、以下の3項分解が得られる：

\[\hat{\tilde{A}} = \hat{A} + \sum_a \sum_{i,j} \vert p_i^a\rangle \left( \langle \phi_i^a \vert \hat{A} \vert \phi_j^a \rangle - \langle \tilde{\phi}_i^a \vert \hat{A} \vert \tilde{\phi}_j^a \rangle \right) \langle p_j^a \vert\]

これにより、いかなる複雑な物理量（例えば、キネティックエネルギー、静電ポテンシャル、磁気モーメント、コア電子の超微細構造など）も、

• 滑らかな擬波動関数のみから評価されるグローバルな寄与（第1項）
• 各原子コア近傍での AE 部分波による局所期待値と擬部分波による局所期待値の差分（第2項・第3項） の和として、極めて高い精度で効率的に計算可能となる。

5.4 全電荷密度の3項分解

物理量演算子の変換の具体的な適用として、系の電子電荷密度を考える。PAW法における総電荷密度 $\rho(\mathbf{r})$ は、以下の3つの成分の代数和として定義される：

\[\rho(\mathbf{r}) = \tilde{\rho}(\mathbf{r}) + \sum_a \left( \rho^a_{\text{AE}}(\mathbf{r}) - \tilde{\rho}^a(\mathbf{r}) \right)\]

1. グローバル擬電荷密度 $\tilde{\rho}(\mathbf{r})$: 空間全体を覆う滑らかなグリッド上で、擬波動関数から直接計算される：

\[\tilde{\rho}(\mathbf{r}) = \sum_n f_n \vert \tilde{\Psi}_n(\mathbf{r})\vert ^2 + \tilde{\rho}_{\text{core}}(\mathbf{r})\]

ここで $f_n$ は占有数、$\tilde{\rho}_{\text{core}}$ は滑らかな擬コア電荷密度である。 2.  **局所全電子原子電荷密度 $\rho^a_{\text{AE}}(\mathbf{r})$**:
増強球 $\Omega_a$ の内部でのみ定義される、本物の急峻な電荷密度：

\[\rho^a_{\text{AE}}(\mathbf{r}) = \sum_{i,j} D_{ij}^a \phi_i^a(\mathbf{r})\phi_j^{a*}(\mathbf{r}) + \rho^a_{\text{core}}(\mathbf{r})\]

ここで $D_{ij}^a = \sum_n f_n \langle \tilde{\Psi}_n \vert  p_j^a \rangle \langle p_i^a \vert  \tilde{\Psi}_n \rangle$ は原子 $a$ における局所密度行列、$\rho^a_{\text{core}}$ は原子核および内核電子による AE コア電荷密度である。 3.  **局所擬原子電荷密度 $\tilde{\rho}^a(\mathbf{r})$**:
同じく増強球 $\Omega_a$ の内部でのみ定義される、滑らかな電荷密度：

\[\tilde{\rho}^a(\mathbf{r}) = \sum_{i,j} D_{ij}^a \tilde{\phi}_i^a(\mathbf{r})\tilde{\phi}_j^{a*}(\mathbf{r}) + \tilde{\rho}^a_{\text{core}}(\mathbf{r})\]

ここで $\tilde{\rho}^a_{\text{core}}$ は増強球内での擬コア電荷密度である。

この分解により、クーロンポテンシャル（ハートリーポテンシャル）を解く際には、グローバルなポアソン方程式は平滑な電荷密度 $\tilde{\rho}(\mathbf{r})$ および補償電荷（Compensation charge）の和に対してのみ高速な FFT グリッド上で解き、原子近傍の急峻な静電相互作用は 1次元の動径格子上で解析的に閉じた形式で解くという、極めて合理的なハイブリッドスキームが実現される。

5.5 USPP と PAW 法の数学的対応：凍結コア・ゼロプロジェクター極限の証明

超ソフト擬ポテンシャル (USPP) 法と PAW 法は、開発の経緯は異なるものの、数学的には極めて密接に関連している。実際、USPP は PAW 法において特定の近似（凍結コア近似およびハミルトニアンの近似）を導入した極限として厳密に導出することができる。

これを証明するため、PAWにおける恒等演算子 $\hat{I}$ に対する期待値（すなわち、波動関数の一般化重なり積分）を考える。 全電子状態の正規直交性 $\langle \Psi_i \vert \Psi_j \rangle = \delta_{ij}$ から出発し、PAW の線形変換式 $\vert \Psi\rangle = \hat{\mathcal{T}} \vert \tilde{\Psi}\rangle$ を代入する：

\[\delta_{ij} = \langle \Psi_i \vert \Psi_j \rangle = \langle \tilde{\Psi}_i \vert \hat{\mathcal{T}}^\dagger \hat{\mathcal{T}} \vert \tilde{\Psi}_j \rangle\]

ここで、重なり演算子を $\hat{S} = \hat{\mathcal{T}}^\dagger \hat{\mathcal{T}}$ と定義する。 演算子 $\hat{\mathcal{T}} = 1 + \sum_a \sum_i \left( \vert \phi_i^a\rangle - \vert \tilde{\phi}_i^a\rangle \right) \langle p_i^a \vert $ を代入して、$\hat{\mathcal{T}}^\dagger \hat{\mathcal{T}}$ を直接計算する。異なる原子サイト間の増強球が重ならない（$\Omega_a \cap \Omega_b = \emptyset$ for $a \ne b$）と仮定すると：

\[\begin{aligned} \hat{S} &= \left( 1 + \sum_a \sum_i \vert p_i^a\rangle \left( \langle \phi_i^a \vert - \langle \tilde{\phi}_i^a \vert \right) \right) \left( 1 + \sum_b \sum_j \left( \vert \phi_j^b\rangle - \vert \tilde{\phi}_j^b\rangle \right) \langle p_j^b \vert \right) \\ &= 1 + \sum_a \sum_{i,j} \vert p_i^a\rangle \left( \langle \phi_i^a \vert \phi_j^a \rangle - \langle \phi_i^a \vert \tilde{\phi}_j^a \rangle - \langle \tilde{\phi}_i^a \vert \phi_j^a \rangle + \langle \tilde{\phi}_i^a \vert \tilde{\phi}_j^a \rangle \right) \langle p_j^a \vert \\ & \quad + \sum_a \sum_i \vert p_i^a\rangle \left( \langle \phi_i^a \vert - \langle \tilde{\phi}_i^a \vert \right) + \sum_a \sum_j \left( \vert \phi_j^a\rangle - \vert \tilde{\phi}_j^a\rangle \right) \langle p_j^a \vert \end{aligned}\]

ここで、PAW法における擬部分波 $\vert \tilde{\phi}_i^a\rangle$ とプロジェクター $\langle p_j^a \vert $ の双直交関係を拡張し、増強球の外部においてプロジェクター $\langle p_j^a \vert $ の外挿に関して、以下の関係式が満たされていると仮定する：

\[\vert \phi_i^a\rangle - \vert \tilde{\phi}_i^a\rangle = \sum_k \left( \vert \phi_k^a\rangle - \vert \tilde{\phi}_k^a\rangle \right) \langle p_k^a \vert \tilde{\phi}_i^a \rangle\]

このとき、上記展開式の中の交差項は相殺され、さらに $q_{ij}^a$ を以下のように定義する：

\[q_{ij}^a = \langle \phi_i^a \vert \phi_j^a \rangle - \langle \tilde{\phi}_i^a \vert \tilde{\phi}_j^a \rangle\]

すると、重なり演算子 $\hat{S}$ は以下のように簡略化される：

\[\hat{S} = 1 + \sum_a \sum_{i,j} q_{ij}^a \vert p_i^a\rangle \langle p_j^a \vert\]

この式は、原子サイトインデックス $a$ を USPP のイオンサイト $I$ に、プロジェクター $\langle p_i^a \vert $ を Vanderbilt プロジェクター $\langle \beta_\mu^I \vert $ に、そして $q_{ij}^a$ を電荷増強強度 $q_{\mu\nu}^I$ にそれぞれ置き換えることで、USPPにおける一般化重なり演算子の定義式と厳密に一致する：

\[\hat{S} = 1 + \sum_{I,\mu\nu} q_{\mu\nu}^I \vert \beta_\mu^I\rangle \langle \beta_\nu^I\vert\]

同様に、PAW の電荷密度再構成式において、全電子コア状態を完全に凍結（Frozen core）し、コア-価電子相互作用を擬ポテンシャルパラメータに吸収させることで、USPPの電荷密度再構成式が導かれる。 したがって、USPP は、PAW法における「全電子波動関数の再生および任意の演算子の厳密な再構成」という一般化された枠組みから、期待値評価の対象を恒等演算子（重なり演算子）および有効ハミルトニアンに限定し、その他の物理量は擬化されたままで評価する「frozen-coreかつ演算子再構成を限定した簡略極限」として位置づけられる。この数学的一致は、両手法の理論的連続性を示すものであり、現代の第一原理計算コード（Quantum ESPRESSO や VASP など）において、両者が同列に扱われる背景となっている。

---

結論：擬ポテンシャル理論の数理的発展の総括

本稿では、密度汎関数理論 (DFT) における平面波基底計算のボトルネックを解消するための擬ポテンシャル理論について、その黎明期から現代に至る数理的発展を体系的に解説した。

1. OPW 法から規格保存型 (NCPP) へ: 内核電子と価電子の物理的役割の分離を出発点とし、OPW 法によるポテンシャルの物理的相殺効果を数学的に検証した。さらに、HSC らの設計指針と対数微分のエネルギー微分に関する恒等式の証明を通じて、規格保存条件がトランスファラビリティを決定づける本質であることを明らかにした。また、KB 変換による $O(N_p)$ スケーリングの達成と、それに伴うゴースト状態の数学的・物理的起源およびその判定法・回避策を示した。
2. 解析的擬ポテンシャル (GTH) の洗練: 数値格子の補間誤差を完全に排除する設計思想として、GTH ポテンシャルを詳述した。ガウス電荷分布のポアソン方程式からの長距離項の厳密な導出、およびガンマ関数を用いた動径プロジェクターの規格化条件の完全な証明を通じて、その解析的整合性を示した。
3. 超ソフト型 (USPP) から PAW 法への統合: 局在性の強い軌道に対する規格保存条件の限界を打破するため、電荷増強強度 $q_{\mu\nu}^I$ を用いた USPP の一般化直交条件および一般化固有値方程式への発展を示した。そして、本物（全電子）の波動関数と擬波動関数を線形変換演算子 $\hat{\mathcal{T}}$ でつなぐ PAW 法へと昇華させ、USPP がその frozen-core 近似および演算子近似の下での厳密な極限であることを数理的に証明した。

この理論的 Progression は、計算の「簡便さと速度」（擬ポテンシャル法）と「物理的厳密さと精度」（全電子法）という、一見相反する二つの要請を数学的トリックによって見事に調和させた結晶である。現在、この数理設計はマテリアルズ・インフォマティクスや大規模量子ダイナミクスなどのフロンティアを支える強固な基盤となっており、その数理的厳密性を理解することは、第一原理計算の結果を正しく解釈し、新たな理論を構築する上で不可欠な道標である。