Bulk modulus for polar covalent crystals.

A microscopic empirical model of bulk modulus based on atomic-scale parameters is proposed. These parameters include the bond length, the effective bonded valence electron (EBVE) number, and the coordination number product of two bonded atoms, etc. The estimated bulk moduli from our model are in good agreement with experimental values for various polar covalent crystals including ionic crystals. Our current work sheds lights on the nature of bulk modulus, provides useful clues for design of crystals with low compressibility, and is applicable to complex crystals such as minerals of geophysical importance.

In the past several decades, great efforts have been made to establish direct correlations between the macroscopic mechanical properties and the microscopic parameters of solid. Such correlations can reveal the crucial factors that determine materials mechanical properties. More importantly, they can be used to estimate mechanical properties of designed crystals and provide insightful information before experimental trials. Up to now, some important progresses have been achieved. For example, microscopic hardness models have been constructed in terms of atomic-scale parameters and used extensively experimentally as well as in superhard materials design1. However, microscopic models for other mechanical properties, such as bulk modulus234, are still in the primary stage.

Bulk modulus measures material's resistance to uniform compression, and is highly related to the chemical composition and crystal structure. Materials with high bulk modulus are potentially superhard materials according to the rough correlation between bulk modulus and hardness5. Moreover, bulk modulus is relevant to geophysics since it is involved in the interpretation of earth's seismic properties. Numerous experimental and theoretical efforts on bulk modulus have thus been undertaken. However, the measurements of the bulk moduli are nontrivial and values for many materials are still unavailable or unreliable. Concurrently, theoretical calculations of bulk moduli are performed mainly with two approaches. The first is based on ab initio techniques. Unlike hardness, bulk modulus is a strictly defined thermodynamic quantity (B = −VdP/dV), and can thus be obtained theoretically by either using the stain–stress method6, or fitting the calculated total energy with Murnaghan equation of state7 and Birch–Murnaghan equation of state8. These ab initio calculations, relying on computational facilities, can also provide details of the bonding, structural, and electronic properties of solids9. The second is empirical approaches based on readily accessible parameters, e.g. chemical bond length, valence, and ionicity. The empirical models can estimate bulk moduli of crystals with the advantage of simplicity and comparable accuracy as ab initio ones, but usually different quantification schemes are being used for different class of crystals.

To put present work in context, a brief review of empirical approaches for bulk modulus estimation is necessary. In 1923, Bridgman proposed an empirical relationship between bulk modulus and molar volume for metals10. Sequentially, researchers revealed the bulk moduli of minerals are controlled primarily by the specific volume per ion pair111213. Jayaraman correlated the bulk modulus with the unit cell volume and the effective covalence product for rocksalt-structured rare-earth compounds14. It was established that, for a given class of materials with identical crystal structure and similar bonding properties, the bulk modulus would scale with the unit cell volume according to , where B0 and p, constants characteristic of selected materials family, are mainly determined by the type of chemical bond (e.g. bond ionicity) and the dependence of the interatomic potential on the interatomic distance, respectively15. Hazen and Finger proposed a more general bulk modulus–volume relationship for cation coordination polyhedral in a variety of structure types including oxides, silicates, halides, sulfides, phosphides, and carbides: , where B, S2, , Z and Z are the polyhedral bulk modulus, an empirical ionicity, the mean cation–anion separation, cation and anion formal charge, respectively16. However, application of this model to rutile-type oxides showed large deviations from experimental values as pointed out by the same authors17.

In his seminal work for zincblende semiconductors, Cohen assumed the bulk modulus of covalent crystals should obey the scaling relation between covalent bond energy and the bond volume change, and proposed a bulk modulus model using bond length and an ionicity parameter2. The different power-law dependent of bulk modulus on the bond length for semiconductors (B ∝ d−3.5) and ionic crystals (B ∝ d−3) were also pointed out2. Most recently, Cohen's work was developed by Kamran3 and Verma18 for diamond-like and zincblende covalent crystals. In both cases, the d−3.5 dependence of B is retained with different descriptions of the ionicity.

Other empirical bulk modulus models have also been proposed based on different considerations, such as the electronegativity4, transition pressure19, Debye temperature20, plasmon energy21, partitioning of the elastic constants with the atomic basin concept22, etc. Details of these models are less relevant to present work. It should be noted that some of the key variables used in these models are undeducible from the crystal structural information and have to be provided otherwise, which limit these models' applicability.

The bulk modulus of a crystal microscopically depends on the nature of its chemical bonds such as bond length and type, as corroborated in the empirical relationships between bulk modulus and specific volume (or bond length). Also, a crystal structure dependence of bulk modulus was revealed42324. To formulate a microscopic bulk modulus model with easy accessibility and broad applicability, input variables should be directly deducible for selected crystals. Parameterizing of the influences of bond type and crystal structure is thus prerequisite. Recently, we suggested a universal quantification of the bond strength in covalent and ionic crystals by using the proposed effectively bonded valence electron (EBVE) number, n, of a chemical bond2526, which can be used as a parameter to characterize bond type (See Method section for the details). For single-bond crystals, Levine pointed out the normalized bond volume, ν = n3/Ω, can be used as a measure of crystal geometry, where n, d, and Ω are the bond number in the unit cell, bond length, and unit cell volume, respectively27. Further analysis shows ν is proportional to the product of coordination numbers for two bonded atoms, and distinguishes zincblende, wurtzite, rocksalt, CsCl, cuprite, fluorite, rutile, and corundum structures (Figure 1). This product can thus be used to characterize the crystal structure. We emphasize that both EBVE number and coordination number product can easily be calculated for a designated crystal and applied to multi-bond crystals.

Figure 1

Plot of ν against the product of coordination numbers of two bonded atoms.

Thus far, three main parameters affecting the bulk modulus, i.e. the bond length (d), EBVE number (n) and coordination number product (p = N), are determined. Previous studies revealed the power-law behavior of the bond length (or unit cell volume) on bulk modulus215. Also, the exponential dependence of the bond strength on EBVE number is emphasized25. These studies highlighted the roles of the bond length and valance electrons on bulk modulus, which can be quantified tentatively in terms of fundamental variables with the formalism, where C, l, m, and k are the fitting constants. In the following, we deduce the empirical expression of bulk modulus for simple ANB8-N type covalent crystals with one coordination number. The generalizations to covalent crystals with diverse coordination numbers and to ionic crystals are presented sequentially, followed by generalization to multi-bond crystals. The application of our model to geophysics interesting minerals is also demonstrated.

Results

ANB8-N type covalent crystals

We first consider single-bond ANB8-N type covalent crystals with zincblende, wurtzite, and rocksalt structures where only one type of chemical bond and one coordination number are presented. With the aid of available experimental bulk modulus data (shown in Table 1), our approach is to decompose the bulk modulus into parameters associated with individual chemical bonds, i.e. d, n, and p, by using Eqn. 1. The constants, C, l, m, and k, can be extrapolated from the experimental bulk modulus data for these covalent crystals with the Levenberg-Marquardt method2829, giving The key variables, experimental bulk moduli341630313233343536, and calculated bulk moduli from our model and from other models34 are listed in Table 1 for selected ANB8-N type covalent crystals. The comparison between calculated bulk moduli and experimental values is emphasized in Figure 2. The overall accuracy of our model is comparable with the electronegativity model4 and better than that of Kamran et al3. Unlike the other two models, where either electronegativity or ionicity has to be determined through other way, all the input parameters in our model can be directly deduced from the designated crystal with simple arithmetic calculations, providing a unique advantage for its application. To further check the performance of Eqn. 2 for other similar crystals not included in Table 1, we list the results for zincblende-structured BAs37 and BSb38, rocksalt-structured PbS39, PbTe40, and AgCl41 in Table 2. The consistence of bulk moduli from the experiments (ab initio calculations for BSb) and our empirical model is satisfying.

Table 1

Parameters related to the calculation of bulk modulus as well as experimental and calculated bulk moduli for simple ANB8-N covalent crystals

Crystal	Stucture	d	nAB	p	BExp	BCal	BCala	BCald
C	A4	1.55	0.707	16	442a	433	418.2	442.1
Si	A4	2.35	0.707	16	100a	103	97.4	95.6
Ge	A4	2.45	0.707	16	78a	89	84.2	81.3
Sn	A4	2.81	0.707	16	55b	55	52.1	52.6
SiC	B3	2.81	0.707	16	211a	222	202.9	203.6
BN	B3	1.57	0.643	16	369a	377	373	398.9
BP	B3	1.97	0.643	16	173a	172	180.4	170.8
AlN	B4	1.9	0.643	16	201a	195	180.9	196.6
AlP	B3	2.36	0.643	16	86a	92	88.3	89.9
AlAs	B3	2.43	0.643	16	77a	83	80.5	79.1
AlSb	B3	2.66	0.643	16	59a	61	56.1	61.3
GaN	B4	1.97	0.643	16	190a	172	164	193.8
GaP	B3	2.36	0.643	16	89a	92	86.6	88.8
GaAs	B3	2.45	0.643	16	76a	81	77.4	78.6
GaSb	B3	2.65	0.643	16	57a	62	59.6	62.1
InN	B4	2.16	0.643	16	137a	125	111.1	129.4
InP	B3	2.54	0.643	16	73a	71	66.1	69.3
InAs	B3	2.61	0.643	16	60a	65	61.2	62.8
InSb	B3	2.81	0.643	16	47a	50	47.7	51.4
BeO	B4	1.64	0.474	16	244c	252
BeS	B3	2.1	0.474	16	105d	107	132.7	112.9
BeSe	B3	2.2	0.474	16	92d	91	113.2	95.1
BeTe	B3	2.4	0.474	16	67d	68	86.5	72.2
MgO	B1	2.1	0.316	36	163e	178		154.8
MgS	B1	2.6	0.316	36	80f	85	72.2	76.6
CaO	B1	2.4	0.316	36	114g	112		98.8
SrO	B1	2.58	0.316	36	88g	87		78.7
BaO	B1	2.77	0.316	36	69h	68		62.5
ZnO	B4	1.97	0.474	16	139d	134		141
ZnS	B3	2.34	0.474	16	75i	74	82.8	76.6
ZnSe	B3	2.45	0.474	16	62d	63	69.3	66.0
ZnTe	B3	2.64	0.474	16	51d	49	55.6	53.4
CdS	B4	2.52	0.474	16	62d	57	62.6	60.7
CdSe	B4	2.62	0.474	16	56d	50	54.4	53.7
CdTe	B3	2.81	0.474	16	42d	39	42.9	44.0
HeSe	B3	2.63	0.474	16	49d	49	55.1	53.0
HgTe	B3	2.78	0.474	16	44d	41	46.1	44.6
MnO	B1	2.22	0.316	36	147j	147
FeO	B1	2.17	0.316	36	174j	159

aRef. 3, b Ref. 30, c Ref. 31, d Ref. 4, e Ref. 32, f Ref. 33, g Ref. 34, h Ref. 16, i Ref. 35, j Ref. 36.

Figure 2

Calculated bulk moduli from different models versus the experimental values for simple ANB8-N type covalent crystals.

The inset shows magnified image between 30 and 135 GPa.

Table 2

Bulk moduli determined from experiments (ab initio calculations for BSb) and Eqn. 2 for zincblende-structured BAs and BSb, rocksalt-structured PbS, PbTe, and AgCl

Crystal	Stucture	d	nAB	p	BExp	BCal
BAs	B3	2.07	0.643	16	148a	145
BSb	B3	2.22	0.643	16	116b	114
PbS	B1	2.97	0.316	36	50c	54
PbTe	B1	3.23	0.316	36	38d	40
AgCl	B1	2.89	0.165	36	47e	47

aRef. 37, b Ref. 38, c Ref. 39, d Ref. 40, e Ref. 41.

AnBm type covalent crystals

For selected AnBm () type crystals (Table 3), bulk moduli calculated from Eqn. 2 are systematically smaller than experimental values1317224243. It should be noted that an average bond length was used for rutile and corundum structures in Table 2. Previously we chose p = N as a characteristic of crystal geometry, which did not explicitly reveal the diversity of the coordination numbers as presented in AnBm crystals. Here we define Nmax = max[N, N] and Nmin = min[N, N]. We choose two schemes to address this issue. The first one is to use N and N individually as input parameters, with the power-law indexes determined by fitting the experimental data. The second one is to introduce an asymmetry parameter, Λ = Nmax/Nmin, for two bonded atoms with different coordination numbers. The bulk modulus formulae are and B(GPa) = CpΛexp(kn), respectively. To keep the consistency with Eqn. 2, only the power-law indexes to Nmax, Nmin, and Λ are adjusted during data fitting, and determined to be 1.201, 0.672, and 0.370, respectively. These two schemes give similar results (Bcal,1 and Bcal,2 in Table 3). Here we choose the second scheme, Thus far, a uniform empirical bulk modulus model for both ANB8-N type and AnBm type covalent crystals is formulated.

Table 3

Parameters as well as experimental and calculated bulk moduli for AnBm covalent crystals

Crystal	Stucture	d	nAB	p	Λ	BExp	BCal, 1	BCal, 2
Ag2O	C3	2.044	0.474	8	2	84a	79	81
Cu2O	C3	1.85	0.474	8	2	112a	111	114
ZrO2	C1	2.195	0.474	32	2	235b	223	222
TiO2	C4	1.96	0.596	18	2	224c	233	235
SnO2	C4	2.05	0.596	18	2	205d	200	201
GeO2	C4	1.878	0.596	18	2	258e	269	271
Al2O3	D51	1.915	0.474	24	1.5	252c	253	246
Fe2O3	D51	2.025	0.474	24	1.5	207c	209	203

aRef. 22, b Ref. 42, c Ref. 13, d Ref. 43, e Ref. 17.

Ionic crystals

In examining the bulk moduli for I-VII rocksalt ionic crystals, Cohen gave an approximate scaling, B = 550d−3, which differs from that for diamond- and zincblende-structured covalent crystals in the prefactor and power-law index2. We use Eqn. 1 to fit the experimental bulk moduli of rocksalt ionic crystals23, using fixed l (0.914) and k (1.485). The power-law index m determined here (−3.43) is very close to the previous value (−3.46). For simplicity, an identical value of −3.46 is used, and the empirical bulk modulus formula for rocksalt-structured ionic crystals is deduced as We apply Eqn. 4 to cesium halides with CsCl-structure and produce bulk moduli in great consistent with experimental ones23, as listed in Table 4.

Table 4

Parameters as well as experimental and calculated bulk moduli for ionic crystals

Crystal	Stucture	d	nAB	p	BExp	BCal
LiF	B1	2.01	0.165	36	76.4a	76.7
LiCl	B1	2.57	0.165	36	32.9a	32.8
LiBr	B1	2.75	0.165	36	26.0a	25.9
LiI	B1	3.01	0.165	36	19.3a	19.0
NaF	B1	2.31	0.165	36	47.1a	47.4
NaCl	B1	2.82	0.165	36	23.9a	23.8
NaBr	B1	2.99	0.165	36	19.6a	19.4
NaI	B1	3.24	0.165	36	14.9a	14.7
KF	B1	2.67	0.165	36	28.7a	28.7
KCl	B1	3.15	0.165	36	16.5a	16.2
KBr	B1	3.3	0.165	36	13.9a	13.8
KI	B1	3.53	0.165	36	11.1a	10.9
RbF	B1	2.83	0.165	36	24.1a	23.5
RbCl	B1	3.29	0.165	36	14.0a	13.9
RbBr	B1	3.44	0.165	36	12.0a	12.0
RbI	B1	3.67	0.165	36	9.6a	9.6
CsF	B1	3.00	0.165	36	19.3a	19.2
CsCl	B2	3.57	0.123	64	17.6a	16.7
CsBr	B2	3.72	0.123	64	15.4a	14.5
CsI	B2	3.96	0.123	64	12.5a	11.7
Li2O	C1	2.01	0.237	32	88b	100
CaF2	C1	2.365	0.247	32	84.0c	85.6
SrF2	C1	2.5	0.247	32	71.3c	70.7
BaF2	C1	2.683	0.247	32	55.5c	55.3

aRef 23, b Ref. 22, c Ref. 44.

Now we consider ionic crystals, such as Li2O, CaF2, SrF2, and BaF2, which possess one type of chemical bond and two distinct coordination numbers. At first glance, the bulk modulus for these crystals should be calculated in the same way as for AnBm type covalent crystals by introducing Λ0.37 term into Eqn. 4. However, the calculated bulk moduli are systematically smaller than experimental values except for Li2O2244. The ratios of Bcal to Bexp for CaF2, SrF2, and BaF2 are approximately the same, indicating a common factor we may miss in the bulk modulus formula for ionic crystals. We note an obvious difference in cationic charge state of Li (+1) and alkali-earth metal atoms (+2) in these crystals. The valence electrons in ionic crystals, unlike those in covalent crystals distributed on the bond, are transferred from cations to anions, which weakens the resistance to compression and results in smaller bulk modulus. It is plausible that cations with higher charge state have a stronger tendency to pull the valence electrons back to the bond region, thus improve the resistibility of chemical bonds to compression and increase the bulk modulus. Therefore, we introduce an additional term, Q′, to Eqn. 4 to account for such effect, where Q is the cationic charge state and r is the power-law index determined to be 0.58 from CaF2, SrF2, and BaF2 data. The bulk modulus formula for ionic crystals can be defined as The key parameters, experimental and calculated bulk moduli are listed in Table 4 for selected ionic crystals. Again, a great consistency is achieved between the calculated bulk moduli and the experimental ones.

Multi-bond crystals

Up to now, we have formulated the bulk modulus for simple polar covalent crystals with parameters deduced from the crystal structure. The generalization of our bulk modulus model to complex multi-bond crystals is straightforward, since all the parameters used in our model are directly correlated to the chemical bond and can easily be applied to complex multi-bond crystals. The only issue is how to weight the contributions to bulk modulus from individual constitutional chemical bonds. Bulk modulus for a multicomponent material can be attributed to the contributions from individual components. From the definition of bulk modulus, B = −VdP/dV, we get where B and v = V/V are the bulk modulus and volume fraction of the μ component, respectively. For a multi-bond crystal, B can be determined by using previous equations, and the volume fraction of distinct chemical bond can be calculated as , where d and N are the bond length and bond multiplicity of the μ-type bond in the unit cell, respectively. Bulk modulus for a multi-bond crystal can then be calculated with Eqn. 6.

Key variables and bulk moduli are listed in Table 5 for representative multi-bond crystals132245464748. The global consistency between calculated bulk modulus and experimentally determined value are remarkably good considering the very simple assessment method and the relatively complex structure in multi-bond crystals. Here we calculate the bulk modulus of spinel Fe3O4 to show how our model works for multi-bond systems. In Fe3O4, all oxygen atoms are 4-coordinated, while one third of Fe atoms, Fe1, are 4-coordinated and the other two thirds of Fe atoms, Fe2, are 6-coordinated. The bond multiplicities of O-Fe1 and O-Fe2 bonds are 4 and 12, respectively. B for O-Fe1 and O-Fe2 components determined from Eqn. 3 are 144 and 188 GPa, respectively. Using Eqn. 6, the bulk modulus of spinel Fe3O4 is calculated to be 177 GPa, in good consistent with the experimental value of 187 GPa13. Also shown in Table 5 are the theoretically proposed β-C3N446 and z-BC2N47 structures. The calculated bulk moduli for these structures with Eqn. 3 and 6 are consistent with the values from first-principles calculations, indicating the prediction power of our simple semi-empirical model.

Table 5

Parameters as well as experimental and calculated bulk moduli for multi-bond crystals

Crystal	Bond	Nμ	d	nAB	p	Λ	Bμ	BCal	BExp
Fe3O4	Fe-O	4	1.927	0.474	16	1	144	178	187a
	Fe-O	12	2.073	0.474	24	1.5	188
SrTiO3	Ti-O	6	1.95	0.555	36	1	327	147	174b
	Sr-O	12	2.76	0.165	72	2	134
MgAl2O4	Mg-O	4	1.88	0.474	16	1	157	210	202a
	Al-O	12	1.95	0.474	24	1.5	233
β-Si3N4	Si-N	6	1.73	0.707	12	1.333	252.9	251	250c
	Si-N	6	1.728	0.707	12	1.333	253.9
	Si-N	6	1.704	0.707	12	1.333	266.5
	Si-N	6	1.767	0.707	12	1.333	235
c-Si3N4	Si-N	1	1.78	0.781	16	1	298	308	300c
	Si-N	3	1.88	0.588	24	1.5	311
β-C3N4	C-N	6	1.457	0.707	12	1.333	458	466	451d
	C-N	6	1.447	0.707	12	1.333	469
	C-N	6	1.449	0.707	12	1.333	467
	C-N	6	1.447	0.707	12	1.333	469
z-BC2N	B-C	4	1.559	0.6	16	1	361	411	403e
	C-C	4	1.527	0.707	16	1	456
	C-C	8	1.537	0.707	16	1	445
	B-N	8	1.546	0.643	16	1	397
	B-N	4	1.579	0.643	16	1	369
	C-N	4	1.542	0.781	16	1	491
SiO2	Si-O	4	1.76	0.596	18	2	340	329	314f
	Si-O	2	1.808	0.596	18	2	310

aRef. 13, b Ref. 22, c Ref. 45, d Ref. 46, e Ref. 47, f Ref. 48.

Discussion

We can rewrite Eqn. 5 as where Q = 0.464Q0.58 (Q is the cationic charge state) is an ionic term which need to be considered for ionic crystals. We thus formulate our empirical bulk modulus for polar crystals with Eqn. 7 and 6. In Figure 3, bulk moduli from first principles calculations and our empirical model are compared with respect to the experimental values for some ANB8-N type covalent crystals. Obviously, the accuracy of our model is comparable to LDA-U results and better than LDA-N, GGA-U, and GGA-N results. The advantage of our empirical model is obvious for saving a lot of calculation endeavor compared with these first-principles calculation methods.

Figure 3

Comparison of bulk moduli deduced from first-principles calculations and our empirical models with respect to the experimental values.

The inset shows the corresponding results of α-, β-, and γ-Mg2SiO4.

It is interesting to further test our model to some geophysical interested materials. These materials usually have complex crystal structures and demonstrate polymorphism, which are difficult to be treated even with state of the art calculation methods. In Table 6 we list results for polymorphous Mg2SiO4, i.e. forsterite (α)49, wadsleyite (β)50, and ringwoodite (γ) Mg2SiO451, based on our model. These materials are the main constituents of the middle and lower crust as well as upper mantle of the earth. The knowledge of the elastic properties of them is of great importance to understand earth's seismic properties. Our bulk modulus estimations, as listed in Table 6, clearly reflex the relative values of different phases of Mg2SiO4 with an acceptable accuracy. We compare the bulk moduli of α-, β-, and γ-Mg2SiO4 from first principles calculations and our empirical model with respect to the experimental values (See the inset to Figure 3). Although the empirical model estimates bulk modulus with accuracy inferior to those from first principles calculations, the divergence is in the range of experimental error. Considering the simple formula of our model, it can thus be used to quickly estimate bulk modulus for complex crystals.

Table 6

Parameters, experimental and calculated bulk modulus for different Mg2SiO4 phases

Crystal	Bond	Nμ	d	nAB	p	Λ	Bμ	BCal	BExp
α-Mg2SiO4	Mg-O	16	2.131	0.325	24	1.5	137	153	129a
	Mg-O	16	2.084	0.325	24	1.5	148
	Mg-O	16	2.068	0.325	24	1.5	152
	Mg-O	16	2.065	0.325	24	1.5	152
	Mg-O	16	2.21	0.325	24	1.5	121
	Mg-O	8	2.177	0.325	24	1.5	127
	Mg-O	8	2.045	0.325	24	1.5	158
	Si-O	16	1.637	0.832	16	1	431
	Si-O	8	1.613	0.832	16	1	453
	Si-O	8	1.655	0.832	16	1	415
β-Mg2SiO4	Mg-O	16	2.046	0.325	24	1.5	158	163	172b
	Mg-O	8	2.115	0.325	24	1.5	141
	Mg-O	16	2.093	0.325	24	1.5	146
	Mg-O	4	2.095	0.325	24	1.5	146
	Mg-O	4	2.035	0.325	24	1.5	161
	Mg-O	16	2.123	0.325	24	1.5	139
	Mg-O	16	2.128	0.325	24	1.5	138
	Mg-O	16	2.016	0.325	24	1.5	166
	Si-O	8	1.638	0.832	16	1	430
	Si-O	16	1.632	0.832	16	1	436
	Si-O	8	1.701	0.832	16	1	378
γ-Mg2SiO4	Mg-O	96	2.066	0.325	24	1.5	153	168	185b
	Si-O	32	1.665	0.832	16	1	407

aRef. 49, b Ref. 50, c Ref. 51.

Some cautions must be exercised during applying our bulk modulus model. First, three dimensional bond network with clearly defined bonds are needed since our model is based on the partitioning of bulk modulus into individual bonds. Our empirical model thus does not work for transition metal (TM) borides where the chemical bonds between TM atoms and boron atoms are not well defined. Second, our model cannot be applied to crystals such as α-SiO2 and anatase TiO2 with an open-packed crystal structure. In these structures, bonds may not deform iso-structurally under compression: bond bending and/or non-centrosymmetric forces play important roles13. These structures can be discriminated by the normalized bond volume (ν). As shown in Figure 1, ν for α-SiO2 and anatase TiO2 are obviously shifted away from the fitting straight line, indicating an open-packed crystal structure. Third, we note substantial deviations between the calculated bulk moduli from our model and experimental values for TM carbides and nitrides. Unlike TM oxides which are dominated with ionic character, TM carbides and nitrides are characterized with covalent bonds. The difficulty to make an effective assessment for TM carbides and nitrides mainly arises from this covalent character, with additional complexity from the diverse experimental values reported even for the same material. It is well known that the bond strength of s-p-d hybridized chemical bond is greater than that of s-p hybridized chemical bond52. In addition, f-electrons would push out the spatial extent of the valence d-orbitals, allow d-electrons to make stronger directional bonds, and enhance bulk modulus53. Both effects are not easy to quantify within our simple model. Finally, it should be noted that our microscopic model is deduced based on the experimental bulk modulus data obtained at room temperature. Accordingly, our model works at room temperature and the temperature effect on bulk modulus is out of the scope of current work5455.

To end this paper, we emphasize that our model only need input variables which can be directly deduced from the crystal structure, providing us a general yet powerful tool for bulk modulus evaluation. Pivotal parameters determining bulk modulus are clearly revealed from our model: Short bond length, high EBVE number, and large coordination number are preferred for achieving high bulk modulus. Moreover, the contributions from individual chemical bonds to bulk modulus can readily be partitioned in our model, which is important for design of materials with low compressibility.

Methods

During compression, it is the disturbed valence electrons that determine the bulk modulus. To establish an effective quantification model, we must find a practical way to estimate the population of these electrons. Considering two atoms, A and B, forming a bond in a crystal, the valence electrons are Z and Z with coordination numbers of N and N, respectively. The EBVE number, n, of A-B bond in terms of n and n as where n = Z/N and n = Z/N are the nominal valence electrons contributed to A-B bond. The EBVE numbers of diamond (0.707) and NaCl (0.163) are in good consistent with the Mulliken overlap population from first-principles calculations. It should be pointed out that lone pair electrons must be excluded during calculations of the nominal valence electrons contributed to A-B bond. For example, in rutile TiO2, oxygen atom is 3-coordinated with two electrons forming a lone pair (not participating bonding). no is thus 4/3 instead of 2.

For the first-principles calculations, crystal structural relaxations and property calculations were performed using density functional theory within both generalized gradient approximation (GGA) and local density approximation (LDA) as implemented in the CASTEP code56. Both ultrasoft (U) and norm-conserving (N) pseudopotentials were employed in the calculations. k-point samplings in the Brillouin zone were chosen by using the Monkhorst-Pack scheme with a resolution of 2π × 0.04 Å. The elastic constants, Cij, were calculated theoretically based on the stress and strain relation (Hooke's law), then the bulk moduli for polycrystals were estimated using the Voigt-Reuss-Hill approximation which is the arithmetic average of the upper (Voigt) and lower (Reuss) bounds for the actual macroscopic effective elastic constants6. The calculated bulk moduli are labeled with GGA-U, GGA-N, LDA-U, and LDA-N according to the selected electron-electron exchange functional and pseudopotential.

Author Contributions

B.X. and Y.J.T. designed and coordinated the overall study. Q.Q.W. performed the first principles calculations. B.X. and Y.J.T. analyzed all data and wrote the paper. All authors contributed to the discussion of the results.