【DFT】DFT-D4モデルの数理と実装：原子部分電荷依存性を導入した次世代分散力補正の全貌

最終更新：2025-12-30

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

序論：幾何構造依存から電子状態依存へ#

密度汎関数法（DFT）、特に一般化勾配近似（GGA）やハイブリッド汎関数は、その計算コストの低さと適度な精度により、現代量子化学計算の主力ツールとなっている。しかし、これらの標準的な汎関数は、電子密度の局所的な性質に基づいているため、長距離の電子相関、すなわち**ロンドン分散力（London dispersion force）**を記述できないという本質的な欠陥を抱えている。

この問題を解決するために、Stefan Grimmeらは半古典的な事後補正項を加えるDFT-Dシリーズ（D2, D3）を開発し、計算化学のデファクトスタンダードとしての地位を確立した。特に2010年に発表されたDFT-D3は、原子の配位数（Coordination Number: CN）に基づいて分散係数 $C_6$ を内挿することで、周囲の化学環境による原子の大きさの変化を模倣することに成功した。

しかし、D3モデルには理論的な限界が存在した。それは、**「幾何構造が同じであれば、電子状態が変化しても分散力パラメータは変化しない」**という点である。例えば、遷移金属錯体において中心金属の酸化数が変化した場合や、有機分子における電荷移動錯体のような系では、原子間の距離（幾何構造）がほとんど変化しないにもかかわらず、電子雲の広がり（分極率）は劇的に変化する。配位数のみに依存するD3モデルでは、この「電子的な呼吸」を捉えることができず、酸化還元反応や帯電系において誤差を生む原因となっていた。

2017年および2019年、Eike CaldeweyherとStefan Grimmeらは、この問題を解決するために原子部分電荷（Atomic Partial Charge）への依存性を明示的に取り入れたDFT-D4モデルを提案した。D4モデルは、幾何構造から原子電荷を高速に算出し、その電荷量に応じて原子の動的分極率をスケーリングさせることで、第一原理計算コストを維持したまま、電子状態の変化に追随する分散力補正を実現した。

本稿では、DFT-D4モデルの数理的構造について、基礎となる電気陰性度等化法（EEQ）から分散係数の導出、エネルギー計算に至るまで、実装可能なレベルで詳細に解説する。

1. 数理的背景：電荷依存的分散係数の導出#

DFT-D4の核心は、原子の動的分極率 $\alpha(i\omega)$ を、その原子の部分電荷 $q$ の関数として定式化した点にある。

1.1 全エネルギー表式#

DFT-D4における分散力補正エネルギー $E_{\text{disp}}$ は、以下の2体項（2-body）と多体項（Many-body）の和で表される。多体項については、D3と同様にAxilrod-Teller-Muto（ATM）項で近似されるが、係数の決定方法が異なる。

E_{\text{disp}}^{\text{D4}} = E^{(2)}_{\text{disp}} + E^{(3)}_{\text{disp}}

2体相互作用項は、原子対 $A, B$ 間の距離 $R_{AB}$ を用いて以下の級数展開で記述される。

E^{(2)}_{\text{disp}} = - \sum_{A>B} \sum_{n=6,8,...} s_n \frac{C_n^{AB}}{R_{AB}^n + f_{\text{damp}}(R_{AB})}

ここで、 $C_n^{AB}$ は $n$ 次の分散係数、 $s_n$ は汎関数ごとのスケーリング因子、 $f_{\text{damp}}$ は短距離での発散を防ぐための減衰関数である。 D4モデルの最大の革新は、この $C_6^{AB}$ 係数が原子 $A, B$ の部分電荷 $q_A, q_B$ に依存して決定される点にある。

1.2 電気陰性度等化法（EEQ）による電荷算出#

D4モデルは「事後補正」としての高速性を維持するため、DFT計算から得られる電子密度（Mulliken電荷など）を使用しない。代わりに、幾何構造のみから原子電荷を算出する電気陰性度等化法（Electronegativity Equilibration: EEQ）、具体的にはGasteiger-Marsili法を拡張したモデルを採用している。

系の全エネルギー $E_{\text{EEQ}}$ を原子電荷 $q_A$ の2次関数として定義する。

E_{\text{EEQ}}(\{q\}) = \sum_A \left( \chi_A q_A + \frac{1}{2} \eta_A q_A^2 \right) + \sum_{A \neq B} q_A q_B J_{AB}

ここで、

$\chi_A$ : 原子 $A$ の電気陰性度
$\eta_A$ : 原子 $A$ の化学的硬度（Chemical Hardness）
$J_{AB}$ : 原子間のクーロン相互作用項

このエネルギーを、全電荷 $Q_{\text{tot}} = \sum q_A$ の制約条件下で最小化する（ラグランジュの未定乗数法）。これにより、以下の線形方程式系が得られる。

\begin{pmatrix} \eta_1 & J_{12} & \dots & J_{1N} & -1 \\ J_{21} & \eta_2 & \dots & J_{2N} & -1 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ J_{N1} & J_{N2} & \dots & \eta_N & -1 \\ 1 & 1 & \dots & 1 & 0 \end{pmatrix} \begin{pmatrix} q_1 \\ q_2 \\ \vdots \\ q_N \\ \mu \end{pmatrix} = \begin{pmatrix} -\chi_1 \\ -\chi_2 \\ \vdots \\ -\chi_N \\ Q_{\text{tot}} \end{pmatrix}

パラメータの詳細設定: D4モデルでは、 $\eta_A$ と $J_{AB}$ に独自の工夫が施されている。

硬度 $\eta_A$ の配位数依存性: 原子の硬度は一定ではなく、周囲の配位数（Fractional Coordination Number: CN）に依存して変化する。
$\eta_A = \eta_A^{\text{ref}} + \Delta \eta_A \times \text{CN}_A$
ここで、 $\text{CN}_A$ はD3モデルと同様の定義を用いる。
$\text{CN}_A = \sum_{B \neq A} \frac{1}{1 + e^{-k_0 (R_{\text{cov}, A} + R_{\text{cov}, B}) / R_{AB} - 1}}$
このCN依存性により、共有結合環境にある原子の硬さの変化を取り込む。
クーロン項 $J_{AB}$ の遮蔽: 短距離でのクーロン発散を防ぐため、Ohno-Klopman型の相互作用ポテンシャルを採用する。
$J_{AB} = \frac{1}{\sqrt{R_{AB}^2 + (\gamma_{AB})^{-2}}}$
ここで $\gamma_{AB}$ は原子対の硬度パラメータに関連する項である。

この連立一次方程式を解くことで、系全体の幾何構造を反映した平衡電荷 $\{q_A\}$ が得られる。

1.3 動的分極率のスケーリングとCasimir-Polder積分#

得られた部分電荷 $q_A$ を用いて、原子 $A$ の虚振動数依存動的分極率 $\alpha_A(i\omega)$ を決定する。D4モデルでは、中性原子の参照分極率 $\alpha_A^{\text{ref}}(i\omega)$ を、電荷に応じてスケーリングする手法をとる。

\alpha_A(i\omega, q_A) = \alpha_A^{\text{ref}}(i\omega) \times \zeta_A(q_A)

スケーリング関数 $\zeta_A(q_A)$ : スケーリング関数は、原子内の有効電子数の変化と、それに伴う分極率の変化を記述する。単純な線形近似ではなく、指数関数的な減衰を考慮した以下の形式が採用されている。

\zeta_A(q_A) = \exp\left( -\frac{3}{2} \ln \left( \frac{Z_{\text{eff}, A} + \delta q_A}{Z_{\text{eff}, A}} \right) \right)

あるいは、論文によっては以下のようなPadé近似形式等価な記述や、より直接的に有効核電荷 $Z_{\text{eff}}$ と部分電荷の関係から導かれる形式で表現される。本質的には、**「正に帯電すると電子雲が収縮して分極率が下がり、負に帯電すると膨張して分極率が上がる」**という物理的描像を再現する。ここで $Z_{\text{eff}, A}$ は原子 $A$ の有効核電荷である。

分散係数 $C_6$ の計算: 原子対 $A, B$ の $C_6$ 係数は、Casimir-Polder積分によって算出される。

C_6^{AB}(q_A, q_B) = \frac{3}{\pi} \int_0^{\infty} \alpha_A(i\omega, q_A) \alpha_B(i\omega, q_B) d\omega

D4モデルでは、TD-DFT等の高精度計算によって事前に算出された各元素の参照分極率データ $\alpha_A^{\text{ref}}(i\omega)$ をデータベースとして保持している。実行時には、EEQで得た電荷を用いてこの分極率をスケーリングし、数値積分を行うことで $C_6$ をオンザフライで決定する。これにより、任意の電荷状態に対応した分散係数が得られる。

1.4 高次項と多体項#

$C_8$ 係数: D3と同様に、再帰的な関係式を用いて $C_6$ から推定する。
$C_8^{AB} = 3 C_6^{AB} \sqrt{Q_A Q_B}$
ここで $Q_A, Q_B$ は原子固有のパラメータである。
ATM 3体項: 3原子 $A, B, C$ 間の分散係数 $C_9^{ABC}$ も、スケーリングされた動的分極率を用いた一般化Casimir-Polder積分により算出される。
$C_9^{ABC} \approx \frac{3}{\pi} \int_0^{\infty} \alpha_A(i\omega) \alpha_B(i\omega) \alpha_C(i\omega) d\omega$
簡易的には、幾何平均近似を用いて $C_6$ から推定されることもあるが、D4の完全な実装では動的分極率の積の積分が基本となる。

1.5 減衰関数 (Damping Function)#

短距離での物理的でない発散を防ぎ、DFTの相関部分との二重カウントを避けるために減衰関数が必要である。D4では主に**Rational Damping（Becke-Johnson型減衰の変種）**が使用される。

E_{\text{disp}}^{(2)} = - \sum_{A>B} s_6 \frac{C_6^{AB}}{R_{AB}^6 + (a_1 R_0^{AB} + a_2)^6} - \sum_{A>B} s_8 \frac{C_8^{AB}}{R_{AB}^8 + (a_1 R_0^{AB} + a_2)^8}

ここで $R_0^{AB} = \sqrt{C_8^{AB}/C_6^{AB}}$ は原子対のカットオフ半径、 $a_1, a_2$ は汎関数ごとに最適化されるパラメータ、 $s_6, s_8$ はスケーリング因子である。

2. 実装ロジック：DFT-D4計算のアルゴリズム#

上記の数理的背景に基づき、DFT-D4エネルギーおよび勾配を計算するプログラムの構造を疑似コード（Python風）で示す。このロジックは、実際のライブラリ（dftd4）の動作原理を反映している。

import numpy as np
from scipy.integrate import simps

class DFT_D4:
    def __init__(self, atoms, coordinates, functional_params):
        """
        atoms: List of atomic numbers
        coordinates: Nx3 array of atomic positions
        functional_params: Dict containing s6, s8, a1, a2 for the chosen DFT functional
        """
        self.atoms = atoms
        self.coords = coordinates
        self.params = functional_params
        
        # Load Reference Data (Pseudo-code)
        # alpha_ref[element] gives frequency-dependent polarizability points
        # param_ref[element] gives chi, eta, Z_eff for EEQ
        self.ref_data = load_d4_reference_data()

    def compute_eeq_charges(self):
        """
        Step 1: Calculate Partial Charges via EEQ Model
        Solves the linear equation system A * x = b
        """
        N = len(self.atoms)
        A_mat = np.zeros((N + 1, N + 1))
        b_vec = np.zeros(N + 1)
        
        # Calculate Coordination Numbers (CN) for hardness scaling
        cn = self.compute_coordination_numbers()

        for i in range(N):
            Zi = self.atoms[i]
            # Hardness scaling: eta = eta_ref + delta_eta * CN
            eta_i = self.ref_data[Zi].eta + self.ref_data[Zi].d_eta * cn[i]
            chi_i = self.ref_data[Zi].chi
            
            A_mat[i, i] = eta_i
            A_mat[i, N] = -1.0
            A_mat[N, i] = 1.0
            b_vec[i] = -chi_i
            
            for j in range(i + 1, N):
                R_ij = np.linalg.norm(self.coords[i] - self.coords[j])
                # Ohno-Klopman Coulomb integral
                gamma_ij = self.compute_gamma(self.atoms[i], self.atoms[j])
                J_ij = 1.0 / np.sqrt(R_ij**2 + gamma_ij**(-2))
                
                A_mat[i, j] = J_ij
                A_mat[j, i] = J_ij

        b_vec[N] = 0.0 # Total charge constraint (assuming neutral system)
        
        # Solve linear system
        x = np.linalg.solve(A_mat, b_vec)
        return x[:N] # Return charges q

    def compute_dispersion_energy(self):
        """
        Step 2 & 3: Interpolate Polarizabilities and Compute Energy
        """
        q = self.compute_eeq_charges()
        N = len(self.atoms)
        E_disp = 0.0
        
        # Frequency integration grid (usually pre-defined Gauss-Legendre)
        omega_grid = self.ref_data['omega_grid'] 
        
        # Pre-compute scaled polarizabilities for each atom
        alpha_scaled = []
        for i in range(N):
            Zi = self.atoms[i]
            alpha_ref_i = self.ref_data[Zi].alpha_iw # Array over omega
            Z_eff = self.ref_data[Zi].Z_eff
            
            # Scaling function zeta(q)
            # Note: Simplification. Actual D4 uses a more specific rational/exp form
            # relating to the number of effective electrons.
            zeta = np.exp(-3.0/2.0 * np.log((Z_eff + q[i]) / Z_eff))
            
            alpha_scaled.append(alpha_ref_i * zeta)

        # Pairwise summation
        for i in range(N):
            for j in range(i + 1, N):
                R_ij = np.linalg.norm(self.coords[i] - self.coords[j])
                
                # Casimir-Polder Integration
                # Integrate[ alpha_i(iw) * alpha_j(iw) ] * (3/pi)
                integrand = alpha_scaled[i] * alpha_scaled[j]
                C6_ij = (3.0 / np.pi) * np.trapz(integrand, omega_grid)
                
                # C8 calculation (Recursive)
                Q_i = self.ref_data[self.atoms[i]].Q_param
                Q_j = self.ref_data[self.atoms[j]].Q_param
                C8_ij = 3.0 * C6_ij * np.sqrt(Q_i * Q_j)
                
                # Damping Parameters
                R0_ij = np.sqrt(C8_ij / C6_ij)
                a1 = self.params['a1']
                a2 = self.params['a2']
                s6 = self.params['s6']
                s8 = self.params['s8']
                
                # Energy Terms with Rational Damping
                term6 = s6 * C6_ij / (R_ij**6 + (a1 * R0_ij + a2)**6)
                term8 = s8 * C8_ij / (R_ij**8 + (a1 * R0_ij + a2)**8)
                
                E_disp -= (term6 + term8)
        
        # (Optional) Add 3-body ATM term here
        
        return E_disp

    def compute_gamma(self, Z1, Z2):
        # Helper for Coulomb integral parameter
        # Typically derived from chemical hardness or atomic radii
        return 1.0 # Placeholder
        
    def compute_coordination_numbers(self):
        # Compute fractional coordination numbers
        # Based on covalent radii and interatomic distances
        return np.zeros(len(self.atoms)) # Placeholder

実装上の注意点#

参照データ: 実装には、各元素ごとの参照分極率 $\alpha^{\text{ref}}(i\omega)$ 、電気陰性度、硬度、有効核電荷などのパラメータセットが不可欠である。これらは論文のSupplementary Informationや公式実装（GitHub上の dftd4）から取得する必要がある。勾

配（力）の計算: 構造最適化を行うためには、エネルギー $E_{\text{disp}}$ の原子座標微分 $\nabla E_{\text{disp}}$ が必要である。D4では電荷 $q$ が座標に依存するため、 $\partial C_6 / \partial q \times \partial q / \partial R$ という項（Charge-Cycling項）が発生する。これは連立方程式の解の微分を含むため、実装の複雑度はD3より格段に高い。

3. 実利的な成果とD3との比較#

D4モデルの導入により、どのような化学的課題が解決されたのか、原著論文のベンチマーク結果を基に解説する。

3.1 酸化還元状態と金属錯体#

D3モデルの最大の弱点は、原子価の変化を認識できないことである。

事例: チタン錯体などの遷移金属系において、酸化数が変化すると中心金属のイオン半径（電子雲の広がり）は大きく変化する。

D3の結果: 配位環境が変わらなければ $C_6$ 係数は一定であり、酸化によるイオン収縮（分極率の低下）を反映できず、分散力を過大評価する傾向があった。

D4の成果: EEQにより中心金属の正電荷の増加を検知し、自動的に分極率をスケーリングダウンさせる。これにより、反応熱や構造パラメータの精度が劇的に向上した。MOR41（Metal-Organic Reactions）ベンチマークセットにおいて、D4はD3と比較して誤差を有意に減少させている。

3.2 立体障害と「有機分子の呼吸」#

巨大な有機分子や超分子系において、電荷移動は局所的な分極率に影響を与える。

事例: S66やGMTKN55などの標準ベンチマークセット。

成果: 中性分子が主体のセットではD3とD4の差は小さいが、双極子モーメントが大きい系や、分子内電荷移動が顕著な系ではD4が優位である。特に、D4はD3が抱えていた「過剰な引力」を適正化する傾向があり、長鎖アルカンのフォールディングやタンパク質の二次構造安定性などの記述において、より物理的に妥当な結果を与える。

###3.3 周期境界条件への適用

2019年の拡張により、D4は周期的境界条件（PBC）下の結晶系にも適用可能となった。EEQモデルにおけるクーロン相互作用の長距離総和（Ewald法など）を実装することで、金属有機構造体（MOF）や分子結晶の格子エネルギー計算において、帯電状態を考慮した高精度な分散力補正が可能となった。

結論#

DFT-D4モデルは、Grimmeの分散力補正スキームにおける「電子状態依存性」の欠如という最後のミッシングリンクを埋めるものである。電気陰性度等化法（EEQ）と動的分極率のスケーリングを組み合わせることで、計算コストをDFT計算本体に比べて無視できるレベル（ $N^3$ スケーリングだが係数が非常に小さい）に抑えつつ、酸化状態や電荷移動に感応する物理的に洗練されたモデルを構築した。数理的には、Casimir-Polder積分をオンザフライで行うアプローチへと進化したことで、パラメータのブラックボックス性が低減し、より第一原理に近い補正となっている。このモデルは、特に遷移金属触媒や光励起状態、帯電した表面吸着系など、電子状態の変化が本質的な役割を果たす複雑な化学系において、標準的なDFT計算の信頼性を底上げする不可欠な要素技術である。

参考文献#

E. Caldeweyher, C. Bannwarth, and S. Grimme, “Extension of the D3 dispersion coefficient model”, J. Chem. Phys. 147, 034112 (2017).
E. Caldeweyher, S. Ehlert, A. Hansen, H. Neugebauer, S. Spicher, C. Bannwarth, and S. Grimme, “A generally applicable atomic-charge dependent London dispersion correction”, J. Chem. Phys. 150, 154122 (2019).

補足: 「有機分子の呼吸（Electronic Breathing）」の表現について#

「有機分子の呼吸（Electronic Breathing）」とは、比喩的な表現であり、学術的には**「電子状態（部分電荷）の変化に伴う、原子の分極率（電子雲の広がりやすさ）の動的な膨張・収縮」**を指します。DFT-D4の文脈において、この概念は非常に重要です。D3とD4の違いを明確にしながら解説します。

物理的な意味：電子雲の「膨張」と「収縮」原子が電子を受け取ったり放出したりすると、その原子の「分散力に対する強さ（分極率）」は変化します。

「吸う」状態（電子過剰 / 負電荷 $\delta^-$ ）:原子が電子を受け取ると、電子間の反発が増え、有効核電荷による束縛が相対的に弱まります。結果として、電子雲は膨張し、外部電場に対して変形しやすくなります（分極率が増大）。分散力（引力）は強くなります。

「吐く」状態（電子不足 / 正電荷 $\delta^+$ ）:原子から電子が引き抜かれると、残った電子は原子核に強く引きつけられます。結果として、電子雲は収縮し、硬くなります（分極率が減少）。分散力（引力）は弱くなります。

このように、分子内の電荷移動に応じて、各原子がまるで呼吸するかのように電子雲のサイズ（分散力の強さ）を変える現象を指しています。

D3モデルの限界（呼吸できないマネキン）

従来のDFT-D3モデルでは、分散係数 $C_6$ は「原子の幾何学的な配位数」だけで決まっていました。問題点: 結合距離が変わらない限り、原子が正に帯電しようが負に帯電しようが、D3モデルは「原子の大きさ（分極率）は同じ」とみなします。

比喩: 服（電荷）を着せ替えても、マネキンの体型（分散力）が変わらないようなものです。

結果: 電荷移動が起きている系（例：プッシュ-プル型の共役系有機分子や、分極したタンパク質）において、電子が抜けて「縮むべき」原子に対しても、中性原子と同じ強い分散力を割り当ててしまい、**引力の過大評価（Overbinding）**を引き起こしていました。

D4モデルの革新（呼吸する生体）

DFT-D4モデルは、原子部分電荷 $q$ に依存して分極率をスケーリングします。

解決策: $C_6(q)$ という関数により、正電荷の原子は分散力を弱め、負電荷の原子は分散力を強めるという補正をオンザフライで行います。

比喩: 息を吸えば胸が膨らみ、吐けば縮むように、電子状態に合わせて原子の「見かけの大きさ」が変化します。

結果: 「立体障害」と表現されるような混み合った領域において、局所的な電荷の偏り（分極）による反発や引力の変化を正確に捉えることができ、過剰な引力が補正され、正しい構造安定性が得られます。

まとめ#

「有機分子の呼吸」とは、幾何構造（骨格）が固定されていても、電子の移動によって原子の分極率がダイナミックに変化する現象のことです。D4モデルはこの「呼吸」を再現できるため、巨大分子のフォールディングや超分子集合体において、D3より物理的に妥当な結果を与えます。