Seven confluence principles: a case study of standardized statistical analysis for 26 methods that assign net atomic charges in molecules

This article studies two kinds of information extracted from statistical correlations between methods for assigning net atomic charges (NACs) in molecules. First, relative charge transfer magnitudes are quantified by performing instant least squares fitting (ILSF) on the NACs reported by Cho et al. (ChemPhysChem, 2020, 21, 688–696) across 26 methods applied to ∼2000 molecules. The Hirshfeld and Voronoi deformation density (VDD) methods had the smallest charge transfer magnitudes, while the quantum theory of atoms in molecules (QTAIM) method had the largest charge transfer magnitude. Methods optimized to reproduce the molecular dipole moment (e.g., ACP, ADCH, CM5) have smaller charge transfer magnitudes than methods optimized to reproduce the molecular electrostatic potential (e.g., CHELPG, HLY, MK, RESP). Several methods had charge transfer magnitudes even larger than the electrostatic potential fitting group. Second, confluence between different charge assignment methods is quantified to identify which charge assignment method produces the best NAC values for predicting via linear correlations the results of 20 charge assignment methods having a complete basis set limit across the dataset of ∼2000 molecules. The DDEC6 NACs were the best such predictor of the entire dataset. Seven confluence principles are introduced explaining why confluent quantitative descriptors offer predictive advantages for modeling a broad range of physical properties and target applications. These confluence principles can be applied in various fields of scientific inquiry. A theory is derived showing confluence is better revealed by standardized statistical analysis (e.g., principal components analysis of the correlation matrix and standardized reversible linear regression) than by unstandardized statistical analysis. These confluence principles were used together with other key principles and the scientific method to make assigning atom-in-material properties non-arbitrary. The N@C60 system provides an unambiguous and non-arbitrary falsifiable test of atomic population analysis methods. The HLY, ISA, MK, and RESP methods failed for this material.


Introduction
Herein, statistical analysis is performed to better understand relationships among the large number of different methods for assigning net atomic charges (NACs) to atoms in molecules. Two related topics are explored. First, how do the relative charge transfer magnitudes of different NAC methods compare? Which NAC methods exhibit relatively small charge transfer magnitudes compared to other methods? Which exhibit relatively large charge transfer magnitudes? Second, which NAC method should be selected if the goal is to model a diverse set of properties related to NACs? For example, which NAC method assigns NACs having the overall strongest linear correlations to various other methods for assigning NACs?
Answering these questions requires an extensive dataset for statistical analysis. Cho et al. computed NACs for $2000 molecules and ions using 26 different charge assignment methods. 1 These charge assignment methods spanned many categories, including: (a) electron density partitioning into overlapping atoms, (b) electron density partitioning into nonoverlapping atoms, (c) NACs optimized to reproduce the molecular electrostatic potential (MEP), molecular dipole moment, or molecular dipole moment derivatives, (d) projection of the rst-order density matrix to give NACs having a complete basis set limit, (e) projection of the rst-order density matrix to give NACs having no complete basis set limit, and (f) various other schemes. The $2000 systems they studied were from the GMTKN55 database, which includes main group molecules and ions. 2 Cho et al.'s quantum chemistry calculations were performed using the PBE0 hybrid functional, 3,4 def2-TZVPP basis set, 5 and using geometries from the online GMTKN55 database 2 without further optimization. Their dataset comprises 29 934 atoms-inmolecules for which NACs were reported. 1 The present article studies the general question of how to design computed quantitative descriptors that are correlated to experimentally observed measured properties, where the computed quantitative descriptor itself is not unambiguously measurable experimentally for most materials. For most materials, the charge of an atom in the material is not itself unambiguously measurable experimentally. 6 Nevertheless, centuries of chemical science history show regarding some atoms in materials as positively charged (aka cations) and others as negatively charged (aka anions) is extremely useful for conceptually explaining chemical properties of materials. 7 Therefore, NAC is a useful computed quantitative descriptor for modeling or explaining experimentally observable properties such as molecular dipole moments, electric eld surrounding molecule, chemical reactivity, spectroscopic properties, etc. that are related to atom-in-material charges.
Is it possible to make any denite statements about how strongly correlated different NAC denitions are to any conceivable experimentally measured chemical property related to atom-in-material charges simply by studying statistical correlations in-between different NAC denitions even without knowing the experimentally measured chemical property to be modeled or explained? Surprisingly, I show herein the answer is yes. I derive a theory of conuence that shows some denitions for assigning NACs are positioned to produce average or better correlations to any and all conceivable properties related to atom-in-material charges. By the same reasoning, a bond order denition can be constructed that exhibits average or better correlations to any and all conceivable chemical properties related to bond orders. Accordingly, assigning properties to atoms in materials is not arbitrary.
More generally, this theory of conuence has transformative implications for all mathematical and physical sciences wherever the goal is to design a computed quantitative descriptor that is itself not a direct experimental observable (at least in most cases) but is correlated to a large number of experimentally observable properties. Conuence means a "joining together". Here, I show many statistical properties that were formerly considered distinct have strict equivalence or nearequivalence that eliminates much of the ambiguity in statistical analyses. Specically, the seven conuence principles explained herein show how to design a broadly applicable quantitative descriptor that exhibits average or better correlations to any and all conceivable related properties. Much like the theory of quantum mechanics that was developed in the twentieth century, this theory of conuence has profound and wide-ranging impacts that force us to interpret the world around us in new ways. This theory of conuence shows that dening quantitative descriptors that are not unambiguously measurable experimentally is still not an arbitrary process, because statistical correlations in-between possible alternative denitions determine which denition exhibits average or better correlations to any and all conceivable related properties.
The rest of this article is organized as follows. Section 2 explains the computational methods and theory behind them. Section 2.1 describes how the source data was checked for consistency to remove a small number of bad data points. Section 2.2 describes the rational and procedure for using a standardized reversible least squares tting called instant least squares tting (ILSF) to compute the relative charge transfer magnitudes of different charge assignment methods. Section 2.3 describes the principal components analysis (PCA) method. Section 2.4 presents mathematical theory governing maximally correlated descriptors. Section 3 presents computational results. Section 3.1 uses ILSF to quantify charge transfer magnitudes and explains atomic population method classication. Section 3.2 identies highly correlated descriptors using the correlation matrix and PCA applied to the NAC database. Section 3.3 presents results on the sensitivity of ranking to the choice of included charge assignment methods. Section 3.4 compares computed AIM populations for a benchmark system having unambiguous experimental values. Section 4 explains seven conuence principles that comprise the theory of conuence. Section 5 explains how these conuence principles work together with other key principles and the scientic method to make assigning atom-in-material properties nonarbitrary. Section 6 concludes. Section 7 contains several mathematical proofs.

Checking the source data for consistency
I checked the source NAC database 1 for consistency as follows. Because the correct net charge of every molecule or ion in the database is integer-valued, the running sum of NACs should reach an integer for the last atom-in-material of every molecule or ion in the database. The database was divided into blocks containing approximately 500 atoms-in-materials per block. Each block contained many molecules/ions, and each molecule/ ion belonged to only one block. (A system containing two molecules or ions spaced far apart (aka 'spatially separated') could be divided into two blocks, with one whole molecule or ion in each block.) For each charge assignment method, the running sum of NACs was computed for each block. For a particular block, the running sum should be equivalent between any two charge assignment methods.
Discrepancies between this expected behavior took three forms. First, some of the methods that computed NACs by numerical real-space integration had small, negligible integration errors; these NACs required no correction. Second, the MBSBickelhaupt NACs were missing for an extremely small number of atoms in materials. This occurred for a spatially separated Li + ion in four places, for which the MBSBickelhaupt NAC was manually set to +1. A [Li$(OH 2 )] + complex was missing MBSBickelhaupt NACs, so this system was entirely removed from the dataset for all charge assignment methods. Third, erroneous quantum theory of atoms in molecules (QTAIM) NACs were reported for a few systems. The spatially separated Li 2 (two occurrences), B 2 , C 2 , and P 2 (three occurrences) QTAIM NACs were manually set to zero, because they were erroneously reported to have large NACs (+0.26 to +0.65). Two systems containing 7 (i.e., H 3 Li 3 C) and 16 (i.e., H 7 BO 2 NaMg 2 Al 2 Cl) atoms were removed from all charge assignment methods, because their erroneously reported QTAIM NACs did not approximately sum to the system's net charge.

Instant least squares tting (ILSF)
Let {a i } and {b i } denote the NAC sets of two methods, where the subscript i runs over all atoms in materials. Standard deviations are computed in the usual manner: where M ¼ (N À 1) for a sample standard deviation and M ¼ N for a population standard deviation. 37 As described in standard statistics textbooks, the population standard deviation is computed from every datapoint in an entire population, while the sample standard deviation is computed when a data subset has been drawn from a larger population. 37 All equations in this article work whether the {s} correspond to sample or population standard deviations, but the same choice must be made for all regressed variables. Herein, the entire population of 29 907 atoms in materials were used to compute s (i.e., M ¼ N ¼ 29 907). The covariance matrix is dened as 38 If L ab ¼ 0, L ab > 0, or L ab < 0, the two variables a and b are said to be uncorrelated, positively correlated, or negatively correlated, respectively. 37 The covariance of a variable with itself is called that variable's variance: 37 The correlation matrix is dened as 37,38 From eqn (4), the covariance and correlation matrices equal each other when all variables have unit standard deviation: This can be achieved by standardizing the variables: 39 where Least-squares regression is a potential way to simultaneously quantify the relative charge transfer magnitudes and correlations between two methods for assigning NACs. Linear models could be constructed as If these two models are equivalent, then solving eqn (9) for b i yields an equation equal to eqn (10). Examining eqn (9) and (10), these two linear models are equivalent if We dene a reversible least-squares tting as one for which tting {a i } to {b i } (eqn (9)) yields a model equivalent to tting {b i } to {a i } (eqn (10)). Because simple least squares tting minimizes the results of simple least squares tting of {a i } to {b i } is not equivalent to tting {b i } to {a i }. For example, simple least squares tting yields the two inequivalent models where VDD are the VDD NACs, and QTAIM are the QTAIM NACs. The contradiction between these two models is obvious. Specically, solving eqn (13) for QTAIM gives QTAIM ¼ 6.0951 Â VDD À 0.0099, which does not even approximately equal eqn (14). The two approaches illustrated in Fig. 1 solve this problem. Both approaches minimize the squared deviations in both w and z variables: corresponds to the line between these two points being perpendicular to the model line.
Orthogonal distance regression was shown to be equivalent to a special case of total least squares regression. 40,41 Moreover, the resulting linear model for orthogonal distance regression corresponds to the major axis in principal components analysis (PCA). 38,41,42 Here, I show that by standardizing the independent variables it is possible to achieve a quadfecta for bivariate linear regression between any two positively correlated quantitative descriptors. Namely, the simultaneous accomplishment of: (1) orthogonal distance regression, (2) total least squares regression with Euclidean metric, (3) PCA regression, and (4) an instantaneous universal bivariate linear model. I now prove this instant least-squares tting (ILSF) can be achieved by standardizing the variables (eqn (6)-(7)), where s a ¼ 1 and s b ¼ sign(L ab ). If L ab ¼ 0, then w i and z i are uncorrelated, and the model collapses to the point (a model Otherwise, w i and z i are positively correlated and the ILSF yields the extremely simple linear model A remarkable property of eqn (16) is this linear model equation is identical for all conceivable pairs (â i ,b i ) of positively correlated real-valued standardized variables. That is, the same model equation describes the ILSF between any conceivable pair of real-valued positively correlated standardized quantitative descriptors in the universe. The name 'instant least squares tting' denotes the amazing result that the ILSF optimized linear model of eqn (16) can be written down instantaneously without having to perform computerized calculations. Section 7.1 below proves this ILSF model simultaneously optimizes the total least squares and orthogonal distance regression of the standardized variables.
ILSF is not the same as Deming regression. In Deming regression, deviations in the x and y variables are normalized by their measurement uncertainties (which approximately equal their root-mean-squared deviations from the model line). 42,43 In ILSF, standardized variables are used which normalize deviations in the x and y variables by the root-mean-squared deviations from their average values. Also, ILSF is not the same as a simple least-squares t on two standardized variables, because simple least squares tting yields irreversible models.

Principal components analysis (PCA)
PCA nds the eigenvalues and eigenvectors of the correlation and/or covariance matrices. 38,39,44 The principal components are sorted from highest to lowest eigenvalue. 38,39,44 The eigenvector having the largest eigenvalue is the rst (aka 'main') principal component. 38,39,44 The PCA eigenvectors are uncorrelated to each other (i.e., the covariance between any two different eigenvectors is zero). 38,39,44 This naturally follows from the fact that eigenvectors of any real symmetric matrix can be represented as an orthonormal basis. 38,45 If no eigenvalue is repeated (i.e., all eigenvalues are distinct), then the orthonormal eigenvectors are uniquely determined. 45 However, if two or more eigenvalues are equal, any rotation of the subspace formed from the corresponding eigenvectors yields new (and equally good) eigenvectors having the same eigenvalue. 38,45 For standardized variables, the correlation and covariance matrices are equal yielding unique results. For unstandardized variables, PCA of the correlation matrix is invariant to rescaling the variables, while PCA of the covariance matrix is not. 38 For example, consider PCA of three variables (A, B, C) compared to PCA of (A, B, D) where D is dened as 2C. PCA of the correlation matrix yields identical results for both variable sets, while PCA of the covariance matrix does not.
For PCA of the covariance matrix, the main principal component is the linear combination that results in the highest possible variance, subject to the normalization constraint where the subscript i represents a datapoint, the superscript (k) denotes which principal component (i.e., rst, second, third, etc.), the superscript j denotes which variable, and V is the total number of variables. 38 Because PCA of the covariance matrix is not scale invariant, it should only be used when the various Fig. 1 Geometry illustrating the error measures used in total least squares (approach 1) and orthogonal distance regression (approach 2). The red line represents the model equation. The green dot represents the measured datapoint. Approach 1 minimizes t 2 , and approach 2 minimizes h 2 .
variables are measured on a similar scale (e.g., all variables have the same measurement units). 38,39 For PCA of the correlation matrix, the eigenvalues sum to the total number of variables. 38 In this case, the eigenvalues represent how many standardized variables worth of variance are explained by each principal component. 38 For example, an eigenvalue of 10.3 means that principal component explains as much variance as 10.3 standardized variables. A principal component with an eigenvalue less than one represents less variance than one standardized variable. The goal of PCA is to reduce the number of variables required to explain the data. For PCA of the correlation matrix, the square root of the variance of standardized variables explained by the k th principal component (PCk) expands as where v a (k) is the coefficient for standardized variable w (a) in the k th eigenvector of the correlation matrix, and l (k) is corresponding eigenvalue. Because the PC's are normalized inserting eqn (19) into eqn (20) gives

Maximally correlated descriptors
This article focuses on conuence principles for a group of mutually positively correlated descriptors. A set of quantitative descriptors is mutually positively correlated if and only if all elements in the correlation matrix are positive and non-zero which is equivalent to all elements in the covariance matrix being positive and non-zero. Although these could represent different experimentally measurable physical properties, the focus in this article is on computed quantitative descriptors that are correlated to many experimentally measurable physical properties but are not themselves uniquely measurable experimentally for most situations. Net atomic charges are a prime example. Centuries of chemical sciences history establish the charges of atoms in materials as a fruitful concept for explaining many chemical phenomena, but various different ways to assign NACs can be conceived. How does one determine the most suitable denitions for broad use? A denition suitable for broad use should be simultaneously correlated to the various physical properties related to that concept. For example, a NAC denition suitable for broad use should be simultaneously correlated to the experimentally measured chemical properties that are related to the concept of charges of atoms in materials. Such a denition would be a superdelegate that captures the essence of the group of mutually positively correlated descriptors. Because the experimentally measured chemical properties closely related to the concept of charges of atoms in materials must be strongly correlated to some particular NAC denition(s), the superdelegate can be chosen by identifying the group member that maximizes the sum of correlations to group members: The average standardized variable at datapoint i is Standardizing this descriptor yieldŝ The sum in eqn (23) can be expanded as where U(a, f) is the correlation between a and f. Hence, the group member that maximizes the sum of correlations to all group members is the group member that is maximally correlated to the average standardized variable. As a further performance characteristic, we can ask how correlated this average is to all group members Inserting eqn (27) into eqn (29), this simplies to Combining eqn (28) and (30) gives the correlation between standardized variableâ and the average standardized variable f: which quanties the relative ability of standardized variableâ to serve as a delegate for the mutually positively correlated descriptors group. Section 7.2 below proves that f maximizes possible summed correlations to the variables {â}. That is, S f $ S s for any conceivable descriptor s that is a linear combination of the standardized variables.
How is f related to the main principal component (MPC) of the correlation matrix? The MPC is the eigenvector with the largest eigenvalue. By denition, a matrix times one of its eigenvectors yields the corresponding eigenvalue (a scalar) times that eigenvector. A common method to nd the principal eigenstate is the identity (32) where (ñ max , l max ) is a principal eigenvector and its eigenvalue.
In eqn (32), p and p À 1 are powers of the matrix. However, eqn (32) only holds if the trial vector is not orthogonal toñ max : Sincef is maximally correlated to the descriptor group's variables, it is a good initial guess forñ max . Substitutingf forñ trial in eqn (32) yields the rst renement Hence, it follows that the coefficient for standardized variableâ in the MPC for PCA of the correlation matrix is approximately proportional to S a (i.e., its summed correlation to all variables in the descriptor group). Since the MPC is normalized, this means each variable's coefficient in the MPC is approximately given by Accordingly, the order of coefficients (largest to smallest) in the MPC of the correlation matrix is approximately the same order as S a (largest to smallest). Repeated renement via eqn (32) could potentially lead to subtle differences between these two orders, but it is highly unlikely that a bottom 25% variable according to the S a criterion would become a top 25% variable according to its coef-cient in the MPC of the correlation matrix, and vice versa.
This analysis clearly reveals a close link between PCA of the correlation matrix, highly correlated descriptors, the average standardized variable f, and the superdelegate. Specically, the superdelegate is the descriptor from the group that has the highest correlation to all group members, and it is likely to have the largest coefficient in the MPC of the correlation matrix. Consequently, this superdelegate will also have relatively high correlation to the MPC of the correlation matrix. Moreover, this superdelegate has high correlation to f, and f has high correlation to the MPC of the correlation matrix.

Charge transfer magnitudes and atomic population method classication
Because the average charge transfer magnitude and the correlation matrix are completely independent of each other, both should be considered when assessing the statistical performance of different charge assignment methods. It is possible to have high statistical correlation between two charge assignment methods even though they predict vastly different charge transfer magnitudes. Theoretically, one of these two charge assignment methods could predict reasonable charge transfer magnitudes while the other might severely under-estimate or over-estimate charge transfer magnitudes. This could occur even if the correlation between the two methods is essentially 1.00. Consider two hypothetical methods that assign NACs directly proportional to each other. For example, A ¼ 5B. The correlation matrix is unchanged if A is swapped for B. In contrast, the average charge transfer magnitude is directly affected by a scaling factor. In this example, method A has ve times the charge transfer magnitude of method B.
The charge transfer magnitude of each NAC method was quantied by its root-mean-squared (rms) deviation from its average value (i.e., s as dened in eqn (1)). Table 1 also lists the average charge, q avg , for each method across the 29 907 atomsin-molecules. The small q avg differences are due to integration imprecisions. There was a factor of 4.9 between the methods with smallest (i.e., Hirshfeld) and largest (QTAIM) charge transfer magnitudes for molecules. To make the results easier to interpret, the fourth column lists s/s DDCE6 as the relative charge transfer magnitude.
The h column of Table 1 indicates whether the NACs have a mathematical limit as the basis set is improved towards completeness. Individual atom-in-material descriptors (e.g., net atomic charges, atomic spin moments (ASMs), bond orders, spdfg populations, polarizabilities, etc.) only have clear chemical and physical meaning when they converge to well-dened values as the basis set is improved (i.e., they have complete basis set limits). Therefore, population analysis methods lacking a complete basis set limit are not useful for computing these properties. Regardless of whether or not an atomic population analysis method has a complete basis set limit, it can still act as a useful basis representation to expand quantum mechanical operators. For example, the electron-electron Coulomb electrostatic energy of a material can be expressed exactly as a polyatomic multipole expansion plus charge overlap terms. This Coulomb energy can be expanded exactly using any population analysis method that reproduces the material's electron distribution, irrespective of whether that population analysis method has a complete basis set limit. However, when the population analysis method lacks a complete basis set limit it is only the computed coulombic energy and not the individual populations that carry any physical meaning. Consequently, individual values of atom-in-material descriptors reported in scientic publications should be computed using methods having a complete basis set limit.
A complete basis set limit is a necessary but not a sufficient condition for computing highly valuable atom-in-material descriptors. Several other criteria are also required: (a) the atom-in-material descriptor values should be highly correlated to many experimentally measured properties, (b) the population analysis method should yield a correct atom-in-material descriptor value for carefully chosen benchmark systems having well-known and unambiguous atom-in-material properties, and (c) the population analysis method should yield atom-in-material descriptor values that are chemically consistent amongst themselves (e.g., the NAC value should be chemically consistent with the ASM value 15 ). These and related criteria are explained more fully in Section 5.
Every quantum chemistry calculation in which the electron density is properly computed from a wavefunction yields a nonnegative total electron density r r $ 0 (Caution: Quantum chemistry algorithms that merely estimate the electron density using response theory may yield rðrÞ\0 for some positionr; this does not correspond to the proper electron density of any wavefunction.) The sixth column of Table 1 indicates whether each method partitions the total electron density into non-negative atom-in-material electron densities having a complete basis set limit. If yes, then partitioning into overlapping versus non-overlapping fr A ðrÞg is indicated. "No" means either that the electron density is not partitioned, that some partitions can have negative density values at some spatial positions, or that the method lacks a complete basis set limit. The electron density partitions in eqn (37) usually correspond to atoms, but the QTAIM method can have some non-nuclear attractors (i.e., one or more electron density partitions that are not atoms). 46,47 Such non-nuclear attractors are a modeling No basis set or quantum chemistry calculation is required to compute EEQ NACs. b Many different charge electronegativity equilibration schemes have been proposed. Many of these are not robust, because they sometimes produce extremely high NAC magnitudes. c The IBO method currently requires the rst-order density matrix to be idempotent. d Whether or not the RESP NACs have a complete basis set limit depends on the type of tting constraints used. If and only if the tting constraints have no basis set dependence or have a complete basis set limit, then the corresponding RESP NACs will have a complete basis set limit. Whether the RESP NACs are robust depends on how the constraints are constructed. e Not rotationally invariant. f Methods that project populations from a quantum chemistry calculation basis set (aka 'source basis set') onto a small basis set (aka 'target basis set') have a basis set limit with respect to improving the source basis set towards completeness, but their results depend on the small target basis set onto which the populations are projected. g QTAIM partitions are robust only when they have been sufficiently smoothed so that noise does not create spurious virial compartments.
advantage for electrides but can be a modeling disadvantage for other materials. 15,48 The seventh column in Table 1 briey summarizes the charge assignment strategy. The Hirshfeld and VDD methods partition the molecule's deformation density into overlapping and non-overlapping partitions, respectively. Methods marked "1PDM projection" project components of the one-particle density matrix (1PDM). Methods marked "dipole intent" were developed to approximately reproduce the molecular dipole moments of reference compounds. "Dipole t" indicates the NACs are optimized to reproduce each molecule's quantummechanically computed dipole moment. The EEQ method requires no quantum chemistry calculation. "MEP t" indicates the NACs minimize some error measure between the quantummechanical molecular electrostatic potential (MEP) and the electrostatic potential of the NAC model; these methods may differ by the grid points and integration weights used to construct the error measure. The DDEC6 method optimizes NACs to simultaneously give small errors across both electrostatic and chemical properties. The APT method optimizes the NACs to reproduce changes in the molecular dipole moment as the atoms vibrate, assuming each NAC is constant as the molecule vibrates. 10 Entries marked "reference orbitals", "reference ions", "reference radii", "spherical averaging", and "Slater functions" indicate a key feature of the charge assignment scheme. The QTAIM method assigns non-overlapping virial compartments. 32,[49][50][51] Cho et al. misclassied the DDEC6 method as an "iterative Hirshfeld variant" (page 694 of ref. 1), which it is not. The Hirshfeld and VDD approaches are based on deformation density partitioning using overlapping and non-overlapping compartments, respectively. 21,36 As shown in Table 1, deformation density approaches yield the lowest average charge transfer magnitudes of all charge assignment methods. The iterative Hirshfeld (aka Hirshfeld-I) method was developed by Bultinck et al. and sets the atomic weighting function equal to a quantummechanically computed reference ion density, where the reference ion's charge is self-consistently updated to match the assigned AIM charge. 25 The earliest DDEC methods used a combination of spherical averaging and charge-compensated reference ions for which the reference ion charges were self-consistently updated to match the assigned AIM charges. 52 Unfortunately, the Hirshfeld-I and early DDEC methods suffer the runaway charges problem in which vastly different NACs are sometimes assigned to symmetry equivalent atoms in materials. 15,53 The DDEC6 method uses a xed sequence of seven charge partitioning steps to solve the runaway charges problem. 15 DDEC6 is the sixth generation improvement of the Density-Derived Electrostatic and Chemical (DDEC) methods. 15,54-56 DDEC6 uses: (a) tail constraints on the atomic weighting functions to prevent them from becoming too diffuse or contracted for buried atom tails, (b) reference ion charges that approximate the number of electrons in the volume dominated by each atom, (c) reference ion smoothing and conditioning to allow the reference ions to expand or contract according to the material's local environment, (d) a weighted spherical average to more accurately reproduce the electrostatic potential surrounding the material, and (e) a xed sequence of seven charge partitioning steps to avoid the runaway charges problem. 15,53 The last column in Table 1 includes comments on specic convergence issues. Methods that can converge to vastly different solutions depending on the initial guess do not have a convex optimization functional for some materials; the Hirshfeld-I and MBIS methods are such examples. 15,27 Methods with a convex optimization landscape that is nearly at for buried atoms can assign buried atom charges that are not chemically meaningful; the CHELPG, HLY, ISA, and MK electrostatic potential tting methods are such examples. 23,33,52 Many different charge electronegativity equilibration schemes have been proposed. [16][17][18]20,[57][58][59][60][61] Many charge electronegativity equilibration schemes sometimes produce extremely high NAC magnitudes. 58,61,62 The ACP and i-ACP NACs are sensitive to the choice of valence electrons for each chemical element; for example, vastly different results might be obtained depending on whether Cs element is considered to have one (i.e., 6s 1 ) or nine (i.e., 5s 2 5p 6 6s 1 ) valence electrons. This unfortunate dependency arises, because the ACP and i-ACP methods are dened to t the entire valence electron population of an atom-in-material using only one Slater exponential decay function. 8,24 [CsO 4 ] + has strong polarcovalent bonding between the Cs and O atoms not purely ionic bonding. 63 In [CsO 4 ] + , the 5s and 5p 'semi-core' electrons are key participants in the polar-covalent bonding, thus acting as valence electrons along with higher subshells. 63 Examining Table 1, the deformation density methods (i.e., Hirshfeld and VDD) had the smallest charge transfer magnitudes, while partitioning based on Virial compartments (i.e. QTAIM) had the largest. Methods designed to approximately (i.e., ACP, CM5, i-ACP) or exactly (i.e., ADCH) reproduce the molecular dipole moment had larger average charge transfer magnitudes than the deformation density group but smaller than the MEP tting group (CHELPG, RESP, MK, HLY). The DDEC6, IBO, Bickelhaupt, and ISA methods had average charge transfer magnitudes similar to the MEP tting group. Many methods (e.g., Hirshfeld-I, MBIS, Becke, APT, etc.) had average charge transfer magnitudes larger than the MEP tting group.
As an illustrative example, Table 2 summarizes selected calculations for the water molecule. Water was chosen for two reasons. First, it participates in many biological, environmental, geological, and chemical processes. Second, its three-atom bent geometry permits NACs to be directly derived from its calculated molecular dipole moment. This corresponds to the ADCH oxygen NAC of À0.693. Larger molecules containing more than two distinct atom types do not have uniquely determined NACs derived only from the molecule's dipole moment, because multiple NAC values could reproduce the same molecular dipole moment. The CM5 oxygen NAC of À0.642 was slightly smaller in magnitude than the ADCH value. All four MEP tting methods (CHELPG, HLY, MK, and RESP) yielded practically identical oxygen NAC of À0.715 to À0.704. Moreover, the oxygen NAC that minimized the RMSE over the 788833 grid points for data listed in Table 2 was also within this same range. The DDEC6 (À0.802) and Hirshfeld-I (À0.900) oxygen NACs were somewhat larger in magnitude than the MEP tting group. As expected, the deformation density (i.e., Hirshfeld and VDD) NACs were too small in magnitude to approximate the molecular dipole moment or MEP. Also as expected, the QTAIM NACs were too large in magnitude to approximate the molecular dipole moment or MEP. When atomic dipoles are included, the molecular dipole moment is reproduced exactly.
The data in Table 2 were computed as follows. The optimized molecular geometry, electron density distribution, and reference electrostatic potential were computed using GAUSSIAN 16 (ref. 64) soware. The dipole moment magnitude of the computed PBE0/def2TZVPP optimized geometry and electron density was 0.765 au, which was used as the reference dipole moment. Using an in-house program, the RRMSE was computed over a uniform grid of 788833 points between 1. The ACP and i-ACP charges were computed using the ACP 8,68 and i-ACP 24,69 programs. Although the ACP and i-ACP methods could potentially be used to compute atomic dipoles, these were not available in the soware versions used.

Identifying highly correlated descriptors
Fig. 2 displays the correlation matrix between all 20 methods having a complete basis set limit. As explained in Section 2.2 above, this also equals the covariance matrix of the standardized variables. Fig. 2 is related to Only the ADCH, APT, Becke, QTAIM, and VDD methods have correlation < 0.9 to the DDEC6 method. Excluding Hirshfeld, the remaining 14 methods connected to DDEC6 form the main block. ADCH is almost connected to the main block through the ADCH-CM5 correlation ¼ 0.8999, but it has no correlation $ 0.9 to any method except self. A small side block containing the deformation density methods (Hirshfeld and VDD) is connected to the main block only through the Hirshfeld-DDEC6 correlation ¼ 0.908. Another small side block containing i-ACP (also part of the main block), APT, and QTAIM is connected to the main block through i-ACP. Within the main block, DDEC6 is the most connected (15 correlations $ 0.9) and IBO and EEQ are the least connected (each having 6 correlations $ 0.9).
Several atomic population analysis methods optimize similarity between atom-in-material electron distributions and those of quantum-mechanically computed reference atoms. These methods require a library of quantum-mechanically computed reference atoms. Among the 26 atomic population analysis methods considered here, these include DDEC6, Hirshfeld, Hirshfeld-I, and IBO. Only the neutral uncharged ground-state reference atoms are required for the Hirshfeld and IBO methods, while ground-state reference ions in various charge states are required for the DDEC6 and Hirshfeld-I methods. The IBO method uses an ingenious projection to represent the molecular orbitals in terms of polarized atom-inmaterial orbitals. 22 Currently, the IBO method is limited to idempotent density matrices. 22 The DDEC methods use chargecompensated reference ions and reference ion conditioning to polarize reference ions by their material environment. 15,52,53 The similar average charge transfer magnitudes (see Table 1) of DDEC6 and IBO NACs is notable. As shown in Fig. 2, the correlation between DDEC6 and IBO NACs is 0.954. Although the Hirshfeld and IBO methods are both based on neutral uncharged reference atoms, their average charge transfer magnitudes differ by a factor of 2.5. The average charge transfer Table 2 Relative root mean squared errors (RRMSE) in electrostatic potential of the water molecule for 20 charge assignment methods having a complete basis set limit. Errors in the predicted molecular dipole moment magnitude are also listed. Methods listed from smallest to largest NAC magnitude on oxygen. For some of the non-negative AIM density partitioning methods, the errors including atomic dipoles are listed in parentheses

Method
Oxygen magnitude of the Hirshfeld-I method is about 1.2 times that of the DDEC6 method. Correlation between the DDEC6 and Hirshfeld-I methods is high at 0.988. In spite of the similar names, the Hirshfeld NACs have slightly less correlation to the Hirshfeld-I NACs (0.879) than to both the IBO (0.893) and DDEC6 (0.908) NACs. Table 3 summarizes PCA of the correlation matrix. All methods had positive coefficients in the MPC. The 14 methods of the main block had the largest MPC coefficients. Comparing columns 2 and 6 of Table 3 shows the approximation of eqn (35) is almost exact. The eigenvalue shows the MPC accounts for 17.158 variables' (85.8%) worth of correlation. This clearly reects the size of the main block (14 methods) plus some contributions from small side blocks weakly connected to the main block. The other principal components (i.e., PC2, PC3, PC4, etc.) account for less than one variable's worth of correlation apiece.
Conuence, which will be more thoroughly explained in Section 4 below, occurs when a quantitative descriptor yields high correlations across a broad group of related descriptors and physical properties. Here, there are three key indicators that conuence occurs among NAC descriptors. The rst, and perhaps most important, is the MPC accounts for the vast majority (i.e., 85.8%) of correlation within this NAC descriptor group. The second is that each of the remaining principal components is extremely weak, accounting for less than one variable's worth of correlation apiece. The third is that at least one individual NAC descriptor (e.g., DDEC6) is highly correlated to a large percentage (e.g., 15/20 ¼ 75%) of NAC descriptors. Of course, the correlation matrix's pronounced main block illustrated within Fig. 2 is a consequence of these three factors.
Since conuence is present within this group of NAC descriptors, it is useful to ask: "Which individual member of the group best represents the group as a whole?" Different criteria can be conceived to determine this: (1a) the member having the largest coefficient in the correlation MPC, (1b) the member having the highest correlation to the correlation MPC, (2a) the member having the highest summed correlation to all group members (i.e., largest S a ), (2b) the member having the highest correlation to the average standardized variable f (i.e., largest U (af) ), or (3) the member having strong correlations to the largest number of other group members. Because criteria (1a) and (1b) are proportional to each other (see conuence principle # 2 of Section 4), they always give identical rankings. Because criteria (2a) and (2b) are proportional to each other (see eqn (28)), they always give identical rankings. By eqn (35), rankings according to criteria (1a) and (2a) will be similar but not necessarily identical. Clearly, a member that has strong correlations to a large number of other group members must also have a relatively high S a ; therefore, criteria (2a) and (3) oen give somewhat similar results. Consequently, in practice the results are oen similar irrespective of which criterion is chosen. Table 3 ranks NAC methods according to criterion (1a). Table  4 ranks NAC methods according to criterion (1b), (2a), (2b), and (3). Rankings for criterion (3) were performed separately using two different thresholds: the numbers of NAC methods having correlation $0.8 and $ 0.9 to each method. The top (DDEC6), the bottom (Becke), the 2 nd (MBIS), the 8 th (RESP), the 9 th (MK), the 11 th (CM5), and the 15 th (Hirshfeld) ranked methods had consistent rankings across all ranking criteria. Across the different ranking criteria, small variations in the placements of other methods were observed.
Cho et al. previously reported MBSBickelhaupt and Hirshfeld-I as having largest correlation to the unstandardized covariance MPC among 16 NAC methods. 1 Because the unstandardized covariance matrix is sensitive to multiplying a variable by a scale factor, PCA of the unstandardized Fig. 2 Correlation matrix between 20 methods having a complete basis set limit for assigning net atomic charges in molecules. Stoplight colors indicate the covariance values: green $ 0.9, 0.8 # yellow < 0.9, red < 0.8. Blue shading marks blocks of values $ 0.9. There are three primary groups: (a) a main group that covers a large number of methods, (b) the i-ACP, APT, and QTAIM group, and (c) the VDD and Hirshfeld group. The DDEC6 method is strongly correlated to all members of group (a) plus the i-ACP method in group (b) and the Hirshfeld method in group (c). No other charge assignment method besides DDEC6 is strongly correlated to some members of all three groups. The ADCH and Becke methods are not strongly correlated to any charge assignment methods besides self. The ADCH-CHELPG entry is red rather than yellow, because its value is 0.7996 which is below the 0.8 cutoff. The ADCH-CM5 entry is yellow rather than green, because its value is 0.8999 which is below the 0.9 cutoff.
covariance matrix tends to favor contributions from variables having larger s, when compared to PCA of the correlation matrix. Because our goals in this paper are to examine the charge transfer magnitudes and correlation properties, we refer readers to Cho et al.'s work 1 for a detailed discussion of PCA of the unstandardized covariance matrix. Table 3 The first four eigenvalues and principal components coefficients for correlation PCA of 20 charge assignment methods having a complete basis set limit. The methods are listed in order from largest to smallest contribution to the MPC. The last column is listed for comparison to the MPC coefficient of column 2 Returning to a discussion of Table 4, it is instructive to ask how well DDEC6 performs compared to the correlation MPC and compared to the average standardized variable f. Among all conceivable descriptors, S max ¼ S f ¼ 18.4806 is the highest possible sum of correlations to the 20 NAC methods. The sum of correlations between the MPC and the NAC methods is S MPC ¼ 18.4795 and almost as high as S f . The S DDEC6 ¼ 18.204 is 0.985 times S f , which also equals the correlation between DDEC6 NAC and f. Correlation between DDEC6 and correlation MPC is almost the same at 0.986. Hence, the DDEC6 NAC captures much of the same information that is captured by f and the correlation MPC.
Examining other high performing methods, the top four ranked methods have S DDEC6 À S a < S max À S DDEC6 , while this inequality does not hold for methods ranked h and beyond. Hence, the top four ranked methods (i.e., DDEC6, MBIS, ISA, and Hirshfeld-I) have relatively small differences between their S a values. Consequently, the average charge transfer magnitudes should also be considered when selecting among these four methods. Among these four methods, the average charge transfer magnitudes from Table 1

Sensitivity of ranking to the choice of included methods
A key question is "How robust are the rankings of the topranked methods to changes in which other methods are included in the dataset?" For example, what happens if the dataset is spammed with trivial variations of one charge assignment method? For example, electrostatic potential tting methods such as CHELPG, HLY, and MK differ only in the choice and weighting of grid points on which the root mean squared error (RMSE) of the electrostatic potential is computed and minimized. With slightly different choices in the grid points and their weightings, one could easily produce a thousand slightly different variations of electrostatic potential tting methods. If these are included in the dataset would they force one of the electrostatic potential tting methods into the topranked position? Somewhat surprisingly, the answer is no. Spamming the database with trivial variations of one method is not sufficient to elevate that method into the top-ranked position if the method being spammed is highly correlated (i.e., U ab > 0.9) to the top-ranked method. For some of the ranking criteria, the top-ranked charge assignment method does not change under such a scenario. Examining Fig. 2 and Table 4, the DDEC6 method is highly correlated (i.e., U ab > 0.9) to 15 charge assignment methods, including the CHELPG method which is highly correlated to 9 charge assignment methods. If 1000 new electrostatic potential tting methods that are trivial variations compared to CHELPG are added to the dataset, this increases the number of methods highly correlated to both DDEC6 and CHELPG by exactly 1000. The new numbers of highly correlated methods (1015 to DDEC6 and 1009 to CHELPG) do not change the relative order of these two methods at all. This new dataset yields rankings identical to the original dataset for each of the top 16 methods according to the U ab > 0.8 ranking criterion and for each of the top 4 methods according to the U ab > 0.9 ranking criterion.
Moreover, such spamming can be easily detected by a ranking abnormality. In the above example, S a for CHELPG would increase from 17.600 in the original dataset to 1017.600 in the modied dataset, while S a for DDEC6 would increase from 18.204 to (18.204 + 1000 Â 0.9252) ¼ 943.4. The better ranking of DDEC6 than CHELPG for number of methods with U ab > 0.9 and U ab > 0.8 but worse ranking of DDEC6 compared to CHELPG for S a in the modied dataset is a clear indication the modied dataset contains a cluster of methods highly similar to CHELPG which are not as conuent as DDEC6 across the entire database. In other words, this ranking abnormality (i.e. different top-ranked method for S a criterion compared to U ab > 0.9 criterion) makes the spamming obvious and easy to detect.
What happens if the method being spammed is low-ranked in the original dataset? For example, if 1000 trivial variations of the Becke method were added to the dataset? Since the Becke method has low correlations to all of the other charge assignment methods in the original dataset, this spamming would force all of these trivial variations of the Becke method into the top-ranked positions of the modied dataset for any of the ranking criteria used in Table 4. However, it would be easy to detect this sham conuence. When genuine conuence occurs, the conuent method exhibits conuence not only across various computed descriptors but also across the physical properties those computational descriptors are intended to describe. Although the Becke method performed well for the water molecule (see Table 2), it gave the wrong sign and magnitude of NAC for Eu in [Eu@C 60 ] + (see Table 7). Specically, the Becke NAC of À4.427 for the Eu atom in [Eu@C 60 ] + is chemically wrong. Also, Table 1 shows the average charge transfer magnitude of the Becke method is relatively high compared to methods optimized to reproduce the MEP.
Another important question is whether the rankings would remain similar if new methods are added to the dataset that are not trivial variations of the already included methods. Moreover, will the rankings be adversely affected if the quality of these newly added methods is dubious? To address this question, the dataset is re-analyzed by adding the six charge assignment methods that do not have a complete basis set limit. Comparing the new rankings listed in Table 5 to the original rankings in Table 4, the DDEC6 and MBIS methods remain in the rst and second spots, respectively, for all of the metrics. The MBSBickelhaupt method (which is one of the newly added methods) is now in the third spot according to all the metrics. Hirshfeld-I now places fourth according to all the metrics, while it originally placed fourth according to all the metrics except one for which it originally placed third. This analysis shows the relative rankings of the methods are only weakly affected by adding a modest number of new methods, even if those new methods are of dubious quality.
Another useful question is whether the data for one particular method has the potential ability to dramatically alter the rankings. A way to frame this question is to ask how the rankings could potentially change if one of the charge assignment methods in the original dataset is swapped for a new charge assignment method having any conceivable properties. Examining Table 4, the number of (U ab > 0.9) ranking criterion is the most robust to this kind of method swap. For charge assignment method A, swapping one of the other charge assignment methods (B) for an arbitrary new one (B 0 ) could affect the number methods having (U ab > 0.9) to method A by: (i) +1 if method B 0 is highly correlated to method A while method B is not, (ii) by À1 if method B is highly correlated to method A while method B 0 is not, and (iii) otherwise this number will be unchanged by the swap. Examining Table 4, a change in AE1 in the number of (U ab > 0.9) for each method would leave DDEC6 and MBIS in either the rst or second spots. Hence, any conceivable change to a single charge assignment method only has a small potential impact on the (U ab > 0.9) ranking criterion.
Finally, consider the grouping of methods into families of related methods. The electrostatic potential tting family includes CHELPG, HLY, MK, and RESP. The deformation density family includes Hirshfeld and VDD. Stockholder partitioning methods include a diverse set that spans a wide variation in average charge transfer magnitudes: Hirshfeld, ACP, i-ACP, DDEC6, ISA, Hirshfeld-I, MBIS, and Becke. Although from a methodology perspective the stockholder partitioning methods form a class, their charge assignment results are diverse. For example, Hirshfeld NACs are highly correlated to VDD NACs (both are based on the deformation density) but not to the Hirshfeld-I NACs. 1 From a statistical perspective, DDEC6 NACs were very highly (>0.95) correlated to MBIS, Hirshfeld-I, ISA, and IBO NACs for molecules, 1 but the DDEC6 average charge transfer magnitude more closely resembled that of the electrostatic potential tting group, IBO, and i-ACP than the average charge transfer magnitudes of MBIS, Hirshfeld-I, and ISA.
The high conuence ranking of DDEC6 cannot be solely attributed to either the presence of other stockholder partitioning methods in the dataset nor to the presence of electrostatic potential tting methods in the dataset. Consider a pared down dataset in which all stockholder partitioning methods except DDEC6 and all electrostatic potential tting methods are removed so that only ADCH, APT, CM5, DDEC6, EEQ, IBO, MBSMulliken, QTAIM, and VDD remain. As shown in Table 6, DDEC6 remains the top-ranked method in this pared down dataset.

An unambiguous scientic test of atomic population analysis methods
Confusion on whether it is possible to apply the scientic method to quantify properties of atoms in materials pertains to the issue of whether atom-in-material properties can be experimentally measured. While it is generally believed that NACs are not directly measurable experimentally, the situation is actually two-fold. For the vast majority of materials NACs are not directly measurable experimentally, but a few carefully chosen materials Table 5 Rank of each charge assignment method according to its amount of correlation to other charge assignment methods. The S a and U(a, f) ranking criteria always give the same order of methods. This provide clear enough experimentally measured atomic population data for falsiable scientic tests. It is obvious that atom-in-material properties for a completely isolated atom are experimentally measurable. For example, the NAC of a completely isolated Na + ion could be denitively measured in an experiment to be +1. However, this is not helpful, because all atomic population analysis methods would yield the correct NAC in this case. The challenge is to come up with more interesting cases where the experimental result is unambiguous and some population analysis methods fail unambiguously.
Here, I show that such situations do indeed occur. In other words, I show it is possible to unambiguously falsify some atomic population analysis methods using the scientic method. By unambiguously, I mean the conclusion is independent of opinions, interpretations, and perspectives.
As an example, consider the endohedral N@C 60 system in which a N atom sits inside a C 60 cage. Electron paramagnetic resonance (EPR) and electron nuclear double resonance (ENDOR) experiments showed the ground spin state is S ¼ 3/2. 70 (The ground spin state of an isolated N atom is also S ¼ 3/2.) These spectra also show the N atom occupies a central position and interacts only weakly with the C 60 cage. [70][71][72][73] ". from the missing nuclear quadrupole interaction a symmetric oncentre equilibrium position of the nitrogen atom can be deduced, implying an isotropic g-matrix." 74 The interaction between the enclosed N atom and C 60 cage is sufficiently weak that at room temperature the cage spins freely around the enclosed N atom leading to a spherically symmetric environment observed in the EPR and ENDOR experiments. 73 How much spin density is transferred between the enclosed N atom and the C 60 cage? ". because of the undetectable 13 C hyperne interaction, the admixture of fullerene molecular orbitals to the central atom wavefunction seems to be extremely small and, as a result, spin rotational interaction can also be neglected. (A 13 C hyperne interaction of the order of 0.05 mT corresponding to approximately 1.5 MHz is expected for a unit spin density on the C 60 shell. The observed 50 kHz linewidth therefore puts an upper limit of 3% to the transferred spin density.)" 74 In other words, the amount of spin transferred from the enclosed N atom to the C 60 cage is small or negligible. How much net charge is transferred from the enclosed N atom to the C 60 cage? "The UV/vis spectrum of N@C 60 is indistinguishable within experimental error from that of C 60 , conrming negligible coupling between nitrogen in its atomic ground state and C 60 cage molecular wave functions." 75 If the C 60 cage in N@C 60 carried a substantial net charge, this would have altered its UV-vis spectrum compared to isolated C 60 . Because the UV-vis spectrum was unaltered, net charge transfer from the enclosed N atom to the C 60 cage is negligible or small in magnitude.
The [Eu@C 60 ] + system exhibits remarkably different behavior than N@C 60 . First, the Eu atom in Eu@C 60 is markedly off-center. 76 Second, there is strong interaction between the Eu atom and the C 60 cage. In contrast to the UV-vis spectrum of N@C 60 which was equivalent to the isolated C 60 spectrum, the Eu@C 60 UV-vis spectrum shows dramatic differences. 77 Comparing the Eu L III -edge XANES spectra of Eu@C 60 to reference compounds showed the Eu atom in Eu@C 60 is in the +II oxidation state. 77 This implies the seven 4f electrons comprising a half-lled subshell remain on the Eu atom, 77 along with potentially part of the 6s electrons. The 4f electrons have a smaller average radius and are more tightly bound than the 6s electrons. "The [isolated] C 60 host has only deeply held paired electrons. 78 (Experiments show C 60 has a rst ionization energy of 6.4-7.9 eV, an electron affinity of approx. 2.6-2.8 eV, and a rst optical transition of approx. 3.2 eV. 79-82 )" 15 Therefore, electrons may be transferred from the Eu atom to the C 60 cage, but would not be transferred from the C 60 cage to the Eu atom. Together, these results show the Eu atom in [Eu@C 60 ] + should have a NAC between approximately 1 and 2 and an ASM between approximately 7 and 8. Table 7 summarizes computed NACs for 20 methods having a complete basis set limit. ASMs are also listed for those methods that compute them. These calculations used the PBE/ def2TZVPP optimized geometries and wavefunctions computed in GAUSSIAN 16. 64 The same soware programs were used to compute the NACs of these systems as were used for the water molecule in Section 3.1. An extremely ne (0.04 bohr) grid was used for the QTAIM method. Default settings were used for all other methods. The Multiwfn defaults for CHELPG, MK, and RESP used vdW radii of 1.5 A for C, 1.5 (MK and RESP) or 1.7 (CHELPG) for N, and 1.4554 for Eu. As recommended in the paper introducing the RESP method, a hyperbolic penalty Table 6 Rankings of nine charge assignment methods in a pared down dataset. The S a and U(a, f) ranking criteria always give the same order of methods function was used with two-stage tting and constants of a ¼ 0.0005 (stage 1), a ¼ 0.001 (stage 2 on selected atoms), and b ¼ 0.1 (both stages). 33 For comparison, Table 7 also shows a onestage RESP tting using the strong constraint (a ¼ 0.001, b ¼ 0.1) on all atoms. Several observations are: (1) Because the nuclear charge of N is +7, its maximum possible NAC of +7 would be achieved if all electrons were removed from this atom. The two-stage RESP NAC of 9.116 for the N atom clearly shows this method assigns a negative number of electrons (i.e., À2.116 electrons) to this atom. The same problem occurred for the HLY and MK analysis of N in N@C 60 . Because the number of electrons cannot properly be negative, these methods are falsied for the N@C 60 system. The one-stage RESP NACs using the strong constraint gave a NAC of 6.553 for the N atom which is much too high even though it is slightly below the atomic number of 7 for N.
(2) The ISA method gave a NAC of À3.082 for the N atom in N@C 60 , which is much too large in magnitude. Therefore, ISA is falsied for this material.
(  60 ] + leaves 9 (valence electrons for neutral Eu) -2.691 ¼ 6.309 valence electrons which are too few to explain the QTAIM ASM of 6.932 for Eu in this material. (If all of these remaining valence electrons were spin polarized they would produce an ASM of 6.309.) Hence, this QTAIM NAC is a bit too high in magnitude.
These results show some atomic population analysis methods are falsied for these materials using the scientic method. This does not necessarily imply those particular methods will not work for other materials, but it indicates those methods may not be reliable across diverse material types.
The observant reader will notice N@C 60 contains a 'buried' nitrogen atom. For comparison, the water molecule studied in Table 2 does not contain any buried atoms. A buried atom is any atom whose shortest distance to the material's van der Waals surface exceeds that atom's van der Waals radius. Materials with buried atoms are plentiful: all liquids, all solids (except one-and two-atom thick materials), and some gasses and plasmas contain buried atoms. Some molecules containing ve or more atoms have buried atoms. As indicated in Table 1 and described in prior literature, the CHELPG, HLY, ISA, and MK methods fail for many materials with buried atoms. 23,33,52 The RESP method was developed with the intention to x this problem, 33 but results for N@C 60 presented here show the RESP method is not reliable for xing this problem in some materials. Changing the form or strength of the RESP constraints could potentially address this problem, but this example clearly demonstrates the extreme challenge associated with trying to nd a RESP constraint that works well across diverse materials. Notably, it is not as easy as just making the constraints stronger or weaker, because a RESP constraint that is too strong for one material (or for one part of a material) may be too weak for another material (or for a different part of the same material).
Although the N@C 60 material contains a buried atom, the presence or absence of buried atoms played no role in the decision to select this material as a benchmark system. N@C 60 was chosen as a benchmark material, because to the best of the author's knowledge published experimental spectroscopic results have characterized its net atomic charges and atomic spin moments more accurately and denitely than for any other known material containing unpaired electron spins and at least two different atom types. As described earlier in this section, these experimental data show unambiguously that there is small or negligible charge and spin transfer from the N atom to the C 60 cage and the system's ground state is a spin quartet.

Seven confluence principles
The word conuence means a coming together, joining, or merging. In the statistical context of this paper, conuence denotes a joining together or merging of statistical characteristics. Two statistical characteristics that are normally thought to be distinct may actually merge to become a single characteristic. Also, various physical or statistical properties may be simultaneously highly correlated to a single quantitative descriptor. An analogy is useful. As illustrated in Fig. 3, consider a group of darts aimed at some target. The dart located in the center of the group never lands the farthest from any conceivable target. This centrally located dart exhibits conuence properties including high correlation to the other individual darts and to the main principal component of the dart group. If the group of darts follows a spherically symmetric distribution, then a centrally located dart lands closer to the target than at least $50% of the darts. In other words, the centrally located dart performs average or better for diverse targets. Other individual darts may land closer to the bullseye for specic targets, but the centrally located dart is best positioned for general-purpose use across diverse targets.
Conuence is the missing link that shows how to dene quantitative descriptors that are not directly experimentally observable (at least in most cases) to achieve high correlations to a host of related physical properties. In this article, we consider the task of assigning properties to atoms in materials. Atoms are the conceptual foundation of chemistry; however, many properties of individual atoms in materials are not directly observable experimentally for most materials. For example, the partial charge (i.e., NAC) of an atom in a material is not a direct experimental observable for most materials. Nevertheless, the concept of charged atoms (i.e., anions and cations) has been crucial to understanding the chemistry of many materials. By using conuence, a NAC descriptor can be constructed that exhibits good correlations to a host of chemical properties related to the partial charges of atoms in materials.
The remainder of this section precisely denes conuence and seven associated conuence principles.

Denition
A quantitative descriptor is dened as conuent among a group of positively correlated quantitative descriptors if this quantitative descriptor has sufficiently high correlation to the group's average standardized variable f. The precise threshold for "sufficiently high" must be (arbitrarily) chosen. Example: as shown in Table 4, the correlation U(DDEC6, f) ¼ 0.985 can be considered "sufficiently high" to label DDEC6 as a conuent descriptor for NAC methods.

Conuence principle #1
For a group of positively correlated quantitative descriptors, the descriptor with the highest correlation to the group's average standardized variable (U(a, f)) also has the highest sum of correlations to the individual group members (S a ). Proof: eqn (31) shows S a ¼ S f U(a, f). Because S f is the same for all group members, the group member with highest U(a, f) also has highest S a . Example: as shown in Table 4, the highest values correspond to S DDEC6 ¼ 18.204 and U(DDEC6, f) ¼ 0.985, which are related by S a /U(a, f) ¼ S f ¼ 18.4806. Implication: the centrally located dart exhibits not only the strongest correlation to the group's average position, but also the highest sum of correlations to all positions of the individual darts in the group.

Conuence principle #2
PCA of the correlation matrix for a group of quantitative descriptors yields coefficients for the k th principal component (PCk) that are directly proportional to each descriptor's correlation to this PC. Proof: the correlation between standardized variableâ and the k th principal component directly expands to give where v a (k) is the coefficient forâ in the k th eigenvector of the correlation matrix, l (k) is the corresponding eigenvalue, and s PC ðkÞ ¼ ffiffiffiffiffiffiffi l ðkÞ p (eqn (21)) is the standard deviation of PCk across the datapoints. Example: as shown in Table 3 Table 4. Implication: aer performing PCA of the correlation matrix, the ranking of variables according to their coefficients in the MPC is identical to the ranking of variables according to their correlation to the MPC.

Conuence principle #3
For a group of positively correlated quantitative descriptors, a quantitative descriptor's correlation to the group's average standardized variable f is similar (though not necessarily equal) to the same descriptor's correlation to the correlation MPC. Proof: see Section 7.3. Example: in Table 4, the largest difference magnitude between a single descriptor variable's correlation to f and MPC is 0.007. Implication: ranking variables according to (i) U(a, f) (equivalent to S a ranking) or (ii) U(a, MPC) (equivalent to MPC coefficient ranking) yields similar (not necessarily equal) results. Fig. 3 In a group of darts aimed at any target, the centrally located dart never lands farthest from the target. If the group of darts follows a spherically symmetric distribution, then a centrally located dart lands closer to the target than at least $50% of the darts. In other words, the centrally located dart performs average or better for diverse targets. This centrally located dart exhibits confluence properties including high correlation to the other individual darts and to the main principal component of the dart group.

Conuence principle #4
Among a group of positively correlated quantitative descriptors, the quantitative descriptor exhibiting conuence to the group's average standardized variable f has predictive advantages across a broad range of target applications. Explanation: here, the term "predictive advantages" refers to the fact that a centrally located dart will not land farthest from any related target (see Fig. 3). If the darts are approximately uniformly distributed over a spherical region, then the center dart lands closer to any target than at least $50% of the darts. This analogy extends to quantitative descriptors where a centrally located descriptor is a descriptor that is highly correlated to f. Example: as an example, DDEC6 NACs (which have high correlation to f) give good performance across both chemical properties and electrostatic properties of molecules.

Conuence principle #5
If a group of positively correlated quantitative descriptors contains two conuent descriptors a and b, then descriptors a and b are somewhat highly correlated to each other. Proof: using standardized variablesâ,b, andf, the correlations are proportional to the dot products over the sample data points: where dot product has the following denition Because the variables are standardized, In other words, a and b must both be approximately parallel to f, which can only occur if they are also approximately parallel to each other. This therefore implies that $(â,b) z M, and thus that U ab z 1. Example: comparing Fig. 2 to Table 4, the 15 descriptors having correlation > 0.9 to the DDEC6 NACs (the most conuent descriptor among the 20 NAC methods) were exactly the same 15 descriptors having highest correlation to f. The Spearman rank correlation between correlation to DDEC6 and correlation to f was 0.90 across the 20 methods, which reveals similar (but not completely identical) rankings according to correlation to DDEC6 NACs and correlation to f. Implication: this principle shows arbitrariness in designing descriptors is dramatically reduced when those descriptors are designed to be conuent. Specically, two different descriptors, each designed to be conuent across the same descriptor group, will be highly correlated to each other and thus not arbitrarily valued.

Conuence principle #6
If quantitative descriptor A is optimized to be conuent among a group of target physical properties, this same quantitative descriptor is expected to be conuent among a group of quantitative descriptors that are individually highly correlated to individual physical properties in this group. Explanation: suppose there are a group of physical properties designated P1, P2, P3, etc. that are experimentally measured across a sample population. Suppose further the quantitative descriptor A has been optimized to give high positive correlations between descriptor A and each individual physical property P1, P2, P3, etc. across this sample population. In other words, the correlation between descriptor A and property P1 is high across this sample population. The correlation between descriptor A and property P2 is also high across this sample population, and so forth for properties P3, etc. Now suppose there is another quantitative descriptor B1 that is optimized to give high positive correlation to physical property P1 across this sample population, but not necessarily high correlation between B1 and physical property P2 or P3 across this sample population. Now suppose there is another quantitative descriptor B2 that is optimized to give high positive correlation to physical property P2 across this sample population, but not necessarily high correlation between B2 and physical property P1 or P3 across this sample population. Now suppose there is another quantitative descriptor B3 that is optimized to give high positive correlation to physical property P3 across this sample population, but not necessarily high correlation between B3 and physical property P1 or P2 across this sample population. Likewise descriptors B4, B5, etc. are highly correlated to physical properties P4, P5, etc., respectively. Since descriptor A is conuent among related physical properties P1, P2, P3, etc., then it will also be conuent among a group of descriptors B1, B2, B3, etc. that are optimized to be highly correlated to physical properties P1, P2, P3, etc. Proof: because the standard deviation of any standardized variable across the sample population equals one, high positive correlations between descriptor A and properties P1, P2, etc. can only occur if for the vast majority of data points in the sample, where the overbar represents the average across the sample population and A i , P1 i , etc. represent the descriptor and property values for the i th datapoint in the sample population. Since descriptor B1 is highly positively correlated to property P1, it follows that for the vast majority of data points in the sample. Similarly, a high positive correlation between descriptor B2 and property P2 means that for the vast majority of data points in the sample. Combining eqn (41)- (43) gives for typical data points in the sample, which completes the proof. Implication: some NAC methods (e.g., ACP, ADCH, CM5, i-ACP) were optimized to reproduce molecular dipole moments.
Others were optimized to reproduce the electrostatic potential surrounding the molecule (e.g., CHELPG, HLY, MK, RESP). Others were optimized to maximize the similarity to quantummechanically computed reference atom densities (e.g., Hirshfeld, Hirshfeld-I) or orbitals (e.g., IBO). Others were optimized to reproduce constrained (e.g., MBIS) or unconstrained (e.g., ISA) spherically averaged AIM distributions. Others were optimized to reproduce electronegativity trends (e.g., EEQ) or number of electrons in the volume dominated by each atom (e.g., QTAIM). There are two different approaches to achieve high correlations across the majority of these descriptors. In approach II, we optimize a quantitative descriptor to be strongly correlated to a collection of many various related quantitative descriptors. In other words, we could optimize a NAC method to give NACs that are strongly correlated to the NACs produced by many various NAC methods. In approach I, we optimize a quantitative descriptor to strongly correlate to many various physical properties (MEP, molecular dipole moments, element electronegativities, etc.). For example, optimizing NACs to reproduce a variety of physical and chemical properties. Regardless of whether approach I or approach II is chosen, the end result is similar: the resulting descriptor will be conuent across this descriptor group and the related physical and chemical properties. Example: the DDEC6 NACs were designed to be conuent across various physical and chemical properties; they were not developed with the goal of giving strong correlations to other NAC assignment methods. 15 Nevertheless, they consequently developed strong correlations to other NAC assignment methods as demonstrated by the data in Table 4 and Fig. 2.

Conuence principle #7
The This means that both f and MPC are maximally correlated to the individual group members, and therefore likely strongly correlated to each other. Example: in agreement with the f and MPC optimization criteria, Table 8 shows f exhibits higher summed correlation compared to MPC, while MPC exhibits higher summed squared correlation compared to f. Because of conuence, the differences are tiny. Expanding the correlation between f and MPC gives  (45) gives U(f, MPC) ¼ 0.99994, which clearly indicates an almost perfect correlation between f and MPC for the NAC methods. Clearly, U(f, MPC) must be less than one, but it is extremely close to one for the group of NAC methods.

5.
How these confluence principles work together with other key principles and the scientific method to make assigning atom-in-material properties non-arbitrary NACs have no complete basis set limit. Consequently, a concept emerged in the early days that NACs are an extremely ill-determined 'arbitrary' concept. This idea of NAC arbitrariness was further encouraged by many poorly performing methods introduced in subsequent decades, oen without a clear understanding of the limitations of various approaches. For example, the Bickelhaupt, MBSBickelhaupt, Stout-Politzer, and Ros-Schuit methods lack rotational invariance; they produce different results when the entire molecule is rotated with respect to the coordinate axes. Because scalar properties like NACs are not vectors or tensors, their values should be independent of coordinate system orientation. Therefore, NAC methods lacking rotational invariance are unphysical. The Lowdin method, which was not included in Cho et al.'s dataset, was another early method that exhibits strong basis set dependence. While the original Lowdin method lacks rotational invariance, 83 the subsequent Davidson-Lowdin 84 method is rotationally invariant but still lacks a complete basis set limit. Some methods also had convergence problems: either failing to converge in some cases or converging to non-unique solutions. Some good methods were also introduced along the way. Methods for assigning atom-in-material properties should preferably work across an extremely wide range of material types. Arguments for a specialized atomic population analysis method that is specically optimized for a narrow material class are intrinsically weak. An atomic population analysis method that is specically optimized to describe one material class (and not other material classes) will be unable to describe systems containing that one material class together with other material types. For example, one person may claim to have developed a new atomic population analysis method that is specically optimized to describe ionic liquids and not other materials, while another person may claim to have developed a different atomic population analysis method that is specically optimized to describe metal-organic frameworks and not other materials. Neither of these methods are capable of describing the behavior of ionic liquids in metal organic frameworks, because they fail to simultaneously describe both material classes. As a second example, even though there are many different kinds of molecules, a charge assignment method that only works for molecules is quite limited, because it cannot even describe systems in which molecules react on solid surfaces. Since the number of possible chemical combinations is innite, this requires a generalpurpose atomic population analysis method that applies across an extremely wide range of material types.
While there exists some exibility in constructing atomic population analysis methods, this exibility should be constrained in several key ways:

Criterion 1
Atom-in-material properties should be mathematically well-dened with a complete basis set limit and rotational invariance. 10

Criterion 2
The method for computing atom-in-material descriptors should be physically well-motivated and derivable from fundamental principles. 10,56,85,86 Criterion 3 If the value of an atom-in-material descriptor corresponds to functional minimization, this functional should be convex to ensure the minimum is unique. 15,53,85 Moreover, nearly at optimization landscapes should be avoided. 33,53 Criterion 4 For the reasons discussed above, methods for assigning atomin-material properties should preferably work across an extremely wide range of material types containing both surface and buried atoms.

Criterion 5
For carefully selected benchmark systems, the computed atomin-material descriptor value should approximately match known reference values. For NACs and ASMs, the N@C 60 system discussed above is one such example. Examples of known bond orders include the H 2 (BO ¼ 1), N 2 (BO ¼ 3), and O 2 (BO ¼ 2) molecules.

Criterion 6
The method for computationally assigning atom-in-material properties should be compatible and consistent across various quantum chemistry methods. For example, it should give consistent results for different basis set types (e.g., plane waves, Gaussian, etc.) as well as for methods having idempotent (e.g., DFT, HF) and non-idempotent (e.g., CCSD, CAS-SCF, SAC-CI, etc.) rst-order density matrices. 10 One way to achieve this is to make the assigned atom-in-material descriptors functionals of the electron density and spin magnetization density distributions. 56

Criterion 7
Each computed atom-in-material descriptor should be designed to achieve conuence across related properties. In other words, it should be strongly correlated to many related experimentally measured and theoretically computed physical and chemical properties.

Criterion 8
An atomic population analysis method should preferably be capable of computing a whole suite of atom-in-material properties (net atomic charges, atomic spin moments, bond orders, spdfg populations, etc.) as opposed to assigning only one atomin-material property.

Criterion 9
The assigned values of atom-in-material properties should be chemically consistent. For example, the number of electrons assigned to an atom should be non-negative. Also, various atom-in-material descriptor values (e.g., NACs, ASMs, bond orders, spdfg populations) should be non-contradictory (i.e., approximately consistent with each other). For example, a hydrogen atom should not be assigned an ASM of 0.9 and a NAC of 0.75, because the former requires at least 0.9 electrons to reside on this atom while the latter requires 0.25 electrons to reside on this atom.

Criterion 10
The assigned atom-in-material descriptors should have good transferability between similar chemical environments. 87,88 Conformational transferability is especially important when parameterizing exible force elds. 27,53,89 Criterion 11 An atomic population analysis method should have reasonable computational costs. 10 (Note: the prior literature contains detailed computational cost studies for a small number of individual atomic population analysis methods, but no study has been published to date that systematically compares computational costs across a wide range of different atomic population analysis methods. 52,55,90-93 )

Criterion 12
The atomic population analysis method should not require the manual adjustment of computational parameters for individual systems; it should work out of the box without requiring system-specic tweaking from a human. 85

Criterion 13
The assigned electron density partitions "should be localized around the atomic nucleus and should not have intricate structures far from their dening nuclear center. This requirement is usual necessary, albeit insufficient, for chemical transferability and conformational stability." 85

Criterion 14
The assigned atom-in-material properties should properly reect the material's symmetry. 10

Criterion 15
For atom-in-material descriptors in which the sum over atoms has a well-dened value, this value should be properly reproduced. For example, the NACs should sum to the unit cell's net charge, 10 for collinear magnetism the ASMs should sum to the number of spin-up minus spin-down electrons in the unit cell, etc. Also, the local values of the electron density partitions should add up to the total electron density at each position in space (eqn (37)).
How does conuence specically relate to assigning atom-inmaterial charges? NACs could be optimized to reproduce the electrostatic potential surrounding a molecule (e.g., CHELPG, HLY, MK methods, etc.), to reproduce the molecular dipole moment (e.g. ADCH, etc.), to reproduce dipole moment derivatives (i.e., APT), to correspond to virial compartments (i.e., QTAIM), to match deformation densities (i.e., Hirshfeld, VDD), to project onto a basis of atomic orbitals (e.g., IBO, MBSMulliken, etc.), to maximize similarity between reference ions and assigned atom-in-material electron distributions (e.g. Hirshfeld-I, etc.), or to satisfy other criteria. But do these optimization criteria require different NAC methods? Could a NAC method be developed that simultaneously provides reasonably good correlations to most of these criteria?
Conuence is the concept of a centrally located method that is ideally positioned to give good correlations to a broad range of related physical properties. In the analogy of a group of darts aimed at targets, the central dart never lands farthest from any target. Conuence is a viable approach to constructing a truly general-purpose atomic population analysis method.
Conuence removes much of the arbitrariness associated with constructing an atomic population analysis method. It may appear arbitrary whether the NACs should be optimized to reproduce (a) the MEP, (b) the molecular dipole moment, or (c) the number of electrons in the local volume dominated by each atom-in-material, etc. However, much of this arbitrariness can be removed by optimizing NACs to achieve conuence across these various physical properties. We may conceivably construct at least two different atomic population analysis methods K and L which are each conuent across these target physical properties. According to conuence principle #5, the results of atomic population analysis methods K and L will be positively correlated to each other. Because each of the target properties (a) to (c) listed above are linear in the charge transfer magnitude, methods K and L must have similar charge transfer magnitudes to be conuent across these properties. Hence, the results of methods K and L must be approximately similar.
Manz and co-workers developed an extremely wide-ranging suite of atomic population analysis tools called the Standard Atoms in Materials Framework (SAMF): (i) ASMs for materials with collinear and non-collinear magnetism, 86 (ii) bond orders, 56 (iii) orbital bond order components that sum to the correct bond orders, 94 (iv) atom-in-material polarizabilities, dispersion coefficients, and quantum Drude oscillator parameters, 95,96 (v) various generations of charge partitioning schemes, 15,52-54 (vi) many linear-scaling computational algorithms, 55,56,86,96 and (vii) a complete library of chargecompensated reference ions for all charge states of chemical elements 1 to 109. 15,52,53,97 This charge-compensated reference ion library and methods to compute ASMs, bond orders, bond order components, polarizabilities, dispersion coefficients, and quantum Drude oscillator parameters can in principle be used with multiple charge assignment methods. However, consistently high accuracy is obtained only when using a high-quality and extremely versatile charge partitioning method such as DDEC6 or similar. 15,54,56,95 When these techniques are used with DDEC6 or potentially similar high-quality partitioning, all 15 criteria described above can be satised.
The DDEC6 method is strongly based on conuence. DDEC6 atomic population analysis was developed to achieve conuence across various observable chemical properties rather than to maximize its correlation to other atomic population analysis methods. The DDEC6 NACs are simultaneously optimized to give strong correlations to: (i) the electrostatic potential surrounding the material, (ii) the number of electrons in the local volume dominated by each atom-in-material, (iii) dipole moments, (iv) element electronegativities, and other properties. 15 The DDEC ASMs are simultaneously optimized to resemble proportional spin partitioning and spherical averaging of the spin magnetization density vectors. 55,86 The Manz bond orders, which oen use DDEC6 charge partitioning, are based on the conuence of atom-in-material exchange propensities. 56 The MCLF atom-in-material polarizabilities and dispersion coefficients (which oen use DDEC6 charge partitioning) are based on m-scaling, conduction limit upper bound, and other scaling principles to achieve accurate results for both surface and buried atoms in diverse materials. 95,96 According to conuence principle #6, this will naturally result in strong correlations between DDEC6 NACs and NACs computed by other methods. This is clearly demonstrated by the data in Table 4 and Fig. 2. According to conuence principle #4, this gives DDEC6 atomic population analysis predictive advantages across a wide range of target applications.
In summary, there is some exibility in designing atomic population analysis approaches, but various approaches that produce a chemical descriptor strongly correlated to many related physical properties must also produce strong correlations in-between these different atomic population analysis approaches. In other words, there may be several paths to achieve similar descriptor values. This is why methods like DDEC6 and IBO that are based on entirely different approaches yield similar average charge transfer magnitudes and highly correlated NACs for small molecules. It is not the values themselves of atom-in-material properties, but various paths to achieve similar values that embodies most of the exibility of constructing general-purpose atomic population analysis methods. Incorporating diverse material classes (e.g., molecules, dense solids, porous solids, conductors, insulators, magnetic materials, multi-reference systems, etc.) and computational approaches (e.g., DFT and various correlated wavefunction methods) can reveal which strategies are broadly applicable and which perform well only for specic material types. Methods that perform well only for limited material types should be replaced by more broadly applicable methods.

Conclusions
Assigning properties to atoms in materials is not arbitrary. It is theoretically possible to develop a method to assign NACs that simultaneously has average or better correlations to any and all physical and chemical properties related to atom-in-material charges. In other words, it is theoretically possible to develop a universally good method to assign NACs and other atom-inmaterial properties. This can theoretically be achieved by centrally locating the atomic population analysis method so that it exhibits strong correlations to other atomic population analysis methods. Among existing atomic population analysis methods, the DDEC6 method currently comes closest to this ideal.
Linear least-squares tting is an extremely common technique. However, simple least squares tting (SLSF) is not reversible: a SLSF of variable y to x yields a linear model that is not mathematically equivalent to a SLSF of variable x to y. 42 A bivariate standardized reversible linear least squares tting was introduced here that solves this problem and has four important properties: (i) it is a total least squares t with Euclidean metric, (ii) it is an orthogonal distance regression, (iii) it is a PCA regression, and (iv) it has a universal model equation that requires no computerized calculations. Because of property (iv), it is called instant least squares tting (ILSF). The ILSF universal linear model equation can be applied to any pair of positively correlated quantitative descriptors; however, it will achieve the best results when those two descriptors are approximately linearly correlated to each other.
As an example, this ILSF was used to compute average charge transfer magnitudes of 26 different methods to assign NACs across $2000 molecular systems. The Hirshfeld and VDD methods (which partition deformation densities) had the smallest average charge transfer magnitudes, while QTAIM (which assigns virial compartments) had the largest average charge transfer magnitude. NACs (e.g., ACP, ADCH, CM5, i-ACP) intended to reproduce the molecule's dipole moment had smaller average charge transfer magnitudes than those NACs optimized to reproduce the electrostatic potential surrounding the molecule (e.g., CHELPG, HLY, MK, RESP). The Bickelhaupt, DDEC6, IBO, and ISA methods gave average charge transfer magnitudes similar to the electrostatic potential tting group.
The correlation matrix between 20 NAC methods having complete basis set limits was extensively analyzed. This correlation matrix had a main block comprised of 14 NAC methods with strong inter-correlations plus two small side blocks weakly connected to the main block. Principal components analysis showed the main (or rst) principal component accounts for 17.16 variables' worth (85.8%) of the correlation in this group. Each of the other principal components accounted for less than one variable's worth of correlation apiece. The NAC methods were ranked according to how strongly correlated they are to all 20 NAC methods. The top (DDEC6), the bottom (Becke), the 2 nd (MBIS), the 8 th (RESP), the 9 th (MK), the 11 th (CM5), and the 15 th (Hirshfeld) ranked methods had consistent rankings across various ranking criteria. The DDEC6 method exhibited correlation >0.9 to 15 of 20 methods, and it had a summed correlation ¼ 18.204. The Becke method exhibited R < 0.7 to all other NAC methods. The DDEC6 NACs had correlation R ¼ 0.986 and 0.985 to the MPC and average standardized variable f, respectively. The MBIS, ISA, and Hirshfeld-I NACs also exhibited high correlations to these descriptors and other NAC methods.
Calculations in Section 3.3 showed the top ranking is relatively stable to the choice of which other charge assignment methods are included in the dataset. For example, the top-ranked method was unchanged when the number of different charge assignment methods was increased to 26 or decreased to 9.
Although NACs are not unambiguously measurable experimentally for most materials, N@C 60 is an interesting benchmark case for which experimental spectroscopy results show negligible or small charge and spin transfer from the N atom to the C 60 cage. Calculations in Section 3.4 for N@C 60 showed that some charge assignment methods give unphysical results for this material while other charge assignment methods performed well. This example demonstrates that it is possible, at least in some cases, to falsiably test atomic population analysis methods using the scientic method.
Seven conuence principles were derived that explain many connections between correlation properties. For example, NAC methods ranked similarly according to the sum of correlations to other methods (S a ), correlation (U(a, f)) to the average standardized variable f, coefficient in the correlation main principal component (MPC), correlation to this MPC, and number of NAC methods to which a NAC method is strongly correlated (e.g., U ab > 0.9). Many relations between these correlation properties were derived and proved.
A quantitative descriptor with strong correlations to many related descriptors has predictive advantages across multiple applications. This can be illustrated via an analogy to a group of darts aimed at a target. A dart near the center of the group lands closer to the bullseye than at least $50% of the darts in all cases, irrespective of the target. This conuence should be used to construct general-purpose atomic population analysis methods. A general-purpose atomic population analysis method should have NACs that are conuent across properties related to atomic charges, bond orders that are conuent across properties related to bond orders, ASMs that are conuent across properties related to the atom-in-material ordering of unpaired electron spins, and so forth.
In addition to achieving conuence, a general-purpose atomic population analysis method should also satisfy many other criteria as described in Section 5. For example, it should give chemically accurate results across diverse material types, give consistent results across different quantum chemistry methods (e.g., various basis sets and exchange-correlation theories), be capable of computing many different atom-in-material descriptors that are approximately chemically consistent with each other, have approximate transferability of descriptor values between similar chemical environments, be computationally efficient, not require manual tweaking for individual materials, have good convergence properties, and so forth.
Finally, the correlations reported in this article for $2000 main group molecules (which contain many surface atoms) should not be extrapolated to dense solids (which contain many buried atoms). It oen occurs that two charge assignment methods give similar results for surface atoms but vastly different results for buried atoms. 33,52 Future studies should consider more diverse material types. The most conuent method identied in this study (i.e., DDEC6) was previously shown to perform well across an extremely broad range of material types: small organic and inorganic molecules, dense solids, porous solids, nanostructured materials, large biomolecules, ionic liquids, polymers, organometallic complexes, heterogenous and homogenous catalysts, metal-organic frameworks, conductors, semi-conductors, and insulators, etc. 15,[54][55][56]98,99 Moreover, DDEC6 yields a wide range of atom-in-material properties: bond orders, 56 net atomic charges and atomic multipoles, 15,54 atomic spin moments, 55,86 polarizabilities and dispersion coefficients and quantum Drude oscillator parameters (when used in conjunction with the MCLF method 95,96 ), electron cloud parameters, 15,100 and bond order components. 94 DDEC6 has been used to construct exible force elds for various materials. [101][102][103][104][105][106][107][108][109][110][111][112] 7. Appendix: mathematical proofs 7.1 Proof that total least squares and orthogonal distance regression yield the same reversible linear model Consider a linear model of the form Because the error measure in eqn (15) is reversible, the model's predicted values are related to the measured values via the equations Using the error measure of approach 1, which rearranges to give The minimum occurs when First, simplies to 0 ¼ 2(1 + z À2 )N(zw avg + h À z avg ) By denition (see eqn (6) and (7)), z avg ¼ w avg ¼ 0. Because w and z are real-valued and positively correlated, z is real-valued and non-zero. Also, N $ 2. Accordingly, eqn (54) simplies to h ¼ 0.
Putting h ¼ 0 into eqn (51) and simplifying gives where the sums have been evaluated in terms of the correlation matrix elements (eqn (2) and (5)). Second, which simplies to which has only two real-valued roots: z ¼ À1 and z ¼ 1. The condition z(1 + z À2 ) ¼ U wz would also yield vL/vz ¼ 0, but cannot be satised for any real value of z for 0 < U wz # 1. Inserting these into eqn (55) yields z ¼ 1 is the global minimum solution, because w and z are positively correlated by construction (i.e., U wz > 0). Next, I show the same solution results from orthogonal distance regression of the standardized variables (approach 2). As shown in Fig. 1, the perpendicular distance is related to distances Dw and Dz that were considered in the total least squares tting. Specically, the area of the blue triangle in Fig. 1 is given by area ¼ 1 2 base Â height ¼ 1 2 (Dw)(Dz) ¼ 1 2 ht. Hence, Substituting eqn (49) into (60) and simplifying gives The orthogonal regression minimizes the sum of squared error (SSE) Hence has the same solution h ¼ 0 as above. Expanding reveals vL (2) /vz has exactly the same z ¼ À1 and z ¼ 1 roots as vL (1) / vz. For these two roots, combining eqn (58), (59) and (62) yield Hence, approach 1 (total least squares with Euclidean metric) and approach 2 (orthogonal distance regression) of the standardized variables produce exactly the same global minimum solution The quality of the t can be quantied by À1 the sum of squared perpendicular errors divided by the sum of squared deviations of one variable from its average value: Hence, the t quality equals the correlation. Because the variables are standardized, the same result occurs if z is used in place of w in the denominator of eqn (68).

Proof that f maximizes summed correlations to the variables {â}
I now prove that f is the descriptor that maximizes S s for any conceivable descriptor s that is a linear combination of the standardized variables in a positively correlated descriptor group: The standard deviation is The objective function to be maximized expands as (71) This is maximized when vS s vK ðbÞ ¼ 0 Differentiating eqn (71) Eqn (31) shows U(a, f) ¼ S a /S f . Comparing eqn (77) to eqn (31) proves U(a, MPC) is approximately proportional to U(a, f). Furthermore, the power-law method to determine the largest eigenvalue (see eqn (32)) implies that Using p ¼ 2 andñ trial ¼1 (i.e., a vector lled with ones), gives Inserting eqn (79) into eqn (77) and comparing to eqn (31) gives the nal result U(a, MPC) z U(a, f)

Proof the MPC maximizes the sum of squared correlations to individual members of a descriptor group
This section proves the following: the MPC of the correlation matrix is the solution to a conuent optimization, where the MPC is a normalized linear combination of members of a descriptor group: (a) the MPC maximizes variance across the dataset and (b) the MPC maximizes the sum of squared correlations to individual members of the descriptor group. Either criterion (a) or (b) could be enforced leading to identical MPC. Enforcing the normalization constraint, the MPC variance is Following criterion (a), the variance is maximized by which is manifestly the eigenstate equation dening correlation MPC. The quantity maximized by criterion (b) is The derivative expands as is an eigenvector of U gb , this derivative simplies to vG vn MPC This maximizes G due to the derivative being zero. Therefore, criterion (b) has the same solution as criterion (a).

Conflicts of interest
There are no conicts of interest to declare.