An accurate description of electrostatics is one of the crucial prerequisites for modeling chemical interactions and developing force fields. Molecular electrostatics can be represented by a multipole expansion of atom-centered moments; however, such atomic moments are not uniquely defined and depend on the underlying fitting or partitioning scheme. Fitting approaches are appealing because they reproduce the electrostatic potential directly, but they are often ill-conditioned and non-unique, making the resulting atomic moments sensitive to the fitting protocol and less transferable across chemical environments. Density-based partitioning offers an attractive alternative by deriving atomic moments directly from the molecular electron density, although existing methods vary in how accurately they reproduce electrostatics. Basis-space approaches such as Gaussian Distributed Multipole Analysis (GDMA) [ 1 ] are computationally efficient but can perform poorly at short range, whereas Hirshfeld-type methods are more robust and transferable, though still strongly dependent on the choice of proatom densities. Bridging the gap between fitting and partitioning approaches, we introduce a new method for constructing atomic densities consistent with charges fitted to the electrostatic potential (ESP), analogous to a recently proposed reverse-engineering approach for Charge Model 5.[2] Specifically, we use the charge-constrained Additive Variational Hirshfeld (c-AVH) framework[3] to obtain atomic densities consistent with ESP-fitted charges. This new approach, called ESP-c-AVH, is applied to diverse datasets of organic molecules, inorganic species, and protein fragments. We compare the multipole moments produced by ESP-c-AVH with those obtained from several other methods, including ESP-fitted dipoles and quadrupoles[4 ], GDMA, and Hirshfeld-based partitioning schemes like Minimal Basis Iterative Stockholder (MBIS)[ 5]. We assess and compare their performance to reproduce molecular multipole moments and electrostatic potential. The results show that ESP-c-AVH moments perform competitively with MBIS and, in most cases, outperform GDMA and ESP-fitted moments.

References
[1] A. J. Stone, J. Chem. Theory Comput., 2005, 1, 1128–1132.
[2] F. Heidar-Zadeh and O. Hosseinzadeh, J. Phys. Chem. Lett., 2026, 17, 2841–2848.
[3] 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.
[4] C. Kramer, T. Bereau, A. Spinn, K. R. Liedl, P. Gedeck and M. Meuwly, J. Chem. Inf. Model., 2013, 53, 3410–3417.
[5] T. Verstraelen, S. Vandenbrande, F. Heidar-Zadeh, L. Vanduyfhuys, V. Van Speybroeck, M. Waroquier and P. W. Ayers, J. Chem. Theory Comput., 2016, 12, 3894–3912.