A new computational approach for evaluating bending rigidity of graphene sheets incorporating disclinations

Abstract

Two-dimensional (2D) materials exhibit remarkable flexibility and can be transformed into various shapes. Graphene sheets (GSs), in particular, can form conical or saddle-like shapes through the introduction of lattice defects known as disclinations, represented by 5- and 7-membered rings, respectively. These rotational-type lattice defects possess relatively large spontaneous curvature and significantly affect the bending rigidity of the GS. Despite increasing interests in exploiting such deformations for material design, evaluating the bending rigidity of GSs with lattice defects remains challenging owing to the complexity introduced by curvature and defect configurations. In this study, we propose a novel computational method that integrates the Helfrich theory of membranes with molecular-dynamics simulations to analyze the effect of curvature and defect patterns on the bending rigidity of GSs. This hybrid approach enables the direct evaluation of bending rigidity from atomic and geometric structures, eliminating the need for experimental bending tests. Using this method, we reveal, for the first time, contrasting trends in bending rigidity between GSs with monopole and dipole disclinations. In the presence of disclination monopoles, the bending rigidity remains independent of the specific structural pattern. Conversely, disclination dipoles, comprising both conical and saddle-shaped surfaces, induce local shape distortions that lead to localized variations in bending rigidity. These findings provide important guidelines for the design of 2D materials with specific bending rigidities, supporting the development of new materials.

---

1. Introduction

Two-dimensional (2D) materials exhibit remarkable flexibility and can be folded into a variety of shapes as in origami.¹ Moreover, lattice defects are known to induce out-of-plane deformation in 2D materials.² For example, introducing 5- or 7-membered rings into the hexagonal lattice of a graphene sheet (GS) can produce conical or saddle-shaped shapes, respectively (Fig. 1). These structures arise from the removal or insertion of equilateral triangles in a GS, forming rotational-type lattice defects termed disclinations (Fig. 1(a)).³ An extracted triangle, resulting in a 5-membered ring, constitutes a positive disclination monopole, while an inserted triangle, forming a 7-membered ring, corresponds to a negative disclination monopole. In both cases, the magnitude of the Frank vector is π/3. When 5- and 7-membered rings are positioned adjacent to each other, the resulting configuration is referred to as a connected disclination dipole, or more generally, a dislocation monopole, which is a type of line defect in a crystalline lattice; this atomic mismatch is characterized by a Burgers vector b.⁴ Although such disclinations induce local curvature and may affect mechanical properties,⁵ a reliable method for calculating the bending rigidity of GS with lattice defects has yet to be established. This is primarily due to the complexity introduced by curvature and defect patterns. Therefore, the aim of this study is to establish a robust and efficient computational method for quantifying the bending rigidity of GSs with disclinations.

Fig. 1 (a) Creation of a kirigami model of GS with disclinations and analytical model, (b) carbon nanohorn, (c) hybrid material of GS and carbon nanotube, and (d) egg-tray graphene.

In recent years, GS-based materials with various shapes have been developed by utilizing the chirality-altering functionality and curvature-stabilizing effects of disclinations. For example, carbon nanohorns⁶ are a combination of fullerene and conical graphene and contain a 5-membered ring near the apex (Fig. 1(b)). These nanohorns have been experimentally observed and considered promising for applications such as chemical sensing⁷ and molecular delivery.⁸ In hybrid materials of GS and carbon nanotubes (CNTs),⁹ 7-membered rings are observed at their linkages (Fig. 1(c)). Wave-shaped GSs, known as egg-tray graphene, feature periodic arrangements of 5- and 7-membered rings (Fig. 1(d)),¹⁰ and have attracted attention as impact-resistant materials as their resistance to impact is strongly influenced by the local curved-surface structure.¹¹ Furthermore, this wave-shaped GS has a negative Poisson's ratio.¹² Conversely, GSs with 7-membered rings can form saddle-shaped surfaces with negative Gaussian curvature, rendering them promising candidates for nanospring applications, where the geometry acts like a strong elastic spring.¹³ Thus, disclinations have a significant impact on the shape and mechanical properties of materials and play an important role in material design. Changes in chirality and spontaneous curvature in 2D materials affect bending rigidity.^14,15 In experimental studies, the bending rigidity, of ideal GSs without disclinations has been measured using atomic-force microscopy (AFM)¹⁶ or by fabricating GSs on atom-sized steps¹⁷ with a range of 1.2–4.13 eV. In simulation studies, the bending rigidity has been assessed using molecular-dynamics (MD) methods,^18–25 density-functional theory (DFT),^26–28 and continuum theory.²⁹ The MD method has been used in the range 0.7–2.27 eV and the DFT in the range 1.4–1.49 eV. These studies have provided valuable insights into the intrinsic properties of flat GS. However, the bending rigidity of GS with lattice defects remains poorly understood owing to the substantial variability introduced by local curvature, defect patterns, thickness, and defect density, which complicate systematic evaluation. Liu et al.³⁰ analyzed the bending rigidity of planar graphene allotropes containing 5- and 7-membered rings and showed that the bending rigidity decreased with increasing defect density. Wang et al.³¹ also analyzed the bending rigidity of egg-tray graphene and showed that the bending rigidity increases linearly with thickness when the effects of curvature and structural pattern are ignored. Nevertheless, these studies provided limited insights into the effect of local curvature and structural patterns on bending rigidity. Furthermore, conducting bending tests that directly measure the bending rigidity of GS with disclinations is technically difficult because the inherent curvature of the GS reduces its accuracy. In addition, theoretical estimates based on continuum plate theory³² give bending rigidity values of approximately 20 eV, significantly higher than experimental and simulation-based results.

To address these issues, we develop a hybrid approach that unifies the Helfrich theory of membranes^33,34 and the MD method. Although originally developed for fluid membranes such as lipid bilayers, the Helfrich theory can be adapted to covalently bonded 2D materials owing their geometric and mechanical similarities. For example, the formation of fullerenes, a type of nanocarbon material, is similar to the vesicle formation process of lipid membranes.³⁵ By applying energy minimization using the MD method, we construct atomic-scale models and quantify the bending response induced by lattice defects. To validate the proposed method, we evaluate the bending rigidity of GSs with disclinations as a representative case. The proposed method enables direct evaluation of bending rigidity from lattice defect-induced out-of-plane deformations, eliminating the need for bending tests.

The results of this study provide a generalizable and robust approach applicable to a wide range of 2D materials, offering fundamental insights for the design of curvature-engineered functional nanomaterials.

2. Analytical models

Different types of disclinations induce varying spontaneous curvatures. These, along with the interactions between adjacent disclinations, lead to diverse arrangements and disclination densities, which in turn result in variations in the curvature of graphene sheet (GS) surfaces. Four types of analytical models were built in this study: a positive disclination monopole (Fig. 2(a)), a negative disclination monopole (Fig. 2(b)), a connected disclination dipole (CDD) (Fig. 2(c)), and a separated disclination dipole (SDD) (Fig. 2(d)). In the monopole-defect models (Fig. 2(a) and (b)), the lattice-defect core was placed at the center of the analytical models. Single positive and negative disclinations were placed in the dipole models (Fig. 2(c) and (d)), and the number of six-membered rings (n = 1, 5, 10–50) between them was varied.

Fig. 2 Analytical models using: (a) positive disclination monopole, (b) negative disclination monopole, (c) connected disclination dipole, (d) separated disclination dipole.

After constructing the atomistic model with reference to the paper-cut model, we conducted structural optimization of the analytical model was performed via MD simulations using the open-source, large-scale atomic/molecular massively parallel simulator (LAMMPS) (version 23Jun2022).^36,37 The optimization was performed using the conjugate-gradient method, and post structural optimization, we computed the potential energy and in-plane stress of optimized structure.

The interactions between carbon atoms were described by adaptive intermolecular reactive bond order (AIREBO) potentials,²⁰ which can describe covalent and non-covalent bonds between carbon atoms, as well as dihedral angles, and both in-plane and out-of-plane deformations of the GS caused by disclinations. The details of AIREBO are provided in ESI Section 1.1 eqn (S1) and (S2).† The cutoff distance was 1.70 Å; the simulation time step was 1 fs; temperature was set to zero, corresponding to a minimum energy state with no thermal fluctuations (i.e., absolute zero). Furthermore, the total number of atoms was set to 120 000–125 000 atoms, for all models, to minimize the effect of size differences on the bending rigidity.³⁸ After structural optimization, the atomic structures were visualized using the Open Visualization Tool (OVITO).³⁹

3. Methods

Disclinations cause large deformations in GSs; hence, the use of thin-membrane theory is appropriate to account for out-of-plane deformations. In general, the bending rigidity of a rigid plate is calculated using κ_b = Et³/12(1 − ν²) where E is the Young's modulus, ν is Poisson's ratio, and the plate thickness is t.³² However, the physical justification for the application of plate theory is not guaranteed in the case of GS, which cannot be regarded as uniformly elastic. Therefore, we use the total strain energy U_total of the GS as expressed by the following equation:⁴⁰

U_total = U_s + U_b + U_edge (1.1)

U_edge = (const.) (1.4)

The first term represents the stretching energy U_s due to in-plane strain. The second term represents the bending energy U_b due to the out-of-plane deformation based on the Helfrich theory of membrane bending. U_egde represents the free surface energy, the deformation considered in this study does not change this energy, and thus it can be neglected. Here, we adopt the Einstein summation convention, where repeated indices imply summation, ε_ij(i, j = x, y) are the in-plane strains, and ε_kk = ε_x + ε_y. In the bending-energy term, H is the mean curvature, K is the Gaussian curvature, and κ_G is the Gaussian rigidity. The spontaneous curvature represents the curvature at which the membrane segment has minimal energy when considered locally. In this study, we assume the spontaneous curvature to be zero. The integration was performed over the entire area with respect to dS. The disclination core is treated as a singularity; hence, eqn (1.1) cannot be applied.

The Helfrich theory of membrane bending can effectively describe the out-of-plane deformation of fluid membranes, such as lipid bilayers. This theory is applicable to GSs due to shared characteristics between 2D thin films, such as large out-of-plane deformations,⁴¹ small values of bending rigidity,⁴² and thermal fluctuations.⁴³ For GSs with interatomic bonds, in-plane stress must also be considered; hence, eqn (1) represents a combination of stretching and bending energies. The generalized Föppl–von Kármán equations, which account for both in-plane and out-of-plane deformations, are used to describe the stress fields of GSs with lattice defects.^40,44 However, these equations are notoriously difficult to solve. In this study, we simplify the analysis by calculating in-plane stress through energy optimization using MD simulations and curvature using differential geometry, directly applying these to eqn (1). The method also has the advantage that bending rigidity can be assessed without the need for bending tests.

3.1. Stretching energy

By rewriting the stretching energy U_s in the first term of eqn (1) into an expression for stress σ_ij(i, j = x, y) using Hooke's law, we obtain the following equation:where the Young's modulus E = 962 GPa and Poisson's ratio ν = 0.143 were independently calculated from tensile tests on a flat GS, as described in ESI Section 1.2.† The Young's modulus in Fig. S1† is consistent with values reported in previous studies.²⁴ For the film thickness, the interlayer distance of graphite, t = 3.34 Å, was used.⁴⁰ The presence of disclinations generates in-plane strain, which in turn generates static in-plane stress. This in-plane stress is calculated for each atom during structural optimization of the analytical model by MD simulation (stress distributions are shown in ESI Section 2 Fig. S3†). Notably, no bending tests were carried out. This method accounts for the effects of out-of-plane deformation on the in-plane stresses.

3.2. Bending energy

The mean curvature H and the Gaussian curvature K are calculated using the principal curvatures k₁ and k₂.³⁴

K = k₁k₂ (4)

when the surface is given by z = f(x,y), the mean curvature and Gaussian curvature are as follows:where the following definitions apply.

The function z = f(x,y) is expressed in Monge's form, which is valid in this case because the deformation is not large. For large deformations, the surface has to be considered as a manifold. In such cases, the coordinates can be locally redefined so that the formulation takes the same form of eqn (1). In this study, an approximate surface approach was adopted to calculate the local mean and Gaussian curvature at each atom of the GS (Fig. 3).

Fig. 3 (a) Surface approximation. (b) Relationship between the distance from the defect core and strain energy and bond length, and (c) strain-energy distribution diagram of the GS with the positive disclination monopole.

The calculation procedure is as follows:

1. The nearest-neighbor atoms A, B, and C of the focus atom, I, are selected.

2. The plane connecting these three points are defined as the ABC plane and its perpendicular unit normal vector considered.

3. An orthonormal basis v₁, v₂, v₃ is constructed with respect to the unit normal vector.

4. An orthogonal transformation of the atomic coordinates is performed using the basis.

5. The aforementioned process is repeated for atoms A, B, and C.

Up to 9 atoms, including neighboring atoms, to the quadratic surface equation.

f(x,y) = a₁ + a₂x + a₃y + a₄x² + a₅xy + a₆y² (8)

Then, an approximate surface was created. The coefficients a₁ − a₆ were calculated using the least squares method as in eqn (8).

a₂ = f_x, a₃ = f_y, 2a₄ = f_xx, a₅ = f_xy, 2a₆ = f_yy (9)

Furthermore, in eqn (1), the Gaussian rigidity is the resistance of the material to topological changes in the absence of edges.⁴⁵ For closed surfaces such as fullerenes, the Gaussian curvature term becomes a constant and vanishes, in accordance with the Gaussian–Bonnet theorem.⁴⁶ Strictly speaking, disclinations are topological changes in the material; because the analytical models used in this study have edges, the bending and Gaussian rigidities must be discussed simultaneously. However, simultaneously determining these properties requires solving for two variables, which is a highly challenging task. Although several attempts have been made to derive the Gaussian rigidity, various issues remain unresolved. In this study, the Gaussian curvature is close to zero (ESI Section 2 Fig. S4†), except close to the defects; therefore, the Gaussian curvature term is neglected in the calculation of eqn (1).

3.3. Bending rigidity

The relationship of the strain energy per atom u_total, the stretching energy u_s, and the bending energy u_b, can be obtained using eqn (10),

u_b = u_total − u_s = 2κ_bH²S (10)

where the per-atom strain energy u_total can be calculated as the absolute value of the difference between the potential energy of a GS with disclinations and that of an ideal GS calculated using MD simulations. The distributions of 2κ_bH²S and u_b is shown in ESI Section 2 Fig. S5 and S6.† The area per atom is calculated using the following formula:where l is the length of the bonds between carbon atoms, calculated from an ideal GS to be l = 1.3968 Å. The relative error between a GS with positive disclination and an ideal GS in the region of calculation is only 0.036%, confirming the validity of this assumption (Fig. 3(b)). Furthermore, as nonlinearity occurs when the curvature is large,⁴⁷ cutoff values were set to eliminate the effects of the free surface and disclination core. Specifically, a range R < 230 Å (R is the distance from the center of the analytical model) was cut off at a sufficient distance from the free surface, and the curvature and energy cutoff value 2H²S = 5.34 × 10⁻³, and u_total = 5.14 × 10⁻³ eV determined from CNTs as shown in ESI Section 1.3 Fig. S2 and Table S2† were applied. Linear fitting was performed using the least-squares method.

4. Results and discussion

4.1 Disclination monopole

We first calculated the bending rigidity of a GS with a positive disclination monopole that formed a conical surface. The atomistic-model diagram in Fig. 4(a) shows the GS with a positive disclination monopole, indicating that the disclination core and free surface were removed. The disclination core was truncated to a circular shape, corresponding to the mean curvature and the bending-energy distribution, and its size was approximately R = 10.7 Å. The bending energy and 2H²S show an approximately linear relationship (Fig. 4(a)). By performing a linear fit, we determined the slope of this relationship, which represents the bending rigidity, to be 0.952 eV. The result of the fitting with errors in the least squares method are shown in ESI Section 1 Tables S3 and S4.†

Fig. 4 Results for positive disclination: (a) Relationship between 2H²S and bending energy and (b) effect of size on bending rigidity. Results of negative disclination: (c) Relationship between 2H²S and bending energy. (d) Fitting only atoms of positive Gaussian curvature. (e) Plot of stretching energy with distance from disclination core. (f) Plot of bond length per atom with distance from disclination core. The black line is the ideal GS bond length (l = 1.397 Å).

The bending rigidity of an ideal GS without disclinations is 0.967 eV in the analysis using the MD method with the AIREBO potential;²⁴ this is greater than 0.952 eV indicating that the bending rigidity is lower and more flexible due to the presence of disclinations. In addition, previous studies have reported that for the same aspect ratio, bending rigidity increases with size and converges when a certain value is reached.³⁸ We investigated the effect of size on GS with disclinations. As shown in ESI Section 2 Fig. S7 and S8,† we selected regions with low energy and stable structure to ensure linearity in all analytical models. The results showed that the bending rigidity of the GS with disclinations increases with increasing size and eventually converges (Fig. 4(b)). This is because increasing the size of the structure increases the number of unaffected atoms and enhances the statistical weight of the effective linear region, which in turn stabilizes the overall structure. Moreover, the analytical model employed in this study is sufficiently large and the effect of size is negligible.

Next, bending rigidity was calculated for a GS with a negative disclination monopole, forming a saddle-shaped surface. The disclination core in the negative disclination was shaped like a quatrefoil, which aligns with the observed mean curvature and bending energy distribution. The presence of negative bending energy is caused by an over-calculation of the stretching energy in eqn (10); this region may be inelastic. The bending rigidity for the GS with negative disclination was calculated to be 0.939 eV, which is slightly lower than that of the positive disclination monopole. Furthermore, the plots for the negative disclination varied and exhibited nonlinear behavior (Fig. 4(c)). This is attributed to the saddle-shaped deformation of the negative disclination, which is not axisymmetric like the positive disclination with core deformation and depends not only on the distance from the core, but also on the azimuth. Assessments of surface shape are generally made by Gaussian curvature, with positive disclinations representing spherical or conical surfaces, negative disclinations representing saddle-shaped surfaces, and zero representing developable surfaces such as cylindrical surfaces. Surfaces of negative disclination are saddle shaped (negative Gaussian curvature) when viewed as a whole; however, when local Gaussian curvature is calculated for each atom on the surface, both positive and negative Gaussian curvatures are present (Fig. 4(d)). This phenomenon is not captured by continuum theory and is unique to GSs with atomic-scale structures, where curvature varies with interatomic distances and bond angles. Atoms with negative Gaussian curvature are shown to exhibit higher stretching energy than those with positive Gaussian curvature (Fig. 4(e)). Furthermore, a comparison with the bond lengths in an ideal GS reveals that tensile and compressive deformations are asymmetrically distributed around these atoms (Fig. 4(f)). These factors are thought to be responsible for the variation in the structure. When only atoms with positive Gaussian curvature were extracted for the negative disclination, the nonlinearity disappeared (Fig. 4(d)). This indicates that the cause of the nonlinearity lies in the atoms with negative Gaussian curvature. The bending rigidity was 0.950 eV. These results demonstrate that in disclination monopoles with low disclination density and simple geometry when nonlinearity is eliminated, the bending rigidity is approximately the same, for both positive and negative disclinations, despite differences in structural patterns. The difference between the behavior of bending rigidity in positive and negative disclinations is related to the strength of the nonlinearity. This indicates that bending rigidity is not a property of the disclination itself, but a local property of the GS due to the local deformation by the disclination.

Comparing the bending rigidity values of the present study with those of previous studies (Fig. 5), we conclude that the results of the present study lie within a reasonable range of the bending rigidity values calculated by the MD method in previous studies, indicating the validity of our results (the 2D Young's modulus Y is calculated by Y = Et). These results indicate that the method used in this study can be applied for the calculation of the bending rigidity of GS with disclinations.

Fig. 5 Comparison of the bending rigidity values obtained by this study with those obtained by previous studies. Reference values and calculated values are shown in Table 1. The numbers in square brackets are reference numbers.

Table 1 Comparison of results with prior studies. Y is the 2D Young's modulus and κ_b is the bending rigidity

Authors	Y (N m^−1)	κ b (eV)	Models	Methods	Ref.
Brenner	236	0.83	Ideal GS	MD	18
Yakobson	363	0.85	Ideal GS	MD	19
Stuart	279	1.56	Ideal GS	MD	20
Brenner	243	1.41	Ideal GS	MD	21
Sears	304	1.63	Ideal GS	MD	22
Wang	337	0.80	Ideal GS	MD	23
Lebedeva	407	1.02	Ideal GS	MD	24
Lebedeva	236	0.80	Ideal GS	MD	24
Lebedeva	243	1.40	Ideal GS	MD	24
Lebedeva	277	0.97	Ideal GS	MD	24
Lebedeva	243	1.40	Ideal GS	MD	24
Lebedeva	315	0.73	Ideal GS	MD	24
Lebedeva	280	2.27	Ideal GS	MD	24
Lebedeva	449	1.60	Ideal GS	MD	24
Lebedeva	265	1.36	Ideal GS	MD	24
Yan	235	0.81	Ideal GS	MD	25
Kudin	345	1.49	Ideal GS	DFT	26
Munoz	394	1.48	Ideal GS	DFT	27
Wei	357	1.44	Ideal GS	DFT	28
Atalaya	220	0.80	Ideal GS	Other methods	48
Wu	235	0.80	Ideal GS	Other methods	49
Tu	358	1.62	Ideal GS	Other methods	50
Zhou	377	1.14	Ideal GS	Other methods	51
Liu	298	1.15	R57	DFT	30
Liu	243	1.18	C-57	DFT	30
Liu	309	1.12	H567	DFT	30
Liu	283	1.42	PSI-graphene	DFT	30
Kunihiro	321	0.952	Positive disclination	MD + Helfrich	This work
Kunihiro	321	0.950	Negative disclination	MD + Helfrich	This work

---

4.2. Connected disclination dipole

Surfaces with CDDs combine conical and saddle shapes due to the presence of adjacent positive and negative disclinations, which influence the surrounding atomic structure. Therefore, the shape of the disclination core is also affected by both curvature type. If the CDD is regarded as a dislocation monopole,⁴ at a distance approximately R = α|b| (α = 3–6) are deleted, which is a valid defect core (the magnitude of the Burgers vector |b| is 2.46 Å). The bending rigidity is highly scattered and very weakly correlated with CDD compared with disclination monopole (Fig. 6); there is minimal linear correlation between 2H²S and the bending energy. This suggests that eqn (10) cannot be directly applied to composites with multiple surfaces curvatures. One factor contributing to the nonlinearity is that the CDD contains a negative disclination (a 7-membered ring). When analyzed in the same manner as the negative disclination monopole, the large stretching energy and negative Gaussian curvature emerge as key contributing factors. In the positive disclination monopoles, the stretching energy was below 3.81 × 10⁻⁵ eV in all ranges. Therefore, we first eliminated the nonlinear region by selecting atoms with stretching energy less than 3.81 × 10⁻⁵ eV and positive Gaussian curvature, and then plotted these atoms. However, the variation remained too large to accurately determine the bending rigidity (Fig. 6(c) and (d)). Next, we refind the analysis by separating the regions strongly affected by each of the positive and negative disclinations. Specifically, the midpoint of the line connecting the two centers of disclination core was set at z = 0, and the region was divided into a region above z ≥ 0 and a region below z < 0. This enables a forced division of the region into positive and negative disclination portions.

Fig. 6 Calculation and regionalization of bending rigidity at CDD: (a) and (b) when the disclination core and free surface are removed, and (c) and (d) when Gaussian curvature and stretching energy are limited. The color maps in the scatter and atomistic diagrams represent the height in the z-direction and correspond in the two diagrams.

We divide the structure into two regions based on the z-height (z ≥ 0 and z < 0), as indicated in Fig. 7(a), which shows a spatial partition of the structure presented in Fig. 6(d), color-coded by z-height. Fig. 7(b) and (c) are plots of regions divided by z-height. Evidently, linearity was observed within each region divided by z-height. The bending rigidity calculated in the z ≥ 0 region was 0.963 eV (Fig. 7(b)), which closely matches the result of the disclination monopole. Conversely, the z < 0 region (Fig. 7(c)) was further divided into two regions. This indicates a localized change in bending rigidity, which was calculated to be 1.916 eV for the upper side (blue) and 0.097 eV for the lower side (green) (Fig. 7(c)). In particular, the bending rigidity on the lower side was very low, indicating almost no resistance to bending, despite the presence of curvature. This phenomenon may be attributed to the simultaneous influence of positive and negative disclination, which may locally relax the in- and out-plane stress in this region, resulting in almost no bending energy. These results indicate that in CDD, both positive and negative disclinations affect the local bending rigidity. In addition to stretching energy, another possible cause of variation in CDD is the effect of the Gaussian curvature term, which is neglected in eqn (1). Gaussian rigidity will need to be considered in this context. These findings must be carefully considered when designing GS-base materials with various geometries.

Fig. 7 (a) Atomic diagram when the regions are divided. Fitting after cutoff due to stretching energy: (b) z ≥ 0 region, and (c) z < 0 region.

4.3. Separated disclination dipole

We then investigated the effect of varying the disclination distance by inserting a six-membered ring between the disclinations in CDD. We observed that the z-height of GS increased with the distance between the disclinations (Fig. 8). Specifically, the height of disclinations increased linearly, while the height of the analytical model initially increased rapidly and then more gradually (Fig. 8(d)). This phenomenon may be attributed to the fact that the curvature of the entire analytical model changed along with the change in height. Previous studies have shown that height and bending rigidity increase linearly when curvature and structural pattern are kept constant.³¹ In this study, simultaneous changes in height and curvature were analyzed using the disclination dipole model. For the bending rigidity calculations, atoms with positive Gaussian curvature and in-plane energy u_s < 3.81 × 10⁻⁵ eV were selected to avoid nonlinearity, as was in the CDD analysis. Then the bending rigidity was evaluated by dividing the region according to z-height. Fitting results for all analytical models are shown in the ESI Section 2 Fig. S9.† In SDD, the behavior of bending rigidity was found to vary with the number of six-member rings n between disclinations. Specifically, for the n = 1 (Fig. 9(a)–(c)) and n = 5 models, two bending rigidities were calculated for both z ≥ 0 and z < 0; for the n = 10 (Fig. 9(d)–(f)), two bending rigidities were calculated for z < 0 only; after the n = 20 (Fig. 9(g)–(i)), one bending rigidity was calculated for both regions, and each value was almost identical. When n = 1, the overall fitting was also divided into two regions (Fig. 9(a)). The results showed that after n = 20, there was no need to separate the regions, and linearity was observed throughout the structure. This suggests that the structure becomes more stable as the distance between disclinations increases, a phenomenon explained by the convergence of the SDD energy of the SDD to that of a disclination monopole when the separation is sufficiently large.⁴ Notably, after n = 20, the SDD begins to exhibit behavior similar to that of a disclination monopole. The differences in behavior between monopoles in CDD and SDD may stem from the exclusion of the Gaussian curvature term in eqn (1), which can lead to inaccuracies in energy estimation when the distance between the disclinations is short. Furthermore, the limitations of applying the Helfrich theory to GS with discrete structures, and the coupling between in-plane and out-of-plane deformations that interact between closely spaced disclinations need to be carefully considered. A more stringent analysis would require an energy model that considers the effects of Gaussian rigidity and interactions.

Fig. 8 SDD models and height: number of six-membered rings between the disclinations is (a) n = 1, (b) n = 10, and (c) n = 20. (d) Relationship between the number of six-membered rings between the disclinations and the height in the z-direction (purple: maximum height, green: defect cores height).

Fig. 9 Fitting diagrams of (a), (d) and (g) all regions, (b), (e) and (h) z ≥ 0 regions and (c), (f) and (i) z < 0 regions for SSD with (a–c) n = 1, (d–f) n = 10, and (g–i) n = 20. Atomic figures show the region where the fitting was done.

In addition, a weighted average was calculated for each n using the number of atoms as weights, to observe the overall trend, using the following formula:

where κ_k is the bending rigidity values in each of the regions, N_k is the total number of atoms in each region and m represents the number of regions present. eqn (12) may be broadly applicable to GSs characterized by localized bending rigidities, and is applied in this study to both CDD and SDDs. When the distances between disclinations are minimal, there are two bending rigidities at z ≥ 0 or z < 0 or both (Fig. 10(a) and (b)). This phenomenon is attributed to the interaction of both positive and negative disclinations and the instability of the structure. This instability can be attributed to the rapid increase in z-height that occurs when a single six-membered ring is inserted between disclinations (Fig. 8(d)). Conversely, as the distance between disclinations increases, both the bending rigidity of each region and the weighted-average bending rigidity converges to a constant value of approximately 0.940 eV (black line in Fig. 10). This is because the energy and curvature changes become smaller and more stable as the distance between disclinations increases (Fig. 10(d)). Interestingly, the bending rigidity decreases from large values in the z ≥ 0 region and increases from small values in the z < 0 region, finally converging to identical values in both regions. This confirms that the structure stabilizes with increasing distance between the disclinations. The convergence value (approximately 0.940 eV) is slightly lower than the bending rigidity of the flat GS (0.967 eV) and the disclination monopole (approximately 0.950 eV). This result indicates that the bending rigidity is lower for higher disclination densities. The trend is consistent with previous studies showing that planar GS allotropes with higher defect density have lower bending rigidity than GSs without defects.³⁰ In this study, bending rigidity was evaluated by introducing one or two disclinations in the center of the analytical model. Slight differences in bending rigidity were observed even for models with very small disclination density, suggesting that disclination density is a strongly influential factor on bending rigidity. The results also showed that the bending rigidity decreased as the height of the analytical model increased. This finding contrasts with the linear relationship between height and bending rigidity reported in previous research.³¹ The discrepancy is likely due to the simultaneous variation of both height and curvature in the model used in this study. These results suggest that, for GSs with disclination dipoles where height and curvature change simultaneously, the bending rigidity is more strongly influenced by the disclination density than by the height. To summarize for SDD, when the distances between disclinations are close, the local bending rigidity trends are different due to the instability of the structure. Furthermore, as the distance between disclinations increases, the structure stabilizes, and the bending rigidity converges to a constant value (0.940 eV). The results also suggest that the bending rigidity of the GS is more strongly affected by the disclination density than by the height. This finding suggests that the control of disclination density is important in the design of GS materials.

Fig. 10 Relationship between distance between disclinations and bending rigidity of SDD: (a) z ≥ 0, and (b) z < 0. (c) Overall fitting bending rigidity and weighted-average bending rigidity. (d) Relationship between distance between disclinations and 2H²S or bending energy per atom. Bars represent average values per atom.

5. Conclusions

In this study, we investigated the bending rigidity of GSs with disclinations using a novel approach that combines MD simulations and the Helfrich theory of membranes. The results revealed a distinct trend in bending rigidity associated with disclination monopoles and dipoles. For the disclination monopole, the bending rigidity exhibited slight variations depending on the structural pattern, primarily due to the nonlinearity associated with negative disclination. For the disclination dipole, the combination of conical and saddle-shaped surfaces caused a local shape change; a corresponding local change in bending rigidity was observed. These findings highlight the potential of designing GSs with specific bending rigidity characteristics, such as saddle-shaped nanosprings that combine flexibility and rigidity, or egg-tray graphene with localized impact resistance. However, the limitations of the Helfrich theory of membrane bending, particularly for surfaces with negative Gaussian curvature, highlight the need for further research to clarify the relationship between defect placement, density, and rigidity. Such efforts are expected to provide a deeper understanding of the deformation mechanisms of GS with disclination and to facilitate the design of GS-based materials with diverse mechanical properties.

Abbreviations

GS	Graphene sheet
2D	Two-dimensional
CNTs	Carbon nanotubes
AFM	Atomic-force microscopy
MD	Molecular-dynamics
DFT	Density-functional theory
CDD	Connected disclination dipole
SDD	Separated disclination dipole
LAMMPS	Large-scale atomic/molecular massively parallel simulator
AIREBO	Adaptive intermolecular reactive bond order
OVITO	Open visualization tool

Author contributions

Yushi Kunihiro: conceptualization, methodology, data curation, writing – original draft. Xiao-Wen Lei: conceptualization, methodology, data curation, supervision, project administration, writing – reviewing and editing, funding acquisition. Takashi Uneyama: methodology, writing – reviewing. Toshiyuki Fujii: writing – reviewing.

Conflicts of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability

All relevant data are within the manuscript and the data supporting this article have been included as part of the ESI.†

Acknowledgements

This work was supported by JST PRESTO (Grant Number JPMJPR2199), JSPS KAKENHI (Grant Number JP23K25994).

References

1. Y. Wang, Y. Zhang, R. Gover, J. Yang and Y. Zhang, Carbon, 2023, 207, 67–76.
2. Y. Wang and Z. Liu, Nanoscale, 2018, 10, 6487–6495.
3. M. Rozhkov, A. Kolesnikova, I. Yasnikov and A. Romanov, Low Temp. Phys., 2018, 44, 918–924.
4. S. Tsuchida, X.-W. Lei, Y. Kunihiro and T. Fujii, Mater. Today Commun., 2024, 41, 110687.
5. T. Lu, X.-W. Lei and T. Fujii, Carbon, 2024, 228, 119287.
6. S. Iijima, M. Yudasaka, R. Yamada, S. Bandow, K. Suenaga, F. Kokai and K. Takahashi, Chem. Phys. Lett., 1999, 309, 165–170.
7. J. Zhang, J. Lei, C. Xu, L. Ding and H. Ju, Anal. Chem., 2010, 82(3), 1117–1122.
8. M. Lahiani, J. Chen, F. Irin, A. Puretzky, M. Green and M. Khodakovskaya, Carbon, 2015, 81, 607–619.
9. G. Dimitrakakis, E. Tylianakis and G. Froudakis, Nano Lett., 2008, 8(10), 3166–3170.
10. W. Liu, L. Zhao, E. Zurek, J. Xia, Y. Zheng, H. Lin, J. Liu and M. Miao, npj Comput. Mater., 2019, 5, 71.
11. Y. Tomioka, T. Natsuki, J. Shi and X.-W. Lei, Nanomaterials, 2022, 12(3), 436.
12. H. Qin, Y. Sun, J. Liu, M. Li and Y. Liu, Nanoscale, 2017, 9(12), 4135–4142.
13. M. Tadayon, S. Amini, A. Masic and A. Miserez, Adv. Funct. Mater., 2015, 25(41), 6437–6447.
14. E. Jomehzadeh, M. Afshar, C. Galiotis, X. Shi and N. Pugno, Int. J. Non Linear Mech., 2013, 56, 123–131.
15. Y. Zhu, P. Wang, S. Xiao, S. He, J. Chen, Y. Jiang, Y. Wang, J. He and Y. Gao, Nanoscale, 2018, 10, 21782–21789.
16. M. Ashino, K. Nishioka, K. Hayashi and R. Wiesendanger, Phys. Rev. Lett., 2021, 126, 146101.
17. E. Han, J. Yu, E. Annevelink, J. Son, D. Kang, K. Watanabe, T. Taniguchi, E. Ertekin, P. Huang and A. van der Zande, Nat. Mater., 2020, 19, 305–309.
18. D. Brenner, Phys. Rev. B: Condens. Matter Mater. Phys., 1990, 42, 9458–9471.
19. B. Yakobson, C. Brabec and J. Bernholc, Phys. Rev. Lett., 1996, 76, 2511–2514.
20. S. Stuart, A. Tutein and J. Harrison, J. Chem. Phys., 2000, 112, 6472–6486.
21. D. Brenner, O. Shenderova, J. Harrison, S. Stuart, B. Ni and S. Sinnott, J. Phys.: Condens. Matter, 2002, 14, 783–802.
22. A. Sears and R. Batra, Phys. Rev. B: Condens. Matter Mater. Phys., 2004, 69, 235406.
23. L. Wang, Q. Zheng, J. Liu and Q. Jiang, Phys. Rev. Lett., 2005, 95, 105501.
24. I. Lebedeva, A. Minkin, A. Popov and A. Knizhnik, Phys. E, 2019, 108, 326–338.
25. J. Yan and W. Zhang, J. Sound Vib., 2021, 514, 116464.
26. K. Kudin, G. Scuseria and B. Yakobson, Phys. Rev. B: Condens. Matter Mater. Phys., 2001, 64, 235406.
27. E. Munoz, A. Singh, M. Ribas, E. Penev and B. Yakobson, Diamond Relat. Mater., 2010, 19, 368–373.
28. Y. Wei, B. Wang, J. Wu, R. Yang and M. Dunn, Nano Lett., 2013, 13(1), 26–30.
29. B. Azizi, M. Shariati, S. Souq and M. Hosseini, Appl. Math. Models, 2023, 114, 466–487.
30. L. Liu, R. Zhu and J. Zhao, Appl. Surf. Sci., 2021, 569, 151048.
31. M. Wang, L. Jiao, R. Zhu, Z. Tan, S. Dai and L. Liu, Mol. Model., 2022, 28, 364.
32. S. Timoshenko and J. N. Goodier, in Theory of Elasticity, McGraw-Hill, 1969.
33. W. Helfrich, Z. Naturforsch., C:J. Biosci., 1973, 28, 693–703.
34. S. A. Safran, in Statistical thermodynamics of surfaces, interfaces, and membranes, Perseus Book, 1994.
35. S. Irle, G. Zheng, Z. Wang and K. Morokuma, J. Phys. Chem. B, 2006, 110(30), 14531–14545.
36. A. Thompson, H. Aktulga, R. Berger, D. Bolintineanu, W. Brown, P. Crozier, P. Veld, A. Kohlmeyer, S. Moore, T. Nguyen, R. Shan, M. Stevens, J. Tranchida, C. Trott and S. Plimpton, Comput. Phys. Commun., 2022, 271, 108171.
37. S. Plimpton, J. Comput. Phys., 1995, 117(1), 1–19.
38. Q. Wang, Appl. Phys. Lett., 2010, 374, 1180–1183.
39. A. Stukowski, Modell. Simul. Mater. Sci. Eng., 2010, 18, 015012.
40. T. Zhang, X. Li and H. Gao, J. Mech. Phys. Solids, 2014, 67, 2–13.
41. F. Ahmadpoor and P. Sharma, Extreme Mech. Lett., 2017, 14, 38–43.
42. J. Hernández-Muñoz, F. Bresme, P. Tarazona and E. Chacón, J. Chem. Theory Comput., 2022, 18, 3151–3163.
43. M. Zelisko, F. Ahmadpoor, H. Gao and P. Sharma, Phys. Rev. Lett., 2017, 119, 068002.
44. H. Seung and D. Nelson, Phys. Rev. A, 1988, 38(2), 1005–1018.
45. M. Rueda-Contreras, A. Gallen, J. Romero-Arias, A. Hernandez-Machado and R. Barrio, Sci. Rep., 2021, 11, 9562.
46. A. Gray, E. Abbena and S. Salamon, Modern differential geometry of curves and surfaces with Mathematica, Champman & Hall/CRC, 3rd edn, 2006.
47. K. Girgin, Z. Girgin and M. Aksoylu, Smart Mater. Struct., 2014, 47, 1115–1130.
48. J. Atalaya, A. Isacsson and J. Kinaret, Nano Lett., 2008, 8, 4196–4200.
49. J. Wu, K. Hwang and Y. Huang, J. Mech. Phys. Solids, 2008, 56, 279–292.
50. Z. Tu and Z. Ou-Yang, J. Comput. Theor. Nanosci., 2008, 5, 1192–1192.
51. X. Zhou, J. Zhou and Z. Ou-Yang, Phys. Rev. B: Condens. Matter Mater. Phys., 2000, 62, 13692–13696.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d5nr01102g

---