Particle-Mesh & FFT-Based QEq

• Particle-mesh and FFT-based QEq is a computational framework that maps atomic charges onto a grid to solve the global charge equilibration problem efficiently.
• The method circumvents expensive dense Coulomb matrix inversion by using FFT-driven Poisson solvers and smooth charge density assignments.
• An iterative projected conjugate-gradient solver ensures rapid convergence and scalability for large atomistic systems, optimizing energy and force evaluations.

Particle-mesh and FFT-based charge equilibration (QEq) is a computational paradigm for solving the global charge assignment problem in atomistic simulations, particularly within fourth-generation machine learning potentials (MLPs). In these frameworks, the necessity to capture long-range charge transfer and electrostatics elevates the QEq step from a local to a global, structure-dependent optimization. The conventional QEq approach, while accurate, becomes a computational bottleneck for large systems due to the explicit calculation and inversion of a dense Coulomb matrix. Particle-mesh and FFT-based reformulations circumvent this limitation through the mapping of atomic charge densities to a regular grid and the use of spectral Poisson solvers, achieving quasi-linear scaling in both energy and force evaluations (Gubler et al., 2024).

1. Charge Equilibration: Energy Functional and Linear System

In the charge equilibration method, the total energy functional for an atomistic system is defined as

$$E(\mathbf{R},\mathbf{Q}) = E_{\text{elec}}(\mathbf{R},\mathbf{Q}) + \sum_{i=1}^N \left[ E_i + \chi_i q_i + \frac{1}{2} J_i q_i^2 \right]$$

where $\chi_i$ and $J_i$ are, respectively, the atomic electronegativity and hardness for atom $i$, and $E_i$ a reference offset. The unknowns $q_i$ are the atomic partial charges. The electrostatic contribution,

$$E_{\text{elec}} = \frac{1}{2} \sum_{i,j} q_i q_j A_{ij}$$

involves matrix elements

$$A_{ij} = \iint \rho_i(|\mathbf{r}-\mathbf{r}_i|)\,\rho_j(|\mathbf{r}'-\mathbf{r}_j|)\frac{1}{|\mathbf{r}-\mathbf{r}'|}\,d\mathbf{r}d\mathbf{r}'$$

where $\rho_i$ describes the spatial distribution of charge for atom $i$. Charge conservation $\sum_i q_i = Q_{\text{tot}}$ is imposed via a Lagrange multiplier. The stationarity of the Lagrangian yields the $(N+1)\times(N+1)$ linear system: $$\begin{pmatrix} M & \mathbf{1} \ \mathbf{1}^T & 0 \end{pmatrix} \begin{pmatrix} \mathbf{q} \ \lambda \end{pmatrix} = \begin{pmatrix} -\boldsymbol{\chi} \ Q_{\text{tot}} \end{pmatrix}$$ where $M_{ij} = A_{ij}$ for $i \ne j$ and $A_{ii} + J_i$ for $i=j$.

2. Reformulation on the Particle–Mesh Grid

Instead of explicitly assembling the dense $A_{ij}$, the formalism employs smooth, analytical charge assignments (e.g., Gaussians, B-splines) to interpolate atomic charges onto a 3D Cartesian mesh,

$$\rho(\mathbf{r}_l) = \sum_i q_i \rho_i(\mathbf{r}_l - \mathbf{r}_i)$$

thus discretizing the overall charge density. The framework then maps the problem of Coulomb interactions to mesh-based density–potential convolutions, allowing for an efficient Poisson equation treatment under periodic boundary conditions.

The electrostatic potential $V(\mathbf{r})$ is obtained by solving

$-\nabla^2 V(\mathbf{r}) = 4\pi\rho(\mathbf{r})$

discretized on the grid. The reciprocal-space (FFT) solution involves

$$\widehat{V}(\mathbf{G}) = \frac{4\pi}{|\mathbf{G}|^2} \widehat{\rho}(\mathbf{G}), \quad \mathbf{G}\ne 0$$

with the $G=0$ term determined by charge neutrality.

3. FFT-Based Solution of the Poisson Problem

The solution process consists of forward and inverse FFTs and pointwise operations in reciprocal space:

• Forward FFT: transform the mesh charge density to obtain $\widehat{\rho}(\mathbf{G})$
• Spectral Poisson Solve: compute $$\widehat{V}(\mathbf{G}) = (4\pi/|\mathbf{G}|^2) \widehat{\rho}(\mathbf{G})$$, except at $\mathbf{G}=0$
• Inverse FFT: transform back to real space to acquire $V(\mathbf{r}_l)$

The mesh-based workflow maintains efficiency by localizing all structure-dependent communication to grid operations and one global spectral solve per mat–vec.

4. Iterative QEq Solver via Projected Conjugate Gradient

The mesh-based approach facilitates the use of an iterative projected conjugate-gradient (CG) method for the reduced QEq linear system. Each mat–vec $M\mathbf{p}$ is realized by:

1. Assigning the trial vector $\mathbf{p}$ onto the mesh as charge densities
2. Computing the mesh potential via the FFT-based Poisson solver
3. Interpolating the resulting $V(\mathbf{r})$ back to the atomic centers by numerical integration
4. Adding the diagonal hardness terms

This eliminates the need to ever assemble or store $A_{ij}$. The system’s symmetry and positive-definiteness (within the charge-neutral hyperplane) ensure rapid CG convergence—empirically, $O(10{-}100)$ iterations—regardless of system size, so the cost per QEq solution is dominated by grid-FFT operations.

5. Forces and Stress: Energy Derivatives Including Implicit Charge Response

Forces and stress tensors require the derivatives of the QEq energy with respect to nuclear coordinates, explicitly incorporating the induced charge response $\partial \mathbf{q}/\partial \mathbf{R}$. For energy functionals depending on the QEq solution, the total force is

$$\frac{dE}{dR_k} = \frac{\partial E}{\partial R_k} + \sum_{i,j} \lambda_i \frac{\partial A_{ij}}{\partial R_k} q_j + \sum_i \lambda_i \frac{\partial\chi_i}{\partial R_k}$$

with $\boldsymbol{\lambda}$ obtained from the adjoint linear system $$M\boldsymbol{\lambda} = -\partial E/\partial \mathbf{Q}$$, subject to $\mathbf{1}^T \boldsymbol{\lambda} = 0$. The double sum can be evaluated efficiently on the mesh as

$$\int [ V^{Q}(\mathbf{r})\,\partial_k\rho^{\lambda}(\mathbf{r}) + V^{\lambda}(\mathbf{r})\,\partial_k\rho^{Q}(\mathbf{r}) ]\,d\mathbf{r}$$

where $V^{Q}$ and $V^{\lambda}$ are mesh potentials from $\mathbf{q}$ and $\boldsymbol{\lambda}$, respectively.

6. Algorithmic Complexity and Scaling Characteristics

Each CG mat–vec comprises two spread (assignment) and interpolate operations plus two FFTs, yielding an operation count of $$O(N_{\text{grid}}\log N_{\text{grid}}) \sim O(N\log N)$$ per iteration; the number of iterations is nearly system-size-independent. Thus, QEq solution, as well as force and stress evaluation, scales quasi-linearly with system size. Benchmarks confirm practical scaling close to $O(N \log^2 N)$, outperforming explicit Coulomb-matrix ($O(N^2{-}N^3)$) approaches for systems larger than several hundred atoms. A plausible implication is that the particle-mesh QEq formalism enables routine treatment of long-range charge fluctuations in large-scale atomistic simulations, previously prohibitive due to cost (Gubler et al., 2024).

7. Generality and Applicability of the Particle–Mesh QEq Framework

The particle-mesh and FFT-based QEq method is formulated for universal applicability, extending beyond fourth-generation MLPs to a range of classical and machine-learned force fields requiring global charge determination. The approach generalizes to any context where the system energy is a quadratic functional of smoothly distributed charges, subject to global constraints. Its main technical prerequisite is the existence of parametric, differentiable atomic density assignment schemes compatible with grid interpolation and FFT Poisson solvers. This suggests broad potential for the method in coupling quantum-accurate electrostatics with scalable molecular modeling and simulation frameworks.