Determination of elastic constants of titanium diboride (TiB2) from first principles using FLAPW implementation of the density functional theory

Computational Materials Science 35 (2006) 134–150 www.elsevier.com/locate/commatsci

Determination of elastic constants of titanium diboride (TiB2) from first principles using FLAPW implementation of the density functional theory K.B. Panda, K.S. Ravi Chandran

*

Department of Metallurgical Engineering, 135 South 1460 East, Room 412, University of Utah, Salt Lake City, UT 84112, United States Received 18 November 2004; received in revised form 18 January 2005; accepted 16 March 2005

Abstract First principles computational calculations of anisotropic elastic constants of titanium diboride, TiB2, were performed using the implementations of the Hartree–Fock (HF) method and the density functional theory (DFT). TiB2 has hexagonal crystal structure, thus five independent elastic constants are needed to completely determine its elastic properties including polycrystalline elastic modulus, PoissonÕs ratio and the elastic anisotropy of the crystal. The HF method employed molecular orbitals constructed from the linear combination of atomic orbitals (LCAO). The DFT calculations were based on the full potential linearized augmented plane wave (FLAPW) method with the generalized gradient approximation (GGA). Five independent elastic distortions of the unit cell were employed to determine the anisotropic elastic constants under the unrelaxed and relaxed configurations of Ti and B atoms in the unit cell. The calculation methods as well as the internal atomic relaxations of the elastic cell distortions were found to have a significant effect on the numerical values of elastic constants. Estimations of polycrystalline elastic constants and their comparison with the experimentally determined values were also performed. The agreement of the DFT (FLAPW) calculations including internal atomic relaxations, with the experimental data is very good. The HF calculations overestimated the elastic constants upto around 20%. Elastic anisotropy, the nature of chemical bonding and the electronic charge transfer between constituent atoms in TiB2 have also been explored to assess the origins of high elastic stiffness of this compound. Ó 2005 Elsevier B.V. All rights reserved.

1. Introduction Titanium and boron react to form an interesting set of compounds, TiB, Ti3B4 and TiB2 at B concentrations (wt.%) of about 18%, 22% and 30% and they have important scientific and technical significance. TiB and Ti3B4, are orthorhombic [1,2], while TiB2 is hexagonal [3,4]. TiB2 is one of the well known metallic-like ceram-

*

Corresponding author. Tel.: +1 801 581 7197; fax: +1 801 581 4937. E-mail address: [email protected] (K.S. Ravi Chandran).

0927-0256/$ - see front matter Ó 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2005.03.012

ics that is electrically conducting and it has high elastic modulus (about 565 GPa) [4], high hardness (about 24 GPa, [4] VickerÕs), high specific strength and wear resistance that are indispensable in specialized applications. Some applications of this material include armor, cutting tools, crucibles, wear-resistant coatings and electrodes in aluminum-extraction cells. The present study illustrates the computational procedure employed to obtain the five independent elastic constants of the hexagonal TiB2 from first principles calculations on the basis of computational implementation of the Hartree–Fock (HF) and the density functional theory (DFT) of electronic structure. One of the objectives is to evaluate the effectiveness of these approaches

K.B. Panda, K.S. Ravi Chandran / Computational Materials Science 35 (2006) 134–150

in providing accurate elastic constants of TiB2. The study, based on ground state electron densities, also explores the nature of bonding between B–B atoms and between the Ti–B and Ti–Ti atoms. The density of states in the valence band and the three-dimensional character of electron charge density distributions along different crystallographic directions have been analyzed to understand the origin of high hardness and elastic modulus in TiB2. 1.1. First principle based computational determination of properties Computational calculation of ground state properties of assemblies of atoms using first principles methods is an effective approach in the determination of structural, magnetic, optical, dielectric and superconducting properties of compounds, as this involves no a priory assumptions about the electronic structure and atomic interactions. Studies of electronic structure of transition metal compounds such as TiC [5] TiCxN1
x [6] and the pressure-dependent nonlinear elasticity and the associated phase transitions in earth-mantle compounds such as MgSiO3 [7] and MgO [8], have recently shown the success of this approach. In particular, the notion that elastic properties are calculable from the ground state crystal energies determined by the equilibrium electron densities of the cluster, molecular unit or the unit cell, is quite appealing, especially when experimental work is difficult. Further, the ground state electron density information can be used to gain insight into the directionality in chemical bonding between constituent atoms and its effect on elastic anisotropy and the bulk properties. In this study, Crystal98 [9], the code that implements HF method, and the DFT implementation, WIEN2K [10], a code that uses full potential linearized augmented plane wave method (FLAPW) to self-consistently determine the ground state crystal energies, were used to determine the electronic structure, anisotropic elastic constants and nature of chemical bonding in TiB2. The FLAPW method is an all electron method and is (within DFT framework) widely applicable to all atoms of the periodic table. It is widely considered to be the most precise electronic structure calculations method in solid state physics. The advantage of the FLAPW method compared with the traditional LMTO method is that within the muffin-tin spheres centered on atomic sites, no shape approximation is made to the potential or charge density, and in the interstitial region, no volume average is taken. The relative effectiveness of the two approaches used in this study (HF versus DFT) in providing accurate elastic constants is assessed. The research not only provides the most accurate values of anisotropic elastic constants of TiB2, but also provides insight into the chemical origins of high elastic modulus and elastic anisotropy in this compound.

135

2. Methodology 2.1. Crystal structure of TiB2 TiB2 crystallizes in C32 structure, which is hexagonal (P6/mmm) [11] belonging to the space group 191 with mmm symmetry. The lattice parameters obtained from the X-ray powder diffraction work are a0 = 3.030 ± ˚ , b0 = 3.03 ± 0.01 A ˚ and c0 = 3.232 ± 0.01 A ˚, 0.01 A respectively [12]. Fig. 1(a) illustrates the crystal structure of TiB2. The fundamental building block of TiB2 is the trigonal prism with one B atom at the center and Ti atoms in corners. The hexagonal unit cell is formed by the close-packed vertical stacking of the trigonal prisms, with each prism sharing all its faces with neighboring prisms. The primitive cell containing the 1 Ti and 2 B atoms is illustrated in Fig. 1(b). There is one B and one Ti in nonequivalent atomic positions of Ti (0, 0, 0) and B (1/3, 2/3, 1/2). 2.2. Calculation of elastic constants Crystals deform almost in a linear elastic manner at small strains and the strain energy density increment of a homogeneous and elastically deformed crystal is given by dE ¼ rij ddij

ð1Þ

where rij are the elements of stress tensor and dij are the elements of strain tensor. In a more generalized form the elastic constants can be expressed as cijkl ¼

o2 E odij odkl

ð2Þ

According to Wallace [13], the internal energy of a crystal under a general strain dij can be expressed by expanding the internal energy E(V) of the deformed crystal with respect to the strain tensor, in terms of the TaylorÕs series, as X V0 X EðV Þ ¼ EðV 0 ; 0Þ þ V 0 rij dij þ cijkl dij dkl þ    2! i;j;k;l i;j ð3Þ where V0 is the volume of the unstrained crystal and E(V0, 0) is the corresponding ground state energy. In the above expression the strain tensors subscripts (ij, kl) are expressed in Voigt notation scheme (11 = 1, 22 = 2, 33 = 3, 23 = 4, 31 = 5 and 12 = 6). It may be noted that Eq. (3) is of the form EðV Þ ¼ M 0 þ M 1 di þ M 2 d2i þ   

ð4Þ

with M0 = E(V0, 0), M1 = V0ri and M2 = (V0cii/2.) Therefore, the second-order coefficient of a polynomial fit to strain energy density with respect to the deforming strain component should give the elastic constant for that strain

136

K.B. Panda, K.S. Ravi Chandran / Computational Materials Science 35 (2006) 134–150

Fig. 1. Crystal structure of TiB2: (a) illustration of the stacking of the trigonal prisms; (b) illustrates the close packing of the prism; (c) shows the unit cell of TiB2 containing one Ti and two B atoms.

2M 2 V0 M2 cii ¼ 2V 0

cii ¼

D will move every bravais lattice points of the undistorted lattice to a new position in the distorted R 0 matrix, which is given by

for the normal strains ði ¼ 1 to 3Þ and for the shear strains ði ¼ 4 to 6Þ

R0 ¼ RD ð5Þ

Eq. (5) implies that other coefficients, M0, M1 and M3 are zero and this will be the case for pure linear elastic strain. However, these coefficients may not be exactly zero in actual calculations, but can be made close to zero with extremely small strains. Therefore, small elastic strains, di, i = 1–6, can be applied to the undeformed unit cell lattice with the equilibrium atomic configuration and elastic constants can be determined from the resulting change in energy. The five independent elastic constants for the hexagonal crystals are c11, c12, c13, c33, c44 and they can be determined by selectively imposing strains either individually or in combination along specific crystallographic directions. Table 1 illustrates the deformations employed in determining (c11 + c12), (c11
c12), c33 and c44, which leads to the determination of all the elastic stiffness constants. 2.3. Unit cell distortions, changes in symmetry and internal atomic displacements Considering the R matrix containing the axis vectors for the undeformed bravais lattice, any distortion matrix

ð6Þ

Here D is the distortion (for pure shear, rotation is ignored) matrix containing the deformation strain. Table 1 describes the distortion matrices employed, the crystal structure and space group after deformation and the relationship between the elastic constant and the second-order coefficient of the polynomial. The three distortions involving (c11 + c22), c33 and (2c11 + 2c12 + c13 + c33) do not change the hexagonal symmetry but are accompanied by a volume change. This is not, however, the case for distortions involving (c11
c12) and c44. In case of (c11
c12), the lattice parameters a1, a2 and c take different values after the deformation, thereby making the distorted lattice belong to a monoclinic system. Similarly in case of c44, the pure shears, in addition to changing symmetry (Table 1) also change the atomic position and the deformed positions have to be found out for the ground state energy calculations. In the unit cell distortions involved in the determinations of (c11 + c22), c33 and (2c11 + 2c12 + 4c13 + c33), the symmetry of the unit cell is preserved but the internal atomic positions change according to the distortion. As an example, Fig. 2(a) and (b) illustrates the c33

K.B. Panda, K.S. Ravi Chandran / Computational Materials Science 35 (2006) 134–150

137

Table 1 Distortion matrix and elastic stiffness relations for hexagonal TiB2 Distortion matrix 0 B @

0 B @

0

1

0

C 0A 1

1þd

0

0

1

0

1
d

0

0

1þd

0

0

1þd

0

0

1 0 B @0 1

0

0 0

1þd

1 0

0

B @0 1 0 d 0 B @

0 0

Crystal structure after deformation

Space group (symmetry) after deformation

No. of symmetry operations

Hexagonal

191 (P6/mmm)

24

Monoclinic

10 (P2/m)

8

Hexagonal

191 (P6/mmm)

24

Triclinic

1 (P1)

4

Hexagonal

191 (P6/mmm)

24

C 0A 1

1 C A

1

C dA 1

1þd

0

0

1þd

0

0

0

1þd

0

1 C A

Relationship between elastic coefficient with second-order polynomial c11 þ c12 ¼

M2 V0

c11
c12 ¼

M2 V0

c33 ¼

2M 2 V0

c44 ¼

M2 2V 0

2c11 þ 2c12 þ 4c13 þ c33 ¼

2M 2 V0

Fig. 2. c33 distortion in a hexagonal TiB2 unit cell.

distortion of the TiB2 hexagonal unit cell. The volume change is illustrated by the gray region (Fig. 2(b)). The ground state crystal energy E(V0, d3) in the deformed condition is calculated using the new lattice parameters and internal atomic coordinates. The pure shear distortions involved in the determination of c44 change the symmetry of the hexagonal crystal to a new space group in the triclinic system (Table 1).

Fig. 3(a) illustrates schematically the c44 shear deformation. Fig. 3(b) shows the three-dimensional view of the atoms in the unit cell. Fig. 3(c) shows the positions of B and Ti atoms in both the deformed and the undeformed crystal with the exaggeration of atomic displacements for clarity. The gray circles represent the atoms in the undeformed crystal. The atoms in the deformed crystal are shown in black. Both the lattice parameters

138

K.B. Panda, K.S. Ravi Chandran / Computational Materials Science 35 (2006) 134–150

Fig. 3. Pure shear distortion d4 involved in c44 in the hexagonal unit cell of TiB2. (a) The schematic of shear distortion, (b) the deformed state of the unit cell as viewed in the a direction, (c) schematic of Ti and B atom displacements (not all atoms in the same plane) and (d) illustration of parameters involved in the definition of the shear strain d4, and new atomic positions are shown in the first quadrant.

ÔbÕ and ÔcÕ and the angle ÔaÕ change as a result of the pure shear c44. The shear deformation changes the crystal symmetry from hexagonal to triclinic with the unit cell parameters of the deformed unit cell as a = 88.0173°, ˚ , b = 3.0304 A ˚, c= b = 90° and c = 120°; a = 3.03 A ˚ 3.23248 A. The shear strain d4 is p ov ow
b ¼ þ ð7Þ d4 ¼ 2eyz ¼ 2 dz dy

symmetrically under pure shear, thus leading to (ov/ dz) = (ow/dy) = eyz with the absence of rotations. Similarly, in case of (c11
c12) deformation, the symmetry of the lattice changes and the approach similar to that used for c44 was followed to determine atomic displacements due to deformation.

The displacements of internal Ti and B atoms after shear into new deformed positions were determined from the undeformed fractional coordinates, fb and fc, of an interior atom (B or Ti) in the b and c directions, respectively. The internal atomic displacements in the Y and Z directions, ov* and ox* due to the shear are then given by

 ov  1  ov ¼ ð8Þ c ¼ eyz c ¼ d4 fc c dz 2

 ow  1 ð9Þ ow ¼ b ¼ eyz b ¼ d4 fb b dy 2

The ground state crystal energies of the undeformed state and that of the deformed states of TiB2 unit cell under different distortions were determined using the computational implementation of the Hartree–Fock (HF) method in Crystal98 [9] code that uses the LCAO scheme of molecular orbital representation as well as the density functional method (DFT) in WIEN2k [10] code. The DFT implementation uses the FLAPW method that allows wave function representation in the periodic crystal in the reciprocal space with a choice of relativistic or nonrelativistic calculations using the generalized gradient approximation (GGA) of the electron density space. In the Hartree–Fock method, the molecular orbitals were constructed from a linear combination of atomic orbitals (LCAO) and were used as the basis for the wave functions. A 86-4113G [14] basis set for Ti and a 6-21G [15] basis set for B were used in the calculations. The

Thus, from lattice parameters, fractional coordinates of nonequivalent atoms and the deformation strain, the internal atom displacements required to impose the shear can be determined. It may be noted that in Eqs. (8) and (9), the crystal was considered to deform

2.4. Computational details

K.B. Panda, K.S. Ravi Chandran / Computational Materials Science 35 (2006) 134–150

-1807.2625

Energy (Rydbergs)

-1807.263 -1807.2635 -1807.264 -1807.2645

-1807.24 -1807.245

Energy (Rydbergs)

atomic basis sets for most of the elements have been developed to serve as a good starting point and the 621G basis set was chosen [16,17]. However, in performing the calculation for a specific crystal/compound, the exponents and contraction coefficients of the outer valence functions need to be optimized. In the present work, the outermost exponents for the orbitals of B were optimized using BILLY, a code developed by Towler at the University of Cambridge [18]. High tolerances in the bi-electronic coulomb and exchange series overlap between the atomic orbitals, 10
8 for coulomb and exchange overlap, the coulomb penetration, the first exchange pseudo-overlap and 10
12 for the second exchange pseudo-overlap were employed in order to achieve high numerical accuracy. The Fock matrix was diagonalized at 133 k-points corresponding to a shrinkage factor of 12 in the Monkhorst net and 793 k-points corresponding to a shrinkage factor of 24 in the Gilat net. With these choices, the total energy was converged to within 10
6 Hartrees per atom within 20–30 SCF iterations. In the DFT calculations, it is important to ensure that the ground state energy is converged within a specified tolerance in terms of the iterations for self-consistency. This can be achieved only by an accurate numerical sampling of electron densities in the Brillouin zone. This is largely dependent on the number of sampling points for numerical integrations of electron densities in the k-space in the irreducible Brillouin zone. A test was first made using the TiB2 unit cell to determine the minimum k-points needed. Fig. 4 illustrates the variation of ground state energies of equilibrium TiB2 unit cell as a function of k-points. At and above 5000 kpoints, the changes in the ground state energies were less than 9 lRy and therefore for 5000 points were used in the present set of calculations.

139

-1807.25 -1807.255 -1807.26 -1807.265 -1807.27 -1807.275 5

6

7

8

9

10

R mt * K max

Fig. 5. Rmt * Kmax optimization for TiB2.

In WIEN2K, the self-consistent ground state energy calculations are performed at two levels. First, within the cores of atoms, atomic-like functions in the basis set are used to determine the core electronic contributions to the energy whereas in the interstitial regions, Fourier-series representations are used to determine the contribution of valence electrons. The radius of the boundary between the two regions is the muffin-tin radius, Rmt and for the present calculations, the minimum of the Rmt values taken was 1.6 Bohr. It is important to check the convergence of the crystal energy as a function the term RmtKmax where Kmax is the magnitude of the largest k vector in the basis set. With a 5000 point kmesh, RmtKmax was varied within the range of 6.0–9.0 with a step size of 0.5. Fig. 5 illustrates the variation of crystal energy as a function of RmtKmax. A value of RmtKmax = 7 or greater was considered to be a reasonable approximation for the basis set in the calculations without being computationally very time consuming. All calculations of crystal energies in undeformed or deformed condition were made using the same set of parameters: Rmt = 1.6 Bohr, RmtKmax = 7, with 5000 k-points. With these choices, the total energy was converged to within 10 lRy per atom within 16–24 iterations for self-consistency of electron densities. Relaxations of internal atoms in the deformed condition were performed using the module ‘‘PORT’’ in WIEN2K, until the interatomic forces between the atoms decreased to less than 2 mRy/a.u.

-1807.265

3. Results and discussion

-1807.2655 -1807.266 0

3.1. Structure optimization 1000

2000

3000

4000

5000

6000

Number of k-points

Fig. 4. K-mesh optimization for TiB2.

7000

8000

The structure parameters for TiB2 were first optimized in WIEN2K to ensure that the calculations were

140

K.B. Panda, K.S. Ravi Chandran / Computational Materials Science 35 (2006) 134–150

Table 2 The experimental and calculated equilibrium structure parameters of TiB2 ˚ ) b (A ˚ ) c (A ˚ ) Volume (A ˚ 3) c/a a (A Experimental Present calculation

3.030 3.038

3.030 3.038

3.232 3.234

25.6972 25.8559

0.016

a

3

Expt. V = 25.697284 A o

0.014

3

Calc. V = 25.8559054 A o

1.0666 1.0645

0.012

0.008 0.006 0.004 0.002 0 23

24

25

26

27

28

29

Volume (A3/cell) 0.001

b

Expt. c/a =1.066666603 0.0008

Calc. c/a = 1.064719945

0.0006

δE (Ry)

based on equilibrium lattice constants. The additional purpose is to check whether or not the crystal energy minimum is achieved for a given set of calculation parameters such as Rmt, k-points as these are the only numerical parameters that affect the accuracy of calculations. The experimental lattice parameters (Table 2) were taken as the basis for doing the optimization by energy minimization for the hexagonal symmetry. The volume was optimized first. In WIEN2K, this is performed by applying hydrostatic strain by changing the lattice parameters while keeping the c/a ratio constant. The volume variations applied for the present study were
10%,
5%, 0%, +5% and +10% of the experimental equilibrium volume. With the optimized volume, the c/a ratio was varied to find the values at which crystal energy minimum was achieved. Fig. 6 illustrates the crystal energy variations with the structure parameters. The values corresponding to the minimum energy is listed in Table 2. There is very good agreement between the experimental values and the values obtained using WIEN2K. The calculated lattice parameters and the c/a ratio differ from the experimental values by not more than 0.3%.

δE (Ry)

0.01

0.0004

0.0002

3.2. TiB2 single crystal elastic constants The five independent elastic constants for TiB2 were determined by imposing five different deformations given in Table 1 on the equilibrium lattice of hexagonal unit cell and determining the dependence of the resulting change in energy on the deformation. Up to six deformations,
2%,
1%,
0.5%, 0.5%, 1% and 2%, were applied to determine each elastic constant. The ground state energies of TiB2 unit cell were determined for all these deformation strains. Fig. 7(a)–(j) illustrates the plots of change in total energy as a function of strain for the five distortions. The HF-method based calculations are presented in Fig. 7(a), (c), (e), (g) and (i). The calculations by the DFT-GGA method are presented in Fig. 7(b), (d), (f), (h) and (j). The lines in the figures represent the third-order polynomial fit (Eq. (4)) to the data from which the values of M2 were determined. The elastic constants were then calculated using equations in Table 1. Among the five independent elastic constants, determination of crystal energy for c44 was the most time consuming because of the triclinic unit cell with a substantial reduction in symmetry (Table 1) resulting from the deformation.

0

-0.0002 1.04

1.05

1.06

1.07

1.08

1.09

c/a Fig. 6. The structure optimization plots using GGA calculations (a) volume optimization and (b) c/a optimization.

Table 3 lists the values of the single crystal elastic constants determined in the present study. Also included are the results obtained in other theoretical works and experimental measurements. The HF results of the present work appear to be similar to that of the HF studies done by others [19]. There is a reasonable agreement between the values determined in both the studies except for c11, c13 and c33. The values of c11 and c13 (consequently c66) are significantly higher and the value of c33 is lower than that of Ref. [19]. The reason behind this is unclear because these two set of calculations were done and on different versions of CRYSTAL code (CRYSTAL98 in this study versus CRYSTAL95 in

K.B. Panda, K.S. Ravi Chandran / Computational Materials Science 35 (2006) 134–150

141

0.004

0.006

a

b

0.005 0.003

δE (Ry)

δE (Ry)

0.004

0.003

0.002

0.002

0.001

0.001 0 0

C -0.001 -0.03

11

+C

C

12

-0.02

-0.01

0

0.01

0.02

11

-0.001 -0.03

0.03

+C

12

-0.02

-0.01

δ

0

0.01

0.02

0.03

δ 0.003

0.0035

c

d

0.003

0.0025

0.0025

0.002

δ E (Ry)

δ E (Ry)

0.002 0.0015

0.0015

0.001

0.001 0.0005

0.0005

0

0

C -C

C -C 11

-0.0005 -0.03

12

-0.02

11

-0.01

0

0.01

0.02

0.03

-0.0005 -0.03

unrelaxed Relaxed

12

-0.02

-0.01

δ 0.0025

0

δ

0.01

0.02

0.0012

e

f

0.002

0.001

0.0015

0.0008

0.001

0.0006

δ E (Ry)

δ E (Ry)

0.03

0.0005

0

0.0004

0.0002

-0.0005

0

C

33

-0.001 -0.03

/2 -0.02

C -0.01

0

0.01

0.02

0.03

δ

-0.0002 -0.03

33

/2 -0.02

-0.01

0

0.01

0.02

0.03

δ

Fig. 7. Strain energy of deformation as a function of deformation parameter for five different deformations of TiB2 crystal using both HF and DFT.

Ref. [19]). It does not seem possible that the numerical implementation differences in the two versions would lead to such a large difference in c33. This aspect may be a subject of further investigation, which is beyond

the scope of the present study. It may be noted, however, that the calculated polycrystalline elastic models (E) and bulk modulus (B) on the basis of HF elastic constants are nearly equal suggesting that increases in some

142

K.B. Panda, K.S. Ravi Chandran / Computational Materials Science 35 (2006) 134–150 0.0035

0.003

Unrelaxed Relaxed B: YZ plane Relaxed B all directions

g 0.003

0.0025

h

0.0025 0.002 0.0015

δ E (Ry)

δE (Ry)

0.002

0.001

0.0015 0.001

0.0005

0.0005

0

0

2C -0.0005 -0.03

2C

44

44

-0.02

-0.01

0

0.01

δ

0.02

0.03

-0.0005 -0.03

0.01

-0.02

-0.01

0

0.01

δ

0.02

0.03

0.006

i

j 0.005

0.008

0.004

δE (Ry)

δ E (Ry)

0.006

0.004

0.003

0.002

0.002 0.001 0

0

C -0.002 -0.03

11

+C

-0.02

12

+2C

-0.01

13

+ C /2

C

33

0

0.01

0.02

0.03

-0.001 -0.03

δ

11

+C

-0.02

12

+2C

-0.01

13

+ C /2 33

0

0.01

0.02

0.03

δ

Fig. 7 (continued)

constants are compensated by the large decrease in c33 in the present study relative to that of Ref. [19]. Nevertheless, c11, c12, c13 and c66 in both studies are significantly higher than the experimental data that it did not seem (in light of the DFT-GGA results to be presented later) that any internal relaxation of atoms would result in a better agreement with the experimental data [20]. Therefore, HF calculations were not pursued further. The DFT calculations in the present study yielded results that agreed better with the experimental data than the HF results. In Table 3, the elastic constants determined in the present study are also compared with the DFT results of a previous study [21]. It should be noted that the previous DFT study [21] employed plane wave pseudo-potential (PWPP) approximation using ultrasoft pseudo-potential, while the present calculations use the full potential linearized augmented plane wave (FLAPW) method, a more recent and a more accurate method that is reported [22,23] to yield accurate electronic structure of compounds. Since a true comparison

can be based only on elastic constants with a fully relaxed atomic structure, we first focus on the effect on relaxation on elastic constants. 3.3. Effect of relaxation of atoms The unrelaxed DFT elastic constants in Table 3 are for the atomic coordinates of the deformed unit cell strictly following the linear elastic deformations imposed. However, this does not necessarily correspond to the equilibrium atomic state of the deformed unit cell, as the B–B, Ti–B and Ti–Ti bond distances were not allowed to equilibrate after deformation. However, the internal B and Ti atoms must be allowed to relax to find a lower energy minimum of the crystal in the relaxation space as schematically illustrated in Fig. 8. The true relaxation of atoms involves adjustments of atomic coordinates in all directions with the deformed unit cell external shape fixed. Additionally, the symmetry of the unit cell should not change from that of the unrelaxed

K.B. Panda, K.S. Ravi Chandran / Computational Materials Science 35 (2006) 134–150 Table 3 The calculated single crystal and polycrystalline elastic constants of hexagonal TiB2 c12

c13

c33

c44

c66

Ed

Kd

td

Method

c11

HF (this work) HF [19] DFT (FLAPW) (this work, unrelaxed) DFT (FLAPW) (this work, relaxed) DFT (PWPP) (unrelaxed) [21]a DFT (PWPP) (relaxed) [21]b Experiment [20]

874 130 137 444 282 372 786 127 87 583 271 330 654 75 99 443 344 290

679 304 0.14 670 299 0.13 625 249 0.09

650

79 100 443 256 285

570 249 0.12

659

62 100 461 260 299

585 251 0.12

655

65

99 461 260 295

582 251 0.12

660

48

93 432 260 306 565c 240 0.11

a Strain amplitude of 2%, no relaxation of internal degrees of freedom. b Strain amplitude of 2%. c Ref. [4]. d Calculated values for polycrystals: E = Elastic modulus, K = bulk modulus, t = PoissonÕs ratio.

Fig. 8. Schematic illustrating the method of relaxation. The dotted curve is the plot of energy versus deformation after relaxation. The lowest value of energy is obtained from the minima in the relaxation space.

and deformed configuration during relaxation as this may push the crystal to some other minimum energy configuration. It should be noted that c11 + c12, c33 and (2c11 + 2c12 + 4c13 + c33) distortions preserve the hexagonal symmetry and relaxations are not possible as the positions of nonequivalent Ti and B atoms are fixed in special positions required by the space group. However, (c11
c12) and c44 deformations involve symmetry reduction (Table 1). A closer look at the crystal unit cell before and after the deformation involving (c11
c12) and c44 clearly indicates that the Ti atoms will always remain in the corner positions and hence there is no relaxation is allowed for the Ti atoms. But this is not the case with B atoms and they may be relaxed in all the

143

three directions. Therefore, relaxations of B atoms were performed for the (c11
c12) and c44 deformations. We used the PORT module in WIEN2K that finds the global minimization of the total energy by minimizing the interatomic forces with the help of the first derivatives of E with respect to the atomic positions. The relaxation of atoms for (c11
c12) and c44 deformations were performed until the interatomic forces were less than 2 mRy/a.u. Table 4 illustrates the forces on atoms before and after relaxation. It is evident that for c44 distortion forces on Ti atoms are large possibly due to the fact that electron density is concentrated long B–B bonds. The c44 distortion results in a triclinic system and as evident from Table 4, relaxation leads to force minimization on all the atoms. For c11
c12, there was force only in the Z-direction due to deformations (tension in X and compression in Y direction) and the minimization is not as dramatic as in c44 probably due to the fact that the symmetry reduction was not as high as in c44. Table 5 illustrates the relaxed atomic positions for the 2% deformed state for c11
c12 and c44. As evident from the table, the crystal energy is lowered by about 700 lRy1 and 30 lRy1 for c44 and c11
c12, respectively, after relaxation. The crystal energies were recalculated on the basis of relaxed atomic position and the results are presented in Fig. 7(d) and (h). The effect of relaxation on the crystal energy for c44 is illustrated by first relaxing the B atom in the YZ plane and keeping the X-coordinates fixed. This lowers crystal energy at all strains slightly (Fig. 7(h)). Further relaxation was done by allowing the B atoms to relax in all the three directions, which resulted in further reduction in energy as seen in Fig. 7(h). The elastic constants were recalculated on the basis of the relaxed states. The relaxed single crystal elastic constants are also given in Table 3. The changes in c11, c12 and c66 are small since the force minimization in c11
c12 deformation is small and c66 depends on c11 and c12. However, the change in c44 is about 88 GPa which is significant. It can be seen that with the exception of c12, all the constants are quite close to the experimental data of Spoor et al. [20]. The present data may be compared with the PWPP calculations of Milman et al. [21] were both elastic constants based on unrelaxed and relaxed atomic positions were determined. However, that work did not report the details of relaxation method and most elastic constants did not show any change beyond the

1 These values can be used to approximately estimate the accuracy of the elastic constants calculated here. The lowering of energies by 700 lRy and 30 lRy corresponding to relaxations of c44 and c11 leads to changes of 88 GPa (c44) and 4 GPa (c11 or c12) in elastic constants. These values correspond to changes of about 0.12 GPa (c44) and 0.13 GPa (c11 or c12) for 1 lRy. The present calculations were performed with an SCF convergence tolerance of 10 lRy. This means an uncertainty of about 1 GPa in the elastic constants.

144

K.B. Panda, K.S. Ravi Chandran / Computational Materials Science 35 (2006) 134–150

Table 4 Forces on atoms before after relaxation for c44 and c11
c12 deformations Parameter

Atom

Forces (unrelaxed) (mRy/a.u.)

Forces (relaxed) (mRy/a.u.)

Fx

Fy

Fz

Fx

Fy

Fz

c44

Ti: B1:

10.882
5.78


7.694 3.51

12.717
2.2


0.887
0.134

1.845
1.451

1.516
1.989

c11
c12

B2: Ti: B


4.999 0 0

4.092 0 0


10.405 0
3.661

1.004 0 0


0.375 0 0

0.517 0
1.107

Table 5 The B atom position before and after relaxation with an applied deformation of 2% for c11
c12 and c44 deformations Deformation

Boron atom positions (unrelaxed)

Boron atom positions (relaxed)

Energy (Ry) (unrelaxed)

Energy (Ry) (relaxed)

c44

B1: (0.3425, 0.66, 0.5053) B2: (0.6757, 0.33, 0.5107)

B1: (0.333, 0.667, 0.5106) B2: (0.6659, 0.334, 0.4999)


1807.262691


1807.263420

c11
c12

B1: (0.33, 0.66, 0.5) B2: (0.66, 0.33, 0.5)

B1: (0.3324, 0.6662, 0.5) B2: (0.66753, 0.3337, 0.5)


1807.263541


1807.263572

reported statistical error. In particular, the elastic constant c44 is identical with and without relaxation, which is strange. Generally, the extent of relaxation is dependent on deformation and symmetry reduction. Distortions involving the largest reduction in symmetry have the greatest potential to relax. That is particularly true in the case of c44 deformation in TiB2 where the symmetry changes from hexagonal to triclinic. The relaxation produced large reductions in interatomic forces (Table 4) in c44 deformation compared to c11
c12. This is also due to the fact that the restriction of special positions of B atoms in hexagonal lattice is no longer a constraint in the space group after deformation. Thus, although Milman et al. show c11
c12 relaxing to a smaller degree, which agrees with our calculations, the fact that they did not see any relaxation for c44, which, in fact is expected to relax significantly, was not explained in that study. Therefore, one has to conclude that the present calculations are more consistent with expectations than that of Milman et al. This trend can also be seen in comparing the present data with the experimental results. With the exception of c12, all constants calculated in this study are quite close to the experiment data. In particular, the value of c33 determined here (443 GPa) is close to the experimental value of 432 GPa whereas Milman et al.Õs value (461 GPa) is significantly higher. The agreement of Milman et al.Õs polycrystalline E and B values with the experimental data is also surprising in view of the lack of any relaxation in c44 deformation and the large values of c33 in their work relative to present data. This may be a consequence of mutual cancellation of errors in individual elastic constants. Therefore, it may be concluded that the present elastic constant calculations are more reliable.

3.4. Elastic properties of TiB2 polycrystalline aggregate The theoretical polycrystalline elastic modulus for TiB2 can be determined from the set of five independent elastic constants. There are two approximation methods to calculate the polycrystalline modulus, namely the Voigt method [24] and the Reuss [25] method. The Voigt and Reuss assumptions result in the theoretical maximum and minimum values of the isotropic elastic modulus, respectively. They are also commonly referred to Voigt–Reuss bounds meaning the extremes within which the polycrystalline elastic modulus should fall. For hexagonal TiB2, the Voigt (KV) and Reuss (KR) bulk moduli are given by 1 K V ¼ ½2ðc11 þ c12 Þ þ c33 þ 4c13

9 KR ¼

ðc11 þ c12 Þc33
2c213 c11 þ c12 þ 2c33
4c13

ð10Þ ð11Þ

Similarly, the upper and the lower bounds for the shear modulus of polycrystalline TiB2 aggregate according to Voigt and Reuss approximations are given by 1 ðc11 þ c12 þ 2c33
4c13 þ 12c44 þ 12c66 Þ ð12Þ 30

2 ðc11 þ c12 Þc33
2c213 c44 c66 5 GR ¼ 2 3K V c44 c66 þ ððc11 þ c12 Þc33
2c213 Þ2 ðc44 þ c66 Þ GV ¼

ð13Þ where 1 c66 ¼ ðc11
c12 Þ 2

ð14Þ

K.B. Panda, K.S. Ravi Chandran / Computational Materials Science 35 (2006) 134–150

The arithmetic average of the Voigt and the Reuss bounds is called the Voigt–Reuss–Hill (VRH) average [26] and is commonly used to estimate elastic moduli of polycrystals. The VRH averages for shear modulus (G) and bulk modulus (K) are 1 1 G ¼ ðGR þ GV Þ; K ¼ ðK R þ K V Þ ð15Þ 2 2 The isotropic bulk modulus and shear modulus values for polycrystalline TiB2 were first calculated using the above equations. The polycrystalline elastic modulus (E) and the PoissonÕs ratio (m) were then computed from these values using the following relationship: E¼

9KG ; 3K þ G

t¼

3K
2G 2ð3K þ GÞ

ð16Þ

Table 3 lists the polycrystalline elastic modulus, bulk modulus and the PoissonÕs ratio of TiB2 calculated in this study. The values computed in the present study based on the HF calculations as well as the DFT calculations both before and after relaxations, can be compared with experimental [20] data and other calculations [19,21]. The E of polycrystalline TiB2 based on HF method is (679 GPa) quite unrealistic. However, the internal relaxation of atoms in DFT reduced the polycrystalline E value from 625 GPa in the unrelaxed condition to 570 GPa after relaxation. The experimental elastic modulus values of polycrystalline TiB2 are in the range of 565–569 GPa [4]. After correction for temperature, these values are likely to be in the range of 570–575 GPa. Thus, the 570 GPa as determined by DFT-GGA with relaxation of interatomic forces in this study is in excellent agreement with the experimental data. This agreement is also better compared to the polycrystalline E value of 582 GPa calculated by Milman et al. 3.5. Elastic anisotropy of TiB2 Most of the crystals exhibit elastic anisotropies of varying degree and there have been different ways to represent the elastic anisotropy of crystals. The compression and the shear anisotropic factors provide measures of the degrees of anisotropy in atomic bonding in different crystallographic planes. For a transversely isotropic material (hexagonal), the anisotropies in compressibility and in shear are given by Acomp ¼

s33 þ 2s13 s11 þ s12 þ s13

ð17Þ

Ashear ¼

2c44 c11
c12

ð18Þ

The shear anisotropic factors calculated in the present work are given in Table 6. A value of unity means that the crystal exhibits isotropic properties and values

145

Table 6 Various measures of anisotropy in TiB2 as determined by calculations of the present study and the experimental data Method

Acomp

Ashear

Acomp (%)

Ashear (%)

Experiment [20] GGA

1.54 1.1467

0.85 0.896

1.27 1.25

1.54 1.0946

other than unity represent varying degree of anisotropy. From the table it can bee seen that TiB2 exhibits very little anisotropy. Another way of measuring the elastic anisotropy is given by the percentage of anisotropy in the compression and shear [27]. They are defined as Acomp ¼

KV
KR  100 KV þ KR

ð19Þ

Ashear ¼

GV
GR  100 GV þ GR

ð20Þ

These values can range from zero (isotropic) to 100% representing the maximum anisotropy. The percentage anisotropy values were computed for TiB2 and are shown in Table 6 along with the values estimated by Spoor et al. [20]. It can be seen that the anisotropy both in shear and compression is almost nonexistent. An illustrative way of describing the elastic anisotropy is a three-dimensional surface representation showing the variation of elastic modulus with direction [28]. This directional dependence of the elastic modulus for hexagonal crystals is given by the following relationship: 1 ¼ ð1
l23 Þ2 s11 þ l43 s33 þ l23 ð1
l23 Þð2s13 þ s44 Þ E

ð21Þ

where s11, s22, etc., are the elastic compliance constants and l1, l2 and l3 are the direction cosines. Fig. 9(a) and (b) illustrate the directional dependence of the TiB2 in the elastic modulus calculated using the compliance constants computed from the present calculation as well as that obtained by experiment by Spoor et al. The projection of the elastic constants on the XY and YZ planes are shown in Fig. 9(c) and (d), respectively. It is evident from the plots that the computed data using GGA in the present work agrees very well with the experimental data. The variation in magnitude of the elastic constant with respect to direction would indicate the overall anisotropy exhibited by TiB2. A deviation from a spherical shape indicates the degree of anisotropy in the system. As can be seen, Fig. 9(a) and (b) show a small deviation from a spherical shape and hence one can conclude that TiB2 does not exhibit a high degree of anisotropy. However, while in-plane anisotropy in XY plane is nonexistent (Fig. 9(c)), a significant in-plane elastic anisotropy in YZ plane is revealed in Fig. 9(d). This is consistent with the hexagonal symmetry of XY plane where Ti–Ti and B–B bonds exist in ‘‘graphite-like’’ layers alternatively stacked in the c-direction and the fact

146

K.B. Panda, K.S. Ravi Chandran / Computational Materials Science 35 (2006) 134–150

Fig. 9. Directional dependence of the elastic modulus in TiB2: (a) Present DFT-FLAPW-GGA calculations using Wien2K; (b) Experimentally obtained by Spoor et al. [20]. The plane projections of the directional dependence of the elastic modulus are shown for comparison: (c) projections on XY-plane and (d) on YZ-plane. The units shown are in GPa.

that weak Ti–B bonding is observed, as will be discussed later. The elastic modulus in the c-direction is significantly higher than in other two directions. The directional dependence of the bulk modulus and the shear modulus can be calculated using the following relationship: 1 ¼ ðs11 þ s12 þ s13 Þ
ðs11 þ s12
s13
s33 Þl23 K  1 s44 ¼ s55 þ s11
s12
ð1
l23 Þ G 2 þ 2ðs11 þ s33
2s13
s44 Þð1
l23 Þl23

agree very well with the data of Spoor et al. (Fig. 10(a) and (b)). The in-plane bulk modulus anisotropy in XY plane is nonexistent (Fig. 10(c)) and it is large in YZ plane (Fig. 10(d)), similar to that observed for the elastic modulus. Similar observations were made with the shear modulus anisotropy and are not presented here.

ð22Þ 3.6. Analysis of chemical bonding

ð23Þ

Fig. 10(a)–(d) illustrates the directional dependence of the bulk modulus by our GGA calculation as well as that by Spoor et al. [20]. It is evident that bulk modulus is nearly isotropic and the present calculations

The minimal degree of anisotropy and the large values of elastic constants exhibited by TiB2 can be understood by examining the nature of chemical bonding between the atoms. The electro-negativity of Ti (1.54) and B (2.04) would suggest a directional bonding between Ti and B. In order to understand the bonding between the B–B, Ti–B and Ti–Ti atoms, the density

K.B. Panda, K.S. Ravi Chandran / Computational Materials Science 35 (2006) 134–150

147

Fig. 10. Directional dependence of the bulk modulus in TiB2: (a) Present DFT-FLAPW-GGA calculations using Wien2K, (b) Experimentally obtained by Spoor et al. [20]. The plane projections of the directional dependence of the bulk modulus are shown for comparison: (c) projections on XY-plane and (d) on YZ-plane. The units shown are in GPa.

ger than between B and Ti, even though some hybridization between B and Ti is evident in Fig. 11 as 4 Total

3.5

Fermi Level (Ef)

Ti-3d B-2p

3

B-2s

DOS ( /eV )

of states (DOS) for the valence bands and electron charge density maps along different crystallographic planes were made from the electron density data of the undeformed crystal. Fig. 11 illustrates the total density of states and the projected density of states for Ti-3d, B-2s and B-2p orbitals. The sharp valley in DOS around the Fermi level (Ef) is called the pseudo-gap, which can be present owing to either a large difference in ionicity [29]. In Fig. 11, the density of states is largely contributed by B-2p levels below Ef and by Ti-3d levels above Ef. The correspondence of B-2p and Ti-3d peaks at E =
3.2 eV as well as at E = 2.8 eV suggests that the hybridization between B-2p and Ti-3d orbitals occur both below and above the Fermi level. Additionally, the fact that B-2p states contribute most to the total DOS below Ef is consistent with the higher electronegativity of B relative to Ti. This indirectly suggests that localization of electrons between B atoms is much stron-

2.5 2 1.5 1 0.5 0

-12

-8

-4

0

Energy (eV) Fig. 11. Total and projected density of states in TiB2.

4

148

K.B. Panda, K.S. Ravi Chandran / Computational Materials Science 35 (2006) 134–150

pointed out before. In fact, a recent experimental investigation [30] showed the existence of broad bands of charge accumulation between B atoms confirming the above observation. These observations suggest a ‘‘graphite-like’’ structure of TiB2, with the B atoms at Z = 0.5 positions forming a hexagonal network of covalent bonds on (0 0 2) plane. However, regarding the nature of intersheet Ti–B bonding two different viewpoints exists in literature [30,31]. The experimental work [30] concluded that there is no charge transfer from Ti to B and claimed to be consistent with the cluster calculations [32] that showed a weak Ti–B bond index of 0.137 (in the scale of 0–1 with 1 being the index of the strongest covalent bond). This suggests that there is no significant intersheet Ti–B bonding. However, the theoretical study [31] concluded that there is significant Ti–B intersheet bonding exists based on B-2p and Ti-3d hybridization. The present calculations, based on DOS data are consistent with the latter study as well as two other studies, one based on tight binding linear muffin tin orbital method [29] and the other based on discrete-variational Xa method [33]. However, the weak Ti–B bonding however small, seems to be sufficient and is essential to distinguish TiB2 from layered hexagonal compounds like

graphite. The most striking feature of the present calculation is the large population of Ti-3d states at about E = 2.8 eV above the Fermi level, which probably explains the electrical conductivity of TiB2. In order to illustrate the three-dimensional electron density topography, we have considered three planes revealing the topology of Ti–Ti, Ti–B and the B–B bonding. Figs. 12–14 illustrate the electron charge density maps in different crystallographic planes. Fig. 12(a) illustrates the two-dimensional electron charge density distribution in the basal (0 0 1) plane, which consists of four Ti atoms. The three-dimensional electron density distribution is shown in Fig. 12(b). From Fig. 12(b) it can be clearly seen that there is no strong bonding between Ti–Ti atoms as there is no accumulation of charges in excess of the background charge density of the plane. It appears that the Ti–Ti bonding is predominantly metallic in TiB2, a state that is created by the overlap of bonding d-orbitals of Ti atoms in the (0 0 1) plane. This is consistent with the observations of heavily populated Ti-3d states above Fermi level and the fact that TiB2 is electrically conductive. The nature of Ti–B and B–B bonding is illustrated in Fig. 13 with the help of the (1010) plane. Fig. 13(a)

Fig. 12. (a) Electron charge density distribution in the (0 0 1) plane in the TiB2 unit cell illustrating the Ti–Ti bond and (b) 3D plot showing the bonding charge density in the (0 0 1) plane.

Fig. 13. (a) Electron charge density distribution in the (1 0  1 0) plane in the TiB2 unit cell illustrating the Ti–B bond and (b) 3D plot showing the bonding charge density in the (1 0 1 0) plane.

K.B. Panda, K.S. Ravi Chandran / Computational Materials Science 35 (2006) 134–150

149

Fig. 14. (a) Electron charge density distribution in the (0 0 2) plane in the TiB2 unit cell illustrating the B–B bond and (b) 3D plot showing the bonding charge density in the (0 0 2) plane.

shows the two-dimensional contour electron density maps. Fig. 13(b) shows the details of the three-dimensional electron charge distribution. It can be seen that there is a small charge accumulation of electron charge between the Ti–B atoms, which is consistent with the idea of B-2p and Ti-3d hybridization noted in Fig. 11. However, the charge accumulation along Ti–B bond is much smaller than that along the B–B bond. This clearly confirms the B–B bond as covalent and Ti–B bond as ionic (because B is more electronegative than Ti and the net charge transfer from Ti to B can be only larger than that from B to Ti). This perhaps explains the high value of elastic constant c11 (650 GPa) and a relatively lower value of elastic constant c33 (443 GPa). This difference also seems to explain the significant in-plane anisotropies of elastic and bulk moduli, shown in Figs. 9(d) and 10(d). Finally, the nature of ‘‘graphite-like’’ hexagonal network of B–B bonding in the (0 0 2) plane is illustrated in Fig. 14(a) in the orientation where the two B atoms having the shortest interatomic distance. The three-dimensional charge density distribution shown in Fig. 14(b) highlights the nature of charge density accumulation between B–B bonds. In the figure, the highly directional nature of bonding and the ‘‘graphite-like’’ B network, as evidenced by the large amount of charge accumulation in the three prongs between the B atoms, can be seen. This indicates that there is a strong covalent interaction between the B–B atoms, supporting the fact that the creation pseudo-gap around the Fermi level is due to the strong hybridization effect of the B(p)–B(p) electron interactions between the B atoms themselves. It can be concluded that the nature of bonding between Ti–B is quite weak compared to that in B–B and its contribution to the pseudo-gap creation is fairly small. Additionally, it may be recalled that the weak metallic bonding dominating Ti–Ti bond is confirmed from Fig. 12(a). Thus the strength of bonding in TiB2 can be ranked in the ascending order as: Ti–Ti (metallic), Ti–B (ionic) and B–B (covalent).

4. Conclusions (1) The HF based calculations consistently yield higher values of anisotropic elastic constants for TiB2 and these values are anywhere from 5% to 20% higher than the DFT-FLAPW-GGA based calculations and experimental values. (2) The DFT-FLAPW-GGA calculations of the present study yielded the single crystal elastic constants that are in very good agreement with the experimental data. The values of the individual elastic constants in GPa are: c11 = 650, c12 = 79, c13 = 100, c33 = 443, c44 = 256 and c66 = 285. The values are estimated to be accurate to about 1 GPa. (3) Internal relaxation of atoms after deformation has been found to have a significant effect on the single crystal elastic constants. The largest change upon relaxation was seen in c44, where the B atom tends to relax moving away from their undeformed equilibrium configuration. We have observed very little relaxation in c11
c12. (4) The calculation also yielded polycrystalline elastic, bulk and shear moduli that are in very good agreement with the experimental values. The spatial variations of moduli are also in good agreement with the experimental data. (5) TiB2 exhibits very little anisotropy. The percentage anisotropy factor observed for the shear and compression modes is around 1% suggesting that the elastic properties are largely isotropic. However, significant in-plane elastic anisotropies exist in YZ-plane, caused by the large difference in bond strength between B–B and Ti–B bonds. (6) The bonding observed in TiB2 has a mixed character. The strength and character of bonding in TiB2 can be ranked in ascending order as: Ti–Ti (metallic), Ti–B (ionic) and B–B (covalent). It is concluded that the highest value of elastic constant

150

K.B. Panda, K.S. Ravi Chandran / Computational Materials Science 35 (2006) 134–150

c11 (650 GPa) is due to B–B covalency while c33 (443 GPa) is relatively lower because Ti–B bonding is ionic and weak.

Acknowledgements The authors would like to thank Dr. Mike Towler at the Cavendish Laboratory, University of Cambridge for help with the basis set optimization, and Dr. P. Blaha, Vienna University of Technology, Austria; Stefaan Cottenier, Instituut voor Kern-en Stralingsfysica, Belgium for providing several clarifications during the calculations. Victor Bazterra of Center for High Performance Computing (CHPC), University of Utah provided invaluable help with respect to calculations using Crystal98 code. References [1] B.F. Decker, R. Kasper, The crystal structure of TiB, Acta Cryst. 7 (1954) 77–80. [2] R.G. Fenish, A new intermediate compound in the titanium– boron system, Ti3B4, Trans. Metall. Soc. AIME 236 (1966) 804. [3] T. Lundstrom, in: V.I. Matkovich (Ed.), Boron and Refractory Borides, Springer-Verlag, 1977. [4] R.G. Munroe, Material properties of titanium diboride, J. Res. Natl. Inst. Stand. Technol. 105 (2000) 709–720. [5] P. Blaha, K. Schwarz, J. Redinger, Electron density of TiC, J. Phys. F: Met. Phys. 15 (1985) 263–266. [6] S. Jhi, J. Ihm, S.G. Louie, M.L. Cohen, Electronic mechanism of hardness enhancement in transition-metal carbonitrides, Nature 399 (1999) 132–134. [7] B.B. Karki, L. Stixrude, S.J. Clark, M.C. Warren, G.J. Ackland, J. Crain, Elastic properties of orthorhombic MgSiO3 perovskite at lower mantle pressure, Am. Mineral. 82 (1997) 635–638. [8] B.B. Karki, L. Sixtrude, R.M. Wentzcovitch, High pressure elastic properties of major materials of earthÕs mantle from first principles, Rev. Geophys. 39 (2001) 507–534. [9] R. Dovesi, C. Pisani, C. Roetti, M. Causa, V.R. Saunders, CRYSTAL 88, An ab-initio all-electron LCAO-Hartree–Fock program for periodic systems, improved and updated version of the original copyrighted CRYSTAL code, QCEP pgm no. 577, Quantum Chemistry Program Exchange, Indiana University, Bloomington, Indiana, 1989. [10] P. Blaha, K. Schwarz, G.K.H. Madsen, D. Kvasnicka, J. Luitz, An augmented plane wave + local orbitals program for calculating crystal properties, Karlheinz Schwarz, Techn. Universitat Wien, Austria, 2001, ISBN 3-9501031-(1-2). [11] J.D.H. Donnay, H.M. Ondik (Eds.), Crystal Data Determinative Tables, II: Inorganic Compounds, third ed., Joint Committee on Powder Diffraction Standards, 1973, pp. H-116–H-117.

[12] B. Post, F.W. Glaser, D. Moskowitz, Transition metal diborides, Acta Metall. 2 (1954) 20. [13] D.C. Wallace, Thermodynamics of Crystals, Wiley, NY, 1972. [14] V.R. Saunders, unpublished work. [15] M. Causa, A. Zupan, Density functional LCAO calculations of periodic systems. A posteriori correction of the Hartree–Fock energy of covalent and ionic crystals, Chem. Phys. Lett. 220 (1994) 145–153. [16] J.S. Binkley, J.A. Pople, W.J. Hehre, Self-consistent molecular orbital methods. 21. Small split-valence basis sets for first-row elements, J. Am. Chem. Soc. 102 (1980) 939–945. [17] M.M. Francl, W.J. Pietro, W.J. Hehre, J.S. Binkley, M.S. Gordon, D.J. DeFrees, J.A. Pople, Self-consistent molecular orbital method: a polarized basis set for second row elements, J. Chem. Phys. 77 (1982) 3654–3665. [18] M. Towler, Available from: . [19] C.A. Perottoni, A.S. Pereira, J.A.H. da Jornada, Periodic Hartree–Fock linear combination of crystalline orbitals calculation of the structure, equation of state and elastic properties of titanium diboride, J. Phys.: Condens. Matter 12 (2000) 7205–7222. [20] P.S. Spoor, J.D. Maynard, M.J. Pan, D.J. Green, J.R. Hellman, T. Tanaka, Elastic constants and crystal anisotropy of titanium diboride, Appl. Phys. Lett. 70 (1997) 1959–1961. [21] V. Milman, M.C. Warren, Elastic properties of TiB2 and MgB2, J. Phys.: Condens. Matter 13 (2001) 5585–5595. [22] W.A. Hofer, J. Redinger, P. Podloucky, Modeling STM tips by single absorbed atoms on W(1 0 0) films: 3d and 4d transitionmetal atoms, Phys. Rev. B 64 (2001) 125108-1–125108-7. [23] N.A. Hill, Density functional studies of multielectron magnetoelectrics, Ann. Rev. Mater. Res. 32 (2002) 1–37. [24] W. Voigt, Lehrbook Der Kristallphysik, second ed., Teubner, Leipsig, 1928. [25] A. Reuss, Berechnung der Fliessgrenze von Mischkristallen auf Grund der Plastizitaettsbediengung fuer EinKristalle, Z. Angew. Math. Mech. 9 (1929) 49–58. [26] R. Hill, The elastic behavior of a crystalline aggregate, Proc. Phys. Soc. London, Sect. A. 65 (1952) 349–354. [27] D.H. Chung, W.R. Buessem, in: F.W. Vahldiek, S.A. Mersol (Eds.), Anisotropy in Single Crystal Refractory Compounds, Plenum, New York, 1968, pp. 217–245. [28] J.F. Nye, Physical Properties of Crystals, Oxford University Press, Oxford, 1985. [29] P. Vajeeston, P. Ravindran, C. Ravi, R. Asokamani, Electronic structure bonding and ground-state properties of AlB2-type transition-meta diborides, Phys. Rev. B 63 (2001) 045115-1– 045115-12. [30] G. Will, Electron deformation density in titanium diboride chemical bonding in TiB2, J. Solid State Chem. 177 (2004) 628– 631. [31] K. Lie, R. Hoier, R. Brydson, Theoretical site-and symmetryresolved density of states and experimental EELS near-edge spectra of AlB2 and TiB2, Phys. Rev. B 61 (2000) 1786–1794. [32] P.G. Perkins, A.V.J. Sweeney, An investigation of the electronic band structures of NaB6, KB6, TiB2 and CrB, J. Less-Common Met. 165 (1977) 165–173. [33] M. Mizuno, I. Tanaka, H. Adachi, Chemical bonding in titaniummetalloid compounds, Phys. Rev. B 59 (1999) 15033–15047.