【計算化学】SCC-DFTB (Self-Consistent-Charge Density Functional Tight-Binding) の数理的定式化と実装論：電荷移動系への拡張と2次摂動の導出

最終更新：2025-12-31

注意: この記事はAIによって自動生成されたものです。正確な数式やパラメータ定義については、必ず引用元の原著論文（M. Elstner et al., Phys. Rev. B 58, 7260 (1998)）をご確認ください。

序論：DFTB1の限界と電荷自己無撞着性の必要性#

密度汎関数タイトバインディング法（DFTB）の初期形式（PorezagらによるNon-SCC-DFTB、通称DFTB1）は、ホモ核共有結合系において高い計算効率と精度を両立させた。しかし、この手法はKohn-Shamエネルギーの電子密度揺らぎに関する展開を1次で打ち切る近似に基づいているため、系内の長距離静電相互作用や電荷移動（Charge Transfer）を記述できないという原理的な限界を有していた。

例えば、極性分子、ヘテロ核結合、あるいは生体分子のような電荷の偏りが機能に直結する系において、DFTB1は双極子モーメントを過小評価し、結合エネルギーの精度が著しく低下する傾向にあった。

これに対処するため、1998年にM. Elstner、D. Porezag、G. Seifertらは、密度揺らぎの2次項までを考慮に入れた拡張手法、SCC-DFTB (Self-Consistent-Charge DFTB) を提唱した。この手法は、原子上の部分電荷（Partial Charges）と静電ポテンシャルを相互に依存させながら反復的に決定する「自己無撞着（Self-Consistent）」な手続きを導入することで、化学反応や分極現象の記述能力を飛躍的に向上させた。本稿では、SCC-DFTBの数理的導出から実装アルゴリズムまでを詳細に解説する。

1. 数理的背景：エネルギー汎関数の2次摂動展開#

SCC-DFTBの理論的骨子は、DFTの総エネルギーを参照密度 $\rho_0$ の周りで2次までTaylor展開することにある。

1.1 エネルギー展開式#

全電子密度 $\rho(\mathbf{r})$ を、参照密度 $\rho_0(\mathbf{r})$ （中性原子密度の和）と、そこからの揺らぎ $\delta\rho(\mathbf{r})$ に分解する。

\rho(\mathbf{r}) = \rho_0(\mathbf{r}) + \delta\rho(\mathbf{r}), \quad \rho_0(\mathbf{r}) = \sum_A \rho_0^A(\mathbf{r})

Kohn-Shamエネルギー $E^{KS}[\rho]$ を展開すると、以下のようになる。

E^{KS}[\rho_0 + \delta\rho] = E^{(0)}[\rho_0] + E^{(1)}[\rho_0, \delta\rho] + E^{(2)}[\rho_0, (\delta\rho)^2] + \dots

0次項と1次項の和は、従来のDFTB1における「バンド構造エネルギー」と「斥力ポテンシャル」に対応する。

E_{band} + E_{rep} = \sum_{i}^{occ} \langle \psi_i | \hat{H}^0 | \psi_i \rangle + \frac{1}{2} \sum_{A \neq B} U_{rep}^{AB}(R_{AB})

ここで $\hat{H}^0$ は参照密度 $\rho_0$ に基づくハミルトニアンである。

1.2 2次項の導出と物理的意味#

SCC-DFTBの核心は、以下の2次項 $E^{(2)}$ の取り扱いにある。

E^{(2)} = \frac{1}{2} \iint \delta\rho(\mathbf{r}) \left( \frac{1}{|\mathbf{r} - \mathbf{r}'|} + \frac{\delta^2 E_{xc}}{\delta\rho(\mathbf{r})\delta\rho(\mathbf{r}')} \right) \delta\rho(\mathbf{r}') d\mathbf{r} d\mathbf{r}'

括弧内の第一項は古典的なCoulomb相互作用（Hartree項）、第二項は交換相関エネルギーの2次応答（XCカーネル）を表す。

1.3 密度揺らぎの単極子近似 (Monopole Approximation)#

積分を効率的に評価するために、密度揺らぎ $\delta\rho$ を各原子 $A$ に局在した球対称な電荷分布 $F_{00}^A$ の重ね合わせとして近似する。これは多重極展開における単極子（Monopole）項のみを残すことに相当する。

\delta\rho(\mathbf{r}) \approx \sum_A \Delta q_A F_{00}^A(|\mathbf{r} - \mathbf{R}_A|)

ここで $\Delta q_A$ は原子 $A$ 上の正味の電荷変動（Mulliken電荷などから算出される中性原子からの偏差）である。これを $E^{(2)}$ の式に代入すると、SCCエネルギー項 $E_{SCC}$ が得られる。

E_{SCC} = \frac{1}{2} \sum_{A} \sum_{B} \Delta q_A \Delta q_B \gamma_{AB}

1.4 $\gamma$ 行列とHubbardパラメータ#

ここで導入された $\gamma_{AB}$ は、原子Aと原子Bの間の有効クーロン相互作用を表す。

\gamma_{AB} = \iint F_{00}^A(\mathbf{r}) \left( \frac{1}{|\mathbf{r} - \mathbf{r}'|} + \frac{\delta^2 E_{xc}}{\delta\rho\delta\rho'} \right) F_{00}^B(\mathbf{r}') d\mathbf{r} d\mathbf{r}'

この $\gamma_{AB}$ の具体的な関数形は、以下の極限挙動を満たすように解析的にモデル化される。

長距離極限 ( $R_{AB} \to \infty$ ): 2つの点電荷間のCoulomb相互作用に漸近する。 $\gamma_{AB} \to \frac{1}{R_{AB}}$
オンサイト極限 ( $R_{AB} \to 0, A=B$ ): 同一原子上の電子間反発を表すHubbardパラメータ $U_A$ に収束する。 $\gamma_{AA} \approx U_A = \text{IE}_A - \text{EA}_A$ ここで IE はイオン化ポテンシャル、EA は電子親和力である（厳密には化学的硬さ $\eta$ に関連する）。

Elstnerらは、電荷分布 $F_{00}^A$ にSlater型関数やGaussian型関数を仮定することで、以下の補間公式を導出した。

\gamma_{AB}(R_{AB}) = \frac{1}{\sqrt{R_{AB}^2 + \frac{1}{4} \left( \frac{1}{U_A} + \frac{1}{U_B} \right)^2 \tau_{AB}}} \quad (\text{一例})

※ 実際の論文では、指数関数を含むより複雑な形式が採用されている場合があるが、本質的には短距離の $U$ と長距離の $1/R$ を滑らかに繋ぐ関数である。

2. 実装論：SCCループと修正ハミルトニアン#

エネルギー最小化条件 $\delta E_{tot} / \delta c_{\mu i} = 0$ から、SCC-DFTBにおける一般化固有値問題が導かれる。

2.1 行列要素の修正#

Kohn-Sham方程式におけるハミルトニアン $H_{\mu\nu}$ は、電荷 $\Delta q$ に依存する項によって修正される。

\sum_{\nu} (H_{\mu\nu} - \epsilon_i S_{\mu\nu}) c_{\nu i} = 0

行列要素は以下のように与えられる。

H_{\mu\nu} = H_{\mu\nu}^0 + H_{\mu\nu}^{SCC}

H_{\mu\nu}^{SCC} = \frac{1}{2} S_{\mu\nu} \sum_{K} (\gamma_{AK} + \gamma_{BK}) \Delta q_K \quad (\mu \in A, \nu \in B)

ここで、 $H_{\mu\nu}^0$ はDFTB1で用いられる参照ハミルトニアン（Slater-Kosterテーブルから取得）、 $S_{\mu\nu}$ は重なり行列である。この式は、原子 $A$ と $B$ の軌道間の相互作用が、周囲の全原子 $K$ が作る静電ポテンシャル $(\sum \gamma \Delta q)$ によってシフトすることを示している。

2.2 マリケン電荷 (Mulliken Charge) の計算#

電荷変動 $\Delta q_A$ は、全電子数分布から中性原子の価電子数 $Z_A^{val}$ を引いたものである。

\Delta q_A = q_A - Z_A^{val}

q_A = \sum_{i}^{occ} n_i \sum_{\mu \in A} \sum_{\nu} c_{\mu i} c_{\nu i} S_{\mu\nu}

ここで $n_i$ は軌道の占有数（スピン縮退を含む場合は2.0等）である。

2.3 力 (Forces) の計算#

構造最適化やMDに必要な原子にかかる力 $\mathbf{F}_A = -\nabla_A E_{tot}$ は、以下の3項から成る。

Hellmann-Feynman力: バンドエネルギーの微分。
斥力項の微分: $\nabla E_{rep}$ 。
SCC力: $\gamma$ 行列の距離依存性に由来する項。 $\mathbf{F}_A^{SCC} = - \sum_{B} \Delta q_A \Delta q_B \frac{\partial \gamma_{AB}}{\partial R_{AB}} \frac{\mathbf{R}_{AB}}{R_{AB}} - \sum_{i}^{occ} n_i \sum_{\mu\nu} c_{\mu i} c_{\nu i} \frac{\partial H_{\mu\nu}^{SCC}}{\partial \mathbf{R}_A}$ 第3項は特に重要であり、電荷分布の変化に伴う力の寄与を表す。

3. アルゴリズムと擬似コード#

SCC-DFTBの実装には、電荷 $\Delta q$ を収束させるための反復（SCF）ループが必要となる。また、収束加速のためにBroyden混合法などが必須となる。

3.1 全体フローチャート#

初期化: Slater-Kosterファイルの読み込み、初期電荷 $\Delta q_A = 0$ の設定。
行列構築 (H0, S): 距離依存のSKテーブルから取得。
SCFループ:
- $\gamma$ 行列の計算。
- シフト項 $H^{SCC}$ の計算と $H = H^0 + H^{SCC}$ の更新。
- 対角化 $HC = SCE$ 。
- 新電荷 $q^{new}$ の計算（Mulliken解析）。
- 電荷混合 (Mixing): $\Delta q^{(n+1)} = (1-\alpha)\Delta q^{(n)} + \alpha \Delta q^{new}$ 。
- 収束判定: $|\Delta q^{new} - \Delta q^{(n)}| < \text{tol}$ なら脱出。
エネルギー・力の計算: 収束した密度行列を用いて最終値を算出。

3.2 Pythonによる擬似コード実装#

import numpy as np
from scipy.linalg import eigh

class SCC_DFTB_Solver:
    def __init__(self, atoms, sk_params, hubbard_u):
        """
        atoms: 原子オブジェクトリスト
        sk_params: Slater-Kosterパラメータ
        hubbard_u: 各元素のHubbard Uパラメータ辞書 {'C': 0.4, 'H': 0.2 ...}
        """
        self.atoms = atoms
        self.sk = sk_params
        self.U = hubbard_u
        self.n_atoms = len(atoms)
        self.charges = np.zeros(self.n_atoms) # Delta q
        
    def calculate_gamma_matrix(self, positions):
        """
        原子間距離とHubbard UからGamma行列を計算
        gamma_AB = 1 / sqrt(R^2 + 1/4 * (1/Ua + 1/Ub)^2 * tau) 近似
        """
        gamma = np.zeros((self.n_atoms, self.n_atoms))
        for i in range(self.n_atoms):
            for j in range(self.n_atoms):
                if i == j:
                    gamma[i, j] = self.U[self.atoms[i].element]
                else:
                    dist = np.linalg.norm(positions[i] - positions[j])
                    Ua = self.U[self.atoms[i].element]
                    Ub = self.U[self.atoms[j].element]
                    # Elstner (1998) Eq. 8 or similar approximation
                    # ここでは簡易的なOhno-Klopman型を示す
                    cutoff = 0.5 * (1.0/Ua + 1.0/Ub)
                    gamma[i, j] = 1.0 / np.sqrt(dist**2 + cutoff**2)
        return gamma

    def run_scf(self, max_iter=100, tol=1e-6, mix_param=0.2):
        # 1. Base Hamiltonian (H0) & Overlap (S) の構築
        H0, S = self.build_base_matrices()
        
        # 原子価電子数
        Z_val = np.array([self.sk.get_valence(a.element) for a in self.atoms])
        
        # 2. SCF Loop
        positions = np.array([a.pos for a in self.atoms])
        
        for iteration in range(max_iter):
            # Gamma行列の更新 (MDでない場合は固定で良い)
            gamma = self.calculate_gamma_matrix(positions)
            
            # SCC項によるHamiltonianの摂動
            # H_mu_nu += 0.5 * S_mu_nu * (pot_A + pot_B)
            # pot_A = sum_K gamma_AK * dq_K
            
            potentials = np.dot(gamma, self.charges) # 原子ごとの静電ポテンシャル
            
            H_scc = np.zeros_like(H0)
            orb_idx = 0
            atom_orb_map = [] # 軌道インデックスと原子の対応
            
            # インデックスの対応付け（前処理推奨）
            for i, atom in enumerate(self.atoms):
                n_orb = self.sk.get_orb_count(atom.element)
                for _ in range(n_orb):
                    atom_orb_map.append(i)
            
            # H_SCC の構築 (ベクトル化可能だが概念的にループ記述)
            for mu in range(H0.shape[0]):
                for nu in range(H0.shape[1]):
                    A = atom_orb_map[mu]
                    B = atom_orb_map[nu]
                    shift = 0.5 * (potentials[A] + potentials[B])
                    H_scc[mu, nu] = S[mu, nu] * shift
            
            H_total = H0 + H_scc
            
            # 対角化
            evals, evecs = eigh(H_total, S)
            
            # Mulliken電荷の計算
            # q_A = sum_occ n_i * sum_mu_in_A (c_mu_i * sum_nu c_nu_i * S_mu_nu)
            new_charges = np.zeros(self.n_atoms)
            
            # 密度行列 P_mu_nu = sum_occ n_i * c_mu_i * c_nu_i
            # 閉殻系を仮定 (2.0 electrons)
            n_occ = int(np.sum(Z_val) / 2)
            P = 2.0 * np.dot(evecs[:, :n_occ], evecs[:, :n_occ].T)
            
            PS = np.dot(P, S)
            
            current_orb = 0
            for i, atom in enumerate(self.atoms):
                n_orb = self.sk.get_orb_count(atom.element)
                # 原子の電子数 (Trace of PS block)
                q_elec = np.trace(PS[current_orb:current_orb+n_orb, :])
                # Delta q = q_core - q_elec (符号定義に注意: 本文数式は q_A - Z)
                # ここでは q_A は電子数ではなく正味の電荷とする
                new_charges[i] = Z_val[i] - q_elec
                current_orb += n_orb
                
            # 収束判定
            diff = np.linalg.norm(new_charges - self.charges)
            if diff < tol:
                print(f"Converged at iter {iteration}")
                self.charges = new_charges
                
                # エネルギー計算
                E_band = np.sum(evals[:n_occ]) * 2.0
                E_scc = 0.5 * np.dot(self.charges, np.dot(gamma, self.charges))
                E_rep = self.calculate_repulsive_energy()
                return E_band + E_scc + E_rep
            
            # Simple Mixing (Broyden推奨)
            self.charges = (1 - mix_param) * self.charges + mix_param * new_charges
            
        raise RuntimeError("SCF did not converge")

    def build_base_matrices(self):
        # DFTB1と同様にSKテーブルからH0, Sを作成
        # ... (前回の記事参照)
        return H0, S

    def calculate_repulsive_energy(self):
        # E_rep の計算
        # ...
        return E_rep

4. 実利的な成果と適用事例#

SCC-DFTBは、電荷移動の記述が可能になったことで、適用範囲が劇的に拡大した。

4.1 生体分子シミュレーション#

ペプチド結合の分極や、DNA塩基対間の水素結合は、静電相互作用が支配的である。SCC-DFTBはこれらの結合エネルギーをDFT（B3LYP等）に近い精度で再現できるため、タンパク質のフォールディングや酵素反応のQM/MM計算におけるQM領域ソルバーとして標準的に利用されるようになった。

4.2 酸化物および表面化学#

酸化チタン ( $TiO_2$ ) やシリカ ( $SiO_2$ ) などの金属酸化物は、金属と酸素間で大きな電荷移動を伴う。DFTB1では記述困難であったこれらの系のバンドギャップや表面吸着エネルギーが、SCC補正によって精度良く計算可能となった。これにより、色素増感太陽電池や触媒反応のシミュレーションが進展した。

4.3 有機反応機構#

極性溶媒中での反応や、電荷分離を伴う遷移状態を持つ反応（例：Sn2反応）において、電荷分布の変化に応じたエネルギー障壁の予測が可能となった。

5. 結論と展望#

SCC-DFTB（DFTB2）は、DFTB1の「高速性」を維持しつつ、密度揺らぎの2次項を取り入れることで「化学的反応性」と「静電相互作用」の記述能力を獲得した。数理的には、多重極展開の単極子近似と、Hubbardパラメータを用いた $\gamma$ 関数の導入が成功の鍵であった。現在では、さらに3次項を含めたDFTB3が開発され、水素結合のプロトン移動や高帯電系における精度が向上しているが、SCC-DFTBはその基礎となる標準モデルとして、物質科学から生命科学まで幅広い分野で利用され続けている。

参考文献#

M. Elstner, D. Porezag, G. Jungnickel, J. Elsner, M. Haugk, Th. Frauenheim, S. Suhai, and G. Seifert, “Self-consistent-charge density-functional tight-binding method for simulations of complex materials properties”, Phys. Rev. B 58, 7260 (1998).
T. Frauenheim, G. Seifert, M. Elstner, Z. Hajnal, G. Jungnickel, D. Porezag, S. Suhai, and R. Scholz, “A self-consistent charge density-functional based tight-binding method for predictive materials simulations in physics, chemistry and biology”, Phys. Status Solidi B 217, 41 (2000).