# From ELF to Compressibility in Solids

Julia Contreras-García, Miriam Marqués, José Manuel Menéndez, José

Manuel Recio

► To cite this version:

Julia Contreras-García, Miriam Marqués, José Manuel Menéndez, José Manuel Recio. From ELF to Compressibility in Solids. International Journal of Molecular Sciences, 2015, 16 (4), pp.8151-8167. 10.3390/ijms16048151 . hal-01224318

HAL Id: hal-01224318

https://hal.sorbonne-universite.fr/hal-01224318

Submitted on 6 Nov 2015

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# From ELF to Compressibility in Solids

Julia Contreras-García <fnref n="1" />,<fnref n="2" />,<fnref n="3" />,*, Miriam Marqués <fnref n="1" />,<fnref n="4" />,<fnref n="5" />, José Manuel Menéndez <fnref n="1" />,<fnref n="4" />
and José Manuel Recio <fnref n="1" />,<fnref n="4" />

Academic Editor: Abhijit Chatterjee

Received: 15 January 2015 / Accepted: 18 March 2015 / Published: 13 April 2015

**Abstract:** Understanding the electronic nature of materials' compressibility has always been a major issue behind tabulation and rationalization of bulk moduli. This is especially because this understanding is one of the main approaches to the design and proposal of new materials with a desired (e.g., ultralow) compressibility. It is well recognized that the softest part of the solid will be the one responsible for its compression at the first place. In chemical terms, this means that the valence will suffer the main consequences of pressurization. It is desirable to understand this response to pressure in terms of the valence properties (charge, volume, etc.). One of the possible approaches is to consider models of electronic separability, such as the bond charge model (BCM), which provides insight into the cohesion of covalent crystals in analogy with the classical ionic model. However, this model relies on empirical parametrization of bond and lone pair properties. In this contribution, we have coupled electron localization function (ELF) *ab initio* data with the bond charge model developed by Parr in order to analyze solid state compressibility from first principles and moreover, to derive general trends and shed light upon superhard behavior.

<footnote n="1">¹ MALTA-Consolider Team, Universidad de Oviedo, E-33006 Oviedo, Spain;
E-Mails: miriam.marques@uva.es (M.M.); jmenendezmontes@gmail.com (J.M.M.);
jmrecio@uniovi.es (J.M.R.)</footnote>

<footnote n="2">² Université Pierre et Marie Curie, Univ Paris 06, UMR 7616, Laboratoire de Chimie Théorique,
F-75005 Paris, France</footnote>

<footnote n="3">³ Centre National de la Recherche Scientifique, Unité Mixte de Recherche 7616,
Laboratoire de Chimie Théorique, F-75005 Paris, France</footnote>

<footnote n="4">⁴ Departamento de Química Física y Analítica, Universidad de Oviedo, E-33006 Oviedo, Spain</footnote>

<footnote n="5">⁵ Departamento de Física Teórica, Universidad de Valladolid, E-47011 Valladolid, Spain</footnote>

<footnote>* Author to whom correspondence should be addressed; E-Mail: contrera@lct.jussieu.fr;
Tel.: +33-1-4427-3879.</footnote>
**Keywords:** ELF; compressibility; energy models

# 1. Introduction

A thorough understanding of many important electronic and mechanical properties of a crystalline material does not require the knowledge of the total energy of the solid, but just the energetics of its bonding network. Many phenomenological approaches have been developed to illustrate the advantage of this view [1]. These so-called bond energy models are useful for properties prediction and have become rather popular in the last years as trusty guides for material design. Specifically to the field of materials strength, Gilman's book [1] collects a number of valuable correlations between a variety of observable elastic moduli of different type of solids and electronic densities coming only from the active electrons of the valence space. Based on the ideas developed by Gilman [2], Gao *et al.* "initialized a link" between macroscopic observables and first principles calculations with the aim to evaluate the indentation hardness of multicomponent crystals [3]. Li *et al.* resorted to bond electronegativity in order to evaluate the bond energy and thus figure out the most stable bond combination for hard materials [4]. In conjunction with evolutionary algorithms [5], Li's model is theoretically generalized by Lyakhov and Oganov making use of the bond-valence model and graph theory [6].

Inherent to all these semi-empirical models, there are still elusive concepts to quantum-mechanical formalisms, as ionicity/covalent scales, electronegativities, covalent radii, *etc*. Efforts from the theoretical chemistry community have been devoted to reformulate these concepts from first principles. As an example, we can cite recent reformulations of electronegativity [7–9]. These chemical concepts do not derive directly from wavefunction analysis, but a posteriori models are needed to extract them. The bondon [10] and the reformulation of chemical bonding as a quantum condensate are some of the formulations that have been recently proposed [11]. Another approach that has gained many followers is the topological formalism. This approach divides three-dimensional space into chemically meaningful disjoint regions following the gradient of a scalar function, such as the electron density [12,13] or the electron localization function (ELF) [14,15]. Since we are interested in valence properties, our function of choice is the latter one. ELF is able to recover Lewis entities, such as bonds and lone pairs and quantify their properties (volume, charge, *etc*.). Due to its direct link with Lewis theory, we have chosen the model by Borkman and Parr [16] as our starting point to merge topology and bond energy models.

Briefly, our study is based on the success of the bond charge model [16–18] in predicting molecular properties, and the ability of the electron localization function [14,15] to provide the required parameters for its free from empirical data application. The bulk modulus (inverse of the compressibility) will be the focus of this first principles-like extension of BCM to the solid state.

The remainder of the paper is organized as follows. Firstly, the original definition of the bond charge model for homo and heteronuclear binary molecules and covalent solids will be presented. Secondly, the computational details for the calculations performed over a broad group of diamond-type and zinc-blende-type solids from the IV, III–V and II–VI groups are given. In Section 3, the basic parameters of this *ab-initio* BCM (hereafter, NEWBCM) are derived from first principles. The ability of
this model to predict solid state properties (compressibility) will be tested against experimental data in Section 4. At the end of the section, this new insight into macroscopic properties will be used to understand macro-micro requirements for superhard materials. The article finishes with the main conclusions and guidelines for prospective work.

## 2. Theoretical Background

In this section, the two approaches to be merged are introduced: the bond charge model and the ELF topological formalism.

### 2.1. The Bond Charge Model: Original Definition

#### 2.1.1. Homo and Heteronuclear Diatomic Molecules

Within the BCM formalism [16], a homonuclear diatomic molecule ($A_2$) can be approximated by two cores holding a $q/2$ positive charge each, and a bond charge in between them, $-q$, that moves freely along a bond length $R_B = \nu R$ which is a fraction ($\nu \le 1$) of the interatomic distance $R$ (see Figure 1 (left)). The energy associated to this model at distances close to the equilibrium ($R \approx R_e$) is given in atomic units by:

$$E(q, R) = E_0 + E_1 + E_2 = 2E_A - \frac{Cq^2}{R} + \frac{D'q}{R_B^2} \quad (1)$$

where $E_0$ represents the core energy (which equals to $2E_A$ for an $A_2$ molecule), $E_1$ refers to the coulombic attractive interactions and $E_2$ is related to the kinetic energy of the bond electrons approximated as particles in a box of length $R_B$. The constants, which are system independent, are then given by (see [16] for details):

$$C = \frac{7}{4}, \quad D' = \frac{\pi^2}{2m} \quad (2)$$

and the parameters $q$ and $\nu$ are inferred from empirical measures of the harmonic force constant, $k_e$, and the equilibrium distance, $R_e$.

In order to adapt this model to heteronuclear diatomic molecules, $AX$, or even to homonuclear molecules in an asymmetric surrounding, two considerations have to be taken into account. On the one hand, the polarity of the bond clearly influences the situation of the bond charge (Figure 1 (right)), being $r_1$ and $r_2$ the distances to its position from atoms $A$ and $X$, respectively. On the other hand, the positive charge lying on each of the two nuclei will be different from $\frac{1}{2}q$ by $\frac{1}{2}\delta$ and $-\frac{1}{2}\delta$ [18]. The changes induced by the electronegativity difference, $\Delta\chi = \chi_A - \chi_X$, in $E_2$ are negligible, but the bond charge shift gives rise to noticeable changes in the coulombic interactions:

$$E_1 = \frac{q^2}{2} \left[ \frac{(1 - \delta^2)}{2R} - \frac{(1 + \delta)}{r_1} - \frac{(1 - \delta)}{r_2} \right] \quad (3)$$

so that the total energy reads as follows:

$$E = E_A(q, \delta) + E_X(q, \delta) + \frac{q^2}{2} \left[ \frac{(1 - \delta^2)}{2R} - \frac{(1 + \delta)}{r_1} - \frac{(1 - \delta)}{r_2} \right] + \frac{D'q}{R_B^2} \quad (4)$$
The determination of the model parameters, $q$, $\delta$, $r_1/r_2$ and $R_B$ has become a matter of choice. Probably, the most extended method has been to take into account the electronegativity of the atoms as shown by Parr *et al.* [19,20].

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Bond charge model for homonuclear (**left**) and heteronuclear (**right**) molecules, $A_2$ and $AX$, respectively. The meaning of the main parameters that are needed in each case is shown in the figures: $q$ and $R_B$ stand for neutral homonuclear molecules as $q$, $\delta$, $r_1/r_2$ and $R_B$ do for heteronuclear ones; B label stands for the bond (see text for more details). $r'_1$ and $r'_2$ are also shown for clarity.</figcaption></figure>

### 2.1.2. Binary Covalent Solids

The model by Parr and Borkman was extended by Martin [21] to describe the energy of covalent crystals. For simplicity, we consider a zinc-blende-type solid, $AX$, with a formula unit $AXB_4$ accounting for the four-fold coordination of the two elements in the unit cell, $B$ represents the bond entity. According to BCM, the energy of the covalent solid per couple of atoms $AX$ is given by

$$E_{AX} = E_A(q_B, \delta) + E_X(q_B, \delta) + E_B(q_B, \delta, R_B) \quad (5)$$

$$E_B = -\frac{Mq_B^2}{R_B} + \frac{4D'q_B}{R_B^2} \quad (6)$$

where $M$ is the Madelung constant for these lattices with charges $2q_B(1 + \delta)$ and $2q_B(1 - \delta)$ at $A$ and $X$ atomic sites, and $-q_B$ at each bond position, as dictated by the multiplicity and the charge neutrality condition. We should recall that $E_A$ and $E_X$ account for the energy of the corresponding core electrons not involved in the bonding, and, in practice, are independent of geometrical parameters. The greatest modification introduced in the application of BCM to solids is the presence of an infinite sum, accounting for all $A$, $X$, and $B$ charged species, which is subsumed into the Madelung constant. It also important to highlight that due to the inclusion of the bond charge in the summation, this Madelung constant is system dependant since the position of the bond charge will not correspond to a special Wyckoff position (invariable) except for homonuclear solids. The Madelung constants have been calculated thanks to the environ code [22] and are collected in Table 1.
<table><caption>**Table 1.** Parameters of the non-empirical new bond charge model (NEWBCM), $[q_B, R_B, r_1/r_2]$, for a broad group of diamond-type and zinc-blende-type solids from the IV, III–V and II–VI groups. $r_1(r_2)$ stands for the distance between the charge, $q$, and the nucleus $A(X)$. $r'_1(r'_2)$ stands for the part of $R_B$ belonging to $r_1(r_2)$, i.e., $r'_1 = r_1 - r_A$ and $r'_2 = r_2 - r_X$, where $r_A$ is the core size of atom A and $r_X$ is the core size of atom X. $M$ and $E_B$ stand for the Madelung constant and the bonding energy of the crystal, respectively. The cell parameter, $a$, and the electronegativity difference [23], $\Delta\chi = \chi_X - \chi_A$, are also shown for reference. Lengths are in Å, charges in electrons, energy in hartrees.</caption><thead><tr><th>AX</th><th>a</th><th>Δχ</th><th>r<sub>1</sub>/r<sub>2</sub></th><th>r'<sub>1</sub>/r'<sub>2</sub></th><th>R<sub>B</sub></th><th>q<sub>B</sub></th><th>M</th><th>E<sub>B</sub></th><th>δ</th></tr></thead><tbody><tr><td>C</td><td>3.557</td><td>0.0</td><td>1.00</td><td>1.00</td><td>0.938</td><td>1.950</td><td>10.856</td><td>-11.644</td><td>0.000</td></tr><tr><td>Si</td><td>5.431</td><td>0.0</td><td>1.00</td><td>1.00</td><td>1.132</td><td>1.950</td><td>8.577</td><td>-7.623</td><td>0.000</td></tr><tr><td>Ge</td><td>5.658</td><td>0.0</td><td>1.00</td><td>1.00</td><td>1.003</td><td>2.100</td><td>7.414</td><td>-8.625</td><td>0.000</td></tr><tr><td>SiC</td><td>4.360</td><td>0.7</td><td>1.17</td><td>0.75</td><td>0.985</td><td>1.950</td><td>9.356</td><td>-9.556</td><td>0.000</td></tr><tr><td>BN</td><td>3.616</td><td>1.0</td><td>1.07</td><td>0.86</td><td>0.936</td><td>1.975</td><td>11.006</td><td>-12.136</td><td>0.241</td></tr><tr><td>BP</td><td>4.538</td><td>0.1</td><td>0.96</td><td>1.21</td><td>1.050</td><td>2.000</td><td>9.662</td><td>-9.739</td><td>0.250</td></tr><tr><td>BAs</td><td>4.777</td><td>0.0</td><td>0.70</td><td>0.82</td><td>1.023</td><td>2.050</td><td>8.711</td><td>-9.468</td><td>0.268</td></tr><tr><td>BSb</td><td>5.156</td><td>0.1</td><td>0.62</td><td>0.88</td><td>0.999</td><td>1.988</td><td>8.554</td><td>-8.954</td><td>0.245</td></tr><tr><td>AlP</td><td>5.464</td><td>0.6</td><td>1.14</td><td>0.97</td><td>1.122</td><td>1.975</td><td>8.925</td><td>-8.206</td><td>0.266</td></tr><tr><td>AlAs</td><td>5.661</td><td>0.5</td><td>0.96</td><td>0.83</td><td>1.080</td><td>2.025</td><td>8.028</td><td>-8.068</td><td>0.284</td></tr><tr><td>AlSb</td><td>6.136</td><td>0.4</td><td>0.81</td><td>0.77</td><td>1.093</td><td>1.963</td><td>7.304</td><td>-6.810</td><td>0.261</td></tr><tr><td>GaAs</td><td>5.653</td><td>0.4</td><td>1.19</td><td>1.18</td><td>0.992</td><td>2.100</td><td>7.662</td><td>-9.008</td><td>0.238</td></tr><tr><td>GaSb</td><td>6.096</td><td>0.3</td><td>1.00</td><td>1.12</td><td>0.986</td><td>2.038</td><td>6.764</td><td>-7.538</td><td>0.215</td></tr><tr><td>InP</td><td>5.869</td><td>0.4</td><td>1.64</td><td>1.38</td><td>1.012</td><td>1.970</td><td>8.494</td><td>-8.615</td><td>0.269</td></tr><tr><td>InAs</td><td>6.058</td><td>0.3</td><td>1.35</td><td>1.11</td><td>0.958</td><td>2.020</td><td>7.268</td><td>-8.187</td><td>0.287</td></tr><tr><td>CdS</td><td>5.818</td><td>0.8</td><td>1.89</td><td>1.48</td><td>0.976</td><td>1.975</td><td>10.016</td><td>-10.512</td><td>0.494</td></tr><tr><td>CdSe</td><td>6.077</td><td>0.7</td><td>1.56</td><td>1.27</td><td>0.938</td><td>2.025</td><td>8.591</td><td>-9.697</td><td>0.506</td></tr><tr><td>ZnS</td><td>5.410</td><td>0.9</td><td>1.67</td><td>1.50</td><td>1.001</td><td>2.025</td><td>10.130</td><td>-10.976</td><td>0.457</td></tr><tr><td>ZnSe</td><td>5.668</td><td>0.8</td><td>1.37</td><td>1.30</td><td>0.970</td><td>2.075</td><td>8.631</td><td>-10.134</td><td>0.470</td></tr><tr><td>ZnTe</td><td>6.104</td><td>0.5</td><td>1.15</td><td>1.25</td><td>0.971</td><td>2.038</td><td>7.415</td><td>-8.389</td><td>0.460</td></tr></tbody></table>

From Equation (6), the bond energy of a zinc-blende-type lattice can be finally expressed within the framework of BCM as:

$$E_B = -\frac{Mq_B^2}{R_B} + \frac{4D'q_B}{R_B^2} \qquad (7)$$

It is interesting to highlight here that this expression is similar to the classical Madelung energy expression for ionic crystals. However, in our case the equation is done in terms of Lewis entities (cores, bonds) instead of ions. The expression will then hold as long as the orthogonality of these entities holds too. In solids such as diamond, where for example, good localized orbitals can be found, Equation (7)
will give an accurate description of the bonding energy in the crystal. However, as the localized picture fails (e.g., down the periodic table), we should expect the fits to worsen. This equation will progressively fail as we go down the periodic table or move towards metallic solids, due to a less obvious separability.

If we reasonably assume that $(\partial E_B / \partial R) \approx 0$ at equilibrium then

$$ \left( \frac{\partial E_B}{\partial R_B} \right)_{R=R_e} = \frac{M q_B^2}{R_B^2} - \frac{8 D' q_B}{R_B^3} = 0 \qquad (8) $$

where we have assumed that the ratio $r_1/r_2$ remains constant upon compression, not introducing changes in $M$ (see discussion in last section and Figure 5). It is very important to notice that all the parameters have to be evaluated at equilibrium. We then have for the constant $M$:

$$ M = \frac{8D'}{q_B R_B} \qquad (9) $$

Hence, the ground state bond energy at equilibrium is given by:

$$ E_{B,eq} = -\frac{4D'q_B}{R_B^2} = -\frac{Mq_B^2}{2R_B} \qquad (10) $$

where the negative value reflects the stabilization introduced by the chemical bond (see Table 1). This equation, just like in the original bond charge model [16,18], is a manifestation of the virial theorem. Upon substitution of Equation (9) in Equation (7), it can be seen that at equilibrium, the electrostatic interaction is twice as big as the kinetic energy term.

## 2.2. ELF Topology

The Electron Localization Function, ELF, was originally designed by Becke and Edgecombe to identify “localized electronic groups in atomic and molecular systems” [14]. Several interpretations have been given to this function. For example, ELF has been generalized and reinterpreted in the light of Markovian processes [24,25]. A very intuitive interpretation was given by Savin [15], according to which the ELF core, $\chi$, can be understood as a local measure of the excess of local kinetic energy of electrons due to the Pauli principle, $t_p$:

$$ \chi_{\sigma}(r) = \frac{t_p(\vec{r})}{t_{HEG}(\vec{r})} \qquad (11) $$

relative to the homogeneous electron gas kinetic energy density, $t_{HEG}(\vec{r}) = \frac{3}{5}(6\pi^2)^{2/3} \rho_{\sigma}^{5/3}(\vec{r})$. $t_p$ is calculated as the difference between the definite positive kinetic energy, $t(\vec{r})$, and the von Weizsäcker kinetic energy functional:

$$ t_p(\vec{r}) = t(\vec{r}) - \frac{1}{4} \frac{|\rho_{\sigma}(\vec{r})|^2}{\rho_{\sigma}(\vec{r})} \qquad (12) $$

ELF values are confined in the [0,1] range by a Lorentzian transformation which facilitates its interpretation:

$$ ELF = \frac{1}{1 + \chi^2} \qquad (13) $$

ELF maxima appear associated with localized electron pairs, such as bonds, lone pairs and atomic shells, thus matching the Valence Shell Electron Pair Repulsion [26] model. These maxima can be used
to partition the system into regions following the gradient flux (this would be like partitioning a set of mountains by the valleys) [12]. This partitioning gives rise to a set of non overlapping regions (or basins), labelled as $\Omega_i$. Their properties can be determined by integration over the region. Most commonly, the volume integration of the density, $\rho$, assigns a basin population, $\bar{N}(\Omega_i)$:

$$ \bar{N}(\Omega_i) = \int_{\Omega_i} \rho(\vec{r}) d\vec{r} \qquad (14) $$

The properties of these basins have a chemical meaning inherited from the partition, so that the integration of the density over a basin associated with a bond maximum, can be understood as the bond charge.

### 3. Computational Details

Electronic structure calculations for a number of IV, III–V and II–VI solids with the diamond and zinc-blende-type structures were performed within the density-functional theory formalism. These compounds were selected to ensure the existence of a basin associated with the bond.

We have used the all-electron full potential linearised augmented-plane wave (FP-LAPW) approach as implemented in the ELK (for electrons in k-space) code [27]. The FP-LAPW method is among the most precise DFT-based methods for crystal structures [28]. It divides space into an interstitial region and non-overlapping muffin-tin (MT) spheres centered at the atomic sites. In the interstitial region, the basis set is described by plane-waves, whereas in the MT spheres, the basis set is described by radial solutions of the one particle Schrödinger equation and their energy derivatives multiplied by spherical harmonics. A high plane-wave cut-off $K_{max}R_{MT} = 9$, where $R_{MT}$ is the smallest muffin-tin radius of the corresponding atomic species and $K_{max}$ is the maximum size of the reciprocal lattice vectors, guaranteed convergence in all the studied systems. Moreover, spherical harmonics within each atomic sphere were expanded up to the maximum angular quantum number $l_{max} = 14$, to avoid discontinuities of the electron density and its derivatives at the muffin boundaries. Lattice parameters were set to the room temperature experimental values of the corresponding structures [29]. We have used both, the Perdew-Wang [30] local density approximation (LDA) and the Perdew-Burke-Ernzerhof [31] generalized gradient approximation (GGA) for the exchange-correlation functional. Brillouin-zone integrations were approximated using 8 × 8 × 8 Monkhorst-Pack grids, and the self-consistent iteration process was repeated until the absolute change in total energy was lower than $1 \times 10^{-5}$ hartrees.

In a second step, the topology of the ELF functions coming from the previous calculations is investigated using CRITIC, a code developed by some of the authors [32,33]. CRITIC is able not only to find ELF maxima and saddle points but it can also construct basins and integrate properties.

The training set has been chosen to lay in the shared-electron range. This characteristic is necessary so that there is a bond basin linkable to each pair of nearest neighbors in the crystal. In other words, it is possible to assign bonding parameters within the ELF framework.

#### 3.1. Bond Length: $R_B, r_1/r_2$

The most common approach in the analysis of a bonding structure by means of the ELF approach is to use its 3D topology. However, the BCM uses a 1D definition of the chemical bond. Thus, we have
devised an alternative route to obtain the BCM parameters from ELF. We have resorted to the 1D-ELF
topology: ELF is analyzed along the straight line between pairs of atoms. In some cases, specially as
going down the periodic table, this line might differ both from the AIM and ELF bonding lines, but it
enables us to stick to the BCM formalism. On the one hand, the position of the bond charge within
this length, r₁/r₂, (shifted from the center by the relative electronegativity of atoms) is determined by
the local ELF maximum along the same direction. On the other hand, the distance between the two
closest minima to the bond maximum along the internuclear direction provide the bond length, R<sub>B</sub> (See
Figure 2).

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Definition of length-related bond charge model parameters, [$R_B$, $r_1/r_2$], for a homonuclear molecule, $N_2$ (**left**), and a heteronuclear molecule, BN (**right**).</figcaption></figure>

**4. Defining the Parameters from First Principles**

ELF profiles along the internuclear line were constructed for each solid [14]. These profiles were
then analyzed in order to find the position of the local maxima and minima, which enable the quantitative
determination of *ab initio* R<sub>B</sub>, r<sub>1</sub> and r<sub>2</sub> values according to the previous definition. Results are collected
in Table 1.

In homonuclear cases the local maximum is located in the center of R<sub>B</sub> and r<sub>1</sub>/r<sub>2</sub> = 1.0 by definition
(see N<sub>2</sub> in Figure 2 (left)). Instead, in heteronuclear cases, the ELF profile is deformed and the maximum
of the bond (or attractor) is no longer placed in the middle of R<sub>B</sub> (see BN in Figure 2 (right)).

A similar value obtained for R<sub>B</sub> is obtained in all solids, with numbers clustering around 1 Å. This
means that the size of the bond (at least along the interbonding direction) depends neither on the lattice
parameter of the crystal (which spreads a range of more than 2.5 Å), nor on the two atoms involved in
the bond.

Secondly, we observe that r₁/r₂ ≈ 1 ratios are obtained for the purely covalent families and those
whose Δχ ≈ 0. There is one remarkable outlier: BN shows r₁/r₂ ≈ 1 with Δχ=1. In order to
understand this fact, an extra column has been added to Table 1 with the displacement of the attractor
within the bond length R<sub>B</sub>, that is, taking away core contributions: r'₁ = r₁ − r<sub>A</sub>, r'₂ = r₂ − r<sub>X</sub>.

This means that we now decompose the bond length into R<sub>B</sub> = r'<sub>1</sub> + r'<sub>2</sub>, so that R = r<sub>A</sub> + r'<sub>1</sub> + r'<sub>2</sub> + r<sub>X</sub>.
The physical meaning of these quantities is depicted in Figure 1 (right). This quantity enables us to
understand that although the charge in BN is nearly at the same distance from B and from N, it is
displaced within $R_B$, as expected from their electronegativity difference (see Figure 2 (right)). This decomposition also helps to get rid of the core size in the comprehension of the system. Let's take the case of BP: since P belongs to the second row and B to the first one, we have $r_X \gg r_A$, so that even though the center of the bond is displaced within the bond towards P ($r'_1 > r'_2$), once the core radii are taken into account, the distance relationship is inversed ($r_2 > r_1$).

### 4.1. Bond Charge: $q_B$

Obtaining the charge associated with ELF entities through the integration over the corresponding basins might be considered in principle a straightforward and accurate method. However, for heavy elements, some ill-defined cases appear for integration, where the low symmetry of the basins and the appearance of spurious critical points at the muffin surfaces may lead to serious difficulties if we seek a quantitative analysis of its topology.

Based on pseudopotential principles, we assume that the core electronic distributions are practically unaltered by the environment, and thus $q_B$ values can be directly obtained by mere difference between the total and core number of electrons. Core values for the elements up to the 4th row involved in the solids under this study have been extracted from the ELF integrations gathered in [34]. Beyond the 4th row, calculations were perfomed for the atoms with the PBE functional [31] and the 6-31G* basis set, using the Gaussian 03 program [35]. Their core/valence charge, which is very stable to the level of calculation, was calculated with the TopMod program [36]. This assumption has the added value of providing extremely fast calculations. It has also been validated by comparison with the charges obtained by the integration over the corresponding 3D-basins with the CRITIC code [37]. For instance, integrated values for $q_B$ of 1.952 and 1.948 for C and BN, respectively, are in good agreement with the corresponding data in Table 1. They are also independent of the exchange-correlation functional used.

It has been shown that ELF charges lie very close to chemical intuition values [38]. This is further seen in Table 1. Nearly 2 electrons are recovered in all cases for the bonds, in agreement with their bond order of 1. This is different from the typical BCM charges obtained from fitting, which in many cases deviate from the expected values [17,18]. Furthermore, this image of $q_B$ correlating to the bond order does not agree with the definition of the $q_B$ parameter in the traditional model. The step forward that ELF provides concerns also the evaluation of the charge transfer ($\delta$). In all cases (see Table 1), $\delta$ also shows reasonable values following the expected trends for the IV, III–V, and II–VI crystal families dictated by their electronegativity [23]. Keep in mind that $\delta$ is calculated with respect to valence configurations, so that for example carbon and silicon, having the same configuration, [34] give rise to $\delta = 0$ in SiC.

This is an important characteristic of this Non Empirical New Bond Charge Model (NEWBCM), parameters are not only obtained *ab initio*, but they are closer to chemical expectations and recover/predict values in good agreement with experience (at difference with the initial BCM parameters).

## 5. Testing the Model: Response to Pressure

Since pressure can introduce significant changes in the stability of phases, it is most useful in the synthesis of novel phases and metastable materials. Pressure allows precise tuning of a fundamental
parameter, the interatomic distance, which controls the electronic structure and virtually all the interatomic interactions that determine materials properties [39–41]. With pressure tuning, properties can often be more rapidly and cleanly optimized than with chemical tuning, which necessitates the synthesis of a large number of different materials and can induce disorder, phase separation, and other undesirable effects. Pressure tuning is therefore a useful tool in the search for new solid state materials with tailored properties.

## 5.1. Bulk Modulus

Among the description of the responses to pressure, compressibility ($\kappa$) or its inverse, the bulk modulus ($B_0$), probably constitute the most widely used parameters. The inverse relationship between bulk modulus and volume has been widely explored [42–45], enabling to split the response of solids under pressure into local contributions (see for example [43,46–48]).

Following our BCM model, the valence charge is distributed along the chemical bonds, so that the compressibility of solids can be approximated to a good extent by the compressibility of the bonds, $\kappa_B$, whereas cores are assumed to remain untouched: $\kappa \simeq \kappa_B$ [43,49]. The bulk modulus of the solid is then given by (see S.I. for the development):

$$B_{0,B} = \left[ V_B \left( \frac{\partial^2 E_B}{\partial V_B^2} \right) \left( \frac{\partial V_B}{\partial V} \right) \right]_0 = \frac{M q_B^2}{32 \sqrt{3} R_B R^3} \quad (15)$$

where $V_B = \frac{1}{4}[V_{AXB_4} - (V_A + V_X)]$ is the bond’s volume within the AXB$_4$ unit.

Equation (15) for our IV, III–V and II–VI compounds is analyzed in Figure 3. It is worth highlighting that a very good agreement is obtained for all data. The slope for the bulk modulus is around 1 (1.14), and the offset is near 9 GPa (8.99 GPa) which departures from zero but it is not too large if we compare it with the range (up to 500 GPa) displayed by $B_0$ in these compounds. The Pearson regression coefficient $R^2$ is 0.9774. A similar agreement was obtained by Gilman et al. [1] in an enlightening analysis from empirical measures where the relationship was analyzed in terms of valence density (i.e., $q_B/V$). The coherency between these two results can be perfectly understood from the analysis of Equation (15) if we take into account that $R = \frac{\sqrt{3}}{4}a$ and $V = a^3$ for these crystals:

$$B_{0,B} = \frac{2M q_B^2}{9R_B V} \quad (16)$$

In his fit, Gilman made use of the cell volume and an average valence occupation $q_B/V$ as would be expected from Lewis theory (i.e., a default bond occupation $q_B = 2$). Whereas the bond occupation has been found in our calculations to be very close to two, the dependence on the cell size can be improved by taking calculated bond lengths into account. Indeed, Equation (16) shows that the dependence on the lenght parameters should be $B_0 \propto (R_B V)^{-1}$ (i.e., $B_0 \propto (R_B R^3)^{-1}$, exponent nearly –4), whereas Gilman took $B_0 \propto V^{-1}$ alone (exponent on $R = -3$). Indeed, a logarithmic fit of the $B_0$ experimental data vs. $R$ gives a slope of 3.8 (±0.1) (regression coefficient $R^2 = 0.9863$) which explains the need to increase the length dependence from $R^{-3}$ to $R^{-3}R_B^{-1}$. As far as the bond charge is concerned, default occupations are very close to the values obtained with ELF, however, the quantitative values provided by electronic structure calculations enable us to obtain a better linear behavior. This is shown in Figure 4, where
experimental and NEWBCM values have been included (in order to represent experimental values, the bond length dependence has been set to $R^4$). It is interesting to note that, although calculations provide a better linear behavior ($R^2$(NEWBCM) = 0.9927, $R^2$(Exptal.) = 0.9830), the use of default or nominal ($q_B$ = 2) values already provides all the physics.

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Experimental vs. non empirical new bond charge model (NEWBCM) bulk modulus (GPa).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Experimental (empty symbols) and NEWBCM (full symbols) bulk moduli (in GPa) for IV, III–V and II–VI compounds. Bulk moduli in GPa, charge in electrons, *R* in atomic units.</figcaption></figure>
## 5.2. From Micro- to Macroscopic: Understanding Compressibility

Understanding the inherent relationship between (quantifiable) microscopic properties and macroscopic properties is the corner stone of material design.

We have shown that the concepts derived from ELF can be used to understand the microscopic origin of solids compressibility. If we analyze the dependencies of Equation (15), we see that the least compressible materials are found within the first period, where smaller interatomic distances, $R \propto V^{1/3}$, and bond lengths, $R_B$, are found. While these results were already known from experience, the microscopic approach can also reveal new information. With this aim in mind, it is interesting to study the evolution of NEWBCM parameters ($r_1, r_2, R_B$ and $q_B$) upon compression.

Let's start with $R_B$. Since $R_B = \nu R$ can be expressed as $R_B = R - r_A - r_X$, and the core radii remain practicably untouched upon compression, a linear relationship necessarily holds for $R_B$ (which virtually subsumes all the changes in the cell parameter) with $R$ for a given material.

Furthermore (and less intuitively), the linear relationship also holds for $r_1$ and $r_2$ with $R$ at constant $\delta$. According to Parr's electronegativity approach [20], the electronegativity in the molecule AX, $\chi_{AX}$, can be expressed in terms of the radii of the isolated atoms, $r_A$ and $r_X$, as: [19]

$$ \chi_{AX} = \frac{\chi_A r_A + \chi_X r_X}{R} \qquad (17) $$

but also in terms of their properties in the molecule AX:

$$ \chi_{AX} = K \frac{q + \delta}{r_1} \qquad (18) $$

where $K$ is a constant.

If we combine Equations (17) and (18) and reorganise, we see that:

$$ r_1 = \frac{K(q + \delta)}{r_A \chi_A + r_X \chi_X} R = \alpha R \qquad (19) $$

where $K$, $q$, $r_A$, $r_X$, $\chi_A$ and $\chi_X$ are constants which only depend on the system, so that $\alpha$ itself is a constant as long as $\delta$ remains constant. In other words, $r_1$ will scale linearly with the interatomic distance, as long as the charge transfer does not change upon pressurization. Of course, similar equations apply to $r_2$, with $r_2 = \beta R$.

In order to check this hypothesis, we have repeated our DFT-ELF procedure upon compression of the zero pressure reference cell (no optimization needed due to the symmetry of the cell). Two prototypical examples of covalent (BAs, $\Delta\chi = 0.14$) and polar (BN, $\Delta\chi = 1$) solids were chosen to compare their corresponding compression with respect to their bonding type. Figure 5 collects the results upon compression for $r_1, r_2$ and $R_B$. Negligible deviations from linearity are found in both cases, leading to three main conclusions:

* $R_B$ is the quantity showing the greatest compression of the model since it represents the valence size, i.e. the chemically sensitive part as pressure is applied [43,49].

* $\alpha, \beta$ and $\nu$ are constant upon compression for the various types of bonding considered. Further analysis shows that $\alpha, \beta$ and $\nu$ can certainly be related to hardness, since they determine which ion compresses the most. Indeed a look at Figure 5 shows that in the case of BN, where the hardness
of atoms is more similar, $r_1$ and $r_2$ compress at similar rates, whereas in BAs, the compression of As overcomes that of B.

* As a consequence of the previous observation and from Equation (19), we can consider that the charge transfer is also constant upon compression for the ranges considered.

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Evolution of the geometrical parameters, $r_1$, $r_2$ and $R_B$, upon compression of the interatomic distance *R* for BN (**left**) and BAs (**right**). Data relative to the corresponding equilibrium values denoted by a 0 superscript. All numbers in Å.</figcaption></figure>

If we analyze these constants, it is observed that the hardest materials, diamond and BN, not only have the smallest *R* and the smallest $R_B$, but also the highest $\nu = R_B/R$ ratio (see Table 1).

The size of the parameters $r_1$ and $r_2$ (or $\alpha$ and $\beta$) also plays an important role. Since these parameters are related to the relative hardness of atoms, the smaller $\alpha$ and $\beta$ are, the less compressible the whole system becomes. What is most important, both parameters should be similarly small. Otherwise, the softest ion will subsume the pressure. Ideally, for a given $r_1$, the compressibility is minimized when the relationship $\alpha = \beta$ is fulfilled. This condition is equivalent to locating the bond position in the middle of both atoms (see Figure 5 (left)).

These dependencies enable us to establish new relationships for materials with low compressibility, directly derived from microscopic conditions. From Table 1 we see that there are two cases which deserve our attention as far as their $r_1/r_2$ ratio is concerned: diamond and BN represent the two different ways of achieving $\alpha = \beta$ (i.e., $r_1 = r_2$):

* Perfectly covalent case (diamond): in this case $A = X$ (or $\Delta\chi = 0$), so that all lengths are the same.

* Polar case (BN): in this case $A \neq X$ ($\Delta\chi \neq 0$). However, a relationship holds between the core and the bond location displacement. It is found that the cores differential size, $\epsilon = r_X - r_A$ for $\epsilon > 0$ is equivalent to the displacement of the bond location from the center toward the smaller atom A, $r'_1 = r'_2 + \epsilon$ (remember that $r_1 = r'_1 + r_A$ and $r_2 = r'_2 + r_X$). In BN, a good combination of atomic sizes and polarity compensate to provide $r_1/r_2 = 1.03$.

**6. Summary and Conclusions**

We have introduced the Electron Localization Function as the source of the necessary parameters for the analysis of the Bond Charge Model in solid state applications, named NEWBCM. The availability
of charges and bond location has enabled us to develop a new BCM-based model in which, parameters are not only obtained *ab initio*, but they are closer to chemical intuition and recover good agreement with experience.

The ability of NEWBCM to describe compressibility trends among these crystal families has been exploited in order to settle the directions for novel superhard materials design. Some of these hints are not related to macroscopic properties but to microscopic ones. Thus, quantitative relationships between electronic structure and macroscopic properties, such as compressibility, can now be proposed.

The micro-macroscopic properties that guide the quest for superhard materials can be summarized as follows: (i) Small atoms; (ii) Small bond; (iii) High $v$; (iv) Hard atoms; and (v) Atoms of similar hardness. Although most of these guidelines are related (e.g., (i) and (iv)) and extremely intuitive, some of them, such as (iii) would not be deduced at first glance, and provide useful quantitative insight into the search and prediction of novel superhard materials from microscopic grounds.

As an example, it has been shown that the hardest compounds within our study set, diamond and BN, both follow these principles, and more specifically, both show an $r_1/r_2$ ratio very close to one in spite of belonging to two very different bonding types (covalent and very polar, respectively). These conditions thus unite their different macroscopic properties into one relevant microscopic condition for superhardness.

Although this model is only applicable to zinc-blende types of structures, work is in progress to extend it to other structural types. This would, for example, enable us to understand the driving forces for chemical changes in phase transitions.

## Acknowledgments

Financial and computational support from the MALTA-Consolider Ingenio-2010 program under project CSD2007-00045 is acknowledged. JCG and MM thank the Spanish Ministry of Education and Principado de Asturias (FICYT), respectively, for postdoctoral grants. MM, JMM, and JMR thank the Spanish Ministry of Education for financial support (CTQ2009-14596-C02-02 and CTQ2012-38599-C02-02).

## Author Contributions

Julia Contreras-García has carried out some of the calculations lead the research and written the article; Miriam Marqués and José Manuel Menéndez have carried out some of the calculations; José Manuel Recio has supervised the work.

## Conflicts of Interest

The authors declare no conflict of interest.

## References

1. Gilman, J.J. *Electronic Basis of the Strength of Materials*; Cambridge University Press, Cambridge, UK, 2003.
2. Gilman, J.J. Why silicon is hard. *Science* **1993**, *261*, 1436–1439.

3. Gao, F.; He, J.; Wu, E.; Liu, S.; Yu, D.; Li, D.; Zhang, S.; Tian, Y. Hardness of covalent crystals. *Phys. Rev. Lett.* **2003**, *91*, 015502.

4. Li, K.; Wang, X.; Zhang, F.; Xue, D. Electronegativity identification of novel superhard materials. *Phys. Rev. Lett.* **2008**, *100*, 235504.

5. Oganov, A.R.; Glass, C.W. Crystal structure prediction using *ab initio* evolutionary techniques: Principles and applications. *J. Chem. Phys.* **2006**, *124*, 244704.

6. Lyakhov, A.O.; Oganov, A.R. Evolutionary search for superhard materials: Methodology and applications to forms of carbon and $TiO_2$. *Phys. Rev. B* **2011**, *84*, 092103-1-4.

7. Parr, R.; Pearson, R.G. Absolute hardness: Companion parameter to absolute electronegativity. *J. Am. Chem. Soc.* **1983**, *105*, 7512–7516.

8. Putz, M.V. Electronegativity: Quantum observable. *Int. J. Quantum Chem* **2009**, *109*, 733–738.

9. Putz, M.V. Chemical action and chemical bonding. *J. Mol. Stru.: Theochem* **2009**, *900*, 64–70.

10. Putz, M.V. The Bondons: The quantum particles of the chemical bond. *Int. J. Mol. Sci.* **2010**, *11*, 4227–4256.

11. Putz, M.V. Density functional theory of bose-einstein condensation: Road to chemical bonding quantum condensate. *Struc. Bond.* **2012**, *149*, 1–50.

12. Bader R.F.W. *Atoms in Molecules, A Quantum Theory*; Oxford University Press: Oxford, UK, 1990.

13. Martín Pendás, A.; Costales, A.; Luaña, V. The topology of the electron density in ionic materials. III. Geometry and ionic radii. *J. Phys. Chem. B* **1998**, *102*, 6937–6948.

14. Becke, A.D.; Edgecombe, K.E. A simple measure of electron localization in atomic and molecular systems. *J. Chem. Phys.* **1990**, *92*, 5397–5403.

15. Savin, A.; Jepsen, O.; Flad, J.; Andersen, L.K.; Preuss, H.; von Schnering, H.G. Electron localization in solid-state structures of the elements: The diamond structure. *Angew. Chem. Int. Ed. Engl.* **1992**, *31*, 187–188.

16. Borkman, R.F.; Parr, R.G. Toward an understanding of potential-energy functions for diatomic Molecules. *J. Chem. Phys.* **1968**, *48*, 1116–1126.

17. Parr, R.G.; Borkman, R.F. Simple bond-charge model for potential-energy curves of homonuclear diatomic molecules. *J. Chem. Phys.* **1968**, *49*, 1055–1058.

18. Borkman, R.F.; Simons, G.; Parr, R. Simple bond-charge model for potential-energy curves of heteronuclear diatomic molecules. *J. Chem. Phys.* **1969**, *50*, 58–65.

19. Parr, R.; Yang, W. *Density Functional Theory of Atoms and Molecules*; Oxford University Press: Oxford, UK, 1989.

20. Ray, N.K.; Samuels, L.; Parr, R.G. Studies of electronegativity equalization. *J. Chem. Phys.* **1979**, *70*, 3680–3684.

21. Martin, R.M. A simple bond charge model for vibrations in covalent crystals. *Chem. Phys. Lett.* **1968**, *2*, 268–270.

22. Oviedo's Quantum Chemistry Group. Available online: http://azufre.quimica.uniovi.es/software.html (accessed on 8 April 2015).

23. Pauling, L. *The Nature of the Chemical Bond*; 3rd ed.; Cornell Univ. Press: Ithaca, NY, USA, 1960.
24. Putz, M.V. Markovian approach of the electron localization functions. *Int. J. Quantum Chem.* **2005**, *105*, 1–11.

25. Putz, M.V. Density functionals of chemical bonding. *Int. J. Mol. Sci.* **2008**, *9*, 1050–1095.

26. Gillespie, R.J.; Nyholm, R.S. Inorganic stereo-chemistry. *Quart. Rev.* **1957**, *11*, 339–380.

27. ELK Code. Available online: http://elk.sourceforge.net/ (accessed on 8 April 2015).

28. Drief, F.; Tadjer, A.; Mesri, D.; Aourag, H. First principles study of structural, electronic, elastic and optical properties of MgS, MgSe and MgTe. *Catal. Today* **2004**, *89*, 343–355.

29. Semiconductors information web-site. Available online: www.semiconductors.co.uk (accessed on 8 April 2015).

30. Perdew J.P.; Wang, Y. Accurate and simple analytic representation of the electron-gas correlation energy. *Phys. Rev. B* **1992**, *45*, 13244–13249.

31. Perdew, J.P.; Burke K.; Ernzerhof, M. Generalized gradient approximation made simple. *Phys. Rev. Lett.* **1996**, *77*, 3865–3868.

32. Contreras-García, J.; Martin Pendás, A.; Silvi, B.; Recio, J.M. Useful applications of the electron localization function in high-pressure crystal chemistry. *J. Phys. Chem. Solids* **2008**, *69*, 2204–2207.

33. Contreras-García, J.; Silvi, B.; Martín Pendás, A.; Recio, J.M. Computation of local and global properties of the electron localization function topology in crystals. *J. Chem. Theory Comput.* **2009**, *5*, 164–173.

34. Kohout, M.; Savin, A. Atomic shell structure and electron numbers. *Int. J. Quant. Chem.* **1996**, *60*, 875–882.

35. Frisch, M.J.; Trucks, G.W.; Schlegel, H.B.; Scuseria, G.E.; Robb, M.A.; Cheeseman, J.R.; Montgomery, J.A., Jr.; Vreven, T.; Kudin, K.N.; Burant, J.C.; et al. Gaussian 03, Revision C.02; Gaussian, Inc.: Wallingford, CT, USA, 2004.

36. Noury, S.; Krokidis, X.; Fuster, F.; Silvi, B. *TopMoD Package*; Universite Pierre et Marie Curie: Paris, France, 1997.

37. Otero-de-la-Roza, A.; Blanco, M.A.; Martín Pendás, A.; Luaña, V. Critic: A new program for the topological analysis of solid-state electron densities. *Comput. Phys. Commun.* **2009**, *180*, 157–166.

38. Kohout, M.; Savin, A.; Preuss, H. Contribution to electron density analysis. I. Shell structure of atoms. *J. Chem. Phys.* **1991**, *95*, 1928–1942.

39. Ackland, G. J.; Marqués, M.; Contreras-García, J.; McMahon, M.I. Origin of incommensurate modulations in the high-pressure phosphorus IV phase. *Phys. Rev. B* **2008**, *78*, 054120.

40. Marqués, M.; Ackland, G.J.; Lundegaard, L.F.; Contreras-García, J.; McMahon, M.I. Potassium under pressure: A pseudobinary ionic compound. *Phys. Rev. Lett.* **2009**, *103*, 115501.

41. Marqués, M.; Santoro, M.; Guillaume, C. L.; Gorelli, F.; Contreras-García, J.; Howie, R.; Goncharov, A.F.; Gregoryanz, E. Optical and electronic properties of dense sodium. *Phys. Rev. B* **2011**, *83*, 184106.

42. Anderson, D.L.; Anderson, O.L.J. The bulk modulus-volume relationship for oxides. *J. Geophys. Res.* **1970**, *75*, 3494–3500.
43. Contreras-García, J.; Mori-Sánchez, P.; Silvi, B.; Recio, J.M.A. Quantum chemical interpretation of compressibility in solids. *J. Chem. Theory Comput.* **2009**, *5*, 2108–2114.

44. Glasser, L. Volume-based thermoelasticity: Compressibility of inorganic solids. *Inorg. Chem.* **2010**, *49*, 3424–3427.

45. Glasser, L. Volume-based thermoelasticity: Compressibility of mineral-structured materials. *J. Phys. Chem. C* **2010**, *114*, 11248–11251.

46. Glasser, L.; Jenkins, H.D.B. Internally consistent ion volumes and their application in volume-based thermodynamics. *Inorg. Chem.* **2008**, *47*, 6195–6202.

47. Martín Pendás, A.; Costales, A.; Blanco, M.A.; Recio, J.M.; Luaña, V. Local compressibilities in crystals. *Phys. Rev. B* **2000**, *62*, 13970–13978.

48. Ouahrani, T; Menendez, J.M.; Marqués, M.; Contreras-García, J.; Baonza, V.G.; Recio, J.M. Local pressures in Zn chalcogenide polymorphs. *Europhys. Lett.* **2012**, *98*, 56002.

49. Contreras-García, J.; Recio, J.M. On bonding in ionic crystals. *J. Phys. Chem. C* **2011**, *115*, 257–263.

© 2015 by the authors; licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution license (http://creativecommons.org/licenses/by/4.0/).