ARTICLE pubs.acs.org/JPCA
Topological Partition of the Elastic Constants of Crystals Alberto Otero-de-la-Roza and Víctor Lua~na* Departamento de Química Física y Analítica, Facultad de Química, Universidad de Oviedo, E-33006 Oviedo, Spain ABSTRACT: We present a partitioning of the elastic constants of a crystal into atomic contributions by using the atomic basin concept inherent to Bader’s Quantum Theory of Atoms in Molecules. The partition is made by following the evolution of the cell volume and the atomic basin volumes under appropriately defined cell deformations. The method is carefully examined, including internal consistency checks. The transferability of atomic contributions between different crystals is determined by obtaining and comparing the oxygen contribution to the elastic constants of a selection of cubic oxides that includes the rock-salt, perovskite, antifluorite, and cuprite crystal families.
1. INTRODUCTION Solid state and molecular physics register a long and fruitful tradition of determining the properties of a compound as the additive result of atomic and chemical group contributions.1,2 Bader’s Quantum Theory of Atoms in Molecules3,4 (QTAIM) provides a natural framework to explain and predict this partitioning. According to QTAIM, electronic systems are divided in real space into nonoverlapping atomic fragments or basins such that the electron density gradient has a zero flow through the surface separating two neighbor basins: 3F(rS) 3 n(rS) = 0, where n(rS) is the normal vector to the surface at the rS point. The electron density, either experimental or theoretical, is the only source of the partition. Even more, the zero flux condition guarantees58 that the same quantum mechanical operators applied to the whole system are locally well-defined within the atomic basins, and that each basin contributes additively to produce the global property. Volumes,9 charges,10,11 magnetic susceptibilities,12 polarizabilities,13 and energies14,15 are examples of the molecular properties whose atomic partition has been discussed and exploited. In the case of crystals, bulk modulus and its inverse, the isotropic compressibility, are the best known examples of thermodynamic property partitioned using QTAIM.16 The uniform compressibility of a series of oxide spinels has been explained as a consequence of the dominant role of the oxide ion, whose local contribution can be transferred among similar crystals.17,18 More recently, the compressibility of ionic crystals has been analyzed in terms of the contribution of core, valence, and lone pair regions.19,20 Until now, however, there had not been a feasible method for partitioning the nonisotropic elastic components into atomic contributions. Such a method is presented and explored in this article for the first time. Given the importance of the elastic constants in determining the stability of a crystalline phase, the new analytical tool could help improve our understanding of the mechanisms of phase transitions. The quest for harder materials could also benefit from a deeper knowledge of the atomic and bonding mechanisms that regulate the shear component of the elastic deformation. The rest of the article is organized as follows: Section 2 describes the methods and conditions used to calculate the electronic structure r 2011 American Chemical Society
and energy of the crystals and the topological techniques used to analyze the electron densities. The partitioning of the bulk modulus and the compressibility, which will serve as a consistency test for the elastic constants, is presented in section 3. Section 4 contains the most important contribution from this work: starting with the main idea that makes the partition possible, we analyze the particular treatment of the diagonal and nondiagonal stress and shear elastic constants, perform a collection of tests to ensure the internal consistency of the procedure, and end analyzing the atomic contributions on a set of 12 crystals of four different oxide families. The article ends with section 5, presenting the most relevant conclusions and a discussion on future work and applications.
2. ELECTRONIC STRUCTURE CALCULATIONS The electronic structure calculations reported in this article have been done using either WIEN2k21,22 or ABINIT23,24 codes. WIEN2k uses the muffin-tin approach and a basis of augmented plane waves to perform all-electron full potential linearized augmented plane waves (FPLAPW) calculations. Our new text interface runwien25 was used to drive the calculations, particularly the determination of the global elastic constants using finite cell deformations. The ABINIT code uses relativistic pseudopotentials and a basis of plane waves to determine the properties of the valence electrons. The pseudopotentials used in this work have been generated using the FHI code,26 starting from the general pseudopotential database of ABINIT. The PBE96 generalized gradient exchange and correlation functional27 has been used. Care has been taken to achieve energy convergence of 1 105 hartree or better in all calculations. All the results discussed in the paper correspond to the ABINIT calculations, Special Issue: Richard F. W. Bader Festschrift Received: May 5, 2011 Revised: August 26, 2011 Published: September 27, 2011 12953
dx.doi.org/10.1021/jp2041718 | J. Phys. Chem. A 2011, 115, 12953–12961
The Journal of Physical Chemistry A
ARTICLE
once we have confirmed their equivalence to the WIEN2k geometries and elastic constants. As a test of the calculation conditions, we have determined the equilibrium properties of all the crystals that will be examined in this work. The comparison between the theoretical and experimental results can be examined in Table 1. It can be observed that the GGA functional tends to produce equilibrium geometries larger than the experiment, and smaller than bulk moduli. There are exceptions, however, such as SrO and Li2O. The MgSiO3 occurs experimentally on a deformed perovskite structure, but we have decided to examine the high symmetry form for comparison with the CaTiO3 and SrTiO3 . Finally, we have found no data for Pb2O, but we have included it for comparison to the other C3 crystals. The topological analysis has been done in all cases using our critic code,28 improved with the QTREE integration technique.29 This involves working analytically with the electron density of WIEN2k and with a grid of density values in the case of ABINIT. In the last case, additionally, the electron valence density produced by the electronic structure code must be complemented with the atomic core densities produced by the pseudopotential generator. The integration of the basin volumes has been done Table 1. Theoretical and Experimental Equilibrium Properties for the Crystals Examined in This Worka this work crystal structure MgO CaO
a
B1 B1
a
B
c11 c12 c44
experimental a
B
c11 c12 c44
with the QTREE algorithm29 to a precision of 1 103 bohr3 or better in all crystals and geometries. Figure 1 shows, in two (2D) and three (3D) dimensions, the atomic basins of the MgO crystal. It is by following the shape and volume shape of those basins with carefully designed cell deformations that we are going to be able to partition the elastic constants.
3. PARTITIONING THE BULK MODULUS AND THE COMPRESSIBILITY Let us start reviewing how the compressibility, k, and bulk modulus, B, can be partitioned into atomic or group contributions.16 We start from the thermodynamic definition k ¼ B1 ¼
B1
5.1159 104 226 42 63 5.1396 91 174 47 56
MgSiO3
E21
3.5252 226 344 167 176
CaTiO3
E21
3.9636 162 314 86 83 3.795 172
SrTiO3
E21
3.9844 164 317 88 96 3.9053 174 317 102 124
Li2O
C1
4.5219 81 219 13 62 4.545 88 217 24 68
Na2O
C1
5.5055 46 108 15 29
K2O Cu2O
C1 C3
6.1472 36 62 24 29 6.436 4.3569 114 130 105 63 4.2696 112 116 105 13
Ag2O
C3
4.9271 65 67 64 30 4.720
Pb2O
C3
5.4244 58 67 53 38
253
5.49
Cell lengths are given in Å, bulk modulus and elastic constants are given in GPa.
ð1Þ
where the derivatives can be defined under adiabatic or isothermal conditions, producing the corresponding kS or kT properties. For convenience, our theoretical calculations will invoke the static approximation, assuming that the temperature is null (so there are no thermal effects) and the vibrational zero point energy can be neglected. Under this static approach, the pressure can be obtained as p = ∂E/∂V, where E and V are the cell energy and volume, respectively. The QTAIM analysis allows us to strictly divide the cell volume into basin contributions: V = ∑ΩVΩ, where VΩ is the volume of the basin Ω. Using this volume partition in eq 1, k¼
4.2262 148 267 88 148 4.2112 160 298 96 155 4.8888 106 217 50 72 4.8105 116 236 56 82
SrO
1 ∂V V ∂p
∑Ω fΩ kΩ
and B1 ¼
∑Ω fΩ B1 Ω
ð2Þ
where fΩ = VΩ/V is the fractional volume occupied by Ω and kΩ ¼ B1 Ω ¼
1 ∂VΩ VΩ ∂p
ð3Þ
is the local compressibility associated with the Ω basin. Notice the formal equivalence between the definitions of k and kΩ. When applied to oxide spinels, for instance, this partition has been used to justify why many spinels have very similar bulk moduli, because the oxide basin dominates the contribution (fO ≈ 0.7 0.8) and its local compressibility is rather transferable among the spinels.17 Some other relevant works include the analysis of oxides,30 nitrides,3133 semiconductors,34 and
Figure 1. Electron density and topology of the MgO crystal. The left image is a contour plot of the electron density in the [001] plane. The electron density gradient is depicted in the middle image, together with the critical points, bond paths, and basin edges on the same plane. The right figure shows a 3D perspective of the Mg (small, green in the online version) and O (large and red) basins. 12954
dx.doi.org/10.1021/jp2041718 |J. Phys. Chem. A 2011, 115, 12953–12961
The Journal of Physical Chemistry A
ARTICLE
constants and compliances are independent: c11, c12, c33, s11, s12, and s33. In the infinitesimal deformation approach, cij and sij can be expressed as derivatives: ! ! ∂τi ∂εi , sij ¼ ð10Þ cij ¼ ∂εj 0 ∂τj 0 ε ,0
where it is essential, for the rest of the report, that there are different constant conditions in the derivatives: constant strain in the case of the cij but constant stress in the case of the sij. 4.1. Partitioning the cij Elastic Constants. A basin partition of the crystal energy, E = ∑ΩEΩ, would provide a direct partitioning of the elastic energy,
Figure 2. Components of the stress tensor.
ϕ¼
35
molecular solids. There have been some attempts to formalize the rules determining the contribution of a particular basin18,36 and the partition analysis has been extended to the basins determined by the ELF function37 and the electron density Laplacian.20
4. PARTITIONING THE ELASTIC CONSTANTS The elastic equation of state38 starts by assuming that, for small deformations around an equilibrium configuration, the crystal energy can be written as E E0 1 1 ¼ ∑ cijkl εij εkl ¼ ∑ sijkl τij τkl ð4Þ ϕ¼ 2 ij, kl 2 ij, kl V0 where V0 and E0 are the cell volume and energy at equilibrium, and the εij are symmetric deformation components " # ∂δxj 1 ∂δxi þ εij ¼ ð5Þ 2 ∂xj ∂xi where δx = x x0 and x ∈ {a,b,c} is any of the cell lengths. The τij are stress components (see Figure 2), defined by ∂ϕ τij ¼ ∂εij
ð6Þ
∑kl cijkl εkl
and εij ¼
∑kl sijkl τkl
ð7Þ
where cij,kl and sij,kl are the elastic constant and elastic compliance components, respectively. The above equations are symmetric with respect to the ij T kl, i T j, and k T l exchanges. This symmetry is exploited in the Voigt notation39 to reduce the bidimensional τ and ε to one index vectors: ½τ11 , τ22 , τ33 , τ23 , τ31 , τ12 w ½τ1 , τ2 , τ3 , τ4 , τ5 , τ6
ð8Þ
and ½ε11 , ε22 , ε33 , 2ε23 , 2ε31 , 2ε12 w ½ε1 , ε2 , ε3 , ε4 , ε5 , ε6
∑Ω
ϕΩ , ϕΩ ¼
EΩ EΩ 0 V0
ð11Þ
and then of the elastic constants: ! ∂2 ϕΩ Ω Ω cij ¼ cij , cij ¼ ∂εi ∂εj 0 Ω
∑
ð12Þ
ε ,0
We are still working on the partition of the cell energy, so this method will not be pursued here. Alternatively, in this article will explore an indirect route, based on the partition of the cell volume V = ∑ΩVΩ: ! ∂τi ∂τi ∂V ∂τi ∂V Ω ¼ ¼ ¼ cΩ ð13Þ cij ¼ ij ∂εj 0 ∂V ∂εj ∂V Ω ∂εj Ω
∑
ε ,0
where cΩ ij
" ¼
∂V Ω ∂εj
!
∂V ∂τi
∑
1 # ð14Þ ε0 , 0
The advantage of this indirect route is that recent developments have allowed us to determine the basin volumes with the required precision,29,40 typically 45 significant figures on each VΩ. To use eq 14, we need to calculate (E, V, VΩ...) for appropriate finite deformations of the unit cell:
In the limit of infinitesimal deformations, stress and strain are related by linear transformations: τij ¼
τ0
0
R f R ¼ Rð1 þ eÞ
ð15Þ
where R represents the cell shape, and e is the deformation or strain matrix. In the case of cubic systems,41,42 1 0 10 2ε6 2ε5 a 0 0 1 þ ε1 0 C B C B CB 1 þ ε2 2ε4 C R ¼B A @ 0 a 0 A@ 2ε6 2ε5 2ε4 1 þ ε3 0 0 a ð16Þ and the volume of the cell is given by V ¼ V0 f1 þ ε1 þ ε2 þ ε3 þ ε1 ε2 þ ε2 ε3 þ ε3 ε1 þ ε1 ε2 ε3
ð9Þ
Voigt notation transforms the four indices c and s tensors into bidimensional symmetric matrices. In this way, the 81 possible elements of c and s can be seen to reduce to a maximum of 21 independent terms, with further reductions determined by the additional symmetry of the crystal cell. In the case of a cubic crystal, like all the phases studied in this report, only three elastic
4ε24 4ε25 4ε26 4ε1 ε24 4ε2 ε25 4ε3 ε26 þ 16ε4 ε5 ε6 g
ð17Þ
Figure 3 shows the deformations of the cubic cell required for the calculation of the c11, c12, and c44 independent elastic constants. Both c11 and c12 require working with a tetragonal deformation of the crystal, and c44 requires working with a monoclinic phase. Even if the equilibrium reference corresponds to a cubic crystal in all cases, this point is better calculated again for each deformation 12955
dx.doi.org/10.1021/jp2041718 |J. Phys. Chem. A 2011, 115, 12953–12961
The Journal of Physical Chemistry A
ARTICLE
Figure 3. Deformations of the cubic cell required for the calculation of the c11 (left), c12 (middle), and c44 (right) elastic constants.
using the same symmetry as the rest of the points. Otherwise, small differences in the calculation conditions, related to the symmetry, could give rise to a discontinuity in the energy that would ruin the calculation of derivatives. 4.2. Partitioning c11. Figure 4 shows the behavior of the MgO crystal under the ε1 deformation. Eleven values of ε1, from 0.05 to +0.05, have been calculated. A low degree polynomial has then been fitted to the elastic energy, ϕ(ε1), and to the cell and basin volumes, V(ε1) and VΩ(ε1). The derivatives required for the calculation of c11 have been obtained from the polynomial, and we have checked that changes in the polynomial degree influence the final results, quoted in Table 2, less than the significant figures reported. We can observe in Figure 4 that a null strain corresponds, as it should be, to a null stress. The stress, τ1(ε1), is not linear because we are doing finite deformations of the crystal. It is only the slope of the τ1(ε1) curve at the ε1 = 0 point that we need for our calculation. The cell volume, on the other hand, remains truly linear, a direct consequence of eq 17: V(ε1,ε0 = 0) = V0(1 + ε1). The basin volumes appear also linear in the deformation. The oxide basin is the main contributor both to the cell volume and to its slope. Table 2 shows that MgO is a rather ionic crystal, where the Mg basin donates 1.71 electrons to the O basin. Accordingly, the oxide basin occupies most of the cell volume (96 of 127 bohr3) and has the leading contribution to c11 (216 of 267 GPa). Whereas the Mg basin is more difficult to compress (see the small ∂VMg/∂ε1 derivative), the larger volume of the oxide basin more than compensates for this. 4.3. Partitioning c12. The calculation of c12 follows the same route as the calculation of c11 described above. The only difference is that now we have to work with energy and volume surfaces: ϕ(ε1,ε2), V(ε1,ε2), and VΩ(ε1,ε2). A grid of 11 11 points was calculated, with a [0.06,+0.06] range for both strains. The null strain line was carefully avoided to prevent a change in the crystal symmetry for different points in the grid. Several numerical methods were tried to determine the derivatives required, with the final results obtained by means of a 2D polynomial fitted to the grid points. The stability of the results with respect to the polynomial degree and the symmetry between c12 and c21 have been carefully checked. The calculation also gives c11 as a subproduct, and we have also checked that the results are equivalent to the values examined in the previous section. Figure 5 shows the ϕ(ε1,ε2) elastic energy of MgO. The energy shape closely resembles a paraboloid, but the cross terms, involving simultaneously ε1 and ε2, are important components of the fitting 2D polynomial. In this case the volume is not required to be linear in the deformations: V ðε1 , ε2 , ε0 ¼ 0Þ ¼ V0 ð1 þ ε1 þ ε2 þ ε1 ε2 Þ
ð18Þ
Figure 4. Elastic energy ϕ(ε1) (top), τ1 stress (center), and cell and basin volumes (bottom) of the MgO for the ε1 deformation.
Table 2. Integrated Charge, Volume, and c11 Atomic Contributions to the MgO Crystal Properties Mg
O 1.71
Cell
Q Ω (e)
1.71
V Ω (bohr3)
31.70
95.65
127.35
cΩ 11 (GPa)
51.50
215.56
267.06
In practice, however, in addition to the cell volume, the basin volumes also remain essentially planar for the examined deformations. The oxide basin is again the main contributor to this elastic O constant: cMg 12 = 16.87, c12 = 71.13, and c12 = 88.09 GPa. 4.4. Partitioning s11 and s12. The partition of the compliances requires doing calculations in which most of the stress 12956
dx.doi.org/10.1021/jp2041718 |J. Phys. Chem. A 2011, 115, 12953–12961
The Journal of Physical Chemistry A
ARTICLE
Table 3. Compliances s11 and s12 Calculated for the MgO Crystal Mg s11Ω s12Ω
1
O
Cell
(TPa )
0.8628
3.6164
4.4792
(TPa1)
0.2132
0.8970
1.1102
Figure 5. Elastic energy of the MgO for the (ε1,ε2) deformation.
Figure 7. Cell and basin volumes of MgO as a function of the ε4 deformation.
Figure 6. Elastic energy of the MgO as a function of the τ1 stress.
components remain fixed to zero. Namely, the definition ! ∂εi ð19Þ sij ¼ ∂τj 0 τ ,0
leads to the volume-defined partition " ! # Ω ∂V ∂V 1 Ω sij ¼ sΩ ij , sij ¼ ∂εi ∂τj Ω
∑
ð20Þ τ0 , 0
We have tested two different approaches. In the first method, we have exploited the capability of the ABINIT code to introduce constant stress conditions into the electronic structure calculation.4345 The second method assumes the linear stress strain relationship, ε = sτ = c1τ, to determine the cell shape that would correspond to a given stress condition. It must be noticed that the first approach does not assume linearity, which is only fulfilled for infinitesimal deformations, as the results in Figure 4 (center) attest. However, the results obtained from both approaches are equivalent. Figure 6 depicts the behavior of the elastic energy with the τ1 stress. The comparison with Figure 4 shows that ϕ(τ1) is close to a scaled version of ϕ(ε1). In a similar way, the ϕ(τ1,τ2) appears to be similar to the ϕ(ε1,ε2) surface in Figure 5. Table 3, on the other hand, reveals that the oxide basin provides the dominant contribution to both s11 and s12, with the Mg contribution being close to 20% of both compliances. It is interesting to notice that s12 is negative, but this does not imply a violation of the Born stability conditions: c11 þ 2c12 > 0, c11 c12 > 0, c44 > 0
ð21Þ
known as the spinodal, tetrahedral, and rhombohedral shear conditions, respectively. 4.5. Partitioning c44 and s44. The calculation of c44 and s44 cannot follow the same procedure used for the other constants. This is caused by the fact that the derivatives required by eqs 14 and 20 are zero, thus giving rise to a 0/0 indeterminacy. Indeed, eq 17 shows that the cell volume depends quadratically on ε4 and then (∂V/∂ε4)ε0 ,0 = 0. The same behavior is inherited by the basin volumes, as Figure 7 clearly reveals. In a similar way, the derivatives of the cell and basin volumes with respect to τ4 are also zero, as the V(τ4) and VΩ(τ4) also behave as a parabola. Fortunately, we can apply the L’H^opital rule to solve the 0/0 indeterminacy. In this way, the basin contribution to an elastic constant can be written as ∂2 V Ω ∂2 V ð22Þ c44 ¼ 2 = Ω ∂ε4 ∂τ 4 ∂ε4 0
∑
ε ,0
and, similarly, s44 ¼
∑Ω
∂2 V Ω ∂2 V = ∂τ4 ∂ε4 ∂ε24
ð23Þ τ0 ¼ 0
The curvatures of V and VΩ are not null, so the above formulas can be used to determine the basin contributions reported in Table 4. We can see, again, that the contribution from the oxide basin dominates both c44 and s44. The most interesting result is, however, that the Mg basin has a negative contribution to both constants. In effect, Figure 7 shows clearly that VMg has a contrary curvature to both the cell and oxide volumes. The calculation of the curvatures demands a very precise determination of the basin volumes that has only been possible because of the new and improved QTREE technique developed for the critic2 code.29,40 On the other hand, the monoclinic stress and strain are proportional, τ4 = c44s4, and the same energy, ϕ(ε4,ε0 = 0), and volume data can be used to obtain and partition both c44 and s44. 12957
dx.doi.org/10.1021/jp2041718 |J. Phys. Chem. A 2011, 115, 12953–12961
The Journal of Physical Chemistry A
ARTICLE
4.6. Internal Consistency of the Partitioning. The topological partition that we have examined is complete and exact, in the sense that the cell cij and sij are exactly reproduced as a sum of the basin contributions. On the other hand, there are some wellknown properties fulfilled by the elastic constants and compliances of the whole cell. First, both matrices are mutually inverse,
c1 ¼ s
However, eqs 24, 25, and 27 are not fulfilled by the basin contributions:
ð24Þ
Second, in the particular case of cubic crystals, B ¼ ðc11 þ 2c12 Þ=3
ð25Þ
k ¼ 3ðs11 þ 2s12 Þ
ð26Þ
c44 ¼ 1=s44
ð27Þ
ð28Þ
Table 4. Basin Contributions to c44 and s44 in the MgO Crystal Mg
O
Cell
c44 (GPa)
25.19
172.79
147.79
s44 (TPa1)
1.149
7.915
6.767
ð29Þ
Ω cΩ 44 6¼ 1=s44
ð30Þ
ðcΩ Þ1 6¼ sΩ
ð31Þ
We can also notice that the topological partition based on the energy (eq 12) would have different properties than the partition based on the volume (eq 14) that we have explored in this work. Namely, in the case of the energy-based partition, the cΩ 11 and Ω Ω Ω cΩ 12 would be additive to produce B , whereas s11 and s12 would not add to give kΩ. Henceforth, the topological partition is not unique, but there are two equally well-defined forms that show different behaviors and properties. Ω 4.7. Transferability of the cΩ ij and sij Contributions. We have examined in the previous sections the partitioning of the cij and sij constants in the MgO crystal. In this section we analyze the results from 12 different crystals that belong to four different oxide families: perovskites (CaTiO3, SrTiO3, and MgSiO3), rock-salt crystals (MgO, CaO, and SrO), inverse fluorites (Li2O, Na2O, and K2O), and cuprite-like crystals (Cu2O, Ag2O, and Pb2O). This selection of crystals is designed to test the behavior of the oxide basin in different environments. The coordination index of the oxide, in particular, can be 2 (perovskites), 4 (inverse fluorite and cuprite), and even 6 (rock-salt). Whereas all the crystals analyzed are cubic by economy, the programs and methods described are fully general and ready to use on any crystal.
It is most reasonable to check whether these properties also hold for the basin counterparts. The analysis of the partition equations (14 and 20) and of the numerical values of the different magnitudes shows that the relationship between k and the sij (eq 26) can also be used on each basin: Ω f Ω kΩ ¼ 3ðsΩ 11 þ 2s12 Þ
Ω BΩ 6¼ ðcΩ 11 þ 2c12 Þ=3
Table 5. Components of the Elastic Properties for the Cubic Perovskite Crystalsa Ω c (%) f (%) VΩ (bohr3)
Ca
Ti
O
CaTiO3
Sr
Ti
O
SrTiO3
Mg
Si
O
MgSiO3
— —
— —
— —
65 4.5
— —
— —
— —
65 4.9
— —
— —
— —
85 11.7
89
50
94
420.2
Ill
51
88
426.9
38
18
80
296.6
Q Ω (e)
1.63
2.17
1.27
0.0
1.60
2.17
1.26
0.0
1.77
3.27
1.68
0.0
fΩ
0.21
0.12
0.67
1.0
0.26
0.12
0.62
1.0
0.13
0.06
0.81
1.0
BΩ (GPa)
209
206
147
1962
212
208
146
164
284
196
221
226
kΩ (TPa1)
4.79
4.85
6.80
3.61
4.73
4.80
6.86
3.58
3.52
5.09
4.52
4.25
52
29
317
65
30
321
35
24
cΩ 11 (GPa) cΩ 12 (GPa) cΩ 44 (GPa)
14 23
8 29
1 sΩ 11 (TPa )
0.60
0.34
1 sΩ 12 (TPa )
0.13
0.07
1 sΩ 44 (TPa )
3.25
4.09
55
(2)
126
(1)
15
(2)
34
(1)
86
30
26
(2) (1)
0.03 2.01
(2) (1)
3.61
0.72
0.34
0.07
(2)
0.78
0.16
0.07
0.43
(1)
3.74
(2) (1)
11.98
25
8
37
5.33
84
18
2.76
58
(2)
111
(1)
16
(2)
30
(1)
88
65
42
(2) (1)
0.43 1.67
(2) (1)
3.58
0.43
0.31
0.09
(2)
0.78
0.14
0.10
0.36
(1)
4.57
(2)
1.88
(1)
10.48
37
12
17
3.30
95
17
1.19
2.13
64
(2)
158
(1)
32
(2)
76
(1)
345 167
98
(2)
9
(1)
0.27 4.06
(2) (1)
4.25
+0.06
(2)
1.39
1.28
(1)
3.14
(2)
5.67
(1)
176
5.67
The c and f values, given as percent ratios, correspond to the covalent and flatness topological indices46 (see text). Other properties are: B and BΩ (bulk modulus); VΩ (basin volume); Q Ω (basin charge); fΩ (basin fraction of the cell volume; see eq 2); kΩ and sΩ ij (basin contribution to the compressibility and compliances); and cijΩ (basin contribution to the elastic constants). a
12958
dx.doi.org/10.1021/jp2041718 |J. Phys. Chem. A 2011, 115, 12953–12961
The Journal of Physical Chemistry A
ARTICLE
Table 6. Components of the Elastic Properties for the Rock Salt Crystalsa Ω
O
MgO
Ca
O
CaO
Sr
O
SrO
Ω
Li
O
Li2O
O
Na2O
K
O
K2O
—
—
86
—
—
74
—
—
73
c (%)
—
f(%) VΩ
— 32
— 96
25 127.3
— 83
— 114
15 197.1
— 114
— 112
13 225.9
f (%) VΩ
— 20
QΩ
1.71 1.71 0.00
1.48 1.48 0.00
1.47 1.47 0.00
QΩ
0.84 1.68
0.25
0.75
1.00
0.42
0.58
1.00
0.51
0.49
1.00
fΩ
0.26
192
137
148
150
90
105
138
82
104
BΩ
107
kΩ
5.22
7.28
6.76
6.65
11.08
9.23
7.26
12.18
9.69
kΩ
9.32 13.32 12.30 16.78 25.85 016.78 24.32 37.55 28.83
cΩ 11 cΩ 12
52 17
216 71
267 88
67 16
155 36
223 52
85 16
139 26
224 42
cΩ 11 cΩ 12
43 2.39
176 218 9.91 12.30
37 5.17
70 9.75
cΩ 44
26
173
147
7
67
74
28
35
63
cΩ 44
29.8
32.4
62.2
1.35
27.2
28.6
sΩ 11
0.86
3.62
4.48
1.52
3.51
5.03
1.78
2.92
4.70
sΩ 11
0.88
3.72
4.60
3.34
6.27
9.61 11.40
B
Ω
sΩ 12 sΩ 44
1.15 7.92
6.77
1.33
12.23 13.56
7.04
8.87
a
We can start by analyzing the ionicity and flatness topological indices46 of all the crystals. The ionicity index is defined as an average for all the basins of the ratio between the actual topological charge and the nominal oxidation state, Ω
1 Q N Ω ¼ 1 OSΩ N
∑
sΩ 44
15.91
Notation and units follow the conventions of Table 5.
c¼
sΩ 12
0.21 0.90 1.11 0.28 0.65 0.93 0.28 0.46 0.74
—
—
79
—
— 63
— 155
9 281.5
— 129
0.00
0.79 1.59
0.00
0.74 1.48
0.00
0.74
1.00
0.45
0.55
1.00
0.66
0.34
1.00
75
83
60
39
46
41
27
36
107 33 14.92 12.50 0.99
—
74
— 6 134 391.9
26 59 9.95 22.45 26.9
27.9
9.09 20.49
0.05 0.21 0.25 0.42 0.79 1.20 3.14 2.51 5.65 7.67
8.41 16.08
1.66 33.31
34.97
1.32 34.54 35.87
Notation and units follow the conventions of Table 5.
Cu
ð32Þ
ð33Þ
where the numerator is the absolute minimum of density in the cell, necessarily a cage critical point, and the denominator is the density at the bond critical point of maximum electron density. The f ratio provides a measurement of how wide the valence electron density actually is, being close to 100% in metallic systems and close to 0% in highly ionic or highly covalent compounds. Together, the c and f indices have been shown to provide a classification of compounds into the three ionicmetallic-covalent classical prototypes.46 All the crystals reported here show a very low value of f, the highest value being the 25% found on MgO, also a record in charge transfer: c = 86%. A few things appear to be clear from the results reported in Tables 48. The oxide basin is always a leading contributor to the cij and sij. Furthermore, considering each crystal family independently, we can see that cO ij tends to increase significantly when its nearest neighbor is a small cation and decrease when it is a large cation. In rock-salt crystals, for instance: cO 11 is 216 GPa
84
— 13 116 156.0
Table 8. Components of the Elastic Properties for the Cupritelike Crystalsa
thus providing a global measurement of the charge transfer for the whole crystal. The value of c shows a clear difference between the crystal families: rock-salt and inverse-fluorites range in the 7386% typical of very ionic compounds, perovskites appear in the middle, with a 65% transfer for the titanates and 85% for the silicate, and cuprites show the smallest charge transfer, in the range 4658%, characteristic of compounds with a significant degree of covalency. Flatness is defined46 as the ratio between the electron density of the crystal at two significant points: Fmin c f ¼ max Fb
—
Na
c(%)
fΩ
a
Mg
Table 7. Components of the Elastic Properties for the Inverse Fluorite Crystalsa
Cu2O
Ag
O
Ag2O
Pb
O
Pb2O
c (%)
—
—
53
—
—
46
—
—
58
f (%)
—
—
5
—
—
4
—
—
5
VΩ QΩ
102 74 279.1 0.53 1.06 0.00
159 85 403.6 0.46 0.93 0.00
226 102 538.6 0.57 1.18 0.00
fΩ
0.73
0.27
1.00
0.79 0.271
1.00
0.82
0.18
1.00
BΩ
113
105
114
65
49
44
58
kΩ
8.88
9.51
8.79
cΩ 11
88
34
122
50
16
66
47
12
cΩ 12
73
29
101
48
15
63
35
9
44
cΩ 44
51
7.8
58
22
1.3
30
25
8.1
37
67
56
15.00 17.80 15.38
20.39 22.62 17.37 59
sΩ 11
0 18.73 9.19
sΩ 12
8.97 3.53 12.50 78.2 25.9 104.1 16.84 4.22 21.06
sΩ 44
14.75 2.36 0.02
a
O
10 5 27.92 140.01 72.85 212.86 37.70 10.21 47.91 17.14
22.9 0.92 32.8 10.85
17.78
5.86
17.78
3.1
Notation and units follow the conventions of Table 5.
(MgO), 155 GPa (SrO), and 139 GPa (SrO), and the same trend O is observed on cO 12 (71, 36, and 26 GPa, respectively) and c44 (173, 67, and 35 GPa, respectively). There are only a few exceptions to this rule, and they are always related to very small elastic constants in crystals that appear to be close to an instability and where the numerical accuracy of our partition could be close to the actual constant: the negative cO 44(Ag2O) = 1.3 GPa is probably the most notable exception, but then c44(Ag2O) = 30 GPa. Unfortunately, it is not clear how the trend of cO ij with the nearest neighbor size can be generalized to include all the crystal families together. We should expect some dependence on the coordination index of the oxide, but it is not apparent. We could also argue about the correlation with the topological charge of the oxide but, again, this is true on each crystal family but it does not hold when we put all the crystal families together. 12959
dx.doi.org/10.1021/jp2041718 |J. Phys. Chem. A 2011, 115, 12953–12961
The Journal of Physical Chemistry A
ARTICLE
Figure 8. Basin plots for the CaTiO3 (left, perovskite prototype), Li2O (middle, inverse fluorite prototype), and Cu2O (right, cuprite).
It is interesting to observe how different basins of the same species can respond very differently to the applied stress or strain when the symmetry allows it. This happens for the oxide basin in the AMO3 perovskites and for the metallic basin in the M2O cuprite families. The O basin in CaTiO3, for instance, has the shape of a cushion (see Figure 8), easier to deform in the direction of the two bond paths with the tetravalent cation Ti than in the direction of the four bond paths with the divalent cation Ca. Correspondingly, there are two small and one large cO ij contribution in this crystal. When we consider the cation contributions and consider the different families apart, there are also recognizable rules. First, the trends are more clear on the larger than on the smaller cij. In particular, c11 appears more regular in its trends than the much smaller c12. This could be simply related to numerical precision issues. Second, larger cations tend to have a larger contribution. Third, the same cation tends to have similar contributions when occurring on different crystals. The only case that we can test cleanly is Ti having almost the same cTi ij in CaTiO3 and SrTiO3, but it is also compelling to note the similar, although not identical, contribution of Mg on MgO and MgSiO3, of Ca on CaO and CaTiO3, and of Sr on SrO and SrTiO3 . This is, perhaps, a route worth pursuing; examining a few more cases, perovskites and spinels come to mind, where crystals of the same family share cations in different ways.
5. CONCLUSIONS AND OUTCOME • The QTAIM partitioning of the cij and sij elastic constants is possible by following the evolution of the cell energy and volume, E and V, and that of the basin volumes, VΩ, under appropriate finite deformations of the unit cell. • Keeping the stresses constant to calculate the sij compliances is possible either by using the linear strainstress relationship (ε = sτ) or by enforcing constant stress conditions in the electronic structure calculations, like the ABINIT code does. • A good precision on the basin volumes is required, particularly in the computation of the delicate shear moduli. This has been possible with the new techniques developed for critic2, particularly the QTREE integration method.29 • The sΩ ij values obtained by partitioning the cell volume into basin contributions fulfill the internal consistency rules. For Ω instance: kΩ = 3(sΩ 11 + 2s12). • Partitioning the cell energy instead of the volume would give a different scheme with different properties. In that case, the Ω cΩ ij ’s would be consistent with BΩ, but the sij ’s would not be consistent with kΩ. • We are still working on the rules that govern the basin contribution of an atom. Where we have been able to recognize some regular trends, there are still many open questions,
particularly related to the changes in the oxide contribution across the different crystal families.
’ AUTHOR INFORMATION Corresponding Author
*E-mail:
[email protected].
’ ACKNOWLEDGMENT We would first like to acknowledge Prof. Richard F. W. Bader, whose relentless pursuit of rigor and beauty in the description and analysis of chemical bonding is a permanent guide and inspiration. This research is financed by the Spanish Ministerio de Ciencia e Innovacion (MICINN), grant CTQ2009-08376, and the Malta/Consolider initiative CSD2007-00045. A.O.R. is indebted to the Spanish Ministerio de Educacion y Ciencia (MEC) for an FPU grant. ’ REFERENCES (1) de Boer, F. R.; Boom, R.; Mattens, W. C. M.; Miedema, A. R.; Niessen, A. K. Cohesion in Metals: Transition Metal Alloys; Elsevier Science Publishers: Amsterdam, 1988; Vol. 1, p 758. (2) Wildman, S. A.; Crippen, G. M. J. Chem. Inf. Comput. Sci. 1999, 39, 868–873. (3) Bader, R. F. W. Atoms in Molecules: A Quantum Theory; Oxford University Press: Oxford, 1990. (4) The Quantum Theory of Atoms in Molecules: From Solid State to DNA and Drug Design; Boyd, R., Matta, C., Eds.; Wiley-VCH: Weinheim, Germany, 2007; pp 527 and xxxviii. (5) Biegler-K€onig, F. W.; Bader, R. F. W.; Duke, A. J.; Nguyen-Dang, T. T.; Tal, Y. J. Phys. B 1981, 14, 2739–2751. (6) Tal, Y.; Anderson, S. G.; Bader, R. F. W.; Nguyen-Dang, T. T.; Ojha, M. J. Chem. Phys. 1981, 74, 5162–5167. (7) Biegler-K€onig, F. W.; Bader, R. F. W.; Tang, T. H. J. Comput. Chem. 1982, 3, 317–328. (8) Bader, R. F. W.; Popelier, P. L. A. Int. J. Quantum Chem. 1993, 45, 189–207. (9) Wiberg, K. B.; Bader, R. F. W.; Lau, C. D. H. J. Am. Chem. Soc. 1987, 109, 1001–1012. (10) Bader, R. F. W.; Zou, P. F. Chem. Phys. Lett. 1992, 191, 54–58. (11) Bader, R. F. W.; Matta, C. F. J. Phys. Chem. A 2004, 109, 8385–8394. (12) Bader, R. F. W.; Keith, T. A. J. Chem. Phys. 1993, 99, 3683–3693. (13) Bader, R. F. W.; Keith, T. A.; Gough, K. M.; Laidig, K. E. Mol. Phys. 1992, 75, 1167–1189. (14) Bader, R. F. W.; Matta, C. F. Int. J. Quantum Chem. 2001, 85, 592–607. (15) Martín Pendas, A.; Blanco, M. A.; Francisco, E. J. Comput. Chem. 2007, 28, 161–184Special issue on Chemical Bonding. (16) Pendas, A. M.; Costales, A.; Blanco, M. A.; Recio, J. M.; Lua~ na, V. Phys. Rev. B 2000, 62, 13970–13978. 12960
dx.doi.org/10.1021/jp2041718 |J. Phys. Chem. A 2011, 115, 12953–12961
The Journal of Physical Chemistry A
ARTICLE
(17) Recio, J. M.; Franco, R.; Pendas, A. M.; Blanco, M. A.; Pueyo, L.; Pandey, R. Phys. Rev. B 2001, 63, 184101-1–7. (18) Taravillo, M.; del Corro, E.; Contreras-García, J.; MartínPendas, A.; Florez, M.; Recio, J. M.; Baonza, V. G. High Pressure Research 2009, 29, 97–102. (19) Contreras-Garca, J.; Pendas, A. M.; Silvi, B.; Recio, J. M. J. Phys. Chem. B 2009, 113, 1068–1073. (20) Otero-de-la Roza, A.; Lua~na, V. J. Chem. Theory Comput. 2010, 6, 3761–3779(Publication Date (Web): November 10, 2010). (21) Schwarz, K.; Blaha, P.; Madsen, G. K. H. Comput. Phys. Commun. 2002, 147, 71–76. (22) Schwarz, K.; Blaha, P. Comput. Mater. Sci. 2003, 28, 259–273. (23) Gonze, X.; et al. Comput. Mater. Sci. 2002, 25, 478–492. (24) Gonze, X.; et al. Z. Kristallogr. 2005, 220, 558–562. (25) Otero-de-la Roza, A.; Lua~na, V. Comput. Phys. Commun. 2009, 180, 800–812(Source code distributed by the CPC program library: http://cpc.cs.qub.ac.uk/summaries/AECM_v1_0.html). (26) Fuchs, M.; Scheffler, M. Comput. Phys. Commun. 1999, 119, 67. (27) Perdew, J.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865–3868. (28) Otero-de-la Roza, A.; Blanco, M. A.; Martín Pendas, A.; Lua~ na, V. Comput. Phys. Commun. 2009, 180, 157–166(Source code distributed by the CPC program library: http://cpc.cs.qub.ac.uk/summaries/AECB_v1_0. html). (29) Otero-de-la Roza, A.; Lua~na, V. J. Comput. Chem. 2011, 32, 291–305. (30) Gracia, L.; Beltran, A.; Andres, J.; Franco, R.; Recio, J. M. Phys. Rev. B 2002, 66, 224114–17. (31) Mori-Sanchez, P.; Franco, R.; Pendas, A. M.; Lua~na, V.; Recio, J. M. Europhys. Lett. 2001, 54, 760–766. (32) Mori-Sanchez, P.; Marques, M.; Beltran, A.; Jiang, J. Z.; Gerward, L.; Recio, J. M. Phys. Rev. B 2003, 68, 064115–15. (33) Marques, M.; Osorio, J.; Ahuja, R.; Florez, M.; Recio, J. M. Phys. Rev. B 2004, 70, 104114–111. (34) Ouahrani, T.; Otero-de-la Roza, A.; Reshak, A. H.; Khenata, R.; Faraoun, H. I.; Amrani, B.; Mebrouki, M.; Lua~na, V. Phys. B: Condens. Matter (Amsterdam, Neth.) 2010, 405, 3658–3664. (35) Gracia, L.; Marques, M.; Beltran, A.; Martín Pendas, A.; Recio, J. M. J. Phys.: Condens. Matter 2004, 16, S1263–S1270. (36) Contreras-Garca, J.; Mori-Sanchez, P.; Silvi, B.; Recio, J. M. J. Chem. Theory Comput. 2009, 5, 2108–2114. (37) Contreras-Garca, J.; Pendas, A. M.; Silvi, B.; Recio, J. M. J. Phys. Chem. Solids 2008, 69, 2204–2207. (38) Grimvall, G. Thermophysical Properties of Materials; North Holland: Amsterdam, 1999; p 444. (39) Nye, J. F. Physical Properties of Crystals; Oxford University Press: Oxford, U.K., 1985; p 329. (40) Otero-de-la Roza, A.; Lua~na, V. J. Comput. Chem. 2011, 32, 291–305. (41) Mehl, M. J.; Osburn, J. E.; Papaconstantopoulos, D. A.; B., K. Phys. Rev. B 1990, 41, 10311–10323. (42) Mehl, M. J. Phys. Rev. B 1993, 47, 2493–2500. (43) Hamann, D. R.; Wu, X.; Rabe, K. M.; Vanderbilt, D. Phys. Rev. B 2005, 71, 035117. (44) Hamann, D. R.; Wu, X.; Rabe, K. M.; Vanderbilt, D. Phys. Rev. B 2005, 72, 079901. (45) Wu, X.; Vanderbilt, D.; Hamann, D. R. Phys. Rev. B 2005, 72, 035105. (46) Mori-Sanchez, P.; Martín Pendas, A.; Lua~na, V. J. Am. Chem. Soc. 2002, 124, 14721–14723.
12961
dx.doi.org/10.1021/jp2041718 |J. Phys. Chem. A 2011, 115, 12953–12961
