Environment-dependent Interatomic Potential Research Articles

Comparison and Assessment of Different Interatomic Potentials for Simulation of Silicon Carbide.

Interatomic potentials play a crucial role in the molecular dynamics (MD) simulation of silicon carbide (SiC). However, the ability of interatomic potentials to accurately describe certain physical properties of SiC has yet to be confirmed, particularly for hexagonal SiC. In this study, the mechanical, thermal, and defect properties of four SiC structures (3C-, 2H-, 4H-, and 6H-SiC) have been calculated with multiple interatomic potentials using the MD method, and then compared with the results obtained from density functional theory and experiments to assess the descriptive capabilities of these interatomic potentials. The results indicate that the T05 potential is suitable for describing the elastic constant and modulus of SiC. Thermal calculations show that the Vashishta, environment-dependent interatomic potential (EDIP), and modified embedded atom method (MEAM) potentials effectively describe the vibrational properties of SiC, and the T90 potential provides a better description of the thermal conductivity of SiC. The EDIP potential has a significant advantage in describing point defect formation energy in hexagonal SiC, and the GW potential is suitable for describing vacancy migration in hexagonal SiC. Furthermore, the T90 and T94 potentials can effectively predict the surface energies of the three low-index surfaces of 3C-SiC, and the Vashishta potential exhibits excellent capabilities in describing stacking fault properties in SiC. This work will be helpful for selecting a potential for SiC simulations.

Transferability of interatomic potentials for germanene (2D germanium)

The capacities of various interatomic potentials available for elemental germanium, with the scope to choose the potential suitable for the modeling of germanene (2D germanium) allotropes,f were investigated. Structural and mechanical properties of the flat, low-buckled, trigonal dumbbell, and large honeycomb dumbbell single-layer germanium (germanene) phases, were obtained using the density functional theory and molecular statics computations with Tersoff, modified embedded atom method, Stillinger–Weber, environment-dependent interatomic potential, ReaxFF, and machine-learning-based interatomic potentials. A systematic quantitative comparative study and discussion of the findings are given.

Comparing empirical interatomic potentials to modeling silicon surface stress

This study compares six widely used classical many-body potentials – Tersoff-T3, Stillinger–Weber, modified Tersoff, optimized modified Tersoff, environment-dependent interatomic potential and modified embedded atom method – to model silicon surface stress. Functional forms for all potentials are first analyzed to determine atomic stress. This is then followed by a set of surface stress analysis on silicon {100} and {110} surfaces. A significant difference between those potentials is observed to calculate the surface stress in silicon crystal. The Tersoff-T3 and environment-dependent interatomic potentials demonstrate a significant surface stress in unreconstructed silicon {100} and {110} surfaces, while the Stillinger–Weber potential estimates a negligible surface stress. The implication of the surface stress is studied on the mechanical response of silicon nanowires. A dramatic contrast between potentials is recognized to calculate the equilibrium elongation due to the surface stress in silicon nanowires along <100> crystal orientation. While the environment-dependent interatomic potential predicts a contraction in the equilibrium elongation in relaxed nanowires, the Tersoff-T3 potential exhibits a considerable uniaxial expansion in silicon nanowires comparable with density functional theories. Findings in this work give more insight into the capability of empirical interatomic potentials to capture the surface stress in silicon nanostructures.

Thermal transport at 6H-SiC/graphene buffer layer/GaN heterogeneous interface

We investigate the epitaxial growth of graphene on SiC and thermal transport at 6H-SiC/graphene buffer layer/GaN heterogeneous interface. Under Tersoff-Erhart-Albe potential (TEA) and environment-dependent interatomic potential (EDIP), graphene buffer layer (GBL) is epitaxially grown by simulated annealing method, and simulations show that TEA potential is superior to EDIP potential by evaluation of nucleation temperature, radial distribution function, average atomic binding energy, atomic structure and bond length. Further, by thermal relaxation method, single-vacancy defects in graphene can reduce interface thermal resistance (ITR) at 6H-SiC/GBL/GaN heterogeneous interface with maximum reduction of 9.6%, meaning that ITR can be reduced by defect engineering. In addition, with the increase of nucleation areas, especially the bonding of different nucleation zones, ITR also has a highest decline of 17.9%, and energy is easier to pass through the interface, illustrating the importance of graphene integrity for interfacial heat transfer. At last, the increase in temperature (300 K to 1100 K) results in a maximum 19.6% decrease in ITR, indicating that ITR is sensitive to temperature. The selection of potential and the analysis of the heat transfer of the three-materials heterogeneous interface contribute to thermal design of graphene-based semiconductor devices.

Thermal conductivity of strained silicon: Molecular dynamics insight and kinetic theory approach

In this work, we investigated the tensile and compression forces effect on the thermal conductivity of silicon. We used the equilibrium molecular dynamics approach for the evaluation of thermal conductivity considering different interatomic potentials. More specifically, we tested Stillinger-Weber, Tersoff, Environment-Dependent Interatomic Potential, and Modified Embedded Atom Method potentials for the description of silicon atom motion under different strain and temperature conditions. It was shown that the Tersoff potential gives a correct trend of the thermal conductivity with the hydrostatic strain, while other potentials fail, especially when the compression strain is applied. Additionally, we extracted phonon density of states and dispersion curves from molecular dynamics simulations. These data were used for direct calculations of the thermal conductivity considering the kinetic theory approach. Comparison of molecular dynamics and kinetic theory simulations results as a function of strain and temperature allowed us to investigate the different factors affecting the thermal conductivity of the strained silicon.

On the transferability of interaction potentials for condensed phases of silicon

A systematic comparison of structural, thermodynamical and kinetic properties of silicon by molecular dynamics simulation was done using five interatomic potentials, namely the Stillinger-Weber (SW) model, the environment-dependent interatomic potential (EDIP), two versions of the modified embedded-atom method (MEAM-1 and MEAM-2) and the angular dependent potential (ADP). This benchmark allows us to determine which interatomic potential is the most adequate to simulate this material. Structural properties and the vibrational density of states (VDOS) for diamond-structure crystal and amorphous phases were studied. We also obtained the diffusion coefficient for the liquid phase. For the models able to describe crystallization from the melt, the interfacial free-energy was also computed. Results show that the SW potential outperforms the remaining models in several aspects including the structure factor and the coordination number distribution of the liquid and the melting temperature. MEAM-2 follows as the second best model due to its ability to reproduce the crystal-liquid interfacial free-energy and the structure factor and the VDOS for the amorphous phase. EDIP, MEAM-1 and ADP models can reproduce some physical properties but they lack transferability, which is a key aspect of a good interatomic potential.

Thermal conductivity of silicon using reverse non-equilibrium molecular dynamics

Simulations are performed using the reverse non-equilibrium molecular dynamics (rNEMD) method and the Stillinger-Weber (SW) potential to determine the input parameters for achieving ±1% convergence of the calculated thermal conductivity of silicon. These parameters are then used to investigate the effects of the interatomic potentials of SW, Tersoff II, Environment Dependent Interatomic Potential (EDIP), Second Nearest Neighbor, Modified Embedded-Atom Method (MEAM), and Highly Optimized Empirical Potential MEAM on determining the bulk thermal conductivity as a function of temperature (400–1000 K). At temperatures &amp;gt; 400 K, data collection and swap periods of 15 ns and 150 fs, system size ≥6 × 6 UC2 and system lengths ≥192 UC are adequate for ±1% convergence with all potentials, regardless of the time step size (0.1–0.5 fs). This is also true at 400 K, except for the SW potential, which requires a data collection period ≥30 ns. The calculated bulk thermal conductivities using the rNEMD method and the EDIP potential are close to, but lower than experimental values. The 10% difference at 400 K increases gradually to 20% at 1000 K.

Mechanical response of bilayer silicene nanoribbons under uniaxial tension.

Understanding the behaviour of nanoscale systems is of great importance to tailor their properties. To this aim, we investigate the Young's modulus (YM) of defect-free and defective armchair bilayer silicene nanoribbons (SNRs), at room temperature, as a function of length and distance between layers. In this study, we perform molecular dynamics simulations using the environment-dependent interatomic potential to describe the interaction of the Si atoms. We show that the Young's modulus of pristine and defective bilayer SNRs increases with the ribbon length exhibiting size dependence. In general, YM of defective bilayer SNRs is smaller than the value obtained for the defect-free case, as a result of the number of missing bonds. In all cases, as the interlayer distance increases YM decreases and the buckling increases. It is shown that the YM exhibits a quadratic interlayer distance dependence. Finally, when only one layer has a mono-vacancy defect, the atomic stress distribution of the pristine layer is affected by the presence of the vacancy. This effect can be considered as a “ghost vacancy” since the deterioration of the pristine layer is similar to that shown by the defective one. These results show that YM of pristine and defective bilayer SNRs could be tailored for a given length and interlayer distance. It is also found that the fracture stress and the fracture strain of defective bilayers are both smaller than those obtained for the defect-free ones.

Characterization of silicon surface elastic constants based on different interatomic potentials

Mechanical properties of materials are an important factor in designing nanoscale systems. Several researches and experiments have shown that the mechanical properties of the nano-scale materials are different from those of bulk. One of the major reasons for this difference is that the ratio of surface to volume increases at the nano-scale, and the effects of free surfaces become very important. In this paper, we have measured the surface elastic constants of silicon crystalline structure using different interatomic potentials. The potentials employed here are EDIP (Environment-Dependent Interatomic Potential), Stillinger-Weber and Tersoff, and also different crystalline orientations are examined. In order to obtain surface elastic constants, we have first determined the lattice constant resulting in minimum energy of bulk model for each potential, then we have measured the bulk and film elastic constants using respective models. It is observed that by increasing film thickness, elastic constants approach to those of bulk. Next, using the results of the previous sections, the surface elastic constants are measured and discussed for different surfaces using different potentials.

A KIM-compliant potfit for fitting sloppy interatomic potentials: application to the EDIP model for silicon

Fitted interatomic potentials are widely used in atomistic simulations thanks to their ability to compute the energy and forces on atoms quickly. However, the simulation results crucially depend on the quality of the potential being used. Force matching is a method aimed at constructing reliable and transferable interatomic potentials by matching the forces computed by the potential as closely as possible, with those obtained from first principles calculations. The potfit program is an implementation of the force-matching method that optimizes the potential parameters using a global minimization algorithm followed by a local minimization polish. We extended potfit in two ways. First, we adapted the code to be compliant with the KIM Application Programming Interface (API) standard (part of the Knowledgebase of Interatomic Models project). This makes it possible to use potfit to fit many KIM potential models, not just those prebuilt into the potfit code. Second, we incorporated the geodesic Levenberg–Marquardt (LM) minimization algorithm into potfit as a new local minimization algorithm. The extended potfit was tested by generating a training set using the KIM environment-dependent interatomic potential (EDIP) model for silicon and using potfit to recover the potential parameters from different initial guesses. The results show that EDIP is a ‘sloppy model’ in the sense that its predictions are insensitive to some of its parameters, which makes fitting more difficult. We find that the geodesic LM algorithm is particularly efficient for this case. The extended potfit code is the first step in developing a KIM-based fitting framework for interatomic potentials for bulk and two-dimensional materials. The code is available for download via https://www.potfit.net.

Continuum models calibrated with atomistic simulations for the transverse vibrations of silicon nanowires

We have conducted Molecular Dynamics (MD) simulations with the Environment-Dependent Interatomic Potential (EDIP) to obtain the natural frequency of ultra-thin Silicon Nanowires (SiNWs) with various crystallographic structures, boundary conditions and dimensions. As expected, results show that the mechanical properties of SiNWs are size-/orientation-dependent. The observed phenomena are ascribed to the surface effects which become dominant due to the large surface-to-volume number of atoms at the investigated range of dimensions. Due to their accuracy, atomistic simulations are widely accepted for the investigations of such nano-scaled systems; however, they suffer from high computational costs. In this regard, continuum models are appreciated due to their efficacy; but they are not accurately applicable in very small nano-scales where the basic assumptions of continuity are undermined. Nevertheless, we aim in this paper to propose and prove that suitably improved continuum models are still applicable for these systems. In order to account for the surface effects, including the surface elasticity and surface stress, we have proposed a shell-core composite continuum model; at which the bulk core is modeled by the Euler–Bernoulli or Timoshenko beams, while the thin bonding surface layers with different elastic parameters are added to incorporate with the surface effects. We have estimated the surface stress and the surface elastic parameters using our MD simulations results; then we have shown that the calibrated continuum models are capable of predicting the vibrational behavior of SiNWs having thicknesses larger than 2 nm.

Threshold displacement energy and damage function in graphite from molecular dynamics

The environment dependent interatomic potential (EDIP) including Ziegler–Biersack–Littermark (ZBL) interactions for close encounters is applied to cascades starting from a host atom and from an interstitial atom. We find the room temperature displacement threshold to be 25 eV, increasing to 30 eV at 900 K. The latter correlates well with the measured threshold for vacancy production. Additionally, divacancy production is found to occur, including interlayer divacancies from around 60 eV. The data suggest a new, continuous damage function applies, where the threshold region depends on the square root of the primary knock-on atom (PKA) energy in excess of the threshold, evolving to a linear dependence on PKA energy.

On the validity of empirical potentials for simulating radiation damage in graphite: a benchmark

In this work, the ability of methods based on empirical potentials to simulate the effects of radiation damage in graphite is examined by comparing results for point defects, found using ab initio calculations based on density functional theory (DFT), with those given by two state of the art potentials: the Environment-Dependent Interatomic Potential (EDIP) and the Adaptive Intermolecular Reactive Empirical Bond Order potential (AIREBO). Formation energies for the interstitial, the vacancy and the Stone–Wales (5775) defect are all reasonably close to DFT values. Both EDIP and AIREBO can thus be suitable for the prompt defects in a cascade, for example. Both potentials suffer from arefacts. One is the pinch defect, where two α-atoms adopt a fourfold-coordinated sp3 configuration, that forms a cross-link between neighbouring graphene sheets. Another, for AIREBO only, is that its ground state vacancy structure is close to the transition state found by DFT for migration. The EDIP fails to reproduce the ground state self-interstitial structure given by DFT, but has nearly the same formation energy. Also, for both potentials, the energy barriers that control diffusion and the evolution of a damage cascade, are not well reproduced. In particular the EDIP gives a barrier to removal of the Stone–Wales defect as 0.9 eV against DFT's 4.5 eV. The suite of defect structures used is provided as supplementary information as a benchmark set for future potentials.

Pressure-induced transformations in amorphous silicon: A computational study

We study the transformations between amorphous phases of Si through molecular simulations using the environment dependent interatomic potential (EDIP) for Si. Our results show that upon pressure, the material undergoes a transformation from the low density amorphous (LDA) Si to the high density amorphous (HDA) Si. This transformation can be reversed by decompressing the material. This process, however, exhibits clear hysteresis, suggesting that the transformation LDA ↔ HDA is first-order like. The HDA phase is predominantly five-fold coordinated, whereas the LDA phase is the normal tetrahedrally bonded amorphous Si. The HDA phase at 400 K and 20 GPa was submitted to an isobaric annealing up to 800 K, resulting in a denser amorphous phase, which is structurally distinct from the HDA phase. Our results also show that the atomic volume and structure of this new amorphous phase are identical to those of the glass obtained by an isobaric quenching of the liquid in equilibrium at 2000 K and 20 GPa down to 400 K. The similarities between our results and those for amorphous ices suggest that this new phase is the very high density amorphous Si.

Inherent structure analysis of defect thermodynamics and melting in silicon

The inherent structure landscapes associated with various types of microstructure in crystalline silicon are probed within the empirical Environment-Dependent Interatomic Potential. Inherent structure sampling is carried out using a combination of molecular dynamics simulation and energy minimisation. Several isolated point defect clusters are analysed first, followed by the more complex phenomenon of crystal melting. In the latter case, both homogeneous and surface-driven melting phenomena are considered, and selected comparisons are made to face-centered cubic aluminium represented by an embedded-atom potential. It is shown that inherent structure theory provides a useful framework for consolidating various aspects of microstructure thermodynamics in crystals, and also naturally allows for the treatment of more complex situations such as applied mechanical fields.

Pressure effects on the transitions between disordered phases in supercooled liquid silicon

We investigate the pressure effects on the transitions between the disordered phases in supercooled liquid silicon through Monte Carlo simulations and efficient methods to compute free energies. Our calculations, using an environment dependent interatomic potential for Si, indicate that at zero pressure the liquid-liquid phase transition, between the high density liquid and the low density liquid, occurs at a temperature 325K below melting. We found that the liquid-liquid transition temperature decreases with increasing pressure, following the liquid-solid coexistence curve. As pressure increases, the liquid-liquid coexistence curve approaches the region where the glass transition between the low density liquid and the low density amorphous takes place. Above 5 GPa, our calculations show that the liquid-liquid transition is suppressed by the glassy dynamics of the system. We also found that above 5 GPa, the glass transition temperature is lower than that at lower pressures, suggesting that under these conditions the glass transition occurs between the high density liquid and the high density amorphous.

Crystallization of amorphous silicon induced by mechanical shear deformations

We have investigated the response of amorphous silicon (a-Si), in particular crystallization, to external mechanical shear deformations using classical molecular dynamics (MD) simulations and the empirical Environment Dependent Inter-atomic Potential (EDIP) [Phys. Rev. B 56, 8542 (1997)]. In agreement with previous results we find that, at low shear velocity and low temperature, shear deformations increase disorder and defect density. At high temperatures, however, the deformations are found to induce crystallization, demonstrating a dynamical transition associated with both shear rate and temperature. The properties of a-Si under shear deformations and the extent at which the system crystallizes are analyzed in terms of the potential energy difference (PED) between the sheared and non-sheared material, as well as the fraction of defects and the number of particles that possess a crystalline environment.

Amorphous silicon under mechanical shear deformations: Shear velocity and temperature effects

Mechanical shear deformations lead, in some cases, to effects similar to those resulting from ion irradiation. Here we characterize the effects of shear velocity and temperature on amorphous silicon ($a$-Si) modeled using classical molecular-dynamics simulations based on the empirical environment-dependent interatomic potential (EDIP). With increasing shear velocity at low temperature, we find a systematic increase in the internal strain leading to the rapid appearance of structural defects (fivefold-coordinated atoms). The impacts of externally applied strain can be almost fully compensated by increasing the temperature, allowing the system to respond more rapidly to the deformation. In particular, we find opposite power-law relations between the temperature and the shear velocity and the deformation energy. The spatial distribution of defects is also found to depend strongly on temperature and strain velocity. For low temperature or high shear velocity, defects are concentrated in a few atomic layers near the center of the cell, while with increasing temperature or decreasing shear velocity, they spread slowly throughout the full simulation cell. This complex behavior can be related to the structure of the energy landscape and the existence of a continuous energy-barrier distribution.

HRMC_1.1: Hybrid Reverse Monte Carlo method with silicon and carbon potentials

The Hybrid Reverse Monte Carlo (HRMC) code models the atomic structure of materials via the use of a combination of constraints including experimental diffraction data and an empirical energy potential. This energy constraint is in the form of either the Environment Dependent Interatomic Potential (EDIP) for carbon and silicon and the original and modified Stillinger–Weber potentials applicable to silicon. In this version, an update is made to correct an error in the EDIP carbon energy calculation routine. New version program summary Program title: HRMC version 1.1 Catalogue identifier: AEAO_v1_1 Program summary URL: http://cpc.cs.qub.ac.uk/summaries/AEAO_v1_1.html Program obtainable from: CPC Program Library, Queen's University, Belfast, N. Ireland Licensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.html No. of lines in distributed program, including test data, etc.: 36 991 No. of bytes in distributed program, including test data, etc.: 907 800 Distribution format: tar.gz Programming language: FORTRAN 77 Computer: Any computer capable of running executables produced by the g77 Fortran compiler. Operating system: Unix, Windows RAM: Depends on the type of empirical potential use, number of atoms and which constraints are employed. Classification: 7.7 Catalogue identifier of previous version: AEAO_v1_0 Journal reference of previous version: Comput. Phys. Comm. 178 (2008) 777 Does the new version supersede the previous version?: Yes Nature of problem: Atomic modelling using empirical potentials and experimental data. Solution method: Monte Carlo Reasons for new version: An error in a term associated with the calculation of energies using the EDIP carbon potential which results in incorrect energies. Summary of revisions: Fix to correct brackets in the two body part of the EDIP carbon potential routine. Additional comments: The code is not standard FORTRAN 77 but includes some additional features and therefore generates errors when compiled using the Nag95 compiler. It does compile successfully with the GNU g77 compiler ( http://www.gnu.org/software/fortran/fortran.html). Running time: Depends on the type of empirical potential use, number of atoms and which constraints are employed. The test included in the distribution took 37 minutes on a DEC Alpha PC.

Graphene buffer layer on Si-terminated SiC studied with an empirical interatomic potential

The atomistic structure of the graphene buffer layer on Si-terminated SiC is investigated using a modified version of the environment-dependent interatomic potential. The determination of the equilibrium state by the conjuguate gradients method suffers from a complex multiple-minima energy surface. The initial configuration is therefore modified to set the system in specific valleys of the energy surface. The solution of minimal energy forms a hexagonal pattern composed of stuck regions separated by unbonded rods that release the misfit with the SiC surface. The structure presents the experimental symmetries and a global agreement with an ab initio calculation. It is therefore expected that the interatomic potential could be used in classical molecular dynamics calculations to study the graphene growth.