Understanding chemistry with the symmetry-decomposed Voronoi deformation density charge analysis

The symmetry-decomposed Voronoi deformation density (VDD) charge analysis is an insightful and robust computational tool to aid the understanding of chemical bonding throughout all fields of chemistry. This method quantifies the atomic charge flow associated with chemical-bond formation and enables decomposition of this charge flow into contributions of (1) orbital interaction types, that is, Pauli repulsive or bonding orbital interactions; (2) per irreducible representation (irrep) of any point-group symmetry of interacting closed-shell molecular fragments; and now also (3) interacting open-shell (i.e., radical) molecular fragments. The symmetry-decomposed VDD charge analysis augments the symmetry-decomposed energy decomposition analysis (EDA) so that the charge flow associated with Pauli repulsion and orbital interactions can be quantified both per atom and per irrep, for example, for σ , π , and δ electrons. This provides detailed insights into fundamental aspects of chemical bonding that are not accessible from EDA.


| INTRODUCTION
The atomic charge is a fundamental and ubiquitous concept in chemistry.
This quantity is widely used as a key factor for understanding and predicting chemical bonding and reactivity. Over the years, chemists have developed numerous methods to quantify atomic charges. [1][2][3][4][5][6][7][8] A wellknown charge analysis scheme is the Voronoi deformation density (VDD) analysis, 1 which uses the spatial integration of the electronic density within an atomic domain, that is, the Voronoi cell of that atom (vide infra). This charge analysis scheme has shown to be insightful for understanding chemical bonding and reactivity in the fields of physical chemistry, inorganic and organic chemistry, as well as in supramolecular and biochemistry, and material sciences. [9][10][11][12][13] The versatile and compelling application of this method follows from the ability to analyze the flow of the electronic charge as a direct consequence of a chemical bond or interaction. VDD atomic charges (Q) are computed by spatial integration of the deformation density over the Voronoi cell of atom A, which is the space defined by the bond midplanes on and perpendicular to all bond axes between this atom and its neighboring atoms Herein, the deformation density Δρ r ð Þ ¼ ½ρ r ð Þ -P i ρ i r ð Þ is the density change going from a superposition of the original atomic densities at the positions of the molecule, coined the promolecule, to the actual density of that molecule. The promolecular density is defined as the sum of the (spherically averaged) ground-state atomic densities P i ρ i r ð Þ. This is the fictitious situation in which the charge density has not been affected by chemical-bonding effects and in which all atoms have a charge zero. In Equation 1, Q represents the amount of charge that, due to chemical bonding within the molecular system, flows from a position closer to another nucleus towards a position closer to nucleus A (Q A < 0), or vice versa (Q A > 0). 1,3a Besides its application for computing these regular VDD atomic charges Q, the VDD method is also a tool for studying atomic charge changes (ΔQ) caused by the bonding between (polyatomic) molecular fragments (Equation 2). 1 This fragment approach has been shown to be essential in studying weak (intermolecular) interactions, such as hydrogen bonds, in which case the charge rearrangements due to the weak interaction in the complex are small compared to the charge shifts resulting from the formation of strong covalent bonds. 1,13 ΔQ ¼ - Equation (2) yields the charge change (ΔQ) in the final density of the overall complex ρ complex r ð Þ relative to the sum of the initial molecular fragment densities P fragments,i ρ i r ð Þ. This reveals how the chemical bond between these polyatomic subsystems modifies the electron density distribution, that is, ΔQ represents the amount of electronic charge that flows into (ΔQ < 0) or out of (ΔQ > 0) the Voronoi cell of an atom due to the interaction between the molecular fragments. A useful aspect of the VDD method is that by using molecular symmetry ΔQ can be partitioned into contributions of the different irreducible representations (irreps) Γ of the point-group symmetry of the fragment-based promolecule (Equation 3) or a sub-point group thereof. For example, in planar (i.e., C s symmetric) molecular systems, this symmetry decomposition of the total ΔQ provides the insightful distinction between charge shifts within the σ (ΔQ σ ) and π (ΔQ π ) electronic systems, corresponding to the A' and A" irreps, respectively.
The symmetry-decomposed VDD charge analysis, implemented in the Amsterdam Density Functional (ADF) software code, 14-16 offers an even more detailed analysis tool when combined with the canonical energy decomposition analysis (EDA) scheme. 17 The EDA splits the total interaction energy (ΔE int ) into components of classical electrostatic interaction ΔV elstat , Pauli repulsion ΔE Pauli , and orbital interaction ΔE oi (Equation 4); note that only the two latter terms are associated with a change in the system's density that contributes to the VDD charge shifts ΔQ ( Figure 1A,B).
In the first step of the EDA, the classical electrostatic interaction  Figure 1C).
We have previously reported the theoretical background and performance assessment of the VDD method. 1 Herein, we extend on that work by presenting a generalization of the symmetry-decomposed

| ATOMIC CHARGE Q AND ELECTROSTATIC INTERACTIONS
The atomic charge Q of an atom in a molecular fragment can, for example, give insight into how it can electrostatically interact with another fragment. In the EDA scheme, ΔV elstat is the electrostatic interaction between the unperturbed densities of the geometrically deformed fragments, that is between ρ A and ρ B ( Figure 1A,B). It is often assumed that the trends in ΔV elstat can be explained by the atomic charges of the frontier atoms. However, this should be taken with care as atoms, and hence molecules, are no point charges but rather complex charge density distributions of positive nuclei and diffuse negative electron clouds. 18 We will demonstrate this by showing an example where atomic charges Q can explain trends in the electrostatic interaction ΔV elstat , and an example where this somewhat oversimplified assumption fails.
In Figure 2, the VDD atomic charges Q of the hydrogenbonded front atoms are shown for the 2-hydroxypyridine (red), 2-aminopyridine (blue), and 2-pyridone (black) isolated monomers in the geometry they adopt within the planar hydrogen-bonded dimers (the energy difference of the NH 2 groups being planar or pyramidal in the 2-aminopyridine dimer is only 0.05 kcal mol À1 , see Table S4). From 2-hydroxypyridine to 2-aminopyridine the hydrogen-bond donor functionality is changed from an OH to an (2) non-frontier atoms also contribute to the electrostatic interaction.
Thus, the electrostatic effects accumulated by the interaction of all atoms can be effective over long interatomic distances.

| ATOMIC CHARGE CHANGES ΔQ
Besides computing absolute atomic charges Q, the VDD method can also be used to study atomic charge changes ΔQ associated with chemical bonding between polyatomic fragments. In this section, we

| Decomposition by symmetry
The symmetry-decomposition of VDD charge changes ΔQ can be very informative to identify differences in reactivity and bonding of molecules. For example, this analysis can be applied to obtain insight into the effect of substituents on the reactivity and regioselectivity of benzene derivatives in electrophilic aromatic substitution reactions. 10 To this end, the electronic effect was investigated of replacing one of the arylic hydrogen atoms of benzene with a fluorine (F) atom, an amino (NH 2 ) group, or a nitro (NO 2 ) group, yielding fluorobenzene, aniline, and nitrobenzene, respectively. The NH 2 substituent is known to be activating (i.e., aniline is more reactive than benzene) and orthoand para-directing in electrophilic aromatic substitution reactions. 19 On the other hand, the NO 2 substituent is known to be deactivating (i.e., nitrobenzene is less reactive than benzene) and meta-directing.
Finally, the fluorine substituent is also deactivating compared to benzene, but is just like NH 2 orthoÀ/para-directing. These substantially different substituent effects can readily be understood from the symmetry-decomposed VDD charge analysis, using open-shell molecular fragments in C 2v symmetry (planarization of the slightly pyramidal NH 2 group in aniline costs only 0.8 kcal mol À1 , see Table S4). As such, we can distinguish between the charge changes ΔQ caused by induction (in the σ system) and resonance (in the π system), which corresponds in C 2v symmetry to the sum of the charge flow within the irreps A 1 + B 1 and A 2 + B 2 , respectively. 20 The analysis of the changes in atomic charges ΔQ resulting from the electron-pair bond formation between the phenyl radical (PhÁ) and the various substituent radicals (ÁNH 2 , ÁH, ÁF, or ÁNO 2 ) is shown in Figure 3. We indeed confirm that the NH 2 substituent is activating by net electron donation to the phenyl group, while F and NO 2 are deactivating by withdrawing electronic charge from the phenyl ring. Next, we separate ΔQ into inductive and resonance effects of the substituents, by the σ and π electrons, respectively. The inductive σ component (ΔQ σ ) shows that in the case of benzene the C H bond formation leads to a net electron flow from the H substituent to the phenyl ring, which is in line with the lower electronegativity of H compared to C. On the other hand, for the F, NO 2 , and NH 2 substituents, ΔQ σ shows an inductive effect where electronic charge flows from the phenyl ring to the more electronegative substituents. For the resonance in the π-electronic system (indicated by ΔQ π ), the effect of the substituents is different: the NH 2 and F substituents donate electronic charge to the π system of the phenyl ring, whereas the NO 2 substituent accepts electronic charge from the phenyl π system ( Figure 3). The H-substituent only has virtual π orbitals (e.g., 2p π ).
Therefore, for benzene, ΔQ π displays only the donation of a relatively small amount of phenyl π-electron density to the H substituent. The analysis of the VDD charge rearrangements ΔQ also gives insight into the regioselectivity in electrophilic aromatic substitution reactions caused by the substituents. Figure 3 shows that the π-charge flow (ΔQ π ) due to C N bond formation in nitrobenzene is responsible for the largest deactivation at the ortho (+35 milli-electrons) and para (+34 milli-electrons) arylic carbon atoms, making the NO 2 substituent meta directing for electrophilic aromatic substitution. On the other hand, the NH 2 substituent in aniline donates π-electronic density to  Figure 4A) and in turn can engage in more stabilizing electrostatic and orbital interactions when interacting with a hydrogen-bond acceptor. 9 In Figure 4B, the atomic charge shifts ΔQ following the covalent C-N bond (σ C N ) formation in the O-, S-, and Se-amides confirms a higher degree of charge depletion on the NH 2 group going from X = O (+177 milli-electrons) to X = S (+244 milli-electrons) to X = Se (+258 milli-electrons). The total charge change (ΔQ) can be partitioned, making use of the C s symmetry of the amides, into contributions of the A' and A" irreps, corresponding to charge shifts within the σ (ΔQ σ ) and π (ΔQ π ) electronic systems, respectively (see Figure 4B). The electron depletion on the NH 2 groups does not originate from the σ system, because ΔQ σ shows for each amide donation fragments enables one to identify the π-electronic system to cause the relative amide hydrogen-bond donor capabilities. Complementary Kohn-Sham molecular orbital analyses could reveal that the steric size of the chalcogen atoms is at the origin of this difference in πelectronic behavior. 9

| Decomposition into interaction energy components
As shown in the previous section, the analysis of symmetrydecomposed VDD charge changes ΔQ yields intuitive σand π-charge flow components. However, the VDD method can be even more informative when combined with the energy decomposition analysis (EDA) scheme, which allows not only quantifying the decomposition terms of the total interaction energy but also computing the charge shifts associated with these terms. Below, we will demonstrate the applicability of these complementary bonding-analysis techniques In the first example, we use the symmetry-decomposed VDD charge analysis and the EDA to study the nature of metal-CO coordination bonds in octahedral M(CO) 6 complexes with M = Fe 2+ , Mn + , Cr, V À , or Ti 2À . 22 As can be seen from Figure (Table S3).
The VDD analysis in Figure 5C shows that upon coordination to M(CO) 5  shows that the total ΔQ on the CO ligand follows from the relative magnitude of charge flow caused by σ donation (i.e., σ CO À! σ* M(CO)5 ) and π-back donation (i.e., π M(CO)5 À! π* CO ), indicated by ΔQ σ oi and ΔQ π oi , respectively, which is in line with the EDA results in Figure 5B. Note that the Pauli repulsion is not dictating ΔQ as the ΔQ Pauli term is relatively small and roughly equal for all complexes.
The application of the symmetry-decomposed VDD charge analysis in combination with the EDA is not limited to the example shown above but can be done for any type of intra-or intermolecular interaction. In the next example, we demonstrate the useful role that this analysis can play in identifying cooperativity and the origin thereof in hydrogen-bonded assemblies (self-assembly is cooperative if the stability of a polymer containing n monomers is substantially higher than (n-1) times the stability of the dimer). It has been established that cooperativity in hydrogen-bonded polymers stems from the charge separation caused by the σ-charge transfer interactions. 11,13c,d,23 The σ-charge separation leads to increasingly stabilized electrostatic and orbital interactions upon adding monomers, in which the latter follows from a reduction of the HOMO-LUMO gap between the σ orbitals involved in hydrogen bonding. For the identification of charge separation effects and cooperativity in supramolecular assemblies, a combination of the EDA and the VDD analysis can be very insightful as we will demonstrate below by examining the formation of a linear, hydrogen-bonded squaramide polymer.
In Table 1, the decomposition of the total interaction energy ΔE int is shown for the hydrogen-bond interaction between two squaramide monomers forming a dimer and between a squaramide monomer and F I G U R E 5 (A) M(CO) 6 complexes (M = Fe 2+ , Mn + , Cr, V À , Ti 2À ) with C O distances (in Å), (B) energy decomposition analysis (EDA) σ donation and πback donation components of ΔE oi between CO and M(CO) 5 , and (C) change in VDD charge ΔQ (in milli-electrons) of the CO ligand, decomposed into Pauli repulsion and orbital interaction components and into σ (A 1 irrep) and π (E irrep) components; computed at ZORA-BLYP/TZ2P in C 4v symmetry (see text).
dimer forming a trimer. The hydrogen-bond interaction strengthens substantially from the dimer (ΔE int = À15.4 kcal mol À1 ) to the trimer (ΔE int = À20.0 kcal mol À1 ) as a result of a more stabilizing electrostatic interaction ΔV elstat and σ-orbital interactions ΔE σ oi . The dispersion (ΔE disp ) and π-orbital interactions (ΔE π oi ) remain constant, whereas the Pauli repulsion (ΔE Pauli ) even becomes more destabilizing for the trimer, and hence are not responsible for the more stabilizing ΔE int .
T A B L E 1 Decomposition of the interaction energy ΔE int (in kcal mol À1 ) of the hydrogen-bond interaction (highlighted in red) in a squaramide dimer and trimer. a and from the σ-(ΔQ σ ) and π-electronic (ΔQ π ) systems ( Figure 6). The discrimination between σ and π electrons is possible because we examine here planar (C 2v symmetric) structures, where the irreps A 1 + B 2 and A 2 + B 1 correspond to the σ and π electrons, respectively. 20 The decomposition of ΔQ in Figure 6, reveals that the charge separation is a consequence of orbital interactions in the hydrogen bonds as ΔQ oi causes like ΔQ a net negative charge on the hydrogen-bond donating monomer and a net positive charge on the hydrogen-bondaccepting monomer. The Pauli repulsion is not responsible for this charge separation as the ΔQ Pauli term leads to charge shifts opposite to the total ΔQ. Now that we established that the charge separation within the squaramide dimer follows from orbital interactions, we decompose ΔQ oi into contributions of the σ-(ΔQ σ oi ) and π-(ΔQ π oi ) orbital interactions ( Figure 6). From this analysis, it becomes evident that the charge separation between the two hydrogen-bonded monomers, and thus the cooperativity (vide supra), originates from the σ-orbital interactions (indicated by ΔQ σ oi ) and not from the π-orbital interactions (indicated by ΔQ π oi ). The π electrons do, however, provide an additional component of stabilization to the hydrogen-bond interaction by polarization of the squaramide monomers by making the carbonylic oxygen atoms more negative and the N-H group more positive (see ΔQ π oi in Figure 6). This stabilization by π polarization is, however, much less important compared to the stabilization by the σ-charge transfer interactions (ΔE π oi À1.1 kcal mol À1 vs. ΔE σ oi = À7.5 kcal mol À1 in Table 1). We emphasize that the division of the charge changes due to orbital interactions (ΔQ oi ) into σ and π components reveals that the charge shifts on the hydrogen-bond front atoms caused by the σcharge transfer interactions (ΔQ σ oi ), are masked by the Pauli repulsion (ΔQ σ Pauli ) and π polarization (ΔQ π oi ) ( Figure 6) which can only be revealed by utilizing the symmetry-decomposition-analysis feature of the VDD method.
Not only the origin but also the consequence of the σ-charge separation in the growing squaramide polymer can be analyzed and understood by the VDD charge analysis. As we showed in Table 1, the self-assembly of squaramides is cooperative because of more stabilizing electrostatic and σ-orbital interactions as the polymer size increases. The more stabilizing σ-orbital interactions in the trimer (ΔE σ oi ¼ À9:5kcalmol À1 ) compared to the dimer (ΔE σ oi ¼ À7:5kcalmol À1 ) can be quantified by the symmetry-decomposed VDD charge analysis (Figure 7). Figure 7 shows that the stronger σ-charge transfer interaction (ΔE σ oi ) in the trimer yields indeed a larger net charge separation (64 milli-electrons) compared to F I G U R E 7 Change in Voronoi deformation density (VDD) atomic charges ΔQ (in milli-electrons) within the squaramide trimer caused by the hydrogen-bond interaction between a dimer and a monomer, decomposed into contributions stemming from Pauli repulsion and orbital interactions (ΔQ = ΔQ Pauli + ΔQ oi ), and into components of the σ (A 1 and B 2 irrep) and π (A 2 and B 1 irrep) electronic systems (ΔQ = ΔQ σ + ΔQ π ). All computed at ZORA-BLYP-D3(BJ)/TZ2P in C 2v symmetry. Color code of the ball-and-stick structures: H = white; C = gray; N = blue; O = red. the dimer (51 milli-electrons, Figure 6) and that this is a direct consequence of the larger flow of σ electrons (indicated by ΔQ σ oi in Figure 7) from the hydrogen-bond accepting monomer to the hydrogen-bond donating monomer.
In the previous examples, we have mainly shown how the symmetry-decomposition feature of the VDD charge analysis can provide intuitive σand π-charge flow components. However, also the charge flow of higher order bonds such as δ-bonds (and in theory also φ-bonds) can be studied using this method. In the final example, we show this for a textbook example of a system containing a quadruple bond: the octachlorodirhenate(III) dianion [Re 2 Cl 8 ] 2À . 24 In this complex, the rhenium atoms form with their d-type orbitals a quadruple covalent bond consisting of one σ-bond, two π-bonds, and one δbond (see Figure 8A). By using the pertinent point-group symmetry (in this case, C 4v ), 25 the VDD analysis using open-shell molecular fragments can quantify the charge flow in the σ (A 1 irrep), π (E irrep), and δ (B 2 irrep) 20 electronic systems that goes with the quadruple bond formation, as is visualized in Figure 8B. The charge flow in these irreps is primarily caused by the covalent bond formation but also involves some residual charge flow components caused by interactions between the ligands, which is partly covered in ΔQ rest , that is, the B 1 irrep. From  Here, we have presented a generalization of the symmetrydecomposed VDD analysis which can from now on be performed for any point-group symmetry (before only for Abelian point groups) of closed-shell, and also for open-shell (i.e., radical) sets of interacting molecular fragments. Furthermore, we demonstrate how the symmetry-decomposed VDD charge analysis augments the symmetry-decomposed energy decomposition analysis (EDA) so that the charge flow associated with Pauli repulsion and with orbital interactions can be quantified per atom and per irrep, for example, for σ, π, and δ electrons.
Altogether, the symmetry-decomposed VDD atomic charge analysis provides fundamental insights into chemical bonding from the perspective of the electronic density distribution. This augments the picture emerging from the corresponding energy decomposition analysis (EDA).

ACKNOWLEDGMENTS
The authors thank Dr. Erik van Lenthe for helpful discussions and software implementations, and the Netherlands Organization for Scientific Research (NWO) for financial support. This work was carried out on the Dutch national e-infrastructure with the support of SURF Cooperative.