【DFT】DFT-D（分散力補正密度汎関数法）の数理と歴史：第一原理パラメータ化によるファンデルワールス相互作用の記述

最終更新：2025-12-30

注意: この記事はAIによって自動生成されたものです。正確な情報については、必ず引用元の原著論文をご確認ください。

序論：局所近似の限界と分散力補正の必要性#

密度汎関数法（DFT）は、その計算コストと精度のバランスから、現代の量子化学計算および固体物理計算において最も広く利用されている電子状態計算手法である。特に、一般化勾配近似（GGA）や混成汎関数（Hybrid Functional）の登場により、共有結合の結合エネルギーや分子構造の予測精度は飛躍的に向上した。

しかし、標準的なKohn-Sham DFT、特に局所密度近似（LDA）やGGA、および通常の混成汎関数には、原理的かつ重大な欠陥が存在する。それは、長距離電子相関（Long-range Electron Correlation）、すなわち**ロンドン分散力（London Dispersion Force）**の記述能力の欠如である。ロンドン分散力は、離れた位置にある電子密度の揺らぎが誘起する瞬間的な双極子相互作用に起因し、原子間距離 $R$ に対して漸近的に $-R^{-6}$ で減衰する引力ポテンシャルを持つ。標準的な交換相関汎関数は、電子密度の局所的な情報（ $\rho(\mathbf{r})$ ）やその勾配（ $\nabla \rho(\mathbf{r})$ ）に依存する「半局所的（Semilocal）」なモデルに基づいているため、空間的に離れた領域間の相関効果を記述することができず、分散相互作用エネルギーは指数関数的に減衰してしまう。

この欠陥は、希ガス二量体、ファンデルワールス錯体、生体分子のスタッキング相互作用、分子結晶の凝集エネルギーなどの計算において、結合エネルギーの過小評価（あるいは非結合）や不適切な平衡構造の予測という形で顕在化する。

この問題に対し、Stefan Grimmeらは、Kohn-Shamエネルギーに対して事後的に（a posteriori）古典的な分散力ポテンシャルを加算するアプローチ、すなわちDFT-D（Dispersion-corrected DFT）法を体系化した。特に2010年に発表されたDFT-D3法は、第一原理計算に基づくパラメータ決定と、原子の化学的環境（配位数）への依存性を導入することで、経験的補正の域を超えた物理的整合性を獲得し、計算化学におけるデファクトスタンダードとなった。

本稿では、Grimmeらの2010年の論文 “A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu” [1] に基づき、DFT-D3法の数理的導出過程、歴史的背景、およびその実利的な成果について詳細に解説する。

1. 歴史的背景：D1からD3への進化#

DFT-D法の基本概念は、全エネルギー $E_{DFT-D}$ を、Kohn-Sham DFTのエネルギー $E_{KS-DFT}$ と分散力補正エネルギー $E_{disp}$ の和として表現することにある。

E_{DFT-D} = E_{KS-DFT} + E_{disp}

この $E_{disp}$ の設計において、Grimmeらは段階的な改良を行ってきた。

1.1 DFT-D1 (2004) と DFT-D2 (2006)#

初期のDFT-D1（2004）およびその改良版であるDFT-D2（2006）では、分散エネルギーは以下のような古典的な対ポテンシャル和として記述された。

E_{disp}^{(D2)} = - s_6 \sum_{i<j} \frac{C_6^{ij}}{R_{ij}^6} f_{damp}(R_{ij})

ここで、 $C_6^{ij}$ は原子対 $i, j$ の分散係数、 $f_{damp}$ は短距離での発散を防ぐための減衰関数である。 DFT-D2の特徴は、原子ごとの分散係数 $C_6^i$ を原子番号のみに依存する定数として定義し、原子対係数を幾何平均 $C_6^{ij} \approx \sqrt{C_6^i C_6^j}$ で近似した点にある。これは実装が容易であり、多くの系で劇的な改善をもたらしたが、**「分子内の原子の状態変化（混成状態や酸化数など）を反映できない」**という物理的な限界があった。例えば、sp3炭素とsp2炭素は分極率が異なるため、本来異なる分散係数を持つべきであるが、D2では同一視されていた。

1.2 DFT-D3 (2010) の革新#

2010年に提案されたDFT-D3は、この問題を解決するために開発された。D3では、原子対の分散係数 $C_6^{ij}$ が、その原子の周囲の幾何学的環境（配位数）に依存する関数として再定義された。また、パラメータの決定プロセスにおいて、経験的なフィッティングを極力排除し、時間依存密度汎関数法（TDDFT）による第一原理計算値をベースとする手法が確立された。これにより、水素（H）からプルトニウム（Pu）までの94元素に対する一貫したパラメータ化が実現された。

2. 数理的背景：DFT-D3の導出#

DFT-D3における分散力エネルギー $E_{disp}^{(D3)}$ は、2体相互作用項 $E^{(2)}$ と3体相互作用項 $E^{(3)}$ の和として定義される。

E_{disp}^{(D3)} = E^{(2)} + E^{(3)}

ここでは、補正の中核をなす2体項の導出とその物理的背景に重点を置いて解説する。

2.1 2体相互作用エネルギーの定式化#

2体相互作用項は、双極子-双極子相互作用（ $R^{-6}$ ）だけでなく、より高次の双極子-四重極子相互作用（ $R^{-8}$ ）を含めた形式で記述される。

E^{(2)} = - \sum_{AB} \sum_{n=6,8} s_n \frac{C_n^{AB}}{r_{AB}^n} f_{damp,n}(r_{AB})

ここで、総和は全ての原子対 $AB$ について取られる。 $s_n$ は汎関数依存のスケーリング因子、 $C_n^{AB}$ は $n$ 次の分散係数、 $r_{AB}$ は原子間距離、 $f_{damp,n}$ は減衰関数である。

2.2 $C_6$ 係数の第一原理的導出 (Casimir-Polderの式)#

DFT-D3の最大の数理的特徴は、 $C_6$ 係数の決定方法にある。原子または分子AとBの間の分散係数 $C_6^{AB}$ は、Casimir-Polderの式を用いて、それぞれの動的分極率 $\alpha(i\omega)$ の虚振動数積分として厳密に定義される。

C_6^{AB} = \frac{3}{\pi} \int_0^{\infty} \alpha^A(i\omega) \alpha^B(i\omega) \, d\omega

Grimmeらは、この積分を分子中のすべての原子対に対して都度計算するのはコスト的に現実的ではないため、以下の手順で効率的な近似法を構築した。

参照系の計算: 各元素に対して、様々な配位状態や酸化状態を模した参照分子（主に水素化物 $A_mH_n$ ）を設定し、TDDFTを用いて動的分極率 $\alpha(i\omega)$ を計算する。ここから、原子ごとの参照係数 $C_{6,ref}^{AB}(CN_A, CN_B)$ を算出する。これは、原子AとBが特定の配位数（Coordination Number: CN）にあるときの分散係数である。
原子環境の補間: 実際の分子中の原子の分散係数 $C_6^{AB}$ は、事前計算された参照係数の線形結合（補間）として表現される。
$C_6^{AB}(CN_A, CN_B) = \frac{L}{k} \sum_{i}^{N_A} \sum_{j}^{N_B} W_i^A W_j^B C_{6,ref}^{AB}(i, j)$
(※ここで $L$ は定数、 $k$ は重みの正規化因子である)
ここで重要となるのが、原子の環境を記述する配位数 (Fractional Coordination Number, $CN$ ) の定義である。DFT-D3では、原子Aの配位数 $CN_A$ を以下の数式で定義する。
$CN_A = \sum_{B \neq A}^{N_{atoms}} \frac{1}{1 + \exp\left( -k_1 (k_2 \frac{R_{A,cov} + R_{B,cov}}{r_{AB}} - 1) \right)}$
ここで、 $R_{cov}$ は共有結合半径、 $k_1, k_2$ は定数である。この関数は、原子間距離 $r_{AB}$ が共有結合半径の和に近い場合に寄与が大きくなるようなカットオフ関数として機能する。この $CN_A$ の値に基づき、その原子がどのような化学結合状態にあるか（例：炭素原子なら、アルカンのような飽和状態か、グラファイトのような不飽和状態か）を判別し、適切な参照 $C_6$ 値を補間する。これにより、**「幾何構造の変化に応じた分散係数の動的な変化」**を低コストで取り入れることに成功した。

2.3 高次項 $C_8$ の再帰的導出#

$R^{-8}$ 項の係数 $C_8^{AB}$ は、別途計算するのではなく、 $C_6^{AB}$ から再帰的に導出される。Grimmeは、多重極展開の原子内平均に基づく関係式を用いた。

C_8^{AB} = 3 C_6^{AB} \sqrt{Q^A Q^B}

ここで、 $Q^A$ は原子Aの特定の期待値（ $s^{2}$ や $s^{4}$ に関連する量）から導かれるパラメータであり、原子ごとに事前に決定されている。これにより、新たなパラメータを増やすことなく高次の相互作用を取り込んでいる。

2.4 減衰関数 (Damping Function) の設計#

分散項 $C_n / r^n$ は $r \to 0$ で発散するため、短距離領域ではゼロに収束させる必要がある。また、短距離領域はDFTの交換相関汎関数が（不完全ながら）ある程度の相関効果を記述している領域でもあるため、分散項との「二重カウント（Double Counting）」を避ける必要がある。 2010年の原著論文で採用された標準的な減衰関数（Zero-damping）は以下の通りである。

f_{damp,n}(r_{AB}) = \frac{1}{1 + 6 (r_{AB} / (s_{r,n} R_0^{AB}))^{-\alpha_n}}

ここで、 $R_0^{AB}$ は原子対のファンデルワールス半径の和から決定されるカットオフ半径、 $s_{r,n}$ は汎関数ごとに調整されるスケーリングパラメータである。この関数は、長距離（ $r \gg R_0$ ）では1に近づき、短距離（ $r \ll R_0$ ）では $r^n$ に比例して減衰するため、エネルギーの発散を防ぐ。

汎関数依存性: DFT-D3法において、物理定数である $C_6, C_8$ は汎関数に依存しない普遍的な値である。一方で、短距離での接続を調整するパラメータ（ $s_6, s_8, s_{r,6}, s_{r,8}$ ）は、組み合わせるDFT汎関数（B3LYP, PBE, TPSSなど）ごとに最適化される。例えば、分散力を全く記述できないB3LYPでは $s_6=1.0$ となるが、ある程度の中距離相関を持つ汎関数では $s_6$ が調整される場合がある（ただしD3では通常 $s_6=1$ 固定で $s_8$ で調整することが多い）。

2.5 3体相互作用項 (ATM項)#

DFT-D3では、Axilrod-Teller-Muto (ATM) の式に基づく3体相互作用項 $E^{(3)}$ もオプションとして導入された。

E^{(3)} = \sum_{ABC} C_9^{ABC} \frac{1 + 3 \cos\theta_a \cos\theta_b \cos\theta_c}{r_{AB}^3 r_{BC}^3 r_{CA}^3} f_{damp}^{(3)}

これは、巨大な系や結晶中において、多体効果による分散力の遮蔽や増強を記述するために重要となるが、総エネルギーへの寄与は2体項に比べて小さい（通常5-10%程度）。

3. 実利的な成果と検証#

Grimmeらは、原子番号1(H)から94(Pu)までの全元素に対するパラメータ化を行い、広範なベンチマークテスト（GMTKN30データベースなど）を通じてDFT-D3の性能を検証した。

3.1 汎用性と精度の向上#

DFT-D3は、既存のDFT-D2と比較して、特に以下の点で顕著な改善を示した。

金属系・有機金属錯体: D2では適用が難しかった遷移金属や重元素を含む系においても、配位数依存性を取り入れたことで適切な分散補正が可能となった。
化学反応: 反応中心の配位数が変化するプロセスにおいて、D3は分散係数の変化を追跡できるため、反応熱や活性化エネルギーの記述精度が向上した。
軽元素分子の異性体: 直鎖アルカンと分岐アルカンのような、原子組成が同じで構造（原子密度）が異なる系のエネルギー差を正確に再現できるようになった。

3.2 標準偏差による評価#

S22データセット（非共有結合相互作用のベンチマーク）において、DFT-D3を適用したB3LYPやPBEなどの汎関数は、CCSD(T)レベルの参照値に対し、平均絶対偏差（MAD）で0.5 kcal/mol以下の精度を達成することが示された。これは、補正なしのDFTが数kcal/mol以上の誤差（あるいは結合の予測失敗）を生じることと比較すると、劇的な改善である。

3.3 計算コスト#

DFT-D3の計算コストは、Kohn-Sham SCF計算のコストと比較して無視できるほど小さい。エネルギー計算に必要な時間は系のアトム数 $N$ に対して $N^2$ でスケールするが、係数は非常に小さく、さらにカットオフ導入により実質 $O(N)$ で計算可能である。この「低コスト・高精度」という特性が、DFT-D3の爆発的な普及を後押しした。

4. プログラム出力：DFT-D3エネルギー計算のアルゴリズム#

以下に、DFT-D3（2体項、Zero-damping）のエネルギー計算プロセスを模したPythonコードを示す。これは実際の第一原理計算コードの一部として組み込まれるロジックの簡略化版である。

import numpy as np

class DFTD3_Mockup:
    """
    DFT-D3 (Zero-damping) エネルギー計算の簡易実装例
    注意: 実際のC6係数は事前計算された巨大なデータベースから
    配位数(CN)に基づいて補間されるが、ここでは簡略化のため
    固定値および仮想的な関数を使用する。
    """

    def __init__(self, functional_params):
        # 汎関数依存のパラメータ (例: B3LYP-D3 用)
        self.s6 = functional_params.get('s6', 1.0)
        self.s8 = functional_params.get('s8', 1.703)
        self.sr6 = functional_params.get('sr6', 1.261)
        self.sr8 = 1.0 # 通常固定
        
        # 原子固有のパラメータ (簡略値)
        # 実際にはGrimmeの論文のSupplementary Informationにある
        # 配位数依存の補間テーブルを使用する。
        self.C6_ref = {
            'C': 30.0, # 仮想値 (a.u.)
            'H': 1.5,
            'O': 15.0
        }
        self.R0_ref = { # vdW半径 (Angstrom -> Bohr変換が必要)
            'C': 1.7,
            'H': 1.2,
            'O': 1.5
        }
        self.k_cov = 1.889725989 # Angstrom to Bohr

    def get_cn(self, atoms, coords):
        """
        配位数 (Coordination Number) の計算
        CN_A = sum_B 1 / (1 + exp(-k1 * (k2 * (Ra+Rb)/r_AB - 1)))
        """
        k1 = 16.0
        k2 = 4.0/3.0
        n_atoms = len(atoms)
        cns = np.zeros(n_atoms)
        
        for i in range(n_atoms):
            for j in range(n_atoms):
                if i == j: continue
                
                r_ij = np.linalg.norm(coords[i] - coords[j])
                # 共有結合半径の和 (簡略化: vdW半径をスケールして代用)
                r_cov_sum = (self.R0_ref[atoms[i]] + self.R0_ref[atoms[j]]) * 0.8 
                
                term = 1.0 / (1.0 + np.exp(-k1 * (k2 * (r_cov_sum / r_ij) - 1.0)))
                cns[i] += term
        return cns

    def interpolate_c6(self, atom_i, atom_j, cn_i, cn_j):
        """
        配位数に基づくC6係数の補間
        実際には参照データ点からのガウス重み付け補間などを行う。
        ここでは配位数によるスケーリングを模したダミー関数とする。
        """
        c6_pure = np.sqrt(self.C6_ref[atom_i] * self.C6_ref[atom_j])
        
        # 配位数が大きいほど分散係数が変化する効果を模擬
        env_factor = 1.0 - 0.05 * (cn_i + cn_j) 
        return c6_pure * env_factor

    def compute_energy(self, atoms, coords):
        """
        D3エネルギーの計算
        E_disp = - sum_ij [ s6 * C6/r^6 * fd6 + s8 * C8/r^8 * fd8 ]
        """
        E_disp = 0.0
        n_atoms = len(atoms)
        
        # 配位数の計算
        cns = self.get_cn(atoms, coords)
        
        for i in range(n_atoms):
            for j in range(i + 1, n_atoms):
                r_ij = np.linalg.norm(coords[i] - coords[j])
                r_ij_bohr = r_ij * self.k_cov
                
                # C6, C8係数の取得
                c6 = self.interpolate_c6(atoms[i], atoms[j], cns[i], cns[j])
                c8 = 3.0 * c6 * np.sqrt(1.0 * 1.0) # Qパラメータを1.0と仮定した簡易式
                
                # カットオフ半径 (R0_AB)
                r0_ab = self.R0_ref[atoms[i]] + self.R0_ref[atoms[j]]
                r0_ab_bohr = r0_ab * self.k_cov

                # 減衰関数 (Zero-damping)
                # fd_n = 1 / (1 + 6 * (r / (sr_n * R0))^alpha)
                # alpha6 = 14, alpha8 = 16 (Typical values)
                
                fd6 = 1.0 / (1.0 + 6.0 * (r_ij_bohr / (self.sr6 * r0_ab_bohr))**14)
                fd8 = 1.0 / (1.0 + 6.0 * (r_ij_bohr / (self.s8 * r0_ab_bohr))**16) # s8 parameter acts as scaling usually
                # Note: In standard D3, s8 is the prefactor, but for cutoff usually sr8 is 1.0.
                # Here we follow the energy formula structure strictly.
                
                term6 = (self.s6 * c6 / (r_ij_bohr**6)) * fd6
                term8 = (self.s8 * c8 / (r_ij_bohr**8)) * fd8 # Note: s8 is strictly energy scaling factor
                
                E_disp -= (term6 + term8)
        
        return E_disp

# サンプル実行
# メタン二量体のような系を想定 (単位: Angstrom)
atoms = ['C', 'H', 'H', 'H', 'H', 'C', 'H', 'H', 'H', 'H']
# 簡略化座標 (適当な距離に配置)
coords = np.array([
    [0.0, 0.0, 0.0], [0.6, 0.6, 0.6], [-0.6, -0.6, 0.6], [-0.6, 0.6, -0.6], [0.6, -0.6, -0.6],
    [4.0, 0.0, 0.0], [4.6, 0.6, 0.6], [3.4, -0.6, 0.6], [3.4, 0.6, -0.6], [4.6, -0.6, -0.6]
])

# B3LYP-D3 パラメータ
params = {'s6': 1.0, 's8': 1.703, 'sr6': 1.261}

d3_calc = DFTD3_Mockup(params)
energy = d3_calc.compute_energy(atoms, coords)

print(f"Computed DFT-D3 Dispersion Energy: {energy:.6f} Hartree")
print(f"Note: This is a demonstration value using mock parameters.")

5. 結論と展望#

Stefan Grimmeらによって開発されたDFT-D3法は、密度汎関数法における長距離相関の問題に対し、物理的に透明かつ実用的な解決策を提供した。その数理的基盤は、Casimir-Polderの式による第一原理的な分散係数の導出と、配位数を用いた動的な環境効果の取り込みにある。これにより、DFT-D3は単なる経験的補正（Empirical Correction）の枠を超え、**「漸近的に正しい（Asymptotically Correct）」**理論としての地位を確立した。

D3法の成功は、その後に続くD4法（帯電状態への依存性を導入）や、非局所相関汎関数（vdW-DF）との統合的な理解へとつながっている。しかし、実装の簡便さ、パラメータの網羅性（94元素対応）、そして計算コストの低さにおいて、DFT-D3は依然として計算化学における分散力補正の「金字塔」であり続けている。

有機分子の結晶構造予測から、巨大なタンパク質のドッキングシミュレーション、さらには表面吸着現象の解析に至るまで、現代の計算化学においてDFT-D3の恩恵を受けていない分野を探すことは困難であると言える。

参考文献#

S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, “A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu”, J. Chem. Phys. 132, 154104 (2010).