Home
Search
Collections
Journals
About
Contact us
My IOPscience
Molecular dynamics simulations of mechanical properties of monolayer MoS2
This content has been downloaded from IOPscience. Please scroll down to see the full text. 2015 Nanotechnology 26 185705 (http://iopscience.iop.org/0957-4484/26/18/185705) View the table of contents for this issue, or go to the journal homepage for more
Download details: IP Address: 150.216.68.200 This content was downloaded on 25/04/2015 at 17:08
Please note that terms and conditions apply.
Nanotechnology Nanotechnology 26 (2015) 185705 (10pp)
doi:10.1088/0957-4484/26/18/185705
Molecular dynamics simulations of mechanical properties of monolayer MoS2 Si Xiong and Guoxin Cao HEDPS, Center for Applied Physics and Technology, Department of Mechanics and Engineering Science, College of Engineering, Peking University, Beijing, 100871, People’s Republic of China E-mail: [email protected] Received 14 February 2015, revised 18 March 2015 Accepted for publication 20 March 2015 Published 16 April 2015 Abstract
Based on the benchmark created by the first principle calculations, the effectiveness of different empirical potential functions, including CVFF1, CVFF2, SW and REBO potentials, on describing the mechanical behavior of single layer MoS2 (SLMoS2) are evaluated. The mechanical properties including the elastic modulus E, Poisson’s ratio ν, nonlinear elastic modulus D, failure stress and ultimate strain are considered. It is found that under a small deformation, the REBO and CVFF2 potentials can provide an effective description of the elastic behavior of SLMoS2; whereas under a large deformation, the SW potential gives a more accurate stress value and the CVFF2 potential can only predict the right value under biaxial tension. After modifying the cut-off distances in the REBO potential, the failure strain can be reasonably predicted by the REBO potential; whereas the failure stress will be overestimated by about 20–40% for biaxial and uniaxial tension. The present study can provide help on accurately understanding the mechanical behavior of SLMoS2. Keywords: mechanical properties, single layer MoS2, molecular dynamics (Some figures may appear in colour only in the online journal) 1. Introduction
between the AFM load and the deformation of the freestanding structure. Besides experimental work, the reports of the mechanical properties of layered MoS2 are also developed theoretically [13–20]. Under applied in-plane tension, the nonlinear elastic properties along with the intrinsic strength/ultimate strain of SLMoS2 are reported from first-principles calculations [15, 18, 20]. Although first-principles calculations are the most accurate theoretical approach, they have the following limitations: (a) they can only simulate several hundred atoms; (b) more complicated mechanical behavior, including indentation response and bending behavior, cannot be investigated. Therefore, molecular dynamics (MD) simulations based on the empirical potentials (force field) are desirable. The accuracy of MD simulations mainly depends on the effectiveness of empirical potential function selected. Wakabayashi et al [21] developed a constant valence force field (CVFF) type potential, in which the bond interactions between Mo and S atoms (including both bond stretching and bond angle variation terms) were described by harmonic potentials. Morita et al [22] and Becker et al [23] provided the revised empirical
SLMoS2 has a direct band gap of 1.8 eV [1], which makes it a good candidate for applications in nanoelectronics, such as field-effect transistors (FET) [2–5], phototransistors [6, 7], nonvolatile memory cells [8] and so on. In addition, MoS2 can be used as a solid lubricant [9]. To effectively design structures for its practical applications, the mechanical properties of single layer MoS2 (SLMoS2) should be well understood. Since SLMoS2 is typically in the atomic-level thick range, its mechanical properties are highly difficult to measure by the conventional experimental methods. Currently, the mechanical properties of 2D materials, such as graphene and layered MoS2, are mainly measured from its free-standing structure under clamped boundary conditions [10–12]. 2D materials are firstly mounted on the substrate with cylindrical wells (holes) (called a free-standing structure), and then deformed by the tip of an atomic force microscope (AFM). The mechanical properties, including elastic modulus and intrinsic strength, can be estimated from the relationship 0957-4484/15/185705+10$33.00
1
© 2015 IOP Publishing Ltd Printed in the UK
Nanotechnology 26 (2015) 185705
S Xiong and G Cao
potential function formula. Both of them used the Morse potential to describe the Mo–S bond stretching interaction and the non-bond interactions between Mo and S atoms were considered by both the Lennard–Jones (LJ) potential and the electrostatic interaction. By comparing with the experimental values of MoS2 lattice parameters, Varshney et al [24] corrected the parameters in the aforementioned empirical potential. Liang et al [25] parameterized an interatomic potential for the Mo–S system based on the second-generation reactive empirical bond-order (REBO) formalism. Jiang et al [13] presented a parameterization of the Stillinger-Weber (SW) potential to describe the interatomic interactions within single-layer MoS2. Based on the developed empirical potentials, more mechanical behaviors of MoS2 were investigated using MD simulations: such as indentation response [16], buckling behavior [19], bending behavior [14] and temperature effect[17]. However, the effectiveness of the aforementioned empirical potentials on investigating the mechanical behavior of layered MoS2 has not been fully evaluated. In the reported results, the experimentally obtained Mo–S lattice parameters [24] and phonon spectrum [13] of bulk MoS2 were selected as the benchmark of the effectiveness of empirical potential. However, it is still not clear whether or not these empirical potentials can effectively describe the deformation of layered MoS2 or which one among these empirical potentials is most effective for studying the mechanical behavior of layered MoS2, especially for the nonlinear elastic behavior or the failure behavior. In the present work, first-principles calculations are firstly used to investigate the tensile properties of SLMoS2, including elastic modulus, Poisson’s ratio, nonlinear elastic modulus, failure stress and ultimate strain. Then, the similar tensile behavior is studied by MD simulations based on different empirical potentials, including CVFF, SW and REBO potentials. Finally, the effectiveness of different empirical potentials is analyzed for describing the mechanical behavior of SLMoS2. The present study can provide a useful guideline to understand the mechanical properties of SLMoS2 and to design the structure of layered MoS2 in its practical applications.
of SLMoS2 determined from our DFT calculations, the spacing between two Mo atoms λMo–Mo = 0.3181 nm, the distance between two S atom planes δ = 0.3114 nm and the Mo–S bond length λMo–S = 0.2408 nm, which are in good agreement with the reported simulation results [15, 18] and the experimental results of bulk MoS2 crystal (λMo–Mo = 0.3162 nm and δ = 0.318 nm in bulk) [12]. The geometric parameters of the equilibrium structure SLMoS2 determined by DFT calculations are displayed in table 1. The uniaxial/biaxial strain tension is applied by increasing the computational cell size along the loading direction.
2.2. MD simulations of SLMoS2
SLMoS2 has a tri-layer structure (figures 1(c) and (d)). Each Mo atom locates in the center of the triangular prism formed by six neighboring S atoms and the S–Mo–S bond is mainly covalent in nature. In MD simulations, different empirical potentials are selected to optimize the equilibrium structure of SLMoS2 as well as its deformed structure under in-plane tension: including the CVFF potential [24], SW potential [13] and REBO potential [16]. The force field parameters for all different empirical potentials used in the present work are the same as those reported in the original literatures, unless specifically stated that they were modified. MD simulations are carried out using the LAMMPS package [29]. A periodic boundary condition is applied in the lateral directions to remove the lateral boundary effect and simulate the intrinsic properties of SLMoS2. For all simulations, the computational cell includes more than 900 atoms and the size is about 4.74 × 5.47 × 9 nm, and the lateral size is slightly varying with the empirical potentials selected. The simulation results are not sensitive to the selected computational cell size, which has been checked using a larger computational cell size. The uniaxial/biaxial strain tension is applied on the computational cell using the NVT ensemble. The strain rate is kept as 109/s for all cases. The system is kept as 1 K and a constant time step of 1 fs is selected.
2.2.1. CVFF potentials. In the CVFF potential, the potential
2. Computational methods
function can be expressed as the summation of the bondstretching energy (Ubond), the bond angle variation energy (Uangle), the LJ non-bond energy (ULJ) and the Coulomb electrostatic energy (UC) [24]:
2.1. First-principles calculations of 2D SLMoS2
In the present work, the mechanical response of SLMoS2 under the uniaxial/biaxial tension is firstly investigated by density functional theory (DFT) calculations, which are performed using the Vienna Ab Initio Simulation Package (VASP) [26, 27]. DFT calculations are performed based on the six-atom (including two Mo atoms and four S atoms) SLMoS2 unit cell, and the unit cell size is 3.1809 × 5.5096 × 61.5 Å (shown in figures 1(a) and (b)). The generalized gradient approximation is applied using the Perdew–Burke–Ernzerhof functional [28]. A Monkhorst–Pack k-point mesh grid 18 × 18 × 3 is selected, and the plane wave energy cut-off is set at 500 eV. For the equilibrium structure
U=
∑Ubond + ∑Uangle + ∑ULJ + ∑Uc
Ubond = K1 × ( rij − r0 )2 (Harmonic bond potential) 2 Ubond = D × ⎡⎣ 1 − e−α ( rij− r0 ) ⎤⎦ (Morse bond potential)
(
)2
Uangle = K2 × θ ijk − θ0 (Harmonic angle potential) ⎡ ⎛ ⎞12 ⎛ σ ⎞6 ⎤ σ ULJ = ε⎢ ⎜⎜ ⎟⎟ − 2 ⎜⎜ ⎟⎟ ⎥ (LJ nonbond potential) ⎢ ⎝ rij ⎠ ⎝ rij ⎠ ⎥⎦ ⎣ 2
Nanotechnology 26 (2015) 185705
S Xiong and G Cao
Figure 1. Atomic structure of SLMoS2, (a) (b) the computational cell used in the DFT calculations; (c) (d) the top and side view of the computational model used in MD simulations. Table 1. Geometric parameters of SLMoS2.
Parameters λMo–Mo (nm) λMo–S (nm) δ (nm) θMo–S–Mo (°) θS–Mo–S (°)
Uc = C
qi q j rij
DFT
CVFF1
CVFF2
SW
REBO
0.3181 0.2408 0.3114 82.688 80.582
0.3267 0.2365 0.2852 87.384 74.189
0.3169 0.2313 0.2831 86.471 75.447
0.3064 0.2398 0.3238 79.417 84.927
0.3146 0.2438 0.3253 80.357 83.687
(Coulomb electrostatic potential)
2.2.2.
SW potential. In the SW potential, the bond interaction is described by both a two-body interaction (U2(i, j)) and a three-body interaction (U3(i, j, k)) [13]:
(1)
where K1, K2, D are the bond potential parameters, σ and ε are the LJ non-bond potential parameters and C is the Coulomb electrostatic potential parameter. rij is the distance between atoms i and j; θijk is the angle between bond ij and bond jk; qi and qj are the partial charge values for atoms i and j; r0 is the equilibrium bond length between atoms i and j; and θ0 is the equilibrium angle between bond ij and bond jk. The bond stretching potential can be expressed by both the harmonic potential and the Morse potential [24], which are defined as CVFF1 and CVFF2, respectively, in the present work. The potential parameters used in equation (1) are taken from Varshney et al [24].
U=
∑U2 (i ,
j) + ⎡
∑U3 (i , −1
U2 (i , j ) = S1e⎣ ρ ( rij− rmax ) U3 (i , j , k ) =
(
⎤ ⎦
(S r
⎡ S3 e⎣ ρ1 ( rij− rij _
× cos θ jik − cos θ0
j, k )
)2
−4 2 ij
max
)−1 + ρ
−1
)
2 ( rjk − rjk _ max
)−1 ⎤⎦
(2)
where the parameters S1, S2, S3, ρ, ρ1, ρ2, rmax, rij_max and rjk_max are taken from Jiang et al [13].
3
Nanotechnology 26 (2015) 185705
S Xiong and G Cao
2.2.3. REBO potential. In the REBO potential, the bond
u=
interaction is described by a many-body function [16, 25]: 1 ∑ fijc ( rij ) ⎡⎣ V R ( rij) − bij V A ( rij ) ⎤⎦ 2 ⎛ Qij ⎞ ⎟⎟ V R ( rij ) = Aij e−αij rij ⎜⎜ 1 + rij ⎠ ⎝
Ub =
(6)
The stress–strain relationship can be calculated by the derivative of u with respect to ε:
V A ( rij ) = Bij e−βij rij bij = ⎡⎣ 1 +
∑
fikc
1 Ciii εi2 + … 2 (unaxial strain tension)
σi = Cii εi +
−1/2 ( rik ) ⋅ G ⎡⎣ cos θ ijk ⎤⎦ + P ( Ni ) ⎤⎦
( )
P ( Ni ) = −a 0 ( Ni − 1) − a1e−a 2 Ni + a 3 Mo
Ni =
(i = 1, 2)
σb = ( C11 + C22 + 2C12 )
S
εi + ( 2C111 − C222 + 3C112) εi2 + … (biaxial strain tension)
∑ fikc ( rik ) + ∑ filc ( ril )
i≠k fijc ( rij )
1 ( C11 + C22 + 2C12 ) εi2 2 1 + ( 2C111 − C222 + 3C112 ) εi3 + … (biaxial tension) 3
(7)
i≠l
⎧1 ⎪ ⎪ ⎧ ⎡ π r − R min ⎪1⎪ ij ij ⎨ = ⎨ 1 + cos ⎢ ⎢ max min ⎪2⎪ ⎣ Rij − Rij ⎪ ⎩ ⎪ ⎩0
(
(
) ⎤⎥ ⎫⎪⎬ ) ⎥⎦ ⎪⎭
rij ⩽
where i = 1, 2 represents the zigzag (zz) and armchair (ac) directions, respectively. Figure 2 shows the calculated values of σi and σb/2, which are very close to the reported results [15, 18]. Similar to graphene, the σ–ε relationships are almost the same along the zz and ac direction when ε < 10%, which means that SLMoS2 can be considered to be isotropic with ε < 10%. Theoretically, the isotropic, linear elastic properties: elastic modulus E and Poisson’s ratio ν can be calculated as:
Rijmin
Rijmin ⩽ rij ⩽ Rijmax Rijmax ⩽ rij (3)
(
)
2 2 E = C11 − C12 / C11 and ν = C12 / C11.
where VR and VA are the pairwise repulsive and attractive interactions, bij is the many-body bond-order function, Ni is the total number of nearest neighbors within the maximum REBO cutoff distance and fijc is the cutoff function that limits the REBO interactions to nearest neighbors. The pairwise parameters α, β, A, B, Q, Rijmin , Rijmax and a0–a1 depend on the chemical species of the interacting atoms i and j, which are taken from Stewart et al [16].
(8)
For graphene, the isotropic, nonlinear elastic modulus D is introduced to describe the nonlinear elastic behavior under a large tension [30, 31], and thus, the stress-strain relationship can be simplified as: σ = Eε + Dε2 1 D= C111 1+3ν 2 − C222 ν 3 + 3ν 2 2
{ (
)
(
+ 3C112 ν 2 − ν
(
)
)}
(9)
The intrinsic strength (σint) and the corresponding strain (εint) will be derived as [31]:
3. Results and discussions 3.1. DFT calculations of the mechanical behavior of SLMoS2
σint = −E 2 /4D at the strain εint = −E /2D .
The nonlinear elastic constitutive behavior of SLMoS2 can be expressed by a Taylor expansion of the strain energy density u in terms of the powers of strain:
(10)
The same model is used for SLMoS2 in the present work: the values of E, ν and D of SLMoS2 are determined from equations (8) and (9) on the basis of DFT calculations, which will be used as the benchmarks to evaluate the effectiveness 1 1 1 u = Cij εi ε j + Cijk εi ε j εk + Cijkl εi ε j εk εl + … (4) of the empirical potentials selected. It should be noted that the 2! 3! 4! values of C11, C22, C12, C111, C222 and C112 determined from the aforementioned equations might be highly sensitive to the where εx (x = i, j, k and l) are the Lagrange strain, and Cij, Cijk order of the polynomial fitting function as well as the fitting and Cijkl are the second-order, third-order and fourth-order range, which has been shown by the results from graelastic constants, respectively. The subscripts of strain (i, j, k phene [30]. and l) are denoted by the Voigt notations: 11 → 1, 22 → 2, In the present work, the linear elastic parameters (C11, 33 → 3, 23 → 4, 31 → 5 and 12 → 6. For uniaxial/biaxial strain C22 and C12) are determined by equations (5)–(7) after tension, the strain energy density is expressed as: neglecting the third and higher order terms, and the fitting range is selected as 2%. After plugging the determined linear 1 1 ui = Cii εi2 + Ciii εi3 + … (i = 1, 2) parameters into equations (5)–(7) (neglecting the fourth and 2! 3! higher order terms), the nonlinear elastic parameters (C111, (uniaxialstrain tension) (5) C222 and C112) can be fitted from the DFT calculations. To 4
Nanotechnology 26 (2015) 185705
S Xiong and G Cao
Typically, materials are broken at the maximum stress in ′ ), which is different to the DFT MD simulations (i.e., εf = εint calculation result of SLMoS2. Therefore, the DFT benchmark of the failure strain (εf) and stress (σf) for MD simulations are ′ ) and the corresponding stress deterselected as min(εf, εint mined by the DFT calculations, as shown in table 2. 3.2. MD simulations of the elastic behavior of SLMoS2
The geometric parameters (the Mo–Mo distance λMo–Mo, the Mo-S bond length λMo–S, the distance between two S planes δ, the bond angles of θMo–S–Mo and θS–Mo–S) of the equilibrated SLMoS2 structure calculated from MD simulations are displayed in table 1. Compared with the DFT results, all of the empirical potentials give a good prediction in the values of λMo–Mo and λMo–S (the maximum error is about 4%), whereas there is a larger difference in the value of δ and the bond angles θMo–S–Mo and θS–Mo–S (the maximum error is about 9%). Overall, the REBO and SW potentials (the maximum error is roughly less than 5%) give a better prediction than two CVFF potentials. Figures 4 and 5 show the strain energy density u and the σ–ε relationships calculated from the MD simulations. For reference purposes, the DFT results are also displayed as the dashed lines in figures 4 and 5. In MD simulations, the true stress along the loading direction can be directly calculated from the following equation,
Figure 2. Stress–strain (σ–ε) relationships of SLMoS2 under the uniaxial/biaxial strain tension determined from the DFT calculations.
show the strain-softening effect, the fitting strain ranges (ε ⩽ 15%) are selected to determine the nonlinear parameters. After determining the values of E and D, the values of intrinsic strength (σint) and the corresponding strain (εint) can be calculated from equation (10). The determined linear/ nonlinear properties from the DFT calculations are list in table 2. ′ of SLMoS2 determined from The maximum stress σint uniaxial strain tension are quite close along the zz and ac ′ ) along the zz directions, but the corresponding strain value (εint direction is twice that along the ac direction, which is different from the result of graphene [30]. The SLMoS2 structure is still ′ , which is the same as the result of graphene stable with ε ≫ εint [30]. Cooper et al [18] considered that the SLMoS2 will be ′ . Peng et al ′ , i.e., the failure strain εf = εint unstable when σ > σint [15] suggested that although the SLMoS2 structure can still ′ from the σ-ε curve, the resist the tension load with ε > εint structure has been failed on the basis of the geometric parameters of structure. Li et al [20] reported that the failure of the SLMoS2 structure could be caused by the out-of-plane softmode phonon instability. We also checked the change of the geometric parameters of SLMoS2 under the uniaxial/biaxial tension. With the increase of ε (uniaxial/biaxial), the values of the distance between two S atom planes δ gradually decreases but the variation trend is sharply changed at ε = 0.26, 0.42, 0.31 for the biaxial tension, uniaxial tension along both the zz and ac directions, respectively, marked by the dashed lines in figure 3, which are considered as the failure strain (εf) in the present ′ for the uniaxial tension along the zz direction, work. εf ≈ εint ′ for the uniaxial tension along the ac direction whereas εf > εint and the biaxial tension. The values of εf determined in the present work are very close to the results reported by Li et al [20] based on phonon instability (εf = 0.22, 0.42, 0.32). The maximum stresses under the uniaxial/biaxial strain tension determined in the present work are also very close to the reported values [15, 18].
⎡ N 1 Σ = − ⎢ ∑m i vi2 + V ⎢⎣ i = 1
⎤
∑rij fij ⎥⎥ i 10%). The reason is that the errors created by E and D are partially cancelled by each other. Due to the large underestimation of D, the values of σint and εint estimated by equation (10) will be greatly overestimated. In the σ–ε curves calculated from the SW potential, the stress values will be suddenly decreased by 20–30% at some
Figure 3. The DFT results of the variation of the distance between
two S atom planes (δ) with the strain (ε) under the uniaxial/biaxial strain tension. The dashed lines denote the critical strain (εf) corresponding to the structure instability.
underestimate C12, and thus, the CVFF2 can give quite a good estimation of the value of E but there is a huge error in describing the value of ν (as shown in table 2). In addition, the value of D is significantly underestimated, and thus, the values of σint and εint estimated by equation (10) will be greatly overestimated. The differences between two CVFF potentials are mainly from the different equilibrium bond length and bond angles, which cause the different bond stretching and bond angle energies as well as the non-bond energy, as shown in figure 6. CVFF2 has a higher equilibrium Mo-S bond length than CVFF1, which is closer to its true value; CVFF2 considers the 6
Nanotechnology 26 (2015) 185705
S Xiong and G Cao
Figure 4. The strain energy density of SLMoS2 calculated from the empirical potentials under (a) uniaxial strain tension along the zz direction; (b) uniaxial strain tension along the ac direction; (c) biaxial tension.
Figure 5. The σ–ε curves of SLMoS2 calculated from the empirical potentials under (a) uniaxial strain tension along the zz direction; (b) uniaxial strain tension along the ac direction; (c) biaxial tension.
7
Nanotechnology 26 (2015) 185705
S Xiong and G Cao
the REBO potential depends on the empirically selected cutoff interaction distance (e.g., bond length). Under the tensile load, the atomic bond will be stretched, and the bond will be considered to be broken (no bond interaction) when the deformed bond length is larger than the cut-off distance. The failure stresses calculated from the REBO potential are quite close to those calculated from the DFT calculations, whereas the failure strains are far less than their DFT benchmarks (as shown in figure 5). In addition, the REBO potential gives a larger value of εf along the ac direction than that along the zz direction, which is also opposite to the DFT results. In the REBO potential, there are three atomic interactions (Mo–Mo, S–S and Mo–S) for SLMoS2 (listed as CASE A in table 3), in which the Mo–Mo cut-off distance is firstly touched under in-plane tension, and thus, it is the main control factor [16]. This leads to a larger εf along the ac direction than that along the zz direction because the Mo–Mo vector is parallel to the ac direction, and there is an angle of 60° between the Mo–Mo vector and the zz direction. This failure mechanism is different than that shown by the DFT calculations of SLMoS2. In DFT calculations, the Mo–Mo distance gradually increases with ε, which will not cause failure of the structure; whereas the structure is failed by the failure of out-of-plane soft-mode phonon instability [20], which can be further supported by the sudden change of the S–S plane distance δ (figure 3). Therefore, we can modify the cut-off distances of Mo–Mo, S–S and Mo–S interactions in the REBO potential on the basis of the failure strains determined by the DFT calculations, as shown in figure 7. The new cut-off distances of Mo–Mo, S–S and Mo–S interactions are listed in table 3, and the REBO potential with the modified cut-off distances is defined as the REBO* potential in the present work. After modification, the Mo–S interaction becomes the control factor of the failure of SLMoS2; the values of εf = 0.24, 0.31 and 0.23 for the biaxial tension, the uniaxial tension along the zz and ac directions, respectively, which are quite close to their DFT counterparts. However, the modification of the cut-off distances also causes the increase of the failure stress (σf = 17.5 N m−1 for the biaxial and the uniaxial tension along the ac and zz directions, respectively), which are about 20–40% higher than their DFT counterparts. Stewart et al investigated the indentation behavior of crystalline MoS2 by using the REBO potential [16]. The present work is the first attempt to employ the REBO potential to study the tensile behavior of SLMoS2. Compared with the results of graphene, the reported σf and εf along the ac and zz directions calculated from the REBO potential and the DFT calculations are: σf = 25.9, 28.1 N m−1 and εf = 0.2, 0.33 [32], σf = 30.3, 31.1 N m−1 and εf = 0.21, 0.30 [33], respectively. The REBO potential underestimates the value of σf by about 20%. In addition, the elastic modulus of graphene is also underestimated by 30% and the Poisson ratio is overestimated by about 200% from the REBO potential [32]. After considering the more complicated structure of SLMoS2 than graphene, we think that the REBO potential provides a pretty good description of the mechanical behavior of SLMoS2. The higher values of σf are caused by the
Figure 6. The bond stretching energy (Ubond), bond angle variation energy (Uangle) and nonbond energy (ULJ + Uc) of two CVFF potentials.
critical value of strain, and then, the stress continually increases with the strain. The sudden change of the value of σ can be considered to be caused by the failure of the structure of SLMoS2, and thus, the SW potential can be used to describe the failure behavior of the SLMoS2. The maximum value of σ before the sudden decrease is defined as the failure stress σf and the corresponding strain is called the failure strain εf. The values of σf are very close to their DFT counterparts, while the values of εf are much lower than the corresponding values determined from the DFT calculations, as shown in table 2. In addition, the value of εf along the ac direction is higher than that along the zz direction calculated by the SW potential, which is opposite to the DFT results. Currently, the reported results of SLMoS2 based on the SW potential are focused on the phonon spectrum calculation [13], and the present work is the first attempt to use this potential to investigate the tensile response of SLMoS2. 3.2.3. REBO potential. The REBO potential gives a similar
elastic behavior of SLMoS2 as CVFF2. It can accurately describe the elastic behavior under a small deformation (e.g., ε < 5%) and the overestimation of σ gradually increases with ε: the error in both E and ν are about 5%, whereas the value of σ will be overestimated by about 20% at ε = 10%. The reason is that the nonlinear behavior is underestimated by the REBO potential, which can be seen from the underestimation of the value of D (50% lower than that calculated from DFT). Due to the great underestimation of D, the values of σint and εint estimated by equation (10) will be greatly overestimated. The REBO potential is not only used to study the elastic behavior of materials but also to simulate the failure (breaking) behavior of materials because it can consider the bond break and reform, which is a significant advantage over other empirical potentials. The failure strain determined by 8
Nanotechnology 26 (2015) 185705
S Xiong and G Cao
Table 3. The cut-off distances selected in the REBO potential.
Cut-off distance (nm)
Mo–Mo (Rmin, Rmax)
S–S (Rmin, Rmax)
Mo–S (Rmin, Rmax)
Aa B Cb D
3.50, 4.05, 4.05, 4.05,
2.30, 2.30, 2.30, 2.30,
2.75 3.05 2.65, 2.70 2.68, 2.73 2.70, 2.75
3.80 4.35 4.35 4.35
3.00 3.00 3.00 3.00
a b
The parameters used in the original REBO potential. The parameters used in the modified REBO potential (REBO*).
underestimated nonlinear behavior (strain-softening effect) in the REBO potential.
4. Conclusions The effectiveness of four empirical potential functions, including CVFF1, CVFF2, SW and REBO potentials, on describing the mechanical behavior of SLMoS2 are evaluated based on the benchmark created by the DFT calculations. To the best of our knowledge, this is the first attempt to employ MD simulations to investigate the tensile response of SLMoS2. Among four empirical potentials, it is found that the CVFF2 and REBO potentials can provide a good description of the elastic behavior of SLMoS2 under small deformation (ε < 5%); CVFF2 can accurately predict the stress value with a very large deformation under biaxial tension; after the cut-off distance modification, the REBO potential gives a good prediction for the failure strain, but the failure stress will be overestimated by about 20–40% under biaxial and uniaxial tension. The REBO potential gives a good prediction of the equilibrium structure of SLMoS2, and the deviations of all geometric parameters are less than 5%, among which there is only 1% difference for both λMo–Mo and λMo–S. The REBO potential can give a good estimation for both the elastic modulus E and the Poisson’s ratio ν of SLMoS2 (about 5% overestimation for both E and ν); it underestimates the strainsoftening effect (D, the nonlinear elastic modulus) of SLMoS2, and thus, it predicts a higher stress value at a given strain ε than the DFT. Compared with the results of graphene, it is concluded that the REBO potential provides a pretty good description of the overall mechanical behavior of SLMoS2 considering the complicated structure of SLMoS2. Both CVFF1 and SW potentials cannot correctly describe the value of E (the former overestimate E by 46% and the latter underestimate E by 40%). However, the SW potential can give a more accurate description of the stress value under the uniaxial or biaxial tension with a large deformation (e.g., ε > 10%) because the errors created in E and D are partially cancelled with each other. In all MD simulations (based on different potentials), the values of σint and εint will be
Figure 7. The effect of the cut-off distances in the REBO potential
on the σ–ε curves of SLMoS2; (a) uniaxial strain tension along the ac direction; (b) uniaxial strain tension along the zz direction; (c) biaxial tension. The parameters of CASEs (A)–(D) are displayed in table 3.
9
Nanotechnology 26 (2015) 185705
S Xiong and G Cao
significantly overestimated by the two-parameter nonlinear elastic model (equation (10)) due to the underestimation of D. Although DFT calculations can provide a more accurate description of the mechanical behavior of SLMoS2, they cannot be used to investigate the structure with a large size or under a complicated loading condition (e.g., free standing indentation). Therefore, the present study will provide a useful help to theoretically investigate the mechanical behavior of SLMoS2 as well as to understand the testing data in experiments.
[14] [15] [16]
[17]
Acknowledgments [18]
I acknowledge the financial support provided by the Ministry of Science and Technology of China (grant no. 2013CB933702) and the National Science Foundation of China (grant no. 11172002).
[19] [20] [21]
References
[22]
[1] Mak K F et al 2010 Atomically thin MoS2: a new direct-gap semiconductor Phys. Rev. Lett. 105 136805 [2] Mak K F et al 2014 The valley hall effect in MoS2 transistors Science 344 1489–92 [3] Radisavljevic B et al 2011 Single-layer MoS2 transistors Nat. Nanotechnology 6 147–50 [4] Deng Y X et al 2014 Black phosphorus-monolayer MoS2 van der waals heterojunction p-n diode ACS Nano 8 8292–9 [5] Liu H et al 2014 Switching mechanism in single-layer molybdenum disulfide transistors: an insight into current flow across schottky barriers ACS Nano 8 1031–8 [6] Mak K F et al 2012 Control of valley polarization in monolayer MoS2 by optical helicity Nat. Nanotechnology 7 494–8 [7] Li Y L et al 2013 Probing symmetry properties of few-layer MoS2 and h-BN by optical second-harmonic generation Nano Lett. 13 3329–33 [8] Bertolazzi S, Krasnozhon D and Kis A 2013 Nonvolatile memory cells based on MoS2/graphene heterostructures ACS Nano 7 3246–52 [9] Chhowalla M and Amaratunga G A J 2000 Thin films of fullerene-like MoS2 nanoparticles with ultra-low friction and wear Nature 407 164–7 [10] Liu K et al 2014 Elastic properties of chemical-vapordeposited monolayer MoS2, WS2, and their bilayer heterostructures Nano Lett. 14 5097–103 [11] Castellanos-Gomez A et al 2012 Elastic properties of freely suspended MoS2 nanosheets Adv. Mater. 24 772 -+ [12] Bertolazzi S, Brivio J and Kis A 2011 Stretching and breaking of ultrathin MoS2 ACS Nano 5 9703–9 [13] Jiang J W, Park H S and Rabczuk T 2013 Molecular dynamics simulations of single-layer molybdenum disulphide (MoS2): stillinger-weber parametrization,
[23] [24] [25] [26] [27] [28] [29] [30] [31] [32]
[33]
10
mechanical properties, and thermal conductivity J. Appl. Phys. 114 064307 Jiang J W, Qi Z N and Park H S 2013 Elastic bending modulus of single-layer molybdenum disulfide (MoS2): finite thickness effect Nanotechnology 24 435705 Peng Q and De S 2013 Outstanding mechanical properties of monolayer MoS2 and its application in elastic energy storage Phys. Chem. Chem. Phys. 15 19427–37 Stewart J A and Spearot D E 2013 Atomistic simulations of nanoindentation on the basal plane of crystalline molybdenum disulfide (MoS2) Modelling Simul. Mater. Sci. Eng. 21 045003 Zhao J H, Jiang J W and Rabczuk T 2013 Temperaturedependent mechanical properties of single-layer molybdenum disulphide: molecular dynamics nanoindentation simulations Appl. Phys. Lett. 103 231913 Cooper R C et al 2013 Nonlinear elastic behavior of twodimensional molybdenum disulfide Phys. Rev. B 87 035423 Jiang J W 2014 The buckling of single-layer MoS2 under uniaxial compression Nanotechnology 25 355402 Li T S 2012 Ideal strength and phonon instability in singlelayer MoS2 Phys. Rev. B 85 235407 Wakabayashi N, Smith H G and Nicklow R M 1975 Latticedynamics of hexagonal Mos2 studied by neutron-scattering Phys. Rev. B 12 659–63 Morita Y et al 2008 Development of a new molecular dynamics method for tribochemical reaction and its application to formation dynamics of MoS(2) tribofilm Appl. Surf. Sci. 254 7618–21 Becker U et al 2003 Metal island growth and dynamics on molybdenite surfaces Geochim. Cosmochim. Acta 67 923–34 Varshney V et al 2010 MD simulations of molybdenum disulphide (MoS2): force-field parameterization and thermal transport behavior Comput. Mater. Sci. 48 101–8 Liang T, Phillpot S R and Sinnott S B 2009 Parametrization of a reactive many-body potential for Mo–S systems Phys. Rev. B 79 245110 Kresse G and Furthmuller J 1996 Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set Phys. Rev. B 54 11169–86 Kresse G and Furthmuller J 1996 Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set Comput. Mater. Sci. 6 15–50 Perdew J P, Burke K and Ernzerhof M 1996 Generalized gradient approximation made simple Phys. Rev. Lett. 77 3865–8 Plimpton S 1995 Fast parallel algorithms for short-range molecular dynamics J. Comput. Phys. 117 1–19 Cao G X 2014 Atomistic studies of mechanical properties of graphene Polymers 6 2404–32 Lee C et al 2008 Measurement of the elastic properties and intrinsic strength of monolayer graphene Science 321 385–8 Lu Q, Gao W and Huang R 2011 Atomistic simulation and continuum modeling of graphene nanoribbons under uniaxial tension Modelling Simul. Mater. Sci. Eng. 19 054006 Liu F, Ming P M and Li J 2007 Ab initio calculation of ideal strength and phonon instability of graphene under tension Phys. Rev. B 76 064120
Search
Collections
Journals
About
Contact us
My IOPscience
Molecular dynamics simulations of mechanical properties of monolayer MoS2
This content has been downloaded from IOPscience. Please scroll down to see the full text. 2015 Nanotechnology 26 185705 (http://iopscience.iop.org/0957-4484/26/18/185705) View the table of contents for this issue, or go to the journal homepage for more
Download details: IP Address: 150.216.68.200 This content was downloaded on 25/04/2015 at 17:08
Please note that terms and conditions apply.
Nanotechnology Nanotechnology 26 (2015) 185705 (10pp)
doi:10.1088/0957-4484/26/18/185705
Molecular dynamics simulations of mechanical properties of monolayer MoS2 Si Xiong and Guoxin Cao HEDPS, Center for Applied Physics and Technology, Department of Mechanics and Engineering Science, College of Engineering, Peking University, Beijing, 100871, People’s Republic of China E-mail: [email protected] Received 14 February 2015, revised 18 March 2015 Accepted for publication 20 March 2015 Published 16 April 2015 Abstract
Based on the benchmark created by the first principle calculations, the effectiveness of different empirical potential functions, including CVFF1, CVFF2, SW and REBO potentials, on describing the mechanical behavior of single layer MoS2 (SLMoS2) are evaluated. The mechanical properties including the elastic modulus E, Poisson’s ratio ν, nonlinear elastic modulus D, failure stress and ultimate strain are considered. It is found that under a small deformation, the REBO and CVFF2 potentials can provide an effective description of the elastic behavior of SLMoS2; whereas under a large deformation, the SW potential gives a more accurate stress value and the CVFF2 potential can only predict the right value under biaxial tension. After modifying the cut-off distances in the REBO potential, the failure strain can be reasonably predicted by the REBO potential; whereas the failure stress will be overestimated by about 20–40% for biaxial and uniaxial tension. The present study can provide help on accurately understanding the mechanical behavior of SLMoS2. Keywords: mechanical properties, single layer MoS2, molecular dynamics (Some figures may appear in colour only in the online journal) 1. Introduction
between the AFM load and the deformation of the freestanding structure. Besides experimental work, the reports of the mechanical properties of layered MoS2 are also developed theoretically [13–20]. Under applied in-plane tension, the nonlinear elastic properties along with the intrinsic strength/ultimate strain of SLMoS2 are reported from first-principles calculations [15, 18, 20]. Although first-principles calculations are the most accurate theoretical approach, they have the following limitations: (a) they can only simulate several hundred atoms; (b) more complicated mechanical behavior, including indentation response and bending behavior, cannot be investigated. Therefore, molecular dynamics (MD) simulations based on the empirical potentials (force field) are desirable. The accuracy of MD simulations mainly depends on the effectiveness of empirical potential function selected. Wakabayashi et al [21] developed a constant valence force field (CVFF) type potential, in which the bond interactions between Mo and S atoms (including both bond stretching and bond angle variation terms) were described by harmonic potentials. Morita et al [22] and Becker et al [23] provided the revised empirical
SLMoS2 has a direct band gap of 1.8 eV [1], which makes it a good candidate for applications in nanoelectronics, such as field-effect transistors (FET) [2–5], phototransistors [6, 7], nonvolatile memory cells [8] and so on. In addition, MoS2 can be used as a solid lubricant [9]. To effectively design structures for its practical applications, the mechanical properties of single layer MoS2 (SLMoS2) should be well understood. Since SLMoS2 is typically in the atomic-level thick range, its mechanical properties are highly difficult to measure by the conventional experimental methods. Currently, the mechanical properties of 2D materials, such as graphene and layered MoS2, are mainly measured from its free-standing structure under clamped boundary conditions [10–12]. 2D materials are firstly mounted on the substrate with cylindrical wells (holes) (called a free-standing structure), and then deformed by the tip of an atomic force microscope (AFM). The mechanical properties, including elastic modulus and intrinsic strength, can be estimated from the relationship 0957-4484/15/185705+10$33.00
1
© 2015 IOP Publishing Ltd Printed in the UK
Nanotechnology 26 (2015) 185705
S Xiong and G Cao
potential function formula. Both of them used the Morse potential to describe the Mo–S bond stretching interaction and the non-bond interactions between Mo and S atoms were considered by both the Lennard–Jones (LJ) potential and the electrostatic interaction. By comparing with the experimental values of MoS2 lattice parameters, Varshney et al [24] corrected the parameters in the aforementioned empirical potential. Liang et al [25] parameterized an interatomic potential for the Mo–S system based on the second-generation reactive empirical bond-order (REBO) formalism. Jiang et al [13] presented a parameterization of the Stillinger-Weber (SW) potential to describe the interatomic interactions within single-layer MoS2. Based on the developed empirical potentials, more mechanical behaviors of MoS2 were investigated using MD simulations: such as indentation response [16], buckling behavior [19], bending behavior [14] and temperature effect[17]. However, the effectiveness of the aforementioned empirical potentials on investigating the mechanical behavior of layered MoS2 has not been fully evaluated. In the reported results, the experimentally obtained Mo–S lattice parameters [24] and phonon spectrum [13] of bulk MoS2 were selected as the benchmark of the effectiveness of empirical potential. However, it is still not clear whether or not these empirical potentials can effectively describe the deformation of layered MoS2 or which one among these empirical potentials is most effective for studying the mechanical behavior of layered MoS2, especially for the nonlinear elastic behavior or the failure behavior. In the present work, first-principles calculations are firstly used to investigate the tensile properties of SLMoS2, including elastic modulus, Poisson’s ratio, nonlinear elastic modulus, failure stress and ultimate strain. Then, the similar tensile behavior is studied by MD simulations based on different empirical potentials, including CVFF, SW and REBO potentials. Finally, the effectiveness of different empirical potentials is analyzed for describing the mechanical behavior of SLMoS2. The present study can provide a useful guideline to understand the mechanical properties of SLMoS2 and to design the structure of layered MoS2 in its practical applications.
of SLMoS2 determined from our DFT calculations, the spacing between two Mo atoms λMo–Mo = 0.3181 nm, the distance between two S atom planes δ = 0.3114 nm and the Mo–S bond length λMo–S = 0.2408 nm, which are in good agreement with the reported simulation results [15, 18] and the experimental results of bulk MoS2 crystal (λMo–Mo = 0.3162 nm and δ = 0.318 nm in bulk) [12]. The geometric parameters of the equilibrium structure SLMoS2 determined by DFT calculations are displayed in table 1. The uniaxial/biaxial strain tension is applied by increasing the computational cell size along the loading direction.
2.2. MD simulations of SLMoS2
SLMoS2 has a tri-layer structure (figures 1(c) and (d)). Each Mo atom locates in the center of the triangular prism formed by six neighboring S atoms and the S–Mo–S bond is mainly covalent in nature. In MD simulations, different empirical potentials are selected to optimize the equilibrium structure of SLMoS2 as well as its deformed structure under in-plane tension: including the CVFF potential [24], SW potential [13] and REBO potential [16]. The force field parameters for all different empirical potentials used in the present work are the same as those reported in the original literatures, unless specifically stated that they were modified. MD simulations are carried out using the LAMMPS package [29]. A periodic boundary condition is applied in the lateral directions to remove the lateral boundary effect and simulate the intrinsic properties of SLMoS2. For all simulations, the computational cell includes more than 900 atoms and the size is about 4.74 × 5.47 × 9 nm, and the lateral size is slightly varying with the empirical potentials selected. The simulation results are not sensitive to the selected computational cell size, which has been checked using a larger computational cell size. The uniaxial/biaxial strain tension is applied on the computational cell using the NVT ensemble. The strain rate is kept as 109/s for all cases. The system is kept as 1 K and a constant time step of 1 fs is selected.
2.2.1. CVFF potentials. In the CVFF potential, the potential
2. Computational methods
function can be expressed as the summation of the bondstretching energy (Ubond), the bond angle variation energy (Uangle), the LJ non-bond energy (ULJ) and the Coulomb electrostatic energy (UC) [24]:
2.1. First-principles calculations of 2D SLMoS2
In the present work, the mechanical response of SLMoS2 under the uniaxial/biaxial tension is firstly investigated by density functional theory (DFT) calculations, which are performed using the Vienna Ab Initio Simulation Package (VASP) [26, 27]. DFT calculations are performed based on the six-atom (including two Mo atoms and four S atoms) SLMoS2 unit cell, and the unit cell size is 3.1809 × 5.5096 × 61.5 Å (shown in figures 1(a) and (b)). The generalized gradient approximation is applied using the Perdew–Burke–Ernzerhof functional [28]. A Monkhorst–Pack k-point mesh grid 18 × 18 × 3 is selected, and the plane wave energy cut-off is set at 500 eV. For the equilibrium structure
U=
∑Ubond + ∑Uangle + ∑ULJ + ∑Uc
Ubond = K1 × ( rij − r0 )2 (Harmonic bond potential) 2 Ubond = D × ⎡⎣ 1 − e−α ( rij− r0 ) ⎤⎦ (Morse bond potential)
(
)2
Uangle = K2 × θ ijk − θ0 (Harmonic angle potential) ⎡ ⎛ ⎞12 ⎛ σ ⎞6 ⎤ σ ULJ = ε⎢ ⎜⎜ ⎟⎟ − 2 ⎜⎜ ⎟⎟ ⎥ (LJ nonbond potential) ⎢ ⎝ rij ⎠ ⎝ rij ⎠ ⎥⎦ ⎣ 2
Nanotechnology 26 (2015) 185705
S Xiong and G Cao
Figure 1. Atomic structure of SLMoS2, (a) (b) the computational cell used in the DFT calculations; (c) (d) the top and side view of the computational model used in MD simulations. Table 1. Geometric parameters of SLMoS2.
Parameters λMo–Mo (nm) λMo–S (nm) δ (nm) θMo–S–Mo (°) θS–Mo–S (°)
Uc = C
qi q j rij
DFT
CVFF1
CVFF2
SW
REBO
0.3181 0.2408 0.3114 82.688 80.582
0.3267 0.2365 0.2852 87.384 74.189
0.3169 0.2313 0.2831 86.471 75.447
0.3064 0.2398 0.3238 79.417 84.927
0.3146 0.2438 0.3253 80.357 83.687
(Coulomb electrostatic potential)
2.2.2.
SW potential. In the SW potential, the bond interaction is described by both a two-body interaction (U2(i, j)) and a three-body interaction (U3(i, j, k)) [13]:
(1)
where K1, K2, D are the bond potential parameters, σ and ε are the LJ non-bond potential parameters and C is the Coulomb electrostatic potential parameter. rij is the distance between atoms i and j; θijk is the angle between bond ij and bond jk; qi and qj are the partial charge values for atoms i and j; r0 is the equilibrium bond length between atoms i and j; and θ0 is the equilibrium angle between bond ij and bond jk. The bond stretching potential can be expressed by both the harmonic potential and the Morse potential [24], which are defined as CVFF1 and CVFF2, respectively, in the present work. The potential parameters used in equation (1) are taken from Varshney et al [24].
U=
∑U2 (i ,
j) + ⎡
∑U3 (i , −1
U2 (i , j ) = S1e⎣ ρ ( rij− rmax ) U3 (i , j , k ) =
(
⎤ ⎦
(S r
⎡ S3 e⎣ ρ1 ( rij− rij _
× cos θ jik − cos θ0
j, k )
)2
−4 2 ij
max
)−1 + ρ
−1
)
2 ( rjk − rjk _ max
)−1 ⎤⎦
(2)
where the parameters S1, S2, S3, ρ, ρ1, ρ2, rmax, rij_max and rjk_max are taken from Jiang et al [13].
3
Nanotechnology 26 (2015) 185705
S Xiong and G Cao
2.2.3. REBO potential. In the REBO potential, the bond
u=
interaction is described by a many-body function [16, 25]: 1 ∑ fijc ( rij ) ⎡⎣ V R ( rij) − bij V A ( rij ) ⎤⎦ 2 ⎛ Qij ⎞ ⎟⎟ V R ( rij ) = Aij e−αij rij ⎜⎜ 1 + rij ⎠ ⎝
Ub =
(6)
The stress–strain relationship can be calculated by the derivative of u with respect to ε:
V A ( rij ) = Bij e−βij rij bij = ⎡⎣ 1 +
∑
fikc
1 Ciii εi2 + … 2 (unaxial strain tension)
σi = Cii εi +
−1/2 ( rik ) ⋅ G ⎡⎣ cos θ ijk ⎤⎦ + P ( Ni ) ⎤⎦
( )
P ( Ni ) = −a 0 ( Ni − 1) − a1e−a 2 Ni + a 3 Mo
Ni =
(i = 1, 2)
σb = ( C11 + C22 + 2C12 )
S
εi + ( 2C111 − C222 + 3C112) εi2 + … (biaxial strain tension)
∑ fikc ( rik ) + ∑ filc ( ril )
i≠k fijc ( rij )
1 ( C11 + C22 + 2C12 ) εi2 2 1 + ( 2C111 − C222 + 3C112 ) εi3 + … (biaxial tension) 3
(7)
i≠l
⎧1 ⎪ ⎪ ⎧ ⎡ π r − R min ⎪1⎪ ij ij ⎨ = ⎨ 1 + cos ⎢ ⎢ max min ⎪2⎪ ⎣ Rij − Rij ⎪ ⎩ ⎪ ⎩0
(
(
) ⎤⎥ ⎫⎪⎬ ) ⎥⎦ ⎪⎭
rij ⩽
where i = 1, 2 represents the zigzag (zz) and armchair (ac) directions, respectively. Figure 2 shows the calculated values of σi and σb/2, which are very close to the reported results [15, 18]. Similar to graphene, the σ–ε relationships are almost the same along the zz and ac direction when ε < 10%, which means that SLMoS2 can be considered to be isotropic with ε < 10%. Theoretically, the isotropic, linear elastic properties: elastic modulus E and Poisson’s ratio ν can be calculated as:
Rijmin
Rijmin ⩽ rij ⩽ Rijmax Rijmax ⩽ rij (3)
(
)
2 2 E = C11 − C12 / C11 and ν = C12 / C11.
where VR and VA are the pairwise repulsive and attractive interactions, bij is the many-body bond-order function, Ni is the total number of nearest neighbors within the maximum REBO cutoff distance and fijc is the cutoff function that limits the REBO interactions to nearest neighbors. The pairwise parameters α, β, A, B, Q, Rijmin , Rijmax and a0–a1 depend on the chemical species of the interacting atoms i and j, which are taken from Stewart et al [16].
(8)
For graphene, the isotropic, nonlinear elastic modulus D is introduced to describe the nonlinear elastic behavior under a large tension [30, 31], and thus, the stress-strain relationship can be simplified as: σ = Eε + Dε2 1 D= C111 1+3ν 2 − C222 ν 3 + 3ν 2 2
{ (
)
(
+ 3C112 ν 2 − ν
(
)
)}
(9)
The intrinsic strength (σint) and the corresponding strain (εint) will be derived as [31]:
3. Results and discussions 3.1. DFT calculations of the mechanical behavior of SLMoS2
σint = −E 2 /4D at the strain εint = −E /2D .
The nonlinear elastic constitutive behavior of SLMoS2 can be expressed by a Taylor expansion of the strain energy density u in terms of the powers of strain:
(10)
The same model is used for SLMoS2 in the present work: the values of E, ν and D of SLMoS2 are determined from equations (8) and (9) on the basis of DFT calculations, which will be used as the benchmarks to evaluate the effectiveness 1 1 1 u = Cij εi ε j + Cijk εi ε j εk + Cijkl εi ε j εk εl + … (4) of the empirical potentials selected. It should be noted that the 2! 3! 4! values of C11, C22, C12, C111, C222 and C112 determined from the aforementioned equations might be highly sensitive to the where εx (x = i, j, k and l) are the Lagrange strain, and Cij, Cijk order of the polynomial fitting function as well as the fitting and Cijkl are the second-order, third-order and fourth-order range, which has been shown by the results from graelastic constants, respectively. The subscripts of strain (i, j, k phene [30]. and l) are denoted by the Voigt notations: 11 → 1, 22 → 2, In the present work, the linear elastic parameters (C11, 33 → 3, 23 → 4, 31 → 5 and 12 → 6. For uniaxial/biaxial strain C22 and C12) are determined by equations (5)–(7) after tension, the strain energy density is expressed as: neglecting the third and higher order terms, and the fitting range is selected as 2%. After plugging the determined linear 1 1 ui = Cii εi2 + Ciii εi3 + … (i = 1, 2) parameters into equations (5)–(7) (neglecting the fourth and 2! 3! higher order terms), the nonlinear elastic parameters (C111, (uniaxialstrain tension) (5) C222 and C112) can be fitted from the DFT calculations. To 4
Nanotechnology 26 (2015) 185705
S Xiong and G Cao
Typically, materials are broken at the maximum stress in ′ ), which is different to the DFT MD simulations (i.e., εf = εint calculation result of SLMoS2. Therefore, the DFT benchmark of the failure strain (εf) and stress (σf) for MD simulations are ′ ) and the corresponding stress deterselected as min(εf, εint mined by the DFT calculations, as shown in table 2. 3.2. MD simulations of the elastic behavior of SLMoS2
The geometric parameters (the Mo–Mo distance λMo–Mo, the Mo-S bond length λMo–S, the distance between two S planes δ, the bond angles of θMo–S–Mo and θS–Mo–S) of the equilibrated SLMoS2 structure calculated from MD simulations are displayed in table 1. Compared with the DFT results, all of the empirical potentials give a good prediction in the values of λMo–Mo and λMo–S (the maximum error is about 4%), whereas there is a larger difference in the value of δ and the bond angles θMo–S–Mo and θS–Mo–S (the maximum error is about 9%). Overall, the REBO and SW potentials (the maximum error is roughly less than 5%) give a better prediction than two CVFF potentials. Figures 4 and 5 show the strain energy density u and the σ–ε relationships calculated from the MD simulations. For reference purposes, the DFT results are also displayed as the dashed lines in figures 4 and 5. In MD simulations, the true stress along the loading direction can be directly calculated from the following equation,
Figure 2. Stress–strain (σ–ε) relationships of SLMoS2 under the uniaxial/biaxial strain tension determined from the DFT calculations.
show the strain-softening effect, the fitting strain ranges (ε ⩽ 15%) are selected to determine the nonlinear parameters. After determining the values of E and D, the values of intrinsic strength (σint) and the corresponding strain (εint) can be calculated from equation (10). The determined linear/ nonlinear properties from the DFT calculations are list in table 2. ′ of SLMoS2 determined from The maximum stress σint uniaxial strain tension are quite close along the zz and ac ′ ) along the zz directions, but the corresponding strain value (εint direction is twice that along the ac direction, which is different from the result of graphene [30]. The SLMoS2 structure is still ′ , which is the same as the result of graphene stable with ε ≫ εint [30]. Cooper et al [18] considered that the SLMoS2 will be ′ . Peng et al ′ , i.e., the failure strain εf = εint unstable when σ > σint [15] suggested that although the SLMoS2 structure can still ′ from the σ-ε curve, the resist the tension load with ε > εint structure has been failed on the basis of the geometric parameters of structure. Li et al [20] reported that the failure of the SLMoS2 structure could be caused by the out-of-plane softmode phonon instability. We also checked the change of the geometric parameters of SLMoS2 under the uniaxial/biaxial tension. With the increase of ε (uniaxial/biaxial), the values of the distance between two S atom planes δ gradually decreases but the variation trend is sharply changed at ε = 0.26, 0.42, 0.31 for the biaxial tension, uniaxial tension along both the zz and ac directions, respectively, marked by the dashed lines in figure 3, which are considered as the failure strain (εf) in the present ′ for the uniaxial tension along the zz direction, work. εf ≈ εint ′ for the uniaxial tension along the ac direction whereas εf > εint and the biaxial tension. The values of εf determined in the present work are very close to the results reported by Li et al [20] based on phonon instability (εf = 0.22, 0.42, 0.32). The maximum stresses under the uniaxial/biaxial strain tension determined in the present work are also very close to the reported values [15, 18].
⎡ N 1 Σ = − ⎢ ∑m i vi2 + V ⎢⎣ i = 1
⎤
∑rij fij ⎥⎥ i 10%). The reason is that the errors created by E and D are partially cancelled by each other. Due to the large underestimation of D, the values of σint and εint estimated by equation (10) will be greatly overestimated. In the σ–ε curves calculated from the SW potential, the stress values will be suddenly decreased by 20–30% at some
Figure 3. The DFT results of the variation of the distance between
two S atom planes (δ) with the strain (ε) under the uniaxial/biaxial strain tension. The dashed lines denote the critical strain (εf) corresponding to the structure instability.
underestimate C12, and thus, the CVFF2 can give quite a good estimation of the value of E but there is a huge error in describing the value of ν (as shown in table 2). In addition, the value of D is significantly underestimated, and thus, the values of σint and εint estimated by equation (10) will be greatly overestimated. The differences between two CVFF potentials are mainly from the different equilibrium bond length and bond angles, which cause the different bond stretching and bond angle energies as well as the non-bond energy, as shown in figure 6. CVFF2 has a higher equilibrium Mo-S bond length than CVFF1, which is closer to its true value; CVFF2 considers the 6
Nanotechnology 26 (2015) 185705
S Xiong and G Cao
Figure 4. The strain energy density of SLMoS2 calculated from the empirical potentials under (a) uniaxial strain tension along the zz direction; (b) uniaxial strain tension along the ac direction; (c) biaxial tension.
Figure 5. The σ–ε curves of SLMoS2 calculated from the empirical potentials under (a) uniaxial strain tension along the zz direction; (b) uniaxial strain tension along the ac direction; (c) biaxial tension.
7
Nanotechnology 26 (2015) 185705
S Xiong and G Cao
the REBO potential depends on the empirically selected cutoff interaction distance (e.g., bond length). Under the tensile load, the atomic bond will be stretched, and the bond will be considered to be broken (no bond interaction) when the deformed bond length is larger than the cut-off distance. The failure stresses calculated from the REBO potential are quite close to those calculated from the DFT calculations, whereas the failure strains are far less than their DFT benchmarks (as shown in figure 5). In addition, the REBO potential gives a larger value of εf along the ac direction than that along the zz direction, which is also opposite to the DFT results. In the REBO potential, there are three atomic interactions (Mo–Mo, S–S and Mo–S) for SLMoS2 (listed as CASE A in table 3), in which the Mo–Mo cut-off distance is firstly touched under in-plane tension, and thus, it is the main control factor [16]. This leads to a larger εf along the ac direction than that along the zz direction because the Mo–Mo vector is parallel to the ac direction, and there is an angle of 60° between the Mo–Mo vector and the zz direction. This failure mechanism is different than that shown by the DFT calculations of SLMoS2. In DFT calculations, the Mo–Mo distance gradually increases with ε, which will not cause failure of the structure; whereas the structure is failed by the failure of out-of-plane soft-mode phonon instability [20], which can be further supported by the sudden change of the S–S plane distance δ (figure 3). Therefore, we can modify the cut-off distances of Mo–Mo, S–S and Mo–S interactions in the REBO potential on the basis of the failure strains determined by the DFT calculations, as shown in figure 7. The new cut-off distances of Mo–Mo, S–S and Mo–S interactions are listed in table 3, and the REBO potential with the modified cut-off distances is defined as the REBO* potential in the present work. After modification, the Mo–S interaction becomes the control factor of the failure of SLMoS2; the values of εf = 0.24, 0.31 and 0.23 for the biaxial tension, the uniaxial tension along the zz and ac directions, respectively, which are quite close to their DFT counterparts. However, the modification of the cut-off distances also causes the increase of the failure stress (σf = 17.5 N m−1 for the biaxial and the uniaxial tension along the ac and zz directions, respectively), which are about 20–40% higher than their DFT counterparts. Stewart et al investigated the indentation behavior of crystalline MoS2 by using the REBO potential [16]. The present work is the first attempt to employ the REBO potential to study the tensile behavior of SLMoS2. Compared with the results of graphene, the reported σf and εf along the ac and zz directions calculated from the REBO potential and the DFT calculations are: σf = 25.9, 28.1 N m−1 and εf = 0.2, 0.33 [32], σf = 30.3, 31.1 N m−1 and εf = 0.21, 0.30 [33], respectively. The REBO potential underestimates the value of σf by about 20%. In addition, the elastic modulus of graphene is also underestimated by 30% and the Poisson ratio is overestimated by about 200% from the REBO potential [32]. After considering the more complicated structure of SLMoS2 than graphene, we think that the REBO potential provides a pretty good description of the mechanical behavior of SLMoS2. The higher values of σf are caused by the
Figure 6. The bond stretching energy (Ubond), bond angle variation energy (Uangle) and nonbond energy (ULJ + Uc) of two CVFF potentials.
critical value of strain, and then, the stress continually increases with the strain. The sudden change of the value of σ can be considered to be caused by the failure of the structure of SLMoS2, and thus, the SW potential can be used to describe the failure behavior of the SLMoS2. The maximum value of σ before the sudden decrease is defined as the failure stress σf and the corresponding strain is called the failure strain εf. The values of σf are very close to their DFT counterparts, while the values of εf are much lower than the corresponding values determined from the DFT calculations, as shown in table 2. In addition, the value of εf along the ac direction is higher than that along the zz direction calculated by the SW potential, which is opposite to the DFT results. Currently, the reported results of SLMoS2 based on the SW potential are focused on the phonon spectrum calculation [13], and the present work is the first attempt to use this potential to investigate the tensile response of SLMoS2. 3.2.3. REBO potential. The REBO potential gives a similar
elastic behavior of SLMoS2 as CVFF2. It can accurately describe the elastic behavior under a small deformation (e.g., ε < 5%) and the overestimation of σ gradually increases with ε: the error in both E and ν are about 5%, whereas the value of σ will be overestimated by about 20% at ε = 10%. The reason is that the nonlinear behavior is underestimated by the REBO potential, which can be seen from the underestimation of the value of D (50% lower than that calculated from DFT). Due to the great underestimation of D, the values of σint and εint estimated by equation (10) will be greatly overestimated. The REBO potential is not only used to study the elastic behavior of materials but also to simulate the failure (breaking) behavior of materials because it can consider the bond break and reform, which is a significant advantage over other empirical potentials. The failure strain determined by 8
Nanotechnology 26 (2015) 185705
S Xiong and G Cao
Table 3. The cut-off distances selected in the REBO potential.
Cut-off distance (nm)
Mo–Mo (Rmin, Rmax)
S–S (Rmin, Rmax)
Mo–S (Rmin, Rmax)
Aa B Cb D
3.50, 4.05, 4.05, 4.05,
2.30, 2.30, 2.30, 2.30,
2.75 3.05 2.65, 2.70 2.68, 2.73 2.70, 2.75
3.80 4.35 4.35 4.35
3.00 3.00 3.00 3.00
a b
The parameters used in the original REBO potential. The parameters used in the modified REBO potential (REBO*).
underestimated nonlinear behavior (strain-softening effect) in the REBO potential.
4. Conclusions The effectiveness of four empirical potential functions, including CVFF1, CVFF2, SW and REBO potentials, on describing the mechanical behavior of SLMoS2 are evaluated based on the benchmark created by the DFT calculations. To the best of our knowledge, this is the first attempt to employ MD simulations to investigate the tensile response of SLMoS2. Among four empirical potentials, it is found that the CVFF2 and REBO potentials can provide a good description of the elastic behavior of SLMoS2 under small deformation (ε < 5%); CVFF2 can accurately predict the stress value with a very large deformation under biaxial tension; after the cut-off distance modification, the REBO potential gives a good prediction for the failure strain, but the failure stress will be overestimated by about 20–40% under biaxial and uniaxial tension. The REBO potential gives a good prediction of the equilibrium structure of SLMoS2, and the deviations of all geometric parameters are less than 5%, among which there is only 1% difference for both λMo–Mo and λMo–S. The REBO potential can give a good estimation for both the elastic modulus E and the Poisson’s ratio ν of SLMoS2 (about 5% overestimation for both E and ν); it underestimates the strainsoftening effect (D, the nonlinear elastic modulus) of SLMoS2, and thus, it predicts a higher stress value at a given strain ε than the DFT. Compared with the results of graphene, it is concluded that the REBO potential provides a pretty good description of the overall mechanical behavior of SLMoS2 considering the complicated structure of SLMoS2. Both CVFF1 and SW potentials cannot correctly describe the value of E (the former overestimate E by 46% and the latter underestimate E by 40%). However, the SW potential can give a more accurate description of the stress value under the uniaxial or biaxial tension with a large deformation (e.g., ε > 10%) because the errors created in E and D are partially cancelled with each other. In all MD simulations (based on different potentials), the values of σint and εint will be
Figure 7. The effect of the cut-off distances in the REBO potential
on the σ–ε curves of SLMoS2; (a) uniaxial strain tension along the ac direction; (b) uniaxial strain tension along the zz direction; (c) biaxial tension. The parameters of CASEs (A)–(D) are displayed in table 3.
9
Nanotechnology 26 (2015) 185705
S Xiong and G Cao
significantly overestimated by the two-parameter nonlinear elastic model (equation (10)) due to the underestimation of D. Although DFT calculations can provide a more accurate description of the mechanical behavior of SLMoS2, they cannot be used to investigate the structure with a large size or under a complicated loading condition (e.g., free standing indentation). Therefore, the present study will provide a useful help to theoretically investigate the mechanical behavior of SLMoS2 as well as to understand the testing data in experiments.
[14] [15] [16]
[17]
Acknowledgments [18]
I acknowledge the financial support provided by the Ministry of Science and Technology of China (grant no. 2013CB933702) and the National Science Foundation of China (grant no. 11172002).
[19] [20] [21]
References
[22]
[1] Mak K F et al 2010 Atomically thin MoS2: a new direct-gap semiconductor Phys. Rev. Lett. 105 136805 [2] Mak K F et al 2014 The valley hall effect in MoS2 transistors Science 344 1489–92 [3] Radisavljevic B et al 2011 Single-layer MoS2 transistors Nat. Nanotechnology 6 147–50 [4] Deng Y X et al 2014 Black phosphorus-monolayer MoS2 van der waals heterojunction p-n diode ACS Nano 8 8292–9 [5] Liu H et al 2014 Switching mechanism in single-layer molybdenum disulfide transistors: an insight into current flow across schottky barriers ACS Nano 8 1031–8 [6] Mak K F et al 2012 Control of valley polarization in monolayer MoS2 by optical helicity Nat. Nanotechnology 7 494–8 [7] Li Y L et al 2013 Probing symmetry properties of few-layer MoS2 and h-BN by optical second-harmonic generation Nano Lett. 13 3329–33 [8] Bertolazzi S, Krasnozhon D and Kis A 2013 Nonvolatile memory cells based on MoS2/graphene heterostructures ACS Nano 7 3246–52 [9] Chhowalla M and Amaratunga G A J 2000 Thin films of fullerene-like MoS2 nanoparticles with ultra-low friction and wear Nature 407 164–7 [10] Liu K et al 2014 Elastic properties of chemical-vapordeposited monolayer MoS2, WS2, and their bilayer heterostructures Nano Lett. 14 5097–103 [11] Castellanos-Gomez A et al 2012 Elastic properties of freely suspended MoS2 nanosheets Adv. Mater. 24 772 -+ [12] Bertolazzi S, Brivio J and Kis A 2011 Stretching and breaking of ultrathin MoS2 ACS Nano 5 9703–9 [13] Jiang J W, Park H S and Rabczuk T 2013 Molecular dynamics simulations of single-layer molybdenum disulphide (MoS2): stillinger-weber parametrization,
[23] [24] [25] [26] [27] [28] [29] [30] [31] [32]
[33]
10
mechanical properties, and thermal conductivity J. Appl. Phys. 114 064307 Jiang J W, Qi Z N and Park H S 2013 Elastic bending modulus of single-layer molybdenum disulfide (MoS2): finite thickness effect Nanotechnology 24 435705 Peng Q and De S 2013 Outstanding mechanical properties of monolayer MoS2 and its application in elastic energy storage Phys. Chem. Chem. Phys. 15 19427–37 Stewart J A and Spearot D E 2013 Atomistic simulations of nanoindentation on the basal plane of crystalline molybdenum disulfide (MoS2) Modelling Simul. Mater. Sci. Eng. 21 045003 Zhao J H, Jiang J W and Rabczuk T 2013 Temperaturedependent mechanical properties of single-layer molybdenum disulphide: molecular dynamics nanoindentation simulations Appl. Phys. Lett. 103 231913 Cooper R C et al 2013 Nonlinear elastic behavior of twodimensional molybdenum disulfide Phys. Rev. B 87 035423 Jiang J W 2014 The buckling of single-layer MoS2 under uniaxial compression Nanotechnology 25 355402 Li T S 2012 Ideal strength and phonon instability in singlelayer MoS2 Phys. Rev. B 85 235407 Wakabayashi N, Smith H G and Nicklow R M 1975 Latticedynamics of hexagonal Mos2 studied by neutron-scattering Phys. Rev. B 12 659–63 Morita Y et al 2008 Development of a new molecular dynamics method for tribochemical reaction and its application to formation dynamics of MoS(2) tribofilm Appl. Surf. Sci. 254 7618–21 Becker U et al 2003 Metal island growth and dynamics on molybdenite surfaces Geochim. Cosmochim. Acta 67 923–34 Varshney V et al 2010 MD simulations of molybdenum disulphide (MoS2): force-field parameterization and thermal transport behavior Comput. Mater. Sci. 48 101–8 Liang T, Phillpot S R and Sinnott S B 2009 Parametrization of a reactive many-body potential for Mo–S systems Phys. Rev. B 79 245110 Kresse G and Furthmuller J 1996 Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set Phys. Rev. B 54 11169–86 Kresse G and Furthmuller J 1996 Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set Comput. Mater. Sci. 6 15–50 Perdew J P, Burke K and Ernzerhof M 1996 Generalized gradient approximation made simple Phys. Rev. Lett. 77 3865–8 Plimpton S 1995 Fast parallel algorithms for short-range molecular dynamics J. Comput. Phys. 117 1–19 Cao G X 2014 Atomistic studies of mechanical properties of graphene Polymers 6 2404–32 Lee C et al 2008 Measurement of the elastic properties and intrinsic strength of monolayer graphene Science 321 385–8 Lu Q, Gao W and Huang R 2011 Atomistic simulation and continuum modeling of graphene nanoribbons under uniaxial tension Modelling Simul. Mater. Sci. Eng. 19 054006 Liu F, Ming P M and Li J 2007 Ab initio calculation of ideal strength and phonon instability of graphene under tension Phys. Rev. B 76 064120
