Reactive potential for the study of phase-change materials: GeTe

Tersoff [4, 5] proposed a three-body short-range potential in which the total energy E is computed as the sum of terms V_ij, whose strength is controlled by the local coordination

$$\begin{equation} E=\frac{1}{2} \sum_{i \neq j} V_{ij}, \end{equation} \tag{ 2 }$$

in which V_ij is defined as follows:

$$\begin{equation} V_{ij} = f_{\rm c}(r_{ij}) [A_{IJ} \mathrm{e}^{-\lambda_{IJ} r_{ij}} - b_{ij}^{IJ}B_{IJ}\,{\rm e}^{-\mu_{IJ}r_{ij}}], \end{equation} \tag{ 3 }$$

$$\begin{equation} b_{ij}^{IJ} = \chi_{IJ} [1+(\beta_I \xi_{ij}^{IJ})^{n_I}]^{\frac{1}{2n_I}}, \end{equation} \tag{ 4 }$$

$$\begin{equation} \xi_{ij}^{IJ} = \sum_{k \neq i,j} f_{\rm c}(r_{ik}) \,{\rm e}_{ijk}^{IJK}t_{ijk}^I, \end{equation} \tag{ 5 }$$

$$\begin{equation} t_{ijk}^I = 1 + \frac{c_I^2}{d_I^2} - \frac{c_I^2}{d_I^2 + [h_I -\cos(\theta_{ijk})]^2}. \end{equation} \tag{ 6 }$$

The indices i and j define the pair of atoms, I and J the atomic species, r_ij is the distance between the atoms i and j, f_c(r_ij) is a continuous function of r_ij, such that f_c(r_ij) = 1 for r_ij ⩽ R_IJ, f_c(r_ij) = 0 for r_ij ⩾ S_IJ and $f_{\mathrm { c}}(r_{ij})=\frac {1}{2} [1+\cos (\pi \frac {r_{ij}-R_{IJ}}{S_{IJ}-R_{IJ}})]$ for values of r_ij between R_IJ and S_IJ, with R_IJ < S_IJ. The parameters A_IJ, B_IJ, λ_IJ and μ_IJ determine the shape of the Morse potential V_ij and depend on both species types. The function b^IJ_ij controls the attractive term, which is reduced as the atomic coordination increases. Different variants of the Tersoff potential have been proposed; here for the term e^IJK_ijk we used the following definition, in which

$$\begin{equation} e_{ijk}^{IJK}=\exp (\mu_{IJ} r_{ij} - \mu_{IK} r_{ik} ). \end{equation} \tag{ 7 }$$

The term t^I_ijk describe the bending energy, which is minimized when the cosine of the angle θ_ijk, with i being the vertex, approaches the value h_I. The two parameters c_I and d_I control the strength of the penalty angular term.

Billeter et al [6] introduced two additional terms in the formulation of the total energy

$$\begin{equation} E=\frac{1}{2}\sum_{i \neq j} V_{ij} + \sum_I N_I E_I^0 + \sum_i E_i^{\rm c}, \end{equation} \tag{ 8 }$$

where the first term is that of equation (2). The second term in equation (8) is the sum of the core energies E⁰_I of each species. This term was introduced to reproduce the energy of configurations containing different numbers of atoms. The last term introduces an additional penalty energy E^c_i to improve the description when there are changes in coordination, such as formation of defects, interfaces and surfaces. E^c_i is defined as

$$\begin{equation} E_i^{\rm c}=c_{I,1} \Delta z_i + c_{I,2} \Delta z_i^2, \end{equation} \tag{ 9 }$$

where Δz²_i is the difference between the current coordination of atom i of type I, z_i and its typical coordination z⁰_I. The coordination z_i is defined using the functions f^IJ_s,ij and b^IJ_ij defined before via $z_i=\sum _{j \neq i} f_{s,ij}^{IJ}b_{ij}^{IJ}$. The deviation Δz_i from the typical coordination is defined via

$$\begin{eqnarray} \fl f_{\rm s}(z_i)={\rm int}(|z_i-z_0^I|) +\left\{ \begin{array}{@{}ll@{}} 0 & {\rm for}\ z\leqslant z_{\mathrm{ T}}-z_{\mathrm{ B}},\\ \frac{1}{2}+\frac{1}{2}\sin(\frac{\pi}{2}\frac{z-z_{\rm T}}{z_{\rm B}}) & {\rm for} \ z_{\mathrm{ T}}-z_{\mathrm{ B}} \leqslant z \leqslant z_{\mathrm{ T}}+z_{\mathrm{ B}},\\ 1 & {\rm for} \ z \geqslant z_{\mathrm{ T}}+z_{\mathrm{ B}}, \end{array}\right. \end{eqnarray} \tag{ 10 }$$

where z = |z_i − z^I₀| − int(|z_i − z^I₀|); Δz²_i is set constant for values of z larger than z_cut, see figure 1; z_T, z_B and z_cut are three adimensional, species-independent parameters.

Zoom In Zoom Out Reset image size

In the present work we used the augmented potential of [6] but modified the definition of the angular term t^I_ijk by choosing different values of the parameters d_I,J,K and h_I,J,K, which now depend on the three species I, J and K forming the angle. Moreover, we added an additional parameter θ_0,I for each species I. The term in t^I_ijk is now defined as

$$\begin{equation} t_{ijk}^I = 1 + \frac{c_I^2}{d_{I,J,K}^2} - \frac{c_I^2}{d_{I,J,K}^2 + [h_{I,J,K} -\cos(\pi \frac{\theta_{ijk}-\theta_{0,I}}{\pi-\theta_{0,I}})]^2}. \end{equation} \tag{ 11 }$$

The modification in equation (11) allowed a better description of the angle distribution in the amorphous phase, see section 3.

The potentials were constructed following the recipe of [6]. To find the optimal set of parameters A_IJ, B_IJ, λ_IJ, μ_IJ, β_I, n_I, c_I, d_I,J,K, h_I,J,K, R_IJ, S_IJ, E₀, Z₀, z_B, z_T, c_I,1, c_I,2 and θ_0,I that minimize the energy and the force mismatch computed via FP and classical potential. The energy mismatch Δ²_E was defined as

$$\begin{equation} \Delta_{E}^2=\frac{1}{N_{f}}\sum_{n=1}^{N_{f}}\left( E_{{\rm cla},n}-E_{{\rm DFT},n} \right)^2, \end{equation} \tag{ 12 }$$

where the summation is carried out over the total number of configurations, N_f. E_cla,n and E_DFT,n are the energies of the n-configuration computed with the classic potential and DFT, respectively. Δ²_E is the mean square deviation of the errors of all frames. Similarly, it is possible to define the force mismatch Δ²_F as

$$\begin{eqnarray} &&\Delta_{F}^2 =\frac{1}{3 \sum_{n=1}^{N_{{f}}} N_{{\rm at},n}} \sum_{n=1}^{N_{f}} \sum_{i=1}^{N_{{\rm at},n}} \sum_{\alpha=x,y,z}\left( F_{{\rm cla},n,i,\alpha}-F_{{\rm DFT},n,i,\alpha} \right)^2, \end{eqnarray} \tag{ 13 }$$

where F_{cla(DFT),n,i,α} represents the classical (DFT) force of atom i along the Cartesian direction α of the n-configuration and N_at,n is the total number of atoms of the n-system.

In [6] the energy mismatch and the force error were optimized simultaneously by minimizing Δ²_W defined as

$$\begin{equation} \Delta_W^2 = w_{E} \Delta_{E}^2+ (1-w_{E}) r_{\rm B} \Delta_{F}^2, \end{equation} \tag{ 14 }$$

where r_B is the Bohr radius and w_E is a weight parameter varying between 0 and 1. By setting w_E = 1, Δ²_W is equal to Δ²_E, which means that only the energy mismatch is minimized; in contrast, w_E = 0 correspond to minimizing the force error. Another way to measure the goodness of the potential is to compare the directions of the classical and DFT forces to measure their alignment. This is done computing the correlation coefficient C_F which defined as follows. The best alignment correspond to C_F  = 1.

$$\begin{eqnarray} &&C_{F}= \frac{F_{\rm scalar}}{\sqrt{F_{\rm cla}}\sqrt{F_{\rm DFT}}}, \end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray} {F}_{\rm scalar}&=&\sum_{n=1}^{N_{f}} \sum_{i=1}^{N_{{\rm at},n}} \sum_{\alpha=x,y,z} F_{{\rm cla},n,i,\alpha} F_{{\rm DFT},n,i,\alpha},\nonumber \\ {F}_{\rm cla}&=&\sum_{n=1}^{N_{f}} \sum_{i=1}^{N_{{\rm at},n}} \sum_{\alpha=x,y,z} F_{{\rm cla},n,i,\alpha}^2,\nonumber \\ {F}_{\rm DFT}&=&\sum_{n=1}^{N_{f}} \sum_{i=1}^{N_{{\rm at},n}} \sum_{\alpha=x,y,z} F_{{\rm DFT},n,i,\alpha}^2. \end{eqnarray} \tag{ 16 }$$

We started by fitting the forces of isolated Ge–Ge, Ge–Te and Te–Te dimers to obtain the starting coefficients A, B, λ, μ and an initial estimate of the cutoff distances R and S. Subsequently, groups of parameters were fitted simultaneously. We cycled hundreds of times over all groups to maximize the correlation coefficient C_F ; the coefficients c_I,1 and c_I,2 of the additional term of equation (9) were fixed to zero. Finally, during the refinement, we kept fixed R and S. We would like to point out that a larger values of S results in a higher correlation coefficient, but the drawback is that the potential becomes less transferable over a wide range of volumes, because at different densities different numbers of atoms can be included in the cutoff range. The parameters defining the term of equation (9) were optimized using a weight w_E between 0.01 and 0.02 using a subset of crystalline structures at different volumes. The optimal parameters are reported in table 1.

We produced 1000 configurations to build up the database with the reference structures using different approaches: (i) most of the structures were obtained via FP MD simulations at different temperatures; (ii) other structures were produced artificially from c-GeTe either by displacing selected groups of atoms along high-symmetry directions, by removing selected atoms to produce local defects or by randomly displacing the atomic positions; and (iii) in the final steps of the construction of the potential, additional FP calculations of energies and forces were carried out on structures produced with our potential via standard MD and replica-exchange (RE) [7, 8] MD simulations.

Our study employs both the Car–Parrinello scheme (CP) [9] and Born–Oppenheimer (BO) FP MD to simulate liquid, amorphous and crystalline GeTe. The FP simulations were performed within the framework of DFT in the local density approximation supplemented by generalized-gradient corrections [10], as implemented in the CPMD code². Local spin-density calculations were carried out to test the presence of spin polarization on selected snapshots. We never found unpaired electrons in these tests. We employed Goedecker–Teter–Hutter pseudopotentials [11] with a plane-wave expansion of the Kohn–Sham orbitals up to a kinetic energy cutoff of 80 Ry. We used only the γ point to sample the Brillouin zone (BZ). GeTe structures were modeled using three-dimensional periodic boundary conditions at different densities; the box size and the number of atoms for each system are reported in table 2.

We used the PW code, which is part of the Quantum-ESPRESSO package [12], only to compute two equations of state (EOS) of c-GeTe, using the same pseudopotentials and the same setup as described above except for the sampling of the BZ. EOS were used to test our FP simulation setup versus experimental and theoretical data available from literature and to validate the capability of the classical potential to reproduce the energies and the structures over a wide range of volumes, see section 3. GeTe has two crystalline phases: the rhombohedral phase, also called α-GeTe, which has R3m symmetry and is stable at standard conditions, and the cubic one, β-GeTe, which has Fm3m symmetry and is stable at 432 °C [13–15]. The α-GeTe phase was modeled using a single GeTe unit in a rhombohedral cell; geometry optimizations were repeated at several volumes and the resulting energies were fitted via a Murnaghan EOS [16]. At each volume, the edge a and the angle α of the rombohedric box were optimized such as to reproduce a diagonal stress tensor with residual anisotropy in the diagonal components of less than 1 kbar; BZ integration was performed over a 12 × 12 × 12 Monkhorst–Pack mesh [17]. The optimized lattice parameter of the box was a = 4.410 Å and the angle α was equal to 57.57°. These values are close to the experimental values: a = 4.2810 Å and α = 58.358° [18]. To model β-GeTe, we used a cubic box containing four GeTe units. For this system, BZ integration was performed over a 12 × 12 × 12 Monkhorst–Pack mesh. The optimized lattice parameter was a = 6.060 Å  to be compared with the experimental one of a = 5.999 Å [13, 14]. The equation of state of pure Ge was computed using a cubic box containing eight atoms and BZ integration was performed over a 12 × 12 × 12 Monkhorst-Pack mesh: in this case the optimized lattice parameter was a = 5.775 Å, to be compared with the experimental one of 5.658 Å [19]. The equilibrium lattice constant of Te was obtained using an hexagonal box containing three atoms, and BZ integration was performed over a 6 × 6 × 6 Monkhorst–Pack mesh: here the optimized lattice parameters were a = 4.497 Å and c/a = 1.328.

The majority of the structures was produced via FP MD simulations using two different systems: the smaller one contained 64 atoms, the larger one 216 atoms. We used a rombohedric box containing 128 and 127 atoms to simulate perfect α-GeTe and α-GeTe containing a germanium vacancy, respectively, structures (a)–(b) in table 2. With the 64-atom system containing 32 GeTe units, we simulated liquid GeTe via BO MD simulations. Liquid GeTe was modeled with a 12.4-Å cubic box corresponding to the experimental density (5.58 g cm⁻³) of the l-GeTe [14, 20], structures (c) in table 2. The temperature was controlled via a Nosé–Hoover (NH) thermostat [21, 22] with a target temperature of 1200 K. The ionic timestep in BO MD was always 2.4 fs, and the total duration of the MD simulation was 17.0 ps. We reduced the temperature to 1000 K and performed 8.4 ps, followed by 9.7 ps at T = 800 K, and 10.2 ps at 700 K. Positions, forces and energies were recorded periodically. We scaled the atomic coordinates obtained at liquid density at T = 1200 K in a simulation box of 12.137 Å, corresponding to the density 5.95 g cm⁻³, and generated other structures (d) via a 6.1 ps long BO MD simulation at T = 1200 K. We cooled the system from T = 1200–300 K via four independent MD simulations lasting 4.4 (e), 14.2 (f), 17.7 and 21.3 ps (h) by rescaling the velocity at each step with a factor between 0.9992 and 0.9998.

Other structures were generated using a larger system containing 108 GeTe units. Structures labeled (i)–(l) in table 2 were created by introducing up to four germanium vacancies, and (m) and (n) by translating either only Te atoms or the Ge and Te atoms along the [111] direction. This particular distortion was inspired by the work of Kobolov et al [23] in which it was shown that changes in the optical properties in GST can occur during the distortion of the crystalline structure mentioned above. To build a potential that is transferable to a broad range of densities, we randomized the positions of c-GeTe with different displacements (always lower than 0.5 Å) at different lattice parameters, see structures (o) in table 2. With the 216-atom system, the atomic trajectories were generated via CP MD simulations at finite temperature; the configurations were then post-processed to compute their ground-state energies and the forces to be used for fitting (training) and testing the classical potential. This two-steps approach was necessary to correct the small error that affects energies and forces when using the CP scheme at finite temperature [24]. In contrast, BO MD simulations were used to evolve the 64-atom system at different temperatures. In the CP MD simulations we used a 0.12 fs timestep and an electron mass equal to 600 au; but to integrate the equation of motions of the ions in BO MD simulations a larger timestep equal to 2.4 fs was used. Constant temperature was imposed on the ions by a massive NH thermostat [21, 22]. As in GeTe at high temperature the energy gap between the highest occupied and lowest unoccupied orbitals was less than 0.1 eV, adiabaticity of the CP MD trajectory was maintained by coupling two separate NH thermostats to the nuclear and electronic subsystems [25]. We melted c-GeTe in a simulation lasting 5.5 ps by heating the system to T = 1200 K using a cubic box with 17.997 Å long edges, structures (p) in table 2. In a larger box, corresponding to a density equal of 5.58 g cm⁻³, see structures (q), l-GeTe was simulated at T = 1200 K in a 9.2 ps long run. Liquid GeTe was then scaled at higher density, 5.95 g cm⁻³, and quenched from 1200 K to room temperature in a 11.9 ps MD. DFT forces and energies of additional a-GeTe structures obtained via RE MD classical simulations were computed and using to refine the potential (s). We performed a CP MD simulation lasting 18 ps at T = 300 − 400 K and at density equal to 5.65 g cm⁻³ of structure (t), which is particularly favorable in energy obtained via classical RE simulations. Structures (u) result from DFT static calculations of a-GeTe structures obtained via classical MD simulations quenching from 875 K to room temperature in 4 ns.

We also used configurations of pure germanium and pure tellurium. Structures containing 512 Ge-atoms in a cubic box were produced via classical MD and RE simulations using the potential proposed by Tersoff [5] at different temperatures between 300 and 3000 K; the box size was 21.537 Å  (v). Forces and energies on selected snapshots were computed via DFT as described above. To improve the fit of the parameters describing the tellurium we performed short DFT MD simulations lasting 1.5 ps at T = 625 K using a box containing 343 Te-atoms at a density of 5.68 g cm⁻³ (w).

We computed the EOS of the rhombohedral α and cubic β phases of c-GeTe using two systems made of 64 and 108 GeTe units, respectively. DFT and classical EOS are compared in figure 2. We fitted the energy versus volume data via Murnaghan EOS (see section 2). The lattice parameters, the bulk modulus and the bond lengths at the equilibrium volume are reported in table 3. The optimized structures of the rombohedric and the cubic phases at the equilibrium volume are shown in figure 3. We also tested the capability of the classical potential to reproduce the Peierls distortions occurring in c-GeTe over a broad range of volumes using the β phase, which is metastable at standard conditions. This bimodal distribution of short and long bonds has been shown even at high temperature by Fons et al [26]. We optimized the geometries with the classical potential at different lattice parameters, see figure 4, and we compared the DFT and the classical distances of the short and the long bonds. Both the DFT and the classical calculations used the same setup as was employed to compute the EOS: the 4 GeTe and 108 GeTe cubic boxes for the DFT and classical simulations, respectively. The Peierls distortion is reproduced correctly over the entire volume range. A very good agreement between DFT and classical distances is achieved at the technologically most interesting densities, which are close to the amorphous and the equilibrium crystalline volume per atom. We computed the formation energy of a Ge vacancy in the rhombohedral (cubic) phase by removal of a Ge atom from the perfect crystal containing 128 (216) atoms. The energy is computed after relaxation of the structure and using as reference the energy of a Ge atom in its bulk. The formation energies are 0.00 and −0.22 eV for the rhombohedral and cubic phases, respectively. The corresponding DFT values are 0.93 and 0.62 eV for the rhombohedral and cubic phases, respectively. The error is about 0.9 eV for both the geometries. For further validation, we performed the test proposed in [3], namely, to compute the energies of α-GeTe structures in which the Ge–Te distances were changed by shifting the Ge-atoms along the [111] direction, see figure 5. The agreement between DFT and classical data is very good.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To compute the EOS of crystalline Ge we used a cubic 512-atom box. The equilibrium lattice parameters (bulk modulus, B₀) obtained via the fit of the EOS for pure Germanium are 5.775 Å (B₀ = 57.3 GPa) and 5.853 Å (B₀ = 46.0 GPa), DFT and classical, respectively. The equilibrium geometry of pure tellurium was computed using an hexagonal box containing 375 Te atoms, a 5 × 5 × 5 supercell. The equilibrium volume is 34.607 Å³ atom⁻¹, the c/a = 1.220, and the bulk modulus B₀ = 115 GPa, to be compared with 34.880 Å³ atom⁻¹, c/a = 1.328 and B₀ = 21.5 GPa obtained via our DFT calculations. These large deviations show the limit of the short of the Te–Te interaction, about 4 Å, too short to properly model pure Tellurium. In the current formulation, the maximum value of the cutoff cannot be larger than the second neighbor distance of Te–Te in c-GeTe. A larger database of structures including liquid and amorphous, like those produced in [27, 28], combined with a long range interaction can improve the description of pure Tellurium.

The structural and dynamic properties of l-GeTe were investigated using two cubic boxes containing 216 and 4096 atoms at a density equal to 5.65 g cm⁻³. We performed two constant temperature MD simulations lasting 100 ps at T = 1150 K using the Berendsen thermostat with a 24 fs coupling constant. We selected different configurations from the two trajectories corresponding to the small and the large systems, and from those configurations we started 100 and 32 independent MD simulations at T = 1150 K using 216 and 4096 atoms, respectively. These simulations were later used to quench the liquid in the amorphous phase. We computed the pair correlation function and the angle distribution from the 100 trajectories containing 216 atoms; the same was done with the 32 trajectories containing 4096 atoms. Coordination numbers are reported in the supplementary data (available from stacks.iop.org/NJP/15/123006/mmedia). In figures 6 and 7, DFT results taken from [3] and classical results on the 216 systems are compared. The overall agreement is good, with the largest deviation between DFT and classical results appearing in the pair correlation function: we obtained a sharper and higher Ge–Te peak. We plotted only the results of the smaller system, because we did not find appreciable differences between the properties computed for the smaller (216 atoms) and the bigger (4096 atoms) systems. This is probably due to the fact that the interactions of our potential are much shorter than the size of the 216-atom cubic box (which s 18.6 Å in size). We equilibrated the 4096-atom system at T = 1000 K with an MD simulation at the liquid experimental density, 5.58 g cm⁻³ [20], then we performed a 180 ps constant volume and energy (NVE) MD simulation to compute the diffusion coefficient from the mean-squared displacement using the Einstein equation

$$\begin{equation} D_I = \lim_{t\to\infty} \frac{\langle \vert R_{I,i}(t)-R_{I,i}(0) \vert^2 \rangle_{i \in I} }{6t}, \end{equation} \tag{ 17 }$$

where t is the time. The average is carried out over the i atoms of the I-species. R_I(t) are the time dependent coordinates. The diffusion coefficients computed are 5.58 and 4.11 ×10⁻⁵cm² s⁻¹ for the germanium and tellurium species, respectively. Our short-range potential overestimates the theoretical diffusion coefficient computed via FP MD in [29] by about 20 (Ge) and 5% (Te): 4.65 ×10⁻⁵ (Ge) and 3.93 × 10⁻⁵ cm² s⁻¹ (Te). Probably, our short-range Tersoff potential produces a smoother potential that results in a more diffusive GeTe. The effect of the cutoff on the diffusion coefficient has been mentioned in [3], in which the authors reported that an increase of the cutoff from 6.0 to 6.88 Å resulted in an improved description of the diffusion coefficient.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Amorphous GeTe configurations were generated by quenching l-GeTe via a series of classical MD simulations, in each of which the temperature was reduced in subsequent 50 K steps from 1150 K to room temperature. The entire duration of the simulation was 100 ps. We chose the same quenching time as in [3] to compare structures obtained with the same protocol. We performed 100 and 32 independent quenching using different initial configurations of l-GeTe containing 216 and 4096 atoms, respectively, i.e. systems described in section 3.2. The box size corresponded to a of 5.65 g cm⁻³, which is close to the experimental values of a-GeTe: 5.6 ± 0.6 and 5.7 g cm⁻³ in [30, 31], respectively. Using each of the 100 a-GeTe structures as starting point, we then carried out 100 MD simulations at T = 300 K. We computed the pair correlation function and the angle distribution averaging over the 100 trajectories. In figure 8, our classical pair correlation function is compared with DFT data from [3]. The Ge–Ge pair correlation functions are in good agreement; the first Ge–Te presents the largest deviation, the classic peak is much narrower, and a shift of about 0.2 Å is observed in the Te–Te peaks. The deviation in the Te–Te interaction can be due to the short-range potential. We could have choosen a longer cutoff parameter S_I,J for Te–Te, but this resulted in a worse description of c-GeTe at high densities, where the second neighbors Te–Te are at shorter distances. The classical and DFT angle distributions are compared in figure 9: the overall agreement is good. Our classical angle distribution compares well with the DFT one. As reported in [3], we also found the formation of strained three-membered rings (3MR), which are responsible for the small shoulder in the angle distribution at about 60°. These highly strained rings are probably absent in the experimental a-GeTe samples, which are usually annealed at a much slower rate than in our simulations. Possibly they would not appear either in DFT if longer simulations were carried out. Our potential produces a much smaller fraction of 3MR than the classical potential proposed in [3] does. We would like to emphasis that the modified term described in equation (11) improved the angle distribution. Specifically, the second peak in the bimodal angle distribution of germanium at about 165° did not appear with the original equation of the Tersoff potential. We computed the vibrational density-of-states (DOS) of the 100 optimized a-GeTe geometries and compared the classical DOS with the DFT [3], see figure 10. The averaged potential energy of the optimized amorphous structures is 0.27 eV atom⁻¹ higher than that energy of the crystalline phase at the equilibrium volume.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The structural properties computed so far do not require a definition of bonds, with the exception of the angle distribution which does not depend dramatically on a particular definition. On the contrary, the analysis of coordination numbers type resolved, q-order parameter [32], and statistic of rings [33] depends strongly on the criteria used to define a bond. We report in the supplementary data (available from stacks.iop.org/NJP/15/123006/mmedia) additional structural analysis to provide data for comparing with pre-existing data using the definition of bonds of [3]; we find this criteria not enough accurate to define properly the bond network. For this reason, although this paper is focused on the classical potential, in the supplementary data, we show that localized Wannier functions, just to mention one, are more appropriate to define the bond network than criteria based on distances.