Permanent electrostatic moments through the lens of atoms: assessing variational Hirshfeld methods

Abstract

Atomic multipoles are key components of many computational models; however, there is no unique method for computing them. This motivates our benchmark study to evaluate how accurately the atomic moments from variational methods approximate the molecular moments, offering insight into their ability to capture the underlying electron density. We show the chemical utility of higher-order atomic multipoles and demonstrate their conformational stability across diverse organic and inorganic molecules and protein fragments. We focus on the recently proposed Additive Variational Hirshfeld (AVH) method and compare its performance to two other popular variational Hirshfeld approaches, the Iterative Hirshfeld (HI) and Minimal Basis Iterative Stockholder (MBIS) methods, as well as the Charge Model 5 (CM5) and electrostatic potential fitted charges. We first show that commonly used integration grids can introduce significant numerical errors, undermining the reliability of quantitative comparisons of atomic dipoles and quadrupoles. We then demonstrate that while AVH charges may not always provide the most accurate approximation of molecular moments on their own, they outperform HI and MBIS in approximating the molecular quadrupole when atomic dipoles are taken into account. More severely, for HI and MBIS, the addition of atomic dipoles can sometimes increase the molecular quadrupole error. These findings suggest that AVH moments offer a more systematic improvement with increasing multipole order.

---

1 Introduction

Molecular dipole, quadrupole, and higher-order moments are important molecular characteristics that determine how a molecule interacts with external electric fields and with other molecules.¹ These molecular moments arise from the electronic charge distribution and influence physical and chemical properties of molecules, like spectroscopic transition intensities and intermolecular forces. Accurate computation of these moments is essential for reliable quantum chemistry models, which is why modern density functional development often uses dipole moments as a key benchmark for accuracy.^2–5 For a molecule with M atomic nuclei {Z_A} positioned at {R_A} and the electron density ρ_mol(r), the molecular moments are calculated by (in atomic units):where q, µ, and [Q]_ij are the molecular charge (L = 0), molecular dipole vector (L = 1), and the ij-th element of traceless molecular quadrupole matrix (L = 2), respectively; here, L represents the order of the molecular multipole moment. In this work, we will not consider molecular multipoles beyond second-order.

Molecular moments are often expressed in terms of their atomic counterparts, known as the atomic multipole moments, which are computed through postprocessing quantum chemistry calculations. This decomposition is chemically intuitive and provides an atom-centered description of the electronic charge distribution within a molecule. In addition, atomic multipole moments are a key component in computational chemistry research, such as force field development,^6–14 electrostatic potential modeling,^15–19 reactivity analysis,^20–32 and machine learning approaches for predicting long-range interactions.^33–39 Given a set of M atomic densities each centered at the position of the corresponding nucleus, {ρ_A(r)}, which collectively reconstruct the molecular electron density, the molecular moments can be expressed in terms of its atomic-centered atomic moments by replacing the molecular density with its sum of atomic density (in atomic units):

where N_A, q_A, µ_A, and [Q_A]_ij denote the atomic population, atomic charge (l = 0), atomic dipole vector (l = 1), and the ij-th element of the traceless atomic quadrupole matrix (l = 2), respectively; here, l represents the order of atomic multipole moment. This shows the contribution of individual atoms to the molecular moments and demonstrates that atomic multipoles of rank l > L do not contribute to molecular multipoles of rank L.

While the mathematical relationships between molecular and atomic moments are well-defined, computing these atomic terms remains a challenge because there is no universally accepted definition for atomic densities {ρ_A(r)} and atomic moments. This has motivated the development of various density-based atoms-in-molecule schemes that exhaustively partition the molecular electron density into atomic contributions.^40–43 Alternatively, one can bypass the need for explicitly defining atomic densities by using approaches that fit atomic moments to reproduce either molecular moments or the electrostatic potential, ensuring consistency with the overall electronic distribution.^44–46 The caveat of these methods is that because there is no atomic density, other atomic properties (like volume) cannot be computed. The existing methods result in atomic moments that can vary both qualitatively and quantitatively, so assessing their relative performance is critical.^14,47 While atomic charges are commonly used to compare different methods, higher-order atomic moments remain less explored.^48–54 Molecular multipole moments are used as a stringent test for assessing atomic moments.^49,55–58 Atomic point charges cannot capture the anisotropic nature of the electron density,^15,16 so even though some may reproduce the total dipole moment reasonably well, they often fail to accurately represent the quadrupole moment.

This study investigates how the recently proposed additive variational Hirshfeld (AVH)⁵⁹ atomic moments approximate molecular dipole and quadrupole moments across a range of organic, inorganic, and protein-fragment systems. We demonstrate the importance of accurate integration grids for obtaining reliable numerical results when computing atomic moments. We compare our results with mathematically related Hirshfeld methods, the Iterative Hirshfeld (HI)^60,61 and Minimal Basis Iterative Stockholder (MBIS)⁶² methods, as well as Charge Model 5 ⁶³  and electrostatic potential derived charges, which are designed to accurately reproduce the molecular dipole moment and electrostatic potential, respectively. This assessment provides a deeper insight into density-derived atomic moments, guiding their development and application in computational modeling.

2 Preliminaries

2.1 Hirshfeld atoms-in-molecules

The Hirshfeld method defines atomic densities, {ρ_A(r)}, that closely resemble reference atom-centered densities, {ρ_A⁰(r)}, known as proatoms.^43,64–69 For atom A, the corresponding Hirshfeld atomic density is defined as,where ρ_mol⁰(r) denotes the promolecular density defined as the sum of atom-centered proatom densities. While this atomic density definition was empirically proposed by Hirshfeld in 1977,⁷⁰ it was later demonstrated that the family of f-divergences is both necessary and sufficient for deriving Hirshfeld atomic densities.^{64–69,71–74} This rigorously established the mathematical framework of Hirshfeld partitioning schemes, highlighting its pervasiveness and formalizing the development of variational methods for selecting proatoms.^43,59,69 Hirshfeld originally selected isolated neutral ground-state atomic densities as proatoms,⁷⁰ but this arbitrary choice naturally results in small atomic charges and is unsuitable for charged molecules and ionic compounds (e.g., LiH).⁷⁵ To address this, numerous Hirshfeld variants introduce system-specific proatoms, leading to improved atomic densities and properties.^{13,56,59,60,62,76–88} In parallel, reference-free approaches have been proposed that bypass proatoms altogether.⁸⁹ Below, we briefly describe the methods used in this study.

The popular iterative Hirshfeld (HI)⁶⁰ partitioning scheme iteratively optimizes the proatom charges to match the corresponding atomic charges, where the proatom density with fractional population is obtained through piecewise linear interpolation. It was recently shown that this self-consistent approach can be reformulated as minimizing the deviation between atomic and proatom populations using an invex function.⁶¹ Although HI is a widely used approach,^{13,51,61,90–95} it often relies on highly charged anions that are diffuse or unbound, making it computationally undesirable. This leads to atomic densities with slow decay and atomic charges that are often larger than expected based on chemical intuition.^43,59,84,96

2.2 Variational Hirshfeld methods

Leveraging the f-divergence framework,⁶⁹ two variational Hirshfeld methods were proposed, minimal basis iterative stockholder (MBIS)⁶² and additive variational Hirshfeld (AVH).⁵⁹ Although they use different proatom models, both methods optimize their proatom parameters to minimize the extended Kullback–Leibler divergence between molecular and promolecular densities. This guarantees that each proatom model produces a promolecular density that closely approximates the molecular density. MBIS represents each proatom as a sum of core and valence densities, each modeled using an s-type Slater functions, where their width and population are optimized. Due to the presence of nonlinear parameters (i.e., width), the optimization does not have a unique solution, and therefore, MBIS requires reasonable initial guesses.⁶² Since its introduction in 2016, MBIS has been rapidly adopted for a wide range of applications, like deriving non-bonded force field parameters,^8–12 calculating hydration free energies,^97,98 improving host–guest binding affinity predictions,¹⁰ and providing reference atomic moments for training machine learning models.^33,35,37,38

On the other hand, AVH utilizes a linear combination of isolated neutral and charged reference densities, where their coefficients are optimized; this results in a convex optimization with a unique solution. The AVH proatom densities are written as:

where K_A is the number of proatom basis functions for atom A, and c_A = [c_A,1,c_A,2,…,c_A,KA] are the positive basis coefficients for each basis function ρ⁰_A,k(r). Setting K = 1 for all atoms results in the minimal AVH (AVH-M), where a single coefficient scales the neutral reference densities. AVH-M closely aligns with the original Hirshfeld method in spirit. Including all the bound charged states of each atom defines the bound AVH (AVH-B). In this approach, cationic and neutral proatomic basis functions are included for all atoms, while anionic basis functions are only incorporated if they are physically bound (e.g., anions are excluded for nitrogen atoms).^59,96

2.3 Datasets

In this study, we considered datasets that include organic and inorganic molecules as well as protein fragments. The organic dataset consists of 214 molecules containing H, C, N, O, F, Si, P, S, Cl, and Br atoms, with a diverse range of hybridizations and molecular topologies. This dataset combines two previously published datasets: one developed to evaluate the performance of iterative Hirshfeld charges in approximating the molecular electrostatic potential,¹⁷ and the other designed to benchmark quantum chemistry methods for computing molecular dipole moments and polarizabilities against experimental data.³ The Silica dataset consists of 248 silica clusters that were curated to evaluate the transferability of atomic charges derived from the Hirshfeld-E partitioning scheme and their accuracy in reproducing the molecular electrostatic potential.⁹⁹ These clusters follow the Si_xO_2x−yH_2y chemical formula, including all possible silicon-oxide structures with two to eight silicon atoms that adhere to specific geometric and bonding criteria to realistically model silica materials.

The alanine–dipeptide dataset includes 10 conformations of alanine dipeptide sampled to evaluate the conformational stability of atomic charges computed by the constrained iterative Hirshfeld method.⁶¹ These conformers cover a range of molecular geometries relevant to biological systems, sampling ψ dihedral angle from −30°to 150° with steps of 20°, each with variationally optimized ϕ dihedral angles. The protein fragments include backbone and sidechain monomers derived from 200 backbone–backbone and 610 sidechain–sidechain dimers sampled from molecular dynamics trajectories to investigate the nature of pi-contact interactions.¹⁰⁰ Compiling the neutral monomers, we obtain the Backbone dataset, consisting of 400 N-methylacetamide conformers, and the Sidechain dataset containing 790 neutral amino acid sidechains, including aromatic (130 PHE, 180 TRP, 70 TYR, and 60 HIS) and hydrophilic (170 ASN and 180 GLN) conformers. These two datasets provide conformational sampling of N-methylacetamide and different sidechains.

Fig. S1 shows the distribution of molecular dipole norms, highlighting the range and variability of the dipole magnitudes for molecules in each dataset.

3 Results and discussion

To qualitatively and quantitatively compare atomic moments, it is essential to ensure they are numerically reliable. Section 3.1 shows that while charges can be computed accurately using relatively small molecular integration grids, dipole and quadrupole moments require substantially finer grids to achieve numerical stability. We then assess the performance of numerically reliable atomic moments obtained using additive variational Hirshfeld (AVH) and compare them with other Hirshfeld-based methods, including Hirshfeld (H), iterative Hirshfeld (HI), and minimal basis iterative stockholder (MBIS). Section 3.2 and 3.3 investigate how well lower-order atomic moments approximate molecular dipole and quadrupole moments, respectively, and assess the computational robustness of these approximations. When using only atomic charges, we also compare our results to Charge Model 5 (CM5)⁶³ and Hu, Lu, and Yang (HLYGAt)⁴⁶ methods. CM5 refines Hirshfeld charges to better reproduce molecular dipole moments, making it a suitable benchmark for evaluating our variational approaches. HLYGAt provides electrostatic-fitted charges that best approximate the electrostatic potential on molecular surfaces. Section 3.4 further examines the conformational stability of atomic moments in reproducing molecular multipole moments.

3.1 Sensitivity to numerical integration grids

Although molecular moments can be computed analytically, calculating atomic moments requires numerical integration, so their accuracy depends on the molecular grid used.^102,103 Using a Becke–Lebedev¹⁰⁴ molecular grid, Fig. 1 shows how changing the number of radial and angular grid points affects the accuracy of numerical integration when calculating molecular moments using eqn (1)–(3). It shows that molecular charges can be accurately reproduced with relatively small integration grids; however, dipole and quadrupole moments converge more slowly with increasing grid size. In general, the number of radial and angular grid points should be increased together to achieve higher accuracy. For example, increasing the angular grid alone does not necessarily reduce the error to zero if the radial grid is too sparse; this is demonstrated when using only 50 radial points in Fig. 1. Also, for a given number of radial grid points, different numbers of angular grid points are required to achieve the desired accuracy for each moment. For example, using 100 radial grid points for the molecule in Fig. 1, at least 170, 434, and 974 Lebedev angular grid points are required to reach a 10⁻² accuracy threshold for charge (in e), dipole components (in D), and quadrupole components (in D Å), respectively. Fig. S2–S5 demonstrate that these observations hold for different molecules sampled from each dataset. Using a different radial grid type (e.g., Gauss–Chebyshev¹⁰³) can lead to faster convergence with respect to the number of radial and angular grid points, especially at lower radial resolutions, but still, similar overall trends were observed.

Fig. 1 Convergence of numerical integration error when computing molecular multipole moments, against the exact moments from the ωB97X-D/Def2-TZVPD level of theory, for molecular charge (left plot), y-components of molecular dipole (middle plot), and yy-component of molecular quadrupole (left plot) of the largest molecule in the Organic dataset. Molecular integrations were performed with the exponential radial grid and the Becke–Lebedev angular grid of HORTON 2,¹⁰¹ using different numbers of radial (nrad label) and angular (x-axis) grid points. Here, dipole and quadrupole components with the highest errors are presented, and other components of the molecular moments are included in Fig. S2. Dotted vertical green lines indicate the minimum number of Lebedev angular grid points required (170 for charge, 434 for y-component of dipole, and 974 for yy-component of quadrupole) to achieve an error below 10⁻² for each multipole moment when using 100 radial grid points.

This shows the importance of using accurate numerical grids when computing atomic dipole and quadrupole moments to minimize numerical error, particularly when quantitatively comparing different schemes. For the calculations performed in the following sections, we ensured that the molecular grids achieved a numerical integration accuracy of 10⁻³ a.u. for molecular dipole components and 10⁻⁴ a.u. for molecular charges; see Section 4 for further details.

3.2 Approximating molecular dipole

Eqn (5) decomposes the molecular dipole into charge-transfer contributions, arising from atomic charges, and polarization contributions, arising from atomic dipole. This means that the molecular dipole (L = 1) can be approximated using only atomic charges (l = 0) by:To assess the importance of polarization contributions (l = 1), Table 1 reports the mean dipole error for each dataset, capturing deviations in both magnitude and direction from the exact dipole. As expected, the approximated molecular dipoles using HLYGAt charges show an almost perfect agreement with the exact molecular dipole, reaching a mean dipole error <0.05 D for all datasets. On the other hand, Hirshfeld (H) charges exhibit the poorest performance, with the highest mean dipole errors (except for the Silica dataset), underscoring the importance of polarization contributions from the H scheme. CM5 was developed to improve H charges for reproducing molecular dipoles, so it has lower mean dipole errors for all datasets except for Silica, where H charges perform better. Across these datasets, MBIS charges perform comparably with CM5 charges in reproducing molecular dipoles, although they were not explicitly parameterized for this purpose. After HLYGAt, HI charges result in the lowest mean dipole errors for the Backbone and Sidechain datasets, but yield higher errors for the Organic and Silica datasets compared to CM5 and MBIS. AVH-M introduces minimal flexibility to the Hirshfeld proatom model, which lowers the mean dipole error; however, its performance still lags behind CM5, HI, and MBIS. Using a more flexible proatom model, AVH-B charges offer an improved description of molecular dipoles. AVH-B charges perform almost perfectly on the Silica dataset (comparable to ESP-fitted charges), are on par with MBIS and CM5 for the Backbone dataset, and match CM5's performance on the Sidechain dataset. However, for the Organic dataset, AVH's performance is close to HI, but falls below that of CM5 and MBIS.

Table 1 Mean dipole error (in Debye) for approximating the molecular dipole using point charges, as defined in eqn (9), compared to the exact molecular dipole computed at the ωB97X-D/Def2-TZVPD level of theory. The error is defined as the average norm of the difference vector |µ_exact − µ_approx| over all molecules in each dataset. Fig. S6 presents the error distributions across each dataset

Schemes	Organic	Silica	Backbone	Sidechain
H	0.679	0.053	1.113	0.771
AVH-M	0.469	0.095	0.488	0.459
AVH-B	0.413	0.099	0.377	0.389
HI	0.366	1.283	0.280	0.242
MBIS	0.269	0.583	0.362	0.238
CM5	0.202	0.460	0.363	0.341
HLYGAt	0.03	0.02	0.048	0.029

---

Taking a closer look, Fig. 2 shows how accurately different charge models reproduce the molecular dipole moment norm, revealing overall trends for different datasets. In addition, Fig. 3 compares the magnitude and direction of µ_approx to µ_exact for each molecule in various datasets. For all datasets, the magnitude and direction of µ_approx by HLYGAt are in excellent agreement with the reference exact value across datasets. This is because ESP-derived charges implicitly capture the effect of higher multipoles. Although the approximated dipole norms from H charges show over 90% correlation with the exact values, they almost always underestimate the dipole magnitude (except for the Silica dataset and the PHE amino acid conformers, which have the smallest dipoles). At the same time, they generally capture the correct direction. These observations are not surprising, given the characteristically small charges of the Hirshfeld method. Although CM5 charges still tend to underestimate dipole magnitudes, their maximum errors are notably smaller. CM5's maximum |µ_exact − µ_approx| errors are 0.83, 1.18, 0.52, and 0.66 Debye for the Organic, Silica, Backbone, and Sidechain datasets, respectively. They also outperform Hirshfeld charges in capturing the dipole direction, except for the Silica dataset and one other exceptions in Fig. 3. The largest CM5 errors are observed for the Silica dataset, which is unsurprising given that these systems were underrepresented in the dataset used to parameterize CM5.

Fig. 2 Comparing the point-charge approximation of molecular dipole given in eqn (9) to the exact value computed at ωB97X-D/Def2-TZVPD level of theory for (a) organic, (b) silica, (c) backbone, and (d) sidechain datasets using Hirshfeld (H), Iterative Hirshfeld (HI), minimal basis iterative stockholder (MBIS), Additive Variational Hirshfeld with Minimal (AVH-M) and Bound (AVH-B) proatom basis, Charge Model 5 (CM5), and Hu, Lu, and Yang (HLYGAt) electrostatic-fitted charges. The dipole norms are reported in Debye (D), and the R-squared is reported for each partitioning method. The dashed gray line is the y = x representing the identity line. Fig. S6 extends this analysis by presenting the error distributions for each dataset.

Fig. 3 Cosine similarity versus norm difference for point-charge-approximated and exact molecular dipoles, computed at the ωB97XD/Def2-TZVPD level of theory for (a) organic, (b) silica, (c) backbone, and (d) sidechain datasets. The cosine similarity of two vectors, defined as measures the alignment of two vectors. A cosine similarity of 1.0 indicates perfect alignment of vector directions, which–together with a norm difference of zero–represents the ideal agreement between the approximate and exact dipole moments. The dashed gray line at x = 0 highlights the boundary between overestimation and underestimation of the molecular dipole, making these deviations more visually apparent.

MBIS charges tend to overestimate dipole magnitudes due to their generally larger point charge values. MBIS's maximum |µ_exact − µ_approx| errors are 1.34, 1.38, 0.57, and 0.55 Debye for the Organic, Silica, Backbone, and Sidechain datasets, respectively. However, MBIS charges are not as good as CM5 in capturing the direction of molecular dipole moments. HI charges also tend to overestimate dipole magnitudes more often than MBIS charges, but also underestimate them sometimes. For HI charges, the maximum |µ_exact − µ_approx| errors are 2.08, 3.03, 0.50, and 0.69 Debye for the Organic, Silica, Backbone, and Sidechain datasets, respectively. Additionally, HI charges more frequently fail to capture the correct direction of the molecular dipole, as illustrated in Fig. 3 for Silica and Sidechain datasets. This suggests that HI charges may better capture dipole behavior in biomolecular fragments but are less reliable for broader chemical diversity or more polar inorganic systems. Overall, HI charges have the worst performance for the Silica dataset. Like H charges, AVH-M tends to underestimate the dipole norm, though to a much lesser extent (especially evident in the Backbone and Sidechain datasets), and offers a slight improvement in accurately capturing dipole direction. AVH-B charges perform almost perfectly on the Silica dataset (comparable to ESP-fitted charges), are on par with MBIS and CM5 for the Backbone dataset, and match CM5's performance on the Sidechain dataset. However, for the Organic dataset, AVH's performance is close to HI, but falls below that of CM5 and MBIS. For AVH-B charges, the maximum |µ_exact − µ_approx| errors are 1.55, 0.25, 0.59, and 1.00 Debye for the Organic, Silica, Backbone, and Sidechain datasets, respectively.

It has been shown that the dipole moment depends on the level of theory used.^2–4 Fig. S7 evaluates the robustness of molecular dipole norms and mean dipole errors for the Organic dataset, based on nine different levels of theory. This demonstrates that the error in the point-charge approximation of the molecular dipole is insensitive to the level of theory for all partitioning schemes considered. This is consistent with earlier reports for the Organic dataset, indicating that the atomic charges remain unchanged when the method or basis set is changed.^43,58,59

3.3 Approximating molecular quadrupole

Based on eqn (6), the molecular quadrupole can be approximated using atomic charges and dipoles through Q⁽⁰⁾_approx + Q⁽¹⁾_approx, where:Adding the atomic quadrupole term, Q⁽²⁾_approx, results in the total molecular quadrupole (up to the numerical integration error) as defined in eqn (3). For each molecule, the quadrupole error is given by its Frobenius norm:where Q_exact represents the exact molecular quadrupole, computed analytically at a given level of theory. Table 2 presents the mean quadrupole error for each dataset using different partitioning schemes. The results show that the AVH methods, particularly AVH-B, consistently outperform both HI and MBIS across all datasets. This is most notable for the Backbone and Sidechain datasets, where the AVH-B average error is < 0.7 Debye Å. HI results in the highest quadrupole errors, even larger than H, which indicates its exaggeration of atomic moments. Surprisingly, MBIS achieves a comparable error to AVH methods for the Silica and Organic datasets, despite having higher errors in other cases. Fig. 4 presents the corresponding quadrupole error distributions for each scheme and dataset. This more clearly demonstrates that the error distribution from AVH-B is consistently several orders of magnitude smaller for all molecules across datasets. This indicates the superior performance of AVH-B when atomic dipoles are included, and suggests that atomic quadrupole moments from the HI and MBIS methods are often significant and should not be neglected.

Table 2 Mean quadrupole error (in Debye⋅Å) for approximating the molecular quadrupole using atomic charges and atomic dipoles, Q⁽⁰⁾_approx + Q⁽¹⁾_approx defined in eqn (10), computed at the ωB97X-D/Def2-TZVPD level of theory. The error is calculated using eqn (11) for each molecule, and then averaged over all molecules in each dataset. CM5 and HLYGAt are purely charge-based models, so they cannot be used for this approximation

Schemes	Organic	Silica	Backbone	Sidechain
H	0.983	0.573	0.774	0.934
AVH-M	1.001	0.468	0.795	1.044
AVH-B	0.667	0.326	0.236	0.382
HI	1.407	3.040	1.296	2.032
MBIS	0.879	0.398	2.100	2.992

---

Fig. 4 Distribution of errors in approximating the molecular quadrupole moment (in Debye Å) using atomic charges and dipoles, Q⁽⁰⁾_approx + Q⁽¹⁾_approx, defined in eqn (10), computed at the ωB97X-D/Def2-TZVPD level of theory for (a) organic, (b) silica, (c) backbone, and (d) sidechain datasets. For each molecule, the error is calculated as the Frobenius norm of the full difference matrix |Q_exact − Q_approx| defined in eqn (11). The mean absolute error (MAE) values provided in the legend match Table 2.

In addition, Table S1 presents the quadrupole approximation errors when only using atomic charges, denoted by Q⁽⁰⁾_approx, for each dataset. In this case, HI and MBIS have the lowest average errors, except for the Silica dataset, where AVH-M and AVH-B achieve errors that are several orders of magnitude smaller. This is consistent with the results in Section 3.2, where atomic charges were used to approximate molecular dipoles, reinforcing that the generally larger magnitudes of HI and MBIS charges lead to more accurate point-charge approximations of higher-order moments. Similarly, ESP-fitted charges exhibit the lowest errors, while CM5 charges perform comparably to MBIS charges. Hirshfeld charges consistently exhibit higher average errors across most datasets, except for the Silica dataset, where they outperform HI and MBIS and are comparable to AVH. Although including atomic dipoles reduces the overall error distribution and lowers its mean, this trend does not consistently hold for individual molecules. For some molecules, the Q⁽⁰⁾_approx + Q⁽¹⁾_approx approximation can produce larger errors than Q⁽⁰⁾_approx alone. For each scheme, Table S2 reports the percentage of molecules in each dataset for which the error increases when the atomic dipole is included. These convergence anomalies are more common for HI and MBIS, while they are significantly less frequent for AVH methods.

Fig. S9 shows that molecular quadrupole moments themselves are to some extent sensitive to the level of theory, especially when small basis sets (e.g., def2-SVPD) or the Hartree–Fock method are used. The mean quadrupole error for Q⁽⁰⁾_approx remains relatively insensitive to the level of theory used. This is also reflected in the computational robustness of atomic charges in Fig. S10. However, the mean quadrupole error for Q⁽⁰⁾_approx + Q⁽¹⁾_approx shows greater variability with respect to the level of theory, especially when small basis sets are used. This variability is most pronounced for the HI scheme, and to a lesser extent for H, AVH-M, and AVH-B. A likely reason for this sensitivity is the use of negatively-charged proatoms in these partitioning schemes, which are computed at the same level of theory as the molecular calculation. The increased sensitivity observed for HI is consistent with the behavior of its atomic dipoles and quadrupoles, which are more affected by the level of theory than atomic charges, as shown in Fig. S10. Specifically, the asymptotic decay of negatively charged proatomic densities, commonly used in HI, is known to exhibit greater variability with the level of theory.⁹⁶ This effect is less pronounced in the other methods, as they either do not rely on negatively charged proatoms (e.g., H and AVH-M) or use only physically bound anions as proatoms (e.g., AVH-B). MBIS shows minimal dependence on the level of theory because it does not rely on proatoms. These observations support the idea of using fixed proatom densities for schemes that depend on proatoms, to improve their computational robustness.

3.4 Conformational stability

The molecular dipole norm for alanine dipeptide conformers varies significantly with changes in the dihedral angle, ranging from over 6 Debye to less than 2 Debye (Fig. S1). Fig. 5 shows that the dipole and quadrupole approximation errors introduced in Sections 3.2 and 3.3 vary smoothly with molecular conformation. This is expected as atomic charges, atomic dipole norms, and atomic quadrupole norms remain largely stable across the 10 conformers for most partitioning schemes (Fig. S10–S12). For all schemes, except HLYGAt, it is notable that dipole approximation errors are larger for conformers with larger molecular dipoles (see Fig. S1). While HLYGAt charges exhibit a pronounced sensitivity to conformational changes (Fig. S10), they yield the lowest and most conformationally stable errors when approximating the molecular dipole (see left plot in Fig. 5). For other schemes, the dipole errors are more sensitive to conformational change, with H, AVH-M, and MBIS in particular showing noticeable variation in dipole approximation error across the dihedral scan. AVH-B and CM5 charges exhibit very similar profiles, surpassing the accuracy of HI beyond ψ = 70°.

Fig. 5 Conformational sensitivity of dipole and quadrupole moment approximations for Alanine–Dipeptide conformers using various atomic partitioning schemes: (Left) Norm of the difference between the molecular dipole and that approximated using atomic charges (|µ_exact − µ_approx|), and (Right) Frobenius norm of the difference between the molecular quadrupole and the approximation including both monopoles and dipoles (|Q_exact − (Q⁽⁰⁾_approx + Q⁽¹⁾_approx)|). Errors are shown as a function of the ψ dihedral angle across ten conformers.

In approximating molecular quadrupoles using atomic charges and dipoles (see right plot in Fig. 5), the HI scheme exhibits the largest errors, while AVH-B produces errors that are smaller by several orders of magnitude. MBIS errors are slightly higher than AVH-B. These are consistent with results in Table 2. AVH-B shows robust performance, achieving the lowest quadrupole errors overall when both charges and dipoles are considered, and maintaining stability across conformational space.

4 Computational Details

The molecular geometries for each dataset were obtained from their respective source articles. Only the structures in the Organic dataset were optimized at the unrestricted ωB97XD/Def2-SVPD level of theory, and the XYZ Cartesian coordinates for all other datasets were used as provided. Single-point gas-phase calculations were performed at (unrestricted) ωB97X-D/def2-TZVPD level of theory for all datasets. Hybrid functionals combined with augmented triple-zeta basis sets have been shown to accurately predict molecular dipole moments in benchmark studies.^3,5,105 To assess the computational robustness of the Organic dataset, additional single-point gas-phase calculations were performed using Hartree–Fock (HF) and two density-functional theory (DFT) approximations (B3LYP^104,106,107 and ωB97XD¹⁰⁸) combined with the def2-SVPD, def2-TZVPD, and def2-QZVPD basis sets¹⁰⁹ All quantum chemistry calculations were performed with Gaussian16 (version C.01)¹¹⁰ to obtain the wavefunction and the HLYGAt and CM5 charges. Symmetry=None keyword was used in the Gaussian input files, which utilizes the (0, 0, 0) origin for the calculation of molecular multipole moments.

The resulting formatted checkpoint files were used to calculate various atomic charges with the free and open-source HORTON 2,¹⁰¹ ChemTools,^111,112 IOData,¹¹³ Grid,¹⁰³ GBasis,^114,115 and DensPart libraries. The proatom densities were computed at the same level of theory as the molecule on an “insane” Becke–Lebedev grid all datasets.^101,103 However, achieving a numerical integration accuracy of 10⁻³ (in atomic units) for molecular dipole components required different integration grid specifications for each dataset. We initially used the “insane” molecular integration grid for all databases, and if this accuracy threshold was not met, the grid size was further increased. For example, this setting was sufficient for the Silica dataset and most molecules in the Organic dataset, except the largest one (shown in Fig. 1), for which an “exp: 5 × 10⁻⁴ : 5 × 10¹ : 125 : 974” molecular grid was used. For the Backbone dataset, 38 (out of 400) N-Methylacetamide monomers required a larger grid; we used the “exp: 5 × 10⁻⁴ : 3 × 10¹ : 150 : 1202” to achieve the desired numerical accuracy. For the Sidechain datasets, only 129 (out of 790) monomers required an increased grid of “exp: 5 × 10⁻⁴ : 3 × 10¹ : 150 : 1202”, and the other 14 needed an “exp: 5 × 10⁻⁴ : 3 × 10¹ : 200 : 2030” specification. When computing atomic properties, all Hirshfeld methods, except MBIS, used the full molecular integration grid, eliminating additional numerical errors. To minimize the extended KL divergence in the AVH algorithm, we used the sequential least squares programming (SLSQP) algorithm (from SciPy¹¹⁶) with a termination tolerance of 1.0 × 10⁻⁶. The initial guess for proatom charges was chosen randomly between 0 and 1: this does not affect our results because the AVH optimization is convex, but making a random choice proves the numerical robustness of our procedure. For all Hirshfeld-based methods, the neutral and charged proatom densities were computed at the same level of theory as the molecule with an unrestricted wavefunction.

5 Conclusions

Capturing the anisotropy of electron density requires going beyond atomic charges to include higher-order atomic moments. Through a systematic evaluation of variational Hirshfeld–based schemes, this work assesses the robustness and practical usefulness of atomic moments for approximating molecular dipoles and quadrupoles across a diverse set of organic, inorganic, and protein fragment systems. We specifically highlight an often-overlooked aspect of computing atomic moments: before comparing them qualitatively or quantitatively, it is essential to verify their numerical accuracy, which depends on the molecular integration grid used.

Our analysis shows that when molecular dipoles are approximated using atomic charges, the AVH-B scheme outperforms other methods, including CM5, for the Silica dataset and exhibits comparable or slightly lower accuracy for the remaining datasets. Using atomic charges and dipole to approximate the molecular quadrupole, AVH-B delivers consistently superior performance, indicating smoother and faster convergence than HI and MBIS. Moreover, for most molecules, adding dipoles decreases the quadrupole error within the AVH framework, while HI and MBIS frequently show the opposite trend. These results position AVH-B as a promising framework for force-field development, enabling the derivation of higher-order electrostatic parameters directly from quantum-chemical calculations rather than through empirical fitting.

The computational robustness of our conclusions was tested using the Organic dataset evaluated at nine levels of theory. Atomic charges and their resulting molecular dipole and quadrupole approximations are largely stable across methods and basis sets, whereas atomic dipoles and quadrupoles show greater variability, with HI being the most sensitive. Consequently, molecular quadrupole approximations that include atomic dipoles exhibit stronger dependence on the level of theory. Conformational stability, assessed using alanine–dipeptide conformers, reveals that although molecular dipole norms change notably with dihedral rotation, atomic moments remain relatively unaffected. Charge-based dipole and quadrupole approximations vary smoothly with conformation, while inclusion of dipoles in the quadrupole approximation yields stable errors for AVH-B and MBIS but increased variability for HI.

Our results challenge the common assumption that the scheme providing the best atomic charges also produces the most reliable higher-order multipole moments; in fact, this is not necessarily the case. Based on the presented results, we propose AVH-B as a more reliable reference scheme for computing atomic moments to develop next-generation force fields and machine-learning models for atomistic simulations.

Author contributions

All authors reviewed and approved the final manuscript. Maximilian van Zyl contributed to methodology, software, data curation, formal analysis, visualization, writing – review & editing, investigation, and project administration; Carlos Castillo-Orellana contributed to data curation, formal analysis, visualization, validation, and writing – review & editing; Leila Pujal contributed to data curation, writing – review & editing; Farnaz Heidar-Zadeh contributed to conceptualization, methodology, formal analysis, visualization, writing – original draft, project administration, funding acquisition, and resources.

Conflicts of interest

There are no conflicts to declare.

Data availability

The “Computational Details” section describes the software and free/open-source code (with appropriate citations) used to build databases and compute atomic moments. All numerical data and input/output files from this study are available upon request from the corresponding author.

Supplementary information (SI): following the structure of the main paper, this includes additional plots that support our discussions, as well as an assessment of the computational robustness of our findings. See DOI: https://doi.org/10.1039/d5ra07866k.

Acknowledgements

M. V. acknowledges support from the Canada Graduate Scholarships awarded by the Natural Sciences and Engineering Research Council (NSERC) of Canada and Ontario Graduate Scholarship (OGS) program. C. CO. acknowledges the Emerging Leaders in the Americas Program (ELAP). F. H. acknowledges support from the Discovery Grant awarded by NSERC and Research Initiation Grant by Queen's University. This research was enabled in part by support provided by the Centre for Advanced Computing (CAC) at Queen's University (cac.queensu.ca), Compute Ontario (computeontario.ca), and the Digital Research Alliance of Canada (alliancecan.ca).

Notes and references

1. A. Stone, The Theory of Intermolecular Forces, Oxford University Press, 2013.
2. D. Hait and M. Head-Gordon, J. Chem. Theory Comput., 2018, 14, 1969–1981.
3. A. L. Hickey and C. N. Rowley, J. Phys. Chem. A, 2014, 118, 3678–3687.
4. F. De Proft, F. Tielens and P. Geerlings, J. Mol. Struct.:THEOCHEM, 2000, 506, 1–8.
5. J. C. Zapata and L. K. McKemmish, J. Phys. Chem. A, 2020, 124, 7538–7548.
6. Z. Jing, C. Liu, S. Y. Cheng, R. Qi, B. D. Walker, J.-P. Piquemal and P. Ren, Annu. Rev. Biophys., 2019, 48, 371–394.
7. S. Musse and R. J. Wheatley, Mol. Phys., 2023, 121, e2111375.
8. S. Vandenbrande, M. Waroquier, V. V. Speybroeck and T. Verstraelen, J. Chem. Theory Comput., 2017, 13, 161–179.
9. K. M. Visscher and D. P. Geerke, J. Chem. Theory Comput., 2019, 15, 1875–1883.
10. D. González, L. Macaya, C. Castillo-Orellana, T. Verstraelen, S. Vogt-Geisse and E. Vöhringer-Martinez, J. Chem. Inf. Model., 2022, 62, 4162–4174.
11. J. Pulido, L. Macaya and E. Vöhringer-Martinez, J. Chem. Eng. Data, 2024, 69, 2917–2926.
12. C. Castillo-Orellana, F. Heidar-Zadeh and E. Vöhringer-Martinez, J. Chem. Theory Comput., 2025, 21, 2043–2054.
13. J. Carmona-Espíndola, V. García-Melgarejo, E. Núñez-Rojas, S. Mendoza, A. García, J. L. Gázquez and J. Alejandre, J. Chem. Phys., 2024, 161, 144103.
14. P. L. A. Popelier and F. M. Aicken, ChemPhysChem, 2003, 4, 824–829.
15. C. Kramer, A. Spinn and K. R. Liedl, J. Chem. Theory Comput., 2014, 10, 4488–4496.
16. A. Spinn, P. H. Handle, J. Kraml, T. S. Hofer and K. R. Liedl, J. Chem. Theory Comput., 2020, 16, 4443–4453.
17. S. Van Damme, P. Bultinck and S. Fias, J. Chem. Theory Comput., 2009, 5(2), 334–340.
18. X. Guan, I. Leven, F. Heidar-Zadeh and T. Head-Gordon, J. Chem. Inf. Model., 2021, 61, 4357–4369.
19. D. S. Kosov and P. L. A. Popelier, J. Phys. Chem. A, 2000, 104, 7339–7345.
20. P. W. Ayers, R. C. Morrison and R. K. Roy, J. Chem. Phys., 2002, 116, 8731–8744.
21. F. De Proft, P. Geerlings, F. Heidar-Zadeh and P. W. Ayers, Conceptual Density Functional Theory, Elsevier, 2024.
22. W. T. Yang and W. J. Mortier, J. Am. Chem. Soc., 1986, 108, 5708–5711.
23. J. Olah, C. Van Alsenoy and A. B. Sannigrahi, J. Phys. Chem. A, 2002, 106, 3885–3890.
24. P. Bultinck, S. Fias, C. Van Alsenoy, P. W. Ayers and R. Carbo-Dorca, J. Chem. Phys., 2007, 127, 11.
25. F. Heidar Zadeh, P. Fuentealba, C. Cárdenas and P. W. Ayers, Phys. Chem. Chem. Phys., 2014, 16, 6019–6026.
26. E. Echegaray, S. Rabi, C. Cardenas, F. H. Zadeh, N. Rabi, S. Lee, J. S. M. Anderson, A. Toro-Labbe and P. W. Ayers, J. Mol. Model., 2014, 20, 2162.
27. E. Echegaray, C. Cardenas, S. Rabi, N. Rabi, S. Lee, F. H. Zadeh, A. Toro-Labbe, J. S. M. Anderson and P. W. Ayers, J. Mol. Model., 2013, 19, 2779–2783.
28. E. Echegaray, A. Toro-Labbe, K. Dikmenli, F. Heidar-Zadeh, N. Rabi, S. Rabi, P. W. Ayers, C. Cardenas, R. G. Parr and J. S. M. Anderson, in Correlations in Condensed Matter under Extreme Conditions: A Tribute to Renato Pucci on the Occasion of His 70th Birthday, ed. A. La Magna and G. G. N. Angilello, Springer International Publishing, Cham, Switzerland, 2017, pp. 269–278.
29. R. A. Miranda-Quintana, F. Heidar-Zadeh, S. Fias, A. E. A. Chapman, S. Liu, C. Morell, T. Gómez, C. Cárdenas and P. W. Ayers, Front. Chem., 2022, 10, 906674.
30. R. A. Miranda-Quintana, F. Heidar-Zadeh, S. Fias, A. E. A. Chapman, S. Liu, C. Morell, T. Gómez, C. Cárdenas and P. W. Ayers, Front. Chem., 2022, 10, 929464.
31. B. Wang, P. Geerlings, F. Heidar-Zadeh, P. W. Ayers and F. De Proft, J. Chem. Theory Comput., 2025, 21, 1695–1708.
32. B. Wang, P. Geerlings, C. Van Alsenoy, F. Heider-Zadeh, P. W. Ayers and F. De Proft, J. Chem. Theory Comput., 2023, 19, 3223–3236.
33. J. A. Semelak, I. Pickering, K. Huddleston, J. Olmos, J. S. Grassano, C. M. Clemente, S. I. Drusin, M. Marti, M. C. Gonzalez Lebrero, A. E. Roitberg and D. A. Estrin, J. Chem. Theory Comput., 2025, 21(10), 5194–5207.
34. M. J. Burn and P. L. A. Popelier, J. Chem. Theory Comput., 2023, 19, 1370–1380.
35. M. Thürlemann, L. Böselt and S. Riniker, J. Chem. Theory Comput., 2022, 18, 1701–1710.
36. M. Veit, D. M. Wilkins, Y. Yang, J. DiStasio, A. Robert and M. Ceriotti, J. Chem. Phys., 2020, 153, 024113.
37. M. T. Lehner, P. Katzberger, N. Maeder, C. C. Schiebroek, J. Teetz, G. A. Landrum and S. Riniker, J. Chem. Inf. Model., 2023, 63, 6014–6028.
38. R. Zubatyuk, J. S. Smith, J. Leszczynski and O. Isayev, Sci. Adv., 2019, 5, eaav6490.
39. M. G. Darley, C. M. Handley and P. L. A. Popelier, J. Chem. Theory Comput., 2008, 4, 1435–1448.
40. R. F. W. Bader, Atoms in Molecules: A Quantum Theory, Oxford University Press, 1990.
41. F. Heidarzadeh and S. Shahbazian, Int. J. Quantum Chem., 2011, 111, 2788–2801.
42. F. H. Zadeh and S. Shahbazian, Theor. Chem. Acc., 2011, 128, 175–181.
43. F. Heidar-Zadeh, P. W. Ayers, T. Verstraelen, I. Vinogradov, E. Vöhringer-Martinez and P. Bultinck, J. Phys. Chem. A, 2018, 122, 4219–4245.
44. A. Stone, Chem. Phys. Lett., 1981, 83, 233–239.
45. D. S. Kosov and P. L. A. Popelier, J. Phys. Chem. A, 2000, 104, 7339–7345.
46. H. Hu, Z. Lu and W. Yang, J. Chem. Theory Comput., 2007, 3, 1004–1013.
47. U. Koch, P. Popelier and A. Stone, Chem. Phys. Lett., 1995, 238, 253–260.
48. K. B. Wiberg and P. R. Rablen, J. Comput. Chem., 1993, 14, 1504–1518.
49. D. M. Elking, L. Perera and L. G. Pedersen, Comput. Phys. Commun., 2012, 183, 390–397.
50. R. F. W. Bader, A. Larouche, C. Gatti, M. T. Carroll, P. J. MacDougall and K. B. Wiberg, J. Chem. Phys., 1987, 87, 1142–1152.
51. B. Han, C. M. Isborn and L. Shi, J. Chem. Theory Comput., 2021, 17, 889–901.
52. Y. Shi, Z. Xia, J. Zhang, R. Best, C. Wu, J. W. Ponder and P. Ren, J. Chem. Theory Comput., 2013, 9, 4046–4063.
53. S. Y. Liem and P. L. A. Popelier, Phys. Chem. Chem. Phys., 2014, 16, 4122–4134.
54. P. Maxwell and P. L. A. Popelier, Mol. Phys., 2016, 114, 1304–1316.
55. K. Jug, Theor. Chim. Acta, 1975, 39, 301–312.
56. T. A. Manz, RSC Adv., 2020, 10, 44121–44148.
57. T. A. Manz, RSC Adv., 2022, 12, 31617–31628.
58. M. Cho, N. Sylvetsky, S. Eshafi, G. Santra, I. Efremenko and J. M. L. Martin, ChemPhysChem, 2020, 21, 688–696.
59. F. Heidar-Zadeh, C. Castillo-Orellana, M. van Zyl, L. Pujal, T. Verstraelen, P. Bultinck, E. Vöhringer-Martinez and P. W. Ayers, J. Chem. Theory Comput., 2024, 20, 9939–9953.
60. P. Bultinck, C. Van Alsenoy, P. W. Ayers and R. Carbó-Dorca, J. Chem. Phys., 2007, 126, 144111.
61. L. Pujal, M. van Zyl, E. Vöhringer-Martinez, T. Verstraelen, P. Bultinck, P. W. Ayers and F. Heidar-Zadeh, J. Chem. Phys., 2022, 156, 194109.
62. T. Verstraelen, S. Vandenbrande, F. Heidar-Zadeh, L. Vanduyfhuys, V. Van Speybroeck, M. Waroquier and P. W. Ayers, J. Chem. Theory Comp., 2016, 12, 3894–3912.
63. A. V. Marenich, S. V. Jerome, C. J. Cramer and D. G. Truhlar, J. Chem. Theory Comput., 2012, 8, 527–541.
64. R. F. Nalewajski and R. G. Parr, Proc. Natl. Acad. Sci. U. S. A., 2000, 97, 8879–8882.
65. R. F. Nalewajski and R. G. Parr, J. Phys. Chem. A, 2001, 105, 7391–7400.
66. P. W. Ayers, J. Chem. Phys., 2000, 113, 10886–10898.
67. R. G. Parr, P. W. Ayers and R. F. Nalewajski, J. Phys. Chem. A, 2005, 109, 3957–3959.
68. P. W. Ayers, Theor. Chem. Acc., 2006, 115, 370–378.
69. F. Heidar-Zadeh and P. W. Ayers, J. Chem. Phys., 2015, 142, 044–107.
70. F. L. Hirshfeld, Theor. Chim. Act., 1977, 44, 129–138.
71. F. Heidar-Zadeh, P. W. Ayers and P. Bultinck, J. Chem. Phys., 2014, 141, 094–103.
72. F. Heidar-Zadeh and P. W. Ayers, Theor. Chem. Acc., 2017, 136, 92.
73. F. Heidar-Zadeh, I. Vinogradov and P. W. Ayers, Theor. Chem. Acc., 2017, 136, 54.
74. F. HEIDAR-ZADEH and P. W. AYERS, Acta Phys.-Chim. Sin., 2018, 34, 514–518.
75. E. R. Davidson and S. Chakravorty, Theor. Chim. Acta, 1992, 83, 319–330.
76. T. C. Lillestolen and R. J. Wheatley, Chem. Commun., 2008, 45, 5909–5911.
77. T. C. Lillestolen and R. J. Wheatley, J. Chem. Phys., 2009, 131, 144101.
78. T. A. Manz and D. S. Sholl, J. Chem. Theory Comp., 2010, 6, 2455–2468.
79. D. Ghillemijn, P. Bultinck, D. Van Neck and P. W. Ayers, J. Comput. Chem., 2011, 32, 1561–1567.
80. T. Verstraelen, P. W. Ayers, V. Van Speybroeck and M. Waroquier, Chem. Phys. Lett., 2012, 545, 138–143.
81. T. A. Manz and D. S. Sholl, J. Chem. Theory Comp., 2012, 8, 2844–2867.
82. T. A. Manz and N. G. Limas, RSC Adv., 2016, 6, 47771–47801.
83. N. G. Limas and T. A. Manz, RSC Adv., 2016, 6, 45727–45747.
84. T. Verstraelen, P. W. Ayers, V. Van Speybroeck and M. Waroquier, J. Chem. Theory Comput., 2013, 9, 2221–2225.
85. J. Carmona-Espíndola and J. L. Gázquez, Comput. Theor. Chem., 2023, 1229, 114335.
86. J. E. S. Mikkelsen and F. Jensen, J. Chem. Theory Comput., 2025, 21, 1179–1193.
87. R. Benda, E. Cancès, V. Ehrlacher and B. Stamm, J. Chem. Phys., 2022, 156, 164107.
88. Y. Cheng, E. Cancès, V. Ehrlacher, A. J. Misquitta and B. Stamm, Multi-center Decomposition of Molecular Densities: A Numerical Perspective, 2024, https://arxiv.org/abs/2405.08455.
89. S. Fias, F. Heidar-Zadeh, J. S. M. Anderson, P. W. Ayers and R. G. Parr, J. Comput. Chem., 2018, 39, 1044–1050.
90. P. Bultinck, P. W. Ayers, S. Fias, K. Tiels and C. Van Alsenoy, Chem. Phys. Lett., 2007, 444, 205–208.
91. P. Bultinck, D. L. Cooper and D. Van Neck, PCCP, 2009, 11, 3424–3429.
92. E. Francisco, A. Martín Pendás, A. Costales and M. García-Revilla, Comput. Theor. Chem., 2011, 975, 2–8.
93. T. Bučko, S. Lebègue, J. G. Ángyán and J. Hafner, J. Chem. Phys., 2014, 141, 034114.
94. D. E. P. Vanpoucke, P. Bultinck and I. Van Driessche, J. Comput. Chem., 2013, 34, 405–417.
95. A. Krishtal, K. Vanommeslaeghe, A. Olasz, T. Veszprémi, C. Van Alsenoy and P. Geerlings, J. Chem. Phys., 2009, 130, 174101.
96. F. Heidar-Zadeh, P. W. Ayers and P. Bultinck, J. Mol. Model., 2017, 23, 348.
97. M. Riquelme, A. Lara, D. L. Mobley, T. Verstraelen, A. R. Matamala and E. Vöhringer-Martinez, J. Chem. Inf. Model., 2018, 58, 1779–1797.
98. M. Riquelme and E. Vöhringer-Martinez, J. Comput.-Aided Mol. Des., 2020, 34, 327–334.
99. T. Verstraelen, P. W. Ayers, V. Van Speybroeck and M. Waroquier, J. Chem. Theory Comp., 2013, 9, 2221–2225.
100. K. Carter-Fenk, M. Liu, L. Pujal, M. Loipersberger, M. Tsanai, R. M. Vernon, J. D. Forman-Kay, M. Head-Gordon, F. Heidar-Zadeh and T. Head-Gordon, J. Am. Chem. Soc., 2023, 145, 24836–24851.
101. M. Chan, T. Verstraelen, A. Tehrani, M. Richer, X. D. Yang, T. D. Kim, E. Vöhringer-Martinez, F. Heidar-Zadeh and P. W. Ayers, J. Chem. Phys., 2024, 160, 162501.
102. F. M. Aicken and P. L. Popelier, Can. J. Chem., 2000, 78, 415–426.
103. A. Tehrani, X. D. Yang, M. Martínez-González, L. Pujal, R. Hernández-Esparza, M. Chan, E. Vöhringer-Martinez, T. Verstraelen, P. W. Ayers and F. Heidar-Zadeh, J. Chem. Phys., 2024, 160, 172503.
104. A. D. Becke, J. Chem. Phys., 1988, 88, 2547–2553.
105. D. Hait and M. Head-Gordon, Phys. Chem. Chem. Phys., 2018, 20, 19800–19810.
106. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B:Condens. Matter Mater. Phys., 1988, 37, 785–789.
107. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
108. J.-D. Chai and M. Head-Gordon, Phys. Chem. Chem. Phys., 2008, 10, 6615–6620.
109. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297–3305.
110. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussiañ16 Revision C.01, Gaussian Inc. Wallingford CT, 2016.
111. F. Heidar-Zadeh, M. Richer, S. Fias, R. A. Miranda-Quintana, M. Chan, M. Franco-Pérez, C. E. González-Espinoza, T. D. Kim, C. Lanssens, A. H. Patel, X. D. Yang, E. Vöhringer-Martinez, C. Cárdenas, T. Verstraelen and P. W. Ayers, Chem. Phys. Lett., 2016, 660, 307–312.
112. L. Pujal, A. Tehrani and F. Heidar-Zadeh, in Conceptual Density Functional Theory: towards a New Chemical Reactivity Theory, ed. S. Liu, Wiley, 1st edn., 2022.
113. T. Verstraelen, W. Adams, L. Pujal, A. Tehrani, B. D. Kelly, L. Macaya, F. Meng, M. Richer, R. Hernández-Esparza, X. D. Yang, M. Chan, T. D. Kim, M. Cools-Ceuppens, V. Chuiko, E. Vöhringer-Martinez, P. W. Ayers and F. Heidar-Zadeh, J. Comput. Chem., 2021, 42, 458–464.
114. T. D. Kim, L. Pujal, M. Richer, M. van Zyl, M. Martínez-González, A. Tehrani, V. Chuiko, G. Sánchez-Díaz, W. Sanchez, W. Adams, X. Huang, B. D. Kelly, E. Vöhringer-Martinez, T. Verstraelen, F. Heidar-Zadeh and P. W. Ayers, J. Chem. Phys., 2024, 161, 042503.
115. A. Tehrani, M. Richer and F. Heidar-Zadeh, J. Chem. Phys., 2024, 161, 072501.
116. P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt and SciPy 1.0 Contributors, Nat. Methods, 2020, 17, 261–272.

---