Analytic Model Charge Density for Ewald Summation

• Model Charge Density is an auxiliary analytic distribution that matches lower multipole moments of the true charge density to enhance convergence in periodic electrostatics.
• It employs a multipole expansion and subtraction scheme whose cancellation of moments up to order L improves error scaling and reduces computational costs.
• The approach is universally applicable in quantum and classical simulations, remaining independent of basis sets and robust for diverse crystalline systems.

A model charge density is an auxiliary distribution, constructed from the true charge density by exact multipole matching, and serves as an analytic device for accelerating convergence in Ewald summations of electrostatics for periodic systems. By subtracting a model density that reproduces the lower multipole moments of the system and is analytically tractable, one ensures that the rapidly converging lattice sums of the residual (difference) density control all higher multipolar contributions, thus facilitating highly efficient electrostatic potential calculations in both quantum and classical condensed-phase simulations (Ribaldone et al., 29 Jan 2026).

1. Mathematical Formulation of the Model Charge Density

Let $n(\mathbf r)$ denote the true crystalline charge density, incorporating both electrons and nuclei. One defines a model charge density, $\rho_{\mathrm{model}}(\mathbf r)$, which is expanded in a finite multipole basis up to order $L$ around a center $\mathbf R$: $$\rho_{\mathrm{model}}(\mathbf r) = \sum_{\ell=0}^L \sum_{m=-\ell}^\ell \eta_\ell^m[n](\mathbf R) \, \kappa_\ell^m(\mathbf r - \mathbf R)$$ where the multipole moments are

$$\eta_\ell^m[n](\mathbf R) = \int_{\mathbb R^3} n(\mathbf r) X_\ell^m(\mathbf r - \mathbf R) \, d^3 r$$

with $X_\ell^m$ the unnormalized solid spherical harmonics, and $\kappa_\ell^m$ appropriate model functions determined by orthonormality (Ribaldone et al., 29 Jan 2026).

The difference density is then

$$\Delta n(\mathbf r) = n(\mathbf r) - \rho_{\mathrm{model}}(\mathbf r)$$

which is constructed to satisfy

$$\int \Delta n(\mathbf r) X_\ell^m(\mathbf r - \mathbf R) d^3 r = 0 \quad \forall~\ell\leq L,~m.$$

This ensures $\Delta n$ is free of net charge, dipole, and all higher multipoles up to order $L$.

2. Integral Constraints and Explicit Construction

The model imposes the integral matching for all moments up to $L$: $$\int r^\ell Y_{\ell m}(\hat r)\, \rho_{\mathrm{model}}(\mathbf r)\, d^3 r = \int r^\ell Y_{\ell m}(\hat r)\, n(\mathbf r)\, d^3 r, \quad (\ell \leq L)$$ This is realized by setting

$$\kappa_\ell^m(\mathbf r - \mathbf R) = \mathcal R_\ell(r) X_\ell^m(\mathbf r - \mathbf R), \quad r = |\mathbf r - \mathbf R|$$

with radial part

$$\mathcal R_\ell(r) = C_\ell^m\, \delta(r)\, r^{-2(\ell+1)},$$

and normalization constant $C_\ell^m$ chosen such that each component is a point multipole of order $\ell$ centered at $\mathbf R$. The final form is: $$\rho_{\mathrm{model}}(\mathbf r) = \sum_{\ell=0}^L \sum_{m=-\ell}^\ell \eta_\ell^m[n](\mathbf R)\, C_\ell^m\, \delta(|\mathbf r - \mathbf R|)\, |\mathbf r - \mathbf R|^{-2(\ell+1)}\, X_\ell^m(\mathbf r - \mathbf R)$$ yielding a model density formed from analytic point multipoles (Ribaldone et al., 29 Jan 2026).

3. Modified Ewald Summation Scheme

The standard Ewald summation is applied to the difference density $\Delta n$, greatly accelerating convergence since all low-order multipoles vanish by construction. The electrostatic potential is decomposed into three terms: $$\Phi[n](\mathbf r) = \sum_{\mathbf g} \int_\Omega \Delta n(\mathbf r') \, \frac{\operatorname{erfc}(\sqrt\kappa |\mathbf r - \mathbf r' - \mathbf g|)}{|\mathbf r - \mathbf r' - \mathbf g|} d^3r' + \frac{4\pi}{V} \sum_{\mathbf G \neq 0} \frac{e^{-|\mathbf G|^2 / (4\kappa)}}{|\mathbf G|^2} \Delta n(\mathbf G) e^{i \mathbf G \cdot \mathbf r} - \frac{2\pi}{3V} \int_\Omega \Delta n(\mathbf r') |\mathbf r'|^2 d^3 r'$$ where reciprocal- and real-space sums are now absolutely convergent for any system (including those with nonzero total charge or dipole), provided the model density is constructed accordingly (Ribaldone et al., 29 Jan 2026).

4. Error Scaling and Convergence Acceleration

If only monopole cancellation ($\ell = 0$) is performed, the real-space tail falls as $1/R_c$ and the reciprocal-space tail as $1/G_c$. Cancelling all multipoles up to order $L$ causes the first non-vanishing moment to be at order $L+1$, and the convergence scales as $1/R_c^{L+1}$ and $1/G_c^{L+1}$. Each increment in $L$ improves convergence by a full power of the cutoff, so typically, quadrupole-level cancellation ($L=2$) reduces computational cost by more than an order of magnitude (Ribaldone et al., 29 Jan 2026).

5. Implementation and Computational Details

In practical implementations (e.g., CRYSTAL), at each self-consistent-field iteration, one evaluates the model multipoles via

$$\eta_\ell^m = \sum_{\mu\nu} P_{\mu\nu} \int \phi_\mu(\mathbf r) \phi_\nu(\mathbf r) X_\ell^m(\mathbf r - \mathbf R) d^3 r$$

with atom-centered Gaussian basis functions and density matrix $P_{\mu\nu}$. The analytic model multipole corrections (notably for $\rho_{\mathrm{model}}$) are computed using internal analytic forms, and are combined with the numerically evaluated difference density. All terms are processed using a single screening parameter $\kappa$ without accuracy loss, leveraging the absolute convergence of all terms (Ribaldone et al., 29 Jan 2026).

6. Independence from Basis Set and Physical Interpretation

The model charge density construction is entirely basis-set independent: it applies to any quantum or classical description where a charge density $n(\mathbf r)$ is available. Physically, the model charge density acts as an auxiliary device, serving only to capture and analytically subtract the problematic (slowly converging) multipole contributions from the true density, leaving a rapidly vanishing residual. The total electrostatic potential and energies remain exact—no approximation other than the imposed multipole truncation is introduced (Ribaldone et al., 29 Jan 2026).

7. Applications and Impact

The model charge density method is crucial for accurate electrostatics in bulk simulations, enabling rapid and robust convergence for arbitrary unit cells (triclinic, monoclinic, etc.), arbitrary charge distributions (delocalized electrons, polar crystals, metallic/nonmetallic), and any basis (plane-wave, Gaussian, real-space grid). Its adoption clarifies and generalizes earlier implementations (notably in CRYSTAL), extends readily to higher-order multipolar systems, and applies universally across density-functional, Hartree–Fock, and classical force-field contexts (Ribaldone et al., 29 Jan 2026).