First-principles study of the lattice dynamics of Sb₂S₃

Abstract

We present a lattice dynamics study of orthorhombic antimony sulphide (Sb₂S₃) obtained using density-functional calculations in conjunction with the supercell force-constant method. The effect of Born effective charges is taken into account using a mixed-space approach, resulting in the splitting of longitudinal and transverse optical (LO-TO) phonon branches near the zone center. Zone-center frequencies agree well with Raman scattering experiments. Due to the slow decay of the interatomic force constants (IFC), a minimal 2 × 4 × 2 supercell (Pnma setting) with 320 atoms is crucial for an accurate determination of the dispersion relations. Smaller supercells result in artificial acoustic phonon softening and unphysical lifting of degeneracies along high symmetry directions. We propose a scheme to investigate the convergence of the IFC with respect to the supercell sizes. The phonon softening can be attributed to the periodic images that affect the accuracy of the force constants, and the truncation of long-ranged forces. The commensuration of the q-vectors with the supercell size is crucial to preserve degeneracies in Sb₂S₃ crystals.

---

1 Introduction

Sb₂S₃ belongs to the group of metal chalcogenides (A₂B₃, A = As, Sb, Bi and B = S, Se, Te) that form an important class of semiconductors with extensive applications in photovoltaics^1,2 and optoelectronics.³ They hold great promise as photovoltaic converters and thermoelectric cooling devices^4,5 due to their small direct bandgaps, high thermoelectric power, and high absorption coefficient in the visible region.^6,7 There is also a surge of interest in using Sb₂S₃ as a solid-state semiconductor-sensitized solar cell to replace the inorganic dye in dye-sensitized solar cells.^7,8 Sb₂S₃ has been synthesized in various nanostructured forms, such as nanowires and nanotubes,^9–11 which exhibit enhanced ferroelectric, piezoelectric, and conductive properties.

So far, much focus has been placed on the synthesis and electronic properties^12–15 of Sb₂S₃. To complement these known aspects, we present here a study of the lattice vibrational properties of Sb₂S₃ using density-functional theory (DFT). Phonon dispersion is one of the fundamental properties of crystals. The behavior of the branches reflects specific features of the crystal structure and the interactions between the constituent atoms. These lattice dynamical properties are indispensable in order to understand the properties of interest for device engineering and design, such as electronic transport and lattice specific heat.

We obtain phonon dispersions of Sb₂S₃ using the supercell force-constant method,^16–20 fully taking into account the effect of Born effective charges. We note that a similar approach has been employed to study Bi₂S₃ in an earlier work.²¹ Due to the slow decay of IFC, we also investigate the effects of supercell sizes on the accuracy of dispersion relations, and propose a scheme to investigate the convergence of the IFC with respect to supercell sizes.

2 Methods

The orthorhombic phase of Sb₂S₃ belongs to the space group Pnma (#62) containing 20 atoms per primitive cell, five of which are inequivalent (Fig. 1). DFT calculations are carried out using the Quantum ESPRESSO²² suite within the local density approximation (LDA). We use pseudopotentials S.pz-bhs.UPF and Sb.pz-bhs.UPF from http://www.quantum-espresso.org. The electronic wavefunctions are expanded in a plane-wave basis set with a kinetic energy cutoff of 75 Ry. The Monkhorst–Pack²³k-point sampling scheme used for Brillouin zone (BZ) integration has divisions of less than 0.03 Å⁻¹. The total energies are tested to converge within 1 meV per atom.

Fig. 1 Ball-and-stick model of a Sb₂S₃ primitive cell containing 12 S (yellow) and 8 Sb (teal) atoms shown in the Pnma setting. Translucent atoms outside the unit cell have been added to complete the four formula units, each of which contains five inequivalent atoms.

We relax the atomic coordinates and cell dimensions using a Broyden–Fletcher–Goldfarb–Shanno quasi-Newton algorithm to obtain optimized structures with residual forces of less than 10⁻³ eV Å⁻¹ and stresses of less than 10⁻⁴ eV Å⁻³. The obtained equilibrium lattice parameters and atomic coordinates of inequivalent atoms are reported in Table 1 which show good agreement with experimental values. Note that the theoretical lattice constants a, b and c tend to be smaller than the experimental values, which is expected from a LDA calculation.^13,14

Table 1 Equilibrium lattice parameters and inequivalent atomic positions for orthorhombic Sb₂S₃ in the Pnma setting. The experimental data²⁴ are included for comparison. Deviations in lattice constants are smaller than 3.8%

Lattice constants (Å)			Atomic positions
	This work	Expt.		This work			Expt.
				x/a	y/b	z/c	x/a	y/b	z/c
a	11.081	11.299	Sb1 (4c)	0.0200	0.25	0.6741	0.0290	0.25	0.6738
b	3.831	3.831	Sb2 (4c)	0.3437	0.25	0.4677	0.3503	0.25	0.4641
c	10.805	11.227	S1 (4c)	0.0508	0.25	0.1282	0.0493	0.25	0.1226
b/a	0.346	0.339	S2 (4c)	0.3738	0.25	0.0567	0.3745	0.25	0.0612
c/a	0.975	0.994	S3 (4c)	0.2133	0.25	0.8068	0.2077	0.25	0.8071

---

To obtain the full phonon dispersion within the BZ, we adopt a direct, supercell force-constant approach.^17–20,25 We note that density-functional perturbation theory (DFPT)^26,27 may also be employed to obtain the phonon dispersions. The relative merits of these two methods have been discussed in ref. 26. In the supercell force-constant approach, the sth atom in the primitive cell is displaced from its equilibrium position in the α direction by ±δ^α_s = ±0.015 Å. The forces acting on the uth atom in the supercell, F^β_u(±δ^α_s) are calculated using the Hellmann–Feynman theorem. α and β denote the three spatial directions. We then use a finite central-difference scheme to evaluate the matrix elements of the IFC ϕ as

The dynamical matrix at a q-point is obtained by summing the contributions from all atoms in the supercellwhere s and t run over all atoms in the primitive cell, M_s is the mass of the sth atom and R is a lattice translation vector. ϕ^αβ_st(R) denotes the IFC between the sth atom in the primitive cell and another atom in the supercell that is the image of the tth atom in the primitive cell under R. The diagonalization of D(q) then yields the phonon frequencies at q.

2.1 Symmetry reduction

In order to reduce the number of static DFT calculations, we use the symmetry properties of the crystals to transform the forces.¹⁷ Only the inequivalent atoms within the primitive cell are displaced to find the forces F_su between the sth atom in the primitive cell and uth atom in the supercell, which can be represented by a 3 × 3 matrix. To obtain the forces between an equivalent atom and all other atoms in the supercell, we use the space group operation S that maps the sth inequivalent atom to the p_sth equivalent atom in the primitive cell. The forces between p_sth and p_uth atoms can be simply calculated as

F_pspu = G(S)F_suG(S⁻¹) (3)

where G(S) represents the point group part of S in Cartesian coordinates. This approach allows us to displace five atoms in the three spatial directions rather than all 20 atoms in the primitive cell, resulting in substantial saving of the total calculation time.

2.2 Non-analytical correction for q → 0

Due to the polar character of Sb₂S₃, the long-range dipole–dipole interaction gives rise to a macroscopic electric field that affects longitudinal optical (LO) phonon modes and not the transverse optical (TO) modes.²⁸ The LO-TO splitting depends on the direction from which one approaches the Γ point in the BZ. This effect is reflected in the non-vanishing Born effective charge tensor Z*, taking the form of a non-analytical contribution D̃^αβ_st to the dynamical matrix^26,27,29,30 in the limit q → 0:where Ω is the volume of the primitive cell, e is the elementary charge, ε^∞ is the high-frequency dielectric tensor and is the Born effective charge tensor for the sth atom. Z* and ε^∞ may be calculated using density-functional perturbation theory (DFPT).^26,27

In order to include the non-analytical correction in the phonon dispersion, we add³¹ a correction factor [small phi, Greek, tilde] to the real space force-constant ϕ:

Φ^αβ_st(R) = ϕ^αβ_st(R) + [small phi, Greek, tilde]^αβ_st (5)

where Φ(R) is the corrected real-space inter-atomic force-constant matrix. We may calculate the correction factor [small phi, Greek, tilde]^αβ_st by imposing the condition thatfrom which we obtainwhere N is the number of primitive cells in the supercell. The corrected force-constant matrix Φ (eqn (5)) is finally used to calculate the phonon frequencies.

3 Results and discussion

3.1 Born effective charges

The Born effective charge is the first derivative of the macroscopic polarization along the α direction with respect to the displacement of the sth atom along the β direction. This quantity is calculated using linear response theory at the zone center (Γ), and the non-zero values are shown in Table 2. The diagonal elements are different for each Sb and S atom, and off-diagonal elements are present, showing considerable anisotropy in the system.

Table 2 Non-zero components of the Born effective charge tensor Z* and high-frequency dielectric tensor ε^∞ of Sb₂S₃. Components of the five inequivalent atoms are shown for Z*. The values of the equivalent atoms are the same up to a sign and can be calculated using the transformation relation according to eqn (3)

	xx	yy	zz	xz	zx
ε ^∞	9.33	18.7	13.0	0	0
Z*(Sb1)	2.89	5.62	7.36	0.07	1.53
Z*(Sb2)	3.33	7.25	4.50	0.28	0.09
Z*(S1)	−2.35	−4.18	−4.07	1.07	0.83
Z*(S2)	−1.83	−4.80	−4.44	−0.45	−0.33
Z*(S3)	−2.03	−3.90	−3.36	−0.20	−1.34

---

The formal valence charges for Sb and S are + 3 and −2. Our calculation shows maximum effective charges of +7.36 and −4.80 for Sb and S respectively. From the study of ferroelectric compounds, it has been suggested that ions with effective charges close to the formal valence charge behave as closed-shell ions. Conversely, the presence of covalent character in the bonds causes a large amount of delocalized charge to flow through the structure during lattice displacements.^28,32,33 In Sb₂S₃, the significantly larger effective charges suggest that there is substantial covalent character in the bonds. Sb atoms are able to donate electrons to S atoms during lattice displacements, and hence increase the magnitude of their respective Born effective charges. These results agree with X-ray photoelectron spectroscopy studies that describe the bonding in Sb₂S₃ as tight covalent.³⁴

We also obtain the high-frequency dielectric tensor ε^∞ which is diagonal as shown in Table 2. The anisotropy of ε^∞ is an indication of the anisotropy of Sb₂S₃ structure.

3.2 Zone-center phonons

Since Pnma is a centrosymmetric space group, the Raman and infra-red (IR) modes of Sb₂S₃ are mutually exclusive, i.e., a mode cannot be simultaneously Raman and IR active. There are 60 phonon modes at Γ that respect the D_2h point group symmetry:

Γ = 3Γ_acoustic + 30Γ_Raman + 22Γ_IR + 5Γ_silent

3 are acoustic phonon modes (Γ_acoustic = B_1u + B_2u + B_3u). Of the optical phonon modes, 30 are Raman active (Γ_Raman = 10A_g + 5B_1g + 10B_2g + 5B_3g), 22 are IR active (Γ_IR = 4B_1u + 9B_2u + 9B_3u), and 5 are optically silent (Γ_silent = 5A_u). The Sb and S atoms have the site symmetry C_s that restricts their motions within the xz plane for the A_g, B_2g, B_1u and B_3u modes, and along the y axis for the B_1g, B_3g, A_u and B_2u modes. This symmetry is also consistent with the anisotropy reflected in the Born effective charges, where the only non-zero off-diagonal terms are due to the coupling in the x and z directions. The restriction means that motion in the xz plane is independent of the motion along the y direction, and hence the displacements of atoms along the y direction will not result in any polarization along x or z directions, and vice versa.

Currently there is no systematic study to assign the experimentally observed Raman modes in Sb₂S₃. Sereni et al.³⁵ recently performed a polarization-dependent Raman scattering study on single crystal samples of Sb₂S₃ in the 90° and 180° geometries. Contributions from B_1g spectrum were found in the experimental data for the A_g measurements. They also reported first-principles calculations of the zone-center phonons, but theoretical Raman scattering intensities were lacking. This makes the distinction between A_g and B_2g modes, and that between B_1g and B_3g modes very challenging, as the effect of microtwinning makes the a and c directions indistinguishable.³⁶ We present here the non-resonant Raman scattering coefficients which are computed from the second order derivative of the electronic density matrix with respect to a uniform electric field as implemented in Quantum ESPRESSO.³⁷ To assign the modes, the phonon frequencies, scattering coefficients and the space group symmetries are taken into consideration. We find that 13 phonon modes are in good agreement with the experimental values, as shown in Table 3. However, experimental studies on additional scattering geometries are needed to provide a comprehensive assignment.

Table 3 Zone-center phonon modes in Sb₂S₃ obtained using the supercell force-constant method, DFPT and Raman scattering spectroscopy. Only phonon modes with Raman data are shown. A complete list of Raman active modes are listed in the ESI

Raman mode	Intensity (Å^4 u^−1)	DFPT (cm^−1)	Supercell (cm^−1)	Expt.^21 (cm^−1)
B1g	267	47.7	47.7	43
B3g	1320	50.8	50.8	52
Ag	475	54.3	54.2	51
B3g	1140	69.1	69.0	60
Ag	1160	74.5	74.2	72
B2g	573	99.1	99.1	91
Ag	187	100.0	100.0	101
B2g	9.3	124.1	125.0	128
Ag	9970	196.7	197.4	192
B1g	1590	208.2	208.2	207
B3g	1390	231.4	231.4	239
Ag	4500	251.0	251.0	256
Ag	19 100	277.9	278.0	283

---

3.3 Phonon dispersions

Periodic images of displaced atoms can exert sizable effects during the calculation of forces, hence decreasing the accuracy of the IFC. This is not an issue for q-vectors that are commensurate with the supercell, since the phonon frequencies calculated at these points are exact using the supercell method. In order to reduce the effects of periodic images on non-commensurate q-points, a huge supercell has to be used. However, large primitive cells such as that of Sb₂S₃ place a computational constraint on the largest supercells that can be used in DFT calculations. As a result, we investigate the effect of supercell sizes on the phonon dispersions by using sizes of 1 × 4 × 1 containing 80 atoms, 2 × 2 × 2 containing 160 atoms and 2 × 4 × 2 containing 320 atoms. We plot the dispersion relations along the high symmetry directions Γ→ X → S → R → T → Z → Γ as shown in Fig. 2, where we have used the q-vector convention in ref. 38 and 39. Fig. 2(c) shows our most accurate results from the 2 × 4 × 2 supercell calculation.

Fig. 2 First-principles phonon dispersions along high symmetry directions for orthorhombic Sb₂S₃ calculated using the supercell force-constant method for 1 × 4 × 1, 2 × 2 × 2 and 2 × 4 × 2 supercells. The selected q-points are Γ = (0,0,0), X = (½,0,0), S = (½,½,0), R = (½,½,½), T = (0,½,½), and Z = (0,0,½). Blue circles are phonon frequencies calculated from DFPT at q-points commensurate with 2 × 4 × 2 supercell, while red circles are frequencies at q-points non-commensurate with any of the supercells. Imaginary frequencies (represented by negative frequencies) are present in the 1 × 4 × 1 supercell. Degeneracies are preserved in the 2 × 2 × 2 and 2 × 4 × 2 supercells but lost in the 1 × 4 × 1 supercell.

There are in general 60 phonon modes present in the dispersions, as evident in the X → Γ and Z → Γ directions. Double degeneracy occurs along the high symmetry lines X → S and R → T → Z, resulting in only 30 distinct phonon frequencies. Along S → R, the phonon frequencies are quadruply degenerate and only 15 frequencies are present. These degeneracies are preserved in the 2 × 2 × 2 and 2 × 4 × 2 results whereby the zone boundary points (X, S, R, T and Z) are commensurate with the supercell sizes, i.e.,

q·L_i = 2πn_i (8)

where n_i is an integer and L_i are the three supercell lattice vectors. It is interesting to note that commensuration at these points helps to ensure that the entire dispersion along the high symmetry directions also have the correct degeneracies. In contrast, degeneracies are lifted in smaller supercell sizes (Fig. 2(a)), as the zone boundary points are not commensurate with the supercell. Commensuration is thus extremely important in low symmetry crystals such as Sb₂S₃ to preserve the correct degeneracies as opposed to high symmetry crystals which may retain degeneracies even when the commensuration criteria are not met. The conclusion is not specific to Sb₂S₃ as we also observed the same behaviour for Bi₂S₃ as shown in the ESI.†

As shown in Fig. 2, the zone-center phonon frequencies of Sb₂S₃ differ in the X → Γ and Z → Γ directions due to macroscopic electric fields in the polar crystals (the so-called LO-TO splitting). D =(½,¼,0) and B = (0,¼,½) are only commensurate with 2 × 4 × 2 supercell, and therefore the agreement between DFPT and supercell results are excellent. There are, however, huge discrepancies in the smaller supercells, showing that the dispersions are not accurate. (½,⅜,0) and (0,⅜,½) are not commensurate with any of the supercell sizes, but the results from the supercell force-constant method and DFPT show a good agreement in Fig. 2(c).

Apart from affecting the degeneracies, small supercell sizes also result in the artificial softening of phonon branches, as can be seen in Fig. 2(a) and (b). The softening becomes so severe in the 1 × 4 × 1 supercell that imaginary frequencies are introduced in the acoustic phonon modes. Soft modes are also displayed in the dispersions of 1 × 1 × 1 and 1 × 2 × 1 supercells (see the ESI†).

One of the biggest challenges in phonon calculations is to determine whether the soft modes are artificial (an artefact of the numerical methods, etc.) or genuine that are associated with unstable lattice structures and phase transitions.^40,41 The occurrence of artificial phonon soft modes found in this work should serve as a caution that convergence with respect to supercell size must be checked carefully when using the supercell force-constant method. We conclude that softening is totally absent in the dispersions obtained with a 2 × 4 × 2 supercell for Sb₂S₃.

The vibrational density of states (vDOS) using the 2 × 4 × 2 supercell result is shown in Fig. 2(d). A dense 10 × 30 × 10 k-point mesh is used to sample the BZ and the effect of LO-TO splitting is included.

3.4 IFC analysis

We now propose a scheme to analyze the IFC in order to understand the origin of the phonon softening. For each pair of atoms, a 3 × 3 force constant matrix is obtained, with elements corresponding to movements of the each atom along the Cartesian directions. As a measure of the strength of this interaction, we use the trace of the IFC tensor, Tr(Φ_su) which has the advantage of being independent of the coordinate system used. The decrease of Tr(Φ_su) with increasing interatomic distance r_su suggests a suitable range for the interatomic forces.⁴²

Fig. 3 shows the decay in Tr(Φ_su) as a function of distances for different supercell sizes. The sth atom is the inequivalent S1 sulphur atom (Table 1), and u runs through all atoms within the supercell. Although not shown here, similar features are observed when the sth atom is replaced by other inequivalent atoms. Although the interactions are dominated by the on-site and the nearest neighbor terms, there are significant non-zero contributions at large distances which show the long ranged interactions in the crystals. The lattice dynamics of Sb₂S₃ thus cannot be approximated by a simple linear chain model considering only a few nearest neighbor interactions. To correctly describe the dispersions, a large supercell is needed to reduce the effect of periodic images and to capture the interactions at large distances.

Fig. 3 Tr(Φ_su) as a function of r_su = |r_s – r_u|, the distance between the sth and uth atom. The sth atom is the inequivalent S1 sulphur atom, and u runs through all the atoms within the supercell. The red triangles are data from the 1 × 4 × 1 supercell, blue diamonds from 2 × 2 × 2, and black squares for 2 × 4 × 2. The insets show significant interactions beyond the on-site and nearest neighbor atoms.

Due to the small sizes of 1 × 4 × 1 and 2 × 2 × 2 supercells, many values of IFC differ from those of 2 × 4 × 2 by a few orders of magnitude, as can be seen in Fig. 4. In addition, many interactions beyond ∼6 Å are not captured in the smaller supercells. We believe that these effects destabilize the crystal structure and result in the softening of the acoustic phonons. The softening is more pronounced in 1 × 4 × 1 supercells and less so in 2 × 2 × 2, due to the larger size of the supercell that reduces the effect of the periodic images and captures more long ranged interactions.

Fig. 4 The logarithm of the absolute value of Tr(Φ_su) as a function of the distance r_su, for the various supercell sizes of Sb₂S₃. The sth atom is the inequivalent S1 sulphur atom and u runs over all atoms within the supercell. The values of the IFC differ by up to a few orders of magnitude between the smaller supercells and 2 × 4 × 2 supercell.

4 Conclusion

To the best of our knowledge, the phonon dispersion of Sb₂S₃ is obtained for the first time through a systematic lattice dynamics study using the supercell force-constant method. The Born effective charges give rise to LO-TO splitting at the zone center and elucidate the covalent character of the bonds. Both the high frequency dielectric tensor and Born effective charges show considerable anisotropy of the crystals. The use of small supercell sizes results in the softening of the phonon modes that is inconsistent with experiments. We attribute this to the effect of the periodic images on the force constants, as well as the truncation of long ranged interactions. We found that a minimal 2 × 4 × 2 supercell (Pnma setting) is required for an accurate determination of the dispersion relations of Sb₂S₃. Our results suggest that when using the supercell force-constant method, the supercell size has to be tested with other parameters such as the kinetic energy cut-off, the Brillouin-zone sampling or the self-consistent convergence criteria especially when dealing with low symmetry systems such as Sb₂S₃.

Acknowledgements

The authors thank Peter Haynes of Imperial College London for useful discussions and pointing out ref. 13. Y.L and K.T.E.C acknowledge the financial support from the NSS programme, Singapore. T.C.S acknowledges the support by the following research grants: NTU start-up grant (M4080514); SPMS collaborative Research Award (M4080536); and the Singapore-Berkeley Research Initiative for Sustainable Energy (SinBeRISE) CREATE Programme. The authors gratefully acknowledge the use of resources at the A*STAR Computational Resource Centre, Singapore.

References

1. D. Aldakov, A. Lefrancois and P. Reiss, J. Mater. Chem. C, 2013, 1, 3756–3776.
2. J. C. Cardoso, C. A. Grimes, X. Feng, X. Zhang, S. Komarneni, M. V. B. Zanoni and N. Bao, Chem. Commun., 2012, 48, 2818–2820.
3. M. Schubert and W. Dollase, Opt. Lett., 2002, 27, 2073–2075.
4. E. Mãrquez, A. M. Bernal-Oliva, J. M. Gonzãlez-Leal, R. Prieto-Alcón and T. Wagner, J. Phys. D: Appl. Phys., 2006, 39, 1793.
5. Z. Mandouh and S. Salama, J. Mater. Sci., 1990, 25, 1715–1718.
6. Y. Itzhaik, O. Niitsoo, M. Page and G. Hodes, J. Phys. Chem. C, 2009, 113, 4254–4256.
7. S.-J. Moon, Y. Itzhaik, J.-H. Yum, S. M. Zakeeruddin, G. Hodes and M. Grätzel, J. Phys. Chem. Lett., 2010, 1, 1524–1527.
8. J. A. Chang, J. H. Rhee, S. H. Im, Y. H. Lee, H.-j. Kim, S. I. Seok, M. K. Nazeeruddin and M. Gratzel, Nano Lett., 2010, 10, 2609–2612.
9. J. Varghese, S. Barth, L. Keeney, R. W. Whatmore and J. D. Holmes, Nano Lett., 2012, 12, 868–872.
10. Z. Deng, M. Mansuripur and A. J. Muscat, Nano Lett., 2009, 9, 2015–2020.
11. U. Jeong, P. H. C. Camargo, Y. H. Lee and Y. Xia, J. Mater. Chem., 2006, 16, 3893–3897.
12. T. B. Nasr, H. Maghraoui-Meherzi, H. B. Abdallah and R. Bennaceur, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 406, 287–292.
13. R. Caracas and X. Gonze, Phys. Chem. Miner., 2005, 32, 295–300.
14. M. R. Filip, C. E. Patrick and F. Giustino, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 87, 205125.
15. C. E. Patrick and F. Giustino, Adv. Funct. Mater., 2011, 21, 4663–4667.
16. G. J. Ackland, M. C. Warren and S. J. Clark, J. Phys.: Condens. Matter, 1997, 9, 7861.
17. G. Kresse, J. Furthmüller and J. Hafner, Europhys. Lett., 1995, 32, 729.
18. W. Frank, C. Elsasser and M. Fahnle, Phys. Rev. Lett., 1995, 74, 1791–1794.
19. C. K. Gan, Y. P. Feng and D. J. Srolovitz, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 235214.
20. C. K. Gan, X. F. Fan and J.-L. Kuo, Comput. Mater. Sci., 2010, 49, S29.
21. Y. Zhao, K. T. E. Chua, C. K. Gan, J. Zhang, B. Peng, Z. Peng and Q. Xiong, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 205330.
22. P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari and R. M. Wentzcovitch, J. Phys.: Condens. Matter, 2009, 21, 395502.
23. H. J. Monkhorst and J. D. Pack, Phys. Rev. B: Solid State, 1976, 13, 5188–5192.
24. L. F. Lundegaard, R. Miletich, T. Balic-Zunic and E. Makovicky, Phys. Chem. Miner., 2003, 30, 463–468.
25. A. van de Walle and G. Ceder, Rev. Mod. Phys., 2002, 74, 11–45.
26. S. Baroni, S. de Gironcoli, A. D. Corso and P. Giannozzi, Rev. Mod. Phys., 2001, 73, 515.
27. X. Gonze, Phys. Rev. B: Condens. Matter Mater. Phys., 1997, 55, 10337.
28. W. Zhong, R. D. King-Smith and D. Vanderbilt, Phys. Rev. Lett., 1994, 72, 3618–3621.
29. W. Cochran and R. Cowley, J. Phys. Chem. Solids, 1962, 23, 447–450.
30. P. Giannozzi, S. de Gironcoli, P. Pavone and S. Baroni, Phys. Rev. B: Condens. Matter Mater. Phys., 1991, 43, 7231–7242.
31. Y. Wang, J. J. Wang, W. Y. Wang, Z. G. Mei, S. L. Shang, L. Q. Chen and Z. K. Liu, J. Phys.: Condens. Matter, 2010, 22, 202201.
32. M. Posternak, R. Resta and A. Baldereschi, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 8911–8914.
33. A. Roy, R. Prasad, S. Auluck and A. Garg, J. Phys.: Condens. Matter, 2010, 22, 165902.
34. J. Grigas, E. Talik and V. Lazauskas, Phase Transitions, 2002, 75, 323–337.
35. P. Sereni, M. Musso, P. Knoll, P. Blaha, K. Schwarz and G. Schmidt, AIP Conf. Proc., 2010, 1267, 1131–1132.
36. M. N. Iliev, M. V. Abrashev, H.-G. Lee, V. N. Popov, Y. Y. Sun, C. Thomsen, R. L. Meng and C. W. Chu, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 57, 2872–2877.
37. M. Lazzeri and F. Mauri, Phys. Rev. Lett., 2003, 90, 036401.
38. W. Setyawan and S. Curtarolo, Comput. Mater. Sci., 2010, 49, 299–312.
39. K. Rao, S. Chaplot, V. Padmanabhan and P. Vijayaraghavan, Pramana, 1982, 19, 593–632.
40. K. Parlinski, Z. Q. Li and Y. Kawazoe, Phys. Rev. Lett., 1997, 78, 4063–4066.
41. Y. Duan and K. Parlinski, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 104113.
42. A. J. E. Foreman and W. M. Lomer, Proc. Phys. Soc., London, Sect. B, 1957, 70, 1143.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c3cp53879f

---