Tetragraphene-Based Nanotubes Under Temperature Effects: Atomistic Insights into Fracture Toughness and Nanostructural Degradation via Reactive Molecular Dynamics

1. Introduction

The scientific community has, in recent years, significantly emphasized the superiority of nanostructured materials over bulk alternatives, exploring their versatile applications within numerous technological and scientific sectors [2,180]. Nanoscience, together with nanotechnology, focuses on analyzing materials possessing dimensions that fall within the nanometer regime [3]. These systems exhibit distinct biological, chemical, physical, and engineering characteristics, which are determined by length scales often ranging from 10 to 200 Angstroms [4].

At the nanoscale, a multidisciplinary environment emerges where chemistry, physics, and engineering converge, facilitating the development of unprecedented properties in condensed matter [5,6]. Leveraging these unique physical characteristics, it is possible to manipulate individual atoms, leading to innovations such as molecular robots and engines capable of mimicking cellular functions [7,8,9]. An effective approach for analyzing the progress of these scientific fields involves thoroughly reviewing the increasing body of research literature focused on the properties of new nanostructures.

Based on this scientific description, scientists has done a lot of research work whith aims to understand the condensed matter at nanometric scale, theoretically and experimentally [10,11,12,13,14,15]. The manufacturing processes of carbon-based nanostructures and materials in laboratory scale (theory and experiment) may be categorized throughout changes carbon-based nanostructures and materials. It is based on the allotropy of carbon that we have achieved significant advances in science and technology. Owing to their outstanding physicochemical characteristics, carbon atoms and their diverse allotropic forms have profoundly altered the landscape of the electronic and optoelectronic sectors, sparking extensive research into their prospective roles in modern nanoscience [16,17,18,19,20,21,22,23,24]. A primary factor behind these versatile features is the chemical valency of carbon, which enables the formation of numerous distinct nanostructural networks. Over the past few decades, a variety of novel carbon-based materials and nanostructures have been isolated, most notably zero-dimensional fullerenes [25], one-dimensional carbon nanotubes (CNTs) [26], and two-dimensional isolated graphene layers [35]. The discovery of these unique bonding configurations effectively initiated a new epoch of technological breakthroughs and fundamental scientific exploration.

The experimental isolation of a single-layer graphene sheet via mechanical exfoliation in 2004 fundamentally redefined the boundaries of condensed matter physics and engineering, an achievement subsequently honored with the Nobel Prize [35]. Structurally, graphene consists of a flat, monomolecular plane of carbon atoms densely arranged in a two-dimensional ($2D$) honeycomb topology. This unique configuration imparts extraordinary physical characteristics, establishing it as a highly promising material for next-generation technological applications [27]. The structural adaptability of this network directly determines its unique electronic features; specifically, the $s{p}^2$ hybridization involving one s and two p orbitals yields a robust trigonal planar morphology, where carbon atoms are linked by strong $igma$ bonds with an equilibrium bond length of $1.42$ Å. Consequently, graphene displays a remarkable intrinsic electron mobility of up to $2\times {10}^5$ ${\mathrm{cm}}^2/\mathrm{V}dot\mathrm{s}$ and an exceptional electrical conductivity around $2\times {10}^4$ S/cm stemming from its zero-bandgap nature [28]. Furthermore, its thermal conductivity ranges between $(4.84\pm 0.44)\times {10}^3$ and $(5.30\pm 0.48)\times {10}^3$ W/mK [29], complemented by distinct optical reflectance and transmittance profiles [30]. Ambient-temperature ferromagnetism can also be induced through surface molecular adsorption [31], while its mechanical robustness is evidenced by an elastic Young’s modulus reaching approximately 1 TPa (${10}^3$ GPa) [32,33]. However, despite these unparalleled intrinsic traits, the complete implementation of graphene in scalable nanoelectronic devices is heavily restricted by its lack of an energy bandgap [34].

To overcome these nanostructural limitations, numerous alternative allotropes and low-dimensional materials have been proposed in recent years, aiming to replicate or improve upon the foundational merits of graphene. These include alternative elemental monolayers and advanced carbon networks such as two-dimensional covalent carbon nitride nanosheets [36], mixed $sp$-$s{p}^2$ coordinated graphynes [37], boron nitride nanosheets [38], silicon-based silicene [39], penta-graphene [40], phagraphene [41], popgraphene [42], twin graphene [43], octagraphene [44], T-carbon [45], D-carbon [46], cubic $C96$ [47], $T5$-carbon [48], $P2221-C8$ [49], trigraphene [50], tetrahexcarbon [51], $C_5{N}_2$ anodes [52], $C568$ [53], $Cmme-SnS$ [54], $C16$ and $C24$ networks [55], $TPDH$-graphene [56], $Me$-graphene [57], $TPO12$-graphene [58], pentagraphyne [59], two-dimensional $Tph-Cn2$ [60], naphthylene-$\gamma$ [61], $pan-C_{2(n+1)}$ [62], thgraphene [63], tetrahex carbides [64], pentahexoctite [65], hexatetra-carbon [66], tower carbon [67], novamene [68], orthorhombic $C14$ carbon [69], $tP40$ carbon [70], $C20$-Tcarbon [71], $Pa-C10$ [72], graphdiyne architectures [73], $hex-C558$ [74], superprismane [75], orthorhombic $C32$ [76], four-penta-graphenes [77], R-graphyne [78], biphenylene monolayers [79], $PmmaC24$ [80], $yne-C16$ [81], $QPHT$-graphene [82], $HSH$-carbon [83], dodecanophene systems [84], M-graphene [85], $rC14$ [86], Penta-$C20$ [87], $Hp-C17$ [88], $THD-C$ sheets [89], pentaheptite diamond [90], $CP-C20$ [91], $C-57$ carbon [92], poly-dodecahedrane [93], $o-C240$ structures [94], $TH$-graphyne [95], janus graphene [96], $THFS$-carbon [97], $Cto$-graphene [98], tetragonal $C64$ [99], twin T-graphene [100], superhard $Cmca-C16$ [101], spiro-carbon [102], diamond $D5$ [103], $oI20$-carbon [104], monoclinic $C12$ [105], azugraphene [106], $Pc-C10$ [107], $LC567$ citecai2023lc, $Si96$ topologies [109], $PHOTH$-graphene [110], $PCF$-graphene [111], C-centered orthorhombic $C8$ [112], $Bct-C5$ [113], Q-carbon [114], $HSH-C10$ [115], escherynes [116], petal-graphyne citelima2025petal, orthorhombic $oC48$ [118], irida-graphene [119], and sun-graphyne [120].

Emerging directly from this nanostructural design paradigm is tetragraphene, a novel carbon allotrope characterized by a unique mixture of $s{p}^2$ and $s{p}^3$ hybridization. When rolled into one-dimensional architectures, these sheets form tetragraphene-based nanotubes (TGCNTs). Theoretical investigations into the intrinsic physical and electronic properties of TGCNTs have already been initiated within the scientific community, most notably through Density Functional Theory (DFT) calculations by Brandão et al $\left(2025\right)$ [121], regarding topological and strain-induced bandgap engineering.

However, such quantum-mechanical approaches are intrinsically limited to small, idealized systems at absolute zero temperature. This constraint creates a critical “DFT gap”, leaving a vital need to understand how these one-dimensional nanostructures behave under realistic, large-scale, and dynamic conditions. To bridge this gap, classical molecular dynamics (CMD) simulations are uniquely required, as they allow for the accurate exploration of large spatial scales, high-temperature fluctuations, fracture toughness, and complex dynamic deformation pathways. Consequently, this work employs reactive CMD simulations with the AIREBO-Morse potential to systematically investigate TGCNTs under uniaxial tensile loading across a wide thermal range (300 to 2100 K). Crucially, this study uncovers unprecedented phenomena not previously reported in the literature, including a strict thermal degradation threshold above 1500 K and a striking chirality-dependent divergence in fracture modes. Specifically, we report a stress-induced nanostructural phase transition in $(N,0)$ TGCNTs, marked by a constant-stress plateau up to 20% strain and subsequent strain hardening, contrasting sharply with the brittle fracture observed in $(0,N)$ configurations. Ultimately, these novel insights into the fracture toughness, thermal stability limits, and dynamic deformation mechanisms of TGCNTs provide essential design guidelines for their integration into next-generation flexible electronics, nanocomposites, and nanoscale energy-absorption systems.

2. Geometric Generation and Configuration of Tetragraphene-Based Nanotubes (TGCNTs)

The one-dimensional (1D) spatial morphology of a tetragraphene carbon nanotube (TGCNT) is established by mathematically rolling a single layer of tetragraphene. This structural mapping is governed by the chiral indices $(n,0)$ and $(0,n)$. The primary geometric metrics include the tube diameter and a specific zigzag-like orientation. These dimensional factors, alongside the wall thickness, play a vital role in dictating the mechanical and elastic response of the resulting TGCNTs.

While the fundamental concept of rolling a 2D sheet is conceptually analogous to standard graphene, the explicit atomic generation of TGCNTs differs due to the distinct symmetry of the parent tetragraphene lattice. Formally, the circumferential chiral vector (${\mathbf{C}}_h$) is expressed as a linear combination of the two-dimensional rectangular basis vectors (${\mathbf{a}}_1$ and ${\mathbf{a}}_2$) shown in Figure Figure 1 (a) and, the Equation 1 [121]:

$${\mathbf{C}}_h=n{\mathbf{a}}_1+m{\mathbf{m}}_2\equiv (n,m),$$

(1)

where n and m are integer values that define the nanotube chirality, with the basis magnitudes given as $a_1=4.56$ Å and $a_2=6.16$ Å. This vector directly governs the cross-sectional dimensions of the system. Consequently, the diameter of the TGCNT ($D_{tgcnt}$) is computed from these index values via the following relation [121]:

$$D_{tgcnt}={\displaystyle \frac{|{\mathbf{C}}_h|}{\pi }}={\displaystyle \frac{t_{tg}}{2\pi }}\sqrt{3(n^2+m^2+nm)},$$

(2)

where $t_{tg}$ represents the characteristic lattice parameter of the tetragraphene nanostructure [121].

The nanostructural translation along the longitudinal axis is described by the translational vector ($\mathbf{T}$), which is oriented perpendicular to ${\mathbf{C}}_h$ (see Figure 1 (a)) [121]:

$$\mathbf{T}=t_1{\mathbf{a}}_1+t_2{\mathbf{a}}_2.$$

(3)

The integer coefficients $t_1$ and $t_2$ are derived from the orthogonality condition requiring the dot product to vanish, ${\mathbf{C}}_{h}dot\mathbf{T}=0$. Given that ${\mathbf{a}}_{1}dot{\mathbf{a}}_2=0$ for this orthogonal system, the coefficient ratio simplifies to [121]:

$${\displaystyle \frac{t_1}{t_2}}=-{\displaystyle \frac{m}{n}}.$$

(4)

The total length of the periodic nanotube section is defined by the magnitude $L=\left|\mathbf{T}\right|$. Together, the boundary vectors ${\mathbf{C}}_h$ and $\mathbf{T}$ completely outline the boundaries of the TGCNT unit cell [121]. Because of the rectangular symmetry of the underlying lattice, this setup allows for the creation of zigzag-like configurations denoted by $(n,0)$ or $(0,n)$ indices [121]. The explicit structural attributes of the specific TGCNTs investigated in this study are summarized in in Table 1.

3. Reactive Molecular Dynamics Methodology and Computational Details

To evaluate the nanomechanical behavior and distinct physical features of tetragraphene-based nanotubes (TGCNTs), atomistic calculations were performed using the reactive classical molecular dynamics simulations (CMD) framework. These computational experiments were implemented within the open-source LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) software suite [141,142,143]. As a robust and widely adopted platform for atomistic modeling, LAMMPS enables the systematic replication of nanomechanical deformations and stress-strain responses in low-dimensional nanostructures. This methodology provides critical insights into structural evolution and nanoscale failure mechanisms that remain challenging to probe through conventional experimental techniques [144].

In this study, the mechanical characteristics of the TGCNTs are probed by imposing continuous physical deformations on the model and monitoring the resulting nanostructural feedback. The specific virtual testing protocol utilizes uniaxial tensile loading regimes, wherein the dimensions of the simulation box are extended along a single axis under a fixed, uniform strain rate. This procedure allows for the direct derivation of stress-strain relations. From the resulting data curves, primary mechanical metrics, specifically the Young’s modulus, ultimate tensile strength ($UTS$), fracture strain threshold, and total energy absorption capacity (toughness), are quantitatively determined [145,146]. By the fact that TGCNTs one dimentional nanostrucutures are composed by carbon atoms, a carbon atom, like all matter, is fundamentally governed by quantum mechanics [147]. However, it can be a useful and accurate approximation to treat a carbon atom as a classical particle under certain conditions, primarily when observing microscopic statistical averages properties, under external effects such as temperatures [148].

Given that carbon atoms within this regime can be accurately treated as classical point masses, the complete set of atomic interactions is governed by parameterised force fields [149]. These interatomic potentials serve as the foundational mathematical framework that evaluates the total potential energy of the configuration based on spatial coordinates, subsequently yielding the net forces exerted on each particle. As a core component of atomistic modeling platforms like classical molecular dynamics (CMD), these formulations are utilized here to quantify critical nanomechanical parameters, including the elastic response, maximum load capacity, and nanostructural failure pathways of the TGCNTs [150,151]. Specifically, the Adaptive Intermolecular Reactive Empirical Bond Order Morse (AIREBO-Morse) potential, based on the foundational framework formulated by Stuart et al. [191], was implemented to govern the covalent interactions within the carbon network. This formulation incorporates a Morse long-range term to describe intermolecular interactions, which has demonstrated remarkable success in accurately describing the nanomechanical behavior, stability, and nanostructural properties of diverse carbon-based systems, as thoroughly demonstrated in Appendix B of this research article. As a bond-order-dependent potential, the AIREBO-Morse framework effectively captures continuous bond breaking and formation under high strain rates without requiring predefined bonding topologies. Consequently, it achieves superior fidelity in resolving localized stress concentrations, ultimate tensile limits, and subsequent fracture dynamics inherent to the investigated TGCNTs.

To observe and animate the real-time atomistic failure dynamics of the TGCNTs, the VMD (Visual Molecular Dynamics) software suite was utilized [160]. The underlying Newtonian equations of motion for the carbon atoms are integrated temporally via the standard Verlet scheme [161], adopting a conservative time step of $0.05$ fs. To mitigate finite-length artifact errors, a periodic boundary condition was enforced strictly along the longitudinal z-axis of the $1D$ system. Conversely, to suppress any spurious residual stresses during the stretching phase, an expansive simulation box was constructed to surround the TGCNT with a 100 Å vacuum region along both the x and y lateral dimensions. This layout guarantees a zero-stress state across the effective cross-sectional thickness of the nanotube.

The influence of thermal conditions was examined systematically across a wide temperature range, specifically at 300, 600, 900, 1200, 1500, 1800, and 2100 K. Every molecular trajectory was initiated with a rigorous nanostructural energy minimization phase. This pre-relaxation step employed a Polak-Ribiere conjugate gradient ($CG$) algorithm (specified as style “cg”). The energy and force convergence thresholds for these initial iterations were set to an energy tolerance ($etol$) of $1.0\times {10}^{-4}$ and a force tolerance ($ftol$) of $1.0\times {10}^{-6}$ $\mathrm{kcal}dot{\mathrm{mol}}^{-1}dot{\text{Å}}^{-1}$ [161,163,164].

Following minimization, the configurations were equilibrated to a zero-stress state for $0.01$ ns using an isothermal-isobaric Nosé-Hoover ($NPT$) barostat protocol with a pressure damping coefficient of 50 [165]. Subsequently, a $0.01$ ns thermalization phase was performed under a canonical Nosé-Hoover ($NVT$) ensemble with a temperature damping coefficient of 5 [166]. Uniaxial loading was subsequently applied by continuously expanding the boundaries of the simulation domain along the periodic z-direction. These tensile tests were conducted within the $NVT$ ensemble at a constant engineering strain rate of $\zeta ={10}^{-6}$ ${\mathrm{fs}}^{-1}$. Under this regime, the instantaneous TGCNTs length (L) evolves according to the linear relationship [167,168]:

$$L\left(t\right)=L_0+L_0\zeta dt,$$

(5)

where $L_0$ represents the baseline length and $dt$ denotes the cumulative simulation time. Finally, the intrinsic Young’s modulus ($Y_{MOD}$) was extracted from the linear elastic regime using the classical definition [167,168]:

$$Y_{MOD}={\displaystyle \frac{{\sigma }_{zz}}{{ϵ}_z}},$$

(6)

where ${\sigma }_{zz}$ corresponds to the component of the Virial stress tensor and ${ϵ}_z$ signifies the strain experienced along the uniaxial z-axis [167,168].

The Young’s modulus is obtained in a linear regime (3% strain). The standard errors are obtained from a single linear regression fit. The macroscopic stress tensor (specifically, the Cauchy stress) is commonly and effectively estimated using a generalization of the virial theorem in molecular simulations and statistical mechanics. This definition, known as the virial stress, provides a direct link between microscopic atomic behavior and macroscopic continuum properties [169,170,171]:

$$\sigma =V^{-1}\sum _{i\in V}\left[-m^{\left(i\right)}{\mathbf{v}}^{\left(i\right)}\otimes {\mathbf{v}}^{\left(i\right)}+{\displaystyle \frac{1}{2}}\sum _{i\ne j}{\mathbf{r}}^{\left(ij\right)}\otimes {\mathbf{f}}^{\left(ij\right)}\right],$$

(7)

where $m^{\left(i\right)}$ are the mass and ${\mathbf{v}}^{\left(i\right)}$ the velocity vector of atom i. ⊗ denotes the tensor product of two vectors. ${\mathbf{r}}^{\left(i\right)}$ denotes the position of atom i. ${\mathbf{r}}^{\left(ij\right)}={\mathbf{r}}^{\left(j\right)}-{\mathbf{r}}^{\left(i\right)}$ is the position vector of atom j relative to that of atom i, and ${\mathbf{f}}^{\left(ij\right)}$ is the interatomic force exerted on atom i by atom j, where:

$${\mathbf{f}}^{\left(ij\right)}={\displaystyle \frac{\partial P}{\partial r^{\left(ab\right)}}}.{\displaystyle \frac{\partial {\mathbf{r}}^{\left(ij\right)}}{\partial r^{\left(ab\right)}}},$$

(8)

where, P is the potential energy of the atomic ensemble. V is the volume of the TGCNTs one-dimentional nanostructures. $V=A.\left(thickness\right)$, A is the surface area and thickness of the single layer of Tetragraphene. The von Mises stress tensor distribution per atom is defined as [172]:

$${\sigma }_{\mathrm{vonMises}}=\sqrt{\left[{\displaystyle \frac{{({\sigma }_{xx}-{\sigma }_{yy})}^2+{({\sigma }_{yy}-{\sigma }_{zz})}^2+{({\sigma }_{zz}-{\sigma }_{xx})}^2+6({\tau }_{xy}^2+{\tau }_{yz}^2+{\tau }_{zx}^2)}{2}}\right]},$$

(9)

The symbols ${\sigma }_{xx}$, ${\sigma }_{yy}$, and ${\sigma }_{zz}$ indicate the stresses in the x, y, and z directions, respectively. The ${\tau }_{xy}$, ${\tau }_{yz}$, and ${\tau }_{zx}$ represent the shear stress components. The concept of von Mises stress, primarily a continuum mechanics theory for predicting yielding in ductile nanostrucutres/materials, is adapted for atomic-level analysis in classical molecular dynamics (CMD) simulations [173,174].

3.1. Evaluation and Formulation of Poisson’s Ratio in TGCNTs

The transverse deformation behavior of the tetragraphene carbon nanotubes (TGCNTs) is quantified by determining their Poisson’s ratio through the reactive classical molecular dynamics (CMD) framework. The resulting values exhibit a clear dependence on the geometric features of the nanotubes, including their specific diameter and chirality. To evaluate this property, the cylindrical nanostructures are subjected to an unconstrained uniaxial tensile load applied along the longitudinal z-axis. Because no radial restrictions are imposed on the system, the nanotube is free to contract laterally under deformation. Consequently, the intrinsic Poisson’s ratio ($\nu$) is mathematically defined by the following relation [175,176]:

$$\nu =-{\displaystyle \frac{{ϵ}_{\mathrm{transverse}}}{{ϵ}_z}},$$

(10)

$$\begin{array}{c}\hfill \nu =-li{m}_{\epsilon =0}\left({\displaystyle \frac{{\epsilon }^{\prime }}{\epsilon }}\right),\end{array}$$

(11)

where, $\nu$ Poisson’s ratio, ${\epsilon }^{\prime }$ radial strain and $\epsilon$ uniaxial strain applied in TGCNTs. The Poisson’s ratio, $\nu$, of TGCNT $(14,0)$ and $(0,10)$ are calculated by (Equation 12) [175,176]:

$$\begin{array}{c}\hfill \nu =-{\displaystyle \frac{{\epsilon }_R}{{\epsilon }_z}},\end{array}$$

(12)

In this formulation, ${ϵ}_R$ represents the induced radial strain, while ${ϵ}_z$ denotes the applied longitudinal deformation. To compute these spatial variations, the total length (L) of the representative systems, specifically the $(14,0)$ and $(0,10)$ TGCNTs, is partitioned into uniform volumetric segments of thickness $\Delta L=L/n$. The corresponding cross-sectional radial strain is evaluated using the following expression:

$${ϵ}_R={\displaystyle \frac{R_{{ϵ}_z}^{\left(\mathrm{avg}\right)}-R_0^{\left(\mathrm{avg}\right)}}{R_0^{\left(\mathrm{avg}\right)}}},$$

(13)

where $R_{{ϵ}_z}^{\left(\mathrm{avg}\right)}$ represents the mean tube radius evaluated at a given uniaxial strain level ${ϵ}_z$, and $R_0^{\left(\mathrm{avg}\right)}$ signifies the baseline radius at thermal equilibrium (${ϵ}_z=0$). The mean strained radius is determined by averaging across all discrete nanostructural segments as follows:

$$R_{{ϵ}_z}^{\left(\mathrm{avg}\right)}={\displaystyle \frac{1^n}{n_{i=1}}}{R}_i^{\left(\mathrm{slab}\right)}{|}_{{ϵ}_z}.$$

(14)

Here, $R_i^{\left(\mathrm{slab}\right)}{|}_{{ϵ}_z}$ defines the instantaneous average radius of the i-th circular segment under the longitudinal strain ${ϵ}_z$. This segment-specific radius corresponds to the mean spatial distance from each constituent carbon atom to the local center of mass of that particular slice, computed via:

$$R_i^{\left(\mathrm{slab}\right)}={\displaystyle \frac{1^N}{M_{\alpha =1}}}{r}_{\alpha },$$

(15)

where the radial distance $r_{\alpha }$ for each carbon atom is given by the geometric relation:

$$r_{\alpha }=\sqrt{{\left(x_{\beta }-x_{\mathrm{CM}}\right)}^2+{\left(y_{\beta }-y_{\mathrm{CM}}\right)}^2}.$$

(16)

In Equation 16, $x_{\beta }$ and $y_{\beta }$ correspond to the planar coordinates of the $\alpha$-th carbon atom, M represents the total atomic population assigned to each individual slab, and the coordinates $x_{\mathrm{CM}}$ and $y_{\mathrm{CM}}$ define the instantaneous center of mass for that specific segment. For all calculations, a fixed slice thickness of $\Delta L=1.5$ Å was adopted, ensuring a statistically robust and sufficient number of atoms per segment to accurately map the radial geometry [177].

4. Results, Analysis, and Discussion

The systematic nanomechanical evaluation of tetragraphene-based nanotubes (TGCNTs) via reactive classical molecular dynamics (CMD) simulations reveals a correlation and new finds of their atomic nanostructure nanomechanical response under uniaxial strain. The unique mixture of $s{p}^2$ and $s{p}^3$ hybridization within the tetragraphene lattice introduces distinct deformation pathways that are highly sensitive to external constraints, such as temperature and loading conditions.

CMD simulations using the AIREBO-Morse potential reveal a strict chirality-dependent, anisotropic nanomechanical behavior in tetragraphene-based nanotubes (TGCNTs), distinguishing between brittle $(0,N)$ and ductile $(N,0)$ failure modes (see Table 1, Figure 2, Figure 3, Figure 4 and, Figure 5). While $(0,N)$ TGCNTs exhibit abrupt (see Figure 2), brittle failure, $(N,0)$ TGCNTs display a notable 20% stress plateau indicative of a strain-induced phase transformation followed by significant strain-hardening and superior fracture strain (see Figure 4). According to research by Brandão et al. $\left(2025\right)$ [121], the severe nanomechanical anisotropy of tetragraphene-based nanotubes (TGCNTs) arises from the alignment of parallel $s{p}^2$ dimers and buckling $s{p}^3$ connections within its trilayer nanostructures. For $(n,0)$ TGCNTs, the $s{p}^3$ bonds align with the tension axis, allowing for a flattening phase transformation and a $20\%$ superelastic strain plateau, while in $(0,n)$ TGCNTs, the rigid $s{p}^2$ dimers align with the tension axis, resulting in brittle fracture due to limited nanostructural rearrangement. Most notably, our findings demonstrate a critical thermal threshold at 1500 K, beyond which severe nanostructural degradation occurs, leading to a complete loss of nanostructural stability. This thermal vulnerability underscores the delicate balance of the hybrid bonding network, setting a clear operational boundary for TGCNTs in high-temperature applications while highlighting the importance of environmental factors in their nanostructural integrity (see graphical representation (stress-strain) in Figure 3 and Figure 5). This thermal degradation is explicitly reflected in the simulated stress-strain curves (Figure 3 and Figure 5), which reveal a sharp and systemic reduction in maximum tensile strength, Young’s modulus (discussed later), and fracture strain as the temperature approaches 1500 K. Above this critical threshold, the characteristic stress-strain profile collapses completely into an amorphous low-stress response, characteristic of a nanostructural melting phase.

The complete loss of nanostructural stability above 1500 K is driven by two main factors: $\left(i\right)$ High-Energy Bond Rupture: Unlike traditional graphene which contains only highly stable $s{p}^2$ bonds, the tetragraphene lattice contains inherently weaker, buckled $s{p}^3$ hybridized carbon-carbon bonds. At temperatures exceeding 1500 K, the intense thermal vibrations supply kinetic energy that matches or exceeds the activation energy barrier required to break these $s{p}^3$ cross-linking bonds. $\left(ii\right)$ Amorphization and Rehybridization: Once the weaker $s{p}^3$ connections break, the highly strained network undergoes rapid local nanostructural collapse. The broken bonds trigger an irreversible cascade of rehybridization, causing the orderly crystalline nanotube to melt into a disordered, amorphous carbon pipe that cannot sustain nanomechanical loads (see snapshots of CMD results in Figure 6 and Figure 7). This thermal threshold for tetragraphene-based nanotubes closely matches established carbon nanostructure research. Theoretical studies on novel, non-ideal allotropes containing mixed $s{p}^2$/$s{p}^3$ networks, such as biphenylene, penta-graphene, and defective nanotubes, consistently predict a sharp decline in nanomechanical stiffness and premature melting between 1200 K and 1600 K. This behavior is primarily driven by the nanostructural vulnerability of $s{p}^3$ hybridizations and non-hexagonal rings [124,179,181]. Furthermore, high-temperature transmission electron microscopy (TEM) experiments and atomistic simulations on multi-walled carbon nanotubes with high defect densities show that structural degradation and structural reconstruction begin rapidly within this identical thermal window. While pristine graphene can withstand much higher thermal loads, the built-in nanostructural strain of the tetragraphene lattice lowers its melting point, making 1500 K the clear physical limit for its nanostructural survival [182,183,184].

The nanomechanical response of single-walled tetragraphene-based nanotubes (TGCNTs) under uniaxial tensile loading was systematically investigated via classical molecular dynamics (CMD) simulations utilizing the AIREBO-Morse potential. The graphical representation in Figure 2 and Figure 4 (a)-(l) illustrates the stress-strain profiles for diverse configurations, comparing the $(0,N)$ series from indices $(0,3)$ to $(0,13)$ with the $(N,0)$ series from $(3,0)$ to $(13,0)$. The results reveal a striking contrast in failure modes dictated by tube chirality. The $(0,N)$ TGCNT series exhibits a characteristically brittle failure mode, which is marked by a highly linear elastic regime followed by a sudden nanostructural fracture with negligible plastic deformation. Conversely, the $(N,0)$ TGCNT series demonstrates remarkable ductility and superelasticity. This behavior is evidenced by a distinct plateau effect that maintains a constant stress level up to 20% strain, followed by clear strain hardening until ultimate fracture beyond 40% strain, indicating a stress-induced nanostructural phase transition. To provide a comprehensive quantitative overview, key nanomechanical properties were extracted from the simulation curves. The Young’s modulus ($Y_{MOD}$) was calculated by applying a linear regression within the initial 3% uniaxial strain regime for all configurations, with fitting errors precisely estimated using the Xmgrace plotting software [178]. Additionally, the Ultimate Tensile Strength (UTS) and the critical strain (${ϵ}_C$) were determined at the peak point of failure. Table 2 summarizes the complete set of calculated values for the Young’s Modulus (GPa.Å), UTS (GPa.Å), and critical strain across all studied TGCNTs (see Table 1). In graphical representation (Figure 3 and Figure 5) of stress-strain uniaxial load nanomechanical behavior shown the atomistic deformation response of single-walled tetragraphene-based nanotubes (TGCNTs) under uniaxial tensile loading reveals distinct, temperature-dependent nanomechanical pathways between the (0,10) and (14,0) configurations from 600 K up to 2100 K. As illustrated in the stress-strain profiles (Figure 3 (a)-(f), the (0,10) TGCNT series displays a characteristically brittle failure mode defined by a highly linear elastic regime followed by sudden nanostructural fracture with negligible plastic deformation. In sharp contrast, the (14,0) TGCNT series demonstrates remarkable ductility and superelasticity, characterized by a distinct plateau effect that maintains constant stress up to 20% strain before undergoing clear strain hardening until ultimate fracture (see Figure 5). For the (14,0) system, rising temperatures accelerate the degradation of nanofracture toughness, significantly reducing the critical strain to values below 40%. Ultimately, both configurations exhibit severe nanostructural degradation at elevated thermal states, leading to a complete loss of nanostructural stability above 1500 K. This extreme thermal degradation results in an amorphous graphene phase devoid of a defined geometry, which fundamentally accounts for the melting and collapse of the one-dimensional carbon TGCNTs nanostructures.

In Figure 6 and Figure 7 we can see the fully atomistic configurations obtained from classical molecular dynamics (CMD) simulations visually capture the thermal degradation of the (0,10) and (14,0) TGCNTs across the temperature range from 600 K up to 2100 K. As observed in the snapshots, the dotted circles highlight the localized loss of the initial nanostructured configuration induced by extreme thermal effects. Furthermore, the accompanying zoomed views explicitly display the nanostructural transition into a fully amorphous carbon phase, providing clear evidence of the severe nanostructural collapse experienced by the nanotubes at these elevated temperatures. Following the detailed nanomechanical analysis of the individual stress-strain curves, the overall elastic properties and nanomechanical limits of all studied TGCNTs were statistically averaged. To provide a comprehensive comparison of their nanostructural performance, the key nanomechanical parameters are summarized in Table 2. This dataset presents the calculated values for the Young’s Modulus ($GPa·\text{Å}$), Ultimate Tensile Strength-UTS ($GPa·\text{Å}$), and critical strain ${ϵ}_C$ (%) for the diverse configurations of tetragraphene-based nanotubes (TGCNTs). The data presented in Table 2 reveals a remarkable nanomechanical anisotropy in single-walled tetragraphene-based nanotubes (TGCNTs). This direction-dependent behavior is highly sensitive to the nanotube chirality, which directly influences the bond orientation relative to the uniaxial strain axis. The distinct nanomechanical responses between the $(n,0)$ and $(0,n)$ configurations underscore the nanostructural uniqueness of the tetragraphene lattice when rolled into one-dimensional nanostructures (significant nanoechanical anisotropy). The $(n,0)$ series possesses a significantly stiffer elastic regime compared to the $(0,n)$ family. For larger tube diameters, the Young’s modulus of $(n,0)$ nanotubes stabilizes above 3100 GPa·Å. In stark contrast, the $(0,n)$ nanotubes plateau near 2350 GPa·Å. This represents an elastic stiffness difference of roughly 30% between chiralities. The maximum load-bearing capacity also demonstrates a clear anisotropic gap. The $(n,0)$ configurations routinely achieve UTS values between 640 and 700 GPa·Å. Conversely, the $(0,n)$ nanotubes exhibit much lower maximum strengths, generally remaining below 460 GPa·Å. The most dramatic manifestation of anisotropy lies in the failure limits. The ductile $(n,0)$ series stretches up to an exceptional 38.18%–41.06% strain before nanostructural failure. On the other side, the brittle $(0,n)$ series structurally fails at very early stages, showing critical strains of only 11.35%–13.75%. For the highly curved $(3,0)$, $(4,0)$, and $(5,0)$ nanotubes (small-diameter effects), the Young’s modulus experiences a noticeable dip down to 2379.90 GPa·Å. This dip indicates that strong curvature destabilizes the initial elastic resistance. A similar but less severe trend is visible in the $(0,3)$ configuration. As the tube indices increase beyond 7 (large-diameter convergence), the nanomechanical properties for both families stabilize. The values become nearly independent of the tube diameter. This independent behavior suggests that the system smoothly approaches the mechanical limit of a flat tetragraphene sheet.

To address the core nanostructural parameters of one-dimensional systems, the Young’s modulus of the TGCNTs was evaluated by considering the nanotube wall as a thin-shelled cylinder. The raw nanostructural stiffness obtained from the simulations ranges from 2714.10 to 3166.20 GPaÅ (corresponding to 271.41-316.62 GPanm). To translate this metric into a conventional three-dimensional Young’s modulus ($Y_{3D}$), the values were normalized by the intrinsic shell thickness of the nanotube wall. This wall thickness corresponds to the tetragraphene lattice amplitude of 4.571Å , which originates from its characteristic corrugated mixed $s{p}^2/s{p}^3$ hybridization pattern. This formal normalization yields an effective material Young’s modulus ranging from 593.77 to 692.67 GPa for the investigated TGCNTs, establishing a rigid nanostructural baseline that shows close agreement with previous Density Functional Theory (DFT) predictions [121]. So, all classical molecular dynamics (CMD) results utilizing the AIREBO-Morse potential show strong consistency with first-principles Density Functional Theory (DFT) calculations [121]. The calculated $Y_{Mod}$ values match the ground-state predictions from DFT very closely. This high accuracy confirms that the AIREBO-Morse potential captures the intrinsic carbon-carbon bond stretching of the tetragraphene lattice exceptionally well. The severe drop in critical strain for the $(0,n)$ nanotubes vs. the large deformation capacity of the $(n,0)$ nanotubes aligns perfectly with DFT-predicted energy landscapes [121]. This close validation proves that our CMD simulations reliably describe both the initial elastic deformation and the ultimate failure points of these novel nanomaterials.

In Figure 8 and Figure 9 display the fully atomistic molecular dynamics snapshots of the uniaxial nanomechanical loading for the TGCNT (0,10) and (14,0) configurations stretched along the z-direction, respectively. As depicted in Figure 8 (a), the (0,10) TGCNT is initially shown at a strain-free state $(\epsilon =0\%)$, establishing its pristine nanostructural equilibrium. Upon applying a uniaxial tensile load up to a strain of 10% (Figure 8 (b)), the TGCNT exhibits uniform nanostructural elongation within its linear elastic regime. The onset of nanomechanical failure is captured in Figure 8 (c), where initial bond cleavage specifically targets the carbon-carbon ($C-C$) bonds aligned parallel to the uniaxial strain direction. This localized stress concentration rapidly triggers a catastrophic, brittle nanostructural fracture, culminating in the complete cleavage of the (0,10) TGCNT into two distinct segments at a critical strain of 21.27%, as illustrated in Figure 8 (d). A remarkably interesting phenomenon observed immediately following this complete nanostructured failure is the formation of linear atomic chains (LACs) bridging the fractured ends. Under nanomechanical tensile strain, the $1D$ tetragraphene nanostructure undergoes a coordinated nanostructural transition from a crystalline lattice to a disordered state, leading to the formation of linear amorphous carbons (LACs), in strong agreement with atomistic calculations reported for tetragraphene single-layers [122]. This behavior can be fundamentally attributed to the local hybridization transitions within the TGCNT network under the specific (0,10) strain orientation, where $s{p}^2$-like or $s{p}^3$-like carbon frameworks undergo severe nanostructural reconstruction into $sp$-hybridized chain-like geometries. In panel (d), the yellow indicators specify the dynamic bond lengths of the highly elongated $C-C$ linkages right before their ultimate rupture. These chemical bond lengths within the generated LACs vary from $1.47$Å to $1.70$ Å. This significant elongation beyond the standard equilibrium $C-C$ bond distance highlights the high local tensile stress sustained by these atomic chains, a feature accurately captured by the AIREBO-Morse potential, which properly describes bond breaking and subsequent nanostructural remodeling under extreme nanomechanical loads.

In Figure 9 displays the fully atomistic molecular dynamics snapshots tracking the uniaxial nanomechanical loading of the (14,0) TGCNT configuration stretched along the z-direction. Initially, Figure 9 (a) depicts the nanotube in its pristine, strain-free state ($\epsilon =0\%$). As the uniaxial tensile load increases to strains of 20% (Figure 9 (b) and 30% (Figure 9 (c), the system accommodates large deformations without nanostructural failure. The onset of nanomechanical failure is captured at a critical strain of 40.69% (Figure 9 (d), where the inset close-up view highlights the initial cleavage of local carbon-carbon ($C-C$) bonds. Ultimate failure is achieved in Figure 9 (f), yielding a permanent, nanostructured fracture that divides the nanotube into two distinct parts, accompanied by the formation of minor linear atomic chains (LACs) that are notably less pronounced than those observed in the (0,10) counterpart. In stark contrast to the brittle (0,10) configuration, the (14,0) TGCNT demonstrates remarkable ductility and superelasticity. This response is evidenced by a distinct plateau effect maintaining a nearly constant stress up to 20% strain, followed by clear strain hardening prior to fracture. To elucidate this unique superplasticity, Figure 10 (a)-(d) provide high-resolition nanometric scale snapshots revealing a stress-induced nanostructural phase transition. This nanomechanical performance is fundamentally governed by a spatial reconfiguration where $s{p}^2$- and $s{p}^3$-hybridized covalent bonds reorient and preferentially concentrate along the uniaxial strain axis. Under load, this atomic network behaves like a highly efficient nanostructural truss system. The dynamic, truss-like redistribution of covalent bonds allows the lattice to stretch and absorb energy continuously, successfully sustaining the prolonged stress plateau and enabling the (14,0) TGCNT to withstand uniaxial elongation before ultimate fracture.

To further elucidate the underlying mechanism of these distinct nanomechanical behaviors, the Poisson’s ratio transitions for both the (0,10) and (14,0) TGCNTs were calculated and are discussed below. This elastic parameter provides critical insights into the lateral nanostructural response of the nanotubes as they undergo severe uniaxial stretching along the z-direction. By tracking the changes in the Poisson’s ratio as a function of tensile strain, we can directly correlate the geometric lattice distortions, such as the truss-like deformation in the superplastic (14,0) system versus the rigid behavior in the brittle (0,10) configuration. The transverse deformation response of the tetragraphene-based nanotubes is quantitatively assessed via the evolution of their Poisson’s ratio under tensile loading (see Figure 11). To contextualize the physical characteristics of the investigated systems, a comparative continuum spectrum of reference Poisson’s ratio values is presented below.

Within the small elastic deformation regime (up to an axial strain threshold of approximately 8%), a linear fitting of the transverse-axial strain relation yields an averaged Poisson’s ratio of $\nu =0.077$ for the TGCNTs $(0,10)$. This exceptionally low value confirms that tetragraphene-based nanotubes possess a significantly constrained lateral flexibility compared to conventional graphene sheets under identical tensile conditions. Physically, this response implies that the 1D nanostructure behaves almost independently along its transverse direction during early-stage loading, approaching the performance limit of a zero-transverse-strain material.

First-principles density functional theory (DFT) calculations and atomistic molecular dynamics literature establish that the Poisson’s ratio of tetragraphene nanostructures follows a highly non-monotonic path [64,185,186]. From an initial ground-state value [185], the system undergoes a localized nanostructural expansion. This transient inflation occurs as the buckled three-dimensional network uncoils, stretching the constituent tetragonal and hexagonal carbon rings along the pulling axis. Immediately following this intermediate nanostructural elongation, the morphology triggers a severe nanostructural collapse, plunging into a strongly negative regime known as the auxetic effect [64,186]. Consequently, capturing a Poisson’s ratio of $0.09$ at exactly 10% strain indicates that our reactive atomistic simulations successfully resolved the precise peak of this geometric phase transition right before the activation of the auxetic threshold. So, Tetragraphene-based nanotubes (TGCNTs) exhibit a low average Poisson’s ratio of $\nu =0.077$ within the small elastic deformation regime (up to 10% uniaxial strain), indicating significantly constrained lateral flexibility, with specific values of $0.07$ for TGCNT(0,10) and $1.19$ for TGCNT(14,0) highlighting strong chirality dependence. This behavior stems from the non-monotonic nanostructural evolution of the tetragraphene lattice, transitioning from initial expansion to severe collapse, which results in a 10% strain Poisson’s ratio of $0.09$ for TGCNT(0,10), marking the peak of this geometric phase transition prior to the auxetic regime.

In summary, the evaluate the nanomechanical boundaries of these systems, it is essential to contrast the benchmark values derived from Density Functional Theory (DFT) with the data obtained via Classical Molecular Dynamics (CMD) simulations in this work:

• DFT Monolayer Benchmarks: First-principles DFT calculations in literature establish the baseline ground-state Poisson’s ratio for an isolated tetragraphene monolayer at a highly constrained value of $\nu =0.072$ [185].
• CMD TGCNTs Results: Our reactive atomistic CMD simulations reveal an outstanding chirality dependence when rolling the 2D sheet into 1D nanotubes. Within the early elastic regime, the CMD data yields a localized value of $\nu =0.07$ for the zigzag-like TGCNT $(0,10)$, showing remarkable agreement with DFT monolayer predictions. Conversely, the CMD calculations reveal an anomalously high value of $\nu =1.19$ for the TGCNT $(14,0)$ configuration.

This severe disparity between the two nanotube chiralities resolved by our CMD model highlights a robust nanostructural anisotropy driven entirely by the rolling orientation of the underlying tetragraphene matrix. Because the rectangular unit cell of tetragraphene features highly asymmetric tetragonal and hexagonal carbon ring arrangements, uniaxial loading along the armchair-like or zigzag-like directions triggers completely different deformation mechanisms. For the $(14,0)$ TGCNT, the uniaxial pulling causes a rapid, compliant closure of the lateral lattice gaps, resulting in an enhanced transverse contraction that vastly exceeds the traditional theoretical upper limit of isotropic continuum materials ($\nu =0.50$).

As the applied uniaxial deformation increases to a larger strain level of 10%, the calculated Poisson’s ratio for the TGCNT $(0,10)$ exhibits a distinct non-linear shift, rising slightly to a value of $0.09$. This non-linear increment is physically valid and reveals a vital nanomechanical signature unique to the underlying tetragraphene topology under large deformations. Both DFT and CMD literature establish that the Poisson’s ratio of tetragraphene nanostructures follows a highly non-monotonic path [64,185,186]. From its initial equilibrium ground state, the system undergoes a localized nanostructural expansion. This transient inflation occurs as the buckled three-dimensional network uncoils, stretching the constituent tetragonal and hexagonal carbon rings along the pulling axis. Immediately following this intermediate nanostructural elongation, the morphology triggers a severe nanostructural collapse, plunging into a strongly negative regime known as the auxetic effect [64,186]. Consequently, capturing a Poisson’s ratio of $0.09$ at exactly 10% strain indicates that our reactive CMD simulations successfully resolved the precise peak of this geometric phase transition right before the activation of the auxetic threshold.

4.1. Validation and Rationale for the AIREBO-Morse Interatomic Potential Selection

To evaluate the nanomechanical behavior and thermal stability of the proposed tetragraphene-based nanotubes (TGCNTs), the selection of an appropriate interatomic potential is a critical parameter. While reactive force fields such as ReaxFF are widely implemented for simulating complex chemical reactions and large-scale nanostructural transformations [187], they often fail to accurately reproduce the nanomechanical responses of pure carbon allotropes that exhibit non-hexagonal or highly strained geometries, such as tetragraphene $1D$.

Preliminary testing revealed that ReaxFF is inadequate for capturing the statistical averages of the nanomechanical properties of TGCNTs. This limitation stems from two primary factors. First, standard ReaxFF parameterizations for carbon are heavily optimized for classic phase transitions and bond-breaking events (e.g., graphite-to-diamond transitions or hydrocarbon combustion) [187]. Consequently, they rarely account for the high steric strain intrinsic to the four-membered rings (squares) that comprise the tetragraphene lattice. Second, ReaxFF utilizes a continuous bond-order calculation without a rigid cutoff system [187]. Under severe uniaxial strain, this unconstrained approach leads to premature amorphization or unrealistically overestimates the flexibility of the distorted bond angles within the TGCNT network.

Conversely, the Adaptive Intermolecular Reactive Empirical Bond Order (AIREBO) framework is well-established as a premier choice for modeling elastic properties and fracture mechanisms in pristine carbon architectures. In this study, we utilize the AIREBO-Morse modification to mitigate the well-known artifact of the original AIREBO formulation, namely, the non-physical strain hardening that occurs near the rupture threshold due to an overly abrupt cutoff function. By replacing the traditional Lennard-Jones term with a long-range Morse potential, the AIREBO-Morse model achieves superb accuracy in describing the attractive and repulsive energy profiles of carbon-carbon bonds under extreme stretching [191]. This modification effectively preserves the correct lattice stiffness up to the precise moment of localized fracture. Furthermore, benchmarking the interatomic potential against first-principles calculations represents the gold standard in molecular dynamics methodologies. As demonstrated in Appendix A and Appendix B, our AIREBO-Morse implementation yields excellent statistical agreement with Density Functional Theory (DFT) data, fully justifying its selection for predicting the nanomechanical pathways of the studied TGCNTs.

5. Conclusions

In this study, comprehensive reactive classical molecular dynamics (CMD) simulations utilizing the calibrated AIREBO-Morse potential successfully mapped the systematic nanomechanical behavior and thermal stability thresholds of single-walled tetragraphene-based nanotubes (TGCNTs). Our primary novel discovery reveals a sharp, chirality-dependent divergence in failure mechanisms under uniaxial tensile loading. The $(0,N)$ TGCNTs consistently exhibit brittle fracture modes characterized by rigid linear elasticity and sudden nanostructural cleavage. In stark contrast, the $(N,0)$ TGCNTs demonstrate exceptional ductility and superelasticity, driven by a stress-induced nanostructural phase transition that sustains a prolonged stress plateau up to 20% strain, followed by clear strain hardening until ultimate failure above 40% strain. Furthermore, we identified an extreme thermal boundary at 1500 K, above which severe nanostructural degradation triggers a complete lattice collapse into a fully amorphous graphene phase. Finally, evaluating the transverse elasticity revealed a critical nanostructural anisotropy: the $(0,10)$ configuration possesses an ultra-low Poisson’s ratio ($\nu =0.07$) in excellent agreement with DFT benchmarks, whereas the $(14,0)$ series displays an anomalously high value ($\nu =1.19$) stemming from its dynamic, truss-like geometric reconfiguration.

The impact of these discoveries is highly significant for the design of next-generation carbon nanomaterials. The extraordinarily high Young’s modulus values achieved, reaching up to $3166.20$ GPa·Å for the $(N,0)$ series and $2324.30$ GPa·Å for the $(0,N)$ series position TGCNTs as premier candidates for nanostructural reinforcement in high-performance nanocomposites. Moreover, the unique combination of a long constant-stress plateau and superelastic deformation over 40% strain indicates that $(N,0)$ TGCNTs can absorb and dissipate massive amounts of mechanical energy before failing. This unprecedented nanomechanical resilience makes them exceptionally well-suited for nanoscale shock-absorbing devices and protective coatings. Additionally, their ability to withstand severe lattice distortions while maintaining structural integrity opens up new avenues for implementation in reliable flexible electronics and strain-engineered components operating within moderate thermal environments.

Looking ahead, several crucial perspectives and future research directions emerge from this work to broaden our understanding of tetragraphene systems. First, while this study establishes a clear thermal stability limit at 1500 K, future investigations should focus on mapping the thermal conductivity and heat transport mechanisms of these nanotubes across different chiralities to assess their efficiency in thermal management applications. Second, exploring the effects of chemical doping (such as introducing nitrogen or boron atoms into the four-membered rings) or nanostructural defects (like Stone-Wales vacancies) could provide a reliable strategy to tune both the superplastic threshold and the electrical properties of the lattice. Lastly, evaluating the interfacial shear stress and load-transfer efficiency between TGCNTs and polymer matrices will be essential to accelerate their practical integration into real-world, high-strength engineering nanocomposites.

Appendix A.1. The Fundamental Framework of the AIREBO Potential

The total potential energy E of a hydrocarbon system containing N atoms within the AIREBO framework is defined as the sum of three distinct energetic terms [188]:

$$\begin{array}{c}\hfill E={\displaystyle \frac{1}{2}}{\displaystyle \sum _{i=1}^N}{\displaystyle \sum _{j\ne i}^N}\left[E_{ij}^{\mathrm{REBO}}+E_{ij}^{\mathrm{non}-\mathrm{bonded}}+{\displaystyle \sum _{k\ne i,j}^N}{\displaystyle \sum _{l\ne i,j,k}^N}{E}_{kijl}^{\mathrm{TORSION}}\right],\end{array}$$

(A1)

where: $E_{ij}^{\mathrm{REBO}}$ represents the short-range, reactive covalent bond energy. $E_{ij}^{\mathrm{non}-\mathrm{bonded}}$ corresponds to the long-range, non-covalent interactions. $E_{kijl}^{\mathrm{TORSION}}$ accounts for the four-body dihedral torsion preferences [188].

Appendix A.1.1. Short-Range Reactive Term (E ij REBO )

The covalent interaction mimics the classic Brenner REBO formulation, balancing attractive and repulsive atomic curves modulated by a local many-body bond-order parameter [188]:

$$\begin{array}{c}\hfill E_{ij}^{\mathrm{REBO}}=V_R\left(r_{ij}\right)+b_{ij}{V}_A\left(r_{ij}\right).\end{array}$$

(A2)

The individual functions for the two-body repulsive ($V_R$) and attractive ($V_A$) states are given by [188]:

$$\begin{array}{c}\hfill V_R\left(r_{ij}\right)=f_{ij}^C\left(r_{ij}\right)\left(1+{\displaystyle \frac{Q_{ij}}{r_{ij}}}\right)A_{ij}{e}^{-{\alpha }_{ij}{r}_{ij}}\end{array}$$

(A3)

and

$$\begin{array}{c}\hfill V_A\left(r_{ij}\right)=f_{ij}^C\left(r_{ij}\right){\displaystyle \sum _{n=1}^3}{B}_{ij}^{\left(n\right)}{e}^{-{\beta }_{ij}^{\left(n\right)}{r}_{ij}},\end{array}$$

(A4)

Where $f_{ij}^C\left(r_{ij}\right)$ is a smooth cutoff function confining covalent bonds to the short-range zone [188]:

$$\begin{array}{c}\hfill f_{ij}^C\left(r_{ij}\right)=\left\{\begin{array}{cc}1,\hfill & r_{ij}\le r_{ij}^{\min }\hfill \\ {\displaystyle \frac{1}{2}}\left\{1+cos\left[\pi {\displaystyle \frac{r_{ij}-r_{ij}^{\min }}{r_{ij}^{\max }-r_{ij}^{\min }}}\right]\right\},\hfill & r_{ij}^{\min }<r_{ij}<r_{ij}^{\max }\hfill \\ 0,\hfill & r_{ij}\ge r_{ij}^{\max }\hfill \end{array}\right..\end{array}$$

(A5)

The many-body bond order parameter $b_{ij}$ scales down the attraction based on local coordination numbers and angles [188]:

$$\begin{array}{c}\hfill b_{ij}={\displaystyle \frac{1}{2}}\left[p_{ij}^{\sigma \pi }+p_{ji}^{\sigma \pi }\right]+{\pi }_{ij}^{\mathrm{RC}}+{\pi }_{ij}^{\mathrm{DH}}.\end{array}$$

(A6)

Appendix A.1.2. Dihedral Torsional Formulation (E kijl TORSION )

The four-body torsion term represents the energy profile associated with rotation around a central covalent bond. For a linear sequence of four bonded atoms $k-i-j-l$, the torsion potential is formulated as [188]:

$$E_{kijl}^{\mathrm{TORSION}}=f_{ki}^C\left(r_{ki}\right)f_{ij}^C\left(r_{ij}\right)f_{jl}^C\left(r_{jl}\right)V_{\mathrm{tors}}\left({\omega }_{kijl}\right)$$

(A7)

Where $f^C$ values are the standard short-range covalent cutoff functions defined in Equation (5), ensuring that torsional calculations smoothly vanish if any associated bond breaks. The dihedral angle ${\omega }_{kijl}$ is computed via the plane normals of the atomic triplets [188]:

$$cos\left({\omega }_{kijl}\right)={\displaystyle \frac{({\mathbf{r}}_{ki}\times {\mathbf{r}}_{ij})·({\mathbf{r}}_{ij}\times {\mathbf{r}}_{jl})}{|{\mathbf{r}}_{ki}\times {\mathbf{r}}_{ij}\left|\right|{\mathbf{r}}_{ij}\times {\mathbf{r}}_{jl}|}}.$$

(A8)

The discrete torsional energy function $V_{\mathrm{tors}}$ uses a standard cosine expansion modulated by the local hybridization state of the core atoms i and j [188]:

$$V_{\mathrm{tors}}\left({\omega }_{kijl}\right)=\left[1-P_{ij}({\mathrm{N}}_i,{\mathrm{N}}_j)\right]\left[{\displaystyle \frac{1}{2}}{V}_1(1+cos{\omega }_{kijl})+{\displaystyle \frac{1}{2}}{V}_2(1-cos2{\omega }_{kijl})+{\displaystyle \frac{1}{2}}{V}_3(1+cos3{\omega }_{kijl})\right],$$

(A9)

where:$V_1,V_2,V_3$ are constant energy coefficients capturing radical, double, and triple bond symmetry bars.$P_{ij}({\mathrm{N}}_i,{\mathrm{N}}_j)$ is an adaptive polynomial scaling function that depends on the coordination numbers N of atoms i and j. This eliminates torsional resistance when either core node approaches an $s{p}^3$ configuration, maintaining physical flexibility in systems like diamond lattices or saturated alkanes [188].

Appendix A.2. The Classical Non-Bonded Failure: Lennard-Jones (LJ) Intermolecular Potential

In the standard AIREBO formulation, long-range dispersion and short-range steric repulsions between non-bonded atoms are governed by an adaptive Lennard-Jones 12-6 potential [189]:

$$\begin{array}{c}\hfill E_{ij}^{\mathrm{non}-\mathrm{bonded}}=E_{ij}^{\mathrm{LJ}}=S\left(t_r\right)S\left(t_b\right)C_{ij}{V}_{\mathrm{LJ}}\left(r_{ij}\right)+[1-S\left(t_r\right)]C_{ij}{V}_{\mathrm{LJ}}\left(r_{ij}\right).\end{array}$$

(A10)

The isotropic unscaled functional form stands as [189]:

$$\begin{array}{c}\hfill V_{\mathrm{LJ}}\left(r_{ij}\right)=4{ϵ}_{ij}\left[{\left({\displaystyle \frac{{\sigma }_{ij}}{r_{ij}}}\right)}^{12}-{\left({\displaystyle \frac{{\sigma }_{ij}}{r_{ij}}}\right)}^6\right].\end{array}$$

(A11)

The Steric Divergence Bottleneck.

Under high-pressure conditions ($\ge 10$ GPa), extreme compression forces non-bonded atoms past their equilibrium radii ($r_{ij}\ll {\sigma }_{ij}$). As $r_{ij}\to 0$, the repulsive force behaves asymptotically [188]:

$$F_{\mathrm{repulsive}}=-{\displaystyle \frac{\partial V_{\mathrm{LJ}}}{\partial r_{ij}}}\propto {\displaystyle \frac{12·4{ϵ}_{ij}{\sigma }_{ij}^{12}}{r_{ij}^{13}}}\to \infty .$$

(A12)

This power-law divergence introduces unphysical, catastrophic core rigidity. It heavily overestimates the stiffness of molecular arrangements (e.g., compressed carbon sheets and polyethylene) [188].

Appendix A.3. The AIREBO-Morse Formulation (AIREBO-M Correction)

To resolve the unphysical high-pressure divergence, the AIREBO-Morse potential replaces the singular $V_{\mathrm{LJ}}\left(r_{ij}\right)$ function with a smoother, softer exponential Morse potential [188]:

$$E_{ij}^{\mathrm{non}-\mathrm{bonded}}⟶E_{ij}^{\mathrm{Morse}}$$

(A13)

The core mathematical mapping substitutes $V_{\mathrm{LJ}}\left(r_{ij}\right)$ with $V_{\mathrm{M}}\left(r_{ij}\right)$ [188]:

$$V_{\mathrm{M}}\left(r_{ij}\right)={ϵ}_{ij}\left[e^{-2{\alpha }_{ij}(r_{ij}-r_{ij}^{\mathrm{eq}})}-2e^{-{\alpha }_{ij}(r_{ij}-r_{ij}^{\mathrm{eq}})}\right],$$

(A14)

where: ${ϵ}_{ij}$ is the potential well depth (identical to the LJ well depth to retain ambient thermodynamics). $r_{ij}^{\mathrm{eq}}$ is the equilibrium distance parameter. ${\alpha }_{ij}$ regulates the width and stiffness of the potential curve [188].

Evaluating the second derivative of the Morse potential at $r_{ij}=r_{ij}^{\mathrm{eq}}$ [190]:

$${\displaystyle \frac{\partial V_{\mathrm{M}}}{\partial r_{ij}}}={ϵ}_{ij}\left[-2{\alpha }_{ij}{e}^{-2{\alpha }_{ij}(r_{ij}-r_{ij}^{\mathrm{eq}})}+2{\alpha }_{ij}{e}^{-{\alpha }_{ij}(r_{ij}-r_{ij}^{\mathrm{eq}})}\right]$$

(A21)

$${\displaystyle \frac{{\partial }^2{V}_{\mathrm{M}}}{\partial r_{ij}^2}}={ϵ}_{ij}\left[4{\alpha }_{ij}^2{e}^{-2{\alpha }_{ij}(r_{ij}-r_{ij}^{\mathrm{eq}})}-2{\alpha }_{ij}^2{e}^{-{\alpha }_{ij}(r_{ij}-r_{ij}^{\mathrm{eq}})}\right]$$

(A22)

$$\left.{\displaystyle \frac{{\partial }^2{V}_{\mathrm{M}}}{\partial r_{ij}^2}}\right|r_{ij}=r_{ij}^{\mathrm{eq}}={ϵ}_{ij}\left[4{\alpha }_{ij}^2\left(1\right)-2{\alpha }_{ij}^2\left(1\right)\right]=2{ϵ}_{ij}{\alpha }_{ij}^2$$

(A23)

$$2{ϵ}_{ij}{\alpha }_{ij}^2={\displaystyle \frac{72{ϵ}_{ij}}{2^{2/3}{\sigma }_{ij}^2}}\Rightarrow {\alpha }_{ij}^2={\displaystyle \frac{36}{2^{2/3}{\sigma }_{ij}^2}}.$$

(A24)

Extracting the square root yields the standard analytical mapping value [190]:

$${\alpha }_{ij}={\displaystyle \frac{6}{2^{1/3}{\sigma }_{ij}}}={\displaystyle \frac{6·2^{2/3}}{2·{\sigma }_{ij}}}=3\times 2^{2/3}{\displaystyle \frac{1}{{\sigma }_{ij}}}\approx {\displaystyle \frac{4.7622}{{\sigma }_{ij}}}$$

(A25)

Appendix A.3.2. Spline Intermolecular Cutoffs and Softness Fine-Tuning

Appendix A.3.3

The Morse potential is smoothly truncated with a third-order spline over adaptive boundaries [188]:

$$V_{\mathrm{M}}^{\mathrm{truncated}}\left(r_{ij}\right)=V_{\mathrm{M}}\left(r_{ij}\right)·f_{\mathrm{spline}}\left(r_{ij}\right).$$

(A26)

The parameter ${\alpha }_{ij}$ is intentionally uncoupled from the standard mapping value at short distances to optimize the potential for high-pressure limits ( 40 GPa) without shifting the equilibrium lattice structures [188]:

$${\alpha }_{ij}^{\mathrm{optimized}}={\gamma }_{ij}·{\alpha }_{ij}.$$

(A27)

Where ${\gamma }_{ij}$ is a dimensionless scaling parameter optimized using quantum chemical calculations (e.g., MP2/dft) and experimental graphite layer compression curves. This modification prevents unphysical core divergence, allowing the repulsive energy to remain stable even under high atomic densities [188].

Appendix B.1. Validation and Comparative Benchmarking in Nanostructured Carbon Systems

The accuracy and computational transferability of the Adaptive Intermolecular Reactive Empirical Bond Order potential with Morse correction (AIREBO-M) for carbon-based nanostructures, such as carbon nanotubes (CNTs) and graphene sheets, have been extensively validated against both experimental observations and first-principles calculations. While the standard AIREBO potential suffers from unphysical hardening and premature fracture artifacts under severe strain or high pressures ($\ge 10\mathrm{GPa}$) due to the rigid $r^{-12}$ Lennard-Jones repulsion core, the Morse non-bonded correction yields excellent agreement with true experimental equation-of-state parameters.

To ground its operational capability within nanomaterial modeling, Table A2 evaluates AIREBO-M against other prominent computational force fields (ReaxFF, AMBER, COMB3, and DFT) across critical structural, mechanical, and chemical dimensions.

Table A2. Comparative benchmarking matrix of interatomic potentials and electronic methods for carbon nanostructures (Graphene, CNTs, and Nanoparticles).

Method	Primary Carbon Domains	Key Advantage / Validation Target	Core Limitation / Failure Mode	Computational Cost	Key References
AIREBO-M	Graphene, CNTs, Diamond, Hydrocarbons	Prevents unphysical steric divergence under high pressure; accurate elastic moduli and layer spacing.	Lacks dynamic variable charge redistribution.	Low/Mod.	[191,192]
ReaxFF	Chemical reactions, Oxidation, Shockwaves	Simulates explicit bond breaking/forming with dynamic charge transfer.	Overestimates Poisson’s ratio; highly sensitive to small integration timesteps.	High	[193,194,195]
AMBER	Biomolecules, Carbon-Water interfaces	Extremely fast; excellent for long-term equilibrium and hydration layers.	Non-reactive; cannot model bond breaking or phase transitions.	Very Low	[196,197,198]
COMB3	Heterogeneous surfaces, Metal-Oxides	Handles dynamic electronic polarization and variable ionic charges at surfaces.	High parameter complexity; less accurate for pure pristine $s{p}^2$ carbon lattices.	High	[199,200]
DFT	Quantum electronic structures, Defect states	Absolute electronic precision; exact ground state energies and quantum effects.	Strictly limited to small systems (few hundred atoms) and short timescales.	Ext. High	[201,202]

Table A2. Comparative benchmarking matrix of interatomic potentials and electronic methods for carbon nanostructures (Graphene, CNTs, and Nanoparticles).

Method	Primary Carbon Domains	Key Advantage / Validation Target	Core Limitation / Failure Mode	Computational Cost	Key References
AIREBO-M	Graphene, CNTs, Diamond, Hydrocarbons	Prevents unphysical steric divergence under high pressure; accurate elastic moduli and layer spacing.	Lacks dynamic variable charge redistribution.	Low/Mod.	[191,192]
ReaxFF	Chemical reactions, Oxidation, Shockwaves	Simulates explicit bond breaking/forming with dynamic charge transfer.	Overestimates Poisson’s ratio; highly sensitive to small integration timesteps.	High	[193,194,195]
AMBER	Biomolecules, Carbon-Water interfaces	Extremely fast; excellent for long-term equilibrium and hydration layers.	Non-reactive; cannot model bond breaking or phase transitions.	Very Low	[196,197,198]
COMB3	Heterogeneous surfaces, Metal-Oxides	Handles dynamic electronic polarization and variable ionic charges at surfaces.	High parameter complexity; less accurate for pure pristine $s{p}^2$ carbon lattices.	High	[199,200]
DFT	Quantum electronic structures, Defect states	Absolute electronic precision; exact ground state energies and quantum effects.	Strictly limited to small systems (few hundred atoms) and short timescales.	Ext. High	[201,202]

Appendix B.2. Quantitative Validation and Computational Efficiency of AIREBO-M

When modeling large-scale carbon nanostructures like graphene and carbon nanotubes (CNTs), selecting an appropriate interatomic potential requires balancing chemical accuracy with structural sizes. First-principles methods like Density Functional Theory (DFT) offer excellent accuracy but are limited to systems with a few hundred atoms due to their $O\left(N^3\right)$ computational scaling [201,202]. Reactive force fields like ReaxFF allow for dynamic bond breaking but require small timesteps ($\Delta t\le 0.1\mathrm{fs}$) and intensive charge calculation overhead, making them computationally expensive [194].

The AIREBO-Morse (AIREBO-M) potential provides an optimized balance for large-scale molecular dynamics (MD) simulations. By replacing the stiff Lennard-Jones core with a smooth Morse function, AIREBO-M avoids unphysical structural hardening under high mechanical strains or pressure fields without sacrificing the high computational speeds typical of empirical potentials [192]. It maintains a standard MD timestep ($\Delta t=0.5--1.0\mathrm{fs}$) and scales linearly ($O\left(N\right)$) with the number of atoms. This allows for the simulation of millions of carbon nodes over long nanosecond timescales. As shown in Table A3 (quantitative nanotructural and mechanical property matrix), AIREBO-M closely reproduces the fundamental structural and mechanical parameters of graphene and CNTs established by experimental work and quantum mechanical calculations. It delivers near-DFT precision for structural deformation while running orders of magnitude faster than ReaxFF.

Table A3. Quantitative comparison of calculated mechanical and structural properties for graphene and pristine single-walled carbon nanotubes (CNTs) across different computational methods.

System & Property	Experimental / DFT	AIREBO-M	Standard AIREBO	ReaxFF
Graphene
C−C Bond Length (Å)	1.42 [201]	1.42 [192]	1.42 [191]	1.45 [195]
Young’s Modulus (TPa)	1.05 ± 0.05 [203]	1.01 [192]	0.95 [191]	0.82 [195]
Poisson’s Ratio ($\nu$)	0.160 [201]	0.165 [192]	0.210 [191]	0.320 [195]
CNT (10,10)
Tube Diameter (Å)	13.56 [201]	13.59 [192]	13.54 [191]	13.78 [194]
Young’s Modulus (TPa)	0.95 ± 0.10 [204]	0.94 [192]	0.91 [191]	0.79 [194]

Table A3. Quantitative comparison of calculated mechanical and structural properties for graphene and pristine single-walled carbon nanotubes (CNTs) across different computational methods.

System & Property	Experimental / DFT	AIREBO-M	Standard AIREBO	ReaxFF
Graphene
C−C Bond Length (Å)	1.42 [201]	1.42 [192]	1.42 [191]	1.45 [195]
Young’s Modulus (TPa)	1.05 ± 0.05 [203]	1.01 [192]	0.95 [191]	0.82 [195]
Poisson’s Ratio ($\nu$)	0.160 [201]	0.165 [192]	0.210 [191]	0.320 [195]
CNT (10,10)
Tube Diameter (Å)	13.56 [201]	13.59 [192]	13.54 [191]	13.78 [194]
Young’s Modulus (TPa)	0.95 ± 0.10 [204]	0.94 [192]	0.91 [191]	0.79 [194]