# An environment-dependent interatomic potential for silicon carbide: calculation of bulk properties, high-pressure phases, point and extended defects, and amorphous structures

G Lucas, M. Bertolus, L Pizzagalli

► To cite this version:

G Lucas, M. Bertolus, L Pizzagalli. An environment-dependent interatomic potential for silicon carbide: calculation of bulk properties, high-pressure phases, point and extended defects, and amorphous structures. Journal of Physics: Condensed Matter, 2010, 22 (3), pp.035802. 10.1088/0953-8984/22/3/035802. cea-05503953

HAL Id: cea-05503953

https://cea.hal.science/cea-05503953v1

Submitted on 10 Feb 2026

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
An environment-dependent interatomic potential for silicon carbide: calculation of bulk properties, high-pressure phases, point and extended defects, and amorphous structures

This article has been downloaded from IOPscience. Please scroll down to see the full text article.

2010 J. Phys.: Condens. Matter 22 035802

(http://iopscience.iop.org/0953-8984/22/3/035802)

View the table of contents for this issue, or go to the journal homepage for more

Download details:

IP Address: 132.169.58.162

The article was downloaded on 13/07/2010 at 14:43

Please note that *terms and conditions apply*.
# An environment-dependent interatomic potential for silicon carbide: calculation of bulk properties, high-pressure phases, point and extended defects, and amorphous structures

G Lucas<fnref n="1" />, M Bertolus<fnref n="2" /> and L Pizzagalli<fnref n="3" />

E-mail: Laurent.Pizzagalli@univ-poitiers.fr

Received 15 October 2009

Published 21 December 2009

Online at stacks.iop.org/JPhysCM/22/035802

## Abstract

An interatomic potential has been developed to describe interactions in silicon, carbon and silicon carbide, based on the environment-dependent interatomic potential (EDIP) (Bazant *et al* 1997 *Phys. Rev. B* **56** 8542). The functional form of the original EDIP has been generalized and two sets of parameters have been proposed. Tests with these two potentials have been performed for many properties of SiC, including bulk properties, high-pressure phases, point and extended defects, and amorphous structures. One parameter set allows us to keep the original EDIP formulation for silicon, and is shown to be well suited for modelling irradiation-induced effects in silicon carbide, with a very good description of point defects and of the disordered phase. The other set, including a new parametrization for silicon, has been shown to be efficient for modelling point and extended defects, as well as high-pressure phases.

(Some figures in this article are in colour only in the electronic version)

## 1. Introduction

Silicon carbide (SiC) is a commonly studied semiconductor or because of its potential industrial and technological applications in electronics or in the nuclear environment [1, 2]. Besides its high temperature semiconductivity, it has excellent thermal and mechanical properties such as extreme hardness or large high temperature thermal conductivity. It also has excellent resistance to chemicals [3].

This material has been extensively studied both experimentally and theoretically. In the latter case, many calculations were based on atomistic simulations, such as molecular dynamics or Monte Carlo calculations, to study, for instance, dislocations, grain boundaries and interfaces, liquid, disordered or amorphous phases, diffusion and growth

processes [4–8]. With these methods, simulations normally involve a large number of atoms and the largest timescale possible. The reliability of these simulations depends for a large part on the quality of the empirical potential used to describe interatomic forces between atoms. In interatomic potentials, the quantum description of the binding is implicitly included in an analytic functional through empirical parameters. The quality of a potential relies obviously on both the functional form and the database used for fitting the parameters. In the case of covalent materials, it is difficult to properly model the oriented and localized nature of the bonds, and situations involving the formation and rupture of these bonds are usually not well described. In the past several potentials have been specifically proposed for semiconductor materials, the most famous ones being Stillinger-Weber [9],

<footnote n="1">¹ Sulzer Metco AG, Rigackerstrasse 16, CH-5610 Wohlen, Switzerland</footnote>

<footnote n="2">² CEA, DEN, DEC/SESC, Centre de Cadarache, F-13108 Saint-Paul-lez-Durance, France</footnote>

<footnote n="3">³ PhyMat, CNRS UMR 6630, Université de Poitiers, Boulevard Marie et Pierre Curie, SP2MI-Téléport 2, BP 30179, F-86962 Futuroscope Chasseneuil Cedex, France</footnote>
Tersoff [10] or Brenner [11] potentials, but none of them has appeared to be really superior to the others.

Concerning silicon carbide, the most widely used empirical potential is probably the bond order potential developed by Tersoff [10] consisting of two- and three-body terms. Several modified versions of Tersoff and Brenner potentials optimized for specific applications have also been published by Devanathan et al [12] and Gao et al [13]. Although relatively accurate in specific situations, their transferability to other configurations remains limited. Recently, an empirical potential based on the Brenner potential has been proposed by Erhart and Albe [14]. This potential appeared successfully to reproduce a wide range of silicon carbide properties, especially the ones involved in phase transitions, but still properties such as formation energies of point defects are far from being accurate in comparison to first-principles calculations.

In this study, we have developed a new empirical potential for silicon carbide based on EDIP (environment-dependent interatomic potential), initially provided for bulk silicon by Bazant et al [15–18] in order to provide a better description of point and extended defects. A modified version of EDIP for carbon-based materials has been proposed by Marks [19], but the modifications of the functional form for describing π-bonding make its use inadequate for our purpose. Here we generalize the original EDIP to silicon carbide and propose two different potentials: one keeping the original parametrization for bulk silicon and another one with a fully new set of parameters. After a description of EDIP functional form and its generalization to SiC, the results of calculations performed with the new potential are reported. Considered properties include bulk phases and elasticity constants, and point and extended defects for silicon, diamond and silicon carbide. Finally, amorphous silicon carbide is investigated and compared to previous results.

## 2. Functional form

### 2.1. Original form of EDIP

In the original EDIP for silicon, developed by Bazant et al [16–18], the total energy of an atomic configuration {$\vec{R}_i$} is expressed as a sum of individual atomic energies, which contain two-body and three-body terms:

$$E_i = \sum_{j \neq i} V_2(r_{ij}, Z_i) + \sum_{j \neq i} \sum_{k \neq i, k > i} V_3(\vec{r}_{ij}, \vec{r}_{ik}, Z_i). \quad (1)$$

$V_2(r_{ij}, Z_i)$ is an interaction between atoms i and j representing pairwise bonds and $V_3(\vec{r}_{ij}, \vec{r}_{ik}, Z_i)$ is an interaction between atoms i, j and k centred on atom i representing angular forces. These two types of interaction depend on the local environment of atom i thanks to an effective coordination number, defined by

$$Z_i = \sum_{m \neq i} f(r_{im}), \quad (2)$$

where $f(r_{im})$ is a cutoff function that measures the contribution of neighbour m to the coordination of atom i according to their separation $r_{im}$. This function is unitary for

$r < c$, decreases smoothly toward 0 when $r_{im}$ increases, and is equal to 0 for $r > a$:

$$f(r) = \begin{cases} 1, & r < c \\ \exp\left(\frac{\alpha}{1 - x^{-3}}\right), & c < r < a \\ 0, & r > a, \end{cases} \quad (3)$$

where $x = (r-c)/(a-c)$. A neighbour of atom i at a distance $r < c$ is considered a full neighbour, while further neighbours contribute partially to $Z$. The cutoff radii $c$ and $a$ are chosen to reproduce the coordination of the diamond lattice ($Z_i$).

The two-body term includes repulsive and attractive interactions:

$$V_2(r_{ij}, Z_i) = A \left[ \left( \frac{B}{r_{ij}} \right)^\rho - p(Z_i) \right] \exp \left( \frac{\sigma}{r_{ij} - a} \right), \quad (4)$$

and vanishes at the cutoff radius $a$ with all derivatives being continuous. It is equivalent to the two-body term in the Stillinger–Weber potential for small distortions of the diamond structure. In EDIP, the bond strength is modified by the local environment. This dependence has been motivated by theoretical calculations [15, 20], which have shown the weakening of the attractive interaction $p(Z)$ and the lengthening of bonds for increasing coordination. This feature can be accurately reproduced with a Gaussian function [16]:

$$p(Z) = \exp(-\beta Z^2). \quad (5)$$

The three-body term in (1) contains radial and angular factors:

$$V_3(\vec{r}_{ij}, \vec{r}_{ik}, Z_i) = g(r_{ij})g(r_{ik})h(l_{ijk}, Z_i), \quad (6)$$

where $l_{ijk} = \cos \theta_{ijk} = \vec{r}_{ij} \cdot \vec{r}_{ik} / r_{ij}r_{ik}$. The Stillinger–Weber form has been chosen for the radial function:

$$g(r) = \exp\left(\frac{\gamma}{r - a}\right), \quad (7)$$

which goes smoothly to zero at the cutoff distance $a$. The angular function $h(l, Z)$ strongly depends on the local coordination through two functions $\tau(Z)$ and $\omega(Z)$, which control the equilibrium angle and the force of the interaction, respectively. The general shape of this function has been postulated following theoretical considerations [16]:

$$h(l, Z) = H\left(\frac{l + \tau(Z)}{\omega(Z)}\right), \quad (8)$$

where $H(x)$ is a function satisfying the following constraints: $H(x) > 0$, $H(0) = 0$, $H'(0) = 0$ and $H''(0) > 0$. Bazant et al have finally proposed the angular function:

$$h(l, Z) = \lambda\left[(1 - \exp(-Q(Z)(l + \tau(Z))^2)) + \eta Q(Z)(l + \tau(Z))^2\right], \quad (9)$$

where $\omega(Z)^{-2} = Q(Z) = Q_0 \exp(-\mu Z)$ controls the force of the angular interactions. The angular term weakens with increasing coordination, allowing us to represent the transition between covalent and metallic bonds. The first contribution to
<table><caption>**Table 1.** Average of the cross-parameters. In the table, *g*, *a* and *wa* stand for geometric, arithmetic and weighted arithmetic means, respectively.</caption><thead><tr><th colspan="2">Two-body parameters</th><th colspan="2">Three-body parameters</th></tr></thead><tbody><tr><td>a:</td><td>p<sub>ij</sub> = (p<sub>i</sub> + p<sub>j</sub>)/2</td><td>ρ, β, α</td><td>wa: p<sub>ijk</sub> = (2p<sub>i</sub> + p<sub>j</sub> + p<sub>k</sub>)/4 λ, η, Q<sub>0</sub>, μ</td></tr><tr><td>g:</td><td>√p<sub>i</sub>p<sub>j</sub></td><td>A, B, a, c, σ, γ</td><td></td></tr></tbody></table>

this term, $H_1(x) \propto 1 - \exp(-x^2)$, is symmetric around the minimum and weak for small angles, and has been used by Mistriotis *et al* [21]. Here, it contains a dependence on the local environment. The second contribution, $H_2(x) \propto x^2$, is introduced to obtain a more asymmetrical shape, as suggested by tight-binding models and exact inversion of *ab initio* cohesive energy curves [15]. It yields stronger interactions for small angles.

The function $\tau(Z) = l_0(Z) = -\cos\theta_0(Z)$ controls the equilibrium angle $\theta_0(Z)$ of the three-body term according to the coordination. This specificity of EDIP allows us to model the hybridization of atoms in different environments. If a Si atom is threefold- or fourfold-coordinated, its bonds would prefer to be hybridized sp² or sp³ with equilibrium angles $\theta_0(Z = 3) = 120^\circ$ and $\theta_0(Z = 4) = 109.471^\circ$, respectively. Twofold and sixfold coordinations are also possible, but at the expense of a higher cost in energy. The function $\tau(Z)$ has been built to smoothly interpolate between these points ($Z = 2, 3, 4, 6$) with the following form:

$$ \tau(Z) = u_1 + u_2(u_3 \exp(-u_4 Z) - \exp(-2u_4 Z)) \quad (10) $$

with $u_1 = -0.165799$, $u_2 = 32.557$, $u_3 = 0.286198$ and $u_4 = 0.66$.

EDIP potential for bulk silicon has 13 adjustable parameters. These parameters have been fitted by Bazant *et al* using a database including experimental and *ab initio* results, such as bulk properties and formation energies of point and extended defects [18].

## 2.2. Generalization of EDIP to silicon carbide

Two potentials (A and B) have been developed, one keeping the original parametrization of EDIP for bulk silicon and the other allowing all parameters to vary during the fitting procedure.

### 2.2.1. Potential A.
The general principle to obtain potential A is to keep the original parametrization of the original EDIP potential for bulk silicon, as it has already been tested and applied in numerous studies. Only parameters for carbon have been optimized, and the cross-parameters for interactions involving Si and C atoms are simply arithmetic or geometric means of the other parameters (see table 1). The parameters $\lambda, \eta, Q_0$ and $\mu$ of the angular function $h(l, Z)$ are obtained by a considered arithmetic mean giving a higher weight to the central atom.

In the original EDIP, the external cutoff radius is about 3.12 Å, which corresponds roughly to the second-neighbour distance (3.08 Å) in SiC. This leads to an unphysical behaviour (discontinuity) in the cohesive energy curve close to the equilibrium distance of SiC. Indeed each Si atom feels the

interaction from its first C neighbour but also from its second Si neighbour, leading to an effective coordination larger than the expected fourfold coordination. To solve this issue, a corrective term has been introduced in the Si–Si interactions when the central atom is surrounded by Si and C atoms. This function reduces the external cutoff radius by a distance $\delta$, when the partial coordination of the central atom $Z_{C-Si}$ is not zero:

$$ \Delta(r) = \begin{cases} 1, & r < c \\ \exp\left(\frac{\alpha}{1-x^{-3}}\right), & c < r < a - \delta \\ 0, & r > a - \delta \end{cases} \quad (11) $$

where $x = (r - c)/(a - \delta - c)$. This correction is applied on the coordination function, the pairwise function $V_2$ and the angular function $V_3$:

$$ f(r) = \exp\left(\frac{\alpha}{1-x^{-3}}\right) \Delta(r) \quad (12) $$

$$ V_2(r_{ij}, Z_i) = A \left[ \left( \frac{B}{r_{ij}} \right)^\rho - p(Z_i) \right] \exp \left( \frac{\sigma}{r_{ij} - a} \right) \Delta(r_{ij}), \quad (13) $$

$$ V_3(\vec{r}_{ij}, \vec{r}_{ik}, Z_i) = g(r_{ij})\Delta(r_{ij})g(r_{ik})\Delta(r_{ik})h(l_{ijk}, Z_i), \quad (14) $$

where $\delta$ has been fixed to 0.3 Å. This value allows us to avoid the above-described discontinuity around the equilibrium distance while keeping the original behaviour of the Si–Si interactions.

### 2.2.2. Potential B.
In potential B, almost all the parameters have been optimized, allowing a greater flexibility during the fitting procedure, and only the parameters of the angular function $h(l, Z)$ are obtained from a weighted arithmetic average of the pure element parameters. As a consequence, a new parametrization for silicon has been obtained. In this case, since the new external cutoff radius for Si–Si interaction is below 2.95 Å, there is no need to introduce a corrective function.

## 2.3. Fitting procedure

Our EDIP potentials have been parametrized in two successive steps. First, a simulated annealing algorithm [22, 23] has been used. It is a stochastic method allowing us to explore a multivariable space defined by an objective function, which in our case includes a weighted least-squares term and a penalty term. The algorithm allows a local minimum to escape when optimizing a set of parameters. In a second step, the simplex algorithm of Nelder and Mead [24], suitable for nonlinear optimization, is used to refine the solution obtained.
<table><caption>**Table 2.** Parameters of the potentials. In potential A, Si parameters remain the same as in the original EDIP potential for bulk silicon [18]. In the table, *g*, *a* and *wa* stand for geometric, arithmetic and weighted arithmetic means, respectively. A corrective function is introduced in potential A ($δ = 0.3 Å$), see the text for details.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Potential A</th><th colspan="3">Potential B</th></tr><tr><th>Si</th><th>SiC</th><th>C</th><th>Si</th><th>SiC</th><th>C</th></tr></thead><tbody><tr><td colspan="7"><i>f(r)</i></td></tr><tr><td>a (Å)</td><td>3.121 382</td><td>g</td><td>2.291 805</td><td>2.941 586</td><td>2.534 972</td><td>2.212 263</td></tr><tr><td>c (Å)</td><td>2.560 910</td><td>g</td><td>1.569 857</td><td>2.540 193</td><td>1.973 974</td><td>1.741 598</td></tr><tr><td>α</td><td>3.108 385</td><td>a</td><td>1.126 809</td><td>3.066 580</td><td>2.507 738</td><td>1.962 090</td></tr><tr><td colspan="7"><i>V<sub>2</sub></i></td></tr><tr><td>A (eV)</td><td>7.982 173</td><td>g</td><td>11.677 335</td><td>5.488 043</td><td>7.535 967</td><td>10.222 599</td></tr><tr><td>B (Å)</td><td>1.507 546</td><td>g</td><td>0.961 202</td><td>1.446 435</td><td>1.177 019</td><td>0.959 814</td></tr><tr><td>ρ</td><td>1.208 520</td><td>a</td><td>2.877 953</td><td>1.343 679</td><td>2.061 835</td><td>2.827 634</td></tr><tr><td>β</td><td>0.007 097</td><td>a</td><td>0.025 596</td><td>0.008 593</td><td>0.015 347</td><td>0.025 661</td></tr><tr><td>σ (Å)</td><td>0.577 411</td><td>g</td><td>0.710 345</td><td>0.298 443</td><td>0.423 863</td><td>0.536 561</td></tr><tr><td colspan="7"><i>V<sub>3</sub></i></td></tr><tr><td>γ (Å)</td><td>1.124 794</td><td>g</td><td>1.225 825</td><td>1.135 256</td><td>1.191 567</td><td>1.084 183</td></tr><tr><td>λ (eV)</td><td>1.453 311</td><td>wa</td><td>1.780 018</td><td>2.417 497</td><td>wa</td><td>3.633 621</td></tr><tr><td>η</td><td>0.252 324</td><td>wa</td><td>0.811 587</td><td>0.589 390</td><td>wa</td><td>0.275 605</td></tr><tr><td>Q<sub>0</sub></td><td>312.134 135</td><td>wa</td><td>417.547 683</td><td>208.924 548</td><td>wa</td><td>289.305 617</td></tr><tr><td>μ</td><td>0.696 633</td><td>wa</td><td>0.612 215</td><td>0.629 131</td><td>wa</td><td>0.594236</td></tr></tbody></table>

The database for the fitting procedure includes experimental data such as cohesive energies, lattice parameters and elastic constants of Si, C and SiC, which ensure that the potentials will be able to accurately reproduce the structural and elastic properties of these materials. Moreover, we have biased the parametrization in order to improve the description of defects in SiC: the formation energies of several defects, i.e. $C_{TSi}$, $CC_{100}$ and $SiSi_{110}$ interstitials, have been added into the database. Note that we have also attempted to include the $Si_{TC}$ interstitial, but it did not lead to a better description of this specific defect configuration. The best sets of parameters obtained for potentials A and B are given in table 2.

## 2.4. Bulk phases

### 2.4.1. Bulk silicon.
Structural and elastic properties of bulk silicon, computed with the molecular dynamics code XMD [25], are reported in table 3. For the diamond phase, results obtained with the parameter set A are obviously identical to the original EDIP [18]. Overall, properties of the diamond phase are well described with all interatomic potentials. This is not so surprising since potential parameters are usually fitted on these properties. We have also investigated high-pressure and exotic silicon phases, for which a higher coordination leads to an increasing metallic behaviour [26]. DFT calculations suggested that the differences in the cohesive energies of the various phases are rather small and are therefore very difficult to reproduce with classical potentials. The new silicon parametrization of EDIP B gives a significant improvement over EDIP A for these phases. For low energy bulk structures, EDIP B provides a close agreement with DFT data. For the high-pressure BC8 phase, $a_0$ and $ΔE$ are found to be 6.67 Å and 0.130 eV with DFT [27], and 6.698 Å and 0.168 eV with EDIP B. For the β-tin phase,

both EDIP yield overestimated cohesive energies. Since the EDIP lattice parameters are in agreement with the DFT value, it is likely that the potential overestimates the transition pressure between diamond and β-tin phases. Finally, for other high energy phases, we found that structural parameters are usually in agreement with reference data. However, cohesive energies are often overestimated, especially with the original parameter set A. The new parametrization B greatly improves the description.

### 2.4.2. Bulk carbon.
Carbon can be found naturally with two different structures: diamond and graphite, which have respectively sp³ and sp² bonds. At room temperature, graphite is thermodynamically more stable than diamond by less than 1 meV. Graphite sheets are held together by long-range van der Waals interactions, a feature difficult to reproduce with short-range classical potentials. The results reported here have then been obtained by fixing the separation between the sheets to the experimental value. As shown in table 4, EDIP potentials reproduce reasonably well the cohesive energies of both structures, the parameter set B yielding the graphitic structure as slightly more stable than diamond. However, like other potentials, the lattice parameter of graphite is somewhat overestimated by EDIP potentials. The elastic constants of the diamond phase are in good agreement with DFT calculations, except for $C_{44}$, which is underestimated. Overall, there are very few differences between EDIP and previous potentials regarding graphite and diamond properties, except for the Brenner potential which poorly reproduces the elastic constants of the diamond structure.

The EDIP potentials have also been used for modelling several high energy phases, with coordination numbers up to $Z = 8$. These results, shown in table 4, are compared to those found by Furthmüller *et al* within DFT [31]. Globally
<table><caption>**Table 3.** Structure and elastic properties of bulk silicon phases obtained with EDIP potentials, and compared to experiments [28], DFT-LDA calculations [18, 27, 29, 30] and empirical calculations using Stillinger-Weber (SW) [9], Tersoff (T89) [10] and Erhart-Albe (EA05) [14] potentials. Energies are given relative to the diamond phase.</caption><thead><tr><th></th><th>Exp.</th><th>DFT-LDA</th><th>EDIP A</th><th>EDIP B</th><th>SW</th><th>T89</th><th>EA05</th></tr></thead><tbody><tr><td colspan="8" style="text-align:center; font-weight:bold;">Diamond (Z = 4)</td></tr><tr><td>a (Å)</td><td>5.43</td><td></td><td>5.430</td><td>5.425</td><td>5.431</td><td>5.432</td><td>5.429</td></tr><tr><td>E<sub>c</sub> (eV)</td><td>-4.65</td><td></td><td>-4.650</td><td>-4.647</td><td>-4.63</td><td>-4.63</td><td>-4.63</td></tr><tr><td>B (GPa)</td><td>99</td><td></td><td>99</td><td>101</td><td>108</td><td>98</td><td>99</td></tr><tr><td>C<sub>11</sub> (GPa)</td><td>167</td><td></td><td>175</td><td>155</td><td>162</td><td>142</td><td>167</td></tr><tr><td>C<sub>12</sub> (GPa)</td><td>65</td><td></td><td>65</td><td>75</td><td>82</td><td>75</td><td>67</td></tr><tr><td>C<sub>44</sub> (GPa)</td><td>81</td><td></td><td>71</td><td>59</td><td>60</td><td>69</td><td>75</td></tr><tr><td>C<sub>44</sub><sup>0</sup> (GPa)</td><td>112</td><td></td><td>112</td><td>111</td><td>117</td><td>119</td><td>105</td></tr><tr><td colspan="8" style="text-align:center; font-weight:bold;">BC8 (Z = 4)</td></tr><tr><td>a (Å)</td><td>6.67</td><td></td><td>6.689</td><td>6.698</td><td>6.591</td><td>6.644</td><td>6.625</td></tr><tr><td>ΔE<sub>c</sub> (eV)</td><td>0.130</td><td></td><td>0.224</td><td>0.168</td><td>0.201</td><td>0.245</td><td>0.527</td></tr><tr><td colspan="8" style="text-align:center; font-weight:bold;">β-tin (Z = 6)</td></tr><tr><td>a (Å)</td><td>4.730</td><td></td><td>4.854</td><td>4.832</td><td>4.969</td><td>4.905</td><td>4.856</td></tr><tr><td>c/a</td><td>0.552</td><td></td><td>0.519</td><td>0.519</td><td>0.561</td><td>0.524</td><td>0.527</td></tr><tr><td>ΔE<sub>c</sub> (eV)</td><td>0.210</td><td></td><td>0.603</td><td>0.576</td><td>0.213</td><td>0.327</td><td>0.412</td></tr><tr><td colspan="8" style="text-align:center; font-weight:bold;">SH (Z = 6)</td></tr><tr><td>a (Å)</td><td>2.639</td><td></td><td>2.875</td><td>2.785</td><td>2.833</td><td>2.699</td><td>2.659</td></tr><tr><td>c/a</td><td>0.94</td><td></td><td>0.857</td><td>0.883</td><td>0.918</td><td>0.967</td><td>0.965</td></tr><tr><td>ΔE<sub>c</sub> (eV)</td><td>0.293</td><td></td><td>1.216</td><td>0.570</td><td>0.403</td><td>0.469</td><td>0.476</td></tr><tr><td colspan="8" style="text-align:center; font-weight:bold;">SC (Z = 6)</td></tr><tr><td>a (Å)</td><td>2.528</td><td></td><td>2.497</td><td>2.493</td><td>2.612</td><td>2.544</td><td>2.525</td></tr><tr><td>ΔE<sub>c</sub> (eV)</td><td>0.348</td><td></td><td>0.521</td><td>0.544</td><td>0.293</td><td>0.318</td><td>0.397</td></tr><tr><td colspan="8" style="text-align:center; font-weight:bold;">BCC (Z = 8)</td></tr><tr><td>a (Å)</td><td>3.088</td><td></td><td>3.236</td><td>3.166</td><td>3.245</td><td>3.084</td><td>3.043</td></tr><tr><td>ΔE<sub>c</sub> (eV)</td><td>0.525</td><td></td><td>1.601</td><td>0.889</td><td>0.300</td><td>0.432</td><td>0.503</td></tr><tr><td colspan="8" style="text-align:center; font-weight:bold;">FCC (Z = 8)</td></tr><tr><td>a (Å)</td><td>3.885</td><td></td><td>4.078</td><td>3.941</td><td>4.147</td><td>3.897</td><td>3.940</td></tr><tr><td>ΔE<sub>c</sub> (eV)</td><td>0.566</td><td></td><td>1.773</td><td>0.688</td><td>0.423</td><td>0.761</td><td>0.587</td></tr></tbody></table>

the EDIP potentials yield a good structural description with lattice parameters close to the ones calculated in DFT, potential B being the best. Regarding stability, the EDIP potentials tend to overestimate the cohesion energies for structures with high coordinations (Z ≥ 6). Nevertheless, the stability order of these high energy phases remains coherent with DFT calculations. The best agreement between DFT and potential calculations of cohesive energies is obtained with the Erhart-Albe potential, likely due to the inclusion of these energies in the database used to fit the potential.

### 2.4.3. Bulk silicon carbide.
The properties of several common SiC polytypes (3C, 4H and 6H), as well as high energy phases (B1 and B2), have been calculated with the developed potentials. Table 5 shows experimental data [32] as well as theoretical results obtained with DFT calculations [34–36], EDIP A and B potentials, and other empirical potentials for SiC (Erhart-Albe [14], Tersoff [10] and Gao-Weber [13]). The lattice constant and cohesive energy of the 3C cubic phase are well described by all potentials. Elastic constants computed with both EDIP

potentials are in good agreement with experimental data, especially with EDIP B. Nevertheless, the C<sub>44</sub> constant is underestimated. We observed a similar behaviour for the diamond phase, which suggests it may be due to the EDIP functional, and not the fitting procedure. Overall only the Erhart-Albe potential gives a better description than EDIP B. In contrast, the Gao-Weber potential leads to a poor description of the elastic constants for the cubic phase.

The cohesive energies of the 3C, 4H and 6H polytypes differ by a few meV [34–36], as they share the same local tetrahedral environment. As EDIP potentials take into account only first-neighbour interactions, they do not allow us to distinguish between polytypes. Here, cohesive energy and bond lengths remain strictly identical for all three polytypes in potential calculations. We have computed the elastic constants of the 4H polytype, which are found in very good agreement with experimental values. Experimentally, it is established that the cubic 3C phase transforms to a B1 structure (NaCl) above 100 GPa [37], even if the transition path is still a matter of discussion [38]. Results obtained with EDIP potentials for B1 (NaCl) and B2 (CsCl) structures are compared to DFT
<table><caption>**Table 4.** Structure and elastic properties of bulk carbon. Results obtained using EDIP potentials are compared to experiments [3], DFT-LDA calculations [31] and empirical calculations using Erhart and Albe (EA05), Brenner (B90) and Tersoff (T89) potentials [14]. Energies are given relative to the graphite phase.</caption><thead><tr><th></th><th>Exp.</th><th>DFT-LDA</th><th>EDIP A</th><th>EDIP B</th><th>EA05</th><th>B90</th><th>T89</th></tr></thead><tbody><tr><td colspan="8">Graphite (Z = 3)</td></tr><tr><td>a (Å)</td><td>2.46</td><td>2.440</td><td>2.570</td><td>2.558</td><td>2.555</td><td>2.513</td><td>2.530</td></tr><tr><td>E<sub>c</sub> (eV)</td><td>-7.374</td><td>-9.030</td><td>-7.382</td><td>-7.371</td><td>-7.374</td><td>-7.376</td><td>-7.396</td></tr><tr><td colspan="8">Diamond (Z = 4)</td></tr><tr><td>a (Å)</td><td>3.567</td><td>3.528</td><td>3.535</td><td>3.561</td><td>3.566</td><td>3.558</td><td>3.566</td></tr><tr><td>ΔE<sub>c</sub> (eV)</td><td></td><td>0.0004</td><td>-0.0001</td><td>0.0008</td><td>0.0009</td><td>0.0524</td><td>0.0250</td></tr><tr><td>B (GPa)</td><td>444</td><td>460</td><td>438</td><td>441</td><td>445</td><td>484</td><td>425</td></tr><tr><td>C<sub>11</sub> (GPa)</td><td>1081</td><td>1100</td><td>1057</td><td>1079</td><td>1082</td><td>621</td><td>1067</td></tr><tr><td>C<sub>12</sub> (GPa)</td><td>125</td><td>143</td><td>127</td><td>121</td><td>127</td><td>441</td><td>104</td></tr><tr><td>C<sub>44</sub> (GPa)</td><td>579</td><td>587</td><td>480</td><td>490</td><td>635</td><td>393</td><td>636</td></tr><tr><td>C<sub>44</sub><sup>0</sup> (GPa)</td><td></td><td></td><td>543</td><td>551</td><td>673</td><td>642</td><td>671</td></tr><tr><td colspan="8">BC8 (Z = 4)</td></tr><tr><td>a (Å)</td><td></td><td>4.419</td><td>4.394</td><td>4.435</td><td>4.429</td><td>4.351</td><td>4.437</td></tr><tr><td>ΔE<sub>c</sub> (eV)</td><td></td><td>0.693</td><td>0.495</td><td>0.503</td><td>0.722</td><td>0.446</td><td>0.775</td></tr><tr><td colspan="8">SC (Z = 6)</td></tr><tr><td>a (Å)</td><td></td><td>1.744</td><td>1.825</td><td>1.733</td><td>1.783</td><td>1.744</td><td>1.802</td></tr><tr><td>ΔE<sub>c</sub> (eV)</td><td></td><td>2.637</td><td>3.432</td><td>3.577</td><td>3.297</td><td>2.133</td><td>2.974</td></tr><tr><td colspan="8">BCC (Z = 8)</td></tr><tr><td>a (Å)</td><td></td><td>2.326</td><td>2.297</td><td>2.308</td><td>2.160</td><td>2.093</td><td>2.152</td></tr><tr><td>ΔE<sub>c</sub> (eV)</td><td></td><td>4.355</td><td>5.059</td><td>5.105</td><td>3.964</td><td>3.037</td><td>3.771</td></tr><tr><td colspan="8">FCC (Z = 8)</td></tr><tr><td>a (Å)</td><td></td><td>3.021</td><td>2.953</td><td>2.905</td><td>2.859</td><td>2.863</td><td>2.728</td></tr><tr><td>ΔE<sub>c</sub> (eV)</td><td></td><td>4.652</td><td>6.242</td><td>5.925</td><td>4.483</td><td>3.713</td><td>4.441</td></tr></tbody></table>

investigations done by Karch *et al* [36]. We found that both lattice parameters and cohesive energies are overestimated, as with other potentials. With EDIP and Tersoff potentials, the equilibrium volume of the B1 phase is even larger than the cubic structure, in clear contradiction with experiments. It seems that this inability to describe high-pressure phases of silicon carbide correctly is a feature shared by all classical potentials.

## 2.5. Point defects

Initially, EDIP was developed in order to improve the description of point and extended defects in silicon. It was therefore important to preserve this feature as much as possible in the new potential. Particular care has thus been taken during the fitting procedure, and the formation energies of three lowest energy point defects in silicon carbide have been included in the fitting database.

### 2.5.1. Silicon.
The properties of several point defects in bulk silicon, i.e. the vacancy and two interstitial configurations ((110) dumbbell and hexagonal), have been investigated. Their formation energies, as well as the activation energy associated with the concerted exchange migration mechanism [39], are reported in table 6. Considering density functional theory calculations [40, 41] as the reference, a good

agreement is obtained for EDIP A and for the Erhart–Albe potential. While the other potentials yield reasonable vacancy formation energies, the interstitial formation energies are either overestimated (SW, T89) or underestimated (EDIP B).

### 2.5.2. Carbon (diamond).
The formation energies of point defects in diamond obtained using various potentials have been compared to DFT calculations of the vacancy [42] and of the interstitials [43] (see table 7). Concerning the vacancy, EDIP A and B potentials yield formation energies of 1.52 eV and 1.58 eV, respectively. These values are considerably lower than the *ab initio* value, but it has not been possible to increase the vacancy formation energy during the fit without deteriorating the accuracy of the potential in other aspects. Other potentials yield closer values, although still relatively different from the DFT energy. It is likely that the formation of a vacancy induces strong changes in the electronic structure, which are difficult to reproduce with an empirical potential. Considering now interstitials, the most stable configuration according to Breuer *et al* is the split interstitial oriented along (100) (denoted $I_S$ in this paper). A configuration with an interstitial centred on a bond was calculated to be 4.05 eV higher in energy, whereas an interstitial initially in a tetrahedral site relaxed into the $I_S$ configuration. Our calculations indicate that $I_S$ was the most stable state with both EDIP potentials, and that other
<table><caption>**Table 5.** Structure and elastic properties of bulk silicon carbide. Results obtained using EDIP potentials are compared to experiments [32, 33], DFT-LDA calculations [34–36] and empirical calculations using Erhart and Albe (EA05), Tersoff (T89) and Gao and Weber (GW02) potentials [14]. Energies are given relative to the 3C-SiC phase.</caption><thead><tr><th></th><th>Exp.</th><th>DFT-LDA</th><th>EDIP A</th><th>EDIP B</th><th>EA05</th><th>T89</th><th>GW02</th></tr></thead><tbody><tr><td colspan="8">3C-SiC</td></tr><tr><td>a (Å)</td><td>4.3596</td><td>4.344</td><td>4.411</td><td>4.364</td><td>4.359</td><td>4.321</td><td>4.360</td></tr><tr><td>E<sub>c</sub> (eV)</td><td>-6.340</td><td>-7.415</td><td>-6.359</td><td>-6.338</td><td>-6.340</td><td>-6.165</td><td>-6.412</td></tr><tr><td>B (GPa)</td><td>225</td><td>222</td><td>224</td><td>226</td><td>224</td><td>224</td><td>235</td></tr><tr><td>C<sub>11</sub> (GPa)</td><td>390</td><td>390</td><td>437</td><td>394</td><td>382</td><td>437</td><td>254</td></tr><tr><td>C<sub>12</sub> (GPa)</td><td>142</td><td>134</td><td>117</td><td>142</td><td>145</td><td>118</td><td>225</td></tr><tr><td>C<sub>44</sub> (GPa)</td><td>256</td><td>253</td><td>195</td><td>168</td><td>240</td><td>311</td><td>66</td></tr><tr><td>C<sub>44</sub><sup>0</sup> (GPa)</td><td></td><td>273</td><td>260</td><td>255</td><td>305</td><td></td><td></td></tr><tr><td colspan="8">4H-SiC</td></tr><tr><td>a (Å)</td><td>3.073</td><td></td><td>3.119</td><td>3.085</td><td></td><td></td><td></td></tr><tr><td>c (Å)</td><td>10.053</td><td></td><td>10.186</td><td>10.074</td><td></td><td></td><td></td></tr><tr><td>B (GPa)</td><td>220</td><td></td><td>223</td><td>227</td><td></td><td></td><td></td></tr><tr><td>C<sub>11</sub> (GPa)</td><td>501</td><td></td><td>536</td><td>521</td><td></td><td></td><td></td></tr><tr><td>C<sub>12</sub> (GPa)</td><td>111</td><td></td><td>87</td><td>105</td><td></td><td></td><td></td></tr><tr><td>C<sub>13</sub> (GPa)</td><td>52</td><td></td><td>54</td><td>62</td><td></td><td></td><td></td></tr><tr><td>C<sub>33</sub> (GPa)</td><td>553</td><td></td><td>568</td><td>567</td><td></td><td></td><td></td></tr><tr><td>C<sub>44</sub> (GPa)</td><td>163</td><td></td><td>194</td><td>170</td><td></td><td></td><td></td></tr><tr><td>ΔE<sub>c</sub> (eV)</td><td></td><td>-0.0024</td><td>0.000</td><td>0.000</td><td></td><td></td><td></td></tr><tr><td colspan="8">6H-SiC</td></tr><tr><td>a (Å)</td><td>3.081</td><td></td><td>3.119</td><td>3.085</td><td></td><td></td><td></td></tr><tr><td>c (Å)</td><td>15.117</td><td></td><td>15.279</td><td>15.112</td><td></td><td></td><td></td></tr><tr><td>ΔE<sub>c</sub> (eV)</td><td></td><td>-0.0017</td><td>0.000</td><td>0.000</td><td></td><td></td><td></td></tr><tr><td colspan="8">B1 (NaCl)</td></tr><tr><td>a (Å)</td><td>4.046</td><td></td><td>4.455</td><td>4.415</td><td>4.244</td><td>4.329</td><td></td></tr><tr><td>ΔE<sub>c</sub> (eV)</td><td></td><td>0.7</td><td>2.06</td><td>1.26</td><td>1.92</td><td>1.49</td><td></td></tr><tr><td colspan="8">B2 (CsCl)</td></tr><tr><td>a (Å)</td><td>2.631</td><td></td><td>2.761</td><td>2.729</td><td>2.668</td><td>2.640</td><td></td></tr><tr><td>ΔE<sub>c</sub> (eV)</td><td></td><td>2.1</td><td>2.73</td><td>2.61</td><td>3.04</td><td>2.49</td><td></td></tr></tbody></table>

<table><caption>**Table 6.** Formation energies (in eV) of point defects and concerted exchange (CE) migration energy in bulk silicon computed using EDIP potentials. Results are compared to DFT calculations [40, 41] and empirical calculations using Stillinger–Weber (SW), Tersoff (T89) and Erhart–Albe (EA05) [14, 18] potentials.</caption><thead><tr><th></th><th>DFT</th><th>EDIP (Å)</th><th>EDIP B</th><th>SW</th><th>T89</th><th>EA05</th></tr></thead><tbody><tr><td>V</td><td>3.17</td><td>3.22</td><td>2.98</td><td>2.82</td><td>3.70</td><td>3.2</td></tr><tr><td>I<sub>(110)</sub></td><td>3.31</td><td>3.35</td><td>1.64</td><td>4.68</td><td>4.68</td><td>3.0</td></tr><tr><td>I<sub>H</sub></td><td>3.31</td><td>4.16</td><td>2.06</td><td>6.95</td><td>4.61</td><td>4.0</td></tr><tr><td>CE</td><td>4.45–4.80</td><td>4.82</td><td>4.40</td><td>4.64</td><td></td><td></td></tr></tbody></table>

configurations convert to $I_S$ after relaxation, or end up with a very high energy.

### 2.5.3. Silicon carbide.
One of the main objectives of this work was an improvement of the point defects description in silicon carbide compared to existing classical potentials. Here, we show the results obtained for the cubic phase. Due to the similarity between local environments of cubic and hexagonal phases, the general conclusions drawn here should remain valid for the other common polytypes such as 4H and 6H phases. The formation energies for vacancies, antisites

<table><caption>**Table 7.** Formation energies (in eV) of several point defects in diamond calculated with EDIP potentials. Results are compared to DFT calculations [42, 43] and empirical calculations using Erhart and Albe (EA05), Brenner (B90) and Tersoff (T89) potentials [14]. Interstitial formation energies are given relative to the split interstitial $I_S$.</caption><thead><tr><th></th><th>DFT</th><th>EDIP A</th><th>EDIP B</th><th>EA05</th><th>B90</th><th>T89</th></tr></thead><tbody><tr><td>V</td><td>6.98</td><td>1.52</td><td>1.58</td><td>5.24</td><td>8.84</td><td>3.42</td></tr><tr><td>I<sub>S</sub></td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>I<sub>B</sub></td><td>4.05</td><td>I<sub>S</sub></td><td>I<sub>S</sub></td><td>5.85</td><td>3.79</td><td>4.53</td></tr><tr><td>I<sub>T</sub></td><td>I<sub>S</sub></td><td>I<sub>S</sub></td><td>13.51</td><td>13.69</td><td>3.72</td><td>9.49</td></tr></tbody></table>

and many configurations of carbon and silicon interstitials, calculated with DFT [44] (for defects in their neutral charge state), and the various potentials are listed in table 8.

Considering first the vacancies, the formation energies computed using both EDIP potentials are lower than DFT values, but differences are smaller than in the case of carbon. Similar trends are obtained for Erhart–Albe and Gao–Weber potentials. With the Tersoff potential, the carbon vacancy has a higher formation energy than the silicon vacancy, in disagreement with DFT calculations. In the case of antisites, a not perfect but nonetheless satisfactory behaviour is noted
<table><caption>**Table 8.** Formation energies (in eV) of point defects in 3C-SiC obtained with our EDIP potentials. Results are compared to DFT calculations [44] and empirical calculations using Erhart and Albe (EA05), Tersoff (T89) and Gao and Weber potentials [14].</caption><thead><tr><th></th><th>DFT</th><th>EDIP A</th><th>EDIP B</th><th>EA05</th><th>T89</th><th>GW02</th></tr></thead><tbody><tr><td colspan="7">Vacancies</td></tr><tr><td>V<sub>C</sub></td><td>3.63</td><td>1.25</td><td>1.45</td><td>1.90</td><td>3.88</td><td>1.39</td></tr><tr><td>V<sub>Si</sub></td><td>7.48</td><td>3.98</td><td>4.18</td><td>4.55</td><td>3.29</td><td>4.67</td></tr><tr><td colspan="7">Antisites</td></tr><tr><td>C<sub>Si</sub></td><td>3.48</td><td>3.01</td><td>2.40</td><td>2.42</td><td>2.20</td><td>4.43</td></tr><tr><td>Si<sub>C</sub></td><td>4.02</td><td>2.83</td><td>2.74</td><td>2.48</td><td>4.50</td><td>5.05</td></tr><tr><td colspan="7">Carbon interstitials</td></tr><tr><td>C<sub>TC</sub></td><td>CC<sub>(100)</sub></td><td>12.15</td><td>CC<sub>(100)</sub></td><td>12.63</td><td>7.21</td><td>6.02</td></tr><tr><td>C<sub>TSi</sub></td><td>CC<sub>(100)</sub></td><td>8.57</td><td>6.40</td><td>9.38</td><td>4.40</td><td>5.69</td></tr><tr><td>CC<sub>(100)</sub></td><td>6.31</td><td>5.42</td><td>4.82</td><td>4.78</td><td>6.50</td><td>4.41</td></tr><tr><td>CC<sub>(110)</sub></td><td>6.65</td><td>CSi<sub>(100)</sub></td><td>5.13</td><td>9.89</td><td>5.68</td><td>4.67</td></tr><tr><td>CSi<sub>(100)</sub></td><td>6.94</td><td>6.05</td><td>4.67</td><td>8.31</td><td>7.69</td><td>4.80</td></tr><tr><td>CSi<sub>(110)</sub></td><td>CC<sub>(100)</sub></td><td>CC<sub>(100)</sub></td><td>CC<sub>(100)</sub></td><td>4.81</td><td>CC<sub>(100)</sub></td><td>5.32</td></tr><tr><td colspan="7">Silicon interstitials</td></tr><tr><td>Si<sub>TC</sub></td><td>7.04</td><td>SiC<sub>(100)</sub></td><td>SiC<sub>(110)</sub></td><td>17.55</td><td>17.67</td><td>2.60</td></tr><tr><td>Si<sub>TSi</sub></td><td>9.23</td><td>SiC<sub>(100)</sub></td><td>SiC<sub>(110)</sub></td><td>17.30</td><td>15.89</td><td>5.40</td></tr><tr><td>SiSi<sub>(100)</sub></td><td>9.32</td><td>8.88</td><td>8.25</td><td>20.90</td><td>12.52</td><td>4.16</td></tr><tr><td>SiSi<sub>(110)</sub></td><td>8.11</td><td>9.68</td><td>SiC<sub>(110)</sub></td><td></td><td>12.11</td><td></td></tr><tr><td>SiC<sub>(100)</sub></td><td>SiSi<sub>(110)</sub></td><td>9.29</td><td>8.90</td><td>14.14</td><td>SiSi<sub>(110)</sub></td><td>6.17</td></tr><tr><td>SiC<sub>(110)</sub></td><td>Si<sub>TC</sub></td><td>9.48</td><td>7.78</td><td></td><td>12.39</td><td></td></tr></tbody></table>

for all potentials, with formation energy differences of about 1–1.5 eV compared to DFT results for all potentials.

Calculations of carbon interstitials using EDIP potentials lead to appropriate results. The most stable configurations are found to be either $CC_{(100)}$ or $CSi_{(100)}$ dumbbell interstitials, and all other dumbbell configurations are very close in energy. Tetrahedral interstitials are also found to be significantly higher in energy. With the Erhart–Albe potential, formation energy differences between all dumbbell configurations are significant, up to 5 eV between dumbbells with orientations $\langle 100 \rangle$ and $\langle 110 \rangle$. The Tersoff potential yields formation energies of dumbbell interstitials close to the DFT results, but the most stable configuration is the $C_{TSi}$ tetrahedral configuration, which is found unstable with DFT. Finally, the Gao–Weber potential provides a satisfactory description of carbon interstitials.

The properties of silicon interstitials are very difficult to reproduce with empirical potentials. DFT calculations showed that the most stable configuration is the $Si_{TC}$ tetrahedral interstitial, and that formation energies of stable configurations range between 7 and 9.3 eV. Only the Gao–Weber potential yields $Si_{TC}$ as the most stable interstitial configuration. Unfortunately, the computed formation energy is rather low compared to DFT results, and also much lower than the calculated values for carbon interstitials. Conversely, both Erhart–Albe and Tersoff potentials lead to very high formation energy values, more than 10 eV with the former. With both EDIP potentials, the $Si_{TC}$ tetrahedral interstitial is found to be unstable, relaxing to SiC split interstitials. However, the formation energies of all other interstitial configurations are

<table><caption>**Table 9.** Core energy differences (in eV/Burgers vector) for dislocations in silicon, computed with EDIP potentials, and compared to DFT calculations [45–48]. Several core configurations have been considered for the screw dislocation (A for shuffle, B for shuffle/glide, $C_1$ for glide single period, $C_2$ for glide double period), the 60° dislocation (S for shuffle, G for glide), the 30° and the 90° partials (SP for single period, DP for double period). For each dislocation, a configuration is chosen as the reference (zero energy).</caption><thead><tr><th></th><th>EDIP A</th><th>EDIP B</th><th>LDA [45]</th><th>GGA [46]</th><th>GGA [47]</th><th>LDA [48]</th></tr></thead><tbody><tr><td colspan="7">Screw (non-dissociated)</td></tr><tr><td>A</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td></td></tr><tr><td>C<sub>2</sub></td><td>-0.07</td><td>-0.06</td><td></td><td>-0.54</td><td></td><td></td></tr><tr><td>C<sub>1</sub></td><td>0.77</td><td>0.86</td><td>0.86</td><td></td><td>0.62</td><td></td></tr><tr><td>B</td><td>0.22</td><td>1.07</td><td>Unstable</td><td></td><td>Unstable</td><td></td></tr><tr><td colspan="7">60° (non-dissociated)</td></tr><tr><td>S</td><td>0.61</td><td>0.55</td><td></td><td></td><td>0.52</td><td></td></tr><tr><td>G</td><td>0</td><td>0</td><td></td><td></td><td>0</td><td></td></tr><tr><td colspan="7">30° (partial)</td></tr><tr><td>SP</td><td>0.33</td><td>0.37</td><td></td><td></td><td>0.52</td><td></td></tr><tr><td>DP</td><td>0</td><td>0</td><td></td><td></td><td>0</td><td></td></tr><tr><td colspan="7">90° (partial)</td></tr><tr><td>SP</td><td>-0.18</td><td>-0.09</td><td></td><td></td><td>0.05</td><td></td></tr><tr><td>DP</td><td>0</td><td>0</td><td></td><td></td><td>0</td><td></td></tr></tbody></table>

close to DFT results, and the average difference between silicon and carbon interstitial formation energies is also in close agreement with reference values.

Overall, the developed EDIP potentials allow for an improved description of point defects in silicon carbide compared to previous potentials. In particular, the formation energies of interstitials are in good agreement with DFT results. The only weak point is the instability of the $Si_{TC}$ configuration. We have tried to bias the potential fit by including the formation energy of this specific configuration into the fitting database, but with no success.

## 2.6. Extended defects

Covalent systems are characterized by the diversity of the possible structures for dislocation cores, which are usually difficult to model using classical interatomic potentials. The original EDIP has been shown to be one of the most accurate potentials for modelling dislocations in silicon [45]. The potential EDIP B has been tested on dislocation cores to determine whether this property was conserved (table 9). In addition, both A and B potentials have been used to compute the properties of dislocations in diamond (table 10) and cubic silicon carbide (table 11), for which very little reference information is available.

### 2.6.1. Silicon.
The structure and stability of screw, 60°, 30° and 90° dislocation cores have been considered in this work. All possible core geometries obtained are in good agreement with previous investigations, except for the screw core B which is found to be stable after relaxation using both EDIP potentials, whereas it is unstable with first-principles
computations. Nevertheless, this configuration has a high energy, especially with EDIP B. The most stable core is the C₂ configuration, in agreement with previous work [46]. However, the energy difference between the configurations C₂ and A, the second lowest energy structure, is much smaller with EDIP. The non-dissociated core of the 60° dislocation has also been investigated. We found a very good agreement with previous calculations, with very small energy differences between the glide and shuffle configurations.

The most studied cores in silicon are the 30° and 90° partial dislocations. The energy gain due to the reconstruction of the single period core into the double period core is computed to be in the range 0.3–0.4 eV with both EDIP potentials, in good agreement with first-principles results. The structures of both cores are also well reproduced with the potentials. The case of the 90° partial dislocation is more tricky, since it has been shown that the two possible core configurations, which exhibit a single (SP) or a double (DP) period, were almost degenerate in energy [48]. Whether the SP or DP core has the lowest energy should depend on the environment in which the dislocation is located [49]. Here, we found that the SP core is slightly favoured compared to the DP core. The energy difference is very small, especially for EDIP B, in very good agreement with previous calculations.

The original EDIP has been designed with the aim of improving the description of extended defects. Our calculations clearly show that this goal has been achieved. More importantly, the new set of parameters (B) leads to an even better description of dislocation cores in silicon.

### 2.6.2. Carbon (diamond).
For dislocations, we have considered the same core configurations as in silicon. It has to be noted that much less information is available for carbon than for silicon, either from experiments or first-principles calculations. First, for the screw dislocation, we found that the most stable configuration is the single-period glide core C₁, in agreement with previous first-principles calculations [4]. The energy difference with other configurations is large, a feature well reproduced by the EDIP potentials. To our knowledge, the double-period glide core C₂ has never been calculated in diamond. Our calculations suggest that it is not favoured in diamond, in contrast to silicon. Considering now the 60° dislocation, we found that the glide core is much more stable than the shuffle core, with larger energy differences than in silicon. The latter, computed with both EDIP potentials, are in close agreement with previous tight-binding calculations [50].

The 30° and 90° partial dislocations have also been studied with EDIP A and B potentials. Unlike in silicon, the single-period core SP is largely favoured compared to the double-period core DP in the case of the 30° partial dislocation. This result requires confirmation from electronic structure calculations. For the 90° partial dislocation, previous works indicate that the most stable configuration is the DP core, with an energy difference of the order of 0.4 eV with the SP core. Instead, our calculations performed with EDIP B potential show that the SP core is more stable than the DP core. With EDIP A, only the DP core is stable. The parameter set A could then be used for modelling 90° partials in diamond.

<table><caption>**Table 10.** Core energy differences (in eV/Burgers vector) for dislocations in diamond, computed with EDIP potentials, and compared to DFT [4, 51] and tight-binding [50] calculations. The same core notation as in table 9 is used.</caption><thead><tr><th></th><th>EDIP A</th><th>EDIP B</th><th>LDA [4]</th><th>DFTB [50]</th><th>LDA [51]</th></tr></thead><tbody><tr><td colspan="6" style="text-align: center;">Screw (non-dissociated)</td></tr><tr><td>A</td><td>1.66</td><td>1.11</td><td>1.29</td><td></td><td></td></tr><tr><td>C<sub>2</sub></td><td>1.34</td><td>1.05</td><td></td><td></td><td></td></tr><tr><td>C<sub>1</sub></td><td>0</td><td>0</td><td>0</td><td></td><td></td></tr><tr><td>B</td><td>1.01</td><td>0.63</td><td>2.51</td><td></td><td></td></tr><tr><td colspan="6" style="text-align: center;">60° (non-dissociated)</td></tr><tr><td>S</td><td>1.16</td><td>0.97</td><td></td><td>1.18</td><td></td></tr><tr><td>G</td><td>0</td><td>0</td><td></td><td>0</td><td></td></tr><tr><td colspan="6" style="text-align: center;">30° (partial)</td></tr><tr><td>SP</td><td>-1.25</td><td>-0.98</td><td></td><td></td><td></td></tr><tr><td>DP</td><td>0</td><td>0</td><td></td><td></td><td></td></tr><tr><td colspan="6" style="text-align: center;">90° (partial)</td></tr><tr><td>SP</td><td>Unstable</td><td>-0.48</td><td></td><td>0.35</td><td>0.3-0.5</td></tr><tr><td>DP</td><td>0</td><td>0</td><td></td><td>0</td><td>0</td></tr></tbody></table>

For dislocations in diamond, it is difficult to test the accuracy of the potentials developed, due to the lack of relevant data. Nevertheless, although non-dissociated perfect dislocations are well described, a less satisfactory agreement is noted in the case of the 90° partial dislocation.

### 2.6.3. Silicon carbide.
Finally, we have investigated the structure and stability of dislocation cores in silicon carbide with the EDIP potential. Compared to the elemental semiconductors previously considered, there are two possible configurations for 60°, 30° and 90°, depending on whether Si or C atoms are located in the core centre. Here we only investigated the zinc blende phase, but the core structures of dislocations located in the {111} planes of the cubic structure and in the basal planes of hexagonal polytypes are equivalent.

Our calculations with both EDIP potentials indicate that the most stable screw dislocation is the glide C₁ core, with small energy differences between A, C₁ and C₂ configurations. This result is in qualitative agreement with first-principles calculations since it has been previously shown that both A and C₁ could be the most stable core, depending on the choice of the exchange-correlation functional [4, 47]. This feature is not well reproduced by the Tersoff potential, which gives a very large energy for the A core relative to the C₁ core. With EDIP potentials, C₂ is found higher in energy than C₁, but there is no reference data to compare this result with. For the 60° dislocation, we found that, in the case of a carbon core, the glide configuration is largely favoured over the shuffle core for both parameter sets. For a silicon core, the glide configuration is also found to be much lower in energy with potential B. However, with potential A, the shuffle core becomes slightly favoured in energy.

No information is available on the 30° partial dislocations. Our EDIP calculations suggest that the single-period SP configuration is more stable for a silicon core. However, for a carbon core, potentials A and B yield different results:
<table><caption>**Table 11.** Core energy differences (in eV/Burgers vector) for dislocations in cubic silicon carbide, computed with EDIP potentials, and compared to Tersoff, DFT [4, 47] and tight-binding [52] calculations. The same core notation as in table 9 is used.</caption><thead><tr><th></th><th>EDIP A</th><th>EDIP B</th><th>Tersoff</th><th>LDA [4]</th><th>GGA [47]</th><th>DFTB [52]</th></tr></thead><tbody><tr><td colspan="7">Screw (non-dissociated)</td></tr><tr><td>A</td><td>0.14</td><td>0.19</td><td>2.33</td><td>-0.29</td><td>0.04</td><td></td></tr><tr><td>C<sub>2</sub></td><td>0.21</td><td>0.15</td><td>0.27</td><td></td><td></td><td></td></tr><tr><td>C<sub>1</sub></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td></td></tr><tr><td>B</td><td>Unstable</td><td>2.35</td><td></td><td>0.14</td><td></td><td></td></tr><tr><td colspan="7">60° (non-dissociated)/Si core</td></tr><tr><td>S</td><td>-0.2</td><td>1.39</td><td></td><td></td><td></td><td></td></tr><tr><td>G</td><td>0</td><td>0</td><td></td><td></td><td></td><td></td></tr><tr><td colspan="7">60° (non-dissociated)/C core</td></tr><tr><td>S</td><td>1.51</td><td>1.36</td><td></td><td></td><td></td><td></td></tr><tr><td>G</td><td>0</td><td>0</td><td></td><td></td><td></td><td></td></tr><tr><td colspan="7">30° (partial)/Si core</td></tr><tr><td>SP</td><td>-0.64</td><td>-0.45</td><td></td><td></td><td></td><td></td></tr><tr><td>DP</td><td>0</td><td>0</td><td></td><td></td><td></td><td></td></tr><tr><td colspan="7">30° (partial)/C core</td></tr><tr><td>SP</td><td>-0.19</td><td>0.26</td><td></td><td></td><td></td><td></td></tr><tr><td>DP</td><td>0</td><td>0</td><td></td><td></td><td></td><td></td></tr><tr><td colspan="7">90° (partial)/Si core</td></tr><tr><td>SP</td><td>Unstable</td><td>Unstable</td><td></td><td></td><td>0.3</td><td></td></tr><tr><td>DP</td><td>0</td><td>0</td><td></td><td></td><td>0</td><td></td></tr><tr><td colspan="7">90° (partial)/C core</td></tr><tr><td>SP</td><td>-0.04</td><td>0.03</td><td></td><td></td><td>0.3</td><td></td></tr><tr><td>DP</td><td>0</td><td>0</td><td></td><td></td><td>0</td><td></td></tr></tbody></table>

A favours the SP configuration whereas DP is more stable with potential B. Finally, 90° partial dislocations have been investigated and compared to tight-binding DFTB results [52]. The latter revealed that, for both Si and C cores, the DP configuration is more stable, with an energy difference of 0.3 eV per Burgers vector compared to the SP configuration. With both EDIP potentials, the DP configuration is found to be stable, and with the correct geometry. In the case of an Si core, the SP configuration is unstable. For a carbon core, the SP core is stable and almost degenerate in energy with the DP configuration.

In conclusion, we observed that the developed EDIP potentials enable a good description of both the structure and energetics of non-dissociated dislocations in silicon carbide. For the partial dislocations, the agreement is less satisfactory, although the most stable core configurations yielded by DFT or tight-binding calculations are found to be stable with low energies in both potentials.

## 2.7. Amorphous silicon carbide

The amorphous silicon carbide structure obtained with the two EDIP potentials has been investigated using classical molecular dynamics simulations at constant volume. The

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Global radial distributions $g(r)$ of amorphous SiC yielded by (a) EDIP A potential, Tersoff potential and *ab initio* calculations, (b) EDIP B potential and *ab initio* calculations.</figcaption></figure>

results have been compared to results from the Tersoff potential (T89) [10] and from *ab initio* molecular dynamics simulations in the DFT framework done by Finocchi *et al* [53]. The time step chosen was 0.7 fs. To obtain an amorphous sample a supercell containing 4096 atoms was melted at 5000 K for EDIP B and Tersoff potentials and at 5500 K for EDIP A during 7 ps, then quenched to 300 K at three different cooling rates: $T_x = 10^{13}$, $10^{14}$ and $10^{15}$ K s$^{-1}$. The system was then maintained in equilibrium for 14 ps and the global and partial radial distribution functions $g(r)$ calculated over 7 additional ps. The distribution functions were then averaged over 10 different amorphous samples and the graphs finally smoothed using a central moving average.

Almost no difference is observed in the global $g(r)$ and the partial distributions yielded by the Tersoff potential at the three quenching rates. Slightly larger differences are observed between the various $T_x$ for EDIP A. The relative intensities of the various peaks depend on the quenching rate: when $T_x$ decreases the numbers of C–C and Si–Si bonds decrease while the number of Si–C bonds increases, as expected from equilibrium considerations. For EDIP B, only the second-neighbour C–C–C peak is affected by the quenching rate: its intensity decreases when $T_x$ decreases.

In the following, we have chosen to analyse the amorphous materials obtained for $T_x = 10^{15}$ K s$^{-1}$. It is the order of magnitude of the quenching rate used in the *ab initio* molecular dynamics calculations [54, 55] and we want the comparison between the potential and DFT results to be as meaningful as possible. We have represented in figure 1 the global $g(r)$
yielded by EDIP A and B potentials, Tersoff potential and *ab initio* calculations.

As can be seen in figure 1(a) the local environment of the amorphous material yielded by the EDIP A and Tersoff potentials are in very good agreement with the DFT results. Differences appear mainly in the relative intensities of the first two peaks: the C–C peak is less intense for the potentials than in *ab initio*. Moreover, it can be seen that the description of the medium-range order is slightly better for EDIP A than for Tersoff. In contrast, significant differences are observed between the *g(r)* obtained using EDIP B potential and the *ab initio* results, as represented in figure 1(b): the C–C and Si–Si peaks are significantly more intense than the Si–C one, the Si–C peak is displaced and the medium-range order is poorly described. In particular, there is a rather intense peak at 2.4 Å, which does not exist in DFT. This shows that the SiC amorphous structure is crudely described by the EDIP B potential.

To go further into detail in the comparison of potentials with *ab initio*, we have analysed the partial distributions *g*(C–C), *g*(Si–C) and *g*(Si–Si). As can be seen from figure 2(a), the *g*(C–C) partial distributions obtained using the three potentials and DFT are relatively similar, although the intensities of the peaks obtained using the EDIP B potential are significantly more intense than for the other methods. A first intense peak is observed at 1.50, 1.52, 1.49 and 1.50 Å for EDIP A, EDIP B, Tersoff potentials and DFT, respectively. These maxima can be compared with the equilibrium distances in graphite and diamond, 1.42 and 1.54 Å, respectively. This seems to indicate that the carbon atoms in amorphous SiC are mainly hybridized *sp*³. The second peak has a smaller maximum, occurring at 2.53, 2.62, 2.55 and 2.51 Å for EDIP A, EDIP B, Tersoff potentials and DFT, respectively. This distance is close to the second-neighbour distance in diamond (2.51 Å). This peak therefore corresponds to C atoms in the C–C–C configuration. A very faint third maximum (best defined for the EDIP A potential and DFT) corresponding to C–Si–C configurations can be seen around 3 Å. The C–C distributions yielded by Tersoff and EDIP A are therefore in very good agreement with the *ab initio* results. The EDIP B potential slightly overestimates the C–C bond lengths and the peak intensities, but the results are still in reasonable agreement.

The Si–Si partial distributions obtained using the various methods show larger differences, as observed in figure 2(b). The first-neighbour peak is at the same position: 2.37, 2.43, 2.41 and 2.32 Å for EDIP A, EDIP B, Tersoff potentials and DFT, respectively, even if the intensities differ. For all three potentials the Si–Si bonds are slightly elongated compared to the bonds in silicon (*d* = 2.35 Å). A second peak is observed at 2.86, 2.86, 3.04 and 2.97 Å for EDIP A, EDIP B, Tersoff potentials and DFT, respectively. For all three potentials this peak is more defined and intense than the first-neighbour peak, especially for EDIP B, and *g*(Si–Si) comes to 0 between the two peaks. These two features, which are probably due to the potential cutoff functions, are responsible for a less than satisfactory agreement between the potentials and DFT. The number of Si–Si bonds, however, is rather small in amorphous SiC, which explains its good overall description by EDIP A and Tersoff.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Partial *g*(C–C), *g*(Si–C) and *g*(Si–Si) distribution functions of amorphous SiC yielded by the EDIP A, EDIP B and Tersoff potentials, as well as DFT calculations.</figcaption></figure>

Finally, the *g*(Si–C) distributions obtained using Tersoff, EDIP A and EDIP B potentials, as well as DFT, are represented in figure 2(c). It can be seen that EDIP A and Tersoff *g*(Si–C) are very similar and in good agreement with the DFT results. They exhibit a first intense peak at 1.88 and 1.90 Å, respectively, a smaller and broader second peak at 2.75 and 2.84 Å, and a very faint maximum around 3.4 Å. This must be compared to the two peaks and the shoulder observed in DFT at 1.90, 2.81 and 3.4 Å. Moreover, the relative intensities of the first two peaks yielded by the two potentials are very similar to the *ab initio* results. The only slight difference is that the second peak is better defined for Tersoff than for EDIP A and DFT. The first peak is located at the first-neighbour distance in *β*-SiC, while the second and third minima correspond to C–C–Si and C–Si–Si configurations. For EDIP B, in contrast,
<table><caption>**Table 12.** Partial coordination numbers and percentage of homonuclear bonds yielded by the integration of the first peaks of the *g* (C–C) and *g* (C – Si) radial distribution functions obtained using EDIP potentials. Results are compared to the results of Tersoff (T89) potential and DFT [53].</caption><thead><tr><th></th><th>EDIP A</th><th>EDIP B</th><th>T89</th><th>DFT</th></tr></thead><tbody><tr><td>n<sub>C-C</sub></td><td>1.6</td><td>2.9</td><td>1.3</td><td>1.8</td></tr><tr><td>n<sub>C-Si</sub></td><td>1.8</td><td>1.5</td><td>2.2</td><td>2.1</td></tr><tr><td>n<sub>C</sub></td><td>3.4</td><td>4.4</td><td>3.5</td><td>3.9</td></tr><tr><td>%<sub>homo</sub></td><td>45</td><td>66</td><td>37</td><td>46</td></tr></tbody></table>

two peaks with similar intensity and closer to each other are observed at 1.98 and 2.36 Å, and there is no peak around 2.8 Å. The first peak corresponds to an elongated Si–C bond. The origin of the second peak is more difficult to ascertain: it could be Si–Si bonds constrained in a Si–C configuration. This feature is certainly at the core of the poor description of the silicon carbide amorphous structure yielded by EDIP B.

We have also determined the partial coordination numbers, i.e. the average numbers of C–C and C–Si bonds formed by C atoms yielded by the various potentials, by integrating the first C–C and Si–C peaks of the *g* (C–C) and *g* (C – Si) distributions. These partial coordinations, as well as the proportion of homonuclear bonds formed by carbon atoms, are shown in table 12. It can be seen that EDIP A describes well the C coordination in the material, and slightly better than the Tersoff potential. EDIP A yields undercoordinated carbon atoms, in agreement with the *ab initio* results, and the proportion of homonuclear C–C bonds is very close to the proportion obtained in DFT. The proportion of Si–C bonds, however, is slightly underestimated by EDIP A and is better described by T89. The EDIP B potential, in contrast, yields overcoordinated C atoms because of the excessive number of C–C bonds formed.

In conclusion, the analysis of the total and partial radial distribution functions, as well as the coordination numbers, shows that the EDIP A potential yields an amorphous silicon carbide with properties in very good agreement with the *ab initio* results. It describes in particular the high proportion of homonuclear bonds formed in this material. The EDIP A potential is in this case slightly better than the Tersoff potential. In contrast, EDIP B is not well suited for an accurate description of the amorphous SiC structure.

### 3. Conclusion

We proposed a generalization of the semi-empirical interatomic potential EDIP, initially designed for silicon [16, 18], for modelling silicon carbide. A force and energy computing routine is available, allowing for an easy integration into existing simulation codes [56]. Two parametrizations have been obtained: the first is an extension of the original parameter set and requires the introduction of an additional cutoff function (set A). In the second one, all the parameters for Si–Si, C–C and Si–C interactions are new (set B). These two potentials have been tested for silicon, carbon (diamond) and silicon carbide for a wide range of properties, including bulk properties,

stability and structure of high-pressure crystalline phases, formation energy and structure of point and extended defects, and structure of amorphous phases. Our results have been compared to data from the literature when available. Note that in a few cases, especially for dislocations, our work also leads to predictions.

We found that the two proposed potentials A and B enable an accurate description of point and extended defects in silicon, diamond and silicon carbide, thus preserving the initial focus of the original EDIP potential. More specifically, an extra feature of the parameter set B is the very good description of high-pressure phases of the different materials. The potential with the parameter set A is well suited for describing the disordered structure of silicon carbide. Therefore, this potential would be fully appropriate for modelling irradiation-induced effects in silicon carbide, such as the generation of structural defects and amorphization. It also has the advantage of keeping the original EDIP parametrization for silicon.

### Acknowledgment

This work was funded by the joint research programme 'ISMIR' between CEA and CNRS.

### References

[1] Willander M, Friesel M, Wahab Q and Straumal B 2006 *J. Mater. Sci. Mater. Electron.* **17** 1

[2] Katoh Y, Snead L L, Henager C H Jr, Hasegawa A, Kohyama A, Riccardi B and Hegeman H 2007 *J. Nucl. Mater.* **367** 659

[3] Madelung O, Rössler U and Schulz M (ed) 2001 *Landolt-Börnstein: Numerical Data and Functional Relationship in Science and Technology (New Series, Group III vol 41A)* (Berlin: Springer)

[4] Pizzagalli L, Beauchamp P and Rabier J 2005 *Europhys. Lett.* **72** 410

[5] Cicero G, Pizzagalli L and Catellani A 2002 *J. Phys.: Condens. Matter* **14** 13031

[6] Kohler C 2002 *Phys. Status Solidi* b **234** 522

[7] Bertolus M, Ribeiro F and Defranceschi M 2007 *Eur. Phys. J. B* **60** 423

[8] Schmidt A A, Kharlamov V S, Safonov K L, Trushin Yu V, Zhurkin E E, Cimalla V, Ambacher O and Pezoldt J 2005 *Comput. Mater. Sci.* **33** 375

[9] Stillinger F and Weber T 1985 *Phys. Rev. B* **31** 5262

[10] Tersoff J 1989 *Phys. Rev. B* **39** 5566

[11] Brenner D W 1990 *Phys. Rev. B* **42** 9458

[12] Devanathan R, de la Rubia T D and Weber W J 1998 *J. Nucl. Mater.* **253** 47

[13] Gao F and Weber W J 2002 *Nucl. Instrum. Methods Phys. Res. B* **191** 504

[14] Erhart P and Albe K 2005 *Phys. Rev. B* **71** 035211

[15] Bazant M and Kaxiras E 1996 *Phys. Rev. Lett.* **77** 4370

[16] Bazant M, Kaxiras E and Justo J 1997 *Phys. Rev. B* **56** 8542

[17] Bazant M 1997 Interatomic forces in covalent solids *PhD Thesis* Harvard University, Cambridge, MA

[18] Justo J, Bazant M, Kaxiras E, Bulatov V and Yip S 1998 *Phys. Rev. B* **58** 2539

[19] Marks N A 2001 *Phys. Rev. B* **63** 035401

[20] Carlsson A E, Fedders P A and Myles C W 1990 *Phys. Rev. B* **41** 1247

[21] Mistriotis A D, Flytzanis N and Farantos S C 1989 *Phys. Rev. B* **39** 1212
[22] Kirkpatrick S D, Gelatt S C D and Vecchi M P 1983 *Science* **220** 671

[23] Cerny V 1985 *J. Opt. Theory Appl.* **45** 41

[24] Nelder J A and Mead R 1965 *Comput. J.* **7** 308

[25] http://xmd.sourceforge.net/

[26] Crain J, Clark S J, Ackland G J, Payne M C, Milman V, Hatton P D and Reid B J 1994 *Phys. Rev.* B **49** 5329

[27] Yin M T 1984 *Phys. Rev.* B **30** 1773

[28] Simmons G and Wang H 1971 *Single Crystal Elastic Constants and Calculated Aggregate Properties: A Handbook* (Cambridge, MA: MIT Press)

[29] Nielsen O H and Martin R M 1985 *Phys. Rev.* B **32** 3792

[30] Balamane H, Halicioglu T and Tiller W A 1992 *Phys. Rev.* B **46** 2250

[31] Furthmüller J, Hafner J and Kresse G 1994 *Phys. Rev.* B **50** 15606

[32] Lambrecht W R L, Segall B, Methfessel M and van Schilfgaarde M 1991 *Phys. Rev.* B **44** 3685

[33] Kamitani K, Grimsditch M, Nipko J C, Loong C-K, Okada M and Kimura I 1991 *J. Appl. Phys.* **82** 3152

[34] Karch K, Pavone P, Windl W, Schütt O and Strauch D 1994 *Phys. Rev.* B **50** 17054

[35] Park C H, Cheong B-H, Lee K-H and Chang K J 1994 *Phys. Rev.* B **49** 4485

[36] Karch K, Bechstedt F, Pavone P and Strauch D 1996 *Phys. Rev.* B **53** 13400

[37] Yoshida M, Onodera A, Ueno M, Takemura K and Shimomura O 1993 *Phys. Rev.* B **48** 10587

[38] Catti M 2001 *Phys. Rev. Lett.* **87** 035504

[39] Pandey K C 1986 *Phys. Rev. Lett.* **57** 2287

[40] Leung W-K, Needs R J, Rajagopal G, Itoh S and Ihara S 1999 *Phys. Rev. Lett.* **83** 2351

[41] Goedecker S, Deutsch T and Billard L 2002 *Phys. Rev. Lett.* **88** 235501

[42] Zywietz A, Furthmüller J and Bechstedt F 1998 *Phys. Status Solidi* b **210** 13

[43] Breuer S J and Briddon P R 1995 *Phys. Rev.* B **51** 6984

[44] Lucas G and Pizzagalli L 2007 *Nucl. Instrum. Methods* B **255** 124

[45] Pizzagalli L, Beauchamp P and Rabier J 2003 *Phil. Mag.* **83** 1191

[46] Wang C-Z, Li J, Ho K-M and Yip S 2006 *Appl. Phys. Lett.* **89** 051910

[47] Pizzagalli L and Beauchamp P 2008 unpublished

[48] Bulatov V V, Justo J F, Cai W, Yip S, Argon A S, Lenosky T, de Koning M and Diaz de la Rubia T 2001 *Phil. Mag.* A **81** 1257

[49] Lehto N and Öberg S 1998 *Phys. Rev. Lett.* **80** 5568

[50] Blumenau A T, Heggie M I, Fall C J, Jones R and Frauenheim Th 2002 *Phys. Rev.* B **65** 205205

[51] Blaze X, Lin K, Canning A, Louie S G and Chrzan D C 2000 *Phys. Rev. Lett.* **84** 5780

[52] Blumenau A T, Fall C J, Jones R, Öberg S, Frauenheim Th and Briddon P R 2003 *Phys. Rev.* B **68** 174108

[53] Finocchi F, Galli G, Parrinello M and Bertoni C M 1992 *Phys. Rev. Lett.* **68** 3044

[54] Finocchi F 2005 private communication

[55] Finocchi F and Galli G 1994 *Phys. Rev.* B **50** 7393

[56] http://theory.phymat.sp2mi.univ-poitiers.fr