test

16 minute read

密度汎関数理論 (DFT) における交換相関汎関数の実装：実空間数値積分、libxc 連携、および Quantum ESPRESSO における実装解説

1. 実空間数値積分とスピン極性の定式化

1.1 実空間グリッド（FFTメッシュ）上での電子密度の離散化と交換相関エネルギーの離散和計算

第一原理密度汎関数理論（DFT）計算において、交換相関エネルギー $E_{xc}$ は、基底状態の電荷密度 $\rho(\mathbf{r})$ に依存する空間積分として以下のように定義される。 \(E_{xc}[\rho] = \int \rho(\mathbf{r}) \varepsilon_{xc}(\rho(\mathbf{r})) d\mathbf{r}\) ここで、一般化勾配近似（GGA）ではエネルギー密度 $\varepsilon_{xc}$ が局所電子密度 $\rho(\mathbf{r})$ に加えてその勾配の絶対値 $|\nabla\rho(\mathbf{r})|$ にも依存し、meta-GGAではさらに局所運動エネルギー密度 $\tau(\mathbf{r})$ にも依存する。

実際の計算機上での実装、特に平面波ポテンシャル法を用いるシステムでは、この空間積分を実空間上に定義された離散的なグリッド（FFTメッシュ）上で実行する。平面波基底を採用するコードでは、波動関数や電子密度は逆空間（$\mathbf{G}$ ベクトル空間）で表されるが、交換相関ポテンシャル $V_{xc}$ や外部ポテンシャル $V_{\text{ext}}$ などの局所ポテンシャルは実空間において対角化される。この双対性を高速に処理するため、高速フーリエ変換（FFT）を利用して二つの空間を行き来する。したがって、実空間上での離散的な積分点は、FFTを実行するメッシュと一致させることが最も効率的である。

単位胞の体積を $\Omega$、実空間グリッド点の総数を $N_g = N_{g1} \times N_{g2} \times N_{g3}$ とし、各グリッド点の座標を $\mathbf{r}p$ とすると、連続空間での積分は以下のように離散和（数値積分）へと近似される。 \(E_{xc} \approx \sum_{p=1}^{N_g} w_p \rho(\mathbf{r}_p) \varepsilon_{xc}\left(\rho(\mathbf{r}_p), \nabla\rho(\mathbf{r}_p), \tau(\mathbf{r}_p), \dots\right)\) ここで、$w_p$ は各グリッド点における積分の数値的重み（体積要素）である。一般的に用いられる等間隔な直交メッシュにおいては、各グリッド点が占める体積は一様であり、以下のように単純なセルの総体積と格子点数の比となる。 \(w_p = \frac{\Omega}{N_g}\) このようにして、複雑な関数形を持つ交換相関エネルギー密度 $\varepsilon{xc}$ の空間積分が、FFTメッシュ上の電子密度 $\rho(\mathbf{r}_p)$ およびその微分量を用いた単純な離散和に還元され、高速な並列計算が可能となる。

1.2 コロニアルスピン極性（LSDA）と非コロニアルスピン（Non-collinear spin）の処理

スピン磁性やスピン軌道相互作用を扱う際、スピン分極の記述方法には「コロニアル（collinear）」と「非コロニアル（non-collinear）」の二通りのアプローチが存在する。

1. コロニアルスピン極性（LSDA）

通常の局所スピン密度近似（LSDA）に代表されるコロニアルスピン系では、系全体でスピン量子化軸が一方向（通常は $z$ 軸）に固定されていると仮定する。このとき、電子密度は空間の各点 $\mathbf{r}$ において、上向きスピン密度 $\rho_\uparrow(\mathbf{r})$ と下向きスピン密度 $\rho_\downarrow(\mathbf{r})$ の二つの独立な実数スカラー量として記述される。

  • 全電荷密度: $\rho(\mathbf{r}) = \rho_\uparrow(\mathbf{r}) + \rho_\downarrow(\mathbf{r})$
  • 局所磁化（スピン極性）: $m(\mathbf{r}) = \rho_\uparrow(\mathbf{r}) - \rho_\downarrow(\mathbf{r})$

この記述においては、各スピンに対する交換相関ポテンシャル $V_{xc,\sigma}(\mathbf{r})$ （$\sigma \in {\uparrow, \downarrow}$）が以下のように独立に評価される。 \(V_{xc,\uparrow}(\mathbf{r}) = \frac{\delta E_{xc}}{\delta \rho_\uparrow(\mathbf{r})}, \quad V_{xc,\downarrow}(\mathbf{r}) = \frac{\delta E_{xc}}{\delta \rho_\downarrow(\mathbf{r})}\)

2. 非コロニアルスピン（Non-collinear spin）

スピン軌道相互作用（Spin-Orbit Coupling: SOC）が無視できない重元素を含む系や、スピンの向きが空間的に連続的に傾く非共線磁性体を扱う場合、量子化軸を一方向に固定するアプローチは破綻する。このような非コロニアル系では、各空間点 $\mathbf{r}$ におけるスピン状態を、スピン空間における $2 \times 2$ エルミート密度行列 $\rho_{\alpha\beta}(\mathbf{r})$ を用いて表現する。 \(\rho_{\alpha\beta}(\mathbf{r}) = \frac{1}{2} \left[ \rho(\mathbf{r})\delta_{\alpha\beta} + \mathbf{m}(\mathbf{r}) \cdot \boldsymbol{\sigma}_{\alpha\beta} \right]\) ここで、

  • $\alpha, \beta \in {\uparrow, \downarrow}$ はスピンの添字である。
  • $\rho(\mathbf{r})$ は全電荷密度であり、密度行列のトレース（対角和）に対応する。 \(\rho(\mathbf{r}) = \text{Tr}[\rho_{\alpha\beta}(\mathbf{r})] = \rho_{\uparrow\uparrow}(\mathbf{r}) + \rho_{\downarrow\downarrow}(\mathbf{r})\)
  • $\mathbf{m}(\mathbf{r}) = \left(m_x(\mathbf{r}), m_y(\mathbf{r}), m_z(\mathbf{r})\right)$ は局所磁化ベクトルである。
  • $\boldsymbol{\sigma} = (\sigma_x, \sigma_y, \sigma_z)$ はスピン空間におけるパウリ行列ベクトルである。 \(\sigma_x = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}, \quad \sigma_y = \begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix}, \quad \sigma_z = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}\)

密度行列 $\rho_{\alpha\beta}(\mathbf{r})$ を成分ごとに具体的に展開すると、以下のようになる。 \(\rho_{\alpha\beta}(\mathbf{r}) = \frac{1}{2} \begin{pmatrix} \rho(\mathbf{r}) + m_z(\mathbf{r}) & m_x(\mathbf{r}) - i m_y(\mathbf{r}) \\ m_x(\mathbf{r}) + i m_y(\mathbf{r}) & \rho(\mathbf{r}) - m_z(\mathbf{r}) \end{pmatrix}\) 局所磁化ベクトル $\mathbf{m}(\mathbf{r})$ の各成分は、密度行列とパウリ行列の積のトレースから次のように抽出される。 \(m_i(\mathbf{r}) = \text{Tr}[\rho_{\alpha\beta}(\mathbf{r}) \sigma_i] = \sum_{\alpha\beta} \rho_{\alpha\beta}(\mathbf{r}) \sigma_{i,\beta\alpha} \quad (i = x, y, z)\) すなわち、 \(m_x(\mathbf{r}) = 2\,\text{Re}[\rho_{\downarrow\uparrow}(\mathbf{r})]\) \(m_y(\mathbf{r}) = 2\,\text{Im}[\rho_{\downarrow\uparrow}(\mathbf{r})]\) \(m_z(\mathbf{r}) = \rho_{\uparrow\uparrow}(\mathbf{r}) - \rho_{\downarrow\downarrow}(\mathbf{r})\)

1.3 密度行列の対角化と局所固有値（局所スピン密度）

非コロニアルスピン計算において、通常の交換相関汎関数（コロニアルスピン用に定式化されているLSDAやGGAなど）の関数形をそのまま利用するため、各実空間グリッド点 $\mathbf{r}$ において $2 \times 2$ 密度行列 $\rho_{\alpha\beta}(\mathbf{r})$ の局所的な対角化を行う。

固有値 $\lambda$ を求めるための特性方程式 $\det\left(\rho_{\alpha\beta}(\mathbf{r}) - \lambda \delta_{\alpha\beta}\right) = 0$ を解く。 \(\det \begin{pmatrix} \frac{1}{2}(\rho + m_z) - \lambda & \frac{1}{2}(m_x - i m_y) \\ \frac{1}{2}(m_x + i m_y) & \frac{1}{2}(\rho - m_z) - \lambda \end{pmatrix} = 0\) \(\left( \frac{1}{2}\rho - \lambda + \frac{1}{2}m_z \right) \left( \frac{1}{2}\rho - \lambda - \frac{1}{2}m_z \right) - \frac{1}{4}(m_x - i m_y)(m_x + i m_y) = 0\) \(\left( \lambda - \frac{1}{2}\rho \right)^2 - \frac{1}{4}m_z^2 - \frac{1}{4}(m_x^2 + m_y^2) = 0\) \(\left( \lambda - \frac{1}{2}\rho \right)^2 = \frac{1}{4}\left( m_x^2 + m_y^2 + m_z^2 \right) = \frac{1}{4}|\mathbf{m}(\mathbf{r})|^2\) この方程式を解くことで、固有値としての局所スピン密度 $\rho_{\uparrow,\downarrow}(\mathbf{r})$ が得られる。 \(\rho_{\uparrow,\downarrow}(\mathbf{r}) = \frac{1}{2} \left[ \rho(\mathbf{r}) \pm |\mathbf{m}(\mathbf{r})| \right]\) ここで $|\mathbf{m}(\mathbf{r})| = \sqrt{m_x(\mathbf{r})^2 + m_y(\mathbf{r})^2 + m_z(\mathbf{r})^2}$ である。 この固有値は、各空間点において局所的に定義される量子化軸（局所磁化ベクトル $\mathbf{m}(\mathbf{r})$ の方向）に対するスピンアップおよびスピンダウンの電子密度に対応している。

これにより、非コロニアル系における交換相関エネルギー $E_{xc}$ は、局所対角化されたスピン密度を用いて以下のように定義される。 \(E_{xc} = \int \rho(\mathbf{r}) \varepsilon_{xc}\left(\rho_\uparrow(\mathbf{r}), \rho_\downarrow(\mathbf{r}), \nabla\rho_\uparrow(\mathbf{r}), \nabla\rho_\downarrow(\mathbf{r}), \dots\right) d\mathbf{r}\) この枠組みを導入することで、非コロニアルスピン軌道計算であっても、既存のコロニアル用交換相関汎関数パッケージ（libxc など）を各点での対角化を通じて完全に再利用することが可能となる。

1.4 交換相関ポテンシャルおよび有効磁場 $\mathbf{B}_{xc}(\mathbf{r})$ への逆変換

エネルギーを計算した後、Kohn-Sham方程式を自己無撞着に解くためには、実空間における交換相関ポテンシャル行列 $V_{xc,\alpha\beta}(\mathbf{r})$ が必要となる。これはエネルギーの密度行列成分に対する関数微分として定義される。 \(V_{xc,\alpha\beta}(\mathbf{r}) = \frac{\delta E_{xc}}{\delta \rho_{\beta\alpha}(\mathbf{r})}\) この $2 \times 2$ ポテンシャル行列は、電荷密度に対するスカラー交換相関ポテンシャル $V_{xc}(\mathbf{r})$ と、スピン磁化ベクトルと結合する有効磁場ベクトル（ベクトルポテンシャル） $\mathbf{B}_{xc}(\mathbf{r})$ を用いて次のように展開できる。 \(V_{xc,\alpha\beta}(\mathbf{r}) = V_{xc}(\mathbf{r})\delta_{\alpha\beta} - \mathbf{B}_{xc}(\mathbf{r}) \cdot \boldsymbol{\sigma}_{\alpha\beta}\)

この展開式におけるスカラーポテンシャル $V_{xc}(\mathbf{r})$ および有効磁場 $\mathbf{B}{xc}(\mathbf{r})$ を、コロニアル用汎関数から得られるスピンごとのポテンシャル $V{xc,\uparrow}(\mathbf{r}) = \frac{\delta E_{xc}}{\delta \rho_\uparrow(\mathbf{r})}$ および $V_{xc,\downarrow}(\mathbf{r}) = \frac{\delta E_{xc}}{\delta \rho_\downarrow(\mathbf{r})}$ から逆変換によって導出する。

合成関数の関数微分法（チェインルール）を用いると、密度行列の共役成分に対する微分は以下のように表される。 \(V_{xc,\alpha\beta} = \frac{\delta E_{xc}}{\delta \rho_{\beta\alpha}} = \frac{\delta E_{xc}}{\delta \rho_\uparrow} \frac{\partial \rho_\uparrow}{\partial \rho_{\beta\alpha}} + \frac{\delta E_{xc}}{\delta \rho_\downarrow} \frac{\partial \rho_\downarrow}{\partial \rho_{\beta\alpha}} = V_{xc,\uparrow} \frac{\partial \rho_\uparrow}{\partial \rho_{\beta\alpha}} + V_{xc,\downarrow} \frac{\partial \rho_\downarrow}{\partial \rho_{\beta\alpha}}\) ここで、局所固有値の定義 $\rho_{\uparrow,\downarrow} = \frac{1}{2}(\rho \pm |\mathbf{m}|)$ を用いて、密度行列 $\rho_{\beta\alpha}$ に対する偏微分を計算する。電荷密度 $\rho = \text{Tr}(\rho) = \sum_\gamma \rho_{\gamma\gamma}$ および磁化ベクトル成分 $m_i = \text{Tr}(\rho \sigma_i) = \sum_{\gamma\delta} \rho_{\gamma\delta}\sigma_{i,\delta\gamma}$ より、 \(\frac{\partial \rho}{\partial \rho_{\beta\alpha}} = \delta_{\alpha\beta}\) \(\frac{\partial m_i}{\partial \rho_{\beta\alpha}} = \sigma_{i,\alpha\beta}\) これらを用いて磁化の絶対値 $|\mathbf{m}| = \sqrt{\sum_i m_i^2}$ の微分を求めると（$|\mathbf{m}| \neq 0$ とする）： \(\frac{\partial |\mathbf{m}|}{\partial \rho_{\beta\alpha}} = \sum_i \frac{\partial |\mathbf{m}|}{\partial m_i} \frac{\partial m_i}{\partial \rho_{\beta\alpha}} = \sum_i \frac{m_i}{|\mathbf{m}|} \sigma_{i,\alpha\beta} = \frac{\mathbf{m}}{|\mathbf{m}|} \cdot \boldsymbol{\sigma}_{\alpha\beta}\) したがって、局所スピン密度 $\rho_\uparrow, \rho_\downarrow$ の密度行列要素に対する微分は以下のようになる。 \(\frac{\partial \rho_\uparrow}{\partial \rho_{\beta\alpha}} = \frac{1}{2} \left( \delta_{\alpha\beta} + \frac{\mathbf{m}}{|\mathbf{m}|} \cdot \boldsymbol{\sigma}_{\alpha\beta} \right)\) \(\frac{\partial \rho_\downarrow}{\partial \rho_{\beta\alpha}} = \frac{1}{2} \left( \delta_{\alpha\beta} - \frac{\mathbf{m}}{|\mathbf{m}|} \cdot \boldsymbol{\sigma}_{\alpha\beta} \right)\)

これらを $V_{xc,\alpha\beta}$ のチェインルール式に代入して整理する。 \(\begin{aligned} V_{xc,\alpha\beta} &= V_{xc,\uparrow} \left[ \frac{1}{2} \left( \delta_{\alpha\beta} + \frac{\mathbf{m}}{|\mathbf{m}|} \cdot \boldsymbol{\sigma}_{\alpha\beta} \right) \right] + V_{xc,\downarrow} \left[ \frac{1}{2} \left( \delta_{\alpha\beta} - \frac{\mathbf{m}}{|\mathbf{m}|} \cdot \boldsymbol{\sigma}_{\alpha\beta} \right) \right] \\ &= \frac{1}{2}(V_{xc,\uparrow} + V_{xc,\downarrow})\delta_{\alpha\beta} + \frac{1}{2}(V_{xc,\uparrow} - V_{xc,\downarrow}) \frac{\mathbf{m}}{|\mathbf{m}|} \cdot \boldsymbol{\sigma}_{\alpha\beta} \end{aligned}\) この結果を、ポテンシャル行列の展開式 $V_{xc,\alpha\beta}(\mathbf{r}) = V_{xc}(\mathbf{r})\delta_{\alpha\beta} - \mathbf{B}{xc}(\mathbf{r}) \cdot \boldsymbol{\sigma}{\alpha\beta}$ と比較することで、スカラー交換相関ポテンシャル $V_{xc}(\mathbf{r})$ および有効磁場 $\mathbf{B}_{xc}(\mathbf{r})$ の具体的な逆変換式が求まる。 \(V_{xc}(\mathbf{r}) = \frac{1}{2} \left( V_{xc,\uparrow}(\mathbf{r}) + V_{xc,\downarrow}(\mathbf{r}) \right)\) \(\mathbf{B}_{xc}(\mathbf{r}) = - \frac{1}{2} \left( V_{xc,\uparrow}(\mathbf{r}) - V_{xc,\downarrow}(\mathbf{r}) \right) \frac{\mathbf{m}(\mathbf{r})}{|\mathbf{m}(\mathbf{r})|}\)

この定式化は、有効磁場 $\mathbf{B}{xc}(\mathbf{r})$ が局所磁化ベクトル $\mathbf{m}(\mathbf{r})$ と平行であり、その強さは局所的な上向き・下向きスピンに対する交換相関ポテンシャルの差に等しいことを意味する。実空間メッシュの各格子点において得られた $V{xc}(\mathbf{r})$ および $\mathbf{B}{xc}(\mathbf{r})$ から $2 \times 2$ のポテンシャル行列 $V{xc,\alpha\beta}(\mathbf{r})$ が再構成され、SCF計算におけるハミルトニアンの更新に供される。

2. 一般化勾配近似（GGA）のポテンシャル変分と空間微分の離散化

2.1 GGAポテンシャル微分の導出と変分法（部分積分）の適用

局所（スピン）密度近似（LDA/LSDA）では、交換相関エネルギー密度は各点における電子密度 $\rho(\mathbf{r})$ の値のみに依存していた。一方、一般化勾配近似（Generalized Gradient Approximation: GGA）では、電子密度の不均一性を考慮するために、エネルギー密度は電子密度 $\rho(\mathbf{r})$ とその空間勾配の大きさ $|\nabla \rho(\mathbf{r})|$ にも依存する。

スピン分極のない（閉殻）系におけるGGA交換相関エネルギー $E_{xc}^{GGA}[\rho]$ は、以下のように定式化される。 \(E_{xc}^{GGA}[\rho] = \int f(\rho(\mathbf{r}), \sigma(\mathbf{r})) d\mathbf{r}\) ここで、微分演算の都合上、勾配の絶対値そのものではなく、その2乗として定義されるスカラー量 $\sigma(\mathbf{r})$ を独立な変数として導入することが一般的である。 \(\sigma(\mathbf{r}) = |\nabla \rho(\mathbf{r})|^2 = \nabla \rho(\mathbf{r}) \cdot \nabla \rho(\mathbf{r})\) Kohn-Sham 方程式において自己無撞着（SCF）計算を行うには、交換相関エネルギーの電子密度に対する関数微分（Kohn-Sham ポテンシャル） $V_{xc}(\mathbf{r})$ を求める必要がある。 \(V_{xc}(\mathbf{r}) = \frac{\delta E_{xc}^{GGA}}{\delta \rho(\mathbf{r})}\) この関数微分を導出するため、電子密度の微小変化 $\delta \rho(\mathbf{r})$ に対するエネルギーの変分 $\delta E_{xc}^{GGA}$ を考える。 \(\delta E_{xc}^{GGA} = E_{xc}^{GGA}[\rho + \delta\rho] - E_{xc}^{GGA}[\rho] = \int \left[ \frac{\partial f}{\partial \rho} \delta\rho + \frac{\partial f}{\partial \sigma} \delta\sigma \right] d\mathbf{r}\) ここで、$\sigma$ の変分 $\delta\sigma$ は、勾配演算子 $\nabla$ と変分演算子 $\delta$ が交換可能であること（$\delta(\nabla\rho) = \nabla(\delta\rho)$）を利用すると、次のように展開できる。 \(\delta\sigma = \delta(\nabla \rho \cdot \nabla \rho) = 2 \nabla \rho \cdot \delta(\nabla \rho) = 2 \nabla \rho \cdot \nabla(\delta\rho)\) これを変分式に代入する。 \(\delta E_{xc}^{GGA} = \int \left[ \frac{\partial f}{\partial \rho} \delta\rho + 2 \frac{\partial f}{\partial \sigma} \nabla \rho \cdot \nabla(\delta\rho) \right] d\mathbf{r}\) 第2項の被積分関数には変分 $\delta\rho$ の勾配 $\nabla(\delta\rho)$ が含まれているため、このままでは密度に対する関数微分の形（$\int V_{xc}(\mathbf{r}) \delta\rho(\mathbf{r}) d\mathbf{r}$）を直接抽出することができない。そこで、空間的な部分積分（またはベクトルの発散定理）を用いて、勾配演算子を $\delta\rho$ から外す操作を行う。

ここで、任意のベクトル場 $\mathbf{A}(\mathbf{r})$ とスカラー場 $\phi(\mathbf{r})$ に対するベクトル解析の公式： \(\nabla \cdot (\phi \mathbf{A}) = \nabla \phi \cdot \mathbf{A} + \phi (\nabla \cdot \mathbf{A})\) を適用する。ここで $\phi = \delta\rho$、$\mathbf{A} = 2 \frac{\partial f}{\partial \sigma} \nabla\rho$ と置くと、 \(\nabla \cdot \left( \delta\rho \cdot 2 \frac{\partial f}{\partial \sigma} \nabla\rho \right) = \nabla(\delta\rho) \cdot \left( 2 \frac{\partial f}{\partial \sigma} \nabla\rho \right) + \delta\rho \left[ \nabla \cdot \left( 2 \frac{\partial f}{\partial \sigma} \nabla\rho \right) \right]\) となる。この式を整理して、変分式第2項の被積分関数を置き換えると以下のようになる。 \(2 \frac{\partial f}{\partial \sigma} \nabla \rho \cdot \nabla(\delta\rho) = \nabla \cdot \left( 2 \frac{\partial f}{\partial \sigma} \nabla\rho \, \delta\rho \right) - \left[ \nabla \cdot \left( 2 \frac{\partial f}{\partial \sigma} \nabla\rho \right) \right] \delta\rho\) これをエネルギー変分式に代入すると、次の表現を得る。 \(\delta E_{xc}^{GGA} = \int \left[ \frac{\partial f}{\partial \rho} - \nabla \cdot \left( 2 \frac{\partial f}{\partial \sigma} \nabla\rho \right) \right] \delta\rho \, d\mathbf{r} + \int \nabla \cdot \left( 2 \frac{\partial f}{\partial \sigma} \nabla\rho \, \delta\rho \right) d\mathbf{r}\) ここで、右辺第2項の全発散項（Total Divergence Term）にガウスの発散定理を適用する。積分領域を単位胞（セル）の体積 $\Omega$ とし、その境界を $\partial \Omega$、境界上の外向き法線ベクトルを $\mathbf{n}$ とすると、 \(\int_{\Omega} \nabla \cdot \left( 2 \frac{\partial f}{\partial \sigma} \nabla\rho \, \delta\rho \right) d\mathbf{r} = \oint_{\partial \Omega} \left( 2 \frac{\partial f}{\partial \sigma} \nabla\rho \, \delta\rho \right) \cdot \mathbf{n} dS\) となる。 周期境界条件（Periodic Boundary Conditions: PBC）が課されるバルク結晶計算においては、境界 $\partial \Omega$ の対向する面同士で物理量（密度、勾配、および変分）が完全に一致するため、面積分は互いに相殺してゼロとなる。また、分子やクラスターなどの孤立系を扱う場合であっても、系の境界を無限遠にとれば、境界上での電子密度およびその変分 $\delta\rho$ は指数関数的にゼロへ減衰するため、この表面項はやはり無視できる。

したがって、境界項が完全に消失し、変分式は以下のように整理される。 \(\delta E_{xc}^{GGA} = \int \left[ \frac{\partial f}{\partial \rho} - \nabla \cdot \left( 2 \frac{\partial f}{\partial \sigma} \nabla\rho \right) \right] \delta\rho(\mathbf{r}) d\mathbf{r}\) 変分 $\delta\rho(\mathbf{r})$ は任意であるため、被積分関数の括弧内が求める関数微分そのものとなり、GGAのKohn-Sham交換相関ポテンシャル $V_{xc}(\mathbf{r})$ が得られる。 \(V_{xc}(\mathbf{r}) = \frac{\partial f}{\partial \rho} - \nabla \cdot \left( 2 \frac{\partial f}{\partial \sigma} \nabla \rho \right)\)

スピン極性（LSDA）を考慮したGGA（一般に $f(\rho_\uparrow, \rho_\downarrow, \sigma_{\uparrow\uparrow}, \sigma_{\uparrow\downarrow}, \sigma_{\downarrow\downarrow})$ で定義される）の場合、勾配変数として $\sigma_{\uparrow\uparrow} = |\nabla\rho_\uparrow|^2$、$\sigma_{\downarrow\downarrow} = |\nabla\rho_\downarrow|^2$、$\sigma_{\uparrow\downarrow} = \nabla\rho_\uparrow \cdot \nabla\rho_\downarrow$ の3つが定義され、それぞれのスピンに対するポテンシャルは同様の変分操作と部分積分によって以下のように一般化される。 \(V_{xc,\uparrow}(\mathbf{r}) = \frac{\partial f}{\partial \rho_\uparrow} - \nabla \cdot \left( 2 \frac{\partial f}{\partial \sigma_{\uparrow\uparrow}} \nabla\rho_\uparrow + \frac{\partial f}{\partial \sigma_{\uparrow\downarrow}} \nabla\rho_\downarrow \right)\) \(V_{xc,\downarrow}(\mathbf{r}) = \frac{\partial f}{\partial \rho_\downarrow} - \nabla \cdot \left( 2 \frac{\partial f}{\partial \sigma_{\downarrow\downarrow}} \nabla\rho_\downarrow + \frac{\partial f}{\partial \sigma_{\uparrow\downarrow}} \nabla\rho_\uparrow \right)\) この導出から明らかなように、GGAポテンシャルの計算には、エネルギー密度関数の偏微分（$\partial f / \partial \rho$, $\partial f / \partial \sigma$）に加えて、電子密度の勾配 $\nabla\rho$、さらにはその勾配ベクトルに偏微分係数を掛けたベクトルに対する発散（Divergence）操作 $\nabla \cdot (\dots)$ という、2階の空間微分が必要となる。これが、実空間グリッド上での実装において数値的な複雑さと不安定性を引き起こす主要な原因である。

2.2 実空間グリッドにおける勾配・発散演算の離散化スキーム

実空間のFFTグリッド上でGGAポテンシャルを計算する際、連続空間での勾配 $\nabla$ および発散 $\nabla \cdot$ の演算を離散化する必要がある。主に用いられるスキームとして、「実空間中心差分法（Central Finite Difference Method）」と「波数空間微分法（Fourier-space Derivative Method）」の2種類が存在し、それぞれ異なる利点と欠点を持っている。

1. 実空間中心差分法（Central Finite Difference Method）

実空間中心差分法では、各グリッド点における微分値を、その点の周囲の近傍グリッド点における密度の値の差分（線形結合）として近似する。 例えば、格子間隔 $h$ の1次元一様グリッドにおいて、点 $x_i$ における1階微分の2次精度中心差分は以下のように表される。 \(\left(\frac{\partial g}{\partial x}\right)_i \approx \frac{g(x_{i+1}) - g(x_{i-1})}{2h}\) 精度を向上させるため、より広い範囲（ステンシル）の点を用いる高次差分スキームが広く使用される。例えば、4次精度の中心差分は以下の形をとる。 \(\left(\frac{\partial g}{\partial x}\right)_i \approx \frac{-g(x_{i+2}) + 8g(x_{i+1}) - 8g(x_{i-1}) + g(x_{i-2})}{12h}\) さらに高次（6次、8次など）の差分を用いることで、高周波のカットオフ誤差を減らすことができる。3次元空間においては、各デカルト座標軸方向（$x, y, z$）に対してこの差分演算を独立に適用し、勾配ベクトル $\nabla\rho$ を構成する。発散演算 $\nabla \cdot \mathbf{A}$ に関しても、ベクトル場 $\mathbf{A} = 2 \frac{\partial f}{\partial \sigma} \nabla\rho$ の各成分 $A_x, A_y, A_z$ に対して同様の差分処理を行い、それらを足し合わせる。

  • 境界条件の処理: 孤立系（Dirichlet/Neumann 境界条件）および周期境界条件の双方に対して、境界点でのステンシルの回り込み処理を定義するだけで容易に実装できる。特に孤立系においてセル端部で関数値が完全にゼロになる場合、境界外の仮想点（Ghost points）の値をゼロとすることで自然に処理可能である。
  • 計算効率と並列化: 非常に局所的なメモリ呼び出しのみで完結するため、キャッシュ効率が極めて高く、メモリバンド幅がボトルネックになりにくい。また、領域分割（Domain Decomposition）による並列化を行う場合、隣接プロセス間での「ハロー領域（Halo/Ghost cells）」の通信（近接通信）のみで微分演算を実行できるため、MPI並列化におけるスケーラビリティが非常に優れている。
  • 安定性と数値ノイズ: 一方で、有限差分法は高周波のノイズ（グリッドレベルで振動する aliasing ノイズ）に弱い。特に、$\partial f/\partial \sigma$ などの汎関数微分は密度の変化に対して非常に非線形な応答を示すため、電子密度に微小なグリッドノイズが存在すると、2階微分（発散演算）の段階でノイズが増幅され、自己無撞着（SCF）計算の収束が著しく悪化することがある。このため、差分を適用する前に電子密度に平滑化フィルタを適用したり、高次のソフトニング技術を組み合わせたりする工夫が必要となる。

2. 波数空間微分法（Fourier-space Derivative Method）

周期境界条件を前提とする計算では、高速フーリエ変換（FFT）を利用して、すべての空間微分演算を reciprocal 空間（波数空間、逆空間）で実行する手法が極めて自然かつ強力である。 周期的な関数 $g(\mathbf{r})$ は、逆格子ベクトル $\mathbf{G}$ を用いて以下のようにフーリエ展開される。 \(g(\mathbf{r}) = \sum_{\mathbf{G}} \tilde{g}(\mathbf{G}) e^{i \mathbf{G} \cdot \mathbf{r}}\) このとき、空間微分（勾配 $\nabla$）は、波数空間における単純な代数積 $i\mathbf{G}$ に対応する。 \(\nabla g(\mathbf{r}) = \sum_{\mathbf{G}} i\mathbf{G} \tilde{g}(\mathbf{G}) e^{i \mathbf{G} \cdot \mathbf{r}} = \text{FFT}^{-1} \left[ i\mathbf{G} \cdot \text{FFT}[g(\mathbf{r})] \right]\) 具体的には、実空間での電子密度 $\rho(\mathbf{r})$ を前向きFFTで波数空間へ写し、$\tilde{\rho}(\mathbf{G})$ に $i\mathbf{G}$ の各成分（$i G_x, i G_y, i G_z$）を乗じた後、3回の逆FFTを実行することで、実空間での正確な勾配成分 $\nabla\rho$ を得る。 同様に、発散演算 $\nabla \cdot \mathbf{A}(\mathbf{r})$ も以下のように計算される。 \(\nabla \cdot \mathbf{A}(\mathbf{r}) = \text{FFT}^{-1} \left[ i\mathbf{G} \cdot \text{FFT}[\mathbf{A}(\mathbf{r})] \right]\) ここで $\text{FFT}[\mathbf{A}(\mathbf{r})]$ はベクトル場 $\mathbf{A}$ の各成分のフーリエ変換であり、波数空間における内積 $i\mathbf{G} \cdot \tilde{\mathbf{A}}(\mathbf{G})$ を計算した後に、1回の逆FFTによって実空間に戻す。

  • 計算精度: 波数空間での $i\mathbf{G}$ による乗算は、グリッドカットオフ内のすべての波長成分に対して数学的に厳密であり、有限差分法のような「空間の離散化に起因する異方性誤差（Grid Anisotropy）」や「高次項の打ち切り誤差」が一切発生しない（実質的に無限次精度の差分に相当）。
  • 境界条件の処理: 空間が3次元的に完全に周期的であるシステム（バルク結晶など）に限定される。孤立系（分子など）を扱う場合は、十分な大きさの真空領域（スーパーセル）を用意し、必要に応じてPoissonソルバーやクーロンカットオフ法を用いて、隣接セル間の人工的な像相互作用を遮蔽する必要がある。
  • 計算効率と並列化: 微分ごとにFFTを実行する必要があるため、計算量は $\mathcal{O}(N_g \log N_g)$ となる。勾配の計算には1回の正FFTと3回の逆FFT（計4回）が必要であり、発散の計算には3回の正FFTと1回の逆FFT（計4回）が必要となるため、実空間の差分法と比較して演算コストが高い。さらに、大規模な並列計算環境においては、3D FFT を実行する際に全プロセス間での大規模なデータ再配分（MPI All-to-all 通信）が発生するため、数万コア規模の超並列計算において通信オーバーヘッドがボトルネックになりやすい。
  • 安定性と数値ノイズ: 原理的に aliasing エラーが回避されるため、滑らかな電荷密度に対しては非常に安定した結果を与える。しかし、電子密度が極めて急峻に変化する領域（例えば、擬ポテンシャルを用いない全電子計算において原子核の近傍など）では、カットオフ制限によりギブス現象（Gibbs oscillation）による数値的な波打ちが生じ、これが非線形な汎関数微分を通じてゴースト状態やSCFの不安定化を招く場合がある。
項目 実空間中心差分法 波数空間微分法 (FFT)
精度 ステンシル数に応じた有限精度（切り捨て誤差あり） カットオフ内で厳密（無限次精度、方向依存性なし）
境界条件 周期境界・孤立系（開境界）の双方に容易に対応 周期境界条件に限定（孤立系にはスーパーセルが必要）
計算コスト (単一ノード) 低い（局所アクセスのみ、$\mathcal{O}(N_g)$） 高い（複数回の3D FFTが必要、$\mathcal{O}(N_g \log N_g)$）
並列スケーラビリティ 高い（ハロー領域の局所通信のみ） 低い（FFTに伴う MPI_Alltoall グローバル通信がボトルネック）
高周波ノイズへの耐性 弱い（高階差分でノイズが増幅されやすく、平滑化が必要） 強い（aliasingはないが、急峻な変化に対してGibbs振動が発生）

実際のDFT実装（例えば Quantum ESPRESSO や VASP など）においては、平面波基底を採用している都合上、すでにFFTグリッドと波数空間のデータ構造が最適化されているため、多くのコードで波数空間（FFT）を利用した微分がデフォルトで採用されている。しかし、実空間グリッド法（Real-space Grid DFT）や、大規模並列化を極限まで追求する実装においては、通信ボトルネックを回避するために、高次の実空間中心差分法が好んで選択される。

3. meta-GGA汎関数の定式化と非局所運動エネルギー演算子の導出

3.1 meta-GGAにおける軌道変分とポテンシャル演算子の定義

局所密度近似（LDA）や一般化勾配近似（GGA）では、交換相関エネルギー $E_{xc}$ は電子密度 $\rho(\mathbf{r})$ およびその勾配 $\nabla\rho(\mathbf{r})$ にのみ依存する、いわゆる局所または半局所汎関数であった。これに対し、meta-GGA汎関数では、局所的な運動エネルギー密度（Kinetic Energy Density） $\tau(\mathbf{r})$ が独立な変数として追加される。 Kohn-Sham軌道 $\phi_i(\mathbf{r})$ を用いると、運動エネルギー密度 $\tau(\mathbf{r})$ は以下のように定義される。 \(\tau(\mathbf{r}) = \frac{1}{2} \sum_{i} f_i |\nabla \phi_i(\mathbf{r})|^2 = \frac{1}{2} \sum_{i} f_i \nabla \phi_i^*(\mathbf{r}) \cdot \nabla \phi_i(\mathbf{r})\) ここで、$f_i$ は軌道 $i$ の占有数であり、和はすべてのKohn-Sham軌道にわたって取られる（簡略化のため、以下では $f_i = 1$ の閉殻系、または $f_i$ を波動関数 $\phi_i$ に内包させて表記する）。

meta-GGAにおける交換相関エネルギーは、以下のように $\rho$, $\sigma = |\nabla \rho|^2$, および $\tau$ の関数として定式化される。 \(E_{xc}^{mGGA}[\{\phi_i\}] = \int f(\rho(\mathbf{r}), \sigma(\mathbf{r}), \tau(\mathbf{r})) d\mathbf{r}\) ここで重要なのは、エネルギー $E_{xc}^{mGGA}$ が、$\tau(\mathbf{r})$ を通じて Kohn-Sham 軌道 $\phi_i(\mathbf{r})$ に直接（陽に）依存している点である。したがって、自己無撞着（SCF）計算を可能にするための交換相関ポテンシャルは、電子密度に対する単純な乗算ポテンシャル $V_{xc}(\mathbf{r})$ ではなく、各軌道に作用するポテンシャル演算子 $\hat{V}_{xc}^{mGGA}$ として定義され、その具体的な作用は軌道共役 $\phi_i^*(\mathbf{r})$ に対する変分から導かれる。 \(\hat{V}_{xc}^{mGGA} \phi_i(\mathbf{r}) = \frac{\delta E_{xc}^{mGGA}}{\delta \phi_i^*(\mathbf{r})}\)

3.2 運動エネルギー密度 $\tau$ の変分と部分積分による非局所演算子の導出

軌道共役の微小変化 $\delta \phi_i^(\mathbf{r})$ に対する交換相関エネルギーの変分 $\delta_{\phi_i^} E_{xc}^{mGGA}$ を考える。合成関数の微分公式（チェインルール）を適用すると、変分は以下のように記述できる。 \(\delta_{\phi_i^*} E_{xc}^{mGGA} = \int \left[ \frac{\partial f}{\partial \rho} \delta_{\phi_i^*}\rho + \frac{\partial f}{\partial \sigma} \delta_{\phi_i^*}\sigma + \frac{\partial f}{\partial \tau} \delta_{\phi_i^*}\tau \right] d\mathbf{r}\)

ここで、それぞれの変数に対する $\delta\phi_i^*$ 由来の変分は以下の通りである。

  1. 電子密度 $\rho$ の変分: \(\delta_{\phi_i^*}\rho = \phi_i \delta\phi_i^*\)
  2. 勾配二乗 $\sigma$ の変分: \(\delta_{\phi_i^*}\sigma = 2 \nabla\rho \cdot \nabla(\delta_{\phi_i^*}\rho) = 2 \nabla\rho \cdot \nabla(\phi_i \delta\phi_i^*)\)
  3. 運動エネルギー密度 $\tau$ の変分: \(\delta_{\phi_i^*}\tau = \frac{1}{2} \nabla(\delta\phi_i^*) \cdot \nabla\phi_i\)

これらをエネルギー変分式に代入する。 \(\delta_{\phi_i^*} E_{xc}^{mGGA} = \int \left[ \frac{\partial f}{\partial \rho} \phi_i \delta\phi_i^* + 2 \frac{\partial f}{\partial \sigma} \nabla\rho \cdot \nabla(\phi_i \delta\phi_i^*) + \frac{1}{2} \frac{\partial f}{\partial \tau} \nabla\phi_i \cdot \nabla(\delta\phi_i^*) \right] d\mathbf{r}\)

右辺第2項および第3項には変分項の勾配（$\nabla(\phi_i \delta\phi_i^)$ および $\nabla(\delta\phi_i^)$）が含まれているため、これらに対して空間的な部分積分を実行し、$\delta\phi_i^*$ が直接掛かる形へと変形する。

1. 第2項（勾配項）の部分積分

第2.1節と同様に、ベクトル解析の恒等式 $\nabla \cdot (\phi \mathbf{A}) = \nabla\phi \cdot \mathbf{A} + \phi(\nabla \cdot \mathbf{A})$ を用いる。ここで $\phi = \phi_i \delta\phi_i^*$、$\mathbf{A} = 2 \frac{\partial f}{\partial \sigma} \nabla\rho$ と置くと、 \(2 \frac{\partial f}{\partial \sigma} \nabla\rho \cdot \nabla(\phi_i \delta\phi_i^*) = \nabla \cdot \left( 2 \frac{\partial f}{\partial \sigma} \nabla\rho \, \phi_i \delta\phi_i^* \right) - \left[ \nabla \cdot \left( 2 \frac{\partial f}{\partial \sigma} \nabla\rho \right) \right] \phi_i \delta\phi_i^*\) 周期境界条件または無限遠での減衰条件の下では、右辺第1項の全発散項の積分はゼロとなるため、 \(\int 2 \frac{\partial f}{\partial \sigma} \nabla\rho \cdot \nabla(\phi_i \delta\phi_i^*) d\mathbf{r} = - \int \left[ \nabla \cdot \left( 2 \frac{\partial f}{\partial \sigma} \nabla\rho \right) \right] \phi_i \delta\phi_i^* d\mathbf{r}\)

2. 第3項（運動エネルギー密度項）の部分積分

同様に、恒等式 $\nabla \cdot (u \mathbf{B}) = \nabla u \cdot \mathbf{B} + u(\nabla \cdot \mathbf{B})$ において、スカラー場 $u = \delta\phi_i^*$、ベクトル場 $\mathbf{B} = \frac{1}{2} \frac{\partial f}{\partial \tau} \nabla\phi_i$ と置くと、 \(\frac{1}{2} \frac{\partial f}{\partial \tau} \nabla\phi_i \cdot \nabla(\delta\phi_i^*) = \nabla \cdot \left( \frac{1}{2} \frac{\partial f}{\partial \tau} \nabla\phi_i \, \delta\phi_i^* \right) - \left[ \nabla \cdot \left( \frac{1}{2} \frac{\partial f}{\partial \tau} \nabla\phi_i \right) \right] \delta\phi_i^*\) ここでも周期境界条件により全発散項（第1項）の積分は完全に消失するため、 \(\int \frac{1}{2} \frac{\partial f}{\partial \tau} \nabla\phi_i \cdot \nabla(\delta\phi_i^*) d\mathbf{r} = - \int \left[ \nabla \cdot \left( \frac{1}{2} \frac{\partial f}{\partial \tau} \nabla\phi_i \right) \right] \delta\phi_i^* d\mathbf{r}\)

3. 式の結合と演算子の同定

以上の部分積分結果を元の変分式に代入し、共通因子 $\delta\phi_i^*(\mathbf{r})$ で整理する。 \(\delta_{\phi_i^*} E_{xc}^{mGGA} = \int \left\{ \left[ \frac{\partial f}{\partial \rho} - \nabla \cdot \left( 2 \frac{\partial f}{\partial \sigma} \nabla\rho \right) \right] \phi_i - \frac{1}{2} \nabla \cdot \left( \frac{\partial f}{\partial \tau} \nabla\phi_i \right) \right\} \delta\phi_i^* d\mathbf{r}\)

軌道変分の定義式 $\delta_{\phi_i^} E_{xc}^{mGGA} = \int \left( \hat{V}_{xc}^{mGGA} \phi_i \right) \delta\phi_i^ d\mathbf{r}$ と比較することにより、meta-GGA交換相関ポテンシャル演算子 $\hat{V}_{xc}^{mGGA}$ の具体的な作用形式が得られる。 \(\hat{V}_{xc}^{mGGA} \phi_i = \left[ \frac{\partial f}{\partial \rho} - \nabla \cdot \left( 2 \frac{\partial f}{\partial \sigma} \nabla \rho \right) \right] \phi_i - \frac{1}{2} \nabla \cdot \left( \frac{\partial f}{\partial \tau} \nabla \phi_i \right)\)

ここで、局所部分（GGAポテンシャルと同様の形式）を $V_{xc,\text{local}}(\mathbf{r})$ と定義し、運動エネルギー密度微分係数を $\nu_3(\mathbf{r}) = \frac{\partial f}{\partial \tau}(\mathbf{r})$ と置くと、演算子は次のように記述される。 \(\hat{V}_{xc}^{mGGA} \phi_i(\mathbf{r}) = V_{xc,\text{local}}(\mathbf{r}) \phi_i(\mathbf{r}) - \frac{1}{2} \nabla \cdot \left( \nu_3(\mathbf{r}) \nabla \phi_i(\mathbf{r}) \right)\) この演算子の第2項 $-\frac{1}{2} \nabla \cdot (\nu_3 \nabla)$ は非局所演算子（または微分演算子）であり、波動関数に直接作用するため、通常の乗算型ポテンシャルとは本質的に異なる実装アプローチが必要となる。

3.3 平面波基底における非局所ハミルトニアン作用の評価とFFT実装

平面波基底ポテンシャルコードにおいて、波動関数 $\phi_i(\mathbf{r})$ は逆空間のフーリエ係数 $c_i(\mathbf{G})$ によって表現される。 \(\phi_i(\mathbf{r}) = \frac{1}{\sqrt{\Omega}} \sum_{\mathbf{G}} c_i(\mathbf{G}) e^{i \mathbf{G} \cdot \mathbf{r}}\)

ハミルトニアンの作用 $\hat{H}\phi_i$ を評価する際、非局所ポテンシャル項 $\hat{V}_{xc,\text{NL}}^{mGGA} \phi_i = -\frac{1}{2} \nabla \cdot \left( \nu_3(\mathbf{r}) \nabla \phi_i(\mathbf{r}) \right)$ の逆空間における表現を得る必要がある。逆空間において、空間微分演算子 $\nabla$ および $\nabla\cdot$ は、それぞれ $i\mathbf{G}$ との積および内積に対応する。

この双対性を利用し、FFTを用いて実空間での積と逆空間での射影を組み合わせることで、効率的に評価を行うことができる。具体的な手順は以下の通りである。

  1. 逆空間での勾配評価: 波動関数の勾配を逆空間で計算する。 \(\tilde{\mathbf{g}}_i(\mathbf{G}) = i\mathbf{G} c_i(\mathbf{G})\) ここで $\tilde{\mathbf{g}}_i(\mathbf{G})$ は3次元ベクトル値のフーリエ係数である。
  2. 実空間への変換 (逆FFT): 勾配の3つのデカルト成分（$x, y, z$）を、それぞれ独立に逆高速フーリエ変換（$\text{FFT}^{-1}$）により実空間へ写す。 \(\mathbf{g}_i(\mathbf{r}_p) = \nabla \phi_i(\mathbf{r}_p) = \text{FFT}^{-1}\left[\tilde{\mathbf{g}}_i(\mathbf{G})\right]\)
  3. 実空間での物理量の積: 各実空間グリッド点 $\mathbf{r}_p$ において、汎関数微分から得られる実数スカラーポテンシャル $\nu_3(\mathbf{r}_p) = \frac{\partial f}{\partial \tau}(\mathbf{r}_p)$ を、勾配ベクトルに乗算する。 \(\mathbf{w}_i(\mathbf{r}_p) = \nu_3(\mathbf{r}_p) \mathbf{g}_i(\mathbf{r}_p) = \nu_3(\mathbf{r}_p) \nabla \phi_i(\mathbf{r}_p)\)
  4. 逆空間への変換 (正FFT): 得られた実空間ベクトル場 $\mathbf{w}_i(\mathbf{r}_p)$ の3つの成分を、それぞれ順方向フーリエ変換（$\text{FFT}$）により逆空間に戻す。 \(\tilde{\mathbf{w}}_i(\mathbf{G}) = \text{FFT}\left[\mathbf{w}_i(\mathbf{r}_p)\right]\)
  5. 逆空間での発散評価と係数更新: 逆空間において発散演算（$i\mathbf{G}$ との内積）を適用し、非局所ポテンシャル成分のフーリエ係数を構成する。 \(\tilde{h}_{i,\text{NL}}(\mathbf{G}) = -\frac{1}{2} \left( i\mathbf{G} \cdot \tilde{\mathbf{w}}_i(\mathbf{G}) \right)\)

最終的に、この $\tilde{h}{i,\text{NL}}(\mathbf{G})$ は、局所ポテンシャル項の作用 $V{xc,\text{local}}\phi_i$ からの寄与 $\tilde{h}{i,\text{local}}(\mathbf{G}) = \text{FFT}[V{xc,\text{local}}(\mathbf{r})\phi_i(\mathbf{r})]$ と足し合わされ、自己無撞着なハミルトニアン作用の評価に用いられる。

このFFTを用いた演算スキームは、各軌道に対して合計6回の3D FFT（逆変換3回、正変換3回）を要求するため、通常の局所ポテンシャル（正・逆各1回）と比較して計算コストが高い。そのため、実装においてはスピンごとの対称性や並列化による並列化効率の最大化が図られる。

3.4 行列要素 $\langle \phi_j | \hat{V}_{xc}^{mGGA} | \phi_i \rangle$ の積分表現と実空間グリッド上の数値的評価

軌道 $\phi_j$ と $\phi_i$ の間における meta-GGA 交換相関ポテンシャル演算子の行列要素 $\langle \phi_j | \hat{V}_{xc}^{mGGA} | \phi_i \rangle$ は、次のように定義される。 \(\langle \phi_j | \hat{V}_{xc}^{mGGA} | \phi_i \rangle = \int \phi_j^*(\mathbf{r}) V_{xc,\text{local}}(\mathbf{r}) \phi_i(\mathbf{r}) d\mathbf{r} - \frac{1}{2} \int \phi_j^*(\mathbf{r}) \left[ \nabla \cdot \left( \nu_3(\mathbf{r}) \nabla \phi_i(\mathbf{r}) \right) \right] d\mathbf{r}\)

第2項の非局所ハミルトニアン寄与に対して部分積分を実行する。ここでも表面項が完全に相殺される条件を適用すると、以下の対称的な積分表現が得られる。 \(- \frac{1}{2} \int \phi_j^*(\mathbf{r}) \left[ \nabla \cdot \left( \nu_3(\mathbf{r}) \nabla \phi_i(\mathbf{r}) \right) \right] d\mathbf{r} = \frac{1}{2} \int \nu_3(\mathbf{r}) \nabla \phi_j^*(\mathbf{r}) \cdot \nabla \phi_i(\mathbf{r}) d\mathbf{r}\)

したがって、行列要素全体は、以下の二つの積分項の和として表現される。 \(\langle \phi_j | \hat{V}_{xc}^{mGGA} | \phi_i \rangle = \int V_{xc,\text{local}}(\mathbf{r}) \phi_j^*(\mathbf{r}) \phi_i(\mathbf{r}) d\mathbf{r} + \frac{1}{2} \int \nu_3(\mathbf{r}) \nabla \phi_j^*(\mathbf{r}) \cdot \nabla \phi_i(\mathbf{r}) d\mathbf{r}\)

この積分表現を実空間グリッド点 $\mathbf{r}_p$ （積分の重み $w_p = \Omega/N_g$）の上で離散化すると、以下の数値積分和へと帰着される。 \(\langle \phi_j | \hat{V}_{xc}^{mGGA} | \phi_i \rangle \approx \sum_{p=1}^{N_g} w_p \left[ V_{xc,\text{local}}(\mathbf{r}_p) \phi_j^*(\mathbf{r}_p) \phi_i(\mathbf{r}_p) + \frac{1}{2} \nu_3(\mathbf{r}_p) \nabla \phi_j^*(\mathbf{r}_p) \cdot \nabla \phi_i(\mathbf{r}_p) \right]\)

数値的実装の特徴と利点

  1. エルミート性（Hermiticity）の厳密な保証: この対称的な積分表現を用いることで、離散化された実空間グリッド上においてもポテンシャル行列のエルミート対称性（$\langle \phi_j | \hat{V}{xc}^{mGGA} | \phi_i \rangle = \langle \phi_i | \hat{V}{xc}^{mGGA} | \phi_j \rangle^$）が数値誤差の範囲内で厳密に保証される。もし発散項 $\nabla \cdot (\nu_3 \nabla \phi_i)$ を直接数値微分で計算して軌道 $\phi_j^$ との積を取る非対称な定式化を行った場合、差分法やFFTグリッドの有限カットオフに起因する非対称な数値誤差が蓄積し、ハミルトニアンの対称性が崩れ、SCFの数値的不安定化や実数固有値の乖離を招く恐れがある。
  2. 2階空間微分の回避による数値的安定性: 部分積分を行うことで、$\nu_3(\mathbf{r})$ の空間微分や $\nu_3 \nabla\phi_i$ の発散といった、ノイズを増幅しやすい2階の空間微分を完全に回避することができる。必要とされるのは波動関数の一階勾配 $\nabla \phi_i$ および $\nabla \phi_j$ のみであり、これらは波数空間で極めて高精度かつ安定に算出できるため、全系でのエネルギーおよび力の計算精度が大幅に向上する。

4. libxc API の設計思想と変数マッピング

実空間グリッド上で交換相関エネルギー密度および各種ポテンシャル（または汎関数微分係数）を評価する際、スクラッチからすべての汎関数の数学的定義やパラメータ、その微分係数をハードコードすることは、実装およびメンテナンスの観点から極めて非効率である。この課題を解決するため、密度汎関数理論（DFT）計算コミュニティにおいてデファクトスタンダードとして広く利用されているのが、交換相関汎関数ライブラリである libxc である。

本節では、libxc の API 設計思想、実空間グリッドデータとの変数マッピング、および実際の C 言語と Fortran による典型的な評価ループの実装例について詳述する。

4.1 物理単位とグリッドデータの配列構造

libxc は、ポータビリティと計算の簡素化のために、一貫して以下の物理単位（原子単位系）を前提として設計されている。

  • 長さの単位: ボーア ($a_0$)
  • エネルギーの単位: ハートリー ($\text{Hartree}$)

したがって、入力される電子密度 $\rho$ は $[{\text{Bohr}}^{-3}]$、その勾配 $\sigma$ は $[{\text{Bohr}}^{-8}]$、運動エネルギー密度 $\tau$ は $[\text{Hartree} \cdot {\text{Bohr}}^{-3}]$ （または $[{\text{Bohr}}^{-5}]$）の単位で正規化されている必要がある。

また、libxc のすべての評価ルーチンは、実空間の空間グリッド情報を 1 次元にフラット化した配列として受け取る。3 次元 FFT メッシュ上の座標点 $(i, j, k)$ は、あらかじめ適切なインデックスマッピングによって 1 次元配列へと変換され、入力データ（rho、sigma、tau など）として受け渡される。

1. スピン極性のない場合（制限付き計算: Restricted / Spin-Unpolarized）

スピン極性を考慮しない計算において、グリッド点数を $N_g$ とすると、各変数は以下のように 1 次元配列に対応する。

  • rho: サイズ $N_g$ の配列。各要素は全電子密度 $\rho(\mathbf{r}_i)$。
  • sigma: サイズ $N_g$ の配列。各要素は密度勾配の自乗 $\sigma(\mathbf{r}_i) = \nabla \rho(\mathbf{r}_i) ^2$。
  • tau: サイズ $N_g$ の配列。各要素は運動エネルギー密度 $\tau(\mathbf{r}_i)$。

2. コロニアルスピン極性がある場合（制限なし計算: Unrestricted / Spin-Polarized）

スピン分極を考慮する場合、上向きスピン（$\uparrow$）と下向きスピン（$\downarrow$）に対応する物理量を隣接するインデックスに格納する。

  • rho: サイズ $2 N_g$ の配列。奇数インデックス（C言語の 2*i）に $\rho_\uparrow(\mathbf{r}i)$、偶数インデックス（2*i+1）に $\rho\downarrow(\mathbf{r}_i)$ を格納する。
  • sigma: サイズ $3 N_g$ の配列。各点 $\mathbf{r}_i$ に対して以下の 3 つの勾配量を順に格納する。
    1. $\sigma_{\uparrow\uparrow}(\mathbf{r}_i) = \nabla \rho_\uparrow(\mathbf{r}_i) ^2$ （インデックス 3*i）
    2. $\sigma_{\uparrow\downarrow}(\mathbf{r}i) = \nabla \rho\uparrow(\mathbf{r}i) \cdot \nabla \rho\downarrow(\mathbf{r}_i)$ （インデックス 3*i+1）
    3. $\sigma_{\downarrow\downarrow}(\mathbf{r}_i) = \nabla \rho_\downarrow(\mathbf{r}_i) ^2$ （インデックス 3*i+2）
  • tau: サイズ $2 N_g$ の配列。奇数インデックス（2*i）に $\tau_\uparrow(\mathbf{r}i)$、偶数インデックス（2*i+1）に $\tau\downarrow(\mathbf{r}_i)$ を格納する。

4.2 libxc が返却する汎関数微分係数（$\nu_1, \nu_2, \nu_3$）と Kohn-Sham ポテンシャルへのマッピング

libxc の主要な役割は、各グリッド点における入力値に対して、交換相関エネルギー密度 $\epsilon_{xc}$ （粒子あたりのエネルギー）と、その各入力変数に対する偏微分係数（汎関数微分係数）を返すことである。

libxc では、全交換相関エネルギーが $E_{xc} = \int \rho \epsilon_{xc} d\mathbf{r}$ と定義されているため、エネルギー密度関数 $f$ は $f(\rho, \sigma, \tau) = \rho \epsilon_{xc}(\rho, \sigma, \tau)$ に相当する。したがって、libxc から出力される偏微分値（vrho、vsigma、vtau）は、本稿のセクション 2 および 3 で導入した汎関数微分係数 $\nu_1, \nu_2, \nu_3$ と以下のように厳密に 1 対 1 で対応する。

  1. $\nu_1$ (変数 vrho): 電荷密度 $\rho$ に関する微分 \(\nu_1(\mathbf{r}) = \frac{\partial (\rho \epsilon_{xc})}{\partial \rho} = \epsilon_{xc} + \rho \frac{\partial \epsilon_{xc}}{\partial \rho}\)
    • LDA, GGA, meta-GGA すべてのレベルで評価される。
    • スピン極性ありの場合、vrho のサイズは $2 N_g$ となり、上向きスピンに対する $\nu_{1,\uparrow} = \frac{\partial f}{\partial \rho_\uparrow}$ と下向きスピンに対する $\nu_{1,\downarrow} = \frac{\partial f}{\partial \rho_\downarrow}$ が交互に格納される。
  2. $\nu_2$ (変数 vsigma): 密度勾配の自乗 $\sigma$ に関する微分 \(\nu_2(\mathbf{r}) = \frac{\partial (\rho \epsilon_{xc})}{\partial \sigma} = \rho \frac{\partial \epsilon_{xc}}{\partial \sigma}\)
    • GGA および meta-GGA のレベルで評価される。
    • スピン極性ありの場合、vsigma のサイズは $3 N_g$ となり、$\nu_{2,\uparrow\uparrow} = \frac{\partial f}{\partial \sigma_{\uparrow\uparrow}}$、$\nu_{2,\uparrow\downarrow} = \frac{\partial f}{\partial \sigma_{\uparrow\downarrow}}$、$\nu_{2,\downarrow\downarrow} = \frac{\partial f}{\partial \sigma_{\downarrow\downarrow}}$ がこの順に格納される。
  3. $\nu_3$ (変数 vtau): 運動エネルギー密度 $\tau$ に関する微分 \(\nu_3(\mathbf{r}) = \frac{\partial (\rho \epsilon_{xc})}{\partial \tau} = \rho \frac{\partial \epsilon_{xc}}{\partial \tau}\)
    • meta-GGA レベルのみで評価される。
    • スピン極性ありの場合、vtau のサイズは $2 N_g$ となり、$\nu_{3,\uparrow} = \frac{\partial f}{\partial \tau_\uparrow}$、$\nu_{3,\downarrow} = \frac{\partial f}{\partial \tau_\downarrow}$ が交互に格納される。

これらの出力変数を用いることで、実空間上の Kohn-Sham ポテンシャル演算子は以下のように完全に組み立てられる。

  • GGA ポテンシャル: \(V_{xc}(\mathbf{r}) = \nu_1(\mathbf{r}) - \nabla \cdot \left( 2 \nu_2(\mathbf{r}) \nabla \rho(\mathbf{r}) \right)\)
  • meta-GGA ハミルトニアン作用素: \(\hat{V}_{xc}^{mGGA} \phi_i(\mathbf{r}) = \left[ \nu_1(\mathbf{r}) - \nabla \cdot \left( 2 \nu_2(\mathbf{r}) \nabla \rho(\mathbf{r}) \right) \right] \phi_i(\mathbf{r}) - \frac{1}{2} \nabla \cdot \left( \nu_3(\mathbf{r}) \nabla \phi_i(\mathbf{r}) \right)\)

4.3 C 言語および Fortran 2003 による評価実装コードスニペット

以下に、実空間グリッド $N_g = 5$ のダミーデータに対して meta-GGA 汎関数（例として SCAN 交換汎関数）を初期化、評価、クリーンアップする典型的なコードを示す。

1. C 言語による実装例 (test_libxc.c)

#include <stdio.h>
#include <stdlib.h>
#include <xc.h>

int main(void) {
    int ngrid = 5;
    
    // 入力データの配列確保 (Bohr, Hartree 単位系を想定)
    double *rho   = (double *)malloc(ngrid * sizeof(double));
    double *sigma = (double *)malloc(ngrid * sizeof(double));
    double *tau   = (double *)malloc(ngrid * sizeof(double));
    
    // 出力データの配列確保
    double *exc    = (double *)malloc(ngrid * sizeof(double));
    double *vrho   = (double *)malloc(ngrid * sizeof(double));
    double *vsigma = (double *)malloc(ngrid * sizeof(double));
    double *vtau   = (double *)malloc(ngrid * sizeof(double));

    // ダミーの入力データの入力 (実際にはDFTコード内部のFFTグリッド上の値が入る)
    for (int i = 0; i < ngrid; i++) {
        rho[i]   = 0.1;   // 電子密度
        sigma[i] = 0.01;  // 密度勾配の自乗
        tau[i]   = 0.05;  // 運動エネルギー密度
    }

    // libxc 汎関数構造体の宣言と初期化
    xc_func_type func;
    // SCAN 交換汎関数 (XC_MGGA_X_SCAN), 制限なし(XC_UNPOLARIZED)で初期化
    int info = xc_func_init(&func, XC_MGGA_X_SCAN, XC_UNPOLARIZED);
    if (info < 0) {
        fprintf(stderr, "Error: SCAN functional initialization failed.\n");
        return EXIT_FAILURE;
    }

    // meta-GGA の評価を実行
    // 引数: 構造体ポインタ, グリッド点数, rho, sigma, laplacian (NULL指定), tau, 
    //      エネルギー密度出力, vrho出力, vsigma出力, vlapl出力 (NULL指定), vtau出力
    xc_mgga_eval(&func, ngrid, rho, sigma, NULL, tau, exc, vrho, vsigma, NULL, vtau);

    // 評価結果の出力
    printf("--- libxc Meta-GGA (SCAN Exchange) Evaluation Results ---\n");
    for (int i = 0; i < ngrid; i++) {
        printf("Grid %d: zk = %9.6f, vrho = %9.6f, vsigma = %9.6f, vtau = %9.6f\n", 
               i, exc[i], vrho[i], vsigma[i], vtau[i]);
    }

    // 汎関数の後処理とメモリの解放
    xc_func_end(&func);
    free(rho); free(sigma); free(tau);
    free(exc); free(vrho); free(vsigma); free(vtau);

    return EXIT_SUCCESS;
}

2. Fortran 2003 による実装例 (test_libxc.f90)

現代の DFT コードでは、C言語で書かれた libxc コアライブラリと Fortran 2003 の iso_c_binding を介したインターフェースを利用することが一般的である。

program test_libxc
  use xc_f03_lib_m
  implicit none

  integer, parameter :: ngrid = 5
  type(xc_f03_func_t) :: func
  integer :: info, i
  
  ! 配列の宣言
  real(8), dimension(:), allocatable :: rho, sigma, tau
  real(8), dimension(:), allocatable :: exc, vrho, vsigma, vtau
  real(8), dimension(:), allocatable :: lapl, vlapl

  ! 配列の割り当て
  allocate(rho(ngrid), sigma(ngrid), tau(ngrid))
  allocate(exc(ngrid), vrho(ngrid), vsigma(ngrid), vtau(ngrid))
  allocate(lapl(ngrid), vlapl(ngrid))

  ! ダミー入力データの設定
  rho(:)   = 0.1d0
  sigma(:) = 0.01d0
  tau(:)   = 0.05d0
  lapl(:)  = 0.0d0

  ! 汎関数の初期化 (SCAN交換汎関数, スピン分極なし)
  info = xc_f03_func_init(func, XC_MGGA_X_SCAN, XC_UNPOLARIZED)
  if (info < 0) then
     write(*,*) "Error: SCAN functional initialization failed."
     stop
  end if

  ! 汎関数の評価
  ! Fortran 2003 バインディングでは、配列の先頭要素を渡すか、インターフェースに従う
  ! laplacian と vlapl にはダミーの配列要素を実引数として渡す
  call xc_f03_mgga_eval(func, ngrid, rho(1), sigma(1), lapl(1), tau(1), &
       exc(1), vrho(1), vsigma(1), vlapl(1), vtau(1))

  ! 結果の出力
  write(*,*) "--- libxc Meta-GGA (SCAN Exchange) Fortran Results ---"
  do i = 1, ngrid
     write(*, '(A,I0,A,F9.6,A,F9.6,A,F9.6,A,F9.6)') &
          "Grid ", i, ": zk = ", exc(i), ", vrho = ", vrho(i), &
          ", vsigma = ", vsigma(i), ", vtau = ", vtau(i)
  end do

  ! 後処理
  call xc_f03_func_end(func)
  deallocate(rho, sigma, tau, exc, vrho, vsigma, vtau, lapl, vlapl)

end program test_libxc

5. Quantum ESPRESSO における交換相関処理の実装とカスタム汎関数の追加手順

オープンソースの第一原理計算パッケージである Quantum ESPRESSO (QE) は、平面波基底ポテンシャル法を採用した DFT コードの代表例である。QE における交換相関汎関数の処理構造を理解することは、物性物理・材料科学の研究において独自の汎関数を実装したり、libxc との連携をカスタマイズしたりする上で非常に有益である。

本節では、QE の交換相関処理に関係する主要ファイルのコードマッピングと、ユーザーが独自の GGA / meta-GGA 汎関数を QE に直接組み込むためのステップバイステップの開発ガイドを提供する。

5.1 主要ファイルとモジュールの役割

Quantum ESPRESSO の PW/ ディレクトリおよび Modules/ ディレクトリには、交換相関エネルギーとポテンシャル計算を担う多くのルーチンが配置されている。その中心的な呼び出し経路と役割は以下の通りである。

graph TD
    v_of_rho[PW/src/v_of_rho.f90] -->|ポテンシャル評価要求| xc_pot[PW/src/xc_pot.f90]
    xc_pot -->|DFTフラグ判定| funct[Modules/funct.f90]
    xc_pot -->|局所ポテンシャル計算| xc[PW/src/xc.f90]
    xc_pot -->|libxc連携| libxc_inter[Modules/funct.f90 / libxc]
    xc -->|GGA/meta-GGA微分処理| funct

1. Modules/funct.f90 (DFT レジストリ)

QE 全体で使用する交換相関汎関数の種類（LDA、GGA、meta-GGA、Hybrid、libxc、DFT-D 等）を決定し、対応するグローバル論理フラグを管理する中心モジュールである。

  • ユーザーが入力ファイル（&SYSTEM カードの input_dft パラメータ）で指定した文字列（例：'pbe', 'scan', 'XC-000101' 等）は、ここで解析され、対応する内部インデックスや dft_is_gga、dft_is_meta といったフラグに変換される。
  • libxc を用いる場合、libxc の ID に対応する変換テーブルが定義されている。

2. PW/src/v_of_rho.f90 (ハミルトニアン構築コア)

SCF ループの各イテレーションにおいて、総電荷密度 $\rho(\mathbf{r})$ から Kohn-Sham ハミルトニアンの有効ポテンシャル $V_{\text{eff}}(\mathbf{r})$ を実空間グリッド上で構築するメインルーチンである。

  • 電荷密度の入力を受け取り、Hartree ポテンシャルや外部ポテンシャルを加算したのち、交換相関ポテンシャル計算のラッパーサブルーチンである xc_pot を呼び出す。

3. PW/src/xc_pot.f90 (実空間ポテンシャルアセンブラ)

実空間グリッドにおける電子密度およびスピン密度情報から、交換相関エネルギー $E_{xc}$ とポテンシャル $V_{xc}$ （非コロニアルの場合はスカラーポテンシャル $V_{xc}$ および有効磁場 $\mathbf{B}_{xc}$）を実際に計算し、組み立てるルーチンである。

  • 各点における電子密度 $\rho(\mathbf{r})$ およびその勾配 $\nabla \rho(\mathbf{r})$ の計算（FFT を用いた波数空間微分、または実空間差分）を実行する。
  • funct.f90 のフラグに基づいて、内蔵 Fortran 汎関数モジュール（xc.f90 等）または libxc の C API インターフェースを呼び出し、エネルギーおよび汎関数微分値を得る。
  • 得られた微分係数から空間発散 $\nabla \cdot (2 \nu_2 \nabla \rho)$ などの 2 階微分寄与を計算し、全ポテンシャル配列に加算する。

4. PW/src/xc.f90 (内蔵交換相関カーネル)

Quantum ESPRESSO がネイティブで実装している LDA および GGA 汎関数（PBE, BLYP, PW91 など）の具体的な数式評価サブルーチン群が定義されているファイルである。

  • exc および vxc を算出する純粋な Fortran サブルーチンが含まれており、汎関数の計算ロジック自体が記述されている。

5.2 QE に独自の GGA / meta-GGA 汎関数を追加するための開発ステップ

研究者が新規の GGA もしくは meta-GGA 汎関数を開発し、それを Quantum ESPRESSO に直接組み込んで SCF 計算を実行したい場合、以下の具体的なステップを踏むことで実装可能である。

ステップ 1: Modules/funct.f90 での登録とフラグ設定

まず、新しい汎関数を登録し、QE がそれを正しく認識できるようにする。

  1. funct.f90 を開き、汎関数を識別する新しいインデックス定数を定義する。 
    ! 例: カスタムGGA交換汎関数のインデックスを 99 と定義
INTEGER, PARAMETER :: MY_GGA_X = 99
  2. set_dft_from_name サブプロシージャ内に、入力文字列（例：'mygga'）とのマッチングを追加し、対応するフラグを有効にする。 
    IF ( matches( "mygga", dft_name ) ) THEN
   ! 内部の交換・相関汎関数インデックスを設定
   iexch = MY_GGA_X
   icorr = 1  ! 相関には既存のPZ-LDAなどを指定可能
   dft_is_gga = .TRUE.
   dft_is_meta = .FALSE.
END IF

ステップ 2: 汎関数カーネルの新規実装 (PW/src/xc.f90 または Modules/funct.f90)

数式に基づいて、$\rho, \sigma, \tau$ からエネルギーと微分係数を計算する実際のロジックを実装する。

  1. PW/src/xc.f90 に、新しく定義したサブルーチンを追加する。以下はスピン極性のない GGA 交換汎関数のテンプレートである。 
    SUBROUTINE my_gga_exchange( ngrid, rho, sigma, ex, vx, vsig )
  USE kinds, ONLY : DP
  IMPLICIT NONE
  INTEGER, INTENT(in) :: ngrid
  REAL(DP), INTENT(in) :: rho(ngrid), sigma(ngrid)
  REAL(DP), INTENT(out):: ex(ngrid), vx(ngrid), vsig(ngrid)
  INTEGER :: i
  REAL(DP) :: r, s, e, v1, v2
      
  DO i = 1, ngrid
     r = rho(i)
     s = sigma(i)
     IF ( r > 1.d-10 ) THEN
        ! ここにカスタムのエネルギー密度（粒子あたり）と微分の数式を記述する
        ! 例: パラメータ化された勾配補正項の導入
        e  = -0.73855896_DP * (r**(1.0_DP/3.0_DP)) - 0.005_DP * s / (r**(4.0_DP/3.0_DP))
        v1 = -0.98474528_DP * (r**(1.0_DP/3.0_DP)) + 0.00166667_DP * s / (r**(4.0_DP/3.0_DP))
        v2 = -0.005_DP / (r**(1.0_DP/3.0_DP))
            
        ex(i)    = e
        vx(i)    = v1  ! vrho に相当
        vsig(i)  = v2  ! vsigma に相当
     ELSE
        ex(i)   = 0.0_DP
        vx(i)   = 0.0_DP
        vsig(i) = 0.0_DP
     END IF
  END DO
END SUBROUTINE my_gga_exchange

ステップ 3: 評価ディスパッチャーへのフック

作成したサブルーチンを、既存のポテンシャル生成ループから呼び出されるディスパッチャーに登録する。

  1. PW/src/xc.f90 もしくは Modules/funct.f90 の中にある、GGAのエネルギーとポテンシャルを計算するディスパッチャーサブルーチン（例：egga または vgga など）の SELECT CASE (iexch) ブロックを探す。
  2. 以下のように、カスタムインデックスに対する case 条件を追加する。 
    SELECT CASE ( iexch )
   ! ... 他の既存汎関数 ...
   CASE ( MY_GGA_X )
      CALL my_gga_exchange( ngrid, rho, sigma, ex, vx, vsig )
END SELECT

ステップ 4: 自動微分とポテンシャル組み立ての確認

ユーザーが独自の汎関数サブルーチンで vrho ($\nu_1$) と vsigma ($\nu_2$) を正しく返却していれば、QE 側の PW/src/xc_pot.f90 がその値を自動的に受け取り、実空間グリッド上で自動的に空間微分と発散項の評価（$\nabla \cdot (2\nu_2 \nabla\rho)$）を実行する。 したがって、空間微分やポテンシャルの発散、運動エネルギー非局所演算子の作用といった複雑な物理処理を、ユーザー自身で新規に実装し直す必要は一切ない。これが、QE の優れた設計思想であり、開発者は純粋に「局所的な数式とその偏微分」の記述だけに集中することができる。

ステップ 5: コンパイルと動作テスト

  1. Quantum ESPRESSO のルートディレクトリで再コンパイルを実行する。 
    make pw
  2. 簡単な入力ファイルを作成し、input_dft = 'mygga' を指定して SCF 計算を実行し、新しい汎関数による全エネルギーと力が正常に計算され、SCF が収束することを確認する。

6. 結論

本稿では、密度汎関数理論（DFT）における交換相関汎関数を、理論的な数式モデルから実際の計算機プログラムへ落とし込むための具体的な実装アプローチについて網羅的に解説した。

実空間数値積分におけるスピン極性の局所対角化および非コロニアル極性磁場への逆変換から始まり、一般化勾配近似（GGA）で必要とされる2階空間微分の離散化スキーム（有限差分 vs FFT）、さらには meta-GGA 汎関数で現れる軌道依存の非局所運動エネルギーポテンシャル演算子の平面波・実空間ハイブリッド処理に至るまで、DFT 実装の物理的・数学的基盤を明らかにした。

また、現代の DFT 開発において不可欠な libxc API との変数・微分係数マッピング、さらには Quantum ESPRESSO のコードアーキテクチャに基づいたカスタム汎関数の統合ガイドを通じ、理論モデルがどのようにコードへと具現化され、研究の現場でカスタマイズされるかを実践的に示した。

電子構造計算コードは一見すると極めて複雑なブラックボックスに見えるが、その深部にある交換相関処理は、変分原理と部分積分という解析力学の基本的な数学ツール、そして FFT や有限差分といった離散数値解析技術の見事な融合によって成り立っている。本解説が、第一原理計算プログラムの内部構造を理解し、より高度な機能拡張や新規手法開発を目指す研究者や開発者にとっての確かなロードマップとなれば幸いである。

You May Also Enjoy