Introducing DDEC6 atomic population analysis: part 2. Computed results for a wide range of periodic and nonperiodic materials

Abstract

Net atomic charges (NACs) are widely used throughout the chemical sciences to concisely summarize key information about charge transfer between atoms in materials. The vast majority of NAC definitions proposed to date are unsuitable for describing the wide range of material types encountered across the chemical sciences. In this article, we show the DDEC6 method reproduces important chemical, theoretical, and experimental properties across an extremely broad range of material types including small and large molecules, organometallics, nanoclusters, porous solids, nonporous solids, and solid surfaces. Some important comparisons we make are: (a) correlations between various NAC models and spectroscopically measured core-electron binding energy shifts for Ti-, Fe-, and Mo-containing solids, (b) comparisons between DDEC6 and experimentally extracted NACs for formamide and natrolite, (c) comparisons of accuracy of different NAC methods for reproducing the electrostatic potential surrounding a material across one and multiple system conformations, (d) comparisons between calculated and chemically expected electron transfer trends for atoms in numerous dense solids, solid surfaces, and molecules, (e) an assessment of NAC transferability between three crystal phases of the diisopropylammonium bromide organic ferroelectric, and (f) comparisons between DDEC6 and polarized neutron diffraction atomic spin moments for the Mn₁₂-acetate single-molecule magnet. We find the DDEC6 NACs are ideally suited for constructing flexible force-fields and give reasonable agreement with force-fields commonly used to simulate biomolecules and water. We find the DDEC6 method is more accurate than the DDEC3 method for analyzing a broad range of materials. This broad applicability to periodic and non-periodic materials irrespective of the basis set type makes the DDEC6 method suited for use as a default atomic population analysis method in quantum chemistry programs.

---

1. Introduction

The concept of atoms in materials is one of the most universal concepts in the chemical sciences. For example, a water molecule is comprised of two hydrogen atoms and one oxygen atom. Materials contain a cloud of electrons surrounding atomic nuclei. Conceptually partitioning a material into constituent atoms is equivalent to assigning fractions of the electron cloud to each atomic nucleus. The net atomic charge (NAC) assigned to each atom in the material is defined by

q_A = z_A − N_A (1)

where z_A is the nuclear charge (i.e., element number) of atom A, N_A is the number of electrons assigned to atom A, and q_A is the NAC of atom A.

There is some flexibility when assigning quantitative properties to the individual atoms in materials. Different definitions for assigning electrons to each atom in a material lead to different NAC values. For example, a method that assigns 8.7 (0.65) electrons to each oxygen (hydrogen) atom in water yields NACs of −0.7 (+0.35) for each oxygen (hydrogen) atom, while a method that assigns 8.8 (0.60) electrons to each oxygen (hydrogen) atom in water yields NACs of −0.8 (+0.4) for each oxygen (hydrogen) atom.

Currently, there is a pressing need to develop an atomic population analysis method suitable for use as a default method in quantum chemistry programs. Because quantum chemistry programs are used to study a broad range of material types, such a method should have extremely broad applicability. Because different quantum chemistry programs use different basis set types (e.g., plane-waves, Gaussian basis functions, Slater-type orbitals, numerical basis sets, etc.), it is preferable for the method to have no explicit basis set dependence. This will ensure the method converges towards a well-defined mathematical limit as the basis set is improved towards completeness. In contrast, the Mulliken¹ and Davidson–Löwdin² population analysis methods currently used as default methods in some quantum chemistry programs do not have any mathematical limits as the basis set is systematically improved towards completeness.³ Several charge partitioning methods with well-defined basis set limits have been developed, but these have other limitations as described in our prior article.⁴ To address these limitations, we introduced a new atomic population analysis method, called DDEC6, that is suitable for use as a default method in quantum chemistry programs.⁴

The theory and computational methods for the DDEC6 method were described in our prior article.⁴ The purpose of the present article is to test the performance of the DDEC6 method across a wider range of material types. A diverse materials set was carefully selected to evaluate the accuracy of our new charge partitioning method. To test whether the DDEC6 method consistently performs better than the DDEC3 method, we include many systems for which the DDEC3 method was originally tested.⁵ In addition, we study many new materials carefully selected for their ability to make falsifiable tests of a charge assignment method's ability to describe electron transfer. One of the most frequent concerns about charge assignment methods is that it is difficult to compare them directly to experimental data. Therefore, we included many materials having strong experimental data. These comparisons to experimental data allow our Results and discussion to be viewed not only as applications of the DDEC6 method but also as performance tests.

A few remarks are appropriate pertaining to the charge assignment methods against which the new DDEC6 charge partitioning is compared. Because there are so many different charge assignment methods, it was impractical to compare all charge assignment methods for each material studied here. Therefore, we adopted the policy to compare against an appropriate subset of charge assignment methods for each material. Since DDEC6 is the successor to DDEC3, we compared DDEC6 to DDEC3 in most cases. In most cases, we included the charge assignment methods one would expect to perform the best for each kind of material. For example, electrostatic potential fitting (ESP^6,7 or REPEAT⁸) NACs were included in most comparisons based on the electrostatic potential root-mean-squared error (RMSE) and relative root-mean-squared error (RRMSE). For dense solids, Bader's quantum chemical topology^9,10 was included, because it is currently the most widely used charge partitioning method for dense solids. We avoided Mulliken and Davidson–Löwdin charges, because these are extremely sensitive to the basis set choice.³ We included comparisons to the Hirshfeld¹¹ (HD) method in many cases, because it is easy to do even though the HD method usually underestimates NAC magnitudes.^5,12–14 We restricted comparisons to Iterative Hirshfeld¹³ (IH) and related methods to previously published results, because the several different variations of these methods and their various reference ion densities is beyond the scope of this article.^13,15–20 (In our prior article, we presented data for three systems proving for the first time that the IH optimization landscape is sometimes non-convex and converges to non-unique solutions.⁴) Because several of the systems studied here were suggested in an article by Wang et al.²¹ focusing on applications of CM5, we compared DDEC6 to CM5 results in those cases and a few others. For a few molecular systems, we also compared results to Natural Population Analysis³ (NPA) and Iterated Stockholder Atoms²² (ISA) charges. None of the dense materials included comparisons to the ISA charges, because these are known to perform poorly for dense solids.²³ We do not include comparisons to Atomic Polar Tensor²⁴ (APT) or Born effective²⁵ charges, because DDEC6 charges quantify a system's static electron distribution while APT and Born effective charges quantify the system's response to a perturbation.^4,24,25 As discussed in our prior article, APT and Born effective charges are not designed to assign core electrons to the host atom.⁴

2. Computational details

We performed periodic quantum chemistry calculations using VASP^26,27 software. Our VASP calculations used the projector augmented wave (PAW) method^28,29 to perform all-electron frozen-core calculations including scalar relativistic effects with a plane-wave basis set cutoff energy of 400 eV. For all systems, the number of k-points times the unit cell volume exceeded 4000 ų. This is enough k-points to converge relevant properties including geometries and atoms-in-materials (AIM) properties (NACs, ASMs, etc.). Except where otherwise specified, geometry optimizations relaxed both the unit cell vectors and ionic positions. The solid surface calculations (see Section S3 of the ESI†) used the Density Functional Theory (DFT) optimized bulk lattice vectors and relaxed the ionic positions. Where noted, experimental crystal structures or other geometries from the published literature were used. A Prec = accurate (∼0.14 bohr) electron density grid spacing was used. Bader NACs were computed using the program of Henkelman and coworkers.³⁰

We performed non-periodic quantum chemistry calculations using GAUSSIAN 09 (ref. 31) software. ESP NACs were computed in GAUSSIAN 09 using the Merz–Singh–Kollman scheme.^6,7

Computational details for the DDEC6 charge partitioning and for the electrostatic potential RMSE are described in our previous article.⁴ Our parallelized code for computing DDEC6 NACs, ASMs, and bond orders is currently available in the CHARGEMOL program distributed via http://ddec.sourceforge.net.³²

3. Results and discussion

3.1 Representing electron transfer between atoms

3.1.1 Metal oxides and sulfides. We now study electron transfer between atoms in the dense solids shown in Fig. 1. While we were testing modifications of the DDEC method for the sodium chloride crystals,⁴ Wang et al.²¹ pointed out a related problem with the DDEC3 method. Specifically, when DDEC3 NACs are compared for a series of transition metal oxide solids with and without Li atoms, the DDEC3 NACs on the transition metal atoms exhibit a trend that does not match chemical expectations. As shown in Table 1, the DDEC3 NAC on the Co atom is lower in crystalline CoO₂ than in LiCoO₂. In contrast, the Bader, CM5, and HD NACs on the Co atom are higher for crystalline CoO₂ than for LiCoO₂. To assess which trend is correct, Wang et al. plotted isosurfaces of the electron density difference between CoO₂ and LiCoO₂ using the M06L³³ functional and found a slight increase in electron density around the Co and O atoms upon Li addition to CoO₂ to create LiCoO₂.²¹ Thus, the Bader, CM5, and HD methods predict the correct charge transfer direction between these two materials, but the DDEC3 method predicts the wrong charge transfer direction between these two materials.²¹ For reasons clearly explained in prior publications, charge transfer magnitudes predicted by the HD method are usually much too small.^5,12–14 In Table 1, we compare NACs computed using the PBE³⁴ functional optimized geometries and electron distributions. For DDEC3, the previously reported M06L results are also listed for comparison.²¹ As shown in Table 1 (PBE results) and Wang et al.²¹ (M06L results), the CM5 and Bader methods predict a decrease of the transition metal NAC upon lithiation for the solids TiS₂ → LiTiS₂, LiTi₂O₄ → LiTiO₂, MnO₂ → LiMn₂O₄, while the DDEC3 method predicts an increase for all except LiTi₂O₄ → LiTiO₂. Assuming these materials behavior similar to the CoO₂ material, a decrease in the transition metal NAC upon lithiation is chemically expected. Thus, we employed these materials as a test set to evaluate the performance of potential modifications to the DDEC method when developing the DDEC6 method. In addition, we studied the Li₃RuO₂ crystal suggested to us by Ayorinde Hassan. Charge partitioning for the Li₃RuO₂ crystal is challenging due to the large proportion of Li atoms and the nearly neutral Ru atoms, because the neutral Li and Ru reference atoms are much more diffuse than the cationic ones leading to large sensitivity of the reference ion densities on the reference ion charges. As shown in Table 1, the DDEC6 algorithm yields reasonable NACs for all of these materials. Only for TiS₂ → LiTiS₂ is there a small increase from 1.32 to 1.38 in the transition metal DDEC6 NAC upon lithiation. For all of these materials, the Bader and DDEC6 methods gave similar Li NACs, while the CM5 and HD methods gave substantially smaller Li NACs.

Fig. 1 Unit cells used to model metal oxide and sulfide solids. The lines mark the unit cell boundaries. Atoms are colored by element: Li (green), O (red), S (yellow), Ti (light blue), Co (dark blue), Mn (magenta), Ru (beige).

Table 1 Average HD, CM5, DDEC3, DDEC6, and Bader charges of Li, transition metal (TM), and nonmetal atoms. NACs shown are for the PBE optimized geometries and electron densities

Crystal	HD			CM5			DDEC3			DDEC6			Bader
	Li	TM	Anion	Li	TM	Anion	Li	TM	Anion	Li	TM	Anion	Li	TM	Anion
a NACs from ref. 21 using M06L optimized geometries and electron distributions.b LiMn2O4 has a spinel structure that undergoes a charge-ordering transition as shown in experiments;^35–37 the PBE functional shows charge disproportionation between the Mn sites (i.e., a charge-ordered phase) while the M06L functional gives equal NACs on all Mn sites^21 (i.e., a high-temperature phase without charge ordering).
LiCoO2	0.11	0.34	−0.23	0.49	0.73	−0.61	1.03 (1.00^a)	1.47 (1.45^a)	−1.25 (−1.23^a)	0.87	1.07	−0.97	0.88	1.22	−1.05
CoO2	—	0.35	−0.18	—	0.80	−0.40	—	1.14 (1.23^a)	−0.57 (−0.62^a)	—	1.12	−0.56	—	1.39	−0.69
LiTiS2	0.07	0.40	−0.23	0.27	0.79	−0.53	0.98 (0.97^a)	1.67 (1.48^a)	−1.33 (−1.23^a)	0.86	1.38	−1.12	0.89	1.48	−1.18
TiS2	—	0.43	−0.21	—	0.86	−0.43	—	1.06 (1.06^a)	−0.53 (−0.53^a)	—	1.32	−0.66	—	1.61	−0.80
LiTiO2	0.11	0.56	−0.34	0.46	1.16	−0.81	1.05 (1.00^a)	2.17 (2.10^a)	−1.61 (−1.55^a)	0.89	1.65	−1.27	0.89	1.57	−1.23
LiTi2O4	0.16	0.64	−0.36	0.48	1.31	−0.78	1.03 (1.00^a)	2.33 (2.32^a)	−1.42 (−1.41^a)	0.90	1.94	−1.19	0.91	1.84	−1.15
LiMn2O4^b	0.17	0.34	−0.21	0.53	0.84	−0.55	0.99 (1.00^a)	1.56 (1.95^a)	−1.03 (−1.23^a)	0.86	1.23	−0.83	0.89	1.59	−1.02
MnO2	—	0.36	−0.18	—	0.88	−0.44	—	1.24 (1.47^a)	−0.62 (−0.73^a)	—	1.25	−0.63	—	1.69	−0.85
Li3RuO2	0.11	0.31	−0.32	0.33	0.58	−0.79	0.83	−0.18	−1.15	0.72	−0.08	−1.04	0.82	0.12	−1.30

---

3.1.2 Palladium-containing crystals. As additional examples of charge transfer in solids, we studied an interstitial H atom in Pd, Pd₃V, Pd₃In, and Pd₃Hf crystals, plus Pd₃V with no interstitial H atom. Manz and Sholl previously studied these materials with the DDEC2, DDEC3, and Bader methods,⁵ and we used their geometries and PW91 electron densities to now compute the HD, CM5, and DDEC6 NACs. (These geometries are representative local energy minima, not necessarily global energy minima.⁵) Interestingly, the HD and CM5 NACs are negative for V, In, and Hf atoms and positive for Pd atoms, even though the Pauling scale electronegativity³⁸ of Pd (2.20) is greater than V (1.63), In (1.78), and Hf (1.3). The Bader, DDEC3, and DDEC6 NACs followed the Pauling scale electronegativity trends with a negative average Pd NAC and positive X NACs following the expected trend Hf > V > In. All methods gave slightly negative to nearly neutral H NACs within the range −0.32 to +0.02 (Table 2).

Table 2 Average NACs for interstitial H in ordered Pd₃X alloys

Material	H charge^a					Pd charge^a					X charge^a
	Bader	DDEC3	DDEC6	HD	CM5	Bader	DDEC3	DDEC6	HD	CM5	Bader	DDEC3	DDEC6	HD	CM5
a Bader and DDEC3 NACs are from ref. 5.
H in Pd	−0.04	−0.25	−0.05	−0.13	−0.12	0.001	0.008	0.001	0.004	0.004	n.a.	n.a.	n.a.	n.a.	n.a.
Pd3V	n.a.	n.a.	n.a.	n.a.	n.a.	−0.35	−0.10	−0.15	0.02	0.02	1.04	0.31	0.44	−0.06	−0.06
H in Pd3V	−0.22	−0.32	−0.14	−0.15	−0.13	−0.34	−0.09	−0.14	0.03	0.02	1.04	0.32	0.44	−0.06	−0.06
H in Pd3In	−0.05	−0.18	−0.05	−0.11	−0.09	−0.21	−0.08	−0.07	0.06	0.16	0.64	0.27	0.22	−0.16	−0.47
H in Pd3Hf	−0.05	0.02	−0.01	−0.10	−0.09	−0.53	−0.31	−0.23	0.08	0.08	1.58	0.92	0.69	−0.22	−0.22

---

Why did the HD and CM5 methods yield negative NACs for the V, In, and Hf atoms? It is well-known that isolated neutral atoms usually become more contracted upon going from left to right within the same subshell of a periodic table row due to the increasing nuclear charge that contracts the subshell. (Deviations from this trend can occur where electron configurations deviate from the Aufbau principle, such as Pd through Cd.) Moreover, atoms usually become slightly more diffuse or remain about the same size down a periodic table column. Accordingly, an isolated neutral Hf atom is more diffuse than an isolated neutral Pd atom. The Pauling scale electronegativity will usually follow the opposite trend, with the Pauling scale electronegativity increasing left to right within the same subshell of a periodic table row and decreasing down a periodic table column except where electron configurations deviate from the Aufbau principle. Because the neutral Hf reference atom is more diffuse than the neutral Pd reference atom, during HD partitioning the Hf atoms steal electrons from the more electronegative Pd atoms. Thus, in this case, the HD method predicts the wrong charge transfer direction. The CM5 method adds a correction to the HD NACs, but this correction is zero between two transition metal atoms.¹² Consequently, the HD and CM5 NACs are identical for Pd₃V. In the other materials, there is a non-zero CM5 correction between the main-group elements H and In and the other elements, which causes the CM5 NACs to slightly differ from the HD NACs.

To avoid this problem, the DDEC3 and DDEC6 methods include a constraint that forces w_A(r_A) for tails of buried atoms to decay at least as fast as exp(−1.75r_A/bohr).⁵ Second, the DDEC6 method sets the reference ion charge for each atom in the material to a weighted average of a stockholder type charge partitioning and a smoothed localized charge partitioning.⁴ This ensures the reference ion charge resembles the charge in the local vicinity of the atom in order to prevent atoms from becoming too diffuse or too contracted. This makes DDEC6 NACs more accurately describe the true charge transfer direction.

3.1.3 Magnesium oxide. Table 3 compares six different charge assignment methods for (MgO)_n molecules (n = 1 to 6) and crystalline MgO. Geometries of the (MgO)_n molecules and their HD, CM5, and DDEC3 NACs and dipole moments were taken from Wang et al.²¹ These geometries (Fig. 2) were built by removing Mg and O atoms from a rigid (MgO)₆ cluster, rather than optimizing the geometries with DFT.²¹ Following Wang et al.,²¹ we computed electron distributions for the (MgO)_n clusters in GAUSSIAN 09 using the M06L functional and def2-TZVP³⁹ basis set. The geometry and electron distribution of crystalline MgO were optimized in VASP using the PBE functional.

Table 3 Comparison of different charge assignment methods for (MgO)_n molecules and crystalline MgO. For the DDEC6 method, M represents NACs only, D represents the inclusion of atomic dipoles, and SCP represents the inclusion of the spherical charge penetration term

Method	System^a
	MgO	(MgO)2	(MgO)3	(MgO)4	(MgO)5	(MgO)6	Bulk MgO
a NACs and dipole moments for the HD, CM5, and DDEC3 methods for the (MgO)n molecules are from ref. 21.b ESP NAC cannot be reported for bulk MgO, because there are no surface atoms.
NAC of central Mg
HD	0.59	0.69	0.56	0.42	0.28	0.15	0.33
CM5	0.79	0.97	0.93	0.88	0.82	0.77	0.77
DDEC6	0.99	1.45	1.57	1.60	1.61	1.61	1.47
Bader	1.18	1.54	1.65	1.62	1.57	1.48	1.70
DDEC3	1.00	1.55	1.72	1.76	1.79	1.84	2.01
ESP	0.89	1.16	1.04	0.96	0.58	−1.72	^b

---

Method		Molecule^a
		MgO	(MgO)2	(MgO)3	(MgO)4	(MgO)5	(MgO)6	MAE
Dipole moment in a.u.
Full density		2.71	2.32	1.39	1.70	0.48	1.97	0.00
HD		1.87	1.97	1.28	2.08	0.88	1.98	0.35
ESP		2.82	2.74	1.87	2.13	1.12	2.43	0.42
CM5		2.51	2.54	1.91	2.47	1.51	2.61	0.56
DDEC6	M	3.14	3.18	2.37	2.49	1.68	2.91	0.87
	D	2.71	2.32	1.39	1.70	0.48	1.97	0.00
DDEC3		3.17	3.31	2.48	2.57	1.79	3.08	0.97
Bader		3.71	3.39	2.46	3.09	1.63	3.77	1.25

---

Method		Molecule^a
		MgO	(MgO)2	(MgO)3	(MgO)4	(MgO)5	(MgO)6	Average RMSE
RMSE in kcal mol^−1 (RRMSE)
ESP		2.81 (0.10)	5.40 (0.23)	4.59 (0.24)	5.16 (0.27)	4.96 (0.27)	3.79 (0.21)	4.45
CM5		3.98 (0.15)	6.50 (0.28)	4.77 (0.25)	5.90 (0.31)	6.12 (0.33)	5.28 (0.29)	5.43
HD		9.27 (0.35)	10.05 (0.43)	6.06 (0.32)	5.99 (0.31)	4.90 (0.26)	4.01 (0.22)	6.71
DDEC6	M	4.25 (0.16)	7.24 (0.31)	6.75 (0.35)	7.70 (0.40)	7.77 (0.42)	6.58 (0.36)	6.71
	M + SCP	4.63 (0.17)	7.40 (0.32)	6.98 (0.36)	7.82 (0.41)	7.98 (0.43)	6.78 (0.37)	6.93
	D	1.31 (0.05)	2.60 (0.11)	2.87 (0.15)	2.86 (0.15)	3.09 (0.17)	3.03 (0.16)	2.63
	D + SCP	0.86 (0.03)	1.55 (0.07)	1.67 (0.09)	1.61 (0.08)	1.74 (0.09)	1.70 (0.09)	1.52
DDEC3		4.44 (0.17)	8.48 (0.36)	7.52 (0.39)	8.48 (0.44)	8.42 (0.45)	7.08 (0.39)	7.40
Bader		9.10 (0.34)	9.26 (0.40)	10.81 (0.56)	12.87 (0.67)	14.20 (0.76)	14.11 (0.77)	11.73

---

Fig. 2 Molecular electrostatic potential (MEP) of the six (MgO)_n molecules (n = 1 to 6) studied. The MEP is shown on the 0.0004 electrons per bohr³ density contour with a MEP scale ranging from −0.78 volts (red) to 0.78 volts (blue). The numbers appearing beside the terminal Mg atoms are their DDEC6 NACs.

Based on the much lower Pauling scale electronegativity of Mg (1.31) than O (3.44), a substantial transfer of electrons from Mg to O is expected. A simple chemical argument suggests that as the central Mg atom is surrounded by more oxygen anions, electrostatic stabilization of the central Mg cation by the oxygen anions should increase, thereby stabilizing more electron transfer from the central Mg atom to the adjacent oxygen atoms. This simple chemical argument predicts an increase in the central Mg atom NAC as the number of adjacent O atoms increases. Examining Table 3, only the DDEC3 method consistently followed this trend. The trend for terminal Mg NACs can be inferred from the electrostatic potential values. As shown in Fig. 2, the electrostatic potential and DDEC6 NACs are most positive near the terminal Mg atoms following the trend MgO > (MgO)₃ > (MgO)₄ > (MgO)₂ > (MgO)₅ > (MgO)₆.

As shown in Table 3, the Mg NAC in bulk MgO followed the trend DDEC3 (2.01) > Bader (1.70) > DDEC6 (1.47) > CM5 (0.77) > HD (0.33). The DDEC3 NAC of 2.01 for bulk MgO is similar to some recent high-resolution diffraction experiments and their interpretations in terms of fully ionized Mg²⁺ and O²⁻ ions.^40,41 However, the situation is not as straightforward as it first appears, because (i) charge partitioning in the experimentally measured electron distribution depends on model definitions used to assign NACs and (ii) the low-order structure factors in simple cubic crystals (e.g., MgO and NaCl) have low sensitivity to the amount of charge transfer.^21,40–42 Zuo et al. used a convergent beam electron diffraction technique to improve the resolution of the low-order structure factors and concluded the crystal's electron distribution is consistent with fully ionized Mg²⁺ and O²⁻ ions,⁴⁰ but this does not rule out other interpretations.

These tests on (MgO)_n and bulk MgO illustrate some possible compromises between matching the chemical state trends on the one hand and the electrostatic potential trends on the other hand. Dipole mean absolute error (MAE) followed the trend HD (0.35) < ESP (0.42) < CM5 (0.56) < DDEC6 (0.87) < DDEC3 (0.97) < Bader (1.25). Electrostatic potential RMSE (kcal mol⁻¹) followed the trend ESP (4.45) < CM5 (5.43) < HD, DDEC6 (6.71) < DDEC3 (7.40) < Bader (11.73). Although the ESP method gave low dipole MAE and electrostatic potential RMSE, we do not recommend the ESP method for assigning NACs, because the ESP NACs of the central Mg atom fluctuated erratically from 1.16 for (MgO)₂ to −1.72 for (MgO)₆. Because the Bader point charges had the highest dipole moment MAE and the highest average electrostatic potential RMSE, we do not recommend Bader NACs for use in force-field point charge models for classical atomistic simulations. Choosing between the remaining four point charge methods (i.e., HD, CM5, DDEC3, and DDEC6) is complicated by the fact that dipole MAE and electrostatic potential RMSE followed a different trend than the central Mg atom NAC. On the basis of the central Mg atom NAC increasing monotonically from MgO molecule to (MgO)₆ molecule to bulk MgO, the DDEC3 method would be preferable, but the DDEC3 method gave the largest dipole MAE and average electrostatic potential RMSE among these four charge assignment methods. The HD and CM5 methods had comparatively low dipole MAE and average electrostatic potential RMSE, but yielded low values of 0.33 (HD) and 0.77 (CM5) for the Mg NAC in bulk MgO. Results for the DDEC6 method were intermediate for central Mg NAC, dipole MAE, and average electrostatic potential RMSE.

Finally, Table 3 investigates effects of atomic dipoles and spherical charge penetration. Including atomic dipoles for any AIM method (e.g., HD, DDEC6, Bader, DDEC3), eliminates the dipole prediction error to within a grid integration tolerance (e.g., ∼0.01). Because the dipole moment of a spherical charge distribution is zero, the spherical charge penetration term has no effect on the computed dipoles. Including DDEC6 atomic dipoles decreased the average RMSE from 6.71 to 2.63 kcal mol⁻¹. Although the spherical charge penetration term slightly increased the average RMSE at the DDEC6 (M + SCP) level, it dramatically reduced the average RMSE to 1.52 kcal mol⁻¹ at the DDEC6 (D + SCP) level. Notably, the DDEC6 (D + SCP) average RMSE was ∼3 times lower than any of the point charge models.

3.2 Comparison to spectroscopic results for various materials

3.2.1 Net atomic charges extracted from high resolution diffraction experiments. Extracting NACs from high-resolution diffraction data is not straightforward, but it can be done using approximations and models. In ‘Kappa refinement’, the high-resolution diffraction data is first fit to a multipolar model^41,43,44 to determine atomic coordinates, thermal parameters, and an electron density map and then refit to a spherical pseudoatom model^45,46 to determine the NACs. In kappa refinement, the spherical pseudoatoms have the form ρ_at(r_A) = ρ^core_A(r_A) + n^val_Aκ³p^val_A(κr_A), where p^val_A(r_A) is the normalized shape function of the valence density of the neutral reference atom.⁴⁵ The two primary limitations of kappa refinement are that the pseudoatom densities do not necessarily sum to the correct total density ρ(r⃑) and the shape functions for the charged atoms are represented as expanded or contracted versions of the neutral atoms.⁵ Here, we revisit two examples for which DDEC3 and experimentally extracted NACs were compared in ref. 5: the formamide and natrolite structures shown Fig. 3. We refer the reader to the earlier publications for a discussion of the experimental details and analysis.^5,46–48 The same geometries and electron distributions are used in this work as in ref. 5.

Fig. 3 Formamide and natrolite structures. The lines in natrolite indicate the unit cell boundaries. Figure reproduced with permission from ref. 5. © American Chemical Society 2012.

As shown in Table 4, both the DDEC3 and DDEC6 NACs follow a trend similar to the experimentally extracted NACs for natrolite, except the DDEC3 NACs on all atoms except Na are significantly higher in magnitude than the experimentally extracted ones. For all atoms, the DDEC6 NACs are slightly lower in magnitude than the DDEC3 NACs, leading to an overall better agreement between the DDEC6 and experimentally extracted NACs. Only for the Na atom, which was fixed to a value of 1.00 in the experimental analysis,⁴⁸ is the experimentally extracted NAC closer to the DDEC3 value than the DDEC6 value.

Table 4 Experimental and theoretical natrolite NACs. DDEC3 and DDEC6 results computed using the PBE-optimized geometries

	High res. XRD^a	DDEC3^b	DDEC6
a High resolution XRD data from ref. 48.b DDEC3 NACs from ref. 5.
Si1	1.84 ± 0.12	2.172	1.772
Si2	1.65 ± 0.10	2.207	1.760
Al	1.51 ± 0.11	2.067	1.762
O1	−0.90 ± 0.05	−1.227	−1.036
O2	−1.21 ± 0.05	−1.318	−1.103
O3	−1.03 ± 0.05	−1.337	−1.094
O4	−1.07 ± 0.05	−1.320	−1.110
O5	−0.87 ± 0.05	−1.113	−0.913
Na	1.00	1.000	0.896
Ow	−0.59 ± 0.03	−0.926	−0.862
H1	0.24 ± 0.03	0.446	0.408
H2	0.36 ± 0.03	0.435	0.405

---

Table 5 summarizes experimentally extracted and computed NACs for formamide. Theoretical charges were computed using the B3LYP^49,50 functional with aug-cc-pvtz⁵¹ basis set. The high-resolution X-ray diffraction (XRD) results were extracted using fully optimized radial factors for all atoms.⁴⁶ Maximum absolute differences from the experimentally extracted NACs are 0.07 (NPA), 0.11 (DDEC3), 0.12 (DDEC6), 0.13 (IH), 0.17 (ESP), 0.22 (ISA), 0.64 (HD), and 0.96 (Bader). The DDEC6, ESP, DDEC3, and ISA point charge dipoles were within ±5% of the full wavefunction value of 1.55. Errors for the other point charge dipoles were −6% (IH), +23% (NPA), and +66% (Bader). Of course, when atomic dipoles are included, all of the AIM methods (Bader, DDEC3, DDEC6, HD, IH, and ISA) yield the exact dipole moment to the integration grid precision. The errors of the point charge models for reproducing the electrostatic potential followed the trend ESP < DDEC6 < DDEC3 < ISA < IH < NPA < HD < Bader. When atomic dipoles were included, the RMSE for the DDEC6 method decreased from 0.74 to 0.49 kcal mol⁻¹. When the spherical charge penetration term was included, the RMSE values for the DDEC6 method were unchanged (within a computational tolerance of 0.01 kcal mol⁻¹) at 0.74 (M + SCP) and 0.49 (D + SCP), indicating a negligible impact of spherical charge penetration over the RMSE grid points. In summary, the DDEC6 NACs did a good job of simultaneously reproducing the electrostatic potential, the molecular dipole moment, and the experimentally extracted NACs.

Table 5 Experimental and theoretical formamide NACs. Dipoles in atomic units. RMSE in kcal mol⁻¹

	High res. XRD^a	Bader^b	DDEC3^b	DDEC6 M (D)^d	ESP^b	HD^b	IH^b	ISA^b	NPA^b
a From ref. 46.b Bader, DDEC3, ESP, HD, IH, ISA, and NPA NACs are from ref. 5.c Dipole moment of the B3LYP/aug-cc-pvtz wavefunction was 1.55.d M denotes point charge (monopole) model; D denotes the inclusion of atomic dipoles.
O	−0.55 ± 0.04	−1.149	−0.557	−0.506	−0.562	−0.304	−0.537	−0.593	−0.605
N	−0.78 ± 0.07	−1.183	−0.788	−0.662	−0.923	−0.136	−0.862	−0.911	−0.808
C	0.51 ± 0.08	1.469	0.624	0.519	0.680	0.139	0.644	0.726	0.534
H1	0.39 ± 0.03	0.411	0.352	0.329	0.389	0.128	0.360	0.389	0.394
H2	0.40 ± 0.03	0.426	0.369	0.313	0.429	0.133	0.377	0.407	0.388
H3	0.03 ± 0.03	0.026	0.000	0.007	−0.012	0.040	0.018	−0.019	0.096
	Dipole moment^c	2.57	1.59	1.53 (1.55)	1.57	1.13	1.46	1.62	1.91
	RMSE	9.85	0.85	0.74 (0.49)	0.58	3.32	1.43	0.99	3.13
	RRMSE	0.89	0.08	0.07 (0.04)	0.05	0.30	0.13	0.09	0.28

---

3.2.2 Correlations between NACs and spectroscopically measured core electron binding energy shifts. The core electron binding energy shift is defined as the binding energy of a particular core orbital level for an atom in a material compared to the same core orbital level for an atom of the same element in a reference compound.^52–57 Core electron binding energies can be measured using X-ray photoelectron spectroscopy (XPS) or X-ray absorption near edge structure (XANES). Several key factors affect core electron binding energy shifts.^52–57 First, a change in the valence electron population of this atom affects its core electron binding energy, because more valence electrons cause electrostatic shielding of the nuclear charge and a decrease in the core electron binding energy.^52–56 Second, the core electron binding energy is directly affected by the electrostatic potential exerted on this atom by the other atoms in the material: lots of anions nearby will decrease the core electron binding energy and lots of cations nearby will increase the core electron binding energy.^52–56 Third, the core electron binding energy is affected by relaxation in which the electrons rearrange to partially fill the hole left by the ejected photoelectron.^53,55–57 Various simple model equations have been developed to correlate core electron binding energy shifts to easily computed chemical descriptors such as (a) the NAC of the atom emitting the photoelectron, (b) the electrostatic potential exerted on the atom emitting the photoelectron by all the other atoms in the material (as computed using electron distributions or simple point charge models), (c) quantum mechanically computed electrostatic potential near the nucleus of the atom emitting the photoelectron, (d) orbital eigenvalues (aka ‘orbital energies’) computed using the Hartree–Fock or other quantum chemistry methods for chemical models of the initial and final states, and (e) two-electron integrals describing exchange and electrostatic interactions between valence and core orbitals.^52–56

Here, we are most interested in correlations between core electron binding energy shifts and NACs that occur for some crystalline materials.^58–62 We now consider a series of Ti, Mo, and Fe compounds as examples. Table 6 summarizes linear correlations between core electron binding energies and NACs. The HD, DDEC6, and Bader methods gave reasonable performance (i.e., R-squared ≥ 0.704) for all three elements, while the DDEC3 and CM5 methods performed poorly (i.e., R-squared ≤ 0.360) for the Ti compounds. Overall, the strength of the correlation between NACs and core electron binding energies followed the trend HD > DDEC6 > Bader > DDEC3 > CM5. Table S1† summarizes details for the Ti-containing solids. NACs were computed using the PBE electron distributions for the experimental geometries defined by the Inorganic Crystal Structure Database (ICSD) codes in Table S1,† except for SrTiO₃ which was geometry optimized. The poor correlation of the DDEC3 method for the Ti-containing solids was primarily due to high NACs for TiO, TiN, and TiB₂ and a lower NAC for TiCl₄ than for TiO.

Table 6 R-Squared correlation coefficients between NACs and spectroscopically measured core electron binding energies. NAC methods ordered from highest to lowest average R-squared correlation coefficient

	Ti compounds	Mo compounds	Fe compounds
HD	0.795	0.987	0.819
DDEC6	0.704	0.978	0.868
Bader	0.727	0.911	0.817
DDEC3	0.360	0.977	0.905
CM5	0.345	0.898	0.747

---

Table S2† summarizes results for the Mo-containing solids. NACs were computed using the PBE electron distributions for the experimental geometries defined by the ICSD codes in Table S2.† For structures listing two ISCD codes, both crystal structures were included in the correlation to the experimental K-edge energy. Li et al. measured these K-edge energies using XANES.⁶² The K-edge energy is correlated to the binding energy of the K-shell (i.e., 1s) core electrons.^55,63 Our analysis and correlation for these Mo-containing solids is identical to that of Li et al.⁶² using DDEC3 and Bader NACs, except we have extended it to DDEC6, CM5, and HD NACs.

Fig. 4 shows linear regression plots between the average Fe DDEC6 NACs and the oxidation state (left panel) and the 2p_3/2 core electron binding energy (right panel). NACs were computed using the PBE electron distributions based on the following geometries: Fe (NAC is zero due to symmetry), Fe₂SiO₄ (PBE-optimized geometry of anti-ferromagnetic spinel phase⁵), Fe₂O₃ (PBE-optimized geometry of anti-ferromagnetic phase⁵), Fe₃O₄ (PBE-optimized geometry of anti-ferrimagnetic phase⁵), and Fe₃Si (experimental crystal structure⁶⁴). Our analysis for Fe-containing solids is similar to that of Manz and Sholl⁵ using DDEC3 NACs, except we have extended it to DDEC6, CM5, HD, and Bader NACs.

Fig. 4 Left: Correlation between Fe oxidation state and average Fe DDEC6 NAC for Fe, Fe₂SiO₄, Fe₂O₃, Fe₃O₄, and Fe₃Si solids. Right: Correlation between 2p_3/2 core electron binding energy (as measured using XPS) and average Fe DDEC6 NAC for these materials.

3.3 Reproducing the electrostatic potential in one system conformation

Because the assigned atomic electron distributions {ρ_A(r⃑_A)} exactly sum to the total electron distribution ρ(r⃑) at each position in space, all AIM methods yield a formally exact representation of the electrostatic potential in the form of a polyatomic multipole expansion with charge-penetration terms.^23,65,66 For conciseness, it is desirable to have this polyatomic multipole expansion converge rapidly with most of the electrostatic potential described by the leading-order terms. Many force-fields used in classical molecular dynamics and Monte Carlo simulations use point-charge models to estimate the electrostatic interaction energies between chemical species.^67–69 These types of force-fields can be parameterized using NACs and optionally atomic multipoles computed via quantum chemistry calculations. To be suitable for this purpose, we desire the DDEC6 NACs to approximately reproduce the electrostatic potential surrounding a material.

In this section, we compare the accuracy of the DDEC3 and DDEC6 NACs for reproducing the electrostatic potential surrounding a single geometric conformation. Table 7 lists 13 materials including small molecules and ion, a large biomolecule, several metal–organic frameworks, a nanosheet, and a nanotube. This represents several kinds of materials often encountered in classical molecular dynamics or Monte Carlo simulations.^67,68,70,71 For each material, the same electrostatic potential grid point files were used to compute the DDEC3 and DDEC6 RMSE values. The DDEC6 NACs reproduced the electrostatic potential better than the DDEC3 NACs in 8 of these systems. This shows the DDEC6 NACs are a slight improvement compared to the DDEC3 NACs for reproducing the electrostatic potential surrounding a material. Including DDEC6 atomic dipoles improved the RSME by > 0.4 kcal mol⁻¹ for the BN nanotube, lp-MIL-53(Al), IRMOF-1 (XRD and DFT geometries), Zn-nicotinate, and H₂PO₄⁻. This shows that overall including atomic dipoles produces a modest improvement in the RMSE accuracy. Adding spherical charge penetration at the monopole or dipole levels (i.e., M + SCP and D + SCP) had no significant effect for these materials.

Table 7 Accuracy of fitting the electrostatic potential. (Values in parentheses include spherical cloud penetration.) The best values for a point charge model are shown in boldface type. Values at M + SCP, D, or D + SCP are shown in boldface type if they are equal to or better than the best point charge model value

Material	Geom	XC	Basis set	RMSE (kcal mol^−1)				RRMSE
				DDEC3^a		DDEC6		DDEC3^a		DDEC6
				M	D	M (M + SCP)	D (D + SCP)	M	D	M (M + SCP)	D (D + SCP)
a DDEC3 data (except BN sheet, formamide, Zn-nicotinate, water, H2PO4^− and DNA) is from ref. 5.b From ref. 23.c From ref. 5.d From ref. 72.e From ref. 73.f From ref. 74.
B4N4	DFT^b	PW91	6-311+G*	0.26	0.33	0.17 (0.19)	0.35 (0.35)	0.08	0.10	0.05 (0.05)	0.10 (0.10)
BN tube	DFT^b	PW91	Planewave	8.81	2.40	5.91 (5.94)	1.11 (1.09)	2.13	0.58	1.43 (1.44)	0.27 (0.26)
BN sheet	DFT	PBE	Planewave	0.07	0.07	0.07 (0.07)	0.07 (0.07)	0.64	0.64	0.64 (0.61)	0.64 (0.61)
Formamide	DFT^c	B3LYP	Aug-cc-pvtz	0.85	0.40	0.74 (0.74)	0.49 (0.49)	0.08	0.04	0.07 (0.07)	0.04 (0.04)
lp-MIL-53(Al)	XRD^d	PW91	Planewave	1.57	0.59	1.46 (1.47)	0.60 (0.58)	0.80	0.30	0.74 (0.75)	0.30 (0.30)
IRMOF-1	XRD^e	PW91	Planewave	0.83	0.44	0.86 (0.86)	0.26 (0.24)	0.39	0.20	0.40 (0.40)	0.12 (0.11)
IRMOF-1	DFT^f	PW91	Planewave	0.65	0.58	0.82 (0.81)	0.28 (0.27)	0.27	0.24	0.33 (0.33)	0.12 (0.11)
ZIF-8	DFT^f	PW91	Planewave	0.88	0.72	0.85 (0.81)	0.79 (0.76)	0.57	0.47	0.56 (0.53)	0.52 (0.50)
ZIF-90	DFT^f	PW91	Planewave	0.81	0.84	1.03 (0.97)	0.93 (0.88)	0.12	0.12	0.15 (0.14)	0.14 (0.13)
Zn-nicotinate	DFT	PBE	Planewave	0.82	0.44	0.90 (0.89)	0.41 (0.40)	0.46	0.25	0.51 (0.51)	0.23 (0.23)
Water	DFT	B3LYP	6-311++G**	1.31	0.80	1.16 (1.16)	0.88 (0.88)	0.14	0.08	0.12 (0.12)	0.09 (0.09)
H2PO4^−	DFT	M06L	Aug-cc-pvtz	2.16	0.41	1.65 (1.65)	0.49 (0.49)	0.17	0.03	0.13 (0.13)	0.04 (0.04)
DNA	DFT	PBE	Planewave	13.91	13.79	12.67 (12.68)	12.77 (12.77)	0.59	0.58	0.54 (0.54)	0.54 (0.54)

---

Water is the most abundant solvent in biology and chemical processing. Because water is vital to life on earth, it plays a key role in nearly all health applications. Water also plays a key role in environmental, weather, and climate change processes. Consequently, water is the most important molecule for molecular modeling in general. Because many classical atomistic molecular dynamics and Monte Carlo simulations will use DDEC6 NACs for non-water molecules combined with a well-established commonly used water model for the water solvent, it is desirable for the DDEC6 NACs for the water molecule to be approximately consistent with those of commonly used water models. Lee et al. computed NACs for large unit cells of simulated bulk water (∼2500 atoms with PBE functional and large psinc basis sets) and showed the DDEC/cc2 (q_H = 0.3915, q_O = −0.783) and DDEC3 (q_H = 0.402, q_O = −0.804) results are similar to common 3-site water models.^79,80 Table 8 lists commonly used 3-site water models that have been optimized to reproduce various properties of bulk water in classical atomistic simulations.^75–78 For comparison, DDEC6 results for the isolated water molecule with B3LYP/6-311++G** optimized geometry and electron distribution are q_H = 0.3953, q_O = −0.7906. A recent study by Farmahini et al. computed DDEC3 NACs to study changes in the hydrophobicity/hydrophilicity of nanoporous silicon carbide-derived carbon upon fluorine doping.⁸¹ These results show the DDEC methods are well-suited for studying water molecules.

Table 8 Commonly used 3-site water models listed in alphabetical order

Model	O–H distance (Å)	H–O–H angle (°)	H NAC	O NAC
a From ref. 75.b From ref. 76.c From ref. 77.d From ref. 78.
SPC^a	1.00	109.47	0.41	−0.82
SPC/E^b	1.00	109.47	0.4238	−0.8476
TIP3P^c	0.9572	104.52	0.417	−0.834
TIPS^d	0.9572	104.52	0.40	−0.80

---

The B-DNA decamer (CCATTAATGG)₂ structure was obtained from a neutron diffraction experiment performed by Arai et al. (PDB ID 1WQZ).⁸⁴ The 25H₂O molecules in the crystal structure are from the solvent and are hydrogen-bonded to the B-DNA as shown in Fig. 5. We added a Na⁺ ion next to each phosphate group, following previous studies to simulate the B-DNA being in a real solution.⁸⁵ Using the PBE functional, we optimized the positions of the Na⁺ ions in VASP while keeping the experimental B-DNA structure fixed. The left panel of Fig. 5 shows the optimized B-DNA decamer including the Na⁺ ions and hydrogen-bonded water molecules. Because the DDEC6 NACs can be used as input for classical molecular dynamics or Monte Carlo simulations of biomolecules, it is instructive to compare the DDEC6 NACs to those of common biomolecular force-fields.^79,80 The right panel of Fig. 5 compares the DDEC6 NACs to the CHARMM27 and AMBER4.1 forcefield NACs for DNA. There is some scatter in the data, but the overall correlation between DDEC6 and force-field NACs is good with R-squared correlation coefficients of 0.93 (CHARMM27) and 0.91 (AMBER4.1). The phosphorus NAC was 1.5 (ref. 82) for CHARMM compared to 1.166 (ref. 83) for AMBER, with the DDEC6 value of 1.38 in-between. Moreover, our DDEC6 NACs for this B-DNA decamer are similar to our previous DDEC3 results.⁸⁶ Finally, we note that recent articles by Lee et al. studied applications of DDEC NACs to atomistic simulations of large biomolecules including a comparison of force-fields based on AMBER and DDEC NACs for several large proteins.^79,80

Fig. 5 Left: B-DNA decamer (CCATTAATGG)₂, the lines mark the unit cell boundaries. Right: Correlation between DDEC6, CHARMM, and AMBER force-field NACs for all atoms excluding the bound water molecules and added Na⁺ atoms. The black line has a slope of 1 and intercept of 0. CHARMM27 NACs from Foloppe and Mackerell.⁸² AMBER4.1 NACs from Cornell et al.⁸³

3.4 Reproducing the electrostatic potential across multiple system conformations for constructing flexible force-fields

3.4.1 Carboxylic acids. In a previous publication, Manz and Sholl studied the accuracy of HD, DDEC3, ISA, IH, NPA, and ESP NACs for reproducing the electrostatic potential across various conformations of the five 4-X-substituted bicyclo[2,2,2]octane-1-carboxylic acids shown in Fig. 6.⁵ They found the ESP NACs reproduce the electrostatic potential as accurately as possible when optimized individually for each conformation, but have low conformational transferability.⁵ When using a conformationally averaged set of NACs to reproduce the electrostatic potential across the various conformations of each molecule, the ESP NACs performed slightly better than the DDEC3 NACs.⁵ When using NACs from the low energy conformation to reproduce the electrostatic potential across the various conformations of each molecule, the DDEC3 NACs performed slightly better than the ESP NACs.⁵ The DDEC3 NACs also had excellent conformational transferability.⁵

Fig. 6 Structures and molecular electrostatic potentials (MEPs) of the low energy conformations (B3LYP/6-311++G** level of theory) of 4-X-substituted bicyclo[2,2,2]octane-1-carboxylic acids: X = (a) –H, (b) –OH, (c) –CO₂C₂H₅, (d) –Br, and (e) –CN. The MEP is shown on the 0.0004 electrons per bohr³ density contour with a MEP scale ranging from −1.6 volts (red) to 1.6 volts (blue). The electrostatic potential is negative near the oxygen, bromine, and nitrogen atoms and positive near the proton of the carboxylate group.

We now show these desirable properties are further improved by the DDEC6 NACs. The B3LYP/6-311++G** optimized geometries and electron distributions from ref. 5 are used. Table 9 summarizes the fragment charges for each of these charge assignment methods, where the weighted sum as defined by Manz and Sholl⁵ is:

where q_A is the NAC for atom A and N_bonds is the number of bonds in the shortest chain connecting the atom to the substituent group. The purpose of this weighted sum is to smooth out the effects of the NACs, where all of the atoms in the substituent group are weighted by q_A and those not in the substituent group receive a diminished weight that tends towards zero as the atom is far removed from the substituent group. Roberts and Moreland determined σ′ substituent constants using experimentally measured acid dissociation constants.⁸⁷ As shown in Table 9, the HD NACs were most closely correlated to the σ′ values, where the R-squared correlation coefficient is that for linear regression: q_frag = a₀ + a₁σ′.⁵ The DDEC6 NACs showed the second strongest correlation to the σ′ values, with an R-squared correlation coefficient of 0.90 for the weighted sum in eqn (2). This shows the DDEC6 NACs captured the important chemical trend among the substituent groups.

Table 9 Fragment charges for the low energy conformation. NAC methods ordered from highest to lowest R-squared correlation coefficient for weighted sum

X	Substituent net charge^b								Weighted sum of eqn (2)^b
	σ′^a	HD	DDEC6	DDEC3	ISA	IH	NPA	ESP	HD	DDEC6	DDEC3	ISA	IH	NPA	ESP
a From ref. 87.b DDEC3, ESP, HD, IH, ISA, and NPA NACs are from ref. 5.
H	0.000	0.03	0.05	0.04	−0.01	0.07	0.21	−0.02	−0.01	−0.04	−0.03	0.02	−0.06	−0.20	−0.01
OH	0.283	−0.08	−0.22	−0.25	−0.31	−0.25	−0.29	−0.33	−0.05	−0.12	−0.12	−0.09	−0.15	−0.34	−0.10
CO2C2H5	0.297	−0.03	−0.02	−0.02	−0.07	0.06	0.02	−0.04	−0.04	−0.10	−0.08	−0.04	−0.09	−0.29	−0.02
Br	0.454	−0.1	−0.23	−0.25	−0.29	−0.11	−0.02	−0.19	−0.07	−0.17	−0.16	−0.15	−0.13	−0.30	−0.18
CN	0.579	−0.16	−0.17	−0.16	−0.25	−0.12	−0.02	−0.31	−0.11	−0.16	−0.14	−0.10	−0.16	−0.32	−0.10
R^2 corr. coef.		0.92	0.54	0.44	0.51	0.26	0.17	0.44	0.93	0.90	0.81	0.71	0.69	0.59	0.47

---

We now consider accuracy of these charge assignment methods for reproducing the electrostatic potential across various system conformations. As shown in Table 10, the conformational transferability of the charge assignment methods from best to worst ordered HD > NPA > DDEC6 > IH > DDEC3 > ISA > ESP. Reasonable conformational transferability of IH charges and poor conformational transferability of ISA charges have also been shown in prior work.^5,88,89 Table 11 compares the electrostatic potential RMSE and RRMSE values averaged across all molecular conformations for each of the point charge models using (a) NACs optimized individually for each conformation, (b) conformation averaged NACs, and (c) the low energy conformation NACs. DDEC6 values at the M + SCP (individually optimized for each conformation and using the low energy conformation), D (individually optimized for each conformation), and D + SCP (individually optimized for each conformation) levels are also shown for comparison. As expected, the ESP NACs reproduced the electrostatic potential most accurately among all point charge models optimized individually for each conformation. Although the DDEC6 NACs gave significantly higher RMSE and RRMSE values than the ESP NACs, the DDEC6 RMSE and RRMSE values including atomic dipoles (e.g., D and D + SCP) were approximately the same as those for ESP NACs optimized individually for each conformation. When using the conformation averaged NACs, ESP NACs still yielded the best overall results with the DDEC3 and DDEC6 NACs not far behind. When using NACs from the low energy conformation, the DDEC6 method provided the best overall results.

Table 10 Assessment of the conformational transferability of different charge assignment methods. NAC methods ordered from highest to lowest conformational transferability

Substituent	Mean unsigned deviation of NACs^a
	H	Br	CN	OH	Ester
a Mean unsigned deviations of NACs for the DDEC3, ESP, HD, IH, ISA, and NPA methods are from ref. 5.
Conformations	4	4	4	8	16
HD	0.002	0.002	0.002	0.003	0.002
NPA	0.004	0.004	0.004	0.006	0.006
DDEC6	0.005	0.005	0.005	0.007	0.006
IH	0.006	0.006	0.006	0.007	0.007
DDEC3	0.007	0.007	0.007	0.010	0.008
ISA	0.015	0.015	0.016	0.025	0.016
ESP	0.038	0.057	0.051	0.045	0.041

---

Table 11 Average RMSE and RRMSE values for charge assignment methods. NAC methods listed in alphabetical order. The best values for a point charge model are shown in boldface type. Values at M + SCP, D, or D + SCP are shown in boldface type if they are equal to or better than the best point charge model

Substituent		Avg. RMSE^a (kcal mol^−1)					Avg. RRMSE^a
		H	Br	CN	OH	Ester	H	Br	CN	OH	Ester
a RMSE and RRMSE for the DDEC3, ESP, HD, IH, ISA, and NPA methods are from ref. 5.
NACs optimized separately for each conformation
DDEC3		0.81	1.15	0.87	0.89	0.77	0.13	0.18	0.10	0.13	0.10
DDEC6	M	0.90	1.17	1.04	1.16	1.02	0.14	0.18	0.12	0.16	0.13
	M + SCP	0.89	1.16	1.03	1.15	1.01	0.14	0.18	0.11	0.16	0.13
	D	0.33	0.87	0.37	0.47	0.38	0.05	0.13	0.04	0.07	0.05
	D + SCP	0.32	0.88	0.37	0.48	0.38	0.05	0.13	0.04	0.07	0.05
ESP		0.49	0.93	0.38	0.48	0.42	0.07	0.14	0.04	0.07	0.04
HD		2.85	3.26	3.67	3.73	3.27	0.41	0.48	0.40	0.51	0.41
IH		1.12	2.49	1.60	1.35	1.05	0.18	0.38	0.18	0.20	0.14
ISA		0.73	1.45	0.74	0.71	0.70	0.11	0.22	0.08	0.10	0.09
NPA		1.71	3.23	2.56	1.94	2.98	0.25	0.49	0.29	0.27	0.31
Conformation averaged NACs
DDEC3		1.27	1.48	1.29	1.40	1.25	0.18	0.22	0.14	0.19	0.16
DDEC6		1.21	1.41	1.31	1.38	1.27	0.18	0.21	0.14	0.19	0.16
ESP		1.10	1.36	1.00	1.37	1.38	0.15	0.20	0.11	0.19	0.14
HD		2.88	3.31	3.70	3.71	3.31	0.41	0.50	0.41	0.50	0.42
IH		1.48	2.65	1.84	1.75	1.39	0.23	0.41	0.21	0.25	0.18
ISA		1.33	1.81	1.31	1.57	1.43	0.18	0.26	0.14	0.21	0.18
NPA		2.12	3.47	2.86	2.42	3.71	0.30	0.52	0.32	0.33	0.38
All conformations use NACs from low energy conformation
DDEC3		1.39	1.61	1.38	1.73	1.44	0.19	0.23	0.15	0.24	0.18
DDEC6	M	1.31	1.49	1.35	1.63	1.34	0.18	0.22	0.15	0.22	0.17
	M + SCP	1.30	1.49	1.44	1.62	1.33	0.18	0.22	0.15	0.22	0.17
ESP		1.49	1.73	1.26	2.12	1.91	0.19	0.24	0.13	0.28	0.20
HD		2.97	3.24	3.68	3.74	3.31	0.42	0.48	0.41	0.51	0.42
IH		1.55	2.61	1.85	1.98	1.48	0.23	0.40	0.21	0.28	0.19
ISA		1.46	1.97	1.41	1.98	1.74	0.19	0.29	0.15	0.27	0.21
NPA		2.23	3.50	2.93	2.57	4.05	0.31	0.52	0.32	0.35	0.42

---

Finally, we considered the 25 conformations of the –OH substituted carboxylic acid generated by the ab initio molecular dynamics (AIMD) calculations of Manz and Sholl at 300 K (Nosé thermostat).⁵ Following AIMD calculations in VASP using the PW91 functional with D2 dispersion corrections, Manz and Sholl computed the electron distributions and electrostatic potentials in GAUSSIAN 09 using the B3LYP/6-311++G** level of theory for each geometry.⁵ We use these same geometries, electron distributions, and electrostatic potentials here. Our purpose here is to see how the DDEC6 NACs perform compared to the previously reported results⁵ for the DDEC3, ESP, HD, IH, ISA, and NPA methods. As shown in Table 12, the DDEC6 NACs had lower electrostatic potential RMSE and RRMSE values across these AIMD conformations than any of the other six charge assignment methods when using either the low energy NACs or the conformation averaged NACs from Table 11. In summary, all of these tests for substituted carboxylic acids show the DDEC6 NACs have desirable properties for constructing flexible force-fields for classical atomistic simulations of materials: (a) reproduce chemical trends, (b) good conformational transferability, and (b) reasonable accuracy for reproducing the electrostatic potential across various system conformations.

Table 12 Average RMSE (kcal mol⁻¹) and RRMSE values for geometries of the –OH substituted carboxylic acid generated using ab initio molecular dynamics. NAC methods listed in alphabetical order. The best values are shown in boldface type

	DDEC3^a	DDEC6	ESP^a	HD^a	IH^a	ISA^a	NPA^a
a From ref. 5.
Using the low energy conformation NACs of Table 11
RMSE	2.51	2.17	3.23	4.38	2.55	3.04	3.65
RRMSE	0.27	0.23	0.35	0.47	0.27	0.33	0.39
Using the conformation averaged NACs of Table 11
RMSE	2.08	1.88	2.11	4.31	2.11	2.44	3.39
RRMSE	0.23	0.20	0.23	0.47	0.23	0.26	0.37

---

3.4.2 Li₂O molecule. Wang et al. compared several charge assignment methods for the Li₂O molecule constrained to bent angles of 90, 100, … 170° with bond lengths and electron distributions at each of these angles optimized using the M06L functional and def2-TZVP basis set.²¹ They found the CM5 NACs closely reproduced the Li₂O dipole moment while the DDEC3 NACs significantly overestimated the Li₂O dipole moment.²¹ In symmetric non-linear conformations of Li₂O, the NACs that exactly reproduce the molecular dipole moment are uniquely defined (aka ‘Dipole charge’).²¹ Here, we revisit this example to study in greater depth relationships between NACs, molecular dipole moments, electrostatic potential RMSE and RRMSE, and atomic dipole moments.

Table 13 summarizes computed NACs for each of the geometries studied by Wang et al.²¹ plus the global low energy conformation using the M06L/def2-TZVP level of theory. The global low energy conformation is a linear molecule corresponding to a 180° angle. The Dipole charge cannot be computed for this low energy conformation, because its molecular dipole is zero irrespective of the NAC. For all of the charge assignment methods except the Bader method, the NAC increased monotonically as the angle increased. (For the Bader method, the increase was almost monotonic.) In order of smallest to largest Li NACs, the charge assignment methods were HD < CM5 < Dipole charge < NPA, ESP, DDEC6, Bader < DDEC3.

Table 13 Li NAC for different Li–O–Li angles in singlet Li₂O molecules using various charge models. NAC methods listed in alphabetical order

Angle	Li net atomic charge^a									(Geom opt)
	90	100	110	120	130	140	150	160	170	180
a Except for the linear molecule NACs for the HD, CM5, DDEC3, ESP, NPA, and Dipole charge methods are from ref. 21.b Cannot be determined because the dipole moment is zero and the molecule is linear and symmetric.
Bader	0.83	0.85	0.85	0.86	0.86	0.86	0.86	0.87	0.86	0.87
CM5	0.55	0.56	0.57	0.58	0.59	0.60	0.60	0.61	0.61	0.62
DDEC3	0.83	0.88	0.92	0.94	0.96	0.97	0.98	0.98	0.98	0.98
DDEC6	0.77	0.80	0.83	0.85	0.87	0.88	0.89	0.90	0.90	0.90
Dipole charge	0.59	0.59	0.60	0.60	0.60	0.61	0.62	0.64	0.65	^b
ESP	0.61	0.63	0.65	0.67	0.70	0.74	0.79	0.83	0.86	0.87
HD	0.39	0.40	0.40	0.41	0.41	0.42	0.42	0.43	0.43	0.44
NPA	0.78	0.79	0.80	0.81	0.83	0.84	0.85	0.85	0.86	0.86

---

To further understand these trends, Table 14 summarizes the electrostatic potential RMSE and RRMSE, dipole moment MAE, and the mean unsigned deviation (MUD) from the conformation averaged NAC. From best to worst conformational transferability, the charge assignment methods ordered Bader > Dipole charge, HD > CM5, NPA > DDEC6, DDEC3 > ESP. From best to worst accuracy in reproducing the dipole moment, the methods ordered Dipole charge > CM5 > ESP > HD > NPA > DDEC6 > Bader > DDEC3. For the RMSE and RRMSE, HD performed the worst of all the charge assignment methods, and DDEC3 performed the second worst. Among the different point charge models, the ESP NACs provided the lowest RMSE and RRMSE when the NACs were optimized separately for each molecular conformation.

Table 14 Average electrostatic potential RMSE (kcal mol⁻¹), RRMSE, dipole moment MAE in atomic units, and conformational transferability of Li₂O for various charge assignment methods. Charge assignment methods listed in alphabetical order

		Conformation averaged NACs		NACs optimized separately for each conformation		NACs from the lowest energy conformation		Dipole moment MAE	NAC conformational transferability (MUD)
		RMSE	RRMSE	RMSE	RRMSE	RMSE	RRMSE
a Not computed.b Not computed, because the variation in the molecular conformation affects the orientation of the atomic dipoles.c Dipole charges cannot be determined, because the dipole moment of the lowest energy (i.e., linear) conformation is zero irrespective of the NAC values.d Since no Dipole charges were available for the linear molecule, these represent values for the nine non-linear conformations.
Bader		6.46	0.25	6.30	0.25	7.18	0.28	0.63	0.01
CM5		6.92	0.29	6.90	0.29	6.31	0.26	0.07	0.03
DDEC3		8.01	0.31	7.55	0.30	9.09	0.36	0.80	0.05
DDEC6	M	6.48	0.26	5.76	0.23	7.17	0.28	0.59	0.05
	M + SCP	^a	^a	5.85	0.23	7.19	0.28	0.59	0.05
	D	^b	^b	1.43	0.06	^b	^b	0.00	0.05
	D + SCP	^b	^b	1.20	0.05	^b	^b	0.00	0.05
Dipole charge^d		6.68	0.29	6.62	0.28	^c	^c	0.00	0.02
ESP		5.55	0.23	4.40	0.18	7.27	0.28	0.19	0.11
HD		11.40	0.48	11.49	0.48	10.72	0.44	0.50	0.02
NPA		6.10	0.24	5.56	0.22	7.03	0.28	0.53	0.03

---

Across all of the accuracy measures listed in Table 14, the following overall trends were observed: (a) the NPA, DDEC6, CM5, and Bader NACs performed better than the DDEC3 NACs, (b) the NPA NACs performed better than the DDEC6 NACs, (c) across the subset of accuracy measures where the Dipole charges were defined, they performed as good as or better than the CM5 and DDEC3 NACs, and (d) all other comparisons between NAC methods yielded mixed results, with better performance for at least one accuracy measure and worse performance for at least one accuracy measure.

For comparison, Table 14 also lists DDEC6 results including spherical charge penetration and atomic dipoles. Adding spherical charge penetration to the point charges had negligible effect on the results. However, adding spherical charge penetration to the DDEC6 point charges plus atomic dipoles decreased the conformation specific RMSE to 1.20 kcal mol⁻¹, which was dramatically better than any of the point charge only models.

What conclusions can be drawn from these results? We can definitely say the HD NACs were too small in magnitude and the DDEC3 NACs were too large in magnitude for this material.²¹ We can also say the Dipole charge is a limited concept, because it cannot be computed for some molecular conformations. Overall, this example illustrates some of the compromises involved in designing a general-purpose charge assignment method: (a) the molecular dipole moment can be reproduced exactly using a point charge plus atomic dipole model, but this makes the model more complicated than a point charge only model. (b) Fitting the electrostatic potential directly to a point charge model for each conformation leads to comparatively low RMSE values, but this degrades the conformational transferability as demonstrated by the ESP results. (c) The electrostatic potential can be more accurately reproduced by a (truncated) multipole model with charge penetration terms (e.g., D + SCP), but this results in more complicated force-field terms.

3.5 Exact for isolated atomic ion limit

The isolated atomic ion limit corresponds to spatially separated and negligibly overlapping atomic ions, such as occurs when each atomic ion is ∼10 Å or more from all other atoms. It is the only situation for which an exact {w_A(r_A)} is uniquely defined. Specifically, in the isolated atomic ion limit, the atomic weighting factors should equal the spherical average of each isolated density: w_A(r_A) = ρ^avg_A(r_A).

Consider the specific example of a periodic cubic array with Na and F ions at alternating vertices. By symmetry, the atomic dipole, quadrupole, and octupole moments are zero. Since atomic hexadecapole moments are the leading non-zero atomic multipoles, the atomic electron distributions are approximately (but not exactly) spherically symmetric. Thus, in this idealized example, the total electron density can be approximated as the sum of individual spherical ion densities:

We consider two limiting cases: (a) a periodic array having a 20 Å distance between nearest Na atoms and (b) the PBE-optimized low energy crystal structure having 2.27 Å between nearest Na atoms.

There are two possible strategies for the IH method: (i) use reference ions for the isolated ions without charge compensation (as done when the IH method was introduced¹³) or (ii) use reference ions that mimic the ion shapes in condensed crystals by including charge compensation effects (as done in later modifications of the IH method^{17,18,20,42,90}). If choice (i) is made, the reference ion shapes will match the ions in example (a). If choice (ii) is made, the reference ion shapes will match the ions in example (b). Due to charge compensation and electrostatic screening effects, anions in the condensed phase are more contracted than their isolated gas-phase counterparts. Therefore, the IH method must choose which of these two limits to reproduce. Moreover, to yield w_A(r_A) = ρ^avg_A(r_A) the IH reference ions would also need to be computed using a similar exchange-correlation theory and basis set as used to study the material of interest.

The DDEC6 method accurately reproduces both limits, because the reference ion densities are conditioned to the material of interest. In the isolated atomic ion limit (structure (a)), w^DDEC6_A(r_A) = ρ^cond_A(r_A) = [small rho, Greek, macron]^ref_A(r_A,q^ref_A)〈ρ(r⃑)/[small rho, Greek, macron]^ref(r⃑)〉_rA= ρ^avg_A(r_A) and additional conditioning steps do not alter w^DDEC6_A(r_A). Because the DDEC6 method derives {w^DDEC6_A(r_A)} from partitions of ρ(r⃑) (i.e., a conditioning process), the DDEC6 method returns w_A(r_A) = ρ^avg_A(r_A) in the isolated atomic ion limit even when the DDEC6 reference ions are computed using a different exchange-correlation theory, different basis sets, and different local chemical environment than used in the system of interest! In the optimized crystal geometry (structure (b)), the symmetry makes the atomic dipole, quadrupole, and octupole moments zero. Consequently, ρ_A(r⃑_A) ≈ ρ^avg_A(r_A) which implies 〈(ρ(r⃑)/W(r⃑))²〉_rA− 〈ρ(r⃑)/W(r⃑)〉_rA² ≈ 0. Since {w^DDEC6_A(r_A)} are computed via conditioning steps that make 〈ρ(r⃑)/W(r⃑)〉_rA ≈ 1, this means the conditioning process combined with the crystal symmetry makes ρ(r⃑)/W(r⃑) ≈ 1. This gives ρ_A(r⃑_A) ≈ w^DDEC6_A(r_A) ≈ ρ^avg_A(r_A). Thus, the DDEC6 method accurately recovers the nearly exact limit for both structures (a) and (b).

Now consider the NACs for structures (a) and (b). Since the atoms are fully separated in structure (a), the computed Bader, DDEC6, HD, and ISA NACs were the same within an integration tolerance. In this fully separated atomic ion limit, the results depend only on the exchange-correlation theory used to generate the electron distribution. For structure (a), the atomic charge magnitudes were 1.00 (Hartree–Fock method), 0.56 (HSE06 (ref. 91) functional), and 0.58 (PBE functional). Now consider the PBE-optimized low-energy crystal structure having 2.27 Å between nearest Na atoms. Analysis of experimental data shows the NaF crystal is mostly ionic.⁹² The computed HD (0.28) and ISA (0.48) atomic charge magnitudes are too small chemically. The Bader (0.86) and DDEC6 (0.85) atomic charge magnitudes are more reasonable. Thus, even when considering these simple NaF structures, some key advantages of DDEC6 over the HD and ISA methods are apparent.

Although the DDEC6 and HD methods returned the same NACs in the isolated atomic ion limit, there is a crucial distinction between these two approaches. Consider a situation intermediate between the isolated atomic ion limit and strongly overlapping atoms. For example, an array of atomic ions in which the adjacent atoms are ∼6 Å apart. In this situation, the adjacent atoms have small but non-zero overlaps. As explained above, w^DDEC6_A(r_A) approaches ρ^avg_A(r_A) in the isolated atomic ion limit, while w^HD_A(r_A) does not (for charged atoms). Thus, the DDEC6 atomic weighting factors approach the correct limit as the atoms are mostly separated while the HD atomic weighting factors do not.

3.6 Diisopropylammonium bromide ferroelectric crystal

Diisopropylammonium bromide (DIPAB) is a remarkable organic ferroelectric compound in numerous aspects.^93,94 It has a high spontaneous polarization (measured by electric field cycling), a high dielectric constant, and a low coercive field.^93,94 These characteristics are desirable for many potential applications.^93,94

In this section, we study several DIPAB crystal phases. The crystal structures were obtained using X-ray diffraction: P2₁ (293 K), P2₁2₁2₁ (293 K), and P2₁/m (438 K).⁹⁴ Phase P2₁ is ferroelectric.^94,95 Phase P2₁2₁2₁ transitions to P2₁ when heated, and phase P2₁ transitions to P2₁/m when heated more.^94,95 P2₁/m transitions to P2₁ when cooled but not to P2₁2₁2₁ if cooled further.^94,95 The P2₁ ferroelectric phase is stable at temperatures from 90 to 425 K.⁹⁵

A ferroelectric phase has a net dipole moment per unit cell. The polarization density is the average dipole moment per unit volume. Fu et al. performed measurements and calculations of the polarization density of DIPAB phase P2₁.⁹⁴ Using a Sawyer–Tower circuit at 25 Hz, they measured a polarization density of approximately 11 μC cm⁻².⁹⁴ Using a pyroelectric technique, they measured a spontaneous polarization density of approximately 23 μC cm⁻².⁹⁴ Using DFT calculations, they computed a polarization density of 24 μC cm⁻²; however, insufficient computational details were provided.⁹⁴ The computed polarization density depends upon the particular ferroelectric motion. Fu et al. did not specify the particular ferroelectric motion associated with their polarization density measurements and calculations.⁹⁴ Jiang et al. measured DIPAB ferroelectric hysteresis scaling behavior as a function of applied electric field and cycling frequency.⁹⁶ They measured saturation polarization densities of around 5 to 11 μC cm⁻².⁹⁶ Alsaad et al. used the Berry phase method to compute a DIPAB polarization of 23 μC cm⁻², but they did not compute the continuous deformation pathway (i.e., ferroelectric motion) associated with this polarization value.⁹⁷

We performed DFT calculations using the PBE exchange-correlation functional on these three phases using VASP. For P2₁ and P2₁2₁2₁, the lattice constants and angles where held at the experimental values and the positions of all atoms were fully relaxed. For P2₁/m, the experimental X-ray structure shows disordered atoms,⁹⁴ so we created a super-cell including both components of the disordered structure and relaxed the atomic positions. Fig. 7 shows the structures of these three phases. As shown in the lower half of Fig. 7, the DDEC6 NACs had good transferability between these three phases.

Fig. 7 DIPAB molecular crystals. Top: colored by element (gray: C; white: H; red: Br; blue: N). Bottom: colored by DDEC6 NAC (red: −0.61 to −0.67; yellow: −0.42 to −0.43; green: −0.26 to −0.33; light blue: 0.02 to 0.06; medium blue: 0.10 to 0.17; dark blue: 0.21 to 0.29).

Table 15 summarizes computed NACs for these three DIPAB phases. NACs for atoms with similar connectivity were averaged. The MUD was computed as the mean unsigned deviation between the NAC of an individual atom and its connectivity-averaged NAC. As shown in Table 15, the MUD was small for the DDEC6, HD, and CM5 methods and larger for the Bader method. The DDEC6 NACs for the three phases were similar. Only the C̲(H)₃(C) NACs exhibited a maximum difference > 0.1e among the three phases. Because the results for all three phases were similar, Table 15 only shows the HD, CM5, and Bader results for the ferroelectric P2₁ phase. (The individual NACs for all three phases are given in the ESI.†) As expected, the HD NACs are lower in magnitude than NACs computed using the other methods. All methods except HD gave a negative NAC for N. All methods gave a positive NAC for the H atoms bound to N. All methods except Bader gave a positive NAC for the H atoms bound to C. All methods gave a positive C̲(H)(C)₂(N) NAC. All methods except Bader gave a negative C̲(H)₃(C) NAC. All methods gave a negative Br NAC. The Bader and CM5 methods gave a more negative NAC for N than for Br, while the DDEC6 method gave a more negative NAC for Br than for N.

Table 15 Computed NACs for three DIPAB phases

Atom	Phase→	P21	P21	P21	P21	P212121	P21/m
	Connectivity	DDEC6	HD	CM5	Bader	DDEC6	DDEC6
N	N(H)2(C)2	−0.326	0.034	−0.628	−1.146	−0.324	−0.261
H	H2N	0.291	0.099	0.359	0.517	0.290	0.254
H	HC	0.048	0.034	0.123	−0.011	0.049	0.058
H	H3C	0.160	0.037	0.111	−0.103	0.160	0.126
C	C̲(H)(C)2(N)	0.266	0.062	0.009	0.334	0.268	0.227
C	C̲(H)3(C)	−0.533	−0.118	−0.322	0.373	−0.535	−0.429
Br	Br	−0.668	−0.394	−0.402	−0.790	−0.674	−0.617
MUD		0.005	0.003	0.003	0.026	0.004	0.004

---

3.7 Collinear and non-collinear magnetic materials

For both the collinear and non-collinear magnetic systems described below, DDEC6 ASMs were computed with the method of Manz and Sholl⁹⁸ using the DDEC6 atomic electron distributions and the recommended value χ_spin = 0.5. Electron and spin density grids with a uniform spacing of ~0.14 bohr were used. Spherical averages were computed over 100 uniformly spaced radial shells up to a cutoff radius of 5 Å.

3.7.1 Collinear magnetism. Fig. 8 compares DDEC6 to DDEC3 ASMs for all of the collinear magnetic materials studied in the article by Manz and Sholl⁵ that introduced the DDEC3 method: [Cr(CN)₆]³⁻ spin quartet, [Cu₂N₁₀C₃₆H₅₂]²⁺ spin triplet, anti-ferromagnetic CuBTC metal–organic framework, anti-ferromagnetic Fe₂O₃ crystal, anti-ferromagnetic Fe₂SiO₄ crystal, anti-ferrimagnetic Fe₃O₄ crystal (PBE functional and PBE+U_eff (U_eff = 4.0 eV) functionals), Fe₃Si crystal, [GdI]²⁺ using both SDD and planewave basis sets, the MgI, MoI, SnI, TeI, and TiI molecules using both SDD and planewave basis sets, and the ozone triplet spin state and the ozone +1 cation doublet spin state using the PW91, B3LYP, and CCSD methods. The same geometries, electron distributions, and spin distributions were used as input for DDEC6 analysis as were previously used for DDEC3 analysis.⁵ As shown in Fig. 8, the DDEC6 and DDEC3 ASMs are essentially identical. This follows the observation that ASMs are usually less sensitive than NACs to the choice of atomic population analysis method.^5,98,99

Fig. 8 Comparison of DDEC3 and DDEC6 atomic spin moments for systems with collinear magnetism. The black line has a slope of 1 and an intercept of 0.

As an additional example, we consider the Mn₁₂-acetate (formula unit Mn₁₂C₃₂H₅₆O₄₈) single molecule magnet illustrated in Fig. 9 (left panel). Mn₁₂-acetate is one of the most widely studied of all single molecule magnets since its synthesis and discovery by Lis.^100–102 We performed calculations in GAUSSIAN 09 using the PBE functional with LANL2DZ¹⁰³ basis sets and in VASP using the PBE/planewave method. Experiments support a conceptual model with an S = 10 and S_Z = 10 ground state.¹⁰⁴ Accordingly, we set S_Z = 10 as a constraint on the GAUSSIAN 09 and VASP electron and spin distributions we computed. In VASP, we optimized the atomic positions and used the experimental lattice parameters of Farrell et al.¹⁰⁵ (Cambridge Structural Database ID BESXAA). In GAUSSIAN 09, we used an isolated molecule and optimized the atomic positions. As shown in Table 16, the DDEC6 ASMs computed using the PBE functional with both LANL2DZ and planewave basis sets were in good agreement with Robinson et al.'s¹⁰⁶ polarized neutron diffraction experiments and Pederson and Khanna's¹⁰⁷ PBE computations. The DDEC6 results are similar to the NACs and ASMs we previously obtained for Mn₁₂-acetate using DDEC3.⁸⁶ ASMs on all atoms except Mn atoms were almost negligible in magnitude (i.e., ≤0.033 (planewave) and ≤0.077 (LANL2DZ)), which agrees with the experimental finding that “there is no evidence for net [magnetic] moments on the oxygen atoms”.¹⁰⁶ The magnetic anisotropy barrier of a single molecule magnet is the energy required to flip the magnetic moment orientation relative to the molecular structure.¹⁰¹ We computed this barrier by performing 62 single-point spin–orbit coupling calculations in VASP, where the electron and spin distributions were kept constant while the magnetic direction was rotated (by varying the SAXIS parameter in VASP). A 1 × 1 × 2 Monkhorst–Pack k-point mesh was used with Fermi smearing (smearing width = 0.05 eV). To ensure the computed gradient field (and hence exchange-correlation energy density) does not break the crystal symmetry, a spherical cutoff was applied to the reciprocal vectors, {G⃑}, such that only reciprocal vectors with |G⃑| < G_cut are included where G_cut defines a sphere enclosed within the parallelepiped defined by the reciprocal space FFT grid.²⁶ The one-center PAW charge densities were stored up to an angular momentum [small script l] = 6. As shown in Fig. 9 (middle panel), the spin–orbit coupling potential energy surface had global energy minima at the poles and a global energy maximum at the equator with no other local energy minima or maxima. This yielded a magnetic anisotropy barrier of 5.13 meV (59.5 K), which is in good agreement with Fort et al.'s¹⁰⁸ experimental value of 60–62 K and Pederson and Khanna's computed value of 55.6–55.8 K.¹⁰⁷ As shown in Fig. 9 (right panel), DDEC6 NACs computed with the LANL2DZ basis set are nearly identical to those computed using the planewave basis set. This shows the DDEC6 NACs are not overly sensitive to the basis set choice.

Fig. 9 Left: Atomic structure of Mn₁₂-acetate single molecule magnet. Mn type 1 (blue), Mn type 2 (red), Mn type 3 (yellow). In the minimum energy conformations, the Mn ASM vectors are perpendicular to the plane of the page. Middle: Computed spin–orbit coupling potential energy surface. Right: Comparison of DDEC6 NACs computed with LANL2DZ and planewave basis sets.

Table 16 Comparison of DDEC6 ASMs for Mn atoms in the Mn₁₂-acetate single molecule magnet to prior experiments and computations. Our computed magnetic anisotropy barrier is also compared to prior experiments and computations

Atom type	DDEC6 PBE planewave	DDEC6 PBE LANL2DZ	Experiments	Pederson Khanna PBE^b
a Polarized neutron diffraction experiments of Robinson et al.^106b Pederson and Khanna using integration of the spin density over spheres of 2.5 bohr radius to compute the ASMs.^107c Fort et al.^108
Atomic spin moment
Mn type 1	−2.80	−2.56	−2.34 ± 0.13^a	−2.6
Mn type 2	3.82	3.63	3.79 ± 0.12^a	3.6
Mn type 3	3.81	3.57	3.69 ± 0.14^a	3.6
Magnetic anisotropy barrier (Kelvin)
	59.5		60–62^c	55.6–55.8

---

3.7.2 Noncollinear magnetism. The left panel of Fig. 10 shows the globally minimized geometry and non-collinear magnetic structure of the Fe₄O₁₂N₄C₄₀H₅₂ noncollinear single molecule magnet. We computed DDEC6 NACs and the atomic spin magnetization vectors {M⃑_A} for this material using the same electron and spin magnetization density files as ref. 98. As shown in the center panel, the DDEC6 NACs followed a similar trend and magnitude as the DDEC3 NACs reported in ref. 5. The atomic spin magnitudes are the magnitudes of the atomic spin magnetization vectors: M_A = |M⃑_A|. As shown in the right panel, the DDEC6 atomic spin magnitudes were virtually identical to the DDEC3 values. The total wall time from CHARGEMOL program start (before input file reading) to end (after output printing finished) was 16.3 minutes for this calculation run on a single processor core in Intel Xeon E5-2680v3 at the Comet supercomputing cluster. This works out to 8.7 seconds per atom. This calculation utilized a volume of 2.9 × 10⁻³ bohr³ per grid point. These results demonstrate that DDEC6 is well-suited for quantifying NACs and atomic spins in non-collinear magnets.

Fig. 10 Left: Fe₄O₁₂N₄C₄₀H₅₂ noncollinear single molecule magnet structure reproduced with permission from ref. 98 (© American Chemical Society 2011). The arrows show the magnitude and direction of the atomic spin magnetization vectors on each atom. The atomic spin magnitudes are small on all atoms except the four iron atoms. Center: Comparison of DDEC3 and DDEC6 net atomic charges. Right: Comparison of DDEC3 and DDEC6 atomic spin magnitudes. The black lines have a slope of 1 and an intercept of 0.

3.8 Summary of systems studied

Table 17 summarizes the types of systems analyzed using the DDEC6 method. The 231 systems studied include those in this article and the previous article⁴ introducing the DDEC6 method. In this context, the same material investigated with different levels of theory counts as different systems. Four different system classifications are summarized in Table 17: (a) according to magnetic property, (b) according to material type, (c) according to exchange-correlation theory, and (d) according to basis set type.

Table 17 Number of systems and elements studied with the DDEC6 method

	Number of systems
According to magnetic property
Non-magnetic	185
Collinear magnetism	44
Noncollinear magnetism	2
According to material type
Small molecules (<100 atoms)	118
Nonporous solids	69
Porous solids	29
Large molecules (≥100 atoms)	6
Surface slabs or sheets	4
Stretched NaF periodic arrays	3
1-D periodic systems	2
According to exchange-correlation theory
PBE	112
B3LYP	71
PW91	19
M06L	17
CCSD	3
CAS-SCF	2
SAC-CI	2
CISD	1
DFT + dispersion	1
DFT+U	1
HF	1
HSE06	1
According to basis set type
Planewave	123
Gaussian function	108
Number of distinct chemical elements studied = 44

---

Even though the Mn₁₂-acetate single molecule magnet exhibits collinear magnetism, we had to compute its fully noncollinear magnetization densities in order to compute the spin–orbit coupling potential energy surface (Fig. 9) and magnetic anisotropy barrier. Therefore, the two noncollinear magnetism systems included the noncollinear single molecule magnet of Section 3.7.2 as well as the PBE/planewave noncollinear magnetism calculation of the Mn₁₂-acetate single molecule magnet. The PBE/LANL2DZ analysis of the Mn₁₂-acetate single molecule magnet was a collinear magnetism calculation.

We showed the DDEC6 method gives good performance across a broad range of material types: (a) small molecules, (b) nonporous solids, (c) porous solids, (d) large molecules, (e) surface slabs or sheets, (f) near the isolated atomic ion limit (e.g., stretched NaF periodic arrays), and (g) 1-D periodic systems. The 1-D periodic systems included the B-DNA decamer and the BN nanotube. Some materials are described in the ESI.† Fig. S1 in Section S2 of the ESI† compares DDEC6 to DDEC3 NACs for 14 different materials comprised almost entirely of surface atoms. The differences between DDEC6 and DDEC3 NAC values was minor for these materials, but some statistically significant differences occur. Fig. S2 in Section S3 of the ESI† compares DDEC6 to DDEC3 NACs for three solid surface slabs: (a) Mo₂C(110) surface with K adatom, (b) NaF(001) surface, and (c) SrTiO₃(100) surface. Atomic dipoles and quadrupoles for the SrTiO₃(100) surface and bulk crystal are listed in Table S3.† These results were included in the ESI,† because the differences between DDEC3 and DDEC6 results was minor. For SrTiO₃ and NaF, a comparison between DDEC3, DDEC6, Bader, and IH NACs for the bulk crystals is shown in Tables S4 and S5,† respectively.

Yang and Manz have previously reported DDEC6 NACs, ASMs, and bond orders for 96 organometallic complexes containing >100 atoms.¹⁰⁹ Including these would have significantly raised the total for large molecules.

The exchange-correlation theories we considered included coupled-cluster methods (e.g., CCSD¹¹⁰ and SAC-CI^111,112), configuration interaction methods (e.g., CAS-SCF¹¹³ and CISD¹¹⁴), GGA functionals (e.g., PBE³⁴ and PW91 (ref. 115)), a meta-GGA functional (e.g., M06L³³), a hybrid functional (e.g., B3LYP^49,50), a range-separated hybrid functional (e.g., HSE06 (ref. 91)), a DFT+U method (e.g., PBE+U_eff¹¹⁶), a dispersion-corrected functional (e.g., PBE + D3 (ref. 117)), and the Hartree–Fock (HF) method. We do not recommend the HF method, because it neglects the correlation energy. Any of the correlated methods, whether wavefunction or DFT-based, can be utilized as long as it provides acceptable accuracy for the system's energy, electron distribution, and spin magnetization distribution. Electron and spin distributions for some of these levels of theory were first computed in earlier papers on the DDEC methods and were re-used here; the reader is referred to those earlier papers for calculation details.^5,23 For example, Fig. S1 of the ESI† comparing DDEC6 to DDEC3 NACs includes CCSD, SAC-CI, and CAS-SCF results for different ozone spin states, and the reader is referred to earlier DDEC papers^5,23 for the quantum chemistry calculation details of these systems. The PBE + D3 dispersion-corrected functional calculation was for the Li₃ cluster calculation.⁴ The CISD calculation was for the H₂ triplet with constrained bond length.⁴ The PBE+U_eff calculation was for the Fe₃O₄​ (magnetite) bulk crystal appearing as a datapoint in Fig. 8. The HSE06 and HF results were for the stretched NaF periodic array discussed in Section 3.5 above.

Because the DDEC6 results are a functional of the electron and spin magnetization distributions, the results have no explicit basis set dependence. As shown Table 17, our DDEC6 calculations were almost evenly divided between planewave and Gaussian function basis sets.

The materials we studied in this article and the previous article⁴ introducing the DDEC6 method included 44 different chemical elements. However, the number of distinct chemical elements analyzed to date using the DDEC6 method (in works beyond these two articles) extends well beyond this number. The important message is that the DDEC6 method works well across a broad range of chemical elements.

4. Conclusions

We tested the performance of the DDEC6 method across a diverse set of periodic and non-periodic materials with numerous comparisons to experimental data. Our computational tests showed the DDEC6 method has substantially better overall performance than the DDEC3 method. Therefore, we recommend the DDEC6 method be used instead of the DDEC3 method. The DDEC6 method's excellent performance across a wide range of important properties and material types makes it suited for use as a default method in quantum chemistry programs. Our computational tests showed the DDEC6 NACs are well-suited both for describing chemical electron transfer trends between atoms in materials and for constructing flexible force-fields for classical atomistic simulations of complex materials.

We first studied electron transfer trends in several solids and clusters. As pointed out by Wang et al.,²¹ the DDEC3 method predicts the incorrect electron transfer sign for the transition metal atom for the delithiation of solid LiCoO₂ to CoO₂. The DDEC6 method fixes this problem. We also showed the DDEC6 method yields reasonable electron transfer trends for the LiTiS₂, TiS₂, LiTiO₂, LiTi₂O₄, LiMn₂O₄, MnO₂, and Li₃RuO₂ solids. For several Pd-containing alloys, we compared the electron transfer direction predicted by element electronegativities to computed NACs: the Bader, DDEC3, and DDEC6 NACs followed the Pauling scale electronegativity trends while the HD and CM5 NACs did not. For (MgO)_n (n = 1 to 6) clusters, we found the DDEC6 method exhibits overall better performance than DDEC3 for reproducing the electrostatic potential and dipole moments.

We then compared NACs to spectroscopic results for various materials. For natrolite, the DDEC6 NACs were smaller in magnitude than the DDEC3 NACs. For this material, the DDEC6 NACs were closer than DDEC3 NACs to the NACs extracted from high-resolution diffraction data using kappa refinement (with the exception of the Na atom which was not refined). DDEC3 and DDEC6 were both in excellent agreement with formamide NACs extracted from high-resolution diffraction data using spherical atom refinement. For a series of Ti-containing compounds, core-electron binding energy shifts were approximately linearly correlated to the DDEC6, HD, and Bader NACs but not to the DDEC3 and CM5 NACs. All five charge assignment methods gave reasonably good correlations between core electron binding shifts and computed NACs for the Mo-containing and Fe-containing compounds.

For 13 materials studied at the low energy conformation, the DDEC6 NACs reproduced the electrostatic potential slightly better than the DDEC3 NACs in 8 of the 13 materials. This shows the DDEC6 NACs are suited for constructing force-fields for materials containing small molecules, porous solids, water, and large biomolecules. A detailed study across various conformations of Li₂O and five carboxylic acids showed the DDEC6 NACs have excellent conformational transferability and are ideally suited for constructing flexible force-fields to approximately reproduce the electrostatic potential across various system conformations.

We then studied a series of materials containing surface atoms. For a series of systems comprised almost entirely of surface atoms (see ESI Section S2†), the DDEC6 and DDEC3 NACs exhibited similar trends with some statistically significant differences in NAC values. Additional tests for three solid surfaces (K adatom on a Mo₂C (110) surface, NaF(001) slab, and SrTiO₃(100) slab) showed the DDEC6 method maintains a consistent treatment of surface and buried atoms (see ESI Section S3†). As explained in Section 3.5 above, the DDEC6 method is asymptotically exact in the isolated atomic ion limit.

As an example of a material with interesting phase change behavior, we studied three DIPAB crystal phases. One of these three phases is ferroelectric. Prior experimental studies have noted the exceptionally good performance of DIPAB for an organic ferroelectric.^93,94 We showed the DDEC6 NACs for this material are chemically reasonable and have good transferability among the three DIPAB crystal phases.

Finally, we examined materials with collinear and non-collinear magnetism and found the DDEC6 atomic spin moments (ASMs) are essentially identical to the DDEC3 ASMs. For the Mn₁₂-acetate single molecule magnet, the computed DDEC6 ASMs were in excellent agreement with previous experiments¹⁰⁶ and computations.¹⁰⁷ We computed the spin–orbit coupling potential energy surface for this material and found the resulting magnetic anisotropy barrier (5.13 meV) to be in excellent agreement with previous experiments¹⁰⁸ and computations.¹⁰⁷ The NACs and ASMs are computed efficiently. For example, it took only 16.3 minutes to compute them on a single processor core when analyzing a Fe₄O₁₂N₄C₄₀H₅₂ single molecule magnet with non-collinear magnetism containing 112 atoms.

Acknowledgements

Supercomputing resources were provided by the Extreme Science and Engineering Discovery Environment (XSEDE). XSEDE is funded by NSF grant ACI-1053575. XSEDE project grant TG-CTS100027 provided allocations on the Trestles and Comet clusters at the San Diego Supercomputing Center (SDSC) and the Stampede cluster at the Texas Advanced Computing Center (TACC).

References

1. R. S. Mulliken, J. Chem. Phys., 1955, 23, 1833–1840.
2. G. Bruhn, E. R. Davidson, I. Mayer and A. E. Clark, Int. J. Quantum Chem., 2006, 106, 2065–2072.
3. A. E. Reed, R. B. Weinstock and F. Weinhold, J. Chem. Phys., 1985, 83, 735–746.
4. T. A. Manz and N. Gabaldon Limas, Introducing DDEC6 Atomic Population Analysis: Part 1. Charge Partitioning Theory and Methodology, RSC Adv., 2016, 6 10.1039/c6ra04656h.
5. T. A. Manz and D. S. Sholl, J. Chem. Theory Comput., 2012, 8, 2844–2867.
6. B. H. Besler, K. M. Merz and P. A. Kollman, J. Comput. Chem., 1990, 11, 431–439.
7. U. C. Singh and P. A. Kollman, J. Comput. Chem., 1984, 5, 129–145.
8. C. Campana, B. Mussard and T. K. Woo, J. Chem. Theory Comput., 2009, 5, 2866–2878.
9. R. F. W. Bader, Acc. Chem. Res., 1975, 8, 34–40.
10. R. F. W. Bader, P. J. MacDougal and C. D. H. Lau, J. Am. Chem. Soc., 1984, 106, 1594–1605.
11. F. L. Hirshfeld, Theor. Chim. Acta, 1977, 44, 129–138.
12. A. V. Marenich, S. V. Jerome, C. J. Cramer and D. G. Truhlar, J. Chem. Theory Comput., 2012, 8, 527–541.
13. P. Bultinck, C. Van Alsenoy, P. W. Ayers and R. Carbo-Dorca, J. Chem. Phys., 2007, 126, 144111.
14. E. R. Davidson and S. Chakravorty, Theor. Chim. Acta, 1992, 83, 319–330.
15. D. Geldof, A. Krishtal, F. Blockhuys and C. Van Alsenoy, J. Chem. Theory Comput., 2011, 7, 1328–1335.
16. T. Verstraelen, P. W. Ayers, V. Van Speybroeck and M. Waroquier, J. Chem. Theory Comput., 2013, 9, 2221–2225.
17. D. E. P. Vanpoucke, P. Bultinck and I. Van Driessche, J. Comput. Chem., 2013, 34, 405–417.
18. T. Bucko, S. Lebegue, J. Hafner and J. G. Angyan, J. Chem. Theory Comput., 2013, 9, 4293–4299.
19. K. Finzel, A. M. Pendas and E. Francisco, J. Chem. Phys., 2015, 143, 084115.
20. T. Bucko, S. Lebegue, J. G. Angyan and J. Hafner, J. Chem. Phys., 2014, 141, 034114.
21. B. Wang, S. H. L. Li and D. G. Truhlar, J. Chem. Theory Comput., 2014, 10, 5640–5650.
22. T. C. Lillestolen and R. J. Wheatley, Chem. Commun., 2008, 5909–5911.
23. T. A. Manz and D. S. Sholl, J. Chem. Theory Comput., 2010, 6, 2455–2468.
24. J. Cioslowski, J. Am. Chem. Soc., 1989, 111, 8333–8336.
25. P. Ghosez, J. P. Michenaud and X. Gonze, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 58, 6224–6240.
26. J. Hafner, J. Comput. Chem., 2008, 29, 2044–2078.
27. G. Kresse and J. Furthmuller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
28. P. E. Blochl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953–17979.
29. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
30. W. Tang, E. Sanville and G. Henkelman, J. Phys.: Condens. Matter, 2009, 21, 084204.
31. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. J. Montgomery, J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, T. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, O. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski and D. J. Fox, Gaussian 09, Revision C.01, Gaussian, Inc., Wallingford CT, 2010.
32. T. A. Manz and N. Gabaldon Limas, Chargemol program for performing DDEC analysis, version 3.4, accessed December 2015, http://ddec.sourceforge.net.
33. Y. Zhao and D. G. Truhlar, J. Chem. Phys., 2006, 125, 194101.
34. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
35. A. Yamada and M. Tanaka, Mater. Res. Bull., 1995, 30, 715–721.
36. J. Rodriguez-Carvajal, G. Rousse, C. Masquelier and M. Hervieu, Phys. Rev. Lett., 1998, 81, 4660–4663.
37. A. Paolone, P. Roy, G. Rousse, C. Masquelier and J. Rodriguez-Carvajal, Solid State Commun., 1999, 111, 453–458.
38. Electronegativity, in CRC Handbook of Chemistry and Physics, ed. D. R. Lide, CRC Press, Boca Raton, Florida, 88th edn, 2007–2008, p. 9.81.
39. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297–3305.
40. J. M. Zuo, M. OKeeffe, P. Rez and J. C. H. Spence, Phys. Rev. Lett., 1997, 78, 4777–4780.
41. J. M. Zuo, Rep. Prog. Phys., 2004, 67, 2053–2103.
42. D. E. P. Vanpoucke, I. Van Driessche and P. Bultinck, J. Comput. Chem., 2013, 34, 422–427.
43. N. K. Hansen and P. Coppens, Acta Crystallogr., Sect. A: Cryst. Phys., Diffr., Theor. Gen. Crystallogr., 1978, 34, 909–921.
44. Y. A. Abramov, A. V. Volkov and P. Coppens, Chem. Phys. Lett., 1999, 311, 81–86.
45. P. Coppens, T. N. Guru-Row, P. Leung, E. D. Stevens, P. J. Becker and Y. W. Yang, Acta Crystallogr., Sect. A: Cryst. Phys., Diffr., Theor. Gen. Crystallogr., 1979, 35, 63–72.
46. E. D. Stevens, Acta Crystallogr., Sect. B: Struct. Crystallogr. Cryst. Chem., 1978, 34, 544–551.
47. E. D. Stevens, J. Rys and P. Coppens, J. Am. Chem. Soc., 1978, 100, 2324–2328.
48. N. E. Ghermani, C. Lecomte and Y. Dusausoy, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 53, 5231–5239.
49. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
50. P. J. Stephens, F. J. Devlin, C. F. Chabalowski and M. J. Frisch, J. Phys. Chem., 1994, 98, 11623–11627.
51. T. H. Dunning, J. Chem. Phys., 1989, 90, 1007–1023.
52. U. Gelius, P. F. Heden, J. Hedman, B. J. Lindberg, R. Manne, R. Nordberg, C. Nordling and K. Siegbahn, Phys. Scr., 1970, 2, 70–80.
53. U. Gelius, Phys. Scr., 1974, 9, 133–147.
54. W. L. Jolly and W. B. Perry, J. Am. Chem. Soc., 1973, 95, 5442–5450.
55. Handbook of X-Ray and Ultraviolet Photoelectron Spectroscopy, ed. D. Briggs, Heyden, London, 1977, pp. 53–77.
56. Photoemission in Solids I. General Principles, ed. M. Cardona and L. Ley, Springer-Verlag, Berlin, 1978, pp. 60–75, 165–195.
57. J. F. Moulder, W. F. Stickle, P. E. Sobol and K. D. Bomben, Handbook of X-ray Photoelectron Spectroscopy, ULVAC-PHI Inc., Chigasaki, Japan, 1995. pp. 1–261.
58. J. C. Carver, G. K. Schweitzer and T. A. Carlson, J. Chem. Phys., 1972, 57, 973–982.
59. J. Wong, F. W. Lytle, R. P. Messmer and D. H. Maylotte, Phys. Rev. B: Condens. Matter Mater. Phys., 1984, 30, 5596–5610.
60. S. P. Cramer, T. K. Eccles, F. W. Kutzler, K. O. Hodgson and L. E. Mortenson, J. Am. Chem. Soc., 1976, 98, 1287–1288.
61. M. J. Guittet, J. P. Crocombette and M. Gautier-Soyer, Phys. Rev. B: Condens. Matter Mater. Phys., 2001, 63, 125117.
62. L. W. Li, M. R. Morrill, H. Shou, D. G. Barton, D. Ferrari, R. J. Davis, P. K. Agrawal, C. W. Jones and D. S. Sholl, J. Phys. Chem. C, 2013, 117, 2769–2773.
63. X-Ray Absorption: Principles, Applications, Techniques of EXAFS, SEXAFS and XANES, ed. D. C. Koningsberger and R. Prins, John Wiley & Sons, New York, 1988, pp. 1–673.
64. Y. Zuxiang, Rock Miner. Anal., 1984, 3, 231–238.
65. S. Cardamone, T. J. Hughes and P. L. A. Popelier, Phys. Chem. Chem. Phys., 2014, 16, 10367–10387.
66. D. M. Elking, L. Perera and L. G. Pedersen, Comput. Phys. Commun., 2012, 183, 390–397.
67. Q. Yang, D. Liu, C. Zhong and J.-R. Li, Chem. Rev., 2013, 113, 8261–8323.
68. B. R. Brooks, C. L. Brooks, A. D. Mackerell, L. Nilsson, R. J. Petrella, B. Roux, Y. Won, G. Archontis, C. Bartels, S. Boresch, A. Caflisch, L. Caves, Q. Cui, A. R. Dinner, M. Feig, S. Fischer, J. Gao, M. Hodoscek, W. Im, K. Kuczera, T. Lazaridis, J. Ma, V. Ovchinnikov, E. Paci, R. W. Pastor, C. B. Post, J. Z. Pu, M. Schaefer, B. Tidor, R. M. Venable, H. L. Woodcock, X. Wu, W. Yang, D. M. York and M. Karplus, J. Comput. Chem., 2009, 30, 1545–1614.
69. J. M. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman and D. A. Case, J. Comput. Chem., 2004, 25, 1157–1174.
70. C. I. Bayly, P. Cieplak, W. D. Cornell and P. A. Kollman, J. Phys. Chem., 1993, 97, 10269–10280.
71. D. Golze, J. Hutter and M. Iannuzzi, Phys. Chem. Chem. Phys., 2015, 17, 14307–14316.
72. T. Loiseau, C. Serre, C. Huguenard, G. Fink, F. Taulelle, M. Henry, T. Bataille and G. Ferey, Chem.–Eur. J., 2004, 10, 1373–1382.
73. H. Li, M. Eddaoudi, M. O'Keeffe and O. M. Yaghi, Nature, 1999, 402, 276–279.
74. T. Watanabe, T. A. Manz and D. S. Sholl, J. Phys. Chem. C, 2011, 115, 4824–4836.
75. H. J. C. Berendsen, J. P. M. Postma, W. F. von Gunsteren and J. Hermans, Intermolecular Forces, Reidel, Dordrecht, Holland, 1981, p. 331.
76. H. J. C. Berendsen, J. R. Grigera and T. P. Straatsma, J. Phys. Chem., 1987, 91, 6269–6271.
77. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, J. Chem. Phys., 1983, 79, 926–935.
78. W. L. Jorgensen, J. Am. Chem. Soc., 1981, 103, 335–340.
79. L. P. Lee, D. J. Cole, C.-K. Skylaris, W. L. Jorgensen and M. C. Payne, J. Chem. Theory Comput., 2013, 9, 2981–2991.
80. L. P. Lee, N. Gabaldon Limas, D. J. Cole, M. C. Payne, C.-K. Skylaris and T. A. Manz, J. Chem. Theory Comput., 2014, 10, 5377–5390.
81. A. H. Farmahini, D. S. Sholl and S. K. Bhatia, J. Am. Chem. Soc., 2015, 137, 5969–5979.
82. N. Foloppe and A. D. MacKerell, J. Comput. Chem., 2000, 21, 86–104.
83. W. D. Cornell, P. Cieplak, C. I. Bayly, I. R. Gould, K. M. Merz, D. M. Ferguson, D. C. Spellmeyer, T. Fox, J. W. Caldwell and P. A. Kollman, J. Am. Chem. Soc., 1995, 117, 5179–5197.
84. S. Arai, T. Chatake, T. Ohhara, K. Kurihara, I. Tanaka, N. Suzuki, Z. Fujimoto, H. Mizuno and N. Niimura, Nucleic Acids Res., 2005, 33, 3017–3024.
85. T. Tsukamoto, Y. Ishikawa, M. J. Vilkas, T. Natsume, K. Dedachi and N. Kurita, Chem. Phys. Lett., 2006, 429, 563–569.
86. N. Gabaldon Limas, MS Thesis, New Mexico State University, 2014.
87. J. D. Roberts and W. T. Moreland, J. Am. Chem. Soc., 1953, 75, 2167–2173.
88. T. Verstraelen, E. Pauwels, F. De Proft, V. Van Speybroeck, P. Geerlings and M. Waroquier, J. Chem. Theory Comput., 2012, 8, 661–676.
89. T. Verstraelen, P. W. Ayers, V. Van Speybroeck and M. Waroquier, Chem. Phys. Lett., 2012, 545, 138–143.
90. T. A. Manz, J. Comput. Chem., 2013, 34, 418–421.
91. A. V. Krukau, O. A. Vydrov, A. F. Izmaylov and G. E. Scuseria, J. Chem. Phys., 2006, 125, 224106.
92. Z. W. Su and P. Coppens, Acta Crystallogr., Sect. A: Found. Crystallogr., 1995, 51, 27–32.
93. D. A. Bonnell, Science, 2013, 339, 401–402.
94. D. W. Fu, H. L. Cai, Y. M. Liu, Q. Ye, W. Zhang, Y. Zhang, X. Y. Chen, G. Giovannetti, M. Capone, J. Y. Li and R. G. Xiong, Science, 2013, 339, 425–428.
95. A. Piecha, A. Gagor, R. Jakubas and P. Szklarz, CrystEngComm, 2013, 15, 940–944.
96. C. L. Jiang, H. C. Lin, C. H. Luo, Y. Y. Zhang, J. Yang, H. Peng and C. G. Duan, J. Cryst. Growth, 2016, 438, 25–30.
97. A. Alsaad, I. A. Qattan, A. A. Ahmad, N. Al-Aqtash and R. F. Sabirianov, IOP Conf. Ser.: Mater. Sci. Eng., 2015, 92, 012017.
98. T. A. Manz and D. S. Sholl, J. Chem. Theory Comput., 2011, 7, 4146–4164.
99. E. Ruiz, J. Cirera and S. Alvarez, Coord. Chem. Rev., 2005, 249, 2649–2660.
100. T. Lis, Acta Crystallogr., Sect. B: Struct. Crystallogr. Cryst. Chem., 1980, 36, 2042–2046.
101. M. N. Leuenberger and D. Loss, Nature, 2001, 410, 789–793.
102. K. M. Mertes, Y. Suzuki, M. P. Sarachik, Y. Myasoedov, H. Shtrikman, E. Zeldov, E. M. Rumberger, D. N. Hendrickson and G. Christou, Solid State Commun., 2003, 127, 131–139.
103. P. J. Hay and W. R. Wadt, J. Chem. Phys., 1985, 82, 299–310.
104. A. Caneschi, D. Gatteschi, R. Sessoli, A. L. Barra, L. C. Brunel and M. Guillot, J. Am. Chem. Soc., 1991, 113, 5873–5874.
105. A. R. Farrell, J. A. Coome, M. R. Probert, A. E. Goeta, J. A. K. Howard, M. H. Lemee-Cailleau, S. Parsons and M. Murrie, CrystEngComm, 2013, 15, 3423–3429.
106. R. A. Robinson, P. J. Brown, D. N. Argyriou, D. N. Hendrickson and S. M. J. Aubin, J. Phys.: Condens. Matter, 2000, 12, 2805–2810.
107. M. R. Pederson and S. N. Khanna, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 60, 9566–9572.
108. A. Fort, A. Rettori, J. Villain, D. Gatteschi and R. Sessoli, Phys. Rev. Lett., 1998, 80, 612–615.
109. B. Yang and T. A. Manz, Theor. Chem. Acc., 2016, 135, 21.
110. G. E. Scuseria, C. L. Janssen and H. F. Schaefer, J. Chem. Phys., 1988, 89, 7382–7387.
111. H. Nakatsuji and K. Hirao, J. Chem. Phys., 1978, 68, 2053–2065.
112. H. Nakatsuji, Chem. Phys. Lett., 1979, 67, 334–342.
113. N. Yamamoto, T. Vreven, M. A. Robb, M. J. Frisch and H. B. Schlegel, Chem. Phys. Lett., 1996, 250, 373–378.
114. J. A. Pople, R. Seeger and R. Krishnan, Int. J. Quantum Chem., 1977, 12, 149–163.
115. J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D. J. Singh and C. Fiolhais, Phys. Rev. B: Condens. Matter Mater. Phys., 1992, 46, 6671–6687.
116. S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys and A. P. Sutton, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 57, 1505–1509.
117. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.

---

Footnote

† Electronic supplementary information (ESI) available: Tables listing spectroscopically measured core electron energies and Bader, CM5, DDEC3, DDEC6, HD NACs for Ti- and Mo-containing solids; plot comparing DDEC6 to DDEC3 NACs for 14 systems comprised almost entirely of surface atoms; DDEC6 and DDEC3 NACs for Mo₂C(110) surface with K adatom and NaF(001) surface; DDEC6 and DDEC3 NACs, atomic dipoles, and atomic quadrupoles for SrTiO₃(100) surface and bulk crystal; DDEC6, DDEC3, Bader, and IH NACs for NaF and SrTiO₃ bulk crystals; xyz files (which can be read using any text editor or the free Jmol visualization program downloadable from jmol.sourceforge.net) containing geometries, NACs, atomic dipoles and quadrupoles, fitted tail decay exponents, and ASMs. See DOI: 10.1039/c6ra05507a

---