A Computational Study On the Mechanical Properties of Pentahexoctite Single-layer: Combining DFT and Classical Molecular Dynamics Simulations

Studies aimed at designing new allotropic forms of carbon have received much attention. Recently, a new 2D graphene-like allotrope named Pentahexoctite was theoretically proposed. Pentahexoctite has a metallic signature, and its structure consists of continuous 5-6-8 rings of carbon atoms with sp2 hybridization. Here, we carried out fully-atomistic computational simulations, combining reactive (ReaxFF) molecular dynamics (MD) and density functional theory (DFT) methods, to study the elastic properties and fracture patterns of Pentahexoctite monolayer. Results revealed a Young's Modulus of 0.74 TPa, smaller than the graphene one (about 1.0 TPa). The Pentahexoctite monolayer, when subjected to a critical strain, goes directly from elastic to completely fractured regimes. This process occurs with no plasticity stages between these two regimes. Importantly, graphene presents a similar fracture process. The elastic properties calculated with both DFT and MD are in good agreement. Keywords: Reactive Molecular Dynamics, DFT, Mechanical Properties, Carbon Allotrope,


Introduction
The accomplishments in producing the single-layer graphene (by mechanical exfoliation) resulted in the Nobel Prize in physics earned by Novoselov and Geim in 2004 [1]. Graphene is a one-carbon-atom-thick layer of graphite exhibiting unique properties, such as light-weight, good transparency, and high mechanical strength, and its advent gave rise to a new era in materials science [2]. As a consequence of its success, several experimental [3,4] and theoretical [5][6][7][8][9][10] studies have been carried out to propose new graphene-like allotropes. One objective of these studies is to propose new structures that can overcome some graphene limitations, such as high reactive to oxygen and null bandgap [11].
Based on that, a new 2D carbon allotrope graphene-like, named Pentahexoctite (PH), was recently proposed as a new metallic nanostructure composed of continuous 5-6-8 rings of carbon atoms with sp 2 hybridizations [12]. It was found that the Pentahexoctite presents a planar sheet with the absence of negative frequencies in the phonon spectra, which confirms its structural stability. They also reported that PH has mechanical strength comparable to graphene. In relation to the PH electronic properties, it exhibits metallic behavior with direction-dependent flat and dispersive bands at the Fermi level ensuring highly anisotropic transport properties. Its stability and electronic properties were well described in references [12,13]. However, its mechanical properties, elastic constants, and fracture patterns have not been fully investigated.
In this work, we carried out extensive fully-atomistic computational simulations using the reactive (ReaxFF) molecular dynamics (MD) and density functional theory methods to investigate the PH mechanical properties and fracture patterns. Our results show Young's Modulus of 247.93 GPa.nm (or 0.74 TPa). This value is comparable to the graphene one. The PH monolayer, when subjected to a critical strain, goes directly from elastic to completely fractured regimes. This trend is also exhibited by graphene. We also studied the impact of temperature on the PH fracture process. The results suggest that Young's Modulus varies between 0.82 to 0.50 TPa for temperatures ranging from 10 up to 1200 K, respectively. Importantly, the elastic properties calculated with both DFT and MD are in good agreement.

MD Method
Here, we carried out extensive fully-atomistic MD/DFT simulations to investigate the PH mechanical properties. The classical reactive (ReaXFF) MD simulations were performed using the LAMMPS code [14]. ReaxFF is a modern force field designed to be a bridge between quantum and empirical chemical methods. It is an interatomic potential that allows the breaking (C − C covalent bond) and formation of chemical bonds. It is appropriate for studying the elastic properties of nanostructures in the analysis of the mechanical failure.
The ReaxFF parameter set used here is the one described in [15,16]. All MD snapshots presented here were rendered with the molecular visual dynamics (VMD) software [17].
PH has a 2D shape composed of 5-6-8 carbon rings, as illustrated in Figure 1. In this monolayer, the lattice vectors are a 1 = 5.85 Åx and a 2 = 3.78 Åŷ. The bond distances are Å. In the MD simulations, it was subjected to uniaxial stresses along the x and y directions (see Figure 1). The temperature value ranges from 10 up to 1200 K. The PH structural model used in the MD simulations has a dimension 97 × 95 Å, being composed of 3536 carbon atoms. Its thickness has approximately 3.34 Å (i.e., one carbon atom thick). The used covalent interaction cutoff distance was 1.6 Å.
In the MD simulations, Newton's equation of motion was integrated by velocity Verlet algorithm [18]. The time step in all simulations is 0.05 fs. Before the stretching dynamics, the structure was thermalized using a Nose-Hoover thermostat during 10000 fs [19]. To eliminate any residual lattice stress before the stretching dynamics, an NPT simulation [20] was also carried out during 10000 fs. The elastic properties and fracture patterns at different temperatures (between 10-1200 K) were analyzed by stretching the PH structure along the x and y directions. We used a constant engineering tensile strain rate of 10 −6 /f s.

DFT Method
The DFT simulations were carried out using the SIESTA code [21]. The calculations were carried out within the scope of the generalized gradient approximation (GGA) with the exchange-correlation potential PBE [22]. The Monkhorst-Pack [23] grid of 8 × 8 × 1 was considered in the self-consistent calculations. The double-ζ polarization (DZP) basis set was used. The standard value of 0.05 Ry for orbital cutoff [24] and the mesh cutoff of 400 Ry were also considered [25]. We performed PH uniaxial stretching also through DFT calculations. A PH supercell was deformed considering a strain that ranged from 1 up to 30%. This strain was applied by changing the lattice vectors corresponding to the stretching direction. The

Results and Discussion
We begin our discussion by presenting the PH bond (Figure 2 (a)) and angle (Figure 2 (b)) distributions after thermalization. As mentioned above, the monolayers were heated up in the range from 10 K to 1200 K. For low temperatures, we observed two peaks in the bond distribution pattern (Figure 2(a)) centered at 1.44 Å and 1.46 Å. This distribution pattern is consistent with the C-C bond length and they are within the range of the original values for the bond lengths shown in Figure 1 and reported in reference [12]. The increase in temperature leads to larger displacements of the carbon atoms that are imposed by thermal vibrations. In this way, the average bond length of the carbon atoms increases, and the two peaks that were formerly seen at 10 K are no longer distinguishable. This trend can be observed in the curves for temperatures ranging from 100 to 1200 K. However, even upon the structure heating, the new peak remains approximately between 1.4-1.5 Å for all temperatures above 10 K. Importantly, we did not observe any bond break or new structural arrangements for temperatures up to 1200 K, as illustrated in Figure 3 (a). This fact suggests that the PH is structurally stable in this temperature range. However, we can observe the loss of the initial planar character (see Figure 3 (b)) The angle distribution for low temperatures reveals five peaks around 106 • , 109 • , 125 • , 142 • , and 145 • . But the heating spread the peaks, and only three of them can be distinguished around 109 • , 125 • and 144 • . These results further support that PH is also dynamically stable, corroborating previous results reported in the literature [12] through phonon spectra calculations.
We discuss now the PH mechanical properties calculated through MD and DFT simulations. Figure 4 shows the stress-strain curves for temperatures varying from 10 up to 1200 K obtained from the MD simulations, considering uniaxial tensile loading applied along the 4(a) x and 4(b) y-directions. The green curves in Figure 4 show the corresponding DFT results.
The first noticeable feature in this figure is that PH has a significant degree of anisotropy regarding the fracture strain (i.e., the critical strain value, ε C ). PH is more resilient to tensile loading along the y-direction than along the x-one. Another result that we can observe  contrasted to the MD ones for both directions. We attribute the observed differences to the distinct methodologies used in the DFT and MD simulations. For example, DFT results are for zero temperatures. Furthermore, the stretching process used in the MD simulations is dynamic, while in the DFT approach, it is static and based on a simulation protocol that involves stretching-relaxing-stretching the sheet to a predetermined strain value. Therefore, we believe that successive PH relaxations to energy minimum states, under a given strain value, have contributed to increasing its U S.
We summarize the results for the elastic constants in We calculated the virial stress σ C along the stretching directions according to the approach described in reference [26]. 600, 900, and 1200 K are close to results obtained with another graphene-like allotrope that has a lattice structure similar to PH, i.e., the popgraphene monolayer [27]). Like Y M , the U S, σ C , and ε C decrease smoothly with increasing temperature. These last quantities are slightly larger for the y-direction than for the x-one. This can be attributed to the different orientations, along the x and y-directions, of the PH C − −C bonds. We can see from Figure   1(b) that there is one bond parallel to the x-direction (C 4 − C 8 ) and another four bonds quasi-parallel to it (C 1 − C 2 , C 2 − C 3 , C 5 − C 6 and C 6 − C 7 ). On the other hand, there are two bonds parallel to the y-direction (C 1 − C 5 and C 3 − C 7 ).
The PH fracture process from the DFT simulations is shown in Figure 5. Figure 5(a) illustrates the monolayer at null stress. In Figure 5(b), We can note that when the uniaxial strain is along the x-direction, the breaking bonds involve the ones that are parallel to the stretching direction (i.e., bond types C 4 -C 8 ). On the other hand, in the case of the y-direction stretching (see Figure 5(c)), the fracture starts from the bonds parallel to this direction (i.e., bond types C 2 − C 6 ). Each curve refers to a bond type presented in the PH unit-cell, as shown in Figure 1(b). The PH bond lengths evolution as a function of the applied uniaxial tensile strain is shown in Figure 6. The bond C 4 -C 8 (parallel to the x-direction strain) experiences the largest variation in length and are related to the fracture at 15% strain (orange vertical traced line in 6(a)). However, the C 1 -C 5 bond (perpendicular to the stretching direction) has practically no variation in length. It is interesting to note that the C 2 -C 6 bond, perpendicular to the axis of the applied strain, has a slight decrease in its length. For the y-direction case ( Figure   6(b)), the bonds C 1 -C 4 , C 2 -C 6 and C 1 -C 5 experience considerable variations in length, the largest being the bond C 1 -C 5 (parallel to the y-direction). The bonds C 1 -C 2 and C 4 -C 8 remain practically constant throughout the tensile process. Bond fracture is seen at 19% strain (dashed vertical orange line in Figure 6(b)).
Finally, we analyzed the PH fracture pattern from the MD simulations. Figure 7 shows representative MD snapshots for the fracture process at 300K. Figures 79(a-e) and 79(f-j) indicate the cases for the stretching along the x and y-directions, respectively. The color scheme in these figures represents the spatial distribution of the von Mises (VM) stress per-atom values [28]. Red and blue colors indicate the high-and low-stress accumulation, respectively. Importantly the VM values provide useful information on the fracture process since they reveal the fracture point or region. Details about the VM calculations can be found in the LAMMPS manuals [29] and in the reference [26].
From Figure 7, we can see that PH when subjected to a critical strain, goes directly from elastic to completely fractured regimes. This process occurs with no plasticity stages between these two regimes. Importantly, graphene presents a similar fracture process [30]. In the case of the tensile loading along the x-direction, the first bonds to break are the ones parallel to the direction of the applied strain (see Figure 7(e)). After that, the lattice structural collapse occurs abruptly nearby one of the edges at 6.85 % strain, as shown in Figure 7(d).
Since the lattice has periodic boundary conditions, one can see the fragmentation and the formation of linear atomic chains (LACs) on the right edge, whereas on the left one, we can see the appearance of vacancies and dangling bonds. In the case of the tensile loading along the y-direction (see Figures 7(f-j)), we can see that the first stages of the lattice induce the formation of a crack (i. e, a fractured region in (see Figures 7(h))) instead of a single bond breaking. The crack propagation is fast, and the lattice becomes separated into two parts at around 11.5% of strain. In this case, we can also see the LAC formation.

Conclusions
In summary, we have investigated the mechanical properties and fracture patterns of GPa.nm), 40.09 TPa (13.43 GPa.nm), 9.75% at 300 K. These values gradually decrease with increasing temperature. An anisotropy for the PH mechanical behavior was observed. For the y-direction, we obtained higher values for the elastic constants when compared to the x-direction. PH when subjected to a critical strain, goes directly from elastic to completely fractured regimes. This process occurs with no plasticity stages between these two regimes.
Importantly, graphene presents a similar fracture process. The crack propagation is fast and the lattice becomes separated into two parts with the presence of linear atomic chains (LACs).

Acknowledgements
This