ACS Publications. Most Trusted. Most Cited. Most Read
My Activity
Recently Viewed
You have not visited any articles yet, Please visit some articles to see contents here.
CONTENT TYPES

Figure 1Loading Img
RETURN TO ISSUEPREVA: New Tools and Met...A: New Tools and Methods in Experiment and TheoryNEXT

Interacting Quantum Atoms Method for Crystalline Solids

Cite this: J. Phys. Chem. A 2021, 125, 40, 9011–9025
Publication Date (Web):October 1, 2021
https://doi.org/10.1021/acs.jpca.1c06574
Copyright © 2021 The Authors. Published by American Chemical Society
ACS AuthorChoiceACS AuthorChoiceCC: Creative CommonsCC: Creative CommonsBY: Credit must be given to the creatorBY: Credit must be given to the creator
Article Views
750
Altmetric
-
Citations
-
LEARN ABOUT THESE METRICS
PDF (2 MB)
Supporting Info (1)»

Abstract

An implementation of the Interacting Quantum Atoms method for crystals is presented. It provides a real space energy decomposition of the energy of crystals in which all energy components are physically meaningful. The new package ChemInt enables one to compute intra-atomic and inter-atomic energies, as well as electron population measures used for quantitative description of chemical bonds in crystals. The implementation is tested and applied to characteristic molecular and crystalline systems with different types of bonding.

Introduction

ARTICLE SECTIONS
Jump To

Nowadays, many interesting materials, e.g., intermetallic compounds, belong to the families of chemical systems, where the usual valence approaches based on pure ionicity, covalency with 2-center 2-electron bonds, or “metallic bonding” are no longer valid. The interplay of all these features challenges our understanding with the emergence of new bonding scenarios beyond the well-accepted concepts. (1−6) Since chemical bonding analysis represents a way to understand inter-atomic interactions, it is also a way to understand material properties and their relations to specific electron counts or partial structures. For this sake, chemical bonding parameters like atomic charges and bond orders must be consistently related to their energetic counterparts.
This is achieved in the Interacting Quantum Atoms (IQA) method. (7−9) It represents an ideal framework to investigate covalent and ionic interactions on an equal footing. The sum of the individual atomic interaction energies forms a part of the recoverable total energy of the system. All routinely used quantum mechanical packages yield the total electronic energy. Usually this energy is given in terms of kinetic and potential energy contributions. Additionally, the potential energy is split into Coulomb and exchange-correlation parts. For deeper insight into the mutual interactions between the atoms, a suitable analysis of the total energy, respectively the wave function and its components, must be performed. A number of tools and methods are available for this task, but only the periodic energy decomposition analysis (pEDA) (10) has been used to decompose the total energy of extended systems. The decomposition is done in Hilbert space. Also the Crystal Orbital Hamilton Populations (COHP) (11) analysis does the decomposition in Hilbert space but decomposes the band structure energy instead. An alternative and appealing method to decompose the total energy components into contributions of and between position-space atomic regions is the approach of IQA. Unique features of IQA are that it is orbital invariant, methodologically independent of the type of basis set, applicable to correlated wave functions, and does not depend on an external reference state.
The IQA scheme is based on the idea that the one- and two-electron energy integrals computed over the total coordinate space can be evaluated separately over any set of nonoverlapping spatial domains. The connection to a sound chemical description is achieved by the utilization of domains defined by the Quantum Theory of Atoms in Molecules (QTAIM), i.e., by the atomic basins determined by the density gradient field. (12) With this, by successive integration of the one-electron components over each atomic basin, the total one-electron energy is decomposed into atomic contributions. The contributions describe the kinetic energy of the atomic domain as well as the attractive energy between the electrons and the nucleus enclosed within the atomic domain. The integration of the two-electron energy components over the atomic basins yields the intra-atomic and inter-atomic contributions that are used in the bonding analysis to judge the ionic and covalent bonding character. The whole procedure was adopted for molecular wave functions based on the Gauss-type of basis sets in the program Promolden. (13)
Until now, the IQA method has not been implemented for solids. Nevertheless, on the basis of careful calibration studies using a zeroth order approximation of the full IQA interaction energies, the so-called point-charge approximation of the covalent interaction energy and the ionic interaction energy (Madelung energy), was employed to understand chemical bonding and site occupation in ternary half-Heusler (MgAgAs-type) phases. (5) Finally, this scheme was employed to predict new phases with this kind of structure. (6) This may show the strength and scientific potential of the solid-state implementation of the exact methodology to be presented below.
The developed ChemInt software is coupled with the DGrid package. (14) This extends the capability to evaluate the molecular Gauss-type orbitals also to the utilization of Slater-type based molecular wave functions, as used in the ADF package, (15) as well as to access the numerical atomic orbitals expanding the wave functions for molecular and crystalline solid state systems computed with the FHI-aims code. (16) ChemInt is also capable to compute the delocalization indices for all those types of basis sets, including crystalline systems. (17) In particular, the bielectronic exchange energy integrals between position-space atomic regions for solid state wave functions are extremely challenging due to the huge number of participating orbitals, which made parallelization of ChemInt a necessary feature.
Theory and implementation of IQA for crystalline solids are outlined, reliability of the code is tested with an extensive set of molecules, and the method is applied to prototype solid-state compounds exhibiting different bonding scenarios.

Methods

ARTICLE SECTIONS
Jump To

The total electronic energy of a molecule or solid with pairwise Coulomb interactions comprises the kinetic energy of the electrons T, electron–nucleus interactions Ene, electron–electron interactions Eee, and nucleus–nucleus interactions Enn
(1)
Given access to the first order density matrix ρ1(r; r′) and the pair density ρ2(r, r′), the energy can be written as
(2)
Suppose that the space is partitioned into a set of domains that exhaustively recover the total volume and are mutually exclusive. Then, the total energy can be decomposed as
(3)
The above energy terms can be written in compact form as the sum of domain self-energies and interdomain energies,
(4)
Domain self-energies are
(5)
(6)
where α ∈ A stands for all nuclei enclosed in domain A, and
(7)
Interdomain energies are
(8)
(9)
and
(10)
Notice that intradomain electron–electron terms are halved to represent total energies whereas interdomain terms represent nonequivalent interaction energies.
The approach to the energy decomposition given by eq 4 considers the crystal to be formed by interacting atoms, molecules, etc. One may, otherwise, take a reference unit formed by one or more domains, embedded in the crystal. We collect, in an additive manner, E = GE(G), all energy contributions E(G) of reference unit G,
(11)
The interaction energy for reference unit G shall be organized by coordination spheres i,
(12)
where there are number of B neighbors in the coordination sphere i of domain A of the reference unit. The last term condenses the sum over the coordination spheres for each domain in the reference unit into a single sum over coordination shells j for all domains in the reference unit and, with interaction contributions within that shell. The multiplicity factor for the shell of those interactions (i.e., Li–Cl(1) or Cl–Li(1)) takes the form
(13)
where
(14)
and nA, nB are the numbers of symmetry equivalent atoms A and B in the reference unit, respectively.
The pair density can be viewed as the product of two quasi-independent densities minus a term accounting for the difference with the true pair distribution ρ2(r, r′) = ρ(r)ρ(r′) – ρxc(r, r′). (18) This leads in eq 8 to the Coulomb energy contribution
(15)
and the exchange-correlation energy contribution
(16)
With the exception of the exchange-correlation energy, all interdomain terms are classical. Thus, the interaction energy can be arranged as a sum of a classical electrostatic energy EclAB and a covalent energy ExcAB,
(17)
In the special case of the monodeterminantal wave function, the exchange-correlation density is directly given by the first order density matrix, ρxc = ρx = |ρ1(r; r′)|2.

Atomic Boundaries

Atomic domains defined by QTAIM theory are delimited by surfaces S that satisfy the zero-flux condition
(18)
where ρ(r) is the density at point r in the surface and n(r) is a normal vector to the surface. The manifold of all steepest ascent paths
(19)
terminating at given attractor (the common destination limtr(t)) constitutes an atomic (QTAIM) basin. In the following, the QTAIM basins are used as the spatial domains over which the energetical decomposition is performed.
We note at this point that a representation of the surface compatible with atomic-centered grids is required to use the multi/bipolar expansions presented in the following sections. To alleviate the expense of the computation, we precomputed a trust sphere using the algorithm described by Rodríguez et al. (19)
Of fundamental relevance for solid state systems is achieving electroneutrality in the unit cell. In contrast to molecules, all atomic basins determined for a solid state system have finite volume. Thus, the atomic basins must recover the unit cell volume, a circumstance that seems to be problematic in case of atom-centered radial grids. We decided to use radial rays with symmetric spherical t-designs. (20)

Energy Integrals in the IQA Framework

The energy integrals for a system with perfect periodicity are performed over complex crystal orbitals ϕn,k for given band index n and reciprocal vector k
with R given by the unit cell vectors and the linear combination ψn,k(rR) = jcn,j(k) χj(raj) of atomic orbitals χj(raj), centered at aj = rj + R, and with coeficients cn,j(k) to be determined. In this paper, electronic structure calculations were done with FHI-aims, utilizing numerical atomic orbitals.
The band and translational symmetry labels are condensed here to a single label i ≡ (n, k). With this, the crystal orbitals ϕi determine the electron density as well as exchange pair density , where are overlap densities.
The density is replaced in eq 15, and the exchange pair density in eq 16,
(20)
with domains being now QTAIM atoms. Recall from eq 7 that B = A integrals have a prefactor.
Intra-atomic Coulomb integrals (A = B) are computed with the Laplace expansion of 1/|rr′| around the position of the atomic nucleus (7)
(21)
where is the angular ( = (θ,ϕ)) integral of the atomic density
(22)
weighted by real spherical harmonic Sl,m. The discriminant , with r< =  min(r, r′) and r> =  max(r, r′), ensures the convergence of the expansion.
The expansion is valid also for the exchange integral (7)
(23)
where . Atomic overlap densities ρi,jA are zero outside of the basin like the atomic density above.
The integral evaluation is more efficiently performed in the basis of overlap densities . The exchange term is diagonal in this basis with coefficients dij = 2(2 – δij) for monodeterminantal wave functions, (9) and no monadic diagonalization is required.
Two-center integrals (AB) are expanded around two poles, one at each atomic position, with a bipolar expansion (21−23)
(24)
where replaces as the discriminant. (7)
The exchange integral between two different basins is
(25)
Once the electronic structure problem is solved by whatever method, the self-energy of each atom (or fragment) in the unit cell is obtained. The interaction energy of each of them, e.g., A, with the rest of atoms of the lattice, appropriately classified by increasing distance, B, is obtained (). This infinite sum contains quickly convergent terms ExAB, which decay exponentially in insulators and polynomially in metals, and a potentially divergent series made up of EclAB electrostatic interactions. By separating Ecl into a short-range, penetration-like exponentially decaying term and a multipolar one that is summed via the Ewald construction, one gets rid of divergences. We have applied the Ewald construction with the QTAIM monopoles (net atomic charges) via the Environ program (24) and have summed directly the higher order multipolar corrections from the IQA integrations up to near neighbors. Even if the QTAIM monopoles are zero, those terms remain and constitute the well-known electrostatic interaction between uncharged atoms or molecules.

Approximations Based on Population Measures

Inter-atomic bielectronic integrals are a particularly expensive step. As seen in eqs 24 and 25, they involve a bipolar expansion. For classic electrostatic energy integrals of distant basins it is equivalent to use a multipolar expansion (25)
(26)
where
(27)
is a real multipole moment of the net charge density ρt(r) = AZAδ(rRA) – ρ(r) in basin A. is an angular factor (cf. eq 25 in ref (7)). Technically, this approximation is invalid if the circumsphere of the basin A (centered at the nuclear position) intersects the circumsphere of basin B (with center at nucleus B). As a remark, solid state systems normally have compact basins. This limits the extent of the spheres, and it is safe to say that classic interactions with neighbor atoms of the second and higher coordination spheres are well approximated with a multipolar expansion in these cases. Also, internuclear distances in the solid state systems are often large, exceeding 2 Å.
On the other hand, exchange integrals were proven numerically to lead to asymptotic convergent multipolar expansions, with low-order terms being a good approximation for second and further nearest neighbors. (26) The appropriate expression has the same structure as for classic interactions
(28)
but in terms of the sum of overlap density multipoles
(29)
Further approximation of bielectronic integrals is to consider only the first term in the multipolar expansion. For the classic electrostatic term it is the same as assuming spherical density distribution in the atomic domain. By Gauss law, a point charge located at the nucleus position (for A and B) produces a physically equivalent classic force. The classic electrostatic energy,
(30)
is proportional to the product of QTAIM net charges for both basins, where is the average electron population in basin A.
Inter-atomic exchange integrals also admit a similar expression (26)
(31)
that instead involves the delocalization index (DI) δAB between the atomic domains A and B
(32)
The delocalization index measures the number of electrons shared between domains (A and B) and is complemented by the localization index (LI)
(33)
Localization and delocalization indices must obey the following sum rule
(34)
The fluctuation of electron population σ2(A) in basin A is the result of electron pair sharing with atoms of the first coordination shell σ12(A), second coordination shell σ22(A), and so on.
Expressions 30 and 31 are the basis for our more generic scaled point-charge approximation (sPCA) for classic ionic and exchange (covalent) inter-atomic interactions. The scaling parameters scl and sx approximate the exact classic and exchange energy, respectively, with a scaled zeroth order multipolar approximation:
(35)

Computational Details

ARTICLE SECTIONS
Jump To

Density functional theory calculations (PBE functional (27)) for crystalline solids were performed with FHI-aims, an all-electron full-potential electronic-structure code. The numeric atom-centered basis sets used were of standard “tight” type. The computed unit cells and Brillouin zone samplings are listed in Table 1. An external output option of FHI-aims (version 200112) is used to generate DGrid compatible wave function information. ChemInt interfaces DGrid to evaluate wave function properties at discrete points in space, determines the atomic QTAIM basins, and performs the IQA analysis with this space partitioning.
Table 1. Crystallographic Information and Brillouin Zone Sampling for the Calculated Materialsa
systemspace group typelattice parameters [Å]calculated cellk-point mesh
β-N2bP63/mmca = 4.050c1 × 1 × 15 × 5 × 3
  c = 6.604  
CO2Paa = 5.4942 (35)1 × 1 × 14 × 4 × 4
diamondFdma = 3.566606 (36)1 × 1 × 14 × 4 × 4
BN (zincblende)F4̅3ma = 3.61 (37)1 × 1 × 14 × 4 × 4
graphiteP63/mmca = 2.464 (38)2 × 2 × 14 × 4 × 3
  c = 6.711  
BN-bP63/mmca = 2.504323 (39)2 × 2 × 14 × 4 × 3
  c = 6.658852  
MgB2P6/mmma = 3.0846 (40)2 × 2 × 13 × 3 × 5
  c = 3.5199  
LiClFmma = 5.12952 (41)1 × 1 × 14 × 4 × 4
NaClFmma = 5.7915 (42)1 × 1 × 14 × 4 × 4
MgOFmma = 4.213 (43)1 × 1 × 14 × 4 × 4
AlFmma = 4.0494 (38)1 × 1 × 14 × 4 × 4
NaIm3̅ma = 4.235 (44)1 × 1 × 14 × 4 × 4
a

Lattice parameters are given for the crystallographic cell, and calculated cells are specified as multiples of the crystallographic cell.

b

The z coordinate for the Wyckoff site is calculated assuming a bond length of 1.108 Å and that the molecule center is located at the symmetry center of the space group P63/mmc, as in ref (45).

c

Lattice parameters taken from ref (46).

The set of molecules taken for the validation were computed with GAMESS, (28) ADF, (15) and FHI-aims, (16) employing HF, LDA, (29−31) PBE, (27) BLYP, (32,33) and B3LYP (34) functionals. The basis set taken within each program is cc-pvtz (GAMESS), TZ2P (ADF), and standard “tight” (FHI-aims). For each molecule, the geometry was optimized before integrating the energy components. The geometry was fixed for molecules used for comparison against the solid phase (N2, CO2, neopentane, and phenalene). Only the four central atoms of phenalene were fixed while the positions of the other atoms were optimized.

Results and Discussion

ARTICLE SECTIONS
Jump To

IQA energy terms are examined for a number of molecular systems to evaluate the energy error. In subsequent exemplary application studies, some prototype compounds (diamond and zincblende structures, honeycomb networks, rocksalt structures, closest packings) covering covalent, ionic, and metallic situations were investigated.

Validation of the Implementation

Three approaches are explored to validate our implementation. First, similar results should be obtained with alternative implementations when comparable (molecular wave functions expanded with GTOs). In this document only the ethane system is shown (Table 2), but the equivalence was tested also with other simple systems reaching always similar accuracy. Second, properties like the volume, exchange-correlation energy, and total energy in the cell were checked against FHI-aims. The total energy requires an exhaustive computation of long- range interactions to be recovered. The electron count sum for the (crystallographic) unit cell provide a third way to check the validity of our implementation. Only the first is examined in this section and the other two in the next sections.
Table 2. Total Energy Components of Ethane Computed with ChemInt, Promolden, and GAMESSa
energyChemIntPromoldenGAMESS
total energy, E–79.731–79.732–79.730
kinetic energy, T79.24379.243–79.244
total potential energy, Ene–158.974–158.975–158.974
electron–electron energy, Eee67.55167.55267.551
Coulomb energy, ECoul80.52880.52980.528
exchange-correlation energy, Exc–12.977–12.977–12.977
a

The evaluated number of electrons per formula unit are 17.9992 (ChemInt) and 18.0001 (Promolden). Ex was re-scaled from the integration of the PBE functional. (47) Energy in Ha units.

For the first check, we take the same molecular wave function and perform an energy decomposition analysis with two independent IQA implementations: our code and the Promolden code. (13) We see that the integrated quantities have similar accuracy to those of Promolden and are consistent, beyond the chemical accuracy of approximately 1 mHa, with accurate full space integrations given by GAMESS (Table 2).
The precision of IQA decomposition is mainly determined by errors in the atomic boundary, the quality of the grid, and the truncation of the Laplace/bipolar expansion. Another possible source of error can be differences in the evaluation of the energy in ChemInt with respect to the SCF program. The energy difference between reconstructing the total energy from IQA components against the SCF energy is a measure of the error of our method. For the evaluated molecular systems, this energy error per atom is always below chemical accuracy (see Figure 1) with a modest choice of integration parameters. The boundary is determined with a precision of more than 5 × 10–4 Å normal to the surface using an ODE integrator that preserves a local relative error of the same magnitude. The grid has 600 radial points and 5780 points distributed on a symmetric spherical t-design grid. A maximum expansion order of l = 4 (l = 6) for regions near (far from) the nucleus is chosen. Increasing the precision of the surface determination does not yield a better integration. As well, the length of the expansion chosen here is enough to approximate bielectronic integrals. From our observations, the quality of the grid is the most determining factor once the rest of the parameters have been fixed like here.

Figure 1

Figure 1. Energy difference (per atom) between the reconstructed IQA energy and the SCF energy (energy error). The white band denotes the median. Boxes delimit the lower and upper quartiles Q3 and Q4, and fences delimit all energy errors from a sample that takes wave functions computed with GAMESS, ADF, or FHI-aims, employing HF, LDA, PBE, BLYP, and B3LYP functionals for every system. Ex was rescaled from the integration of the corresponding functional. (47) All errors fall below 1 mHa.

Nevertheless, it is important to realize that periodic systems require a significantly larger number of orbitals than molecular systems. Computations with maximum resolution are not always feasible, and one must find a balance of precision and speed. Therefore, we are interested in the minimum required computation needed to achieve a reasonable accuracy level. For all the systems tested, we found that a surface represented with approximately 6000 points is enough to integrate the total volume in the crystallographic unit cell with a volume error ΔVcell = AVAVcell below 0.200 Å3 if we combine it with a radial grid of at least 200 points (Table 3). The symmetrical t-design quadrature offers a more robust integration, devoid of occasional (large) volume integration errors present with Lebedev–Laikov grids. In particular, for half-Heusler LiMgN, the cell volume error would be 2.027 Å3 using a 5810-point Lebedev grid. Instead, it is 0.010 Å3 with a 5780-point t-design grid. The distribution of points offered by the new grids is more homogeneous, thus providing a better overall representation of basin surfaces. Besides this, volume errors take place at the boundary of the basins so the integration of the charge can, and this is the case from our observations, still be within chemically reasonable accuracy error (<0.01).
Table 3. Reconstructing the Crystallographic Unit Cell Volume from QTAIM Basin Volumes: Volume Error, ΔVcell, and Cell Charge, Qcell
compoundstructure typeΔVcell3]Qcell [e]
LiClNaCl0.104–0.005
NaClNaCl0.1680.009
MgONaCl0.0390.002
CsClCsCl0.004–0.000
Cdiamond–0.0040.007
BNzincbende0.0200.000
Cgraphite0.0170.001
BNBN-b0.0090.001
Naα-W (bcc)–0.0280.001
AlCu (fcc)0.0100.000
MgB2AlB20.0140.002
LiMgNMgAgAs0.0100.006
N2β-N2–0.0070.008
CO2CO2–0.0730.000

Molecular Crystals

Molecular crystals are a suitable example to examine the type of strong bonds that one finds in molecules before shifting to the highly diverse bonds present in crystalline materials. For that matter, we analyze β-N2, and CO2 crystals. A priori, N–N and C–O bonds are only expected to differ slightly from a gas phase molecule. Indeed, covalent intramolecular energies in the solid phase bear a striking resemblance to bonds of equal internuclear distance between atoms in gas phase, as shown in Table 4. Intermolecular interactions in the solid barely debilitate the covalent stabilization provided by the N2 triple bonds. As well, the classic interaction EclNN between two nearest nitrogen atoms is almost unaltered. In this regard, the triple bond exhibited by a N2 molecule is transferable to a triple bond in β-N2.
Table 4. Comparison of IQA Two-Center Terms for Covalent Bonds in Molecular and Solid Phases β-N2 and CO2a
systemABRABδABEnnABEneABEneBAECoulABEclABExAB
N2 (mol.)N–N1.1083.04123.399–21.944–21.94420.7160.227–0.906
N2 (solid)N–N(1)1.1083.00023.399–21.952–21.95320.7260.220–0.903
CO2 (mol.)C–O1.1491.39622.114–25.443–13.41315.453–1.289–0.444
CO2 (solid)C–O(1)1.1491.32322.114–25.646–13.15215.269–1.415–0.426
CO2 (mol.)O–O2.2970.43314.743–16.981–16.98119.5630.343–0.054
CO2 (solid)O–O(2)2.2970.40714.743–17.081–17.08119.7930.373–0.052
a

Nearest-neighbor interactions are denoted with AB(1), second neighbors with AB(2), and so on. Energies in Ha; distances in Å.

On the other hand, while the C–O bond is similar in both environments, the strong polarization of this bond permits ionic interactions with neighbor molecules, with O–O being the most intense interaction. The C–O bond is mainly altered by its classic ionic component, which varies by 0.126 Ha in comparison to 0.018 Ha for the covalent energy ExCO. In both scenarios, a strong contribution from dipolar and quadrupolar components to the classic energy is found. The same effect is observed for intramolecular O–O interactions but is more attenuated.

Diamond and Zincblende Structure

Convergence of Bielectronic Integrals

A steady convergence with increasing multipolar order is observed for intra-atomic bielectronic integrals in diamond and zincblende BN (Figures 2 and 3). The Coulomb energy ECoulAA shows a dominating monopolar term (l = 0) followed by zero l = 1, 2 terms as demanded by symmetry. (48) The l = 3, 4 terms are nonzero due to nonspherical distribution of electrons in the basin. The l = 5 term vanishes again due to symmetry. These observations for diamond atoms are in close agreement with those of a carbon atom in methane or the quaternary atom in neopentane (cf. Supporting Information, Figure S2 and Figure S3).

Figure 2

Figure 2. Diamond: Convergence of bielectronic integrals inside the basin of the carbon atom, with increasing multipolar order l, and . Only their magnitude is plotted. Dashed lines indicate the trend of convergence for symmetry allowed terms. Disconnected dots are not allowed by symmetry and are nonzero due to numerical errors. Terms below 10–5 Ha are represented as dots at the bottom.

Figure 3

Figure 3. BN (zincblende): Convergence of bielectronic intrabasin integrals with increasing multipolar order l. For point and lines explanations, see Figure 2.

Convergence with l of nearest-neighbor inter-atomic electron–electron terms in diamond is slower than for intra-atomic terms (see Figures 2 and 4). The convergence of Coulomb and exchange components with the order of the multipolar expansion was examined to assess the possibility of truncating the series. Storage of higher-order terms contributes to a substantial increase of computation time, specially due to the increasing number of averaged overlap densities Rl,mi,j. Coulomb and exchange terms show the expected convergent behavior. However, terms of order l = 3, 4 are still important for the Coulomb integral. This forces us to take several orders in the expansion. This behavior could, as well, vary for other systems. Therefore, it is only a preliminary observation.

Figure 4

Figure 4. Diamond: Convergence of bielectronic interbasin integrals with increasing bipolar order L = lA + lB. Labels in the top right indicate an interaction of atom A with an atom B of the ith coordination sphere as AB(i).

The Coulomb integral ECoulAB between a boron atom and a nearest N atom is well approximated from the (lA, lB) = (0, 0) term (Figure 5). Higher order terms L = lA + lB > 0 clearly reflect the characteristics of the electron distribution in their respective basins. Terms with order L = 1, 2 result from the combinations {(0, 1), (1, 0), (1, 1), (2, 0), (0, 2)}. As we have seen above, since the symmetry in B and N sites is tetrahedral, l = 1, 2 terms are zero (Figure 3). The same argument justifies the absence of L = 5 terms. Thus, only with L = 0, 3, 4, and 6 are the corrresponding bicentric Coulomb terms nonzero. Only due to numerical accuracy are they nonzero. More distant interactions like the closest N–N and B–B exhibit the same pattern, but their contribution from higher L order terms is already small, especially for B–B interactions.

Figure 5

Figure 5. BN (zincblende): Convergence of bielectronic interbasin integrals with increasing bipolar order L = lA + lB. For point and lines explanations, see Figure 2.

Self-Energies and Near Neighbor Interactions

Having a generic method to decompose atomic energies allows us to compare molecular and solid systems on the same footing. Diamond is a prototype of a covalent network with bonded quaternary carbon atoms. A carbon atom with the same coordination in a molecule is the quaternary carbon of neopentane.
Table 5 compares the population and intra-atomic energy components of this atom in the two environments. The average electron population of the quaternary C atom in neopentane is slightly lower than that in diamond carbon atoms. Both intra-atomic classic EclCC and exchange energies ExCC are larger in diamond due to the higher population of the atomic basins.
Table 5. IQA Monocentric Integrals Inside Basins with Td Point Symmetrya
systemANeAλAσ2(A)TAEneAAECoulAAEclAAExAA
neopentaneC5.9293.7752.15437.781–89.69419.154–70.540–4.566
diamondC5.9993.8202.18037.917–90.12819.522–70.606–4.602
BNB2.8392.0750.76523.700–51.7537.923–43.830–3.008
 N9.1617.5171.64355.650–138.16236.638–101.524–7.019
a

Energies in Ha units.

Carbon atoms have a similar inter-atomic distance in diamond and neopentane, approximately 1.545 Å. To facilitate the discussion, the C – C distance in neopentane is fixed to be exactly the same as in diamond. This way, the nucleus–nucleus repulsion is identical. Full relaxation of the geometry does not change our conclusions. Inter-atomic interaction energies are presented in Table 6. As for intra-atomic terms ECoulCC, electron–electron Coulomb energies are larger in the solid phase. This is counterbalanced by EneCC to finally yield a smaller classic energy. In diamond, electron–nucleus interactions are symmetrical, EneAB = EneBA. The sum of all classic electrostatic terms nearly cancels to yield EclCC = 0.019 Ha for neopentane and EclCC = 0.014 Ha for diamond because the electroneutrality of the atom makes the L = 0 term vanish. The classical L > 0 energy term Ecl,L≠0AB has to be positive for neutral nonoverlapping charge distributions. (8)
Table 6. IQA Two-Center Terms for Relevant Interactions in Neopentane, Diamond, and Zincblende BNa
systemABmRABδABEnnABEneABEneBAECoulABEclABExAB
neopentaneC–C(1)41.5440.95712.335–12.058–12.09411.8360.019–0.288
diamondC–C(1)41.5440.91412.335–12.217–12.21712.1130.014–0.284
 C–C(2)122.5220.0427.554–7.551–7.5517.5490.001–0.006
BNB–N(1)41.5630.35711.848–15.502–6.6978.765–1.586–0.106
 N–N(2)62.5530.15210.158–13.290–13.29117.3900.967–0.030
 B–B(2)62.5530.0025.183–2.943–2.9431.6710.968–0.000
a

is the number of equivalent interactions AB per reference unit (G = C1C2 or G = B1N1) (eq 13). Energies in Ha; distances in Å units.

Exchange energies ExCC are almost identical in neopentane and diamond. Conversely, the nearest-neighbor DI δCC in diamond is smaller due to higher delocalization of the wave function (more distant neighbors) in the solid.
Diamond presents a DI of δCC = 0.914 between nearest-neighbor atoms. The delocalization goes down to δCC = 0.044 for second nearest-neighbors. The variation of DIs for the third coordination shell obtained for equivalent interactions is too large to prevent a detailed numerical examination. Up to first and second nearest neighbors the sum rule (34) is shown in Table 7.
Table 7. Cumulative Electron Populations of Atomic Basins in Diamond and Boron Nitride (Zincblende) from Successive Inclusion of Delocalization Shells α with DI δαa
αcoordination sphereδαNelAcumα
0C3.820
1C–C(1)0.9145.648
2C–C(2)0.0445.912
0B2.075
1B–N(1)0.3572.789
2B–B(2)0.0022.801
0N7.517
1N–B(1)0.3578.231
2N–N(2)0.1559.161
a

The sum rule of eq 34 approaches ⟨NelA⟩=6 as delocalization with more distant atoms is included. Boron has ⟨NelA⟩ = 2.839 electrons, and nitrogen has ⟨NelA⟩ = 9.161 electrons.

Boron nitride represents an example of a diamond-type structure with polar bonds. While the nearest-neighbor electrostatic interaction EclCC=0.014 Ha in diamond represents a destabilization, in B–N it is a stabilizing classic interaction between opposite charged atoms with an energy contribution of EclBN=–1.586 Ha (Table 6).
The nearest-neighbor exchange interaction ExBN is only one-third of the analogous nearest neighbors interaction ExCC in diamond. Note that this decrease can be estimated (31) from the strong decrease of the delocalization index from δCC = 0.957 to δBN = 0.357.
For BN, in overall, classic electrostatic energies between nearest neighbors are much larger than the exchange ones. This fact reveals the high importance of classic interactions in the short-range regime.
The closest B–N interactions involve 0.36 electron pair sharing (Table 7) that parallels the value of 0.37 obtained before. (49) Nitrogen atoms share 0.15 electron pairs with any other nearest nitrogen atom (to be compared with 0.35 electron pairs). The sum-rule for N is approached including interactions up to the second coordination shell (see Table 7). On the other hand, boron atoms have no significant electron sharing with other boron atoms. There are, however, 0.038 electrons still missing to recover the average number of electrons in the boron basin.
The localization indices are in agreement with results obtained before by some of us (Table 6).

Additive Interaction Energies

The additive interaction energy (cf. equation 11) for a 2-atom reference unit G = {C1, C2} ≡ C1C2 in diamond can be approximately obtained as follows. The covalent interaction energies between nearest neighbors (−0.284 × 4 = −1.136 Ha) and between second nearest neighbors (− 0.006 × 12 = −0.072 Ha) sum to a total covalent interaction energy of Ex(C1C2) = −1.208 Ha. Higher neighbors are omitted due to a very small value of remaining bond fluctuations . A kind of lower bound of the neglected covalent interaction can be estimated with the zeroth order approximation if we assign all remaining electrons to the third coordination shell, Ha. The stabilizing covalent interactions are counterbalanced by electrostatic interactions Ecl(C1C2) = 0.014 × 4 + 0.001 × 12 = 0.068 Ha (in diamond) from contributions between first and second nearest neighbors, respectively. The total interactions considered amount to Eint(C1C2) = −1.140 Ha.
For zincblende, again, a diatomic reference unit G = {B1, N1} ≡ B1N1 is chosen. The approximate covalent interaction energy up to second nearest neighbor terms sums nearest neighbor B–N(1) interactions (− 0.106 × 4 = −0.424 Ha), second nearest neighbors N–N(2) interactions (− 0.030 × 6 = −0.180 Ha), and second nearest neighbor B–B(2) interactions (0.000 × 6 = 0 Ha) giving Ex(B1N1) = −0.604 Ha. While the electron sum rule is already well satisfied for N, some 0.038 electrons are missing on B. If we use this number to approximate the maximum additional covalent energy, the third nearest neighbor B–N(3) interaction with a distance may be estimated to contribute maximally . It is interesting that a non-negligible covalent N–N(2) interaction has been found, which makes up 30% of the covalent bond energy. It can be considered as the residual bonding of the incompletely occupied N2.161– atoms, which is a necessary consequence of the polar bonding B–N(1) leading to the small DI δBN = 0.357.
The procedure applied for the covalent part is not feasible for the electrostatic part, because of the long-range of the electrostatic interactions and the well-known resulting nonconvergence of the simple coordination series energy in real space. The electrostatic part of the interaction per reference unit can only be approximated by the point charge approximation using the Madelung constant of zincblende fM,R = 1.63806. The validity of the PCA for this procedure can be estimated by determining the PCA scale factor scl for the computed first and second nearest neighbor Ecl(B1N1) energies using Q atomic net charges and the corresponding distances , and . It turns out that the PCA energy values are rather close to the exact IQA ones (, ), such that the PCA is well justified. The Madelung energy per unit cell G = B1N1 is calculated according to Ha, where R is the first-neighbor B–N distance. Although l = 1, 2 (dipoles, quadrupoles) are forbidden by symmetry, there are contributions from classic electrostatic terms with L > 0 due to the large atomic net charges present in the system. Owing to their absence in the Madelung formula, they must be accounted at least for nearest neighbors. We compute them with Ecl(G) – Ecl,L=0(G) = −0.026 Ha, where the last term is just the point-charge contribution (eq 30). The total electrostatic energy is then Ecl = −2.617 Ha. From these values, the total interaction energy amounts to Eint(B1N1) = −3.221 Ha. The fraction of covalent bond energy per G = B1N1 with respect to the total interaction energy yields a 19% contribution of covalent interaction. Based on this fraction, the interaction in BN (zincblende) is characterized by dominating 81% ionic interactions. Still, the covalency is non-negligible, and it represents the local type of bonding, while the electrostatic Coulomb interaction is a collective type of bonding as discussed previously. (50) While the covalent bonding A–B is a local property of A and B, the electrostatic interaction A–B has long-range consequences, alternating stabilizing and destabilizing terms. These two different types of bonding cannot be discussed on the same footing, i.e., in a local picture. For this reason, the energies per reference unit have been calculated for characterization.

Honeycomb Networks

Honeycomb networks can be found in graphite, hexagonal boron nitride, and the intermetallic superconductor MgB2. Table 8 shows nearest-neighbor interactions in the three hexagonal systems. In addition to the σ bonds discussed in the previous section, they might exhibit a certain degree of π double bond character. In graphite, it has been shown, (51) based on delocalization indices, that the shortest C–C bonds with a DI δCC = 1.20 have partial double bond character, whereas in MgB2 the number of shared electron pairs is the same as for a single bond δBB = 1.00. For comparison, similar values of δCC = 1.202 in graphite and δBB = 0.983 in MgB2 are obtained.
Table 8. IQA Two-Center Terms for Relevant Interactions in Hexagonal Graphite, BN, and MgB2a
ABmRABδABEnnABEneABEneBAECoulABEclABExAB
C–C(1)ab31.4231.20213.391–13.098–13.10212.8460.037–0.372
C–C(2)ab62.4640.0547.731–7.692–7.6927.6550.002–0.007
C–C(3)ab32.8450.0376.696–6.679–6.6816.6650.001–0.004
C–C(4)c13.3560.0195.677–5.709–5.7095.7460.004–0.003
C–C(5)c6 + 33.6450.0065.227–5.245–5.2475.2660.001–0.001
B–N(1)ab31.4460.45212.810–16.751–7.0659.244–1.762–0.135
B–B(2)ab32.5040.0045.283–2.936–2.9351.6321.043–0.000
N–N(2)ab32.5040.21210.354–13.604–13.60517.8771.021–0.040
B–N(3)ab32.8920.0086.405–8.419–3.5644.685–0.893–0.001
B–N(4)c23.3290.0065.563–7.329–3.1104.101–0.776–0.001
N–N(5)c33.6300.0247.144–9.411–9.41012.4030.726–0.003
B–B(1)ab31.7810.9837.429–8.392–8.3929.5220.166–0.249
B–B(3)ab63.0850.0524.289–4.953–4.9535.7220.106–0.006
B–B(5)ab33.5620.0173.714–4.303–4.3034.9870.095–0.002
Mg–B(2)c122.5040.06112.681–14.830–10.96412.824–0.290–0.013
Mg–Mg(3)ab33.0850.00324.704–21.364–21.36418.4750.450–0.000
a

The reference units corresponding to are G = C1C2, G = B1N1, and G = MgB1B2. Energies in Ha; distances in Å units. Key: ab, intralayer interaction; c, interlayer interaction

Graphite nearest-neighbor C–C bonds display higher ExAB energies (covalent bond energies) than in the other systems with honeycomb networks. The nearest-neighbor B–N covalent bond is particularly weak. In the honeycomb layer, the Ex proportion of the B–N covalent bond compared to the C–C covalent bond, , parallels the proportion seen in diamond and zincblende BN, . Conversely, classic electrostatic interactions are different for each system studied here. BN is largely stabilized from collective classical interactions between ions. This finding provides conclusive support for interpreting B–N as a very polar bond. However, destabilizing second-neighbor ionic interactions overcome the stabilization achieved by nearest neighbors. In the long-distance limit, those classic ionic interactions would tend to partially cancel. (50) Graphite only has high order multipole contributions to the classic electrostatic energy. Thus, already for first neighbors, a low value of Ecl = 0.037 Ha is obtained. In MgB2, nearest B–B interactions display a considerable classic (interionic) destabilization that is compensated by Mg–B classic electrostatic stabilization.
As in zincblende BN (ExNN = −0.030 Ha), among second-neighbor interactions only N–N still has a notable exchange contribution ExNN = −0.040 Ha in the hexagonal phase.

Phase Stability: Cubic versus Hexagonal Structures

Diamond and graphite are the main allotropes of carbon. The experimental difference in enthalpies is known to be quite small, ΔH298 K = +1.895 kJ/mol, (52) requiring an exquisite degree of accuracy that challenges DFT functionals. Albeit the relative stability problem is already solved, (53) the internal forces leading to a greater stability of graphite are not completely known. The additional stability might result from the delocalized π bonding framework. However, it is not clear if the higher coordinated carbon atoms (via σ bonds) in diamond should be less stable. The origin of graphite stability could as well be attributed to a particular local property like a different hybridization. In the spirit of a recent publication that includes those contributions in an analytical model, (54) we perform a numerical exploration to elucidate why nature chooses graphite over diamond.
Similar to the case of diamond discussed in the Diamond and Zincblende Structure section, the case of graphite involves a combination of stabilizing covalent bond contributions – 0.004 × 3 = −1.170 Ha for ortho, meta, and para neighbors for the chosen C1C2 reference unit. The covalent part of the C–C interaction in the c direction amounts to = Ha to be added to yield final Ex(C1C2) = −1.182 Ha. The destabilizing electrostatic interactions amount to = = +0.126 Ha (intralayer) and Ha (interlayer) summing up to Ecl(C1C2) = +0.139 Ha. A total interaction energy of Eint(C1C2) = −1.043 Ha is obtained. Comparison of the total interaction energy to that of the diamond modification Eint(C1C2) = −1.140 Ha (Ex(C1C2) = −1.208 Ha, Ecl(C1C2) = +0.068 Ha) reveals that the covalent interactions up to third nearest neighbors are energetically smaller, and the destabilizing electrostatic ones are larger. This yields an interaction energy preference for the diamond structure of ΔEint(C1C2) = 0.097 Ha.
Thus, the interaction energies do not explain the obtained total energy difference ΔEFHI-aims(C1C2) = 0.009 Ha between both modifications giving a preference to the graphite modification. This finding suggests that, at least for the current calculation, the preference of the graphite modification is caused by intra-atomic electron reorganization terms. Indeed, inclusion of the intra-atomic energies using ΔE(C1C2) = E(C1C2, diamond) – E(C1C2, graphite) = ΔEself(C1C2) + ΔEint(C1C2) resolves this seeming discrepancy and yields the energetical preference of the hexagonal phases (Table 9). In the framework of Valence Bond Theory, one would argue, the sp3 hybridization of carbon in diamond requires more energy than the sp2 hybridization in graphite. Just as in the solid phase, the 4-connected central carbon atom in neopentane and the 3-connected one in phenalene display a consistent self-energy difference of ΔEself(2 × C1) = 0.068 Ha, thereby demonstrating that the hybridization of the carbon atom is related to the stability of graphite over diamond. This indicates that the relative stability of carbon allotropes is not ruled by nonlocal contributions (π delocalization) but by intra-atomic electron redistribution energies often coined as “hybridization energies”.
Table 9. Relative Phase Stability of Cubic and Hexagonal Phases for C and BN Compoundsa
phaseEx(G)Ecl(G)ΔEint(G)ΔEself(G)bΔE(G)ΔEFHI-aims(G)
C (cubic)–1.2080.06800.1460.0490.009
C (hex.)–1.1820.1390.097000
BN (cubic)–0.604–2.6170B: 00.0120.005
    N: 0.120  
BN (hex.)–0.539–2.6450.037B: 0.07100
    N: 0  
a

Values for a reference unit: G = C1C2 and G = B1N1. ΔE(G) values are referred to the most stable phase. Energies in Ha units.

b

Due to numerical difficulties integrating the total kinetic energy with mHa accuracy, atomic kinetic energies TA were scaled to recover the total kinetic energy from FHI-aims.

The covalent interaction part per diatomic unit G = B1N1 is calculated according to = = −0.528 Ha for ortho (3 B–N(1)), meta (3 B–B(2) + 3 N–N(2)), and para neighbors (3 B–N(3)) respectively. The B–N(4) and N–N(5) interactions in the c direction amount to Ha to yield a total covalent bond energy Ex(B1N1) = −0.539 Ha. Since the Madelung constant for this structure is not fixed by symmetry, it cannot be found tabulated in the literature. Therefore, the ionic electrostatic interaction energy Ha was computed with the Environ program using the QTAIM net charges. This value is slightly larger than Ha obtained for cubic phase. The low symmetry of hexagonal BN allows classic contributions not present in the cubic phase Ecl,L>0(B1N1) = 0.123 Ha that in sum counteract the point-charge Madelung energy. Thus, classic energy components account for Ecl(B1N1) = −2.645 Ha. The portion of covalent interactions of 17% is smaller than 19% obtained for the cubic variant. The total interaction energy obtained amounts to Eint(B1N1) = −3.184 Ha, which is of smaller size than Eint(B1N1) = −3.221 Ha obtained for cubic BN by 37 mHa. In contrast, the total energy difference per BN unit computed with FHI-aims is ΔEFHI-aims(B1N1) = 0.005 Ha energy preference to hexagonal BN. So, the total energies obtained from FHI-aims yield a preference for the hexagonal variants for both cases: carbon, and BN. In contrast, for both cases, the summed interaction energies yield a preference of the cubic structures. This seeming discrepancy is resolved again taking the intra-atomic energies into account as well. According to this finding it is the intra-atomic energies and not the interaction ones that are responsible for the preferred stability of the hexagonal phases (Table 9). For BN, an interesting competition between B and N species is found. Intra-atomic energies indicate a preference of B species for the zincblende environment, while N species prefer the hexagonal environment. In the end the higher stabilization of the N species in the hexagonal environment dominates, and from ΔE(B1N1) = E(B1N1, cubic) – E(B1N1, hex.), the preference for the hexagonal phase is obtained.

Rocksalt-Type Structures

Ionic crystals should present situations with extremely polar bonds. As a case in point we consider LiCl, NaCl, and MgO.
Table 10 shows short-range interactions in the rocksalt-type structure for various compounds. Even nearest-neighbor inter-atomic interactions have a small energy contribution from electron sharing. One can foresee from the delocalization indices that this is the case. Classical interactions take over the exchange component for those systems. As expected from their large QTAIM charges, the classic energy only converges when long-range interactions are also included.
Table 10. IQA Two-Center Terms for Relevant Interactions in LiCl, NaCl, and MgOa
ABmRABδABEnnABEneABEneBAECoulABEclABExAB
Li–Cl(1)62.5650.04410.523–11.076–7.3747.762–0.165–0.007
Li–Li(2)63.6270.0001.313–0.920–0.9200.6450.118–0.000
Cl–Cl(2)63.6270.08442.164–44.393–44.39346.7410.118–0.012
Na–Cl(1)62.8960.06334.173–35.919–31.45233.061–0.138–0.010
Na–Na(2)64.0950.00015.635–14.391–14.39113.2460.099–0.000
Cl–Cl(2)64.0950.04437.344–39.274–39.27241.3010.100–0.005
Mg–O(1)62.1060.12524.116–29.218–20.67425.050–0.725–0.029
Mg–Mg(2)62.9790.00125.579–21.931–21.93118.8030.520–0.000
O–O(2)62.9790.09011.368–13.803–13.80316.7590.522–0.014
a

The net charges of QTAIM atoms are ±0.90e, ± 0.88e, and ±1.71e in LiCl, NaCl, and MgO, respectively. Energies in Ha; distances in Å units.

The additive energy for a reference unit G = LiCl is computed as follows. The covalent interaction energy from interactions with the first and second coordination sphere is Ex(LiCl) = −0.007 × 6 + 0.0 × 6–0.012 × 6 = −0.114 Ha where six Li–Cl(1), six Li–Li(2), and six Cl–Cl(2) interactions are considered. The lattice electrostatic energy assuming point charges is Ha based on a Madelung summation. The effects of nonsphericity of the charge distribution around ions has an effect that is appreciated only for nearest neighbor interactions, as seen from their scaling factors of scl(1)=0.98 and scl(2) = 1.00. To correct the classic energy for first-neighbors, Li–Cl(1), higher order multipolar contributions, taken from the IQA inter-atomic energy as Ha, are included. The total electrostatic energy is therefore Ecl(LiCl) = −0.2848 Ha. Thus, the covalent part contributes to the interaction energy Eint(LiCl) = −0.398 Ha.
Similarly, for a G = NaCl reference unit, covalent interactions contribute with Ex(NaCl) = −0.09 Ha. The classic energy, assuming point charges, is Ha, that is corrected by nonspherical terms Ha. The corrected classic energy Ecl(NaCl) = −0.230 Ha entails a 72% percent of the interaction energy Eint(NaCl) = −0.320 Ha.
MgO exhibits a stronger covalent stabilization than either LiCl or NaCl, Ex(MgO) = −0.258 Ha. However, this increase is paralleled by stronger classical interactions Ha generated by interacting point charges of ±1.71 e. The correction from nonspherical terms is also enhanced, Ha. In overall, the total classic energy is Ecl(MgO) = −1.215 Ha. Thus, the covalent part contributes less, , to the interaction energy Eint(MgO) = −1.473 Ha than the previous systems with rocksalt-type structure.
With respect to interaction energies, the sequence of increasing covalent character is MgO (17%) < NaCl (28%) < LiCl (29%).

Scaled Point-Charge Approximation (sPCA)

The kind of approximation discussed here has been used previously to predict the stability of phases with highly symmetric MgAgAs structure (6) and to discuss Fe–Fe bonding in FeGa3. (4)
The validity of those approximations is tested against formally exact integrals from IQA method. Table 11 compares the zero order (point charge) approximation for interactions where largest deviations are expected, that is, between nearest neighbors.
Table 11. Scaled Point-Charge Approximation (sPCA) (See Eq 35) of Inter-Atomic Bielectronic Integrals for Nearest-Neighbor Interactionsa
SystemRABQAQBδABEclABsclExABsx
β-N21.1080.0030.0032.9990.220b–0.9030.630
CO21.1492.280–1.1401.323–1.4151.182–0.4260.699
Graphite1.4220.000–0.0011.2020.037b–0.3720.832
BN (hex.)1.4462.214–2.2130.452–1.7620.982–0.1350.816
Diamond1.5450.0010.0010.9140.014b–0.2840.907
BN (cubic)1.5612.160–2.1600.357–1.5871.004–0.1060.881
MgB21.782–0.813–0.8130.9830.1660.845–0.2490.852
MgO2.1061.711–1.7120.129–0.7250.984–0.0290.884
LiCl2.5660.897–0.8990.044–0.1650.991–0.0070.819
Al2.8630.0000.0000.2730.0013b–0.0541.071
NaCl2.8960.875–0.8780.064–0.13760.980–0.0100.866
CsCl3.5400.826–0.8270.126–0.10190.996–0.0170.915
Na3.6670.0000.0000.1080.0008b–0.0130.836
a

Energies in Ha, distances in Å, and charges in e units.

b

For interactions between (nearly) noncharged atomic species, the scaling parameters scl become quite large, because the electrostatic interaction is then no longer dominated by a monopolar term. Nevertheless, the absolute error of this assumption is typically small, because the interactions are weak.

Molecular crystals like β-N2 or CO2 present rather short interactions that can not be modeled with the point-charge approximation. Rather like molecules, both solids show very strong covalent bonds as evidenced by their DIs and inter-atomic ExAB. This is in agreement with chemical wisdom. Covalent bonds break the spherical symmetry of the density leading to larger contributions from higher order multipoles to the classic energy. Even more, in N2 the atomic charge is zero but there is still a 0.220 Ha classic destabilization.
Graphite and diamond feature fractional double bonds and covalent single bonds, respectively, which display a corresponding increase of inter-atomic distances that leads to a corresponding decrease of EclCC.
An inflection point is seen for zincblende BN. Despite having also covalent single bonds like diamond, the classic energy closely follows the point-charge approximation scl ≈ 1, indicating that the electron density is close to the spherical one.
Ionic and metallic crystals typically display longer inter-atomic distances and conceptual bond orders <1. Therefore, their covalent and ionic interactions are expected to be accurately approximated with the simple scaled point-charge equations applying an intermediate scaling factor sx ≈ 0.75–0.85.

Conclusions

ARTICLE SECTIONS
Jump To

The new software package ChemInt was developed to decompose the energy of crystalline solids according to the Interacting Quantum Atoms approach. The QTAIM basins are chosen as domains for the energy decomposition. Other possible spatial partitions (i.e., Becke fuzzy atoms, Hirshfeld atoms, ...) deserve future research and will assess the robustness of our conclusions. A number of key issues were addressed which limited the accuracy of the integration. The package ChemInt was used for the examination of some prototype crystalline solids. For the ionic solids studied, covalent interactions contribute 15–30% to the interaction energy. The scale parameters of our scaled point-charge approximation were obtained from the formally exact integration of covalent and classic inter-atomic interaction energies. While certain trends for these parameters are already visible, further investigations are needed to faithfully predict them for all kinds of bonding situations. The phase stability of graphite over diamond and hexagonal BN over zincblende BN was obtained in terms of chemically meaningful energy components. Future explorations along this line are a promising area of research. The obvious limitations of the application to very large systems can be significantly reduced in the future with the exploitation of crystal symmetries and treating only valence electrons. The new program may contribute toward a better understanding of interactions in crystalline materials, in particular in those which do not follow traditional valence rules. Extending the IQA method to planewave bases is perfectly possible, and it can be implemented on any electron structure package because it does not depend on the representation of the wave function in terms of atomic orbitals.

Supporting Information

ARTICLE SECTIONS
Jump To

The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.jpca.1c06574.

  • Computational details of ChemInt calculations and supporting IQA results (PDF)

Terms & Conditions

Most electronic Supporting Information files are available without a subscription to ACS Web Editions. Such files may be downloaded by article for research use (if there is a public use license linked to the relevant article, that license may permit other uses). Permission may be obtained from ACS for other uses through requests via the RightsLink permission system: http://pubs.acs.org/page/copyright/permissions.html.

Author Information

ARTICLE SECTIONS
Jump To

Acknowledgments

ARTICLE SECTIONS
Jump To

D.M.C. thanks J. L. Casals-Sainz for suggesting spherical t-design grids as an alternative to Lebedev–Laikov grids. A.M.P. and E.F. thank the Spanish MICINN, Project PGC2018-095953-B-I00, for funding.

References

ARTICLE SECTIONS
Jump To

This article references 54 other publications.

  1. 1
    Hume-Rothery, W. Materials Science and Engineering; McGraw Hill: 1967; pp 323.
  2. 2
    Fredrickson, D. C.; Lee, S.; Hoffmann, R. The Nowotny Chimney Ladder Phases: Whence the 14 Electron Rule?. Inorg. Chem. 2004, 43, 61596167,  DOI: 10.1021/ic049897h
  3. 3
    Yannello, V. J.; Fredrickson, D. C. Generality of the 18- n Rule: Intermetallic Structural Chemistry Explained through Isolobal Analogies to Transition Metal Complexes. Inorg. Chem. 2015, 54, 1138511398,  DOI: 10.1021/acs.inorgchem.5b02016
  4. 4
    Wagner, F. R.; Cardoso-Gil, R.; Boucher, B.; Wagner-Reetz, M.; Sichelschmidt, J.; Gille, P.; Baenitz, M.; Grin, Y. On Fe–Fe Dumbbells in the Ideal and Real Structures of FeGa3. Inorg. Chem. 2018, 57, 1290812919,  DOI: 10.1021/acs.inorgchem.8b02094
  5. 5
    Bende, D.; Grin, Y.; Wagner, F. R. Covalence and Ionicity in MgAgAs-type Compounds. Chem. - Eur. J. 2014, 20, 97029708,  DOI: 10.1002/chem.201400299
  6. 6
    Bende, D.; Wagner, F. R.; Sichevych, O.; Grin, Y. Chemical Bonding Analysis as a Guide for the Preparation of New Compounds: The Case of VIrGe and HfPtGe. Angew. Chem., Int. Ed. 2017, 56, 13131318,  DOI: 10.1002/anie.201610029
  7. 7
    Martín Pendás, Á.; Blanco, M. A.; Francisco, E. Two-electron Integrations in the Quantum Theory of Atoms in Molecules. J. Chem. Phys. 2004, 120, 45814592,  DOI: 10.1063/1.1645788
  8. 8
    Blanco, M. A.; Martín Pendás, Á.; Francisco, E. Interacting Quantum Atoms: A Correlated Energy Decomposition Scheme Based on the Quantum Theory of Atoms in Molecules. J. Chem. Theory Comput. 2005, 1, 10961109,  DOI: 10.1021/ct0501093
  9. 9
    Martín Pendás, Á.; Francisco, E.; Blanco, M. A. Two-electron Integrations in the Quantum Theory of Atoms in Molecules with Correlated Wave Functions. J. Comput. Chem. 2005, 26, 344351,  DOI: 10.1002/jcc.20173
  10. 10
    Raupach, M.; Tonner, R. A Periodic Energy Decomposition Analysis Method for the Investigation of Chemical Bonding in Extended Systems. J. Chem. Phys. 2015, 142, 194105,  DOI: 10.1063/1.4919943
  11. 11
    Dronskowski, R. Computational Chemistry of Solid State Materials; John Wiley & Sons, Ltd.: 2005.
  12. 12
    Bader, R. F. W. Atoms in Molecules: A Quantum Theory (International Series of Monographs on Chemistry); Clarendon Press: 1990.
  13. 13
    Martín Pendás, Á.; Francisco, E. Promolden: A QTAIM/IQA code. Available from the authors upon request by writing to [email protected]
  14. 14
    Kohout, M. DGrid, ver. 5.2; Dresden: 2021.
  15. 15
    te Velde, G.; Bickelhaupt, F. M.; Baerends, E. J.; Fonseca Guerra, C.; van Gisbergen, S. J. A.; Snijders, J. G.; Ziegler, T. Chemistry with ADF. J. Comput. Chem. 2001, 22, 931967,  DOI: 10.1002/jcc.1056
  16. 16
    Blum, V.; Gehrke, R.; Hanke, F.; Havu, P.; Havu, V.; Ren, X.; Reuter, K.; Scheffler, M. Ab Initio Molecular Simulations with Numeric Atom-Centered Orbitals. Comput. Phys. Commun. 2009, 180, 21752196,  DOI: 10.1016/j.cpc.2009.06.022
  17. 17
    Baranov, A. I.; Kohout, M. Electron Localization and Delocalization Indices for Solids. J. Comput. Chem. 2011, 32, 20642076,  DOI: 10.1002/jcc.21784
  18. 18
    Chen, H.; Friesecke, G. Pair Densities in Density Functional Theory. Multiscale Model. Simul. 2015, 13, 12591289,  DOI: 10.1137/15M1014024
  19. 19
    Rodríguez, J. I.; Köster, A. M.; Ayers, P. W.; Santos-Valle, A.; Vela, A.; Merino, G. An Efficient Grid-Based Scheme to Compute QTAIM Atomic Properties without Explicit Calculation of Zero-Flux Surfaces. J. Comput. Chem. 2009, 30, 10821092,  DOI: 10.1002/jcc.21134
  20. 20
    Womersley, R. S. Contemporary Computational Mathematics - A Celebration of the 80th Birthday of Ian Sloan; Springer International Publishing: Cham, Switzerland, 2018; pp 12431285.
  21. 21
    Blanco, M. A. Métodos Cuánticos Locales para la Simulación de Materiales Iónicos. Fundamentos, algoritmos y aplicaciones. Ph.D. thesis, University of Oviedo: 1997.
  22. 22
    Kay, K. G.; Todd, H. D.; Silverstone, H. J. Bipolar Expansion for r12nYlm1212). J. Chem. Phys. 1969, 51, 23632367,  DOI: 10.1063/1.1672353
  23. 23
    Buehler, R. J.; Hirschfelder, J. O. Bipolar Expansion of Coulombic Potentials. Phys. Rev. 1951, 83, 628633,  DOI: 10.1103/PhysRev.83.628
  24. 24
    Luaña, Víctor. Environ: local environment determination of arbitrary positions in a crystal; Oviedo: 1992.
  25. 25
    Kosov, D. S.; Popelier, P. L. A. Atomic Partitioning of Molecular Electrostatic Potentials. J. Phys. Chem. A 2000, 104, 73397345,  DOI: 10.1021/jp0003407
  26. 26
    Francisco, E.; Menéndez Crespo, D.; Costales, A.; Martín Pendás, Á. A Multipolar Approach to the Interatomic Covalent Interaction Energy. J. Comput. Chem. 2017, 38, 816829,  DOI: 10.1002/jcc.24758
  27. 27
    Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 38653868,  DOI: 10.1103/PhysRevLett.77.3865
  28. 28
    Barca, G. M. J.; Bertoni, C.; Carrington, L.; Datta, D.; De Silva, N.; Deustua, J. E.; Fedorov, D. G.; Gour, J. R.; Gunina, A. O.; Guidez, E. Recent Developments in the General Atomic and Molecular Electronic Structure System. J. Chem. Phys. 2020, 152, 154102,  DOI: 10.1063/5.0005188
  29. 29
    Dirac, P. A. M. Note on Exchange Phenomena in the Thomas Atom. Math. Proc. Cambridge Philos. Soc. 1930, 26, 376385,  DOI: 10.1017/S0305004100016108
  30. 30
    Bloch, F. Bemerkung zur Elektronentheorie des Ferromagnetismus und der elektrischen Leitfähigkeit. Eur. Phys. J. A 1929, 57, 545555,  DOI: 10.1007/BF01340281
  31. 31
    Vosko, S. H.; Wilk, L.; Nusair, M. Accurate Spin-Dependent Electron Liquid Correlation Energies for Local Spin Density Calculations: a Critical Analysis. Can. J. Phys. 1980, 58, 12001211,  DOI: 10.1139/p80-159
  32. 32
    Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785789,  DOI: 10.1103/PhysRevB.37.785
  33. 33
    Becke, A. D. Correlation Energy of an Inhomogeneous Electron Gas: A Coordinate-Space Model. J. Chem. Phys. 1988, 88, 10531062,  DOI: 10.1063/1.454274
  34. 34
    Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab Initio Calculation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Force Fields. J. Phys. Chem. 1994, 98, 1162311627,  DOI: 10.1021/j100096a001
  35. 35
    Downs, R. T.; Somayazulu, M. S. Carbon Dioxide at 1.0 GPa. Acta Crystallogr., Sect. C: Cryst. Struct. Commun. 1998, 54, 897898,  DOI: 10.1107/S0108270198001140
  36. 36
    Bindzus, N.; Straasø, T.; Wahlberg, N.; Becker, J.; Bjerg, L.; Lock, N.; Dippel, A.-C.; Iversen, B. B. Experimental Determination of Core Electron Deformation in Diamond. Acta Crystallogr., Sect. A: Found. Adv. 2014, 70, 3948,  DOI: 10.1107/S2053273313026600
  37. 37
    Dobrzhinetskaya, L. F.; Wirth, R.; Yang, J.; Green, H. W.; Hutcheon, I. D.; Weber, P. K.; Grew, E. S. Qingsongite, Natural Cubic Boron Nitride: The First Boron Mineral from the Earth’s Mantle. Am. Mineral. 2014, 99, 764772,  DOI: 10.2138/am.2014.4714
  38. 38
    Didier, C.; Pang, W. K.; Guo, Z.; Schmid, S.; Peterson, V. K. Phase Evolution and Intermittent Disorder in Electrochemically Lithiated Graphite Determined Using in Operando Neutron Diffraction. Chem. Mater. 2020, 32, 25182531,  DOI: 10.1021/acs.chemmater.9b05145
  39. 39
    Li, M.-R.; Deng, Z.; Lapidus, S. H.; Stephens, P. W.; Segre, C. U.; Croft, M.; Paria Sena, R.; Hadermann, J.; Walker, D.; Greenblatt, M. Ba3(Cr0.97(1)Te0.03(1))2TeO9: in Search of Jahn-Teller Distorted Cr(II) Oxide. Inorg. Chem. 2016, 55, 1013510142,  DOI: 10.1021/acs.inorgchem.6b01047
  40. 40
    Tsirelson, V.; Stash, A.; Kohout, M.; Rosner, H.; Mori, H.; Sato, S.; Lee, S.; Yamamoto, A.; Tajima, S.; Grin, Y. Features of the Electron Density in Magnesium Diboride: Reconstruction from X-ray Diffraction Data and Comparison with TB-LMTO and FPLO Calculations. Acta Crystallogr., Sect. B: Struct. Sci. 2003, 59, 575583,  DOI: 10.1107/S0108768103012072
  41. 41
    Ievinǎ, A.; Straumanis, M.; Karlsons, K. Präzisionsbestimmung von Gitterkonstanten hygroskopischer Verbindungen (LiCl, NaBr). Z. Phys. Chem. 1938, 40B, 146150,  DOI: 10.1515/zpch-1938-4009
  42. 42
    Walker, D.; Verma, P. K.; Cranswick, L. M.; Jones, R. L.; Clark, S. M.; Buhre, S. Halite-sylvite Thermoelasticity. Am. Mineral. 2004, 89, 204210,  DOI: 10.2138/am-2004-0124
  43. 43
    Ewais, E. M.; El-Amir, A. A.; Besisa, D. H.; Esmat, M.; El-Anadouli, B. E. Synthesis of Nanocrystalline MgO/MgAl2O4 Spinel Powders from Industrial Wastes. J. Alloys Compd. 2017, 691, 822833,  DOI: 10.1016/j.jallcom.2016.08.279
  44. 44
    Barrett, C. S. X-ray Study of the Alkali Metals at Low Temperatures. Acta Crystallogr. 1956, 9, 671677,  DOI: 10.1107/S0365110X56001790
  45. 45
    Streib, W. E.; Jordan, T. H.; Lipscomb, W. N. Single-crystal X-Ray Diffraction Study of β Nitrogen. J. Chem. Phys. 1962, 37, 29622965,  DOI: 10.1063/1.1733125
  46. 46
    Schuch, A. F.; Mills, R. L. Crystal Structures of the Three Modifications of Nitrogen 14 and Nitrogen 15 at High Pressure. J. Chem. Phys. 1970, 52, 60006008,  DOI: 10.1063/1.1672899
  47. 47
    Francisco, E.; Casals-Sainz, J. L.; Rocha-Rinza, T.; Martín Pendás, Á. Partitioning the DFT Exchange-Correlation Energy in Line with the Interacting Quantum Atoms Approach. Theor. Chem. Acc. 2016, 135, 170,  DOI: 10.1007/s00214-016-1921-x
  48. 48
    Gelessus, A.; Thiel, W.; Weber, W. Multipoles and Symmetry. J. Chem. Educ. 1995, 72, 505508,  DOI: 10.1021/ed072p505
  49. 49
    Bende, D. Chemical Bonding Models and Their Implications for Bonding-Property Relations in MgAgAs-Type and Related Compounds.; Ph.D. Dissertation; p 97.
  50. 50
    Martín Pendás, Á.; Casals-Sainz, J. L.; Francisco, E. On Electrostatics, Covalency, and Chemical Dashes: Physical Interactions versus Chemical Bonds. Chem. - Eur. J. 2019, 25, 309314,  DOI: 10.1002/chem.201804160
  51. 51
    Wagner, F. R.; Baranov, A. I.; Grin, Y.; Kohout, M. A Position-Space View on Chemical Bonding in Metal Diborides with AlB2 Type of Crystal Structure. Z. Anorg. Allg. Chem. 2013, 639, 20252035,  DOI: 10.1002/zaac.201200523
  52. 52
    Barin, I. Thermochemical Data of Pure Substances; Wiley: 1995; Vol. I.
  53. 53
    Grochala, W. Diamond: Electronic Ground State of Carbon at Temperatures Approaching 0 K. Angew. Chem., Int. Ed. 2014, 53, 36803683,  DOI: 10.1002/anie.201400131
  54. 54
    Popov, I. V.; Görne, A. L.; Tchougréeff, A. L.; Dronskowski, R. Relative Stability of Diamond and Graphite as Seen Through Bonds and Hybridizations. Phys. Chem. Chem. Phys. 2019, 21, 1096110969,  DOI: 10.1039/C8CP07592A

Cited By

This article has not yet been cited by other publications.

  • Abstract

    Figure 1

    Figure 1. Energy difference (per atom) between the reconstructed IQA energy and the SCF energy (energy error). The white band denotes the median. Boxes delimit the lower and upper quartiles Q3 and Q4, and fences delimit all energy errors from a sample that takes wave functions computed with GAMESS, ADF, or FHI-aims, employing HF, LDA, PBE, BLYP, and B3LYP functionals for every system. Ex was rescaled from the integration of the corresponding functional. (47) All errors fall below 1 mHa.

    Figure 2

    Figure 2. Diamond: Convergence of bielectronic integrals inside the basin of the carbon atom, with increasing multipolar order l, and . Only their magnitude is plotted. Dashed lines indicate the trend of convergence for symmetry allowed terms. Disconnected dots are not allowed by symmetry and are nonzero due to numerical errors. Terms below 10–5 Ha are represented as dots at the bottom.

    Figure 3

    Figure 3. BN (zincblende): Convergence of bielectronic intrabasin integrals with increasing multipolar order l. For point and lines explanations, see Figure 2.

    Figure 4

    Figure 4. Diamond: Convergence of bielectronic interbasin integrals with increasing bipolar order L = lA + lB. Labels in the top right indicate an interaction of atom A with an atom B of the ith coordination sphere as AB(i).

    Figure 5

    Figure 5. BN (zincblende): Convergence of bielectronic interbasin integrals with increasing bipolar order L = lA + lB. For point and lines explanations, see Figure 2.

  • References

    ARTICLE SECTIONS
    Jump To

    This article references 54 other publications.

    1. 1
      Hume-Rothery, W. Materials Science and Engineering; McGraw Hill: 1967; pp 323.
    2. 2
      Fredrickson, D. C.; Lee, S.; Hoffmann, R. The Nowotny Chimney Ladder Phases: Whence the 14 Electron Rule?. Inorg. Chem. 2004, 43, 61596167,  DOI: 10.1021/ic049897h
    3. 3
      Yannello, V. J.; Fredrickson, D. C. Generality of the 18- n Rule: Intermetallic Structural Chemistry Explained through Isolobal Analogies to Transition Metal Complexes. Inorg. Chem. 2015, 54, 1138511398,  DOI: 10.1021/acs.inorgchem.5b02016
    4. 4
      Wagner, F. R.; Cardoso-Gil, R.; Boucher, B.; Wagner-Reetz, M.; Sichelschmidt, J.; Gille, P.; Baenitz, M.; Grin, Y. On Fe–Fe Dumbbells in the Ideal and Real Structures of FeGa3. Inorg. Chem. 2018, 57, 1290812919,  DOI: 10.1021/acs.inorgchem.8b02094
    5. 5
      Bende, D.; Grin, Y.; Wagner, F. R. Covalence and Ionicity in MgAgAs-type Compounds. Chem. - Eur. J. 2014, 20, 97029708,  DOI: 10.1002/chem.201400299
    6. 6
      Bende, D.; Wagner, F. R.; Sichevych, O.; Grin, Y. Chemical Bonding Analysis as a Guide for the Preparation of New Compounds: The Case of VIrGe and HfPtGe. Angew. Chem., Int. Ed. 2017, 56, 13131318,  DOI: 10.1002/anie.201610029
    7. 7
      Martín Pendás, Á.; Blanco, M. A.; Francisco, E. Two-electron Integrations in the Quantum Theory of Atoms in Molecules. J. Chem. Phys. 2004, 120, 45814592,  DOI: 10.1063/1.1645788
    8. 8
      Blanco, M. A.; Martín Pendás, Á.; Francisco, E. Interacting Quantum Atoms: A Correlated Energy Decomposition Scheme Based on the Quantum Theory of Atoms in Molecules. J. Chem. Theory Comput. 2005, 1, 10961109,  DOI: 10.1021/ct0501093
    9. 9
      Martín Pendás, Á.; Francisco, E.; Blanco, M. A. Two-electron Integrations in the Quantum Theory of Atoms in Molecules with Correlated Wave Functions. J. Comput. Chem. 2005, 26, 344351,  DOI: 10.1002/jcc.20173
    10. 10
      Raupach, M.; Tonner, R. A Periodic Energy Decomposition Analysis Method for the Investigation of Chemical Bonding in Extended Systems. J. Chem. Phys. 2015, 142, 194105,  DOI: 10.1063/1.4919943
    11. 11
      Dronskowski, R. Computational Chemistry of Solid State Materials; John Wiley & Sons, Ltd.: 2005.
    12. 12
      Bader, R. F. W. Atoms in Molecules: A Quantum Theory (International Series of Monographs on Chemistry); Clarendon Press: 1990.
    13. 13
      Martín Pendás, Á.; Francisco, E. Promolden: A QTAIM/IQA code. Available from the authors upon request by writing to [email protected]
    14. 14
      Kohout, M. DGrid, ver. 5.2; Dresden: 2021.
    15. 15
      te Velde, G.; Bickelhaupt, F. M.; Baerends, E. J.; Fonseca Guerra, C.; van Gisbergen, S. J. A.; Snijders, J. G.; Ziegler, T. Chemistry with ADF. J. Comput. Chem. 2001, 22, 931967,  DOI: 10.1002/jcc.1056
    16. 16
      Blum, V.; Gehrke, R.; Hanke, F.; Havu, P.; Havu, V.; Ren, X.; Reuter, K.; Scheffler, M. Ab Initio Molecular Simulations with Numeric Atom-Centered Orbitals. Comput. Phys. Commun. 2009, 180, 21752196,  DOI: 10.1016/j.cpc.2009.06.022
    17. 17
      Baranov, A. I.; Kohout, M. Electron Localization and Delocalization Indices for Solids. J. Comput. Chem. 2011, 32, 20642076,  DOI: 10.1002/jcc.21784
    18. 18
      Chen, H.; Friesecke, G. Pair Densities in Density Functional Theory. Multiscale Model. Simul. 2015, 13, 12591289,  DOI: 10.1137/15M1014024
    19. 19
      Rodríguez, J. I.; Köster, A. M.; Ayers, P. W.; Santos-Valle, A.; Vela, A.; Merino, G. An Efficient Grid-Based Scheme to Compute QTAIM Atomic Properties without Explicit Calculation of Zero-Flux Surfaces. J. Comput. Chem. 2009, 30, 10821092,  DOI: 10.1002/jcc.21134
    20. 20
      Womersley, R. S. Contemporary Computational Mathematics - A Celebration of the 80th Birthday of Ian Sloan; Springer International Publishing: Cham, Switzerland, 2018; pp 12431285.
    21. 21
      Blanco, M. A. Métodos Cuánticos Locales para la Simulación de Materiales Iónicos. Fundamentos, algoritmos y aplicaciones. Ph.D. thesis, University of Oviedo: 1997.
    22. 22
      Kay, K. G.; Todd, H. D.; Silverstone, H. J. Bipolar Expansion for r12nYlm1212). J. Chem. Phys. 1969, 51, 23632367,  DOI: 10.1063/1.1672353
    23. 23
      Buehler, R. J.; Hirschfelder, J. O. Bipolar Expansion of Coulombic Potentials. Phys. Rev. 1951, 83, 628633,  DOI: 10.1103/PhysRev.83.628
    24. 24
      Luaña, Víctor. Environ: local environment determination of arbitrary positions in a crystal; Oviedo: 1992.
    25. 25
      Kosov, D. S.; Popelier, P. L. A. Atomic Partitioning of Molecular Electrostatic Potentials. J. Phys. Chem. A 2000, 104, 73397345,  DOI: 10.1021/jp0003407
    26. 26
      Francisco, E.; Menéndez Crespo, D.; Costales, A.; Martín Pendás, Á. A Multipolar Approach to the Interatomic Covalent Interaction Energy. J. Comput. Chem. 2017, 38, 816829,  DOI: 10.1002/jcc.24758
    27. 27
      Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 38653868,  DOI: 10.1103/PhysRevLett.77.3865
    28. 28
      Barca, G. M. J.; Bertoni, C.; Carrington, L.; Datta, D.; De Silva, N.; Deustua, J. E.; Fedorov, D. G.; Gour, J. R.; Gunina, A. O.; Guidez, E. Recent Developments in the General Atomic and Molecular Electronic Structure System. J. Chem. Phys. 2020, 152, 154102,  DOI: 10.1063/5.0005188
    29. 29
      Dirac, P. A. M. Note on Exchange Phenomena in the Thomas Atom. Math. Proc. Cambridge Philos. Soc. 1930, 26, 376385,  DOI: 10.1017/S0305004100016108
    30. 30
      Bloch, F. Bemerkung zur Elektronentheorie des Ferromagnetismus und der elektrischen Leitfähigkeit. Eur. Phys. J. A 1929, 57, 545555,  DOI: 10.1007/BF01340281
    31. 31
      Vosko, S. H.; Wilk, L.; Nusair, M. Accurate Spin-Dependent Electron Liquid Correlation Energies for Local Spin Density Calculations: a Critical Analysis. Can. J. Phys. 1980, 58, 12001211,  DOI: 10.1139/p80-159
    32. 32
      Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785789,  DOI: 10.1103/PhysRevB.37.785
    33. 33
      Becke, A. D. Correlation Energy of an Inhomogeneous Electron Gas: A Coordinate-Space Model. J. Chem. Phys. 1988, 88, 10531062,  DOI: 10.1063/1.454274
    34. 34
      Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab Initio Calculation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Force Fields. J. Phys. Chem. 1994, 98, 1162311627,  DOI: 10.1021/j100096a001
    35. 35
      Downs, R. T.; Somayazulu, M. S. Carbon Dioxide at 1.0 GPa. Acta Crystallogr., Sect. C: Cryst. Struct. Commun. 1998, 54, 897898,  DOI: 10.1107/S0108270198001140
    36. 36
      Bindzus, N.; Straasø, T.; Wahlberg, N.; Becker, J.; Bjerg, L.; Lock, N.; Dippel, A.-C.; Iversen, B. B. Experimental Determination of Core Electron Deformation in Diamond. Acta Crystallogr., Sect. A: Found. Adv. 2014, 70, 3948,  DOI: 10.1107/S2053273313026600
    37. 37
      Dobrzhinetskaya, L. F.; Wirth, R.; Yang, J.; Green, H. W.; Hutcheon, I. D.; Weber, P. K.; Grew, E. S. Qingsongite, Natural Cubic Boron Nitride: The First Boron Mineral from the Earth’s Mantle. Am. Mineral. 2014, 99, 764772,  DOI: 10.2138/am.2014.4714
    38. 38
      Didier, C.; Pang, W. K.; Guo, Z.; Schmid, S.; Peterson, V. K. Phase Evolution and Intermittent Disorder in Electrochemically Lithiated Graphite Determined Using in Operando Neutron Diffraction. Chem. Mater. 2020, 32, 25182531,  DOI: 10.1021/acs.chemmater.9b05145
    39. 39
      Li, M.-R.; Deng, Z.; Lapidus, S. H.; Stephens, P. W.; Segre, C. U.; Croft, M.; Paria Sena, R.; Hadermann, J.; Walker, D.; Greenblatt, M. Ba3(Cr0.97(1)Te0.03(1))2TeO9: in Search of Jahn-Teller Distorted Cr(II) Oxide. Inorg. Chem. 2016, 55, 1013510142,  DOI: 10.1021/acs.inorgchem.6b01047
    40. 40
      Tsirelson, V.; Stash, A.; Kohout, M.; Rosner, H.; Mori, H.; Sato, S.; Lee, S.; Yamamoto, A.; Tajima, S.; Grin, Y. Features of the Electron Density in Magnesium Diboride: Reconstruction from X-ray Diffraction Data and Comparison with TB-LMTO and FPLO Calculations. Acta Crystallogr., Sect. B: Struct. Sci. 2003, 59, 575583,  DOI: 10.1107/S0108768103012072
    41. 41
      Ievinǎ, A.; Straumanis, M.; Karlsons, K. Präzisionsbestimmung von Gitterkonstanten hygroskopischer Verbindungen (LiCl, NaBr). Z. Phys. Chem. 1938, 40B, 146150,  DOI: 10.1515/zpch-1938-4009
    42. 42
      Walker, D.; Verma, P. K.; Cranswick, L. M.; Jones, R. L.; Clark, S. M.; Buhre, S. Halite-sylvite Thermoelasticity. Am. Mineral. 2004, 89, 204210,  DOI: 10.2138/am-2004-0124
    43. 43
      Ewais, E. M.; El-Amir, A. A.; Besisa, D. H.; Esmat, M.; El-Anadouli, B. E. Synthesis of Nanocrystalline MgO/MgAl2O4 Spinel Powders from Industrial Wastes. J. Alloys Compd. 2017, 691, 822833,  DOI: 10.1016/j.jallcom.2016.08.279
    44. 44
      Barrett, C. S. X-ray Study of the Alkali Metals at Low Temperatures. Acta Crystallogr. 1956, 9, 671677,  DOI: 10.1107/S0365110X56001790
    45. 45
      Streib, W. E.; Jordan, T. H.; Lipscomb, W. N. Single-crystal X-Ray Diffraction Study of β Nitrogen. J. Chem. Phys. 1962, 37, 29622965,  DOI: 10.1063/1.1733125
    46. 46
      Schuch, A. F.; Mills, R. L. Crystal Structures of the Three Modifications of Nitrogen 14 and Nitrogen 15 at High Pressure. J. Chem. Phys. 1970, 52, 60006008,  DOI: 10.1063/1.1672899
    47. 47
      Francisco, E.; Casals-Sainz, J. L.; Rocha-Rinza, T.; Martín Pendás, Á. Partitioning the DFT Exchange-Correlation Energy in Line with the Interacting Quantum Atoms Approach. Theor. Chem. Acc. 2016, 135, 170,  DOI: 10.1007/s00214-016-1921-x
    48. 48
      Gelessus, A.; Thiel, W.; Weber, W. Multipoles and Symmetry. J. Chem. Educ. 1995, 72, 505508,  DOI: 10.1021/ed072p505
    49. 49
      Bende, D. Chemical Bonding Models and Their Implications for Bonding-Property Relations in MgAgAs-Type and Related Compounds.; Ph.D. Dissertation; p 97.
    50. 50
      Martín Pendás, Á.; Casals-Sainz, J. L.; Francisco, E. On Electrostatics, Covalency, and Chemical Dashes: Physical Interactions versus Chemical Bonds. Chem. - Eur. J. 2019, 25, 309314,  DOI: 10.1002/chem.201804160
    51. 51
      Wagner, F. R.; Baranov, A. I.; Grin, Y.; Kohout, M. A Position-Space View on Chemical Bonding in Metal Diborides with AlB2 Type of Crystal Structure. Z. Anorg. Allg. Chem. 2013, 639, 20252035,  DOI: 10.1002/zaac.201200523
    52. 52
      Barin, I. Thermochemical Data of Pure Substances; Wiley: 1995; Vol. I.
    53. 53
      Grochala, W. Diamond: Electronic Ground State of Carbon at Temperatures Approaching 0 K. Angew. Chem., Int. Ed. 2014, 53, 36803683,  DOI: 10.1002/anie.201400131
    54. 54
      Popov, I. V.; Görne, A. L.; Tchougréeff, A. L.; Dronskowski, R. Relative Stability of Diamond and Graphite as Seen Through Bonds and Hybridizations. Phys. Chem. Chem. Phys. 2019, 21, 1096110969,  DOI: 10.1039/C8CP07592A
  • Supporting Information

    Supporting Information

    ARTICLE SECTIONS
    Jump To

    The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.jpca.1c06574.

    • Computational details of ChemInt calculations and supporting IQA results (PDF)


    Terms & Conditions

    Most electronic Supporting Information files are available without a subscription to ACS Web Editions. Such files may be downloaded by article for research use (if there is a public use license linked to the relevant article, that license may permit other uses). Permission may be obtained from ACS for other uses through requests via the RightsLink permission system: http://pubs.acs.org/page/copyright/permissions.html.

Pair your accounts.

Export articles to Mendeley

Get article recommendations from ACS based on references in your Mendeley library.

Pair your accounts.

Export articles to Mendeley

Get article recommendations from ACS based on references in your Mendeley library.

You’ve supercharged your research process with ACS and Mendeley!

STEP 1:
Click to create an ACS ID

Please note: If you switch to a different device, you may be asked to login again with only your ACS ID.

Please note: If you switch to a different device, you may be asked to login again with only your ACS ID.

Please note: If you switch to a different device, you may be asked to login again with only your ACS ID.

MENDELEY PAIRING EXPIRED
Your Mendeley pairing has expired. Please reconnect

This website uses cookies to improve your user experience. By continuing to use the site, you are accepting our use of cookies. Read the ACS privacy policy.

CONTINUE