Density functional theory study of the mechanical behavior of silicene and development of a Tersoff interatomic potential model tailored for elastic behavior

The supercell of silicene is indicated by the dotted line in figure 1(a). The supercell contains four Si atoms and its x and y axes correspond to the zigzag (ZZ) and armchair (AC) directions, respectively. Two metastable structures were reported in Cahangirov et al [25]: buckled (Hex.2b) and planar (Hex.2p) structures the side views of which are presented in figures 1(b) and (c), respectively. According to DFT calculations, Hex.2b is more stable than Hex.2p in terms of cohesive energy. The DFT row in table 1 lists the lattice constant a, buckling parameter Δ, and cohesive energy of each relaxed structure.

Zoom In Zoom Out Reset image size

DFT calculations were carried out with the Vienna Ab initio Simulation Package (VASP) [32–35] in which a projector augmented-wave method was used as a plane wave basis set. We used the generalized-gradient approximation (GGA) for the exchange-correlation functionals to compare the results between the two. A 25 × 25 × 1 k-point mesh for the Brillouin zone was sampled using the Monkhorst–Pack scheme. The length of the supercell in the z direction was chosen to be 20 Å (or L_z  = 20 Å), which was long enough to minimize the image interaction. The energy cutoff was extended to 500 eV for precise calculations and the convergence criterion was selected to be 10⁻⁵ eV for electronic relaxation and 10⁻⁴ eV for ionic relaxation between relaxation steps. Every structure was relaxed with an energy minimization method, e.g. the Blocked–Davidson scheme [36].

We conducted three types of mechanical test simulation with a buckled silicene structure to calculate the elastic properties: Young's modulus (Y^2D), Poisson's ratio (ν), and bulk modulus (B^2D). The details of the calculation method are described in the electronic supplementary materials S1 (available online at stacks.iop.org/NANO/32/295702/mmedia). The bottom row of table 2 lists the 2D Young's modulus (Y_2D), Poisson's ratio ($\nu$), and 2D bulk modulus (B_2D) of buckled silicene from a DFT calculation of mechanical tests. Y_2D in the ZZ and AC directions are 61.7 (N m⁻¹) and 61.8 (N m⁻¹), respectively, and $\nu$ is obtained as 0.31 regardless of the direction. The elastic constants are similar to previous DFT calculations [5, 7, 15, 16, 37]. Since the 2D Young's moduli along the ZZ and AC directions are almost identical, the isotropic behavior of silicene is expected under small deformations. For isotropic 2D material, B_2D is related to Y_2D and $\nu$ as in equation (1) [38].

$$\begin{eqnarray}&&{B}_{2D}=\displaystyle \frac{{Y}_{2D}}{2(1-\nu )}\end{eqnarray} \tag{ 1 }$$

Table 2. 2D Young's modulus Y_2D, Poisson 's ratio ${\boldsymbol{\nu }},$ and 2D bulk modulus B_2D of buckled silicene. In the references with the asterisks (*), isotropy is assumed.

		Y2D (N m−1)		${\boldsymbol{\nu }}$
	Potentials	ZZ axis	AC axis	ZZ axis	AC axis	B2D (N m−1)	$\tfrac{{{\boldsymbol{Y}}}_{2{\boldsymbol{D}}}/{{\boldsymbol{B}}}_{2{\boldsymbol{D}}}}{2(1-{\boldsymbol{\nu }})}$
MS	MA2013 [29]	52.6	52.6	0.03	0.03	27.4	0.993
	MEAM [30]	26.7	26.8	0.22	0.22	17.3	0.988
	SW92 [31]	37.7	37.7	0.17	0.17	22.7	1.004
	SW2014 [24]	39.6	39.6	0.11	0.11	22.3	0.996
	ReaxFF [17]	41.6	41.6	0.26	0.26	28.2	0.996
	This study	59.8	59.8	0.17	0.17	36.0	0.997
DFT	Qin et al [5]*	63.0	63.0	0.31	0.31	·	·
	Zhao [15]	60.06	63.51	0.41	0.37	·	·
	Peng et al [16]*	63.8	63.8	0.33	0.33	·	·
	Goli et al [37]	66.72	67.63	·	·	40.11	·
	This study—GGA	61.7	61.8	0.31	0.31	45.0	0.994

Thus, we checked the ratio of ${Y}_{2{\rm{D}}}/2\left(1-\nu \right)$ and B_2D and the value in the bottom-right cell of table 2 is 0.994, which is close to 1 and thus indicates that silicene is isotropic for small strain.

Uniaxial tensile tests were carried out for silicene of buckled and planar structures up to ε = 30% in both the ZZ and AC directions. Figure 2 presents the per-atom energy curves and the buckling parameter (Δ). When the buckled silicene is under tension in the ZZ direction, a structural transition to the planar one is observed at ε₁₁ = 21% (green arrow in figure 2(a)), after which the per-atom energy curve (solid black line) of the buckled structure is overlapped with that of the planar one (solid red line). The transition strain of ε₁₁ = 21% is clearly captured in figure 2(b) as the value of the buckling parameter is less than 10% of the initial value (Δ₀ = 0.45 Å). This result is different from that reported by Peng et al [16], according to which the buckling parameter is maintained as about 75.6% of Δ₀ even at ε₁₁ = 18.3%. Similarly, Mortazavi et al also reported that the buckling parameter decreases only by 17.0% at ε₁₁ = 9.5% [39]. The secret to the difference between our results and the previous studies lies in the loading condition. In Peng et al [16] and Mortazavi et al [39] the silicene was under a uniaxial stretch condition, while it is under a uniaxial tensile condition in this study. We were able to reproduce previous results when we repeated the simulations under a uniaxial stretch condition. At ε₁₁ = 25%, the solid lines in figure 2(a) abruptly drop and increase again, implying that another structural change may occur. Indeed, the hexagonal lattice structure transforms into a rectangular lattice structure at ε₁₁ = 25% in figure 3(a). On the other hand, when silicene is under tensile loading along the AC direction, the buckling parameter only decreases by about 16% until ε₂₂ = 21% (see figure 2(b)) and the structure remains planar at a large strain of ε₂₂ = 40% in figure 3(b). Evidently, in figure 2(a), the per-atom energy of the planar structure (dashed red line) is always larger than that of the buckled structure (dashed black line) in the range of ε₂₂ ≤ 30%. In figure 3(b), the atomic bond, parallel to the AC tensile direction, is stretched from l_bond = 2.778 Å to 3.348 Å at ε₂₂ = 40%, while other bonds experience a slight increase of 0.005 Å. This separation tears silicene into several 1D zigzag Si chains. Similar fracture behavior for hexagonal boron nitride was reported in the tensile tests along the AC direction [40].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

There are 13 adjustable parameters (A, B, λ₁, λ₂, λ₃, α, β, n, c, d, h, R, and D) in the Tersoff model [27] but we keep α, R, and D and only consider the remaining 10 parameters during optimizations (the equations of the Tersoff model are in the supplementary S2). Following the philosophy of the original Tersoff model, the bond-order function b_ij in the attractive term is determined relative to a_ij in the repulsive term and hence α is fixed at 0 so that a_ij is always 1. We used the same cut-off function as in the original work and R and D were fixed to 3.0 and 0.2, respectively.

For potential parameter optimization, we consider four different metastable 2D Si structures: buckled (Hex.2b), planar (Hex.2p), hexagonal (Hex.1p), and rectangular (Tetra.1p). Figure 4 shows the per-atom energies of the four structures from the DFT-GGA calculation of the biaxial stretch test in terms of per-atom occupation area or (area)/(number of atoms in the unit cell), and the buckled structure is again the most stable of the four. Parameter optimization is performed by matching these energy curves so that the model can predict to the 2D bulk modulus (B_2D) of buckled silicene (Hex.2b) and other elastic constants (Y_2D and ν) to a reasonable extent, even though we do not explicitly include these elastic constants during the optimization process.

Zoom In Zoom Out Reset image size

For this purpose, the authors proposed the objective function f_obj in equation (2), where X is an array of ten fitting parameters and thus a design variable in the optimization process.

$$\begin{eqnarray}&&{f}_{{\rm{obj}}}\left({\bf{X}}\right)\equiv {f}_{A}\left({\bf{X}}\right)+{f}_{B}({\bf{X}})\end{eqnarray} \tag{ 2 }$$

The first term f_A is a penalty term that was introduced to match the overall shape of the energy curves for all four structural types in MS as closely as possible to those obtained from DFT calculations.

$$\begin{eqnarray}&&{f}_{A}\left({\bf{X}}\right)\equiv \displaystyle \displaystyle \sum _{\alpha }\displaystyle \sum _{m=1}^{{N}_{\alpha }}{\left({\rm{\Delta }}{E}_{\alpha ,m}^{MS}\left({\bf{X}}\right)-{\rm{\Delta }}{E}_{\alpha ,m}^{DFT}\right)}^{2}\end{eqnarray} \tag{ 3 }$$

where α ∈ {Hex.2b, Hex.2p, Tetra.1p, Hex.1p} and N_α is the number of sampled points (denoted as dots) from each curve in figure 4. We have sampled more points (N_α=Hex.2b = 11) from the Hex.2b curve (black solid line) than the other three because the Hex.2b structure is the most stable. In total, 26 configurations are sampled within a ±10% change of equilibrium lattice constant (a₀) from four different types of structure for optimization. The energy difference ΔE in equation (4) is defined in equations (4) and (5).

$$\begin{eqnarray}&&{\rm{\Delta }}{E}_{\alpha ,m}^{MS}\left({\bf{X}}\right)={E}^{MS}\left({\bf{X}};{\left\{{\bf{x}}\right\}}_{\alpha ,m}\right)-{E}^{MS}\left({\bf{X}};\{{{\bf{x}}}_{0}\}\right)\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{\rm{\Delta }}{E}_{\alpha ,m}^{DFT}={E}^{DFT}({\left\{{\bf{x}}\right\}}_{\alpha ,m})-{E}^{DFT}\left(\{{{\bf{x}}}_{0}\}\right)\end{eqnarray} \tag{ 5 }$$

where {x}_α,m is the atomic configuration of the mth sampled point in the αth structural type and {x₀} is the atomic configuration at the energy minimum point of the Hex.2b structure obtained from DFT. The second term f_B is the correction term to shift the energy point evaluated at {x₀} from MS to be as close as possible to the one in DFT (see equation (6)).

$$\begin{eqnarray}&&{f}_{B}\left({\bf{X}}\right)\equiv {\left({E}^{MS}\left({\bf{X}};\{{{\bf{x}}}_{0}\}\right)-{E}^{DFT}(\{{{\bf{x}}}_{0}\})\right)}^{2}\end{eqnarray} \tag{ 6 }$$

The optimal parameter set is searched by minimizing the objective function f_obj with a constraint that the gradient of the interatomic potential energy would be zero at {x₀} or ${(\partial {E}^{MS}/\partial {\bf{x}})}_{{{\bf{x}}}_{0}}=0.$ Since f_obj has multiple local minima and is not a simple quadratic function of X, the search results are dependent on the parameters' initial value. Therefore, multiple initial sets X₀ are produced by setting up a suitable grid in the design space. The preset range of each parameter is equally divided and combinations of initial parameters are made from each grid point. Minimizing f_obj is done by the MATLAB command fmincon, which uses the interior-point algorithm for constrained nonlinear optimization [41, 42]. The number of maximum iterations is fixed at 10⁶ and the criterion of the stopping residual is set to 10⁻⁶.

Table 3 lists the optimized Tersoff parameters and their fitting range. With the optimized parameters, we evaluate buckled (Hex.2b) silicene to obtain the lattice constant of a = 3.7 Å, the buckling parameter Δ = 0.45 Å and the cohesive energy of 4.855 eV in table 1. It is unsurprising that our model outperforms the other potential models: EDIP [26], Tersoff [27], Murty and Atwater (MA95 [28], MA2013 [29]), Stillinger and Weber (SW92 [31], SW2014 [24]), MEAM [30], and ReaxFF [17] because the buckled structure at the energy-minimum point is included during optimization. Among the tested potentials, EDIP, Tersoff, and MA95 produce a planar (Hex.2p) structure as the most stable one in contrast to the DFT prediction and we exclude these three models from the subsequent benchmark test in section 4. Surprisingly, our parameter set predicts that not only Hex.2b but also Hex.2p, Tetra.1p and Hex.1p exist as metastable structures. This polytype feature is attributed to the bond-order nature of the Tersoff potential. Table 1 includes the lattice constants (a) and cohesive energies (E_c ) of Hex.2b, Tetra.1p, and Hex.1p as well.

The reference energy curves in figure 4 were reproduced with the optimized Tersoff model and plotted in figure 5. Figure 5(a) shows the energy curves of the buckled structure (Hex.2b) and planar structure (Hex.2p). The dashed curves are from the Tersoff model and they follow the trend of the DFT-GGA energy curves (solid lines) relatively well, although the entire curves are shifted downward by about 0.2 eV. A good factor of this parameter set is that the crossing between Hex.2b and Hex.2p energy curves is captured, which is crucial for predicting the structural change of silicene under extensional loading. Figure 5(b) compares the results of DFT-GGA and the optimized model for rectangular structure (Tetra.1p) and hexagonal structure (Hex.1p). Note that, contrary to the DFT-GGA results (solid lines), the hexagonal structure (Hex.1p) has lower energy than the rectangular structure (Tetra.1p) in the Tersoff model (dashed lines). While this opposite behavior initially looks disappointing, these two structures will not be observed in MD simulations in ambient conditions because both the Tetra.1p and Hex.1p structures have much higher energy than the buckled (Hex.2b) and planar (Hex.2p) structures. Therefore, we take the parameter set in table 3 as our best.

Zoom In Zoom Out Reset image size

As a benchmark test, we compared the elastic properties predicted from five different interatomic potential models (MA2013 [29], SW92 [31], SW2014 [24], MEAM [30], and ReaxFF [17]) together with ours. The crystal orientation of silicene used in MS was set to be the same as that used in DFT calculations (ZZ direction along the x axis and AC direction along the y axis). The initial structures were created using the conjugated-gradient (CG) method at 0 K. The tensile test along the ZZ direction used a silicene strip of 5120 atoms the area of which was about 32 nm × 11 nm. The test along the AC direction used a silicene strip of 4928 atoms the area of which was 30 nm × 11 nm. Periodic boundary conditions were applied in all directions and the simulation cell size in the silicene-plane normal direction was set to 8 nm, which was long enough to prevent image interaction. During the test, CG relaxation with the Parrinello–Rahman scheme was applied to release stress components except the tensile component and to find a minimum-energy structure at each strain increment. All MS simulations were performed using the Large-scale Atomic/Molecular Massively Parallel Simulator [43] and visualized with OVITO [44] and Atomeye [45]. Among the tested potentials, only ReaxFF and SW2014 were developed specifically for silicene. SW2014 was intended to better predict the phonon curve of silicene than the original SW (SW85 [46]). ReaxFF was developed to find the elastic constants and study the fracture pattern and edge reconstruction of silicene during uniaxial tensile tests.

In table 2, elastic constants (Y_2D, ν, and B_2D) from our model are compared with those from other models. All tested potentials predict the isotropic behavior of silicene. Y_2D from our study is 59.8 (N m⁻¹) in both the ZZ and AC directions and is closest to the reference DFT value (61.4 N m⁻¹). The 2D bulk modulus B_2D is 36 (N m⁻¹), which is 80% of the reference value and still better than any other model. The Poisson ratio ν is 0.17, which is only 55% of the reference but reasonably fair, considering the others' performance. To better visualize the performance of our model at a glance, the authors present figure 6, which is a scatter plot of the elastic constants from the tested parameter sets in the (Y_2D, B_2D) plane. The point (in diamond red marker) from the optimized parameter set is closest to that of DFT-GGA marked with ×. This scatter plot shows that our model outperforms other potential models in that it better predicts silicene's elastic behavior overall.

Zoom In Zoom Out Reset image size

The uniaxial tensile test was carried out quasi-statically with six potential models until ε ≤ 30%. Figures 7(a) and (b) show the stress and strain curves in the ZZ and AC directions, respectively. Our optimized model (red line) closely follows the DFT curve (black line) up to ε₁₁ = 8.1% in the ZZ direction and about ε₂₂ = 10% in the AC direction. All other potential models show much softer behavior across the entire strain range. In figure 7(a), our model starts to deviate at ε₁₁ = 8.1% where the silicene structure changes from Hex.2b to Hex.2p, while DFT predicts that such a transformation occurs at ε₁₁ = 21%. This early transformation is explained by the early crossing of the energy curves of Hex.2b and Hex.2p in our model in figure 5(a). Our model also shows a cusp in the stress curve at the transformation strain and abrupt stiffened behavior after the cusp, which was not observed from the DFT curve within our stress accuracy. As loading beyond the transformation strain (ε_T ) and unloading produce hysteresis in the stress curve, we can say that ε_T is the elastic limit. Interestingly, this stiffening after a structural change is universally observed from all tested interatomic potentials. The Hex.2b structure of SW92, MA2013, MEAM and ReaxFF transformed into Hex.2p at ε₁₁ = 24.3%, 25.6%, 25.8% and 23.1%, respectively. Only the SW2014 model could predict a similar transformation strain of ε₁₁ = 8.45% to our study. We speculate that the transformational stiffening is related to the higher elastic modulus of the planar silicene than that of the buckled, regardless of the tested potentials.

Zoom In Zoom Out Reset image size

Figure 7(b) shows tensile tests along the AC direction that show how the Hex.2b structure is also altered into a Hex.2p structure for all potential models, while the DFT calculation predicts such transition at a much larger strain of ε₂₂ = 40%. The dependency of structural transition strain (ε_T ) on the tensile loading direction implies the development of anisotropic behavior in silicene. Indeed, by comparing the stress curves between figures 7(a) and (b), it can be noted that the DFT results only show isotropic behavior in small strain regimes (<4%) but strong anisotropy for ε > 4% while no potential models can capture the appearance of anisotropy at this early stage of strain. The anisotropy is confirmed if σ₁₁ and σ₂₂ differ by more than 0.1 N m⁻¹, which is approximately of the order of 0.1 GPa in 3D, a stress tolerance used in a DFT study of anisotropy in graphene [47]. Based on the same criteria, the anisotropy appears at ε = 6.4%, 12.4%, 7.1%, 6.8%, 11.5% and 7.1% for MA2013, MEAM, SW92, SW2014, ReaxFF, and our model, respectively. The authors expect that the puckered shape of the buckled structure causes the early development of anisotropy in silicene. This explanation is plausible if we compare silicene to graphene. Graphene, which is analogous to silicene but completely flat, maintains isotropic behavior until 15% strain according to DFT studies [47–49]. Although our optimized Tersoff model does not completely reproduce every feature reported in DFT calculations, it is the most reliable model so far, at least in the elastic regime.

It is noticeable that a discontinuous slope increase appears in the stress–strain curve when the structure changes into Hex.2p in figure 7. This abrupt stiffening is associated with the geometric constraint that is inherent in the Hex.2p structure, where the sum of all three triplet angles between an atom and its three neighbor atoms must be 360°. That means that one angle (θ₁) is dependent on the other (θ₂)—as shown in inlet figure 8(a)—once the silicene becomes planar. To investigate the effect of angle function g(θ) on the potential energy, the per-atom energy depending on the triplet angle is calculated via DFT and the MS-Tersoff model with the optimized parameter set (see figure 8(a)). In this calculation, the interatomic distance is fixed at 2.278 Å, the bond length at the most stable configuration of Hex.2b. The difference between the results of DFT and MS becomes noticeable for θ > 130°. The energy increase is steeper in MS-Tersoff than DFT as θ increased, which can be attributed to the steep increase in the original angle penalty function g(θ) in the Tersoff model with the triplet angle (θ) in figure 8(b). We need a mildly increasing penalty function to better describe the angle effect.

Zoom In Zoom Out Reset image size

A fourth-order polynomial function is selected as the new angle penalty function g*(θ) (equation (7)) that is padded to the original function after 120°. In other words, we use the original function g(θ) for θ ≤ 120° and the new function g*(θ) for θ > 120°. Thus, the elastic modulus of Hex.2b (θ = 116.2°) is unaffected by a new angle penalty function. Additional optimization of g*(θ) is conducted using the method described in section 3.1, with six parameters (p₀, p₁, p₂, p₃, p₄, and h₂). The reference configurations used in this second round of optimization are denoted as dots in figure 8(a), or α ∈ {Hex.2b, Hex.2p} in equations (4)–(6). Three nonlinear constraints in equations (8)–(10) are imposed during optimization to preserve C2 continuity at 120°, which is crucial for stress continuity.

$$\begin{eqnarray}&&\begin{array}{l}g* \left(\theta \right)={p}_{0}+{p}_{1}\left({h}_{2}-\,\cos \,\theta \right)+{p}_{2}{\left({h}_{2}-\cos \theta \right)}^{2}\\ \,+\,{p}_{3}{\left({h}_{2}-\cos \theta \right)}^{3}+{p}_{4}{\left({h}_{2}-\cos \theta \right)}^{4}\end{array}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{\theta \to 120^\circ -}g\left(\theta \right)=\mathop{\mathrm{lim}}\limits_{\theta \to 120^\circ +}g* \left(\theta \right)\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{\theta \to 120^\circ -}\displaystyle \frac{dg\left(\theta \right)}{d\theta }=\mathop{\mathrm{lim}}\limits_{\theta \to 120^\circ +}\displaystyle \frac{dg* \left(\theta \right)}{d\theta }\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{\theta \to 120^\circ -}\displaystyle \frac{{d}^{2}g\left(\theta \right)}{d{\theta }^{2}}=\mathop{\mathrm{lim}}\limits_{\theta \to 120^\circ +}\displaystyle \frac{{d}^{2}g* \left(\theta \right)}{d{\theta }^{2}}\end{eqnarray} \tag{ 10 }$$

The lower part of table 3 summarizes the optimized parameters for g*(θ), which is plotted in figure 8(b) together with g(θ). With g*(θ), the per-atom energy curve clearly follows the trend of the DFT curve after 130° in figure 8(a). The constant shift between MS and DFT remains because we use the original function g(θ) for θ ≤ 120° and the per-atom energy curve must be continuous at θ = 120° even though it is not smooth. Uniaxial tensile tests were performed to calculate the elastic moduli and draw stress–strain curves to confirm the effect of g*(θ). The elastic moduli, including Y_2D and ν, were maintained regardless of the direction, as intended. Figure 9 shows that the slope after transformation to the Hex.2p structure decreased in uniaxial tensile tests along both the ZZ and AC directions.

Zoom In Zoom Out Reset image size