Theory

PyMPFIT implements the MPFIT algorithm introduced by Ferenczy to calculate partial atomic charges that reproduce the electrostatic potential of a distributed multipole analysis. This is in contrast with traditional (e.g., RESP) partial charge methods that reproduce the molecular electrostatic potential for a given number of atomic sites. In these latter methods, site selection, molecular orientation dependence, or charge assignment of symmetrically related centers seemingly trouble such calculations.

Potential derived charges are, of course, methods in which the charges are not obtained from the wave function, but rather, indirectly via electrostatic potentials. In lieu of computing the entire wave function, constraining some lower moments of the fitted charges to molecular multipole moments serves as a better approach to utilizing information inherent to the wave function. In practice, charges based on distributed multipole potentials and electrostatic potentials are equivalent.

Gaussian Distributed Multipole Analysis

The Gaussian Distributed Multipole Analysis (GDMA) introduced by Stone provides a rigorous framework for decomposing a molecular charge density into a series of multipole moments localized on atomic centers. The electrostatic potential at a given point in space, \(\mathbf{r}\), can be expressed as a multipole expansion,

(1)\[ V(\mathbf{r}) = \frac{1}{4\pi\epsilon_0}\sum_{n=0}^{\infty}\frac{1}{r^{(n+1)}}\int (r')^n P_n(\cos{\alpha})\rho(\mathbf{r'})d\tau', \]

where \(P_n\) are Legendre polynomials and \(\rho(\mathbf{r})\) is an arbitrary, localized charge distribution. Molecular charge distributions are hence defined as

(2)\[ \rho(\mathbf{r}) = \sum_{ij} D_{ij}\chi_i(\mathbf{r})\chi_j(\mathbf{r}), \]

where \(D_{ij}\) is an element of the density matrix, and \(\chi(\mathbf{r})\) is a normalized basis function, expressed as a linear combination of Gaussian primitive functions, i.e., local multipoles

(3)\[ \chi_i(\mathbf{r}) = N_i x_i^{a_i} y_i^{b_i} z_i^{c_i} \exp{[-\zeta_i(\mathbf{r_i})^2]}, \]

where \(a_i + b_i + c_i\) are equal to the angular momentum quantum number, \(l\), and \(N_i\) is the basis function coefficient. Evaluating the higher moments of the overlap density distribution (i.e., \(\chi_i\chi_j\)) yields the multipole moment of order \(lm\),

(4)\[ Q_{lm} = \int R_{lm}(\mathbf{r})\rho(\mathbf{r})d^3\mathbf{r}, \]

where \(R_{lm}\) is a regular solid harmonic. GDMA expands Equation (4) by distributing such moments across atoms and bonds (or even arbitrary virtual sites), serving as a physically interpretable bridge between QM charge density and classical force field electrostatics.

Minimal Basis Iterative Stockholder

Minimal Basis Iterative Stockholder (MBIS) introduced by Verstraelen et al. is a variant of the Hirshfeld atoms-in-molecules method that derives atomic charges and multipoles by fitting a pro-density — an atom-centered expansion in \(s\)-type Slater density functions — to the molecular electron density. The pro-density of each atom \(A\) is expanded as

(5)\[ \rho^0_A(\mathbf{r}) = \sum_{i=1}^{m_A} \frac{N_{A,i}}{8\pi\sigma_{A,i}^3}\exp\!\left(-\frac{|\mathbf{r}-\mathbf{R}_A|}{\sigma_{A,i}}\right), \]

where the populations \(N_{A,i}\) and widths \(\sigma_{A,i}\) of each Slater shell are free parameters. These are refined iteratively to minimize the Kullback–Leibler divergence between the pro-density and a reference electron density \(\rho(\mathbf{r})\),

(6)\[ \Delta S = \sum_A \int \rho_A(\mathbf{r}) \ln\!\frac{\rho_A(\mathbf{r})}{\rho^0_A(\mathbf{r})} \, d\mathbf{r}, \]

yielding atomic populations and multipoles that scale linearly with system size and require no empirical radii or precomputed pro-atoms.

MPFIT

The MPFIT (multipole fitting) procedure originally outlined by Ferenczy represents the traditional approach to deriving atom-centered partial charges from GDMA output. Formally, this is achieved by minimizing the difference between the potential created by the distributed multipole moments, \(V^{\text{GDMA}}(\mathbf{r})\), and that generated by the fitted charges, \(V^{Q}(\mathbf{r})\),

(7)\[ f(\mathbf{r}) = V^{\text{GDMA}}(\mathbf{r}) - V^{Q}(\mathbf{r}) = \sum_{a}\sum_{l,m} Q_{lm}^{a}I_{lm}^{a}(\mathbf{r}) - \sum_{i} q_i I_{00}^{i}(\mathbf{r}), \]

where \(Q_{lm}^{a}\) denotes the \(m\)th component of the rank-\(l\) multipole moment centered at site \(a\), \(q_i\) are the point charges, and \(I_{lm}^{a}(\mathbf{r}) = r_a^{-(l+1)}C_{lm}(\theta,\phi)\) are irregular solid harmonics. The optimal charges are obtained by minimizing the integrated squared error,

(8)\[ \frac{\delta}{\delta q_j^a} \int [f(\mathbf{r})]^2 r^2 \sin\theta \, dr \, d\theta \, d\phi = 0, \]

which yields a linear system of the form \(Aq^a = b\), where the elements of \(A\) and \(b\) depend on regular and irregular solid harmonics across atomic centers and multipole ranks. Solving for \(q^a = A^{-1}b\) yields the set of least-squares charges that best reproduce the distributed multipole potential.

The explicit form of the design matrix problem requires additional mathematics. Suppose the distance between the multipole center and the point where the potential is to be calculated (\(r-r_a\)) is greater than the distance between multipoles and the charge site (\(r_{ia}\)), then

(9)\[ I_{00}^{i}(r) = \sum_{lm} R_{lm}(r_{ia}) I_{lm}^a(r) = \sum_{lm} R_{lm}^a(r_i) I^{a}_{lm}(r) \]

where \(R_{lm}(r)\) is a regular spherical harmonic defined as

(10)\[ R_{lm}(r) = r^{2l+1}I_{lm}(r) \]

This allows for simplification of the original equation to,

(11)\[ f(r) = \sum_a\sum_{l,m}I_{lm}^a[r](Q_{lm}^a - \sum_{i}q_{i}^{a}R_{lm}^{a}(r_i)) = \sum_{a}f^{a}(r) \]

We can eliminate \(I_{lm}^{a}(r)\) by integration. Namely, the integration of \([f(r)]^2\) yields the optimum set of net charges based on an appropriate integral bound in polar coordinates. Via chain rule, it can be shown that the integrand becomes:

(12)\[ \frac{\delta [f(r)]^2}{\delta q_{j}^b} = 2 \sum_a\sum_{l,m}\sum_{l',m'} I_{lm}^{a}(r) I_{l'm'}^{b}(r) R_{l'm'}(r_j) \times \left[\sum_i q_i^a R_{lm}^a(r_i) - Q_{lm}^a\right] \]

which, after plugging back into the integral, can be formulated as \(Aq = b\) where

(13)\[ A_{ij}^{ab} = \sum_{l,m}\sum_{l',m'} K_{lm,l'm'}^{ab} R_{lm}^a (r_i) R_{l'm'}^b (r_j) \]

and

(14)\[ b_j^b = \sum_a\sum_{l,m}\sum_{l',m'} K_{lm,l'm'}^{ab} R_{l'm'}^b(r_j) Q_{lm}^a \]

where

(15)\[ K_{lm,l'm'}^{ab} (\rho_1,\rho_2) = \int_{\rho_1}^{\rho_2}\int_{\theta_1}^{\theta_2}\int_{\phi_1}^{\phi_2} I_{lm}^a(r) I_{l'm'}^b(r) r^2 \sin\theta \, dr \, d\theta \, d\phi \]

Computing the off-diagonal components of this approach would be cumbersome (\(a\neq b\)), especially in the case of having some grid construction involved in molecular electrostatic potential- or field-derived methods. To this end, recall that the optimization can be broken up into the sum of \(f(r)\) at each spherical layer around point \(a\):

(16)\[ F^a(\rho_1,\rho_2) = \sum_{l,m} \frac{4\pi}{2l + 1} W_{\rho_1,\rho_2,l} \left[Q_{lm}^a - \sum_iq_i^aR_{lm}^a(r_i)\right]^2 \]

where the \(W_{\rho_1,\rho_2,l}\) factors

(17)\[ W_{\rho_1,\rho_2,l} = \int_{\rho_1}^{\rho_2} r^2r^{-2(l+1)} dr = \frac{1}{1-2l}\left(\rho_2^{1-2l} - \rho_1^{1-2l}\right) \]

weight the importance of the multipoles of rank \(l\). Now, we just need to solve for

(18)\[ \frac{\delta F^a}{\delta q_j^a} = 2 \sum_{l,m} \frac{4\pi}{2l+1} W_{\rho_1,\rho_2,l} \left[ Q_{lm}^a - \sum_i q_iR_{lm}^a(r_i)\right]R_{lm}^a(r_j) = 0 \]

where the new \(Aq^a = b\) matrix equation is solved by creating the following \(A\) and \(b\) matrices:

(19)\[ A_{ij}^a = \sum_{lm}\frac{1}{2l+1}R_{lm}^a(r_i) R_{l'm'}^b(r_j) W_{\rho_1,\rho_2,l} \]

(20)\[ b_j^a = \sum_{lm}\frac{1}{2l+1}R_{lm}^a(r_i) Q_{lm}^a W_{\rho_1,\rho_2,l} \]

where the partial charges are determined via

(21)\[ q^a = A^{-1}b \]

where

(22)\[ q_i = \sum_a q_i^a = \sum_a A_a^{-1}b^a \]

Constrained MPFIT

The unconstrained MPFIT objective \(F = \sum_a F^a\) can be augmented with two classes of constraints for transferable force field development: atom-type equivalence and per-molecule charge conservation. When fitting across multiple molecules simultaneously, these constraints enable transferable charge sets where chemically equivalent atoms in different molecules carry identical charges.

In the constrained formulation, each atom \(i\) contributes a charge \(q_i^a\) at each multipole site \(a\) for which it is within the cutoff radius (the quse mask). The total charge on atom \(i\) is

(23)\[ q_i = \sum_a q_i^a \]

Rather than optimizing all \(q_i^a\) directly, we work with a reduced parameter vector \(\mathbf{p}\) that implicitly enforces atom-type equivalence.

Atom-Type Equivalence

Atoms sharing the same type label are constrained to have equal total charges. Let \(\mathcal{T}\) be a set of atoms sharing the same type, and let \(i_1\) be the first atom in \(\mathcal{T}\) (by index order). The reference total charge for the type is defined as the sum of \(i_1\)’s per-site charges, where each \(q_{i_1}^a\) creates a free parameter in \(\mathbf{p}\) for every contributing site (i.e., where \(\text{quse}_{a,i_1} = 1\)):

(24)\[ q^{\text{total}}_{\mathcal{T}} = \sum_a q_{i_1}^a \]

For every subsequent atom \(i \in \mathcal{T}\) (\(i \neq i_1\)), each contributing site also creates a free parameter in \(\mathbf{p}\), except the last contributing site \(a^*\) (in index order), which absorbs the difference needed to match the reference total:

(25)\[ q_i^{a^*} = q^{\text{total}}_{\mathcal{T}} - \sum_{a \neq a^*} q_i^a \]

The per-site charge distributions may differ between atoms in \(\mathcal{T}\); only the totals are forced equal. If atom \(i\) contributes to only one site, no free parameters are created and \(q^{\text{total}}_{\mathcal{T}}\) is copied directly. This reduces the number of free parameters by one per subsequent atom in \(\mathcal{T}\). The mapping from \(\mathbf{p}\) to the full per-site charge matrix and total charges is performed by expandcharge.

Constrained Objective Function

The full constrained objective augments the unconstrained MPFIT objective (Equation (16)) with a per-molecule charge conservation penalty:

(26)\[ F_{\text{total}} = \underbrace{\sum_a F^a(\rho_1, \rho_2)}_{\text{multipole fit}} + \underbrace{\sum_{\text{mol}} \lambda \left(\sum_{i \in \text{mol}} q_i - Q_{\text{mol}}\right)^2}_{\text{charge conservation}} \]

where \(F^a(\rho_1, \rho_2)\) is defined in Equation (16). The parameter \(\lambda\) (conchg) controls the strength of the charge conservation penalty. Larger values of \(\lambda\) enforce stricter conservation at the cost of a slightly worse multipole fit. \(Q_{\text{mol}}\) is the target total charge for each molecule (e.g., +1 for a cation, 0 for a neutral).

The gradient of \(F_{\text{total}}\) with respect to the full per-site charges has two contributions. The multipole fitting gradient \(\partial F^a / \partial q_j^a\) follows directly from Equation (18). The charge conservation gradient for atom \(j\) in molecule \(m\) is

(27)\[ \frac{\partial F_{\text{con}}}{\partial q_j} = 2\lambda\left(\sum_{i \in m} q_i - Q_m\right) \]

The full-space gradient is then projected into the reduced parameter space via the chain rule through expandcharge, yielding the gradient \(\nabla_{\mathbf{p}} F_{\text{total}}\) passed to the optimizer.

Optimization

The constrained objective is minimized using scipy.optimize.minimize (L-BFGS-B by default) with the analytic Jacobian. The optimizer receives the reduced parameter vector \(\mathbf{p}\), which is expanded to full charges for each function and gradient evaluation. After convergence, the optimal \(\mathbf{p}\) is expanded one final time to obtain the total charge per atom \(q_i\).

Virtual Sites

Coming soon…