【計算化学】非電荷自己無撞着密度汎関数タイトバインディング法（Non-SCC-DFTB）の数理的定式化と実装論：第一原理に基づく強束縛近似の導出

最終更新：2025-12-31

注意: この記事はAIによって自動生成されたものです。正確な情報やパラメータの詳細については、必ず引用元の原著論文（Porezag et al., Phys. Rev. B 51, 12947 (1995) 等）をご確認ください。

序論：第一原理と経験論の架け橋としてのDFTB#

物質の電子状態を計算する手法において、密度汎関数理論（Density Functional Theory: DFT）は現代の標準的なツールとしての地位を確立している。しかし、DFT（特にKohn-Sham法）の計算コストは系に含まれる原子数 $N$ に対して $O(N^3)$ で増大し、プレファクター（係数）も大きいため、数千原子を超える系や長時間の分子動力学（MD）シミュレーションへの適用には依然として困難が伴う。一方で、経験的パラメータを用いる古典的分子動力学法は高速であるが、化学反応（結合の組み換え）や量子力学的効果を記述することは原理的に不可能である。

この間隙を埋める手法として発展してきたのが、密度汎関数タイトバインディング法（Density Functional Tight-Binding: DFTB） である。DFTBは、DFTの総エネルギー表現を出発点とし、電子密度の揺らぎに関する系統的な展開近似を導入することで導出される。これにより、DFTの物理的厳密さを一定程度保持しつつ、半経験的手法に匹敵する計算速度を実現している。

本稿では、DFTBの最も基本的な形式である DFTB1、あるいは Non-SCC-DFTB（Non-Self-Consistent-Charge DFTB）と呼ばれる手法に焦点を当てる。これは、系内の電荷移動に伴う静電相互作用の自己無撞着な緩和を考慮しない（あるいは0次近似として扱う）モデルである。後に開発されたSCC-DFTB（DFTB2）やDFTB3の基礎となるこの手法について、その数理的導出、ハミルトニアンの構成法、および斥力項の取り扱いを、プログラム実装が可能なレベルまで詳細に掘り下げて解説する。

1. 歴史的背景：経験的パラメータから第一原理ベースへ#

1.1 タイトバインディング（TB）モデルの変遷#

タイトバインディング（強束縛）近似は、固体の電子状態を記述する最も初期のモデルの一つである。BlochやSlater、Kosterらによって定式化された古典的なTBモデル（Slater-Koster法）は、ハミルトニアン行列要素を実験バンド構造へのフィッティングパラメータとして扱っていた。これは簡便で強力であったが、パラメータの移送性（Transferability）—すなわち、ある環境で決定したパラメータが別の化学環境で通用するかどうか—に問題を抱えていた。

1.2 ab initio Tight-Bindingへの希求#

1980年代後半から1990年代にかけて、アンダーソン局在や欠陥構造などの複雑な系を扱うために、パラメータを実験値ではなく第一原理計算（DFT）から決定しようとする試みがなされた。SankeyやNiklewskiらは、擬ポテンシャルと局在基底を用いた手法を提案したが、これらは依然として多中心積分を含むなど計算コストが高かった。

1.3 PorezagらによるDFTBの確立 (1995年)#

1995年、D. Porezag、Th. Frauenheim、Th. Köhler、G. Seifert、R. Kaschnerらは、DFTに基づきながら、二中心近似（Two-Center Approximation）を徹底することで、ハミルトニアンおよび重なり行列要素を事前に計算・テーブル化する手法を確立した。これが現在の標準的なDFTB（Non-SCC）の誕生である。彼らの手法は、原子間の相互作用ポテンシャル（斥力項）をDFTの総エネルギーとバンド構造エネルギーの差として定義し、これをペアポテンシャルとしてフィッティングすることで、構造最適化やMDシミュレーションにおいて極めて高い精度と安定性を実現した。特に炭素材料（フラーレン、ナノチューブ、ダイヤモンド、アモルファスカーボン）において顕著な成功を収め、その後のナノテクノロジー研究を支える基盤技術の一つとなった。

2. 数理的背景：DFTからの展開とDFTB1の導出#

DFTBの理論的骨格は、DFTの総エネルギー汎関数を、参照密度 $\rho_0$ 周りで2次までTaylor展開することによって得られる。ここではその導出過程をステップごとに追う。

2.1 DFT総エネルギーの展開#

Kohn-Sham DFTにおける総エネルギー $E^{KS}$ は、電子密度 $\rho(\mathbf{r})$ の汎関数として以下のように書ける。

E^{KS}[\rho] = T_s[\rho] + E_{ext}[\rho] + E_H[\rho] + E_{xc}[\rho] + E_{nn}

ここで、

$T_s[\rho]$ : 相互作用のない運動エネルギー
$E_{ext}[\rho]$ : 外部ポテンシャル（原子核-電子相互作用）エネルギー $\int V_{ext}(\mathbf{r})\rho(\mathbf{r}) d\mathbf{r}$
$E_H[\rho]$ : Hartree（古典的クーロン）エネルギー $\frac{1}{2} \iint \frac{\rho(\mathbf{r})\rho(\mathbf{r}')}{|\mathbf{r}-\mathbf{r}'|} d\mathbf{r} d\mathbf{r}'$
$E_{xc}[\rho]$ : 交換相関エネルギー
$E_{nn}$ : 原子核間反発エネルギー

ここで、真の電子密度 $\rho(\mathbf{r})$ を、参照密度 $\rho_0(\mathbf{r})$ とそこからの微小な揺らぎ $\delta\rho(\mathbf{r})$ の和として表現する。

\rho(\mathbf{r}) = \rho_0(\mathbf{r}) + \delta\rho(\mathbf{r})

通常、参照密度 $\rho_0(\mathbf{r})$ は、系を構成する中性原子の密度 $\rho_0^A(\mathbf{r})$ の重ね合わせ（Superposition of Atomic Densities）として定義される。

\rho_0(\mathbf{r}) = \sum_{A} \rho_0^A(\mathbf{r} - \mathbf{R}_A)

この分解を用いてエネルギー汎関数を展開すると、以下のようになる。

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

2.2 各項の物理的解釈と近似#

0次項と1次項の変形#

展開の0次項と1次項を整理し、Kohn-Sham軌道 $\psi_i$ （占有数 $n_i$ ）を用いて書き下すと、DFTBのエネルギー式が得られる。ここで運動エネルギー項などの取り扱いに注意が必要だが、最終的にエネルギーは「バンド構造エネルギー（Band Structure Energy）」と「斥力エネルギー（Repulsive Energy）」の和として近似される。

E_{total}^{DFTB1} \approx \sum_{i}^{occ} n_i \epsilon_i + E_{rep}[\rho_0]

ここで、 $\epsilon_i$ は以下の有効ハミルトニアン $\hat{H}_{eff}[\rho_0]$ の固有値である。

\hat{H}_{eff}[\rho_0] = -\frac{1}{2}\nabla^2 + V_{eff}[\rho_0] = -\frac{1}{2}\nabla^2 + V_{ext} + V_H[\rho_0] + V_{xc}[\rho_0]

重要な点は、このハミルトニアンが参照密度 $\rho_0$ のみに依存しており、自己無撞着に求める必要がないことである。これが計算高速化の鍵である。

2次項の無視（DFTB1の定義）#

展開の2次項 $E^2$ は、密度の揺らぎ $\delta\rho$ に起因するCoulomb相互作用や交換相関の2次応答を含んでいる。

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

DFTB1（Non-SCC-DFTB）では、この2次項を無視する、あるいは短距離の寄与として斥力項 $E_{rep}$ に含めてパラメータ化するという近似を行う。これにより、電荷の反復計算（SCCループ）が不要となり、一度の対角化のみで電子状態が決定される。この近似は、系内の電荷移動が小さい場合（ホモ核系や極性の低い共有結合系）において正当化される。

2.3 基底関数の構成とLCAO近似#

波動関数 $\psi_i$ を原子軌道（Atomic Orbitals: AO） $\phi_\mu$ の線形結合（LCAO）で展開する。

\psi_i(\mathbf{r}) = \sum_{\mu} c_{\mu i} \phi_\mu(\mathbf{r})

これにより、Kohn-Sham方程式は一般化固有値問題（Roothaan-Hall方程式類似）に帰着される。

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

最小基底（Minimal Basis）と閉じ込めポテンシャル#

DFTBでは、計算効率を最大化するために最小基底系（水素は1s、炭素は2s, 2pのみ等）を用いる。しかし、自由原子の軌道は空間的に広がりすぎており、固体内での局在性を表現するのに適さない。そこでPorezagらは、原子DFT計算において閉じ込めポテンシャル（Confinement Potential） $V_{conf}(r)$ を追加した修正Kohn-Sham方程式を解くことで、基底関数を縮退（収縮）させる手法を提案した。

\left[ -\frac{1}{2}\nabla^2 + V_{eff}^{atom}[\rho_0^A] + V_{conf}(r) \right] \phi_\mu = \epsilon_\mu^{atom} \phi_\mu

典型的には、以下の形式のポテンシャルが用いられる。

V_{conf}(r) = \left( \frac{r}{r_0} \right)^N

ここで $r_0$ は閉じ込め半径（通常原子共有結合半径の2倍程度）、 $N$ は整数（例：2など）である。この操作により、軌道および密度が空間的にコンパクトになり、多中心積分の無視が正当化されやすくなる。

2.4 二中心近似による行列要素の計算#

ハミルトニアン行列要素 $H_{\mu\nu}$ と重なり行列要素 $S_{\mu\nu}$ は以下のように計算される。

S_{\mu\nu} = \langle \phi_\mu | \phi_\nu \rangle

H_{\mu\nu} = \langle \phi_\mu | \hat{T} + V_{eff}[\rho_0] | \phi_\nu \rangle

ここで、有効ポテンシャル $V_{eff}[\rho_0]$ を原子ごとのポテンシャルの和で近似する（Three-center termsの無視）。

V_{eff}[\rho_0] \approx \sum_{K} V_{eff}^K[\rho_0^K]

さらに、 $\phi_\mu$ が原子A、 $\phi_\nu$ が原子Bに属する場合、ポテンシャルの和のうち原子AとB由来のもののみを残す**二中心近似（Two-Center Approximation）**を行う。

H_{\mu\nu} \approx \langle \phi_\mu | \hat{T} + V_{eff}^A + V_{eff}^B | \phi_\nu \rangle \quad (\mu \in A, \nu \in B)

この近似により、 $H_{\mu\nu}$ と $S_{\mu\nu}$ は原子間距離 $R_{AB}$ のみの関数となり、事前に計算してテーブル化（Slater-Kosterファイル）することが可能となる。これがDFTBの高速性の根源である。

3. 斥力ポテンシャル（Repulsive Potential）の定義と構築#

DFTB1の総エネルギーは以下のように表される。

E_{tot} = E_{band} + E_{rep} = \sum_{i}^{occ} \epsilon_i + \sum_{A<B} V_{rep}^{AB}(R_{AB})

ここで $E_{band}$ は行列対角化から得られる引力的な項であり、 $E_{rep}$ は二重カウントされたHartreeエネルギーや交換相関エネルギー、イオン間反発などを補正するための斥力項である。

3.1 斥力項の導出#

理論的には $E_{rep}$ は以下のように定義される。

E_{rep} = E_{tot}^{DFT} - E_{band}^{DFTB}

すなわち、DFTBで再現したいターゲットとなる高精度なDFT計算（通常はLDAやPBE汎関数を用いた全電子計算）を行い、その総エネルギー $E_{tot}^{DFT}$ から、DFTBのバンドエネルギー $E_{band}^{DFTB}$ を差し引いた残差として定義される。

3.2 ペアポテンシャル近似#

Non-SCC-DFTBでは、この $E_{rep}$ を短距離のペアポテンシャルの和として近似する。

E_{rep} \approx \sum_{A<B} U_{AB}(R_{AB})

この $U_{AB}(R)$ は、参照となる結晶構造や分子構造に対して結合距離 $R$ を変化させながらDFT計算とDFTBバンド計算を行い、その差分をスプライン関数や多項式でフィッティングすることで作成される。このプロセスにより、DFTBは参照としたDFTの構造特性（結合長、結合角、弾性定数など）を非常に高い精度で再現する。

4. 実装のためのアルゴリズム詳細#

ここでは、DFTB1を実装するために必要な具体的なステップと、主要部分のPython擬似コードを示す。

4.1 プログラムの全体フロー#

初期化:
- Slater-Koster (SK) テーブル（.skf ファイル等）の読み込み。
- 原子座標の読み込み。
Hamiltonian/Overlap行列の構築:
- 各原子ペアについて距離ベクトルを計算。
- SKテーブルからスプライン補間により $H_{\mu\nu}, S_{\mu\nu}$ を取得。
- Slater-Koster変換則に従い、軌道の向き（ $p_x, p_y, p_z$ 等）を考慮して行列要素を回転・充填。
対角化:
- 一般化固有値問題 $HC = SCE$ を解く（Cholesky分解等で標準固有値問題へ変換）。
- 固有値 $\epsilon_i$ を昇順にソートし、Fermi分布または占有数に従って $E_{band}$ を計算。
斥力エネルギー計算:
- 原子間距離に基づき、SKテーブル内の斥力スプラインから $E_{rep}$ を計算。
総エネルギーと力の出力:
- $E_{tot} = E_{band} + E_{rep}$ 。
- 原子にかかる力（Forces）の計算（Hellmann-Feynman力 + 基底関数微分項 + 斥力項微分）。

4.2 Slater-Koster変換の実装#

SKテーブルには、 $\sigma, \pi, \delta$ 結合に対応する基本積分値（例： $pp\sigma, pp\pi$ ）が距離の関数として格納されている。これを実際の座標系での行列要素に変換するには、方向余弦 $(l, m, n)$ を用いた変換式（Slater-Koster rules）を用いる。

変換式の例（s-p軌道間）#

原子A（位置 $\mathbf{0}$ ）のs軌道と、原子B（位置 $\mathbf{R}=(Rx, Ry, Rz)$ ）のp軌道間の行列要素 $E_{s, px}$ は以下のように表される。

E_{s, px} = l \cdot E_{sp\sigma}

ここで $l = Rx / |\mathbf{R}|$ である。

4.3 Pythonによる実装コード例（主要部）#

以下に、Hamiltonian行列を構築するクラスの概念的な実装を示す。

import numpy as np
from scipy.linalg import eigh

class DFTB1Solver:
    def __init__(self, atoms, sk_params):
        """
        atoms: 原子情報（座標、元素種）
        sk_params: Slater-Kosterパラメータを管理するオブジェクト
        """
        self.atoms = atoms
        self.sk_params = sk_params
        self.n_orbitals = sum([sk_params.get_orbital_count(a.element) for a in atoms])
        
    def build_matrices(self):
        """
        Hamiltonian (H) と Overlap (S) 行列を構築する
        """
        n = self.n_orbitals
        H = np.zeros((n, n))
        S = np.zeros((n, n))
        
        # 軌道インデックスの管理用
        orbital_indices = self._get_orbital_indices()
        
        for i, atom_i in enumerate(self.atoms):
            for j, atom_j in enumerate(self.atoms):
                # 対角ブロック（原子内）
                if i == j:
                    idx = orbital_indices[i]
                    # 原子の軌道エネルギー（onsite terms）をセット
                    # SKパラメータから取得（例: es, ep, ed）
                    energies = self.sk_params.get_onsite_energies(atom_i.element)
                    for k, e in enumerate(energies):
                        H[idx[k], idx[k]] = e
                        S[idx[k], idx[k]] = 1.0
                    continue
                
                # 非対角ブロック（原子間）
                vec_r = atom_j.pos - atom_i.pos
                dist = np.linalg.norm(vec_r)
                
                # カットオフ判定
                if dist > self.sk_params.cutoff:
                    continue
                
                # 方向余弦の計算
                l, m, n_dir = vec_r / dist
                
                # SKテーブルから基本積分値を取得 (Spline補間)
                # integrals = {'ss_sigma': val, 'sp_sigma': val, ...}
                sk_vals_H = self.sk_params.get_hamiltonian_integrals(atom_i.element, atom_j.element, dist)
                sk_vals_S = self.sk_params.get_overlap_integrals(atom_i.element, atom_j.element, dist)
                
                # Slater-Koster変換行列の生成
                block_H = self._slater_koster_transform(sk_vals_H, l, m, n_dir, atom_i.element, atom_j.element)
                block_S = self._slater_koster_transform(sk_vals_S, l, m, n_dir, atom_i.element, atom_j.element)
                
                # 行列への充填
                idx_i = orbital_indices[i]
                idx_j = orbital_indices[j]
                
                # numpyのブロードキャスト等を使って埋める
                # H[idx_i_start:idx_i_end, idx_j_start:idx_j_end] = block_H
                for r in range(len(idx_i)):
                    for c in range(len(idx_j)):
                        H[idx_i[r], idx_j[c]] = block_H[r, c]
                        S[idx_i[r], idx_j[c]] = block_S[r, c]
                        
        return H, S

    def _slater_koster_transform(self, vals, l, m, n, elem1, elem2):
        """
        SK積分値と方向余弦から行列ブロックを作成する
        例: s-s, s-p, p-p 相互作用の展開
        """
        # 軌道タイプに応じたサイズ（s:1, p:3, d:5）
        # ここでは簡単のため s, p 軌道のみのケースを示す
        # 行列サイズは (1+3) x (1+3) = 4x4 などを想定
        
        # 行列ブロックの初期化
        mat = np.zeros((4, 4)) 
        
        # s-s sigma
        mat[0, 0] = vals['ss_sigma']
        
        # s-p sigma (原子iのsと原子jのp)
        # E_s_px = l * E_sp_sigma
        mat[0, 1] = l * vals['sp_sigma']
        mat[0, 2] = m * vals['sp_sigma']
        mat[0, 3] = n * vals['sp_sigma']
        
        # p-s sigma (原子iのpと原子jのs)
        # 対称性に注意（符号反転など）
        mat[1, 0] = -l * vals['sp_sigma']
        mat[2, 0] = -m * vals['sp_sigma']
        mat[3, 0] = -n * vals['sp_sigma']
        
        # p-p ブロック
        # E_xx = l^2 * pp_sigma + (1-l^2) * pp_pi
        # E_xy = l*m * pp_sigma - l*m * pp_pi
        # ... (以下同様に展開)
        
        pp_sigma = vals['pp_sigma']
        pp_pi = vals['pp_pi']
        
        mat[1, 1] = l*l * pp_sigma + (1 - l*l) * pp_pi
        mat[1, 2] = l*m * pp_sigma - l*m * pp_pi
        # ... 他の要素も同様に記述
        
        return mat

    def calculate_energy(self):
        H, S = self.build_matrices()
        
        # 一般化固有値問題を解く
        eigvals, eigvecs = eigh(H, S)
        
        # バンドエネルギー計算（占有軌道の和）
        # 簡単のため電子数をカウントして下から詰める
        n_electrons = sum([self.sk_params.get_valence_electrons(a.element) for a in self.atoms])
        n_occ = n_electrons // 2 # スピン非考慮または閉殻系
        
        E_band = 2.0 * np.sum(eigvals[:n_occ])
        
        # 斥力エネルギー
        E_rep = 0.0
        for i, atom_i in enumerate(self.atoms):
            for j in range(i+1, len(self.atoms)):
                dist = np.linalg.norm(atom_i.pos - atom_j.pos)
                E_rep += self.sk_params.get_repulsive_potential(atom_i.element, atom_j.element, dist)
                
        return E_band + E_rep

4.4 斥力項のフィッティング戦略#

実装において最も職人芸的な部分が斥力項の作成である。Porezagらの手法では、以下の手順で作成される。様々な結合距離を持つ参照系（例：炭素ならダイヤモンド、グラファイト、C2分子など）について、DFT計算で $E_{tot}^{DFT}(R)$ を求める。同じ構造について、斥力項なしのDFTBバンドエネルギー $E_{band}^{DFTB}(R)$ を計算する。差分 $E_{rep}(R) = E_{tot}^{DFT}(R) - E_{band}^{DFTB}(R)$ をプロットする。この差分曲線が、短距離では強い正の値（斥力）、長距離ではゼロに収束する滑らかな関数になることを確認し、これを多項式スプライン（通常はカットオフ距離を設けた3次または5次スプライン）でフィッティングする。

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

DFTB1（Non-SCC）は、特にホモ核共有結合系において、驚異的な計算効率と良好な精度を両立させた。

5.1 炭素材料科学におけるブレイクスルー#

1990年代後半、フラーレンやカーボンナノチューブの研究が爆発的に進展していた時期において、DFTBは最適なツールであった。フラーレンの生成過程: 数千個の炭素原子を含む高温MDシミュレーションが可能となり、フラーレンの形成メカニズム（アニーリング過程）の解明に貢献した。アモルファスカーボン: ダイヤモンドライクカーボン（DLC）などの密度やsp2/sp3比率のシミュレーションにおいて、実験値を良く再現した。これは古典ポテンシャル（TersoffやBrenner）では記述が難しく、かつab initioでは計算コストが高すぎる領域であった。

5.2 大規模生体分子への適用#

電荷の偏りが少ない疎水性タンパク質や脂質膜などの系において、DFTBは構造最適化の初期構造生成や、長時間MDのドライビングフォースとして利用された。

5.3 計算効率#

Non-SCC-DFTBの計算コストは、行列対角化（ $O(N^3)$ ）が支配的であるが、SCCループがないため、SCC-DFTBに比べて通常5〜10倍高速である。さらに、行列要素の計算（ $O(N^2)$ ）や斥力項計算（ $O(N)$ ）は極めて高速であるため、並列化効率も高い。数万原子の系に対しても、線形スケーリング法（Order-N法）を組み合わせることで適用可能である。

6. 結論と展望D#

FTB1（Non-SCC-DFTB）は、第一原理計算であるDFTから出発し、系統的な近似（2次項の無視、二中心近似、最小基底）を経ることで導出される、極めて合理的な電子状態計算手法である。そのハミルトニアンは経験的フィッティングではなく、原子のDFT計算から直接的に構築されるため、物理的な透明性が高い。本手法は、電荷移動が顕著な系（イオン結晶や極性分子反応）には不向きであるという明確な限界を持つが、その限界は後に2次項を取り入れた SCC-DFTB（DFTB2）や、3次項を含む DFTB3 によって克服されていく。しかし、炭素材料や共有結合性結晶などの電荷揺らぎの小さい系においては、DFTB1は現在でも「高速かつそれなりな精度を持つ量子力学的計算手法」として、材料設計の第一線で活用され続けている。

参考文献#

D. Porezag, Th. Frauenheim, Th. Köhler, G. Seifert, and R. Kaschner, “Construction of tight-binding-like potentials on the basis of density-functional theory: Application to carbon”, Phys. Rev. B 51, 12947 (1995).
G. Seifert, D. Porezag, and Th. Frauenheim, “Calculations of molecules, clusters, and solids with a simplified LCAO-DFT-LDA scheme”, Int. J. Quantum Chem. 58, 185 (1996).
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).