test

6 minute read

導入

前稿「密度汎関数理論 (DFT) における擬ポテンシャルの基礎と定式化：規格保存型、超ソフト型、PAW 法、および解析的 GTH ポテンシャルの数理」では、擬ポテンシャルの物理的起源（直交化平面波法：OPW）から、各種定式化の数学的詳細について論じた。

しかし、これらの数理モデルを実際のコンピュータ上で第一原理計算パッケージとして実装し、あるいは独自のポテンシャルを作成・検証するためには、プログラム側のデータ構造や処理アルゴリズムに踏み込む必要がある。本稿では、第一原理計算コード（特に Quantum ESPRESSO）で広く採用されている標準擬ポテンシャルフォーマット UPF (Unified Pseudopotential Format) v2 の XML 仕様と物理変数のマッピング、平面波コード内部でのデータフローと主要サブルーチン（プロジェクターのフーリエ変換、電荷密度の再構成、非局所ハミルトニアン作用素の計算、オーバーラップ演算子の作用）、そして実務としての ld1.x コードを用いた擬ポテンシャルの生成と検証手順について徹底的に解説する。

1. UPF 仕様の全体像と対数径格子の数理

1.1 対数径格子（Logarithmic Radial Grid）の必要性とその定式化

第一原理計算において、電子状態の計算は原子核近傍（$r \to 0$）の急峻な変化と、原子間・真空領域（$r \to \infty$）の平滑な変化の両方を正確に表現する必要がある。全電子計算はもちろん、擬ポテンシャルを用いた計算であっても、価電子の擬波動関数やプロジェクター、ポテンシャルの動的変化を実空間の球座標系で表現する際、一様な径方向格子（等間隔格子：$r_i = i \cdot dr$）を使用することは極めて非効率である。 なぜなら、原子核近傍のクーロン場やコアの節構造（またはプロジェクターの急峻な立ち上がり）を解像するために非常に小さな格子間隔 $dr$ を採用すると、外殻領域まで膨大な数の格子点が必要となり、メモリ使用量と計算コストが発散するためである。

この課題を解決するために、原点付近で格子点を密に配置し、遠方に向かうにつれて指数関数的に格子間隔を広げる対数径格子（Logarithmic Radial Grid）が導入される。代表的な格子点定義として、以下の2つの形式が広く用いられている。

1) 指数関数形式 (Exponential Form)

\[r_i = r_0 e^{i \cdot dx}\]

ここで、$r_0$ は原点付近の最小半径を決定するパラメータ、$dx$ は対数空間における格子間隔の増分、$i$ は格子点のインデックス（$i = 0, 1, \dots, N_{\text{grid}}-1$）である。この形式では、$i = 0$ において $r_0 > 0$ となり、厳密な原点 $r = 0$ は定義されないが、原子核の有限サイズや実質的な最小カットオフを十分に表現できる。

2) シフト付き指数関数形式 (Shifted Exponential Form)

\[r_i = r_0 \left( e^{i \cdot dx} - 1 \right)\]

この形式では、$i = 0$ のときに $r_i = 0$ となるため、厳密な原点（$r = 0$）を含めることができるという数学的利点を持つ。

1.2 格子積分重み（微分体積因子）の導出と径方向積分への適用

対数径空間上の関数 $f(r)$ に対して、径方向の積分

\[I = \int_{0}^{\infty} f(r) dr\]

を実行することを考える。連続変数 $r$ を対数格子インデックス $i$（連続変数とみなす）によって変数変換する。 まず、指数関数形式 $r(i) = r_0 e^{i \cdot dx}$ について考える。$i$ に関する微分は、

\[\frac{dr}{di} = r_0 e^{i \cdot dx} \cdot dx = r_i \, dx\]

となる。したがって、微小な格子間隔に対応する径方向の微分体積因子（積分重み） $dr_i$ は、以下のように導出される。

\[dr_i = \left( \frac{dr_i}{di} \right) di = r_i \, dx \, di\]

これを用いて、積分を $i$ についての離散和に書き換える。ここで離散化ステップを $\Delta i = 1$ とすると、離散積分は以下のように定式化される。

\[\int_{0}^{\infty} f(r) dr \approx \sum_{i} f(r_i) \frac{dr_i}{di} \Delta i = \sum_{i} f(r_i) r_i \, dx\]

この $r_i \, dx$ は、UPFフォーマット内では径方向積分の重み配列 rab(i) として格納されている：

\[\text{rab}(i) \equiv \frac{dr_i}{di} = r_i \, dx\]

したがって、径方向積分は単に $\sum_i f(r_i) \text{rab}(i)$ として簡潔に評価できる。

同様に、シフト付き指数関数形式 $r(i) = r_0 (e^{i \cdot dx} - 1)$ の場合、微分は、

\[\frac{dr}{di} = r_0 e^{i \cdot dx} \cdot dx = (r_i + r_0) dx\]

となる。よって、微分体積因子は、

\[dr_i = (r_i + r_0) dx \, di\]

となり、径方向積分の表現は以下のようになる。

\[\int_{0}^{\infty} f(r) dr \approx \sum_{i} f(r_i) (r_i + r_0) dx\]

この場合の積分重み配列は $\text{rab}(i) = (r_i + r_0) dx$ となる。

なお、3次元空間の球対称な関数についての積分、すなわち径方向体積積分

\[\int g(\mathbf{r}) d\mathbf{r} = \int_{0}^{\infty} g(r) \cdot 4\pi r^2 dr\]

を評価する場合、被積分関数は $f(r) = 4\pi r^2 g(r)$ となり、離散化された形式は、

\[\int g(\mathbf{r}) d\mathbf{r} \approx 4\pi \sum_i g(r_i) r_i^2 \text{rab}(i)\]

となる。指数関数形式のグリッドでは $\text{rab}(i) = r_i \, dx$ であるため、

\[\int g(\mathbf{r}) d\mathbf{r} \approx 4\pi \sum_i g(r_i) r_i^3 \, dx\]

という極めてシンプルな代数計算に帰着される。

1.3 UPF v2 XMLタグの物理的意味とマッピング

Quantum ESPRESSOなどで用いられる統一擬ポテンシャルフォーマット (Unified Pseudopotential Format: UPF) v2 は、擬ポテンシャルの全情報を保持する構造化XMLフォーマットである。以下に、主要なXMLタグとそこに格納されている物理データの数理的な対応関係を詳細に解説する。

1) <PP_GRID>

  • 物理的意味: 擬ポテンシャル計算で用いられる対数径格子および積分重みの定義を格納する。
  • 構造とデータ:
    • 格子点数 $N_{\text{grid}}$、パラメータ $dx$ などの基本属性。
    • 格子点の座標値を表す配列 $r_i$（<PP_R> タグ内）。
    • 各点での積分重み（微分因子） $\text{rab}(i)$（<PP_RAB> タグ内）。
  • 役割: 平面波基底と実空間径格子の相互変換（フーリエ変換）や、実空間における重なり積分を実行する際の基盤となる。

2) <PP_BETA>

  • 物理的意味: 超ソフト擬ポテンシャル (USPP) や Projector Augmented Wave (PAW) 法において重要な、各方位量子数 $l$ チャネルに対するVanderbiltプロジェクター（射影演算子） $\beta_{i}^l(r_p)$ を格納する。

  • 数理的表現: 非局所ハミルトニアン作用素におけるプロジェクター $\vert \beta_i \rangle$ の実空間表現 $\beta_i^l(r_p)$。これは、あらかじめ定義されたカットオフ半径 $r_c^l$ の外側で厳密にゼロとなる：

    \[\beta_i^l(r_p) = 0 \quad (\text{for } r_p > r_c^l)\]
  • 役割: 擬波動関数 $\vert \tilde{\psi} \rangle$ を局在した原子軌道成分へ射影し、非局所エネルギーや電荷密度の補正項を算出するために使用される。

3) <PP_DIAG>

  • 物理的意味: プロジェクター間の非局所結合マトリックス係数 $D_{ij}$（D行列）を格納する。

  • 数理的表現: ハミルトニアンの非局所部分 $\hat{V}{NL} = \sum{ij} D_{ij} \vert \beta_i \rangle \langle \beta_j \vert$ における係数行列 $D_{ij}$ である。これは、イオン対数ポテンシャルと局所ポテンシャル・全電子固有値の差分から決定される：

    \[D_{ij} = D_{ij}^{0} + \int V_{\text{loc}}(r) Q_{ij}(r) dr\]

    （あるいはPAW法における対応するエネルギー補正項）。

  • 役割: プロジェクターを作用させた後の結合強度を表し、非局所ポテンシャルエネルギー $\langle \tilde{\psi} \vert \hat{V}_{NL} \vert \tilde{\psi} \rangle$ を決定する。

4) <PP_QIJ> / <PP_QIJL>

  • 物理的意味: USPPおよびPAWにおいて、規格保存条件の緩和に伴って生じるコア領域の電荷密度の過不足（電荷不足）を補うための補強電荷（Augmentation Charge）密度関数 $Q_{ij}(r_p)$ またはその多重極展開成分を格納する。

  • 数理的表現: 真の全電子波動関数 $\psi_{i}^{\text{AE}}$ と擬波動関数 $\psi_{i}^{\text{PS}}$ のノルムの差を補うように定義される：

    \[Q_{ij}(\mathbf{r}) = \psi_i^{\text{AE}*}(\mathbf{r})\psi_j^{\text{AE}}(\mathbf{r}) - \psi_i^{\text{PS}*}(\mathbf{r})\psi_j^{\text{PS}}(\mathbf{r})\]

    <PP_QIJ> はその径方向の関数 $Q_{ij}(r_p)$、<PP_QIJL> は角度方向を球面調和関数で展開した際の $L$ チャネルごとの分量 $Q_{ij}^L(r_p)$ を保持する。

  • 役割: ポアソン方程式を解くための正しい全電荷密度 $\rho(\mathbf{r})$ の再構成に使用される。

5) <PP_PSWFC>

  • 物理的意味: 各原子軌道 $(n, l)$ に対する擬原子波動関数（Pseudo Atomic Wavefunction） $\chi_{nl}(r_p)$ を格納する。
  • 数理的表現: 対応する孤立原子の全電子原子軌道 $\phi_{nl}(r)$ から生成された、コア領域で節を持たない滑らかな径方向波動関数 $\chi_{nl}(r_p)$。
  • 役割:
    1. 自給自足（SCF）計算の開始時における初期試行波動関数の構築（LCAOベースの初期化）。
    2. ポピュレーション解析（Mulliken解析や Löwdin 電荷解析）におけるプロジェクター。
    3. DFT+U（Hubbard $U$）計算における局在軌道への射影基準。

6) <PP_RCLOC>

  • 物理的意味: 各チャネル $l$ におけるカットオフ半径 $r_c^l$ を格納する。

  • 数理的表現: この半径 $r_c^l$ より外側では、擬ポテンシャルおよび擬波動関数は、全電子ポテンシャルおよび全電子波動関数と厳密に一致する：

    \[V_{\text{PS}}(r) = V_{\text{AE}}(r), \quad \chi^{\text{PS}}(r) = \phi^{\text{AE}}(r) \quad (\text{for } r > r_c^l)\]
  • 役割: 擬ポテンシャルの非局所性が機能する範囲（コア領域）を規定し、計算アルゴリズム上の実空間カットオフの判定などに用いられる。

2. Quantum ESPRESSO 内部構造と呼び出しシーケンス・データフロー

Quantum ESPRESSO（以下、QE）などの平面波基底DFTコード内では、擬ポテンシャルから読み込まれた各種の径方向物理量や非局所演算子が、平面波基底および逆格子空間（$G$空間）へと展開され、セルフコンシステントフィールド（SCF）ループ内で動的に処理される。

ここでは、QEにおける擬ポテンシャル関連サブルーチンの呼び出しシーケンス、データフロー、および個々のコードが果たす役割について数学的定式化とともに詳述する。

2.1 呼び出しシーケンスの全体像

以下の図は、SCF計算の初期化フェーズから、波動関数の操作、ハミルトニアンの作用、一般化オーバーラップ演算子の適用、および電荷密度の再構成に至る呼び出しシーケンスと主要データの流れを示す。

flowchart TD
    subgraph Initialization ["初期化フェーズ (Initialization)"]
        read_pp["read_pp.f90<br>(UPF解析 & 配列初期化: betar, dion, qq 等)"]
        init_us_1["init_us_1.f90<br>(プロジェクターβの実空間からG空間への変換)"]
        read_pp -->|径方向グリッド・プロジェクター| init_us_1
    end

    subgraph SCF_Loop ["SCFループ (Self-Consistent Field Loop)"]
        direction TB
        subgraph Density_Update ["電荷密度更新"]
            addusdens["addusdens.f90<br>(超ソフト電荷補強を含む電荷密度再構成)"]
        end

        subgraph H_S_Operations ["ハミルトニアン & S行列作用"]
            calbec["calbec.f90<br>(並列プロジェクター重なり積分評価: becp = &lt;β\vertψ&gt;)"]
            add_vuspsi["add_vuspsi.f90<br>(非局所ハミルトニアン作用: V_NL ψ)"]
            s_psi["s_psi.f90<br>(一般化オーバーラップ演算子作用: S ψ)"]
            
            calbec -->|becp| add_vuspsi
            calbec -->|becp| s_psi
        end

        Density_Update -->|有効ポテンシャルの構築| H_S_Operations
        H_S_Operations -->|収束判定・新波動関数 ψ| Density_Update
    end

    init_us_1 -->|G空間プロジェクター| H_S_Operations

2.2 各モジュール・サブルーチンの役割と数理的詳細

2.2.1 read_pp.f90：UPFファイルのパースと内部配列の初期化

read_pp（および関連モジュール read_upf_v2）は、XML形式で記述された UPF v2 ファイルをパースし、以下の物理量を表す内部配列群をメモリ上に確保・格納する。

  • betar(r, nb): 実空間対数径格子 r 上で定義された、各擬原子軌道（プロジェクターインデックス nb）に対応する Vanderbilt 非局所プロジェクター $\beta_i(r)$。
  • dion(nb, nb): 原子固有の非局所エネルギーパラメータ行列 $D_{jk}^0$。ハミルトニアンの非局所ポテンシャル部分のベースとなる。
  • qq(nb, nb): 超ソフト擬ポテンシャル（USPP）における規格保存条件の緩和度を規定する $q_{jk}$ 行列（電荷不足量）。
  • qfunc(r, nb, nb): 実空間径格子 r 上での補強電荷密度関数 $Q_{jk}(r)$。

これらの配列は、プログラム内でグローバルに（あるいはモジュール変数として）保持され、後続のフーリエ変換やハミルトニアン演算の基盤となる。

2.2.2 init_us_1.f90：プロジェクターのフーリエ変換

実空間の対数径格子上で定義されている球対称プロジェクター $\beta_i^l(r)$ は、平面波基底（逆格子空間、すなわち $G$ 空間）での演算を行うため、以下の球ベッセル変換（フーリエ変換の動的成分）によって $G$ 空間のプロジェクター $\beta_i^l(G)$ へと変換される。

\[\beta_i^l(G) = 4\pi i^l \int_0^\infty \beta_i^l(r) j_l(Gr) r^2 dr\]

ここで、$l$ は方位量子数、$j_l(Gr)$ は第1種球ベッセル関数である。対数径格子 $r_p$ および積分重み配列 rab(p) を用いた離散化表現は以下のようになる。

\[\beta_i^l(G) \approx 4\pi i^l \sum_{p} \beta_i^l(r_p) j_l(G r_p) r_p^2 \text{rab}(p)\]

計算された $\beta_i^l(G)$ は、平面波のカットオフエネルギーに基づいて逆格子ベクトル $\mathbf{G}$ 上でテーブル化（tab または tab_dft 配列）され、SCFループ中の内挿計算を高速化する。

2.2.3 addusdens.f90：電荷密度の再構成

USPPおよびPAW法では、擬波動関数がコア領域でノルム保存条件を満たさないため、単に滑らかな波動関数の自乗和をとるだけでは物理的な電荷密度を再現できない。そのため、addusdens は滑らかな電荷密度 $\rho_{\text{smooth}}(\mathbf{r})$ にコア領域での補強電荷を加算し、以下の式に沿って全電荷密度 $\rho(\mathbf{r})$ を再構成する。

\[\rho(\mathbf{r}) = \rho_{\text{smooth}}(\mathbf{r}) + \sum_{I, \mu\nu} \rho_{\mu\nu}^I Q_{\mu\nu}^I(\mathbf{r})\]

ここで、$I$ は原子サイトのインデックス、$\rho_{\mu\nu}^I$ は原子 $I$ における局所的な占有数行列（一般化密度行列）であり、波動関数 $\psi_i$ とプロジェクターの重なり $becp$ を用いて以下のように定義される。

\[\rho_{\mu\nu}^I = \sum_{i} f_i \langle \psi_i \vert \beta_\mu^I \rangle \langle \beta_\nu^I \vert \psi_i \rangle = \sum_{i} f_i becp_{i,\mu}^{I*} becp_{i,\nu}^{I}\]

ここで $f_i$ はバンド $i$ の占有数である。$Q_{\mu\nu}^I(\mathbf{r})$ は原子 $I$ において実空間で急峻な変化を持つ補強電荷密度（Augmentation Charge Density）であり、二重グリッド（dense grid と smooth grid）間でのスプライン補間や高速フーリエ変換（FFT）を組み合わせて高精度に加算される。

2.2.4 calbec.f90：プロジェクターと波動関数の重なり積分

calbec は、すべてのバンド $i$ の波動関数 $\psi_i$ とプロジェクター $\beta_j$ との重なり積分 $becp_{i,j}$ を並列評価する。

\[becp_{i,j} = \langle \beta_j \vert \psi_i \rangle = \sum_{\mathbf{G}} \beta_j^*(\mathbf{k}+\mathbf{G}) \psi_i(\mathbf{k}+\mathbf{G})\]

この処理は、ハミルトニアンの適用（add_vuspsi）とオーバーラップ演算子の適用（s_psi）の両方で前提となる。QE内部では、この巨大な行列積を高性能な並列線形代数ライブラリ（BLASの DGEMM/ZGEMM、または並列化された scalapack などのルーチン）を用いて並列実行し、通信・計算コストを最小限に抑えるように設計されている。

2.2.5 add_vuspsi.f90：非局所擬ポテンシャルの適用

add_vuspsi は、波動関数 $\psi_i$ に対してハミルトニアンの非局所擬ポテンシャル演算子 $\hat{V}_{NL}$ を作用させ、その結果をハミルトニアン作用後の波動関数配列（$H\psi$）に蓄積する。

\[\hat{V}_{NL} \psi_i = \sum_{jk} D_{jk} \lvert \beta_j \rangle \langle \beta_k \vert \psi_i \rangle\]

あらかじめ calbec により $becp_{i,k} = \langle \beta_k \vert \psi_i \rangle$ が計算されているため、上式は以下のように簡略化して評価される。

\[\hat{V}_{NL} \psi_i = \sum_{jk} D_{jk} becp_{i,k} \lvert \beta_j \rangle\]

ここで、$D_{jk}$ は read_pp で読み込まれたイオン由来の定数 $D_{jk}^0$ に対して、電子密度の再構成（addusdens）に伴うスクリーニングポテンシャルの寄与を加えた実効的な $D$ 行列（deeq 配列）である。

2.2.6 s_psi.f90：一般化オーバーラップ演算子 $S$ の適用

超ソフト擬ポテンシャル（および PAW 法）を用いる場合、Kohn-Sham 方程式は一般化固有値問題 $H\psi = E S\psi$ へと拡張される。s_psi は、この一般化オーバーラップ演算子 $\hat{S}$ を任意の波動関数 $\psi_i$ に作用させる。

\[\hat{S} \psi_i = \psi_i + \sum_{jk} q_{jk} \lvert \beta_j \rangle \langle \beta_k \vert \psi_i \rangle\]

$becp$ 行列を利用することで、この作用は以下のように実行される。

\[\hat{S} \psi_i = \psi_i + \sum_{jk} q_{jk} becp_{i,k} \lvert \beta_j \rangle\]

ここで $q_{jk}$ は read_pp で初期化された qq 行列である。この演算は、波動関数の規格化、対角化処理中の残差 $\vert r_i \rangle = \vert H\psi_i - E_i S\psi_i \rangle$ の算出、および直交化手続き（グラムシュミット直交化等）において極めて重要な役割を果たす。

3. ld1.x による擬ポテンシャル生成実務と検証プロトコル

3.1 ld1.x の入力パラメータと構成カードの設計

Quantum ESPRESSO に同梱されている原子全電子・擬ポテンシャル計算コードである ld1.x は、任意の原子種の全電子計算を実行すると同時に、指定されたパラメータに従って規格保存型（NC）、超ソフト型（USPP）、あるいは PAW 型の擬ポテンシャルを生成し、さらにはその検証計算までを一貫して実行できる強力なツールである。

ld1.x の入力ファイルは主に3つのネームリスト（&input, &inputp, &test）と、それに続く構成カードから構成される。主要なパラメータとその物理的な意味は以下の通りである。

ネームリスト &input

  • title: 計算全体を識別するためのタイトル文字列。
  • zed: 対象原子の原子番号 $Z$（実数型または整数型）。
  • isxt または dft: 交換相関機能（DFT Functional）の指定。isxt（Exchange-Correlation Functional Index）は数値インデックス（1: LDA, 3: PBE など）を用い、dft は文字列（’LDA’, ‘PBE’, ‘PBEsol’ など）を用いて指定する。生成された擬ポテンシャルファイルを実用する際には、生成時と同じ機能を使用することが一貫性の観点から極めて重要である。
  • rel: 相対論効果の取り扱い。
    • rel=0: 非相対論的シュレーディンガー方程式。
    • rel=1: スカラー相対論的シュレーディンガー方程式（スピン軌道相互作用を平均化したもので、中重元素に必須）。
    • rel=2: 全相対論的ディラック方程式（スピン軌道分裂を完全に解く。フル相対論的擬ポテンシャル生成に使用）。
  • config: 原子全電子状態の初期電子配置を指定する文字列。たとえばケイ素（Si, $Z=14$）であれば '[Ne] 3s2 3p2' のように記述する。これにより、占有数の自動設定が行われる。

ネームリスト &inputp

  • pseudotype: 生成する擬ポテンシャルの数理形式を決定するフラグ。
    • pseudotype=1: 規格保存型（Norm-Conserving, NC）擬ポテンシャル。
    • pseudotype=2: 超ソフト型（Ultrasoft, USPP）擬ポテンシャル。
    • pseudotype=3: Projector Augmented Wave (PAW) 法用データセット。
  • file_pseudopotential: 生成された擬ポテンシャルを出力するファイル名。通常は拡張子 .UPF または .upf を指定する。
  • lloc: 局所ポテンシャル（Local Part）として採用する角運動量チャネル $l_{\text{loc}}$。通常は最外殻のチャネルのいずれか、または十分に高い $l$ チャネルを選択する。例えば lloc=-1 と指定すると、原子核近傍を滑らかにした動的ポテンシャル（超ソフト擬ポテンシャルなどの局所部分）が自動生成される。
  • rcore: 非線形コア補正（Non-linear Core Correction, NLCC）を適用する際のコアカットオフ半径 $r_c^{\text{core}}$。価電子密度の広がりとコア電子密度の広がりが重なる領域において、コア電子によるスピン極性や分極効果を局所密度近似の非線形性を通じて補正する。

構成カード（Configuration Cards）

ネームリスト &inputp の直後に、擬ポテンシャルのフィッティングおよびカットオフ半径をチャネルごとに詳細に規定する「価電子状態カード」を配置する。各行のパラメータは以下の構造を持つ：

n  l  occupation  energy  r_c  r_c_inner
  • n, l: 擬波動関数を作成する全電子軌道の主量子数 $n$ と方位量子数 $l$。
  • occupation: その軌道の基準占有数。
  • energy: 束縛状態ではない高エネルギー領域にプロジェクターを配置する場合に指定する参照エネルギー（通常は空白のままで、束縛固有値が自動的に用いられる）。
  • r_c: 全電子波動関数と擬波動関数をマッチングさせるカットオフ半径 $r_c^l$。これより外側で全電子状態と擬状態が一致する。
  • r_c_inner: USPP / PAW における補強電荷や局所波動関数のマッチングを行うための内側カットオフ半径 $r_{\text{in}}^l$。

3.2 擬ポテンシャル生成テンプレート

ここでは、実際に ld1.x で擬ポテンシャルを生成するための具体的な入力ファイルテンプレートを示す。

3.2.1 規格保存型 (Norm-Conserving) / GTH 擬ポテンシャル生成テンプレート

Troullier-Martins (TM) 法に基づく規格保存型擬ポテンシャルを生成するためのテンプレートである。非線形コア補正（NLCC）も有効にしている。

&input
    title='Silicon (Si) Norm-Conserving Generation',
    zed=14.0,
    rel=1,
    config='[Ne] 3s2 3p2',
    dft='PBE',
    isxt=3
/
&inputp
    pseudotype=1,
    file_pseudopotential='Si.pbe-nc.UPF',
    lloc=2,
    rcore=1.4,
    author='Antigravity'
/
3
3s  0  2.00  0.00  1.70  1.70
3p  1  2.00  0.00  1.70  1.70
3d  2  0.00  0.10  1.70  1.70

[!NOTE] ここで、末尾の 3 は擬ポテンシャルの生成チャネル数を表し、続く3行が各チャネルの定義（$3s, 3p, 3d$）である。$3d$ チャネルは空軌道であるが、分極効果や高エネルギー状態の記述能力（転移性）を向上させるために、局所的プロジェクターとして $1.70\,\text{a.u.}$ のカットオフで導入されている。

3.2.2 超ソフト (Ultrasoft) / PAW 擬ポテンシャル生成テンプレート

以下は、同じくケイ素（Si）に対して、Vanderbilt 形式の超ソフト擬ポテンシャル（USPP）を多重プロジェクター（各チャネル2個のプロジェクター）構成で生成するための入力テンプレートである。

&input
    title='Silicon (Si) Ultrasoft (USPP) Generation',
    zed=14.0,
    rel=1,
    config='[Ne] 3s2 3p2',
    dft='PBE'
/
&inputp
    pseudotype=2,
    file_pseudopotential='Si.pbe-us.UPF',
    lloc=-1,
    rcore=1.3,
    author='Antigravity'
/
4
3s  0  2.00  0.00  1.90  1.30
3s  0  0.00  0.80  1.90  1.30
3p  1  2.00  0.00  2.00  1.40
3p  1  0.00  0.80  2.00  1.40

[!TIP] 各チャネルに対して束縛状態（エネルギー 0.00）と高エネルギー状態（エネルギー 0.80 Ry）の2つのプロジェクターを定義している。マッチング半径 $r_c$ よりも小さい $r_{\text{in}}$（r_c_inner、それぞれ 1.30 および 1.40）を指定することで、コア領域における波動関数の柔軟性を確保し、カットオフエネルギー（平面波基底の数）を大幅に削減することが可能となる。

3.3 擬ポテンシャルの検証プロトコル（Verification Protocols）

生成された擬ポテンシャルは、計算の高速化のために導入された近似モデルであるため、実際のバルク結晶や分子のシミュレーションに投入する前に、その物理的妥当性と転移性を厳密に検証しなければならない。主要な検証手続きとして以下の2つが標準的である。

3.3.1 対数微分（Logarithmic Derivative）プロットとゴースト状態の検出

非局所擬ポテンシャルを導入した際、最も警戒すべき人工的欠陥がゴースト状態（Ghost States）である。ゴースト状態とは、非局所プロジェクターの Kleinman-Bylander (KB) 局所化形式への変換（ポテンシャルの非局所部 $\hat{V}_{\text{NL}}$ をセパラブルな形式へ展開する操作）に伴って発生する、全電子系には存在しない不自然な束縛状態（通常は価電子帯よりはるか下に位置する深い束縛状態、または価電子領域に現れる擬似的な励起状態）である。

これを検出するために、価電子の関与するエネルギー領域 $\epsilon$ において、原子核から十分に離れた一定のカットオフ半径 $R_{\text{cut}}$（通常は最外殻のマッチング半径の直上）における動的波動関数の対数微分（Logarithmic Derivative） $D_l(\epsilon)$ を計算する。

\[D_l(\epsilon) = \left. \frac{d}{dr} \ln \chi_l(r, \epsilon) \right\vert_{r=R_{\text{cut}}} = \left. \frac{\chi_l'(r, \epsilon)}{\chi_l(r, \epsilon)} \right\vert_{r=R_{\text{cut}}}\]

この $D_l(\epsilon)$ を、全電子（AE）計算と擬ポテンシャル（PS）計算の双方についてエネルギー $\epsilon$ の関数としてプロットし、比較する。

  1. 一致度の評価: 正確な擬ポテンシャルであれば、全電子の対数微分プロット $D_l^{\text{AE}}(\epsilon)$ と擬ポテンシャルの対数微分プロット $D_l^{\text{PS}}(\epsilon)$ は、価電子帯および低エネルギー伝導帯の領域（通常 $-2\,\text{Ry} < \epsilon < 2\,\text{Ry}$）で精密に重なり合う。
  2. ゴースト状態の同定: もし全電子計算には存在しないエネルギーにおいて、擬ポテンシャルの対数微分 $D_l^{\text{PS}}(\epsilon)$ に不連続な極（ポール）や異常な特異点が出現した場合、それは非局所ハミルトニアンの固有値が unphysical な状態（ゴースト）を作り出している決定的な証拠である。ゴースト状態が存在すると、バルク計算において全エネルギーが不自然に低下したり、SCFが収束しなくなったりする。
  3. 対策: ゴースト状態が検出された場合、局所ポテンシャルとして採用する角運動量チャネル $l_{\text{loc}}$ を変更するか、マッチング半径 $r_c^l$ を調整して擬ポテンシャルを再生成する必要がある。

3.3.2 遷移性（Transferability）チェック：全電子計算とのエネルギー差比較

擬ポテンシャルは特定の孤立原子の基準電子配置（Reference Configuration）において全電子結果を再現するように作られる。しかし、実際の計算対象（分子や結晶）の中では、原子は多様な化学結合やチャージ状態に置かれる。したがって、基準配置から大きく外れた異なる電子配置（励起状態やイオン化状態）に対しても、全電子計算と同様の応答を示すか（転移性）を検証しなければならない。

検証では、ld1.x の &test ネームリストを使用し、複数の電子配置（例えば、Si において $3s^2 3p^2$ の基底状態から、イオン化状態 $3s^2 3p^1 3d^1$、あるいは励起状態 $3s^1 3p^3$ など）に対して、全電子計算と擬ポテンシャル計算のそれぞれから全エネルギー $E_{\text{AE}}$ および $E_{\text{PS}}$ を算出する。 基準構成に対する各電子配置の励起エネルギー差（Configuration Energy Difference）

\[\Delta E_{\text{excitation}} = E(\text{Configuration}) - E(\text{Reference})\]

を計算し、全電子（AE）と擬ポテンシャル（PS）におけるその差のズレ（エラー）

\[\delta \Delta E = \Delta E_{\text{excitation}}^{\text{PS}} - \Delta E_{\text{excitation}}^{\text{AE}}\]

を評価する。この誤差 $\delta \Delta E$ が全てのテスト配置において数 $\text{mRy}$（約 $0.05\,\text{eV}$）以下に収まることが、高い転移性と信頼性を保証する基準となる。

結論

本稿では、密度汎関数理論（DFT）における擬ポテンシャル実装の深部を解説した。

  1. UPF v2 フォーマットの XML タグ体系と、原点付近の急峻なポテンシャル勾配を効率的に記述する対数径格子 r(i) および微分体積因子 rab(i) の関係を定式化した。
  2. Quantum ESPRESSO などの平面波コードの内部構造に踏み込み、実空間のプロジェクターを逆格子空間へ射影する init_us_1.f90 から、一般化オーバーラップ演算子 $S$ の作用を評価する s_psi.f90 に至るデータフローを解き明かした。
  3. 実務的な ld1.x を用いた擬ポテンシャル生成手順を示し、ゴースト状態を防ぐための対数微分プロットの重要性と、遷移性チェックによる精度検証プロトコルを確立した。

擬ポテンシャルは単なる平面波計算の「高速化ツール」ではなく、電子状態計算の精度の上限を規定する極めて重要な物理モデルである。その内部仕様とデータフローを精確に把握することは、高度な第一原理計算を実行し、あるいは新たな機能開発（DFT+U や非局所相関関数の実装など）を行う上で不可欠な土台である。

You May Also Enjoy