# 有効内殻ポテンシャルの理論

重元素は、明示的な全電子計算を 2 つの独立した意味で高価にします:単純に相関/収束させるべき電子が
多くなること、そして最も重い元素については、化学がほぼ完全に価電子殻で起こるにもかかわらず内殻電子
への相対論効果が無視できなくなることです。**有効内殻ポテンシャル(ECP)**(擬ポテンシャルとも呼ば
れる)は、化学的に不活性な内殻電子を、固定的でパラメータ化された 1 電子演算子で置き換え、電子数と
基底サイズの両方を削減しつつ、(良くフィットされた ECP では)内殻電子が本来必要とする明示的な取り
扱いなしにスカラー相対論効果を暗黙に折り込みます。本章は qc-rs が実装するセミローカル ECP 演算子、
その分かりにくい動径べき乗規約、そして重く強く縮約された原始関数に対して数値的に頑健にする半数値
評価法を導出します — `.design/24qc.scf-ecp.md` と qc-rs の `AGENTS.md` qc-ecp 規則に基づきます。

## 価電子のみの Hamiltonian

核電荷 $Z_A$ と ECP によって除去される内殻電子数 $N_A^{\text{core}}$ から、**有効核電荷**は
$Z_A^{\text{eff}}=Z_A-N_A^{\text{core}}$ であり、非相対論的な価電子のみの電子ハミルトニアンは

$$
\hat H = \sum_i\left[-\tfrac12\nabla_i^2 - \sum_A\frac{Z_A^{\text{eff}}}{r_{iA}} +
\sum_A\hat U_A^{\text{ECP}}(i)\right] + \sum_{i<j}\frac{1}{r_{ij}} + E_{NN}^{\text{eff}}, \qquad
E_{NN}^{\text{eff}} = \sum_{A<B}\frac{Z_A^{\text{eff}}Z_B^{\text{eff}}}{R_{AB}}.
$$

となります。この中のすべての部分 — 縮小された核引力、縮小された核反発、そして*新しい*
$\hat U_A^{\text{ECP}}$ 演算子自体 — は単一の置換から流れ出ます:全電子原子は単に
$N_A^{\text{core}}=0$、$\hat U_A^{\text{ECP}}=0$ の特殊ケースです。価電子数は
$N_e=\sum_AZ_A^{\text{eff}}-Q$($Q$ は分子電荷)から直接従い、これが qc-rs のスピン/多重度のパリティ
検査が使う電子数です — ECP 計算の電子数計算ロジックは全電子ロジックに取って付けた特殊ケースでは
なく、$Z_A$ の代わりに $Z_A^{\text{eff}}$ に適用された全電子ロジック*そのもの*です。

## セミローカルスカラー ECP 演算子

化学的に重要な部分は $\hat U_A^{\text{ECP}}$ 自体です。**セミローカル**ECP(qc-rs が実装する
scalar/AREP — 平均化相対論有効ポテンシャル — の種類)は、純粋な動径**局所チャネル**と、それぞれ
自身の動径補正を運ぶ**角運動量射影子**の和に分かれます:

$$
\hat U_A^{\text{scalar}} = U_{A,L}^{\text{loc}}(r_A) + \sum_{l=0}^{L-1}\sum_{m=-l}^{l}
|Y_{lm}^A\rangle\,\Delta U_{A,l}(r_A)\,\langle Y_{lm}^A|, \qquad
\Delta U_{A,l}(r) = U_{A,l}(r) - U_{A,L}^{\text{loc}}(r),
$$

$Y_{lm}^A$ は原子 $A$ を中心とする実球面調和関数、$L$ は ECP に明示的に表現される最高角運動量です
(角運動量 $l\ge L$ の電子は局所チャネルだけを見ます — ある角運動量を超えると、内殻はもはやチャネル
固有の補正を必要としないという物理的描像です)。射影子 $|Y_{lm}^A\rangle\langle Y_{lm}^A|$ は、中心
$A$ に対する角運動量 $l$ の電子が $U_{A,L}^{\text{loc}}+\Delta U_{A,l}$ を経験し、他のどの電子も
局所チャネルだけを経験することを意味します — これはまさに、擬ポテンシャルのチャネル依存性が、除去
された内殻軌道への直交性(Pauli 排他律)が異なる角運動量の電子にどの程度影響するかの違いを表すと
いう物理的言明です。公開されている ECP データは局所チャネルを、消費側が生のチャネルごとの値から
$\Delta U$ を計算する必要なく、直接 `ul` として表にしていることがあるため、qc-rs は入力段階で
静かに差ポテンシャルへ変換せず — データベースが提供するものを厳密に保存し、(必要なら)評価時に
のみ減算を適用します。

各動径チャネル自体は原始 Gauss 型動径関数の和です:

$$
U_{A,l}(r) = \sum_k C_{A,lk}\,r^{n_{A,lk}-2}\,e^{-\alpha_{A,lk}r^2}.
$$

**動径べき乗規約は装飾的な細部ではなく、文書化された落とし穴です。** 保存される整数 $n_{A,lk}$ —
NWChem/PySCF/libecpint のファイル規約で `RADI_POWER`、qc-mol で `EcpPrim::n` と呼ばれる — は、
上式が明らかに $r^{n-2}$ 依存性を持つにもかかわらず、qc-rs 自身のコードでは`-2` 補正を適用せず、
**基底ファイルから読んだままの値**が保持され渡されます。`-2` は既にファイル規約自体に焼き込まれて
います(PySCF の `ECPscalar` に対して検証済み);qc-rs 内でもう一度 2 を引いて「修正」すると、
補正が二重に適用され間違った動径べきを生みます。これはまさに [基底と AO の章](../20-guide/basis-and-ao.md)
がより一般的に警告する種類の規約の不一致です:保存された整数の意味は、それを生んだファイル形式に
よって定義されるのであって、教科書の式の指数がたまたまどう見えるかによってではありません。

完全な AO 行列要素は、すべての ECP 中心にわたる素朴な 1 電子期待値です:

$$
V^{\text{ECP}}_{\mu\nu} = \sum_A\langle\chi_\mu|\hat U_A^{\text{ECP}}|\chi_\nu\rangle.
$$

$U_A^{\text{ECP}}$ の射影子構造は真にセミローカルな演算子(単なる $1/r$ や恒等演算子に掛かる Gauss
関数ではなく、*中心 $A$ に対する電子の角運動量に依存*し、これは他所を中心とする任意の AO 基底で
対角ではありません)なので、libcint の演算子代数の完全に外側にあります — これが qc-ecp が、qc-rs
の Hamiltonian 組み立ての他の部分が築かれている本番の libcint 1 電子/2 電子積分ドライバを再利用
するのではなく、独自のカーネルを必要とする理由です。

## 原始関数の正規化

PySCF の `ECPscalar` に対して検証された ECP 動径原始関数の正規化は、

$$
N(l,\zeta) = \sqrt{\frac{2\,(2\zeta)^{l+3/2}}{\Gamma(l+3/2)}},
$$

原始関数に対する単位動径自己重なりに一致します — 構造的には通常の Gauss 型 AO 原始関数の正規化と
同じ発想です([基底と AO の章で導出](../20-guide/basis-and-ao.md))、ただし通常の GTO 角運動量指数
ではなく上記の ECP 動径べき乗規約に特化しています。完全な AO 行列は(ECP 専用の特注 AO 表現ではなく)
`qc-grid` の `GaussianShell` と球面値評価器を再利用するので、正規化・一般縮約・AO 順序はすべて
自動的に libcint 自身の規約と一致します — ECP 行列は直接 libcint AO 順序で構築され、下流での別途
の並べ替えパスは不要です。

## 評価:解析的角度因子、半数値動径求積

セミローカル射影子の角度依存性には、真の閉形式解析的取り扱い — 実球面調和角度因子、Gaunt 係数で
因子化された(球面調和関数を掛け合わせる際に角運動量を結合する同じ Gaunt/Clebsch-Gordan 機構)— が
ありますが、2 つの(場合によっては大きく離れた、非中心の)Gauss 型 AO 原始関数の積と ECP 自身の動径
Gauss 展開にわたる**動径**積分には、一般の非中心ケースについて同様にきれいな閉形式がありません。
qc-rs の本番経路は明示的に**半数値**です:厳密な角度因子、数値動径求積 — これは qc-rs 固有の妥協
ではなく PySCF/libecpint が使うのと同じ一般戦略です。

本番ドライバ `qc-ecp::grid_radial::ecp_matrix_grid` は PySCF の `nr_ecp.c` の直接アルゴリズム移植
です(Rust で独立に再実装され、その数値を再現します):**殻対ごとに適応的に精密化される入れ子の
Gauss-Chebyshev 動径グリッド**。積分被積分関数がどれほど鋭く尖っているかによらずすべての殻対に
1 つの動径点数を固定するのではなく、グリッドは粗く始まり(31 点、「LEVEL0」)、与えられた殻対/
チャネルブロックについて、そのブロックの動径積分が厳しい許容誤差($10^{-12}$)に収束するまでのみ
**倍加**します、最大 2047 点の master グリッド(「LEVEL_MAX」)まで:

:::{prf:algorithm} 適応的動径精密化 ECP 行列要素
:label: alg-ecp-adaptive-radial-ja

**入力:** 殻対/チャネルブロックの動径被積分関数(Gauss 動径関数の積と ECP 自身の動径展開)、
収束許容誤差 $\tau=10^{-12}$。
**出力:** その ECP 行列ブロックへの収束済み動径寄与。

1. 最も粗い精密化レベル $L_0$(31 個の master グリッド点)から始める。
2. 動径推定値 $E_L = w_{\text{scale},L}\sum_{i\in\text{subset}_L}w_i\,f(r_i)$ を評価する。これは
   $(0,\infty)$ 上の単一 master Gauss-Chebyshev グリッドのストライド部分サンプルであり、レベル依存
   の再スケーリング $w_{\text{scale},L}=2^{L_{\max}-L}$ を伴い、すべてのレベルが*同じ*積分を
   推定するようにする。
3. $|E_L-E_{L-1}|<\tau$ なら(または $L=L_0$ で比較すべき前の推定値がなければ)$E_L$ を受け入れる。
4. さもなくば点数を倍にし(レベル $L+1$ へ進む、*同じ* master グリッドのより密なストライド部分
   集合 — 新しいグリッドは構築されない)、繰り返す。
:::

この適応性が、この方法を頑健かつ高速に同時にするものです:**なめらかで、よく離れた**殻対
(化学的にはよくあるケースです — 大きな分子のほとんどの AO 対は、与えられた ECP 中心から遠い関数を
少なくとも 1 つ含む)は最も粗いレベルで収束し、$\sim31$ 個の動径点しかコストがかかりません、一方
**鋭く、強く縮約され、互いに近接した**対 — 実際に高い動径分解能を必要とする対 — だけが高価な
2047 点精密化のコストを払います。固定グリッド方式は、難しいケースについて安全であるために*至る
ところで*高価なグリッドを使わなければならないか、それらについて静かな不精確さのリスクを負うか
でしょう;適応精密化は、物理が実際にそれを要求する場所でのみ密グリッドの精度を得ます。

**閉形式解析的動径経路の代わりに、それと並んでこれが存在する理由。** qc-rs はまた、原理的には厳密
だが重元素の強く縮約された原始関数に対して文書化された数値的失敗様式を持つ、解析的な非中心動径
経路(`matrix::ecp_matrix_analytic`、修正球面 Bessel 関数と同じ Gaunt 因子化を経由)も持っています:
閉形式の式はそれ自身が発散しうる前置因子($e^{k^2/4p}\to\infty$、大きな指数について)と、消える
項($e^{-\zeta d^2}\to0$、大きな分離依存の指数について)の積に分かれ、浮動小数点演算は
$\infty\times0$ を、2 つの因子の組み合わせが解析的に与えるはずの有限の積を回復するのではなく
`NaN` と評価します。グリッドベースの動径形式はこれを回避します、その被積分関数
$e^{-\zeta(r-d)^2}\cdot M_\lambda(2\zeta dr)$(別々に発散/消滅する 2 つの因子からではなく、
*結合された*、既に有限な指数から直接構築される修正球面 Bessel 関数の引数)が、全体を通じて数値的
に有界のままだからです — これが、一般ケースでは解析経路ではなく `ecp_matrix_grid` が qc-rs の
本番既定である理由です、解析経路が交差検証には利用可能で有用なまま残っているにもかかわらず。

角度部分は、どちらの動径法が使われるかにかかわらず解析経路から変更なく引き継がれます — それは
動径戦略とは独立に一度検証され、再利用されます:「Type 1」(単一中心、同一原子の AO 対)と
「Type 2」(AO 対の中心が互いに、あるいは ECP 中心と、あるいはその両方と異なる)の角度構造の
どちらについても、Cartesian 単項式ベースの経路と球面調和/Gaunt 経路の両方が実装され、互いに
$\sim10^{-12}$〜$10^{-13}$ まで交差検証されており、各ケースでより速い方が既定として選ばれます
(高角運動量の Type 1 では Cartesian、その中心ごとのキャッシュが有利な Type 2 では Gaunt)。

## 検証済み例:ECP 基底を持つ重元素

```python
import qc
m = qc.chk.new(atom="I 0 0 0; I 0 0 2.666", ao="lanl2dz", unit="angstrom").scf(ref="r").run()
m.scf.energy, m.scf.converged
# (-22.307234220033017, True)
m.scf.energy_components
# {'core': -56.253929996038075, 'coulomb': 28.188845752921235, 'exchange': -3.968212743262761}
```

I₂ は原子あたり $Z=53$ を持ちます(全電子なら 106 電子)が、LANL2DZ の ECP は各ヨウ素の内殻の
大部分を除去し、SCF に明示的に現れるのは価電子($Z^{\text{eff}}$ ごと、53 よりはるかに少ない)
だけです。結果として得られる全電子エネルギー、$\sim-22.3\ E_h$ は、真の全電子ヨウ素二量体のもの
よりはるかに小さい大きさです — まさに ECP 計算の期待される特徴です:深く化学的に不活性な内殻
エネルギーは固定された $\hat U_A^{\text{ECP}}$ 演算子に折り込まれ、明示的なエネルギー寄与として
現れることは決してありません。上記の `core` 成分は既に、縮小された核引力*と* ECP 演算子自身の寄与
の両方を含みます — qc-rs は別個の「ECP エネルギー」項を報告しません、$\hat U_A^{\text{ECP}}$ は
アーキテクチャ上、単に価電子のみのコア Hamiltonian $T+V_{ne}[Z^{\text{eff}}]+\hat U^{\text{ECP}}$
のもう 1 つの部分であって、別途計算・報告される取って付けた補正ではないからです。

:::{exercise}
:label: ex-ecp-theory-ja

1. 動径べき乗規約は、動径式が明らかに $r^{n-2}$ を含むにもかかわらず、`RADI_POWER` を基底ファイル
   から読んだままに保ち、qc-rs 自身のコードでは `-2` 補正を適用しません。コードの指数を教科書の式
   に文字通り一致させるために qc-rs 内で「親切に」もう一度 2 を引くことが、明確化ではなくむしろ
   バグになる理由を一文で説明しなさい。
2. 適応的動径グリッドは、まだ収束していない殻対/チャネルブロックについてのみ点数を倍にします、
   最も細かい(2047 点)グリッドを一様に使うのではなく。*一様に粗い*グリッドを危険にする殻対の
   具体的なクラスは何で、*一様に細かい*グリッドを無駄にするものは何ですか。
3. 閉形式解析的 ECP 動径積分は、真の数学的値が有限であるにもかかわらず、重元素の強く縮約された
   原始関数について文字通り `NaN` と評価されえます。その原因となる 2 つの衝突する極限(1 つは
   発散、1 つは消滅)を特定し、同じ物理をグリッドベースの被積分関数に再編成することが、クランプ
   や特殊ケースで単に回避するのではなく問題を回避する理由を説明しなさい。
:::

:::{solution} ex-ecp-theory-ja
:class: dropdown

1. 公開されている NWChem/PySCF/libecpint のファイル規約の `RADI_POWER` は、その整数が*意味する*
   ものに既に `-2` が焼き込まれています — 保存された値は、文字通り読めば、動径式が必要とする指数
   そのものです。qc-rs 内でもう一度 2 を引くと補正が二重に適用され、この方法で読み込まれるすべての
   ECP のすべての原始関数について静かに間違った動径べきを生みます — クラッシュではなく微妙に
   間違った ECP エネルギーとしてしか現れない系統的で静かな数値誤差であり、無害な単純化ではなく
   真に危険な「修正」となります。
2. 一様に*粗い*グリッドは、**鋭く、強く縮約され、互いに近接した**殻対(大きな指数、小さな分離)
   で危険です — まさに動径被積分関数が最も強く尖り、正確に積分するのが最も難しい対であり、31 点
   グリッドはそれらを静かに分解不足にし、収束したように見えるが数値的に間違った ECP 行列要素を
   生みかねません。一様に*細かい*(2047 点)グリッドは無駄です、妥当なサイズの分子のほとんどの
   殻対はなめらかでよく離れており(与えられた ECP 中心から遠い関数を少なくとも 1 つ持つ)既に
   最も粗いレベルで収束するからです — すべての対について無条件に高価なグリッドコストを払うと、
   簡単な大多数に何の精度利益もなく全体の ECP 評価コストをおよそ $2047/31\approx66$ 倍に
   膨らませます。
3. 解析的動径積分は、結合された指数 $p$ と運動量様パラメータ $k$ が大きいと発散する前置因子
   $e^{k^2/4p}$ と、分離依存の指数が大きいと消える因子 $e^{-\zeta d^2}$ の積に分かれます — 真の
   積は有限です(真の解析的相殺)が、IEEE 浮動小数点演算は 2 つの因子が*別々に*計算されてから
   掛け合わされるため、$\infty\times0$ をその相殺を回復するのではなく `NaN` と評価します。
   グリッドベースの被積分関数は代わりに*結合された*指数引数を直接、単一の数値的に有界な式として
   構築します(求積内の $e^{-\zeta(r-d)^2}$、独立に計算される発散項と消滅項に決して分割されない)
   — これは同じ不安定な式に取って付けた回避策ではなく、同一の背後の物理を真に異なる(数値的に
   安定な)方法で組織化したものです、だからこそ有限のままであるためにクランプや特殊ケースを必要
   としません。
:::

次章、[陰的溶媒和について](solvation-theory.md) は、qc-rs がサポートする 2 つ目の Hamiltonian 修正
を扱います — しかし ECP の固定演算子とは異なり、真に密度依存であり、したがって SCF 自己無撞着
ループに、そして [線形応答の章](linear-response-cphf.md) に従えば応答カーネルにも参加しなければ
ならないものです。