xnn.common.models.d4#

DFT-D4 London dispersion: a charge-dependent correction for any xnn model.

Implements the D4 model of Caldeweyher et al., J. Chem. Phys. 150, 154122 (2019) (doi:10.1063/1.5090222) in pure PyTorch, as an additive energy term that combines with every xnn model family – a GNN, SchNet, a descriptor network or a classical force field – and deploys through the same channels (PyTorch, TorchScript, ASE, LAMMPS). In the paper’s terms this is the default D4 model: EEQ partial charges, BJ (rational) damping for the two-body term and the approximate Axilrod-Teller-Muto three-body term (bj-eeq-atm).

The model, equation by equation (atomic units; the paper’s numbering)#

  • Coordination number (eq 6): an error-function count of neighbors, weighted by a Pauling-electronegativity factor,

    CN_A = sum_B delta^EN_AB/2 (1 + erf(-k0 (R_AB - R^cov_AB) / R^cov_AB)), delta^EN_AB = k1 exp(-(abs(EN_A - EN_B) + k2)^2 / k3),

    with k0 = 7.5, k1 = 4.10451, k2 = 19.08857, k3 = 2 * 11.28174^2 and R^cov_AB the sum of the (4/3-scaled) Pyykko covalent radii.

  • EEQ partial charges (eqs 11-16): the electronegativity-equilibration charges of Gaussian charge densities of width a_A, obtained from the linear system [[A, 1], [1^T, 0]] [q, lambda] = [X, q_tot] with A_AA = J_A + 2 gamma_AA / sqrt(pi) (gamma_AB = (a_A^2 + a_B^2)^-1/2), A_AB = erf(gamma_AB R_AB) / R_AB and X_A = -EN_A + kappa_A sqrt(mCN_A), where mCN is the plain (electronegativity-free) error-function CN of eq 14, softly capped at 8. For periodic structures the 1/r matrix is Ewald-summed.

  • Charge scaling (eqs 2-4): every reference polarizability is scaled by zeta(z, z_ref) = exp(beta1 [1 - exp(gamma_A [1 - z_ref / z])]) with z = Z_eff + q (beta1 = 3, gamma_A twice the element’s chemical hardness).

  • Reference polarizabilities (eq 5): the atom-in-molecule dynamic polarizabilities are partitioned out of the TD-DFT (PBE38/daug-def2-QZVP) polarizabilities of the reference systems A_m X_n by subtracting the charge-scaled contribution of the X_n atoms.

  • Gaussian CN weighting (eqs 7-8): alpha_A(i omega) = sum_ref W_A,ref alpha_A,ref(i omega) with normalized weights sum_j^{N^s} exp(-beta2 j (CN_A - CN_A,ref)^2) (beta2 = 6; N^s grows where reference CNs cluster, fig 4).

  • Casimir-Polder integration (eqs 1, 9): C6^AB = 3/pi int alpha_A(i w) alpha_B(i w) dw on the fixed 23-point trapezoid grid.

  • Two-body energy (eqs 18-21): E = -sum_AB sum_n=6,8 s_n C_n^AB / (R_AB^n + R_0^n) with R_0 = a1 sqrt(C8/C6) + a2 and C8 = 3 C6 sqrt(Q_A Q_B) (Q the <r^4>/<r^2> factors).

  • Three-body ATM energy (eqs 22-27): E = s9 sum_ABC C9^ABC (3 cos cos cos + 1) / (R_AB R_BC R_CA)^3 / (1 + 6 (R_0^ABC / R_ABC)^16) with C9 = sqrt(C6 C6 C6) built from neutral (q = 0) polarizabilities.

Fidelity#

An independent implementation, written from the paper and from the observed behavior of the reference code dftd4 (not copied, and never imported here). The conventions that the paper leaves to the code are reproduced so that energies, forces, virials, charges, coordination numbers, C6 coefficients and polarizabilities match dftd4 to floating-point precision (tests/test_d4.py and examples/fidelity_checks/d4_verification.ipynb): the N^s bookkeeping of eq 8 and its fallback to the highest-CN reference when every Gaussian weight underflows, the soft CN cap in the EEQ model (and its absence in the D4 CN), the effective nuclear charges of the reference systems, the q = 0 C6 coefficients of the ATM term, the real-space cutoffs (60 / 40 / 30 / 25 bohr for pairs / triples / CN / EEQ-CN) with the optional quintic switching windows, and the Ewald conventions of the periodic EEQ (the automatic splitting parameter, the fixed real- and reciprocal-space windows, the Wigner-Seitz image averaging).

The element and reference data (d4_reference.npz next to this file) are the numerical values published with the D4 method (TD-DFT polarizabilities, reference coordination numbers and charges, element constants), extracted from the dftd4 / multicharge / mctc-lib sources by tools/build_d4_reference.py.

Units and cutoffs#

Like every xnn model the public interfaces take Angstrom and return eV; internally the model works in atomic units with CODATA-2018 conversion factors. The D4 cutoffs (default: the upstream 60 / 40 / 30 / 25 bohr) set the neighbor-list radius of D4Dispersion (model.cutoff), so a wrapped short-range model receives only the edges within its own cutoff. For condensed-phase MLIP training, shorter D4 cutoffs (10-15 Angstrom) with a switching window are the sensible choice; the defaults reproduce dftd4 exactly.

Functions

c6_matrix(alpha_iw)

Pairwise C6 from dynamic polarizabilities (paper eqs 1, 9).

d4_options_from_extra(extra)

Pick the DFTD4 keyword arguments out of a config extra dict.

Classes

D4Dispersion([model])

DFT-D4 dispersion as an xnn potential, standalone or wrapped around a model.

DFTD4([s6, s8, a1, a2, s9, alp, ga, gc, wf, ...])

The D4 dispersion model as a TorchScript-compatible energy evaluator.