Molecular Dynamic Simulations of Pristine and Defective Graphene Nanoribbons Under Strain

RESEARCH ARTICLE Copyright © 2013 American Scientific Publishers All rights reserved Printed in the United States of America Journal of Computational and Theoretical Nanoscience Vol. 10, 470–480, 2013 Molecular Dynamic Simulations of Pristine and Defective Graphene Nanoribbons Under Strain Burcu Tüzün1 and Şakir Erkoç1 2 ∗ 1 Department of Micro and Nano Technology, Middle East Technical University, Ankara 06531, Turkey 2 Department of Physics, Middle East Technical University, Ankara 06531, Turkey Structural properties of pristine and defective graphene nanoribbons have been investigated by stretching them under 5 percent and 10 percent uniaxial strain until the nanoribbons fracture. The stretching process have been carried out by performing molecular dynamics simulations at 1 K and 300 K to determine the temperature effect on the structure of the graphene nanoribbons. Results of the simulations indicated that the conformation of the initial graphene nanoribbon model, temperature, and stretching speed have a considerable effect on the structural properties, however they have a slight effect on the strain value. The maximum strain at which fracture occurs is found to be 46.41 percent for zigzag 8 layer pristine graphene nanoribbon at 1 K and fast stretching process. On the other hand, the defect formation energy is strongly affected from temperature and nanoribbon type. Stone-Wales formation energy is calculated to be 1.60 eV at 1 K whereas 30.13 eV at 300 K for armchair graphene nanoribbon. Keywords: Monovacancy, Divacancy, Stone Wales Defect, Carbon Nanosheet, Graphene, Nanoribbon, Molecular Dynamics Simulations, Atomistic Potential. Delivered by Ingenta to: Rice University IP: 178.57.65.157 On: Tue, 14 Jun 2016 15:30:28 Copyright: American Scientific Publishers 1. INTRODUCTION Carbon is an important element in nanoscale because of its capability of forming different low dimensional structures such as graphene (2D), carbon nanotubes (CNT) (1D), and fullerenes (0D) which have superior physical properties.1 A graphene layer is composed of sp2 hybridized carbon atoms.2 Until 2004, it was thought that films with thickness at nanometer scale are thermodynamically unstable. However, in 2004, Novoselov and coworkers3 were able to prepare single layer graphene. Since the isolation of either single or a few atom thickness graphene layers, it and its applications arouse interest among scientists because of its novel properties.4 2D graphene layer can be considered as a building block of different carbon structures. For example, by rolling the graphene layer and in turn, creating carbon nanotubes (CNT),1 by wrapping graphene and forming spherical buckyballs and by putting them one after the other, graphite can be obtained.4 The popularity of 2D graphene layer comes from its unique electronic properties, and there are lots of studies about the electrical properties of them in literature.5–8 However, mechanical strength of graphene is also marvelous.9 The study carried ∗ Author to whom correspondence should be addressed. 470 J. Comput. Theor. Nanosci. 2013, Vol. 10, No. 2 by Lee et al. showed that the breaking strength of graphene is 200 times larger than that of steel.10 It is obvious that the strength of graphene depends on the number and types of defects and also edge termination. There are different types of defects such as dopinginduced defects caused by substitution of non-carbon atoms in graphene sheet, structural defects, grain boundary defects, bond rotation defect and non-sp2 carbon defects resulted by dangling bonds.1 Apart from mechanical strength, defects also play a crucial role in physicochemical properties, chemical activity, structure, topology, cavity,1 optical, electronic and thermal properties.2 Moreover, conductance and ductility of graphene based devices are controlled by defect characteristic of graphene. In other words, structural defects in graphene layer appear during processing or growth has a great effect on the performance of the graphene based devices. Because of this structural defects, graphene based devices may lose their functional quantity or gain superior or new functionalities depending on type, location or size of the defects. For example, graphene is generally used to enhance and improve the strength of composites. However, research carried out by Bartolucci and coworkers11 showed that the strength of graphene-aluminum composites is lower than that of pure aluminum and aluminum-carbon nanotube composites. The lower mechanical strength is explained by the 1546-1955/2013/10/470/011 doi:10.1166/jctn.2013.2721 Tüzün and Erkoç Molecular Dynamic Simulations of Pristine and Defective Graphene Nanoribbons Under Strain Fig. 1. Zigzag and armchair graphene nanoribbon models. The number of layers (L) define the width of the nanoribbon. J. Comput. Theor. Nanosci. 10, 470–480, 2013 471 RESEARCH ARTICLE defects which become reaction sites with aluminum and result in significant amount of aluminum carbide (Al4 C3 ) as points of weakness.11 In some cases, defected graphene is favorable such as earth metal dispersion on graphene layer and hydrogen adsorption. Hydrogen has very limited interaction with graphene12 which results in using metals such as lithium or palladium to adsorb or decompose hydrogen for sensing.12 13 Metal atoms such as lithium binding on graphene vacancies with energy higher then cohesive and cluster forming energy are the active side for hydrogen adsorption.13 The structure is also affected from defect due to different hybridization ability of carbon atoms which cause different stable structure such as carbyne or graphite. Moreover, the flat structure of graphene layer may be deformed due to non-hexagonal carbon arrangements when polygons can not satisfy symmetry rules. This local curvature allow reactivity and Fig. 2. Types of defects in graphene nanoribbons. adsorption.2 The importance of graphene comes from its large applienergy generation devices due to its electro-catalytic activcation areas. For example, graphene is a very good choice ity, enzymatic binding capacity, unique electrical propfor fabrication of ultrafast transistors because of its abilerties and low production cost. It can also be possible ity of outstanding electron transport at room temperature. to detect the bio-molecules, drugs, heavy metals or catGraphene is also preferred in the biosensors and chemical alytic surfaces by the help of the presence of the defects sensors because of its low capacitance resulting in low sigor dopants. However, knowing information concerning nal to noise ratio, and rapid electron transfer resulting in defects is crucial to control chemical, electronic, mechanaccurate and selective detection of biomolecules.14 Some ical and magnetic properties of graphene.1 robust polymer composites (poly-vinyl-alcohol/graphene), Delivered by Ingenta to:InRice the University present study, 2D graphene nanoribbon (GNR) in high electrical conducting (polystyrene/graphene), high Tue, 14 Jun IP: 178.57.65.157 On: 2016 15:30:28 two different geometries, which are zigzag and armchair thermal conducting polymers (epoxy resins/graphene comCopyright: American Scientific Publishers models, have been studied. The GNR models considered posites) and more thermally stable polymers can be proare shown in Figure 1. duced by using graphene as fillers during fabrication.15 Although the electronic properties of GNR have been Moreover, using graphene instead of CNT decrease the studied extensively, there are limited information in the production cost and difficulties that stem from reliabilliterature about the mechanical properties of GNRs. The ity problem due to chirality of CNT.12 15 Graphene is structural defects considered in the present work, which suitable to be used as electrode material in capacitors16 are monovacancy, divacancy and Stone-Wales defects, are and lithium-ion batteries17 because its wide surface area presented in Figure 2. The other property of defects is (2630 m2 g−1 ), high electrical conductivity (64 mS cm−1 ), mobility. Defects are not stationary, they have a certain stability over vast range of temperature and high elecmobility parallel to the graphene plane. The characteristic tron mobility (2 × 105 cm2 V−1 s−1 at electron densities of migration depends on types of defects and temperature.2 of 2 × 1011 cm−2 for 150 nm graphene nano sheet above In the present simulations the migration of defects have Si/SiO2 gate) make the graphene a good choice for energy storage.4 Apart from these, graphene is preferred also in been also investigated. Molecular Dynamic Simulations of Pristine and Defective Graphene Nanoribbons Under Strain RESEARCH ARTICLE 2. METHOD OF CALCULATIONS Tüzün and Erkoç 3. RESULTS AND DISCUSSION Molecular dynamics (MD) method, known since 1960, is The calculated formation energies Efe of the monovacancy, a very powerful technique to study the materials propdivacancy and Stone-Wales defects for zigzag and armchair erties at nanoscale. Many interparticle interactions are GNRs are given in Tables I and II, respectively. Defect fortaken into account by MD technique which can not be mation energy Efe is defined by the differences between the ignored at nanoscale as done at macroscale.18 In the total energies of the pristine and defective GNR with the present study classical molecular dynamics simulations condition that they have the same width, length and edge (MDS) have been performed in a canonical ensemble form. Namely Efe = ETp − ETd , where ETp is total energy system to investigate the structural properties of GNRs. of pristine GNR, and ETd is the total energy of defective According to canonical ensemble, the temperature is taken GNR. The formation energies give the information about as constant and the number of particles is also conthe tendency of the defect formation on the GNRs and stastant. The Verlet algorithm is used for time integration bility. Defects with low formation energies are more stable and velocity calculations. The interaction between partithan defects with high formation energies.20 DFT calcucles are calculated according to Tersoff potential energy lation results for defect formation energies for CNTs and function parameterized for carbon.19 MDSs have been graphene sheets are available in the literature.21–23 The carried out at 1 and 300 K in order to understand the present results may not be compared with the literature effect of temperature on the defect mobility and fracture values because of the definitions of the defect formation mechanism. energies and also the different method of calculations. Before the stretching process, an appropriate working As seen from Tables I and II the defect formation enercell has been generated by defining the initial position of gies are different for each defect. Defect formation enercarbon atoms. Periodic boundary conditions are applied in gies at 1 K for zigzag GNRs are calculated as 21.47, 26.94 the direction which the stretching process is carried out. and 31.87 eV for monovacancy, Stone-Wales and divaMD relaxation was carried out until the total energy of cancy defects, respectively. From these values, at 1 K for the system becomes constant with respect to time. This zigzag GNR, the defect type can be put in order from is when the system gets to the state of equilibrium. 105 more stable to less stable as monovacancy, Stone-Wales MD steps were enough to reach equilibrium for the models and divacancy. However, Efe at 300 K for zigzag GNR is considered in the present study. After Delivered the systemby reaches different than that of the one at 1 K which affects also Ingenta to: Rice University 178.57.65.157 On: Tue,stability 14 Jun order 2016 15:30:28 equilibrium the stretching processIP: starts with the extension between defects. At 300 K, Efe of divaPublishers of the geometry by 5 percent. At eachCopyright: stretching American process, Scientific cancy (32.58 eV) changes slowly whereas Efe of monothe system was relaxed until equilibrium was reached. vacancy (44.00 eV) and Efe of Stone-Wales (50.59 eV) Stretching and relaxation processes have been continued drastically increase by 104.94 percent and 87.79 percent until some of the bonds are broken and the nanoribbon respectively. The change in percentage is calculated as starts to fracture. Efe300 K − Efe1 K /Efe1 K × 100, where Efe300 K and Efe1 K are In the present study GNRs with two different chirality defect formation energies at 300 K and 1 K, respectively. (zigzag and armchair) and with different layer (the width Therefore, at 300 K, divacancy defect on zigzag GNR is of the nanoribbon) have been studied. To see the effect more stable than monovacancy and Stone-Wales defects of the width on the strength of the zigzag model GNRs which is less stable among them. Defect formation enerwith 4, 5, 6, 7, and 8 layers; and armchair model GNRs gies are not only affected from temperature but they are with 6, 7, 8, 9, 10, 11, and 12 layers were generated. also affected from chirality. At 1 K for armchair GNR, Because of the limited space in the journal papers we have formation energy of Stone-Wales defect is calculated as not included the results of all the models generated. We 1.60 eV which is a much smaller value among all defect just included the results of the largest models considered, formation energies. namely 8 layer model for zigzag GNR and 12 layer model As seen from Table II at 1 K the most stable defect on for armchair GNR. The width of the armchair and zigzag armchair GNR is the Stone-Wales defects which is folGNRs was tried to be kept close to each other by taklowed by monovacancy defect (14.58 eV) and divacancy ing the models at different layers in order to understand defect (23.52 eV). At 300 K the Efe of Stone-Wales defect the chirality effect better. For example, the width of the on arm-chair GNR increase drastically by 1783.2 percent, zigzag model with 8 layers is close to the width of the armchair model with 12 layers as seen from Figure 1. In Table I. Defect formation energies (Efe , in eV) of zigzag 8 layer the present study, 8 layer zigzag GNR with 24.60 Å initial graphene nanoribbons. length and 12 layer armchair GNR with 42.60 Å initial Defect type Efe (at 1 K) Efe (at 300 K) length were compared to investigate chirality effect on the Divacancy 31.87 32.58 strength. In order to eliminate the differences in the initial Monovacancy 21.47 44.00 length, strain values were compared instead of final length Stone-Wales 26.94 50.59 at which fracture occurs. 472 J. Comput. Theor. Nanosci. 10, 470–480, 2013 Tüzün and Erkoç Molecular Dynamic Simulations of Pristine and Defective Graphene Nanoribbons Under Strain Table II. Defect formation energies (Efe , in eV) of armchair 12 layer graphene nanoribbons. J. Comput. Theor. Nanosci. 10, 470–480, 2013 473 RESEARCH ARTICLE Figure 4 shows the results of the simulations for the defective armchair GNR models. These results were obtained at 300 K with slow stretching process. As seen Defect type Efe (at 1 K) Efe (at 300 K) from the figures, the number of runs are the same for monoDivacancy 2352 55.44 vacancy and divacancy defective armchair GNR. Moreover, Monovacancy 1458 62.34 the topology or shape of the GNRs is also similar for monoStone-Wales 160 30.13 vacancy and divacancy defective armchair GNR. This is not expected because the absence of two atoms should have a but it is still the most stable structure among them. At more significant effect than that of the absence of one atom. 300 K, Efe of monovacancy (62.34 eV) and Efe of divaOn the other hand, Stone-Wales defective armchair cancy (55.44 eV) increase by 327.57 percent and 135.71 GNR model is stronger than the other defective modpercent, respectively. At 300 K, monovacancy defect is els. It is an expected result since in Stone-Wales defect more stable than divacancy defect as in the case of zigzag the reason of the defect is not the absence of the atoms GNR at 300 K. Strain values for zigzag and armchair but the non-hexagonal rings surrounded by hexagonal carGNRs with monovacancy defect and pristine are listed in bon rings. Moreover, the non-hexagonal ring changes and Table III. As seen from this table, the maximum strain becomes a hexagonal ring during the stretching process. value is obtained as 46.41 for pristine zigzag GNR at 1 K As a result of stretching a Stone-Wales defective GNR under fast stretching process, while the minimum value is changes to a defect-free GNR, namely Stone-Wales defect obtained as 21.00 for monovacancy defected zigzag GNR disappears during streching process. See Figure 4. at 300 K under fast stretching process. As demonstrated Figure 5 shows the monovacancy defect structure on from Table III strain values are affected from temperature. armchair GNRs under stretch at 1 and 300 K. As seen from For example, for zigzag GNRs strain values obtained at the figure, the motion of the monovacancy defect at 300 K 300 K are smaller than that of at 1 K, which is valid also is more remarkable than that of at 1 K. At 300 K before for armchair model. However, the temperature effect seems the stretching process, defects can not be seen clearly from to be make no difference in the strain values of armchair figure because defect moves toward the edge (a sort of GNRs for fast stretching process. diffusion). During stretching, movement of defect continThe stretching results of the defect free armchair model ues until bond breaking takes place and the nanoribbon Delivered Ingenta to: Rice University at 1 and 300 K can be seen in Figure 3. The by number fracture. The side view of the GNR at 300 K looks more IP: 178.57.65.157 On: Tue,wavy 14 Jun of runs, namely a complete stretching process, are differthan2016 that 15:30:28 of GNR at 1 K as a result of temperature Copyright: American Scientific Publishers ent for slow and fast stretching processes. For example, effect. However, the side view of the GNR at 1 K is almost for 12 layer armchair model at 1 K, the stretching prosmooth. cess continues until 6 runs for slow stretching and 3 runs Figure 6 shows the results of the simulations for armfor fast stretching. It is an expected result since extending chair GNR with divacancy defects under stretch. The interis 5 percent for slow stretching while 10 percent for fast esting and more remarkable point in this figure is that the growing the size of the defect at 1 K is more rapidly and stretching. From Figure 3 the effect of temperature is also at the end of the stretching process, the final size of the seen. At 300 K the topology of the GNRs also changes defect at 1 K is larger than the size of the defect at 300 K. during the stretching process. Figures 7 and 8 display the strain energy, or delta energy (E) versus strain graph of the Stone-Wales defected armTable III. Calculated strain values for various graphene nanoribbons. chair GNRs. E is defined as the difference between final Graphene model and and initial configuration energies, namely E = Ef − Ei . defect type Temperature (K) Stretching type Strain Strain is defined as the percent of the length differences, zz8L-pristine 1 Slow 40.71 namely Strain = 100 × Li − L0 /L0 , where L0 is initial zz8L-monovacancy 1 Slow 34.00 length, whereas Li is the length of the nanoribbon at ith zz8L-pristine 1 Fast 46.41 stretching step or run. zz8L-monovacancy 1 Fast 33.10 The difference in Figures 7 and 8 is very interesting, zz8L-pristine 300 Slow 27.63 because the characteristics of E versus strain show oppozz8L-monovacancy 300 Slow 21.55 zz8L-pristine 300 Fast 33.10 site behaviour at low (1 K) and room temperature (300 K) zz8L-monovacancy 300 Fast 21.00 cases. That is E versus strain changes rapidly for fast ac12L-pristine 1 Slow 34.01 stretching case than the slow stretching case at 1 K. On ac12L-monovacancy 1 Slow 34.01 the other hand, E versus strain changes rapidly for slow ac12L-pristine 1 Fast 33.10 stretching case than the fast stretching case at 300 K. ac12L-monovacancy 1 Fast 33.10 The results of the zigzag GNR stretching simulations are ac12L-pristine 300 Slow 27.62 ac12L-monovacancy 300 Slow 21.55 presented in the Figures 9–13. As seen from the Figure 9, ac12L-pristine 300 Fast 33.10 the pristine zigzag GNR keeps its structure without deforac12L-monovacancy 300 Fast 33.10 mation and fracture sudenly, then we stop stretching at this Tüzün and Erkoç RESEARCH ARTICLE Molecular Dynamic Simulations of Pristine and Defective Graphene Nanoribbons Under Strain Fig. 3. Structures of armchair 12 layer defect free graphene nanoribbons under stretch. Delivered by Ingenta to: Rice University IP: 178.57.65.157 On: Tue, 14 Jun 2016 15:30:28 Copyright: American Scientific Publishers Fig. 4. Structures of defected 12 layer armchair graphene nanoribbons. 474 J. Comput. Theor. Nanosci. 10, 470–480, 2013 Tüzün and Erkoç Molecular Dynamic Simulations of Pristine and Defective Graphene Nanoribbons Under Strain RESEARCH ARTICLE Fig. 5. Structures of 12 layer armchair graphene nanoribbons with monovacancy defect under stretch. point. The same situation does also hold in the armchair edge atoms at 1 K and bond breaking takes place. On the other hand, at 300 K, deformation and bond breaking starts pristine GNR. Delivered by Ingenta to: Rice University the 2016 defect.15:30:28 It is interesting to note that the deforFigure 10 shows the zigzag GNRs with monovacancy. IP: 178.57.65.157 On: Tue,near 14 Jun mation on the edge at 1 K starts on the edge which is As seen from this figure the deformation starts from the Scientific Copyright: American Publishers Fig. 6. Structures of 12 layer armchair graphene nanoribbons with divacancy defect under stretch. J. Comput. Theor. Nanosci. 10, 470–480, 2013 475 Tüzün and Erkoç RESEARCH ARTICLE Molecular Dynamic Simulations of Pristine and Defective Graphene Nanoribbons Under Strain Fig. 7. Structures of armchair graphene nanoribbons with Stone-Wales defect under stretch of at 1 K. The graph shows the strain energy (E) versus strain for the corresponding structures. parallel to the side of the triangular type monovacancy. the fixed edges of the GNR for monovacancy, whereas the On the contrary, the side view of the GNR with monoroughness and/or the wavy character are localized close to vacancy at 300 K is different from the side view of the divacancy position. GNR with divacancy at 300 K. As seen from Figure 10 Figures 12 and 13 display the E versus strain graph Delivered by Ingenta to: Rice University and 11 the wavy structure and roughness are localized at of Stone-Wales defected zigzag GNRs. At 1 K the E IP: 178.57.65.157 On: Tue, 14the Jun 2016 15:30:28 Copyright: American Scientific Publishers Fig. 8. Structures of armchair graphene nanoribbons with Stone-Wales defect under stretch of at 300 K. The graph shows the strain energy (E) versus strain for the corresponding structures. 476 J. Comput. Theor. Nanosci. 10, 470–480, 2013 Tüzün and Erkoç Molecular Dynamic Simulations of Pristine and Defective Graphene Nanoribbons Under Strain RESEARCH ARTICLE Fig. 9. Structures of pristine zigzag graphene nanoribbons under stretch of at 1 K and 300 K. Delivered by Ingenta to: Rice University IP: 178.57.65.157 On: Tue, 14 Jun 2016 15:30:28 Copyright: American Scientific Publishers Fig. 10. Structures of zigzag graphene nanoribbons with monovacancy defect under stretch of at 1 K and 300 K. J. Comput. Theor. Nanosci. 10, 470–480, 2013 477 Tüzün and Erkoç RESEARCH ARTICLE Molecular Dynamic Simulations of Pristine and Defective Graphene Nanoribbons Under Strain Fig. 11. Structures of zigzag graphene nanoribbons with divacancy defect under stretch of at 1 K and 300 K. Fig. 12. Structures of zigzag graphene nanoribbons with Stone-Wales defect under stretch of at 1 K and 300 K. 478 Delivered by Ingenta to: Rice University IP: 178.57.65.157 On: Tue, 14 Jun 2016 15:30:28 Copyright: American Scientific Publishers J. Comput. Theor. Nanosci. 10, 470–480, 2013 Tüzün and Erkoç Molecular Dynamic Simulations of Pristine and Defective Graphene Nanoribbons Under Strain RESEARCH ARTICLE Fig. 13. Structures of zigzag graphene nanoribbons with Stone-Wales defect under stretch of at 300 K. The graph shows the strain energy (E) Delivered by Ingenta to: Rice University versus strain for the corresponding structures. IP: 178.57.65.157 On: Tue, 14 Jun 2016 15:30:28 Copyright: American Scientific Publishers no significant effect on the strain value, whereas it has versus strain graph for both fast and slow stretching look great effect on the strain energy. However, the strain value similar. However, they show different charactheristics at has been affected from both defect and temperature. The 300 K; the E versus strain is larger for slow stretching number of runs taken until the fracture takes place for than that of fast stretching. zigzag GNR is slightly larger than that of the runs for armchair GNR which can be concluded as zigzag GNR is 4. CONCLUSIONS more stable (or stronger) than armchair GNR under strain. The strength and stability of armchair and zigzag GNRs Acknowledgment: The authors would like to thank at two different temperatures have been investigated by Alper Ince for critically reading the manuscript. MDS. The result of the simulations were used to calculate strain values and defect formation energies. From the results of the simulations we can conclude that, at References nanoscale, defect formation energies are affected strongly 1. M. Terrones, A. R. Botello-Méndez, J. Campos-Delgado, F. Lópezfrom temperature and chirality. At room temperature the Urías, Y. I. Vega-Cantú, F. J. Rodríguez-Macías, A. L. Elías, mobility of defects can be seen clearly. Moreover, the E. Muñoz-Sandoval, A. G. Cano-Márquez, J.-C. Charlier, and topology and the final size of the defects are also affected H. Terrones, Nano Today 5, 351 (2010). 2. F. Banhart, J. Kotakoski, and A. V. Krasheninnikov, ACS Nano 5, 26 from temperature. The differences in the strength and in (2011). the stability of zigzag and armchair GNR can be seen from 3. K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, the differences in defects formation energies and the numS. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Science 306, 666 ber of runs in the stretching process which takes place until (2004). the fracture occurs. The maximum strain value was cal4. D. A. C. Brownson, D. K. Kampouris, and C. E. Banks, J. Power Sources 196, 4873 (2011). culated as 46.41 percent for zigzag pristine model during 5. Z. Chen, Y.-M. Lin, M. J. Rooks, and P. Avouris, Physica E 40, 228 fast stretching at 1 K, while the minimum strain value was (2007). calculated as 21.00 percent for zigzag monovacancy GNR 6. J. Bai and Y. Huang, Mat. Sci. Eng. R 70, 341 (2010). during fast stretching at 300 K. From these values, it can 7. S. Banerjee and D. Bhattacharyya, Comp. Mater. Sci. 44, 41 (2008). be concluded that the fast and slow stretching process has 8. N. M. R. Peres, Vacuum 83, 1248 (2009). J. Comput. Theor. Nanosci. 10, 470–480, 2013 479 RESEARCH ARTICLE Molecular Dynamic Simulations of Pristine and Defective Graphene Nanoribbons Under Strain 9. H. Bu, Y. Chen, M. Zou, H. Yi, K. Bi, and Z. Ni, Phys. Lett. A 373, 3359 (2009). 10. C. Lee, X. Wei, J. W. Kysar, and J. Hone, Science 321, 385 (2008). 11. S. F. Bartolucci, J. Paras, M. A. Rafiee, J. Rafiee, S. Lee, D. Kapoor, and N. Koratkar, Mat. Sci. Eng. A-Struct. 528, 7933 (2011). 12. W. Wu, Z. Liu, L. A. Jauregui, Q. Yu, R. Pillai, H. Cao, J. Bao, Y. P. Chen, and S.-S. Pei, Sensor Actuat. B-Chem. 150, 296 (2010). 13. E. Rangel, J. R. Ramŕres-Arellano, I. Carrillo, and L. F. Magana, Int. J. Hydrogen Energ. 36, 13657 (2011). 14. Y. Shao, J. Wang, H. Wu, J. Liu, I. A. Aksay, and Y. Lin, Electroanal 22, 1027 (2010). 15. T. Kuilla, S. Bhadra, D. Yao, N. H. Kim, S. Bose, and J. H. Lee, Prog. Polym. Sci. 35, 1350 (2010). Tüzün and Erkoç 16. X. Du, P. Guo, H. Song, and X. Chen, Electrochim. Acta 55, 4812 (2010). 17. G. Wang, X. Shen, J. Yao, and J. Park, Carbon 47, 2049 (2009). 18. D. C. Rapaport, The Art of Molecular Dynamics Simulation, Cambridge University Press, Cambridge (1995). 19. J. Tersoff, Phys. Rev. Lett. 61, 2879 (1988). 20. J. E. Padilha, R. G. Amorim, A. R. Rocha, A. J. R. daSilva, and A. Fazzio, Solid State Commun. 151, 482 (2011). 21. J. Ma, D. Alfè, A. Michaelides, and E. Wang, Phys. Rev. B 80, 033407 (2009). 22. G. E. G. Samsonidze, G. G. Samsonidze, and B. I. Yakobson, Comp. Mater. Sci. 23, 62 (2002). 23. A. V. Krasheninnikov, P. O. Lethinen, A. S. Foster, and R. M. Nieminen, Chem. Phys. Lett. 418, 132 (2006). Received: 19 December 2011. Accepted: 5 January 2012. Delivered by Ingenta to: Rice University IP: 178.57.65.157 On: Tue, 14 Jun 2016 15:30:28 Copyright: American Scientific Publishers 480 J. Comput. Theor. Nanosci. 10, 470–480, 2013