A new parametrization of the Stillinger–Weber potential for an improved description of defects and plasticity of silicon

With the SW potential [2], the total energy E is obtained as a combination of two-body and three-body interactions:

$$\begin{eqnarray} E&=&\mathop{\sum }\limits _{i\lt j}{\Phi }_{2}(i,j)+\mathop{\sum }\limits _{{i\not = j\atop j\lt k}}{\Phi }_{3}(i,j,k) \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} &=&\mathop{\sum }\limits _{i\lt j}\varepsilon {\phi }_{2}({r}_{i j})+\mathop{\sum }\limits _{{i\not = j\atop j\lt k}}\varepsilon {\phi }_{3}({r}_{i j},{r}_{i k},{\theta }_{i j k}). \end{eqnarray} \tag{ 2 }$$

The two- and three-body ϕ functions are

$$\begin{eqnarray} &&\fl {\phi }_{2}({r}_{i j})=A\left [B\left (\frac{\sigma }{{r}_{i j}}\right )^{p}-\left (\frac{\sigma }{{r}_{i j}}\right )^{q}\right ]\exp \left (\frac{1}{{r}_{i j}/\sigma -a}\right ) \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} &&\fl {\phi }_{3}({r}_{i j},{r}_{i k},{\theta }_{i j k})=\lambda \left (\cos {\theta }_{i j k}+\frac{1}{3}\right )^{2}\exp \left (\frac{\gamma }{{r}_{i j}/\sigma -a}\right )\nonumber\\ &&\fl \phantom {\phi _3\:(r_{ij},r_{ik},\theta _{ijk}) =} \times \ \exp \left (\frac{\gamma }{{r}_{i k}/\sigma -a}\right ). \end{eqnarray} \tag{ 4 }$$

The choice of the two-body function is not really justified in the original paper from Stillinger and Weber, only that its derivatives exhibit no discontinuities at the cutoff a [2]. The latter is also true for the three-body function, which is built in such a way that it is positive and vanishes for a perfect cubic diamond crystal when bond angles are equal to the tetrahedral angle defined as cosθ = 1/3. ε and σ are energy and distance scaling parameters, respectively, which set the cohesive energy and lattice constant of the cubic diamond silicon. The other remaining parameters, A,B,p,q,λ and γ, allow us to tune two-body and three-body interactions as a function of interatomic distances and angles.

Since the three-body interactions cancel for the perfect crystal, properties conserving the cubic diamond lattice do not depend on the value of λ and γ. This is the case for the bulk modulus, in addition to the cohesive energy and the lattice constant. This also allows for increasing or decreasing the strength of three-body interactions (typically by changing λ) with respect to two-body interactions, without altering the cubic diamond ground state. This property has been exploited to improve the fracture behaviour [21], and to investigate how the strength of the three-body interactions is related to the plasticity of amorphous silicon [20]. Overall, it seems that increasing λ leads to better results regarding several properties related to fracture [21] and amorphous systems [19, 20]. However, the excess energy of any system with respect to the perfect cubic diamond crystal, such as the formation energy of defects, is linearly proportional to λ. This has obvious consequences, for instance increasing λ from 21 to 42 (all other parameters being kept constant), as in the work of Holland and Marder [21], implies that all excess energies will be twice the original SW potential.

The original set of parameters is reported in table 1, along with the parameters set from Vink and co-workers [22] and the one proposed in this work (SW_m). The latter has been obtained after a thorough analysis of the influence of λ on many properties such as the stability of point and extended defects, and an analysis of the previous investigations mentioned above [21, 22]. We found that λ should be greater than or equal to 31 in order to obtain the correct energy ordering of dislocation configurations [15, 16, 24]. At the same time, λ should be kept as low as possible to minimize the increase of excess energies. Then we set λ = 31 before optimization. We also kept the dimensionless cutoff parameter a equal to its original value. The remaining parameters for SW_m are obtained using a least-square fitting algorithm and a database including elastic constants and the formation energies of silicon interstitials and a monovacancy. We also determined the scaling parameters ε and σ in order to recover the experimental cohesive energy and lattice constant of cubic diamond silicon. We first tried to keep the original exponent values p = 4 and q = 0, which are computationally very convenient. Nevertheless, it was not possible to obtain a good agreement in that case, and this constraint was lifted. Finally, the best results were obtained for the set of parameters SW_m reported in table 1.

In table 2, the results of basic properties computed with the original SW, SW_m and VBWM potentials are shown. The experimental lattice constant is obtained for all sets. Regarding the cohesive energy for the original SW potential, it is slightly higher than the reference data, since ε has been fitted to reproduce the melting temperature of the crystal in that case [2]. An even higher value is obtained with the VBWM set, due to the low value of ε. Instead, our potential yields the correct cohesive energy by construction. Elastic constants are obtained in fair agreement with experiments for the SW and SW_m potentials, although C₁₂ is too high and C₄₄ is too low. In the case of the VBWM potential, the computed elastic constants are all quite low, leading to a silicon crystal softer than expected.

The three potentials have then been tested for reproducing the formation energies of the following point defects: mono- and di-vacancy, and the most stable configurations of the self-interstitial. For the reference values, we considered experiments or state of the art calculations when the former are nonavailable. The results are reported in table 3. In the case of vacancies, the best agreement is obtained for the SW_m and VBWM potentials, while the original SW potential yields too low formation energies. For self-interstitials, I〈110〉 is found as the most stable configuration by all potentials, in agreement with recent quantum Monte Carlo (QMC) calculations [32]. However, in these calculations, all self-interstitial configurations appear almost degenerate, which is not the case here with an energy range of about 2 eV for all potentials. The best result is obtained for the original SW, followed by SW_m. The VBWM potential yields the highest formation energies, due to the fast increase of the excess energy as a function of λ.

An important quantity in the context of irradiation is the threshold displacement energy (TDE), which is the minimum kinetic energy to be transferred to a lattice atom for displacing it from its original site to a new metastable defect position. It allows us to quantify the lattice resistance to the creation of point defects, and is a required input for large-scale calculations of damage during irradiation. TDE is a highly anisotropic quantity in covalent materials, and a well-constructed average can be used [31]. Although the determination of TDE is computationally intensive, accurate first-principles calculations have been recently achieved [31, 33, 34].

In table 4 are reported computed TDE values, averaged over 1000 uniformly randomly chosen lattice directions or corresponding to high symmetry orientations, determined using the different potentials considered here, and the reference values of Holmström et al [31]. Considering first the averaged data, the new potential proposed here clearly outperforms the original SW potential, yielding too low values, and the VBWM potential, yielding too high values. Regarding now specific directions, the most important one is 〈111〉(B), since it corresponds to the absolute minimum of possible TDE values, i.e. it is the direction of the easiest defect creation in the silicon lattice. It appears that all potentials lead to values greater than the first-principles reference, the best agreement being obtained for SW_m and SW.

Table 4.  Threshold displacement energies (eV): average values for the creation of Frenkel pairs ($\overline{{E}_{d}}(\mathrm{FP})$), and for the creation of either Frenkel pairs or bond defects ($\overline{{E}_{d}}(\mathrm{FP}/\mathrm{BD})$), as well as values along two high symmetry orientations (there are two values for 〈111〉 because the two possible directions are not equivalent in the cubic diamond structure).

	$\overline{{E}_{d}}(\mathrm{FP})$	$\overline{{E}_{d}}(\mathrm{FP}/\mathrm{BD})$	〈100〉	〈111〉(A)	〈111〉(B)
Ref. [31]	36.0 ± 2.0	24.0 ± 2.0	20.5 ± 1.5	14.5 ± 1.5	12.5 ± 1.5
SW	30.3 ± 0.3	20.6 ± 0.3	23.5 ± 0.5	20.5 ± 0.5	17.5 ± 0.5
VBWM	41.0 ± 0.4	27.4 ± 0.2	23.5 ± 0.5	26.5 ± 0.5	20.5 ± 0.5
SWm	37.7 ± 0.4	24.5 ± 0.2	22.5 ± 0.5	22.5 ± 0.5	17.5 ± 0.5

We now focus on the performance of the potentials for reproducing plasticity-related properties. We first determined the relative stability of the possible dislocation core configurations. This is a stringent test since interatomic potentials are usually hardly able to reproduce with accuracy the energies associated with the highly distorted core geometries. Table 5 shows the core energy differences for different dislocation cores for the different potentials and first-principles calculations [15, 35–39]. These have been obtained by relaxing a single dislocation using cylinder-like boundary conditions and a large system including 3360 atoms in order to avoid technical issues associated with small cells [40]. In a few cases, additional calculations using quadrupole dislocations and periodic boundary conditions were performed, with no noticeable changes.

The geometries of all stable dislocation cores are shown in figure 1. First, we focus on the non-dissociated screw dislocation, which is important for the low temperature plasticity of silicon [16, 41–46]. All potentials found the sessile core C₂ (figure 1) as the most stable one, albeit with half energy differences compared to the first-principles reference. This core exhibits a double-period reconstruction along the dislocation line, seemingly close to the double-period DP core of the 90° partial dislocation. A, the second stable core, is slightly higher in energy. A known issue of the original SW potential is that the core B [47], unstable with first-principles calculations, is more stable than the shuffle core A [15, 48, 49]. With the new potential SW_m and VBWM, the B core becomes unstable as in first-principles calculations.

Zoom In Zoom Out Reset image size

Another important dislocation is at 60°, which is expected to be first nucleated during the relaxation of epitaxied thin films [50–53]. Recent investigations have revealed that the core S₁, proposed in the pioneering work of Hornstra [54], was not stable, and that the most stable configuration was a sessile G core [36]. No potentials were able to reproduce the first result. For the second one, only the VBWM potential yields the right behaviour, since with both SW and SW_m another sessile configuration, S₃, corresponds to the lowest core energy. This last configuration is quite complex, but ensures the complete passivation of dangling bonds in the dislocation core, as for the G core.

Finally, we have investigated the energetics of 30° and 90° partial dislocations, which govern the plasticity of silicon at high temperatures [26]. The 30° is known to be reconstructed along the dislocation line, while it has been shown that two configurations, with either a single- (SP) or a double-period (DP) along the dislocation line, are quasidegenerate in energy for the 90° partial dislocation [40]. The SW_m potential appears as the best candidate in that case, since it yields the correct reconstruction energy of the 30° partial dislocation, and the two possible core configurations for the 90° partial dislocation are quasidegenerate in energy, in excellent agreement with first-principles calculations.

Another criterion associated with an appropriate modelling of plasticity is the ability to reproduce the behaviour of silicon bulk under high shear stress. This property is, for instance, important for dislocation nucleation [50, 51]. It has been shown in earlier investigations that, despite some limitations, the original SW potential was one of the most efficient ones [14, 55]. In order to check the performance of the new SW_m potential, the critical stress associated with the homogeneous shear of {111} planes has been computed for three high symmetry directions, taking into account the relaxation of the volume. The results are reported in table 6 for the three potentials, and compared to first-principles DFT-LDA calculations [56] made with the SIESTA code [57].

For the easiest directions, 〈110〉 and 〈112〉-T, all potentials tend to slightly overestimate the critical shear stress, and the best quantitative agreement is obtained for the original SW potential. Unfortunately, it also favours an easy shearing along the 〈112〉-AT direction, which is clearly in discrepancy with DFT-LDA results. The correct behaviour, i.e. the lowest stress value for the 〈112〉-T direction followed by 〈110〉 and 〈112〉-AT, is recovered for the SW_m and VBWM potentials. Overall, there is a qualitative improvement using the new parametrization SW_m, or the VBWM potential, compared to the original SW.

Although they cannot be directly related to a real physical mechanism, generalized stacking faults (GSFs) energy surfaces [59] are often computed because they can be combined with Peierls–Nabarro models [60, 61] to predict dislocation core properties [62–64]. Their correct modelling is then a good criterion to determine the aptitude of a given potential for dislocation-related properties. The original SW potential has been shown to be well suited for calculating GSF energy surfaces in good agreement with more accurate calculations. In particular, the SW potential leads to smooth energy variations as in first-principles simulations, a feature which is not well reproduced by more recent potentials [14, 55].

Figure 2 represents the GSF energy curves along the two most important slip systems. The smooth aspect of the curves for the three potentials confirms that this characteristic is due to the simple analytical formulation of the SW potential [14]. For the shuffle 〈110〉{111}, the original SW potential clearly underestimates the energy variation. A clear improvement is obtained for both the SW_m and VBWM potentials, although the computed maxima remain slightly lower than first-principles results. For the second slip system, all potentials overestimate the reference data. The best agreement is obtained with the VBWM potential, followed by the SW and SW_m potentials. Overall, there is a bare gain in performance using the new SW_m potential compared to the original parametrization, but it worth recalling that the latter is already superior to other potentials [14].

Zoom In Zoom Out Reset image size

The ability of the different potentials to describe fracture has been studied by computing the propagation of $(1 1 1)[\bar {2}1 1]$ cracks. $\{ 1 1 1\} \langle \bar {2}1 1\rangle$ is an easy cleavage system in silicon, for which fracture is stable in a wide range of energy release rates G (2 J m⁻² ≤ G ≤ 16 J m⁻²) [65]. Cleaved surfaces are atomically flat with few defects which are predominantly steps with a specific symmetry [66]. In pure tensile mode and at low crack velocity (<800 m s⁻¹), a systematic deviation from the original fracture plane is observed [67]. These experiments have been successfully reproduced by hybrid quantum/classical calculations [68], allowing for physical explanations of the observed behaviours [69, 70].

We describe here the procedure used since fracture simulations are not standard. Our set up consists in a quasi-2D system of 703 872 atoms with dimensions ≃138 nm along $[\bar {2}1 1]\times$, 132 nm along [111], and 0.7681 nm along $[0 1\bar {1}]$. Periodic boundary conditions are applied along the last dimension. A $(1 1 1)[\bar {2}1 1]$ notch was introduced in the sample with an initial crack tip position at almost 1/3 of the total length. This system was successively loaded by an homogeneous rescaling of the coordinates in pure tensile mode (mode I), the load being conserved by fixing positions of atoms belonging to top and bottom planes in the [111] direction. Starting from an initial value of 1.8%, the tensile strain was increased in steps of 0.1%, alternating with atomic position relaxation using a damped molecular dynamics algorithm and a maximum force criterion of 3 × 10⁻⁴ eV Å⁻¹. When crack tip instability occurred, a non-damped molecular dynamics was performed to monitor the crack propagation. With this procedure we do not aim to reproduce the experimental minimal energy release rate for fracture since our simulations include periodicity along the crack front, temperature close to 0 K, and short simulation times (≃15 ps). However, it allows for a rigorous comparison of fracture behaviour for the three potentials in the same conditions.

The SW potential required the largest load (4.10% ± 0.05% strain or energy release rate G = 17.35 J m⁻² ± 0.45 J m⁻²) for crack propagation, due to an efficient spreading of stress at the crack tip following three successive reconstruction stages. The tip structure just before crack opening is characterized by the presence of heptagonal rings, forming elongated atomic chains along the [111] direction (figure 3(a)). The propagating crack quickly becomes highly disordered, with a rather round and amorphous tip structure and the eventual formation of atomic chains along the [111] direction. The cleaved surfaces are also highly disordered and far from being atomically flat (figure 3(d)).

Zoom In Zoom Out Reset image size

A significantly different behaviour is obtained using the VBWM potential (figures 3(b), (e)). In that case, the crack tip starts to propagate when the loading strain is 2.15% ± 0.05% (or G = 4.40 J m⁻² ± 0.2 J m⁻²). For lower values, no tip reconstructions occur, in disagreement with electronic structure calculations [69, 70]. Nonetheless, the crack propagation with this potential is fully brittle, leaving atomically flat surfaces. At higher loads (2.70% strain), one surface remains non-reconstructed whereas the other exhibits the alternating five-atom and seven-atom rings of the Pandey reconstruction [70].

Finally, using the SW_m potential, we found that the crack tip reconstructs upon loading (figure 3(c)), but to a lesser extent than the SW potential, and leading to a structure quite similar to electronic structure results [69]. Crack tip propagation then occurs for a strain of 3.65% ± 0.05% (G = 14.40 J m⁻² ± 0.45 J m⁻²). During the very first stages, a seemingly brittle fracture is obtained, with formation of Pandey-like reconstructions on one surface as in electronic calculations, and atomically flat cleaved surfaces. In this initial dynamic stage, some defects or reconstructions occasionally occur on one surface and among them some structures characteristic of the rotated Pandey reconstructions which were also observed with more precise simulation schemes. However, after 6–9 ps, the fracture process becomes increasingly less brittle, with the blunting of the crack tip and the formation of large defected zones including steps of several atomic planes (figure 3(f)).

The comparison of the three potentials reveals that there is a relation between the maximal load attained and the propagation behaviour of cracks. In fact, high strain values are required for the SW and SW_m potentials, which results in a large strain energy released during crack propagation and the dynamical formation of many defects at the crack tip after a short time. Conversely, propagation is obtained for a much lower strain value in the case of the VBWM potential, allowing for a fully brittle crack propagation. Compared with experimental observations of cleaved surfaces, it appears that this potential is the more appropriate for describing crack propagation. Concerning the structure of the crack tip, electronic structure calculations revealed a specific crack tip reconstruction before propagation. This feature is best reproduced by the SW_m potential, but missed by the VBWM potential. Then none of the tested potentials is completely satisfactory, which is hardly surprising given the complexity of the fracture mechanism.

It is known that the transition pressure needed to transform the cubic diamond structure to β-tin computed using DFT-LDA/GGA ranges from 7 to 9 GPa, i.e. lower than experimental values of 10–12 GPa [72, 73]. Conversely, transition pressures of 15–16 GPa are predicted by QMC techniques [71, 73]. The computed β-tin volume is about 15 Å³/atom with all approaches, the pressure discrepancy coming from the energy difference between the two phases (about 0.5 eV/atom for QMC and 0.25 eV/atom for DFT-LDA, see table 7) [71, 73].

With the original SW potential, the energy difference is comparable to the DFT-LDA result, but the volume is too large (obtained as a²c from values in table 7). The transition pressure (only roughly estimated here as the ratio between energy and volume differences) is likely to be in fair agreement with experiments. For the other potentials, the energy difference is much higher (due to the larger λ), close to the QMC results, but the computed volume is much larger, around 18 Å³/atom. Then the pressure transition will be much higher than in experiments, between 30 and 40 GPa. This is clearly a failure of these potentials which are not able to reproduce the compact metallic character of the high pressure phase. Nevertheless, it can be considered an advantage for plasticity-oriented simulations, where one wants to use high stresses (to speed up MD simulations), without triggering phase transition.

The melting temperature was obtained by molecular dynamics calculations. We searched for the temperature at which the cubic diamond and the liquid phases coexist, i.e. the temperature at which the free energies of solid and liquid phases are equal. Our computations were performed at fixed volume, whereas they should ideally be performed at constant pressure. This may lead to a slight overestimation of the melting temperature.

We found a melting temperature of 1720 K for the original SW potential, in excellent agreement with the experimental value of 1687 K. This is not surprising since the ε parameter has been scaled to reproduce this property [2]. For the SW_m potential, we found that the crystal was stable until 2610 K. Finally, an even higher value of 2850 K was determined using the VBWM potential. For both potentials, the high value of the λ parameter greatly enhances the stability of the cubic diamond phase compared to the liquid. It is similar to what we obtained for the high pressure phase transition.

Obviously, it is difficult to recommend using the SW_m potential for solid–liquid transitions. Nevertheless, a too high melting temperature could be a blessing in disguise. In fact, computational limitations severely restrict the duration of molecular dynamics simulations, which often prevents the occurrence of thermally activated events. A common trick is to use temperatures higher than in experiments to increase occurrence probabilities, assuming that there are no mechanism changes as a function of temperature. A high melting temperature then allows us to do this trick, without the risk of melting the crystalline structure.

The ability of the SW_m potential to describe the amorphous state of silicon has also been investigated. We have prepared an a-Si sample by using the following procedure: (1) crystalline silicon is melted at a temperature of 4500 K in the NVE ensemble; (2) after equilibration, the liquid state is quenched to 0 K in the NPT ensemble at a rate of 10¹¹ K s⁻¹, typical of laser experiments [19]; (3) finally a conjugate gradient minimization is performed while keeping the density constant. The radial distribution function (RDF) of the final system is reported in figure 4, and compared to measurements and results for the original SW and VBWM potentials. In the case of the SW potential, the RDF curve is characterized by a non-zero minimum between the first and second neighbour peaks, as well as an additional peak at approximately 3.3 Å. These features are not present in a-Si models obtained with both SW_m and VBWM potentials, in agreement with experiments.

Zoom In Zoom Out Reset image size

The a-Si sample obtained with SW_m is characterized by an energy of 0.20 eV/atom relative to the crystalline phase. Measurements of enthalpy differences between amorphous and crystal silicon are about half this value [75, 76], but it has been recently shown that first-principles calculations were required to reach a better agreement between simulations and experiments [77]. Nevertheless, our computed energy difference is among the lowest reported in investigations using classical potentials [19, 78, 79].

It is also known that experimentally prepared a-Si samples are about 2% less dense than the crystalline phase [74, 75], a property which is not well reproduced by the original SW potential, which yields an a-Si sample more dense by 3%. Despite being specifically designed for a-Si modelling, the VBWM tends to produce an amorphous model with a density 6% lower than the crystal [22]. Using our modified potential, we found instead that our a-Si sample is 2% less dense, in excellent agreement with experiments.

Although a good description of the amorphous phase of silicon was not a primary goal during our investigations, it clearly appears that the SW_m potential is well suited for this task. It corrects the spurious behaviour of the original SW potential, like VBWM, but in addition it also provides an a-Si model that is in excellent agreement with experiments regarding density and energy.