Modeling of the interaction between polypropylene and monolayer sheets: a quantum mechanical study

Abstract

In this paper, we performed quantum mechanical calculations to determine the best monolayer sheet for preparing polypropylene nanocomposites. For this purpose, six types of monolayer sheets were selected (graphene, silicon-carbide monolayers (SiC-m), boron-nitride monolayers (BN-m), aluminum-nitride monolayers (AlN-m), gallium-nitride monolayers (GaN-m) and silicene). First, the mechanical properties of the different monolayers were studied; the results indicated that the Young's modulus of graphene is six times greater than that of silicene. Next, the interaction energies between the propylene monomer and monolayer sheets were calculated. The results indicated that among the different monolayer sheets, the Young's modulus of graphene is higher than that of the other materials but that silicene has an eightfold stronger interaction with the propylene monomer. Then, we calculated the interaction energy and equilibrium distance between polypropylene and silicene in nanocomposites with increasing propylene monomer contents on the silicene. Finally, Lennard–Jones parameters (σ and ε) were calculated to model the interaction between polypropylene and silicene with a linear spring.

---

1. Introduction

Polymer nanocomposites, used as nanoscale fillers at low loading, have been widely investigated in academia and industry because of their exceptional multifunctional properties compared to customary polymer composites.^1–3 However, in the composites, maximal mechanical improvement is obtained when the exterior load is efficiently transferred from the polymer to the filler via a strong interaction at the interface between the matrix and nanofiller.^4,5 Therefore, the choice of a nanofiller with excellent mechanical properties and strong adhesion with the matrix is an important issue in the processing of composites.

The study of two-dimensional (2D) nanomaterials with potential for use in next-generation electronic devices and the application of nanocomposite materials has seen marvelous progress in the past few years.^6,7 An example of such nanostructures is graphene, which is composed of carbons arrayed in a 2D honeycomb monolayer fashion.^4,8 In 2004, Novoselov et al.⁹ found a way to experimentally separate individual graphene nanoribbon planes. Among 2D nanomaterials, over 20 hexagonal monolayer structures have been theoretically predicted to be stable analogs of graphene, each of which exhibits a honeycomb lattice structure.¹⁰

Other group IV elements, such as silicon (Si), reveal chemistry similar to that of carbon in certain aspects and have structures with stable honeycomb monolayers (namely silicene). In recent years, silicene has been the subject of many experimental and theoretical investigations. Silicene was first predicted to exist based on ab initio calculations in 1994 and has been recently synthesized over silver (111) and (110) substrates.^11–13 Silicene is a zero-gap semimetallic material, and the electronic structure of silicene is quite similar to that of graphene.^14,15 Unlike graphene, silicene is not stable as a perfectly planar sheet and has chair-like distortions in the rings that lead to ordered surface ripples and enhance the reactivity with the surfaces of other materials. In silicene, sp³ hybridization is energetically more favorable, whereas sp² hybridization is more stable in graphene due to the smaller ionic radius of carbon. Although it is difficult to construct stable Silicon structures in silicene, it is possible to form a single-layered structure by mixing C and Si. Unlike graphene, silicon-carbide monolayers (SiC-m) are polar materials and, therefore, may demonstrate unusual physical properties that graphene may not display. For example, it is well known that graphene is a zero-band-gap semimetal, whereas SiC-m exhibit finite band-gap semiconductor properties.¹⁶

Substituting the C in graphene with group III elements and N atoms has been shown to be a promising way to improve the semiconductor properties of graphene.^17,18 For example, honeycomb boron-nitride monolayer (BN-m) and graphene have similar 2D lattice structures but very different physical properties. BN-m is insulator with a direct band gap of 5.97 eV.¹⁹ Due to the strong covalent sp² bonds within the plane of BN-m, the in-plane mechanical strength and thermal conductivity of BN-m have been reported to be close to those of graphene. BN-m have an even greater chemical stability than graphene; BN-m can be steady in air up to 1000 °C (in contrast, for graphene, the corresponding temperature is 600 °C).²⁰ Due to these notable properties, BN-m have found widespread applications in micro- and nano-devices, for example, as insulators with high thermal conductivity in electronic devices, as ultraviolet-light emitters in optoelectronics, and as nano-fillers in high-strength and thermal-conductive nanocomposites.

The modeling and simulation approaches can be divided into two main categories: (I) atomistic simulations at the nanoscale; and (II) continuum and structural mechanics modeling.²¹ The atomic simulations are mainly based on classical molecular dynamics (MD) with empirical interatomic potentials, the semi-empirical methods (tight-binding) and the ab initio quantum mechanics (QM) calculations. However, due to the challenges involved in devising experiments to investigate the interfacial interaction mechanisms involved in polymer-nanofiller composite systems, MD and QM simulations have become increasingly popular. The QM calculations based on density functional theory (DFT) can clearly calculate the interactions involving electrons to provide the most accurate description of the geometry, interaction mechanism and mechanical behavior of systems.²²

To the best of our knowledge, there have been no published investigations regarding the properties of the interactions (interaction energy and Lennard–Jones (L–J) parameters) of carbon- and non-carbon-based monolayer sheets to which propylene monomers have been adsorbed. Therefore, in the present study, we used the DFT method to investigate the interaction energy of propylene monomers with carbon- and non-carbon-based monolayer sheets (graphene, SiC-m, BN-m, aluminum-nitride monolayers (AlN-m), gallium-nitride monolayers (GaN-m) and silicene) to determine a suitable type of monolayer sheet for constructing polypropylene nanocomposites. Furthermore, the mechanical properties of the aforementioned monolayer sheets were investigated to determine the better choice for a nanofiller. We also evaluated the effect of the monomer content on the interaction energy of polypropylene nanocomposites with the best expected monolayer filler (silicene). Finally, the L–J parameters were determined according to the L–J 6–12 potential function between polypropylene and silicene, and the interaction was then modeled as a classical spring.

2. Computational method

The atomic geometry, interaction energy and mechanical properties of carbon- and non-carbon-based monolayer sheets in the presence and absence of propylene monomers were calculated within the framework of the ab initio DFT method^23,24 and performed by the Spanish Initiative for Electronic Simulations with Thousands of Atoms (SIESTA) code.^25,26 A significant limitation of the LDA is its overbinding of extended solids. Therefore, the lattice parameters are typically underpredicted, while the cohesive energies, phonon frequencies, and elastic moduli are typically overpredicted. Due to this limitation of LDA, we used the generalized-gradient approximation (GGA) function to treat the electronic exchange and correlation effects, as described by Perdew–Burke–Ernzerhof (PBE). In all of the procedures, a split-valence double-ζ basis set of localized numerical atomic orbitals was used, including the polarization functions (DZP) with an energy shift of 50 meV and a split norm of 0.25.^26,27 This method offers a reasonable balance between accuracy and computational resources. A 5 × 5 × 1 Monkhorst–Pack grid for the k-point sampling of the Brillouin zone was set, and the atomic positions were relaxed until the residual forces on each atom were lower than 0.03 eV Å⁻¹. An energy cutoff of 120 Ry for the grid integration was chosen to depict the charge density.²⁶

3. Results and discussion

3.1. Mechanical properties of monolayer sheets

Due to the role of adhesion in the exterior load transfer from the matrix to the nanofiller, we first calculate the mechanical properties of different carbon- and non-carbon-based monolayer sheets (graphene, SiC-m, BN-m, AlN-m, GaN-m and silicene) and then consider the interaction between propylene and these monolayer sheets.

In this section, DFT calculations were first performed to obtain the completely optimized geometry for 50-atom monolayer sheets of graphene, SiC-m, BN-m, AlN-m, GaN-m and silicene. After full structural optimization, we found that the bond lengths of graphene, SiC-m, BN-m, AlN-m, GaN-m and silicene monolayer sheets are approximately 1.42, 1.81, 1.45, 1.83, 1.86, and 2.27 Å, respectively; these results are in good agreement with previous theoretical studies.^{10,18,28–31} The bond lengths of these monolayers after optimization are listed in Table 1. As shown in this table, the Ga–N bond length in the GaN-m monolayer sheets is approximately 1.86 Å that is in good agreement with the result obtained by Peng et al. (1.85 Å).³⁰

Table 1 Bond lengths of monolayer sheets after optimization, compared with the literature

Sheets	Element 1	Element 2	Bond length (Å) (this work)	Bond length (Å) (in literatures)	Ref.
Graphene	C	C	1.42	1.41	10
SiC-m	Si	C	1.81	1.79	28
BN-m	B	N	1.45	1.44	29
AlN-m	Al	N	1.83	1.81	18
GaN-m	Ga	N	1.86	1.85	30
Silicene	Si	Si	2.27	2.29	31

---

From a continuum perspective, the elastic properties of a material are determined from the elastic strain energy density, Φ, which is defined as below.¹⁴

E_s = E_tot − E₀ (2)

where E_s is the strain energy, E_tot is the total energy of the strained system, E₀ is the total energy of the strain-free system and V is the volume of the undeformed supercell. Nonlinear elastic constitutive behavior is established by expanding Φ in a Taylor series in terms of powers of strain η, as follows:³²where η_ij is the Lagrangian elastic strain. The summation convention is employed for repeating indices; the lower-case subscripts range from 1 to 3. Herein, C denotes each higher-order elastic modulus tensor; the rank of each tensor corresponds to the number of subscripts. The second-order elastic constants, C_ijkl, third-order elastic constants, C_ijklmn, and fourth-order elastic constants, C_ijklmnop, are given by the components of the fourth-, sixth-, and eighth-rank tensors, respectively. We employed the conventional Voigt notation for subscripts: 11 → 1, 22 → 2, 33 → 3, 23 → 4, 31 → 5, and 12 → 6.³³ Furthermore, note that for strain, we have η₄ = 2η₂₃, η₅ = 2η₃₁, and η₆ = 2η₁₂. The summation convention for upper-case subscripts runs from 1 to 6. Eqn (3) can be rewritten as follows:

The symmetric second Piola–Kirchhoff stress tensor, π_ij, is

Using the Voigt notation for subscripts, eqn (5) can be rewritten as

In this study, we are interested only in the in-plane elasticity. The second-order elastic constants model the linear elastic response. Thus, the elasticity can be expressed in a matrix form as eqn (7).

The strains of ±0.03 are applied along two directions, x (uniaxial) and xy (biaxial). The value of the elastic constants can also be calculated based on the uniaxial/biaxial strain:³²

where A is the area of the undeformed supercell. Variations in the strain energies of graphene with respect to the uniaxial/biaxial strain are presented in Fig. 1. Moreover, the results for the elastic constants for the monolayer sheets are listed in Table 2. Finally, the in-plane Young's modulus (or stiffness), Y_s, and Poisson's ratio, ν, may be obtained from the following relationships:³²

Fig. 1 Strain energy versus strain for uniaxial and biaxial strain in graphene.

Table 2 Elastic constant, in-plane stiffness, and Poisson's ratio of the monolayer sheets, compared with the literatures

Sheets	C11 (N m^−1)	C12 (N m^−1)	ν	Ys (N m^−1) (this work)	Ys (N m^−1) (in literatures)	Ref.
Graphene	354.1	57.0	0.161	345.1	347	17
SiC-m	172.3	47.7	0.277	159.2	166	34
BN-m	281.0	51.4	0.183	273.8	278	30
AlN-m	160.3	63.0	0.393	135.5	136	18
GaN-m	126.5	58.1	0.460	99.7	109	30
Silicene	63.6	18.0	0.283	58.2	62	34

---

The values of Y_s for graphene, SiC-m, BN-m, AlN-m, GaN-m and silicene monolayer sheets calculated according to the present method are equal to 345.1, 159.2, 273.8, 135.5, 99.7, and 58.2 N m⁻¹, respectively. In addition, the calculated Poisson's ratios for graphene, SiC-m, BN-m, AlN-m, GaN-m and silicene monolayer sheets are 0.161, 0.276, 0.183, 0.393, 0.460, and 0.283, respectively. Our results of Y_s and ν for the six monolayer sheets are reported in Table 2. For all of the monolayer sheets, our calculated elastic parameters are in good agreement with the experimental values and previous calculations.^30,33,34 For example, our calculated value of Y_s for graphene (345.1 N m⁻¹) is in good agreement with the experimental value (340 ± 50 N m⁻¹)³³ and the DFT theoretical predictions (347.2 N m⁻¹).¹⁷ Based on these results, the order of the stiffness of the monolayer sheets is as follows: graphene > BN-m > SiC-m > AlN-m > GaN-m > silicene. Furthermore, we found that the in-plane Young's modulus of graphene (345.1 N m⁻¹) is six times greater than that of silicene.

3.2. Interaction between monomer and sheets

Now, we consider the interaction between propylene monomer (–CH₃–CH–CH₂–) and different types of monolayer sheets. For this purpose, we first approach the propylene monomer to the surface of the sheets at different possible sites to find the most favorable configuration of propylene on the monolayer sheets. The interaction energy between the propylene monomer and the monolayer sheet was calculated as follows:³⁵

E_int = (E_pp/g − E_g − E_pp)/n (12)

where E_pp/g is the total energy of the combined system of the propylene monomer and monolayer sheet, E_g is the total energy of the bare sheet, E_pp is the total energy of the free propylene monomer, and n is the number of monomers. Therefore, a negative binding energy corresponds to a stable adsorption structure. After a full structural optimization for all of the systems, we found that the interaction energy between propylene and graphene (propylene/graphene) is approximately −0.11 eV in the favorable configuration. The equilibrium distance between the closest atom of the propylene and monolayer sheet (r_eq) is measured at approximately 2.94 Å. The optimized configurations and the geometrical parameters of the propylene monomer on the surface of graphene in the favorable configuration are depicted in Fig. 2. From Fig. 2(a), we also found that the propylene monomer is attached with the small incline (approximately 1°) on the surface of graphene. Therefore, we can assume that propylene monomer is parallel to the graphene surface. To evaluate the interaction energy between propylene monomers and group III (B, Al, and Ga)-nitride monolayer sheets, we placed the monomer at different positions and orientations above the sheets to find their most favorable configurations. We first computed the interactions between the propylene and the BN-m sheets (propylene/BN-m). We found that interaction energy and the equilibrium distance in this system are equal to −0.20 eV and 2.76 Å, respectively. Fully optimized geometrical structures and the equilibrium distance for this system in the most favorable configuration are illustrated in Fig. 3. Comparison of the calculated binding energy and equilibrium distance for the propylene monomer on the surface of graphene and BN-m sheet reveals that the interaction energy for the propylene on the BN-m sheet is higher than that on graphene, and the equilibrium distance between the monomer and the BN-m sheet is shorter than the corresponding value determined for graphene. The above results reveal that the interaction of propylene with the BN-m sheet is stronger than the interaction with graphene. Finally, we evaluated the interaction energies between propylene and AlN-m and GaN-m sheets (propylene/AlN-m and GaN-m) and found that these values are equal to −0.04 and −0.02 eV, respectively. Therefore, from the obtained results, we can conclude that the interaction energies between the propylene and the monolayer sheets are in the following order: GaN-m < AlN-m < graphene < BN-m. These results indicate that the adhesion of the propylene monomer is dependent on the type of monolayer sheet.

Fig. 2 (a) Side-view and (b) top-view of the optimized structure and geometrical parameters of a propylene monomer on graphene in the most stable configuration.

Fig. 3 Optimized structures and geometrical parameters of a propylene monomer on BN-m in the most stable configuration.

We also calculate the interaction energy between propylene and the other monolayer sheets, namely SiC-m and silicene. The interaction of propylene with SiC-m (propylene/SiC-m) is found to be equal to 0.11 eV, which is lower than that of graphene. Finally, using silicene, the interaction energy between propylene and the monolayer sheet (propylene/silicene) was calculated to be 0.8 eV, which represents an increase of approximately eightfold more than that of the propylene/graphene system. Fig. 4 is a schematic representation of the optimized geometric structure for the propylene/silicene system. As seen in this figure, the equilibrium distance in this system, in agreement with the calculated interaction energy, is equal to 2.35 Å and lower than the equilibrium distance in the propylene/graphene system. According to the above results, we calculate that the interaction energy between the propylene monomer and the silicene sheet is greater than the propylene monomer interaction energies with the other discussed sheets. In addition, it is obvious that the adsorption of the propylene monomer onto the silicene is eight times as strong as the adsorption of the same monomer onto the graphene. Therefore, we propose that using silicene as a nanofiller in a polypropylene matrix can lead to an increase in the mechanical properties of the nanocomposite due to strong interfacial interactions.

Fig. 4 Optimized structures and geometrical parameters of a propylene monomer on silicene in the most stable configuration.

In this section, we study the dependence of the interaction energy on the number of propylene monomers. For this purpose, isotactic chains of propylene monomers [(–CH3–CH–CH2–)_n] with n = 1–7 were placed on the silicene surface with different possible configurations to find the most energetically favorable configuration and the interaction energy. The optimized structures and geometrical parameters of the systems with three, five and seven monomers of propylene on a silicene sheets are presented in Fig. 5. According to our results, with an increasing number of monomers on the surface of silicene, the interaction energy initially decreases rapidly, then slowly decreases, and finally approaches a constant value of −0.19 eV/monomer when seven propylene monomers are adsorbed onto the surface of silicene. The results obtained for the interaction energy between propylene and silicene as a function of the amount of monomer are presented in Fig. 6. Furthermore, as shown in Fig. 5, in agreement with the interaction energy results, we find that the equilibrium distances increase from 2.35 Å for propylene/silicene to 2.85 Å for a seven monomer propylene/silicene interaction. Thus, in modeling the polypropylene (PP)/silicene nanocomposites, the equilibrium distance between the PP polymer and silicene can be defined as equal to 2.85 Å.

Fig. 5 Optimized structures and geometrical parameters of (a) three, (b) five and (c) seven monomers of propylene on a silicene sheet.

Fig. 6 Variation in the interaction energy with the number of monomers adsorbed onto the silicene surface.

Because the monolayer sheets are modeled at the atomistic scale and the polymeric matrix is treated as a continuum, the modeling of the monolayer sheet/polymer interface is rather difficult. According to the above results, van der Waals interfacial interactions exist between the silicene and the PP matrix, which can be modeled using the general form of the L–J potential as below:³⁶

where U(r) is the potential energy between a pair of molecules, r is the distance between the pair of molecules, ε is the strength of interaction between molecules, and σ is the van der Waals separation distance. According to the L–J potential equation, eqn (13), the equilibrium distance between two molecules (r_eq) can occur at the minimum potential energy. Thus, the first derivative is zero. The expression for dU(r)/dr is obtained by differentiating the expression for the L–J potential, and then r_eq is determined by setting the expression to zero

Thus, we can calculate the value of L–J parameter σ according to:

σ = 2^−1/6r_eq (15)

Furthermore, the value of L–J parameter ε can be calculated by obtaining the potential energy at the minimum, as below:

U(r_eq) = −ε (16)

To determine the L–J parameters for the interaction between the silicene and the PP polymer matrix, the distance between the sheet and the molecule is varied, and the potential energy between the propylene monomers and the silicene as a function of the separation is plotted (as shown in Fig. 7). After curve fitting, the values of σ and ε are found to be 2.54 Å and 33 kJ mol⁻¹, respectively. Recently, Androulidakis et al.³⁷ have used MP2 calculations to predict the values of σ and ε for graphene and a methyl methacrylate monomer. They found values of σ and ε equal to 2.81 Å and 18.98 kJ mol⁻¹, respectively.

Fig. 7 Potential energy of a PP/silicene system as a function of the separation distance.

We are interested in comparing the effect of the L–J interaction with the classical spring model. We can Taylor-expand the potential eqn (13) around r_eq to find the effective value of the spring constant. For a sufficiently small displacement, we can neglect the terms beyond the third term in the power series. Furthermore, because U(r) is a minimum at r_eq, . To obtain an expression for the spring constant in terms of the potential parameters σ and ε, we take the second derivative:

This result is a potential energy function for a harmonic oscillator with a spring constant k. Consequently, the spring constant replaced by the van der Waals interaction between the PP polymer matrix and the silicene monolayer sheet is equal to 15.72 N m⁻¹ (94.78 kJ mol⁻¹ Ų). In addition, we calculate the spring constant per unit area and find that this value is equal to 12.5 Pa m⁻¹ (0.75 kJ mol⁻¹ Å⁴).

4. Conclusions

In this work, we first calculated the mechanical properties of different carbon- and non-carbon-based monolayer sheets (graphene, SiC-m, BN-m, AlN-m, GaN-m and silicene) and then considered the interaction between a propylene monomer and these monolayer sheets using the first-principles DFT approach to find the best nanofiller for reinforcing the polypropylene. We found that silicene has a weaker in-plane Young's modulus than the other monolayers (Y_s = 58.2 N m⁻¹) and that graphene has the strongest modulus (Y_s = 345.1 N m⁻¹). Furthermore, the strongest interaction between the propylene monomer and monolayer sheets occurred with silicene, with the interaction energy of 0.8 eV, which is eightfold stronger than that of graphene. We also analyzed the effect of the quantity of propylene monomer on the intermolecular binding between polypropylene and silicene, which was the best proposed nanofiller. The interaction energy between propylene monomers and silicene was found to decrease as the number of monomers increases and finally reach a constant value. Moreover, the L–J parameters (σ and ε) were calculated to model the interaction between PP and silicene using a linear spring. We found that σ, ε, and the spring constant per unit area are equal to 0.254 nm, 33 kJ mol⁻¹, and 12.5 Pa m⁻¹, respectively.

References

1. M. Razavi-Nouri, M. Ghorbanzadeh-Ahangari, A. Fereidoon and M. Jahanshahi, Polym. Test., 2009, 28, 46–52.
2. K. Hu, D. D. Kulkarni, I. Choi and V. V. Tsukruk, Prog. Polym. Sci., 2014, 39, 1934–1972.
3. N. Bitinis, M. Hernández, R. Verdejo, J. M. Kenny and M. A. Lopez-Manchado, Adv. Mater., 2011, 23, 5229–5236.
4. B. Shen, W. Zhai, M. Tao, D. Lu and W. Zheng, Compos. Sci. Technol., 2013, 77, 87–94.
5. A. Fereidoon, M. G. Ahangari and S. Saedodin, J. Macromol. Sci., Part B: Phys., 2009, 48, 196–211.
6. G. Lalwani, A. M. Henslee, B. Farshid, L. Lin, F. K. Kasper, Y. X. Qin, A. G. Mikos and B. Sitharaman, Biomacromolecules, 2013, 14, 900–909.
7. S. Z. Butler, S. M. Hollen, L. Cao, Y. Cui, J. A. Gupta, H. R. Gutierrez, T. F. Heinz and S. S. Hong, et al., ACS Nano, 2013, 7, 2898–2926.
8. L. Yang, W. A. Yee, S. L. Phua, J. Kong, H. Ding, J. W. Cheah and X. Lu, RSC Adv., 2012, 2, 2208–2210.
9. K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva and A. A. Firsov, Science, 2004, 306, 666–669.
10. M. Q. Le, Int. J. Mech. Mater. Des., 2014, 11, 15–24.
11. K. Takeda and K. Shiraishi, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 14916–14922.
12. C. Leandri, G. L. Lay, B. Aufray, C. Girardeaux, J. Avila, M. E. Dávila, M. C. Asensio, C. Ottaviani and A. Cricenti, Surf. Sci., 2005, 574, L9–L15.
13. B. Lalmi, H. Oughaddou, H. Enriquez, A. Kara, S. Vizzini, B. Ealet and B. Aufray, Appl. Phys. Lett., 2010, 97, 223109.
14. P. Xu, Z. Yu, C. Yang, P. Lu, Y. Liu, H. Ye and T. Gao, Superlattices Microstruct., 2014, 75, 647–656.
15. T. Botari, E. Perim, P. A. S. Autreto, A. C. van Duin, R. Paupitz and D. S. Galvao, Phys. Chem. Chem. Phys., 2014, 16, 19417–19423.
16. E. Bekaroglu, M. Topsakal, S. Cahangirov and S. Ciraci, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 075433.
17. Q. Peng, A. R. Zamiri, W. Ji and S. De, Acta Mech. Sin., 2012, 223, 2591–2596.
18. Q. Peng, X. J. Chen, S. Liu and S. De, RSC Adv., 2013, 3, 7083–7092.
19. Y. Kubota, K. Watanabe, O. Tsuda and T. Taniguchi, Science, 2007, 317, 932–934.
20. K. K. Kim, A. Hsu, X. Jia, S. M. Kim, Y. Shi, M. Hofmann, D. Nezich, J. F. Rodriguez-Nieva, M. Dresselhaus, T. Palacios and J. Kong, Nano Lett., 2011, 12, 161–166.
21. M. D. Ganji, A. Fereidoon, M. Jahanshahi and M. G. Ahangari, Solid State Commun., 2012, 152, 1526–1530.
22. M. G. Ahangari, A. Fereidoon, M. Jahanshahi and M. D. Ganji, Phys. E, 2013, 48, 148–156.
23. P. Hohenberg and W. Kohn, Phys. Rev., 1964, 136, B864–B871.
24. W. Kohn, Rev. Mod. Phys., 1999, 71, 1253–1266.
25. P. Ordejon, E. Artacho and J. M. Soler, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 53, R10441–R10444.
26. J. M. Soler, E. Artacho, J. D. Gale, A. Garcia, J. Junquera, P. Ordejon and D. Sanchez-Portal, J. Phys.: Condens. Matter, 2002, 14, 2745–2779.
27. J. Junquera, O. Paz, D. Sanchez-Portal and E. Artacho, Phys. Rev. B: Condens. Matter Mater. Phys., 2001, 64, 235111.
28. E. Bekaroglu, M. Topsakal, S. Cahangirov and S. Ciraci, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 075433.
29. X. Zhong, Y. K. Yap and R. Pandey, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 83, 193403.
30. Q. Peng, C. Liang, W. Ji and S. De, Appl. Phys. A, 2013, 113, 483–490.
31. T. H. Osborn, A. A. Farajian, O. V. Pupysheva, R. S. Aga and L. L. Y. Voon, Chem. Phys. Lett., 2011, 511, 101–105.
32. X. Wei, B. Fragneaud, C. A. Marianetti and J. W. Kysar, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 80, 205407.
33. C. Lee, X. Wei, J. W. Kysar and J. Hone, Science, 2008, 321, 385–388.
34. M. Topsakal, S. Cahangirov and S. Ciraci, Appl. Phys. Lett., 2010, 96, 091912.
35. M. G. Ahangari, A. Fereidoon and M. D. Ganji, J. Mol. Model., 2013, 19, 3127–3134.
36. R. Retkute and J. P. Gleeson, Dynamics of two coupled particles: comparison of Lennard–Jones and spring forces, in SPIE Third International Symposium on Fluctuations and Noise, International Society for Optics and Photonics, 2005, pp. 61–71.
37. C. Androulidakis, E. N. Koukaras, O. Frank, G. Tsoukleri, D. Sfyris, J. Parthenios, N. Pugno, K. Papagelis, K. S. Novoselov and C. Galiotis, Sci. Rep., 2014, 4, 5271.

---