# How DFT-D3 incorporates coordination number (CN) into C6

I'm reading Grimme's DFT-D3 paper but really puzzled by how the $\mathrm{CN}$ is considered in the calculation of $C_\mathrm{6}$ coefficients. (Stefan Grimme, Jens Antony, Stephan Ehrlich, and Helge Krieg. J. Chem. Phys. 2010, 132 (15), 154104.)

According to Fig. 5 in the paper, it seems the $\ce{AB}$ pair is put in some reference system and then based on the precomputed $C_\mathrm{6}$ values in the reference system we could get the $C_\mathrm{6}$ in the environment.

Especially, I'm confused by the fact in the figure that $\mathrm{CN}^\ce{A}$ and $\mathrm{CN}^\ce{B}$ are in different dimension, what does it mean? I will quote the paper quite a bit, but I'll try and summarize a bit after the quotes and equations. You might want to start at the bottom and work backward, a lot of this is just for later reference.

First, some background. Equation 3:

$$E^{(2)} = \sum_{AB} \sum_{n=6,8,10,\dots} s_{n} \frac{C_{n}^{AB}}{r_{AB}^{n}} f_{d,n}(r_{AB})$$

Here, the first sum is over all atom pairs in the system, $C_{n}^{AB}$ denotes the averaged (isotropic) $n$th order dispersion coefficient (orders $n=6,8,10,\dots$) for atom pair $AB$ (...)

So $A$ and $B$ are nuclear indices, and together specify an arbitrary pair.

Although the $C_{6}^{AB}$ values can be compute easily for any pair of free atoms by using Eq. (5) [the Casimir-Polder integral] in principle, this would lead to a rather inconsistent treatment of dispersion in and between molecules.

For reference, here is equation 5:

$$C_{6}^{AB} = \frac{3}{\pi} \int_{0}^{\infty} \alpha^{A}(i\omega) \alpha^{B}(i\omega) \, d\omega$$

Instead they calculate the integrals over pairs of hydrides, subtracting the effect of molecular hydrogen (equation 10), in order to capture the effect of being bonded:

$$C_{6}^{AB} = \frac{3}{\pi} \int_{0}^{\infty} d\omega \frac{1}{m} \left[ \alpha^{A_{m}H_{n}}(i\omega) - \frac{n}{2} \alpha^{H_{2}}(i\omega) \right] \times \frac{1}{k} \left[ \alpha^{A_{k}H_{l}}(i\omega) - \frac{l}{2} \alpha^{H_{2}}(i\omega) \right]$$

At first sight this new approach seems to be a disadvantage because it leads to reference molecule dependent (ambiguous) coefficients. However, as will be discussed in detail below (...), it opens a route to system [coordination number (CN)] dependent "atomic" $C_6$(CN) coefficients. (...) We propose here a radically new "geometric" approach that is based on the concept of a fractional coordination number. (...) We base our ansatz on a fractional CN for an atom $A$, that is a sum of a "counting" function over all atoms in the system,

Equation 15:

$$\text{CN}^{A} = \sum_{B\neq A}^{N_{at}} \frac{1}{1 + e^{ -k_1 ( k_2 (R_{A,\text{cov}} + R_{B,\text{cov}}) / r_{AB} - 1) }}$$

where $R_{A,\text{cov}}$ is a scaled covalent (single-bond) radius of atom $A$ (and analogously for $R_{B,\text{cov}}$).

So the CN$^A$ is some bond distance-dependent scaling factor for the specific atom $A$ under the influence of all other atoms $\{B\}$. Here is how they incorporate their coordination number scaling into the dispersion coefficients. Following emphasis mine:

Turning back to the dispersion coefficients we propose the following general approach. For each element in the Periodic Table, at least one reference molecule is used to compute $\alpha(i\omega)$ values in Eq. (10). (...) If an element is commonly investigated, and/or often found in different chemical environments, we suggest to use several representatives for which values are computed. In the case of carbon for example, ethyne, ethene, and ethane can be used for CNs between 2 and 4. Values computed for $\ce{C-H}$ and the carbon atom complete the set. For the atom of interest in any of the reference systems also the CN is calculated and stored. Using Eq. (10) the $C_{6,\text{ref}}^{AB}$(CN$^A$,CN$^B$) values are computed for this pair with the two atoms in their specific environments. These serve as supporting points in an interpolation procedure to derive the coefficient for any combination of fractional CN values. Because we base our approach on pair-specific coefficients $C_{6}^{AB}$, a two-dimensional (2D) interpolation scheme must be used. After extensive testing we propose a simple Gaussian-distance ($L$) weighted average,

Equation 16:

$$C_{6}^{AB}(\text{CN}^{A},\text{CN}^{B}) = \frac{ \sum_{i}^{N_A} \sum_{j}^{N_B} C_{6,\text{ref}}^{AB}(\text{CN}_{i}^{A},\text{CN}_{j}^{B}) \times L_{ij} }{ \sum_{i}^{N_A} \sum_{j}^{N_B} L_{ij} }, \\ L_{ij} = e^{ -k_3 [ (\text{CN}^A - \text{CN}_i^A)^2 + (\text{CN}^B - \text{CN}_j^B)^2 ] },$$

where $N_A$ and $N_B$ are the number of supporting points (= number of reference molecules) for atoms $A$ and $B$, respectively. The CN$^A$ and CN$^B$ are coordination numbers for the atom pair $AB$ in the system of interest, and the CN$_i^A$ and CN$_j^B$ are those for the two reference systems $i$ and $j$, for which $C_{6,\text{ref}}^{AB}(\text{CN}_{i}^{A},\text{CN}_{j}^{B})$ is the precomputed value. (...) For most of the elements only two reference values are necessary to cover the typical bonding situations. For "important" elements such as carbon, however, more values (5 in this case) seem appropriate. Currently we have computed in total 227 reference systems for the elements up to $Z=94$ which leads to about $2.6\times10^4$ difference values for $C_{6,\text{ref}}^{AB}(\text{CN}_{i}^{A},\text{CN}_{j}^{B})$. They are computed once and read at the beginning of any computation.

Conceptually, here is how a calculation would work:

1. Your SCF has converged and it's time to calculate the -D3 correction. This table is read in from a file into some data structure so that all the $C_{6,\text{ref}}^{AB}(\text{CN}_{i}^{A},\text{CN}_{j}^{B})$ can be looked up easily.

2. Loop over the atom pairs $AB$ of your system as part of equation 3. For atom $A$, calculate its coordination number CN$^A$ using equation 15. Do the same thing for atom $B$. In the equation, $B\neq A$ means all atoms other than $A$, not atom $B$ in the pair!

3. For this pair, the $C_{n}^{AB}$ in equation 3 is really $C_{6}^{AB}(\text{CN}^{A},\text{CN}^{B})$, which requires evaluating equation 16. What equation 16 does is interpolate between those precomputed $C_{6,\text{ref}}^{AB}(\text{CN}_{i}^{A},\text{CN}_{j}^{B})$ for a series of coordination numbers. Because there are two CNs, the interpolation is on a 2D grid, rather than a line or some higher dimension. The interpolation function is chosen to be the $L_{ij}$. From the above quote, each $i,j$ sum is only about 2-5 points, so this is very fast; even the $\ce{C-C}$ grid will only total $5\times5=25$ points.

This is exactly correct. It is done by interpolation on a grid of reference values, where the coordinates/points on the grid are the coordination numbers, and there is a grid for every possible $AB$ pair.
$\tiny{\text{Phew.}}$