Day 37: DDEC6 -- Theory

Now that the PRs for QTAIM are pretty much finalized, I am starting to digest the theories related to DDEC6 so that I can start implementing this soon. Calculating DDEC6 charges involve seven steps which include very specific procedures and coefficients that were optimized based on the method developers’ parameter survey.

Before going into the details mentioned in the original paper (doi:10.1039/c6ra04656h), the paper sets up a few mathematical conventions. It first lists each atom that constitute a molecule as A, and describes the unit lattice cell with L₁, L₂, and L₃. Lattice vectors associated with the setup is v₁, v₂, and v₃, whereas the location of the grid relative to the center of the atom is set as R_A. Here, methods like Bader’s QTAIM which partitions based on the charge density generally has a form of θ(r) = ρ(r) - Σ_{A, L} ρ_A(r_A). Here, A and L in the summation indicates each constituent atoms and grid respectively.

In this setup, the aim is to optimize ρ_A(r_A)/ρ(r) in a certain way. In vectorized charge partitioning, this quantity is equated with spherically symmetric atomic weighing factor w_A(r_A) / W(r), where W(r) = Σ_{A, L}w_A. In ESP charges, from the fact that the atomic charges can be approximated as point charges, W(r) = ρ(r), ρ_A(r_A) = w_A(r_A).

In Hirshfeld method, w is derived from reference densities of neutral atom. This is represented into an equation w_A^HD(r_A) = ρ_A^ref(r_A, q_A^ref=0).

Iterative Hirshfeld improves on Hirshfeld method and assigns q_A^ref = q_A in a self-consistent manner.

Iterative Stockholder method starts the partition from ρ(r) and satisfies the equation w_A^ISA(r_A) = ρ_A^avg(r_A).

In DDEC, first concept that is introduced is ‘conditioned charge densities’. This quantity is defined as ρ_A^conditioned(r_A) = ρ_A^generic(r_A) * <ρ(r) / ρ^generic(r)>_rA. Here, <>_rA denotes spherical averaging and generic refers to any charge densities that were obtained/prepared for use in AIM methods.

In the paper, charge densities that were conditioned one time from reference ion/atoms are indicated with the symbol Y_A^avg(r_A). Then, in DDEC6, w_A(r_A) = (ρ_A^ref(r_A))^1/5(ρ_A^avg(r_A))^4/5.

DDEC6 charges are calculated in largely 7 steps. In the first and second step, ion charge values are obtained. In the first step, following calculations are made:

1. Start from the charge densities of the molecule and set w_A^{1, Stock}(r_A) = ρ_A^ref(r_A, 0). 0 in the parenthesis indicates that the no ionic assumptions are made yet.
2. Then calculate N_A^{1, Stock} = &oint; (w_A^{1, Stock}(r_A) / Σ_{A, L} w_A^{1, Stock}(r))
3. Set w_A^{1, Loc}(r_A) = ρ_A^ref(r_A, 0)⁴ and carry out step 2.
4. Calculate q^Stock = Z_A - N_A^{1, Stock} and likewise for localized weight.
5. Set q_A^{1, ref} = (1/3)q_A^{1, Stock} + (2/3)q_A^{1, Loc}.

In the second step of DDEC6 routine, first step (which is listed above) is repeated simply with q_A^{1, ref} instead of zero.

As seen in this post, calculating DDEC6 charges involve quite a lot of conditioning and optimizing steps based on the numerical behavior that the researchers were able to find through parameter sweeps. I will cover the rest of the steps involved in DDEC6 charge calculation (steps 3-7) in a later post. Happy weekends!