Simulation methodology

First-principles DFT calculations were performed using projector-augmented wave pseudopotentials as implemented in the electronic structure Vienna Ab initio Simulation Package^65,66. Gamma-surface calculations were performed using 6-quintuplet slabs, with changes in total energy obtained after translating the upper 3 quintuplets parallel to the hexagonal basal plane. The translation vector t sampled the lateral area of the conventional unit cell on a uniform 10 × 10 mesh. After each translation, the atomic positions were relaxed in the direction normal to the fault. The lateral lattice vectors (in Cartesian coordinates) a₁ = a(1, 0, 0) and a2=(a∕2)(3,1,0) and the initial atomic positions were obtained for the lattice parameter a computed by relaxing forces and stresses for bulk Bi₂Te₃. Convergence with respect to the Brillouin zone sampling was achieved employing uniform Monkhorst–Pack k-point meshes with sizes up to 7 × 7 × 1. A smooth function γ(t) was obtained by interpolating between the measured points using a multiquadric radial basis function. Several exchange-correlation functionals were employed, namely, LDA⁴⁰, optPBE-vdW and optB88-vdW^41–44, and DFT-D2⁴⁵. The lattice parameters and elastic constants obtained with these functionals are summarized in Supplementary Table ¹ and are in good agreement with previous DFT calculations⁶⁷.

In the SDVPN model used here, the disregistry function δ_i(x) was discretized on a mesh {x^α}_{α = 1,...,N} of N = 300 points with the spacing Δx=a3∕3. The total dislocation energy is E_total = E_elastic + E_misfit, where the elastic strain energy E_elastic depends on the energy coefficients K_ij, which in turn depend on the elastic constants C_ij. The misfit energy Emisfit= ∑α=1Nγδi(xα)Δx was obtained from the DFT gamma-surface. The equilibrium disregistry δ_i(x^α) was found by numerical minimization of E_total under appropriate boundary conditions. As the initial guess, we used the arctangent disregistry function predicted by the classical Peierls–Nabarro model. A continuous disregistry curve was obtained by smooth interpolation between the x^α points. The calculations were performed separately for each of the DFT functionals. Further technical details of the SDVPN calculations can be found in the Supplementary Note ¹.