A Coupled Refined Model of Atomistic and Continuum Parameters of Diatomic Covalent Bonds

1. Introduction

Carbon nanostructures, such as nanotubes and graphene, possess unique mechanical properties — high stiffness, strength, and resistance to deformation. Their behavior at the nanoscale is typically described using molecular dynamics (MD) methods, which allow the explicit consideration of atomic interactions through empirical potentials, such as Morse or Lennard-Jones.

However, the direct application of MD models for describing large-scale systems is associated with significant computational costs, due to the need to simulate a representative volume of material containing a relatively small number of molecules or atoms, as well as the limitation to modeling processes over very short time intervals (from picoseconds to nanoseconds) [1]. This significantly complicates the macroscopic modeling of nanomodified materials and structures derived from them, and in some cases renders such modeling impossible.

In this context, there is a need to develop specialized multiscale models that allow for a consistent transition from atomistic mechanics parameters (bond stiffness under tension, bending, and torsion, potential parameters) to phenomenological parameters of solid mechanics (elastic modulus and Poisson’s ratio, shear modulus, plasticity, etc.). Such a transition requires the development of specialized models in which an individual chemical bond is interpreted as a mechanical element with angular, radial, and torsional stiffness. In this regard, the question arises regarding the correct choice of geometric parameters, in particular the effective diameter of the C–C chemical bond, which allows the reduction of the interatomic interaction Equations to the classical Equations of continuum mechanics.

The main difficulty lies in the fact that the formulas of structural mechanics include the Poisson’s ratio, whereas at the level of an individual chemical bond this parameter is not explicitly defined. As a result, traditional methods for determining the energetically equivalent mechanical properties of an effective rod depend on an arbitrarily chosen value of the Poisson’s ratio ν. Moreover, the direct application of classical beam theory Equations to carbon nanostructures frequently results in nonphysical values of the Poisson’s ratio, exceeding the upper limit of ν = 0.5 permitted for isotropic continuum materials.

Thus, a contradiction arises when directly transferring the molecular dynamics bond stiffnesses of the C–C bond into the classical beam theory Equations. From the well-known relations that energetically link the force constant of molecular mechanics $k_r, { k}_{\theta } и k_{\tau }$ with the physico-mechanical characteristics of a rod with a circular cross-section of diameter d within the straight-normal (Euler–Bernoulli) hypothesis [2] (${\mathrm{k}}_{\mathrm{r}}={\displaystyle \frac{\mathrm{E}\mathsf{\pi }{\mathrm{d}}^2}{4\mathrm{L}}}$, ${\mathrm{k}}_{\mathsf{\theta }}={\displaystyle \frac{\mathrm{E}\mathsf{\pi }{\mathrm{d}}^4}{64\mathrm{L}}}$, and ${\mathrm{k}}_{\mathsf{\tau }}={\displaystyle \frac{\mathrm{G}\mathsf{\pi }{\mathrm{d}}^4}{32\mathrm{L}}}$), it follows that $\mathrm{d}=4\sqrt{{\displaystyle \frac{{\mathrm{k}}_{\mathsf{\theta }}}{{\mathrm{k}}_{\mathrm{r}}}}}$. After substituting the known values ${\mathrm{k}}_{\mathrm{r}}=651.694{\displaystyle \frac{\mathrm{n}\mathrm{N}}{\mathrm{n}\mathrm{m}}},{\mathrm{k}}_{\mathsf{\tau }}=0.278{\displaystyle \frac{\mathrm{n}\mathrm{N}\mathrm{n}\mathrm{m}}{{\mathrm{r}\mathrm{a}\mathrm{d}}^2}}{,\mathrm{k}}_{\mathsf{\theta }}=0.901{\displaystyle \frac{\mathrm{n}\mathrm{N}\mathrm{n}\mathrm{m}}{{\mathrm{r}\mathrm{a}\mathrm{d}}^2}}$ for the C–C covalent bond, we obtain $\mathrm{d}\approx 0.147\mathrm{n}\mathrm{m}$, $\mathrm{E}\approx 5.49\mathrm{T}\mathrm{P}\mathrm{a}$, and $\mathrm{G}\approx 0.871\mathrm{T}\mathrm{P}\mathrm{a}$. At the same time, if the Poisson’s ratio ν is determined from the known relations for an isotropic body, as $\mathsf{\nu }={\displaystyle \frac{\mathrm{E}}{2\mathrm{G}}}-1$, we obtain $\mathsf{\nu }={\displaystyle \frac{5.49}{2\mathrm{*}0.871}}-1\approx 2.15>0.5$, which is unphysical, since it violates the thermodynamic stability conditions for the continuum mechanics model of isotropic bodies (volumetric compression leads to a negative increment energy).

Thus, the use of classical formulas for tension, bending, and torsion of a beam to determine the effective physical parameters of the C–C covalent bond can lead to nonphysical values of the Poisson’s ratio, thereby highlighting the need to address the problem of a consistent and physically justified coupling between atomistic and continuum models.

In this work, we show that there exists a specific range of variation of the effective C–C bond diameter within which the conditions of thermodynamic stability are satisfied when transitioning from molecular-mechanics Equations to phenomenological Equations of solid mechanics. This eliminates the problem of obtaining nonphysical values of the Poisson’s ratio and allows the mean diameter within this range to be considered as a universal bond characteristic, ensuring a consistent transition from the atomistic description of the C–C bond to its continuum analogue.

2. Materials and Methods

2.1. State of the Art in Modeling the Effective Mechanical C–C Bond Properties

The carbon–carbon covalent (C–C) bond is a key element in the structure of carbon materials. Its electron density determines the spatial configuration, which in turn governs the mechanical and physical properties of nanostructures, including graphene, carbon nanotubes, and diamond-like materials. Research aimed at determining the macroscopic characteristics of C–C bonds from molecular mechanics parameters began to develop actively in the late 1980s and received particular attention after the discovery of carbon nanotubes [3]. The main focus was on establishing the relationship between empirical interatomic potentials and effective elastic constants.

To date, several key research directions have emerged in the literature, such as:

• the development and refinement of interatomic interaction potentials [4,5,6,7,8,9];
• the coupling molecular dynamics models with continuum models (beam and shell analogues of bonds) [10,11,12,13,14,15,16];
• the experimental validation (measurements of strength and elasticity of nanotubes and graphene) [4,17,18,19];
• the application of quantum-chemical and density functional theory (DFT) methods for refining C–C bond parameters [20,21,22];
• the development of multiscale theories integrating data from different approaches [23,24,25,26].

The determination of macroscopic mechanical properties of carbon nanostructures, such as carbon nanotubes and graphene, is directly related to the accurate description of the properties and geometric parameters of the C–C covalent bond. To date, one of the key challenges remains the transition from atomistic molecular mechanics parameters (bond stiffnesses, potential functions) to macroscopic continuum mechanics parameters (Young’s modulus, Poisson’s ratio, shear modulus).

However, the literature shows significant discrepancies in the values of the so-called “effective diameter” of the C–C bond, which lead to substantial variations in the calculated effective elastic modulus and Poisson’s ratio of representative C–C bond volumes. Early works [5,6,27], laid the foundation for the theory of interatomic interaction potentials describing atomic interactions.

Based on these studies, the REBO (Reactive Empirical Bond-Order) potentials [7] and their AIREBO [9] modifications were developed, which are widely used in molecular dynamics simulations [8]. These models allow the determination of microscopic bond stiffnesses; however, transitioning to macroscopic characteristics requires the introduction of additional geometric parameters, including the effective bond diameter. In particular, the effective values of Young’s modulus calculated based on REBO and AIREBO potential range from 0.9 to 1.3 TPa [28], while the Poisson’s ratios range from 0.15 to 0.45. At the same time, the use of structural-mechanical analogies sometimes leads to nonphysical values of ν > 0.5 [29,30,31].

Recent studies [30,32] show that when choosing an effective diameter d for the C–C bond in the range of 0.04 nm to 0.10 nm, significant discrepancies can arise in the calculation of the Poisson’s ratio ν. For example, when using d$\ge$0.1 nm, the values of ν can exceed 0.5, which is unphysical from the perspective of structural mechanics, as it violates the condition of thermodynamic stability. With a more precise consideration of molecular-mechanical parameters and a refined value range of d ≈ 0.076 nm – 0.09 nm, the resulting values of ν return to the physically permissible range (0 < ν < 0.5).

Studies addressing this problem have shown that different potentials (REBO, AIREBO, TERSOFF, REAXFF) [4] yield significantly varying values of elastic moduli if the bond diameter correction is not taken into account [23]. Particular attention in the literature is paid to the fact that standard beam and bar theory Equations, when directly applied at the level of a covalent bond, can lead to unphysical results if a diameter value of around 0.1 nm is used without correction. Therefore, the commonly accepted bond diameter d = 0.147 nm should not be used without appropriate corrections [17,30].

Thus, the correct choice of the effective C–C bond diameter is central to achieving a consistent transition from molecular to structural mechanics.

2.2. Selection of calculation schemes for C-C bond

The method proposed in this work for determining the physics-mechanical characteristics of a C–C covalent bond between two carbon atoms consists of the following steps:

• Determination of the effective C–C bond diameter, taking shear deformation into account, using a refined Timoshenko beam model;
• Determination of the effective Poisson’s ratio;
• Determination of the effective Young’s modulus based on the Morse potential;
• Determination of the effective shear modulus of the C–C bond.

Let us consider each of these points separately.

Accounting for shear deformation according to Timoshenko in modeling the C–C covalent bond as a beam model is necessary due to the fundamental difference between the assumptions of classical Euler–Bernoulli bending theory and the physical reality of interatomic interactions at the nanoscale.

Classical Euler–Bernoulli beam theory is based on two key assumptions:

• Plane section hypothesis: It is assumed that the beam’s cross-sections remain plane and perpendicular to the neutral axis during bending.
• Neglect of transverse shear: This model assumes infinite material stiffness with respect to transverse shear ($\mathrm{G}\to \mathrm{\infty }$).

These assumptions work well only for “very long” and “very thin” beams, where the primary contribution to deformation comes from pure bending, and shear deformation is negligibly small. In the case of modeling an individual C–C covalent bond as a beam in a continuum model (for example, to calculate the elastic properties of carbon nanotubes), the C–C bond has a length of approximately 0.142 nm, and the effective “bond diameter,” determined according to the plane section hypothesis (${\mathrm{d}}_{\mathsf{\theta }}\approx 4\sqrt{{\displaystyle \frac{{\mathrm{k}}_{\mathsf{\theta }}}{{\mathrm{k}}_{\mathrm{r}}}}}=0.147\mathrm{n}\mathrm{m}$), is not negligibly small compared to the bond length.

In this context, the C–C bond, according to the classification of beam theories, can be considered a geometrically “short” or “thick” beam. Therefore, the use of Euler–Bernoulli theory in this case leads to an artificially overestimated bond diameter and an underestimated bond elastic modulus as result of compensating for the neglected shear deformation. In addition, when real atomic structures (e.g. molecules or nanotubes) bend or vibrate, the angle between bonds can change significantly, which further necessitates to accounting for transverse shear deformation in the continuum model.

Therefore, in this work, the Timoshenko model is adopted as the fundamental hypothesis for the transition from the atomic model to the continuum, whose advantages can be summarized as follows:

• Timoshenko theory, by introducing the shear angle as an independent parameter, removes the assumption of cross-section perpendicularity to the neutral axis and accounts for both bending deformation and transverse shear deformation independently.
• The total potential deformation energy in the Timoshenko model includes not only the bending energy but also the shear energy. In this case, the Timoshenko model provides a more physically adequate representation of the behavior of continuum bond models when simulating the actual behavior of atoms and bonds, especially when analyzing local bending and high-frequency vibrations (when the wavelengths are comparable to interatomic distances).
• Furthermore, accounting for shear strain plays a significant role in the analysis of high-frequency vibrations (for example, when determining the natural modes and frequencies of vibrations of nanostructures).

In addition, the Timoshenko model allows for a more adequate description of local stiffness variations caused by defects (e.g., Stone–Wales defects) or changes in the orientation of C–C bonds, which is important for a more accurate determination of the strength characteristics of nanostructures.

Thus, accounting for shear according to Timoshenko is not strictly necessary for determining the intrinsic diameter of the C–C bond, but is required for the correct application of continuum models when calculating the mechanical and dynamic properties of systems composed of these bonds (e.g., carbon nanotubes), where classical simplified Euler–Bernoulli theories produce significant errors due to the considerable influence of transverse shear and can lead to unphysical results when obtaining effective physics-mechanical and geometric parameters.

2.3. C-C bond shear effect

The simplest method for accounting for the shear effect was proposed in [33], where a shear coefficient was introduced into the bending Equations of the Timoshenko beam.

There are several methods for determining the shear coefficient $\mathrm{K}$. In general, it should satisfy the following condition:

$${\int }_{\mathrm{S}}{\mathsf{\sigma }}_{\mathsf{\tau }}\mathrm{d}\mathrm{A}=\mathrm{K}\mathrm{A}\mathrm{G}\mathsf{\phi },$$

(1)

where ${\mathsf{\sigma }}_{\mathsf{\tau }}$ is the transverse shear stress acting in the beam, $\mathrm{A}$ is the cross-sectional area of the beam, $\mathrm{G}$ is the shear modulus of the beam material, and $\mathsf{\phi }$ is the rotation angle of the cross-section normal during deformation.

According to [34], the shear coefficient can be expressed by an Equation depending on the Poisson’s ratio and the geometric parameters of the beam’s cross-section as follows:

$$\mathrm{K}={\displaystyle \frac{2\left(1+\mathsf{\nu }\right){\mathrm{I}}_{\mathrm{x}}}{\frac{\mathsf{\nu }({\mathrm{I}}_{\mathrm{y}}-{\mathrm{I}}_{\mathrm{x}})}{2}-\frac{\mathrm{A}}{{\mathrm{I}}_{\mathrm{x}}}\iint \mathrm{x}\left(\mathsf{\chi }+\mathrm{x}{\mathrm{y}}^2\right)\mathrm{d}\mathrm{x}\mathrm{d}\mathrm{y}}},$$

(2)

where $\mathsf{\chi }(\mathrm{x},\mathrm{y},\mathsf{\nu })$ is a harmonic function depending on the coordinates $\mathrm{x},\mathrm{y}$ of the beam’s cross-section and the Poisson’s ratio $\mathsf{\nu }$ [34], ${\mathrm{I}}_{\mathrm{x}}=\iint {\mathrm{y}}^2\mathrm{d}\mathrm{x}\mathrm{d}\mathrm{y}$ is the moment of inertia of the cross-section about the $y$-axis, and ${\mathrm{I}}_{\mathrm{y}}=\iint {\mathrm{x}}^2\mathrm{d}\mathrm{x}\mathrm{d}\mathrm{y}$ is the moment of inertia of the cross-section about the $x$-axis.

For a solid circular cross-section beam, the shear coefficient $\mathrm{K}$ is equal to

$$\mathrm{K}={\displaystyle \frac{6(1+\mathsf{\nu })}{(7+6\mathsf{\nu })}}$$

(3)

and for a solid rectangular cross-section, it is equal to

$$\mathrm{K}={\displaystyle \frac{10(1+\mathsf{\nu })}{12+11\mathsf{\nu }}}$$

(4)

In this case, it does not depend on the ratio of the beam’s width to its height [34].

The total steric potential energy of a diatomic C–C system is associated with the spatial arrangement of the atoms. For small linear elastic deformations, neglecting electrostatic interactions and the energy associated with nonbonded van der Waals forces, it can be expressed as [2,11]:

$${\mathrm{U}}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}={\mathrm{U}}_{\mathrm{r}}+{\mathrm{U}}_{\mathsf{\theta }}+{\mathrm{U}}_{\mathsf{\phi }}+{\mathrm{U}}_{\mathsf{\omega }}={\mathrm{U}}_{\mathrm{r}}+{\mathrm{U}}_{\mathsf{\theta }}+{\mathrm{U}}_{\mathsf{\tau }}$$

(5)

where:

${\mathrm{U}}_{\mathrm{r}}={\displaystyle \frac{1}{2}}{\mathrm{k}}_{\mathrm{r}}{\left(\mathsf{\Delta }\mathrm{r}\right)}^2$ — energy associated with bond tension;

${\mathrm{U}}_{\mathsf{\theta }}={\displaystyle \frac{1}{2}}{\mathrm{k}}_{\mathsf{\theta }}{\left(\mathsf{\Delta }\mathsf{\theta }\right)}^2$ — energy associated with bond angle change during bending;

${\mathrm{U}}_{\mathsf{\tau }}={\mathrm{U}}_{\mathsf{\phi }}+{\mathrm{U}}_{\mathsf{\omega }}={\displaystyle \frac{1}{2}}{\mathrm{k}}_{\mathsf{\tau }}{\left(\mathsf{\Delta }\mathsf{\phi }\right)}^2$ — energy accounting for total torsion under the dihedral angle ${\mathrm{U}}_{\mathsf{\phi }}$ and out-of-plane torsion energy ${\mathrm{U}}_{\mathsf{\omega }}$.

The force constants ${\mathrm{k}}_{\mathrm{r}}$, ${\mathrm{k}}_{\mathsf{\theta }}$ and ${\mathrm{k}}_{\mathsf{\tau }}$ represent the bond’s resistance to tension, bending, and torsion, respectively. $\mathsf{\Delta }\mathrm{r}$, $\mathsf{\Delta }\mathsf{\theta }\mathrm{a}\mathrm{n}\mathrm{d}\mathsf{\Delta }\mathsf{\phi }$ correspond to the increments of longitudinal displacement under tension, angle changes during bending, and bond torsion, respectively. A schematic representation of the interatomic interactions is shown in Figure 1(a).

Following the procedure described in [2], we write the expressions for determining the potential energy of the beam in the cases of tension, bending, and torsion (Figure 1(b)):

$${\mathrm{U}}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}={\mathrm{U}}_{\mathrm{r}}+{\mathrm{U}}_{\mathsf{\theta }}+{\mathrm{U}}_{\mathsf{\phi }}+{\mathrm{U}}_{\mathsf{\omega }}={\mathrm{U}}_{\mathrm{r}}+{\mathrm{U}}_{\mathsf{\theta }}+{\mathrm{U}}_{\mathsf{\tau }}$$

(6)

$${\mathrm{U}}_{\mathrm{M}}={\displaystyle \frac{1}{2}}{\int }_0^{\mathrm{L}}{\displaystyle \frac{{\mathrm{M}}^2}{\mathrm{E}\mathrm{I}}}\mathrm{d}\mathrm{L}={\displaystyle \frac{1}{2}}{\displaystyle \frac{\mathrm{E}\mathrm{I}}{\mathrm{L}}}{\left(2\mathsf{\alpha }\right)}^2-\mathrm{b}\mathrm{e}\mathrm{n}\mathrm{d}\mathrm{i}\mathrm{n}\mathrm{g},$$

(7)

$${\mathrm{U}}_{\mathrm{T}}={\displaystyle \frac{1}{2}}{\int }_0^{\mathrm{L}}{\displaystyle \frac{{\mathrm{T}}^2}{\mathrm{G}\mathrm{J}}}\mathrm{d}\mathrm{L}={\displaystyle \frac{1}{2}}{\displaystyle \frac{\mathrm{G}\mathrm{J}}{\mathrm{L}}}{\left(\mathsf{\Delta }\mathsf{\beta }\right)}^2-\mathrm{t}\mathrm{o}\mathrm{r}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n},$$

(8)

where $\mathrm{E}$ — Young’s modulus; $\mathrm{G}$ — shear modulus; $\mathrm{A}$ — cross-sectional area; $\mathrm{L}$ — beam length; $\mathrm{I}={\displaystyle \frac{\mathsf{\pi }{\mathrm{d}}^4}{64}}$ — moment of inertia; $\mathrm{J}={\displaystyle \frac{\mathsf{\pi }{\mathrm{d}}^4}{32}}$ — polar moment of inertia of a circular cross-section beam; ∆L, α, and ∆β — increment of length, bending angle, and increment of torsion angle, respectively; $\mathrm{d}$ — effective bond beam diameter, nm.

By comparing the molecular mechanics Equations with Equations (6)–(8) of structural mechanics, we obtain the well-known relations [2]:

$${\displaystyle \frac{\mathrm{E}\mathrm{A}}{\mathrm{L}}}={\mathrm{k}}_{\mathrm{r}},{\displaystyle \frac{\mathrm{E}\mathrm{I}}{\mathrm{L}}}={\mathrm{k}}_{\mathsf{\theta }},{\displaystyle \frac{\mathrm{G}\mathrm{J}}{\mathrm{L}}}={\mathrm{k}}_{\mathsf{\tau }}$$

(9)

from which it follows that

$$\mathrm{E}={\displaystyle \frac{4{\mathrm{k}}_{\mathrm{r}}\mathrm{L}}{\mathsf{\pi }{\mathrm{d}}^2}}\mathrm{and} \mathrm{G}={\displaystyle \frac{32{\mathrm{k}}_{\mathsf{\tau }}\mathrm{L}}{\mathsf{\pi }{\mathrm{d}}^4}}$$

(10)

According to [37], the bending stiffness of the effective representative beam volume, accounting for the correction for the shear coefficient K in Equation (3), is given by [33]:

$${\mathrm{K}}_{\mathrm{b}\mathrm{e}\mathrm{n}\mathrm{d}}={\displaystyle \frac{\mathrm{E}\mathrm{I}}{\mathrm{L}}}{\displaystyle \frac{4+\mathsf{\phi }}{1+\mathsf{\phi }}},$$

(11)

$$\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\mathsf{\phi }={\displaystyle \frac{12\mathrm{E}\mathrm{I}}{\mathrm{G}{\mathrm{A}}_{\mathrm{k}}{\mathrm{L}}^2}}-\mathrm{s}\mathrm{h}\mathrm{e}\mathrm{a}\mathrm{r} \mathrm{d}\mathrm{e}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n} \mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{t} \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}$$

(12)

$${\mathrm{A}}_{\mathrm{k}}={\displaystyle \frac{\mathrm{A}}{\mathrm{K}}}=\mathrm{A}{\displaystyle \frac{(7+6\mathsf{\nu })}{6(1+\mathsf{\nu })}}.$$

(13)

Taking into account relations (10), we obtain the expression for the shear deformation constant $\mathsf{\phi }$:

$$\phi ={\displaystyle \frac{9 {\mathrm{k}}_{\mathrm{r}} {\mathrm{d}}^4 (1+\nu )}{16{{\mathrm{k}}_{\mathsf{\tau }}\mathrm{L}}^2\left(7+6\nu \right) }}.$$

(14)

By equating the bending stiffness of the effective representative bond volume (11) to the molecular mechanics force constant for bond bending (${\mathrm{K}}_{\mathrm{b}\mathrm{e}\mathrm{n}\mathrm{d}}= {\mathrm{k}}_{\mathsf{\theta }}$), and substituting (14) and (10) into (11), we obtain the following relation, which establishes the connection between the effective diameter $\mathrm{d}$ and the Poisson’s ratio $\nu$ through the molecular mechanics force constants:

$${\mathrm{k}}_{\mathsf{\theta }}= {\displaystyle \frac{ {\mathrm{k}}_{\mathrm{r}}{\mathrm{d}}^2}{16}} {\displaystyle \frac{ 448 {{\mathrm{k}}_{\mathsf{\tau }} L}^2 +9{\mathrm{k}}_{\mathrm{r}}{\mathrm{d}}^4 +384 {{\mathrm{k}}_{\mathsf{\tau }} L}^2 \nu +9 {\mathrm{k}}_{\mathrm{r}} {\mathrm{d}}^4\nu }{ {112 {\mathrm{k}}_{\mathsf{\tau }}L}^2 + {9\mathrm{k}}_{\mathrm{r}}{\mathrm{d}}^4 +96{ {\mathrm{k}}_{\mathsf{\tau }}L}^2 \nu + {9\mathrm{k}}_{\mathrm{r}} {\mathrm{d}}^4\nu }}.$$

(15)

After combining like terms in (15), we obtain a sixth-order algebraic Equation containing nonzero coefficients only for even powers of $\mathrm{d}$, or, in other words, a bicubic Equation with respect to the sought diameter $\mathrm{d}$, namely:

$${\mathrm{F}\left(\mathsf{\nu },\mathrm{d}\right)=9{{\mathrm{k}}_{\mathrm{r}}\mathrm{k}}_{\mathrm{r}}\left(1+\mathsf{\nu }\right){\mathrm{d}}^6-{144\mathrm{k}}_{\mathrm{r}}{\mathrm{k}}_{\mathsf{\theta }}\left(1+\mathsf{\nu }\right){\mathrm{d}}^4+64\mathrm{k}}_{\mathrm{r}}{{\mathrm{k}}_{\mathsf{\tau }}\mathrm{L}}^2\left(7+6\mathsf{\nu }\right){\mathrm{d}}^2-256{\mathrm{k}}_{\mathsf{\tau }}{\mathrm{k}}_{\mathsf{\theta }}{\mathrm{L}}^2\left(7+6\mathsf{\nu }\right)=0.$$

(16)

Similarly, we obtain the Equation for calculating the Poisson’s ratio $\mathsf{\nu }$:

$$\mathsf{\nu }={\displaystyle \frac{-9{{\mathrm{k}}_{\mathrm{r}}\mathrm{k}}_{\mathrm{r}}{\mathrm{d}}^6+{144\mathrm{k}}_{\mathrm{r}}{\mathrm{k}}_{\mathsf{\theta }}{\mathrm{d}}^4-448\left[{\mathrm{k}}_{\mathrm{r}}{{\mathrm{k}}_{\mathsf{\tau }}\mathrm{L}}^2{\mathrm{d}}^2-4{\mathrm{k}}_{\mathsf{\tau }}{\mathrm{k}}_{\mathsf{\theta }}{\mathrm{L}}^2\right]}{9{{\mathrm{k}}_{\mathrm{r}}\mathrm{k}}_{\mathrm{r}}{\mathrm{d}}^6-{144\mathrm{k}}_{\mathrm{r}}{\mathrm{k}}_{\mathsf{\theta }}{\mathrm{d}}^4+384 [{\mathrm{k}}_{\mathrm{r}}{{\mathrm{k}}_{\mathsf{\tau }}\mathrm{L}}^2{\mathrm{d}}^2-4{\mathrm{k}}_{\mathsf{\tau }}{\mathrm{k}}_{\mathsf{\theta }}{\mathrm{L}}^2]}}.$$

(17)

In [33], based on the same methodology, an expression for ${\mathrm{k}}_{\mathsf{\tau }}$ was obtained, similar to expression (15), namely:

$${\mathrm{k}}_{\mathsf{\tau }}= {\displaystyle \frac{ {\mathsf{\pi }\mathrm{k}}_{\mathrm{r}}{\mathrm{d}}^2}{16\mathsf{\pi }}} {\displaystyle \frac{ 448 {{\mathrm{k}}_{\mathsf{\theta }} L}^2{\mathsf{\pi }}^2+384 L^2\mathsf{\pi }{\mathrm{k}}_{\mathsf{\theta }}\nu +9{\mathrm{k}}_{\mathrm{r}}{\mathsf{\pi }}^2{d}^4 + 9 {\mathrm{k}}_{\mathrm{r}}{\mathsf{\pi }}^2{\mathrm{d}}^4\nu }{ {112{\mathrm{k}}_{\mathsf{\theta }} L}^2{\mathsf{\pi }}^2+96 L^2{{\mathsf{\pi }}^2\mathrm{k}}_{\mathsf{\theta }}\nu + {9 \mathrm{k}}_{\mathrm{r}}{\mathrm{d}}^4{\mathsf{\pi }}^2 + {9\mathrm{k}}_{\mathrm{r}}{\mathrm{d}}^4{\mathsf{\pi }}^2\nu }}.$$

(18)

Equation (16), after a change of variables, is reduced to an algebraic cubic Equation. Solving it using Cardano’s method with the force constants corresponding to the C–C bond [2] (${\mathrm{k}}_{\mathrm{r}}=651.694{\displaystyle \frac{\mathrm{n}\mathrm{N}}{\mathrm{n}\mathrm{m}}},{\mathrm{k}}_{\mathsf{\theta }}=0.875{\displaystyle \frac{\mathrm{n}\mathrm{N}\mathrm{n}\mathrm{m}}{{\mathrm{r}\mathrm{a}\mathrm{d}}^2}},{\mathrm{k}}_{\mathsf{\tau }}=0.278{\displaystyle \frac{\mathrm{n}\mathrm{N}\mathrm{n}\mathrm{m}}{{\mathrm{r}\mathrm{a}\mathrm{d}}^2}}$), we obtain the curve of the C–C bond diameter ${\mathrm{d}}_{c-c}$ as a function of the Poisson’s ratio $\mathsf{\nu }$, shown in Figure 2.

Comparison of the obtained results with the data presented in [33] demonstrates their close agreement (Figure 2).

on the Poisson’s ratio $\mathsf{\nu }$.

3. Results

The results of the function $\mathrm{F}\left(\mathsf{\nu },\mathrm{d}\right)$in Equation (16), with the Poisson’s ratio $\mathsf{\nu }$ varied in the range from 0 to 0.5, which corresponds to physically admissible values of this parameter for continuum materials, are shown in Figure 3 and Figure 4.

for different values of Poisson’s ratio $\mathsf{\nu }$ in the range from 0 to 0.5.

depending on the Poisson’s ratio $\mathsf{\nu }$ in the range from 0 to 0.5.

As can be seen from the graphs in Figure 3 and Figure 4, for Poisson’s ratio values within this range, the C–C bond diameter $\mathrm{d}$ varies from 0.08393 nm at $\mathsf{\nu }=0$ to 0.08474 nm at $\mathsf{\nu }=0.5$. The mean C–C bond diameter ${\mathrm{d}}_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}$corresponds approximately to 0.0844 nm at a Poisson’s ratio of $\mathsf{\nu }=0.32$.

The obtained numerical range of $\mathrm{d}$ for the C-C bond indicates that, when the Poisson’s ratio is confined within physically admissible limits for continuum media, the effective diameter of the C–C bond remains nearly constant and exhibits only a weak dependence on ν. From the standpoint of structural mechanics, such behavior is expected, because Poisson’s ratio primarily affects the redistribution of transverse shear stiffness, rather than the longitudinal stiffness of the representative volume.

In this context, the obtained diameter can be interpreted as the lower-bound estimate of the effective transverse cross-section of the bond, which ensures equivalence between the molecular bending stiffness ${\mathrm{k}}_{\mathsf{\theta }}$ and the bending stiffness of a representative continuum element modeled within the Timoshenko framework. In contrast to the widely accepted nominal diameter of $\mathrm{d}$ = 0.147 nm, which is directly associated with the equilibrium spatial configuration of atomic orbitals, the present estimate represents a mechanically consistent characteristic derived strictly from energy-equivalence considerations. The presented correspondence between continuum and atomistic descriptions provides a fundamentally justified measure of the bond diameter, which prevents overestimation of flexibility under bending and avoids unphysical Poisson’s ratio values (e.g., ν > 0.5) that often arise when using simplified Euler–Bernoulli formulas.

Thus, for engineering calculations of nanostructures, such as carbon nanotubes, it is recommended to use the effective C–C bond diameter ${\mathrm{d}}_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}=0.0844 nm$, which, according to Equation (17), corresponds to a Poisson’s ratio of $\mathsf{\nu }=0.32$. This value of ${\mathrm{d}}_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}$ may therefore serve as a reference geometric parameter in subsequent derivations of homogenized physical and mechanical characteristics for periodic carbon-containing nanostructures (CNT arrays, graphene sheets, graphite laminates) under both static deformation and dynamic loading, including the propagation of waves at characteristic nanoscale lengths.

Let us consider the physical meaning of these results.

As demonstrated in [21], most chemical bonding concepts can be adequately described using the Electron Localization Function (ELF). The distribution of electron density is a fundamental characteristic of a molecule, determining its physical and chemical properties. The effective overlap region of the C–C σ-bond can be evaluated either by performing Density Functional Theory (DFT) calculations [20] or from molecular visualization moiré patterns [21].

The work [21] presents studies on the ELF distribution for various molecules and notes that the central positions of C atoms within molecular structures act as invariant attractors. This indicates that the topology of structurally invariant attractors (the centroids of structural elements) is always preserved regardless of the chemical correlation between neighboring groups (changes in bond angles or unit cell dimensions), as long as the number of primary neighbors around each invariant attractor remains unchanged [21]. Based on this, assuming that the distance between two carbon atoms is 0.142 nm, and by appropriately scaling the data shown in Figure 5, the effective diameter of the C–C bond’s electron density cross-section is found to vary in the range of approximately 0.08 - 0.09 nm. These estimates of the effective transverse width of the C–C σ-bond overlap are in good agreement with the results obtained from DFT calculations, which lie in the range of approximately 0.07 -

0.09 nm [20], as well as with the results presented above.

molecule (extended Hückel method) [21].

ELF calculations and experimental data indicate that the electron density in the C–C σ-bond region is highly concentrated around the bond axis, and the “characteristic bond width,” where the density remains high—e.g., at least 95% of the maximum—is approximately 0.1 nm in diameter. These findings are corroborated by numerous studies on ELF analysis for graphene and diamond [22].

To describe the transverse distribution of electron density, a smooth normal distribution is used, which adequately captures the symmetry and shape of the σ-bond electron cloud. The adopted model allows for the analytical calculation of electron density fractions as a function of the distance from the bond axis.

Assuming a Gaussian distribution for the function $\mathsf{\rho }\left(\mathrm{r}\right)$ along the radius $\mathrm{r}$ perpendicular to the C–C bond axis, we can write:

$$\mathsf{\rho }\left(\mathrm{r}\right)={\mathsf{\rho }}_{\mathrm{m}\mathrm{a}\mathrm{x}}{\mathrm{e}}^{(-{\displaystyle \frac{{\mathrm{r}}^2}{2{\mathsf{\xi }}^2}})}$$

(19)

where $\mathrm{r}$ is the distance from the bond axis (transverse direction) in nm; ${\mathsf{\rho }}_{\mathrm{m}\mathrm{a}\mathrm{x}}$ describes the maximum density at the center, and $\mathsf{\xi }=0.015$ nm is the effective width parameter of the distribution.

Based on Equation (19), a normalized curve of the Gaussian electron density distribution (ELF), ${\mathrm{F}}_{\mathrm{E}\mathrm{L}\mathrm{F}}={\displaystyle \frac{\mathsf{\rho }\left(\mathrm{r}\right)}{{\mathsf{\rho }}_{\mathrm{m}\mathrm{a}\mathrm{x}}}}={\mathrm{e}}^{(-{\displaystyle \frac{{\mathrm{r}}^2}{2{\mathsf{\xi }}^2}})}$, along the radius (perpendicular to the bond axis) is obtained, as shown in Figure 6.

To estimate the radius of the region containing a specified percentage of the electron density, an integral calculation of the cumulative density is performed, assuming a cylindrical geometry for the bond. Based on the chosen profile, the radii corresponding to 20%, 50%, and 90% of the cumulative electron density are determined. These radii uniquely define the transverse size of the region of maximum electron localization. A 2D map of the electron density distribution in the transverse section of the C–C bond, corresponding to the Gaussian distribution given by Equation (19), is shown in Figure 7. Calculations show that the transverse radius of the region containing from 90% to 99% of the electron density increases from 0.04 nm to 0.05 nm. The difference between these cumulative density levels reflects the exponential decay of electron density along the transverse radius of the bond.

Determining the transverse diameter of the bond’s electron cloud is crucial for refining the distribution of mechanical loads in nanostructures, as the thickness of the electron density influences the calculation of bending and torsional stiffnesses when these values are incorporated into multiscale models for an accurate description of the mechanical behavior of carbon materials. Analysis of the electron density distribution curve indicates that the effective radius of the C–C bond lies within the range $\mathrm{r}\approx 0.04\mathrm{n}\mathrm{m}÷0.045\mathrm{n}\mathrm{m}$, while the corresponding effective diameter falls within $\mathrm{d}\approx 0.08\mathrm{n}\mathrm{m}÷0.09\mathrm{n}\mathrm{m}$.

The obtained value of the average C–C bond diameter, according to Equation (16), ${\mathrm{d}}_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}=0.0844\mathrm{n}\mathrm{m}$, is fully consistent with these data. The deviation from the average diameter ranges from –0.514% to 0.43%, which lies within the accuracy limits of the molecular mechanics force constants ${\mathrm{k}}_{\mathrm{r}},{\mathrm{k}}_{\mathsf{\theta }},{\mathrm{k}}_{\mathsf{\tau }}$, and the adopted shear coefficient $\mathrm{K}={\displaystyle \frac{6(1+\mathsf{\nu })}{(7+6\mathsf{\nu })}}$.

Given the known average value of the C–C bond diameter, the effective Young’s modulus of the C–C bond can be determined based on the generalized Morse potential $\mathrm{U}\left(\mathrm{r}\right)$, which represents a semi-empirical approximation of the true interatomic interaction potential and is expressed as [14]:

$$\mathrm{U}\left(\mathrm{r}\right)={\mathrm{D}}_{\mathrm{e}}\left[{ (1-{\mathrm{e}}^{-\mathsf{\alpha }\left(\mathrm{r}-{\mathrm{r}}_0\right)})}^2-1\right]={\mathrm{D}}_{\mathrm{e}}({\mathrm{e}}^{-2\mathsf{\alpha }(\mathrm{r}-{\mathrm{r}}_0)}-{2\mathrm{e}}^{-\mathsf{\alpha }(\mathrm{r}-{\mathrm{r}}_0)}),$$

(20)

where $\mathrm{U}\left(\mathrm{r}\right)$ is the specific potential energy of the bond, ${\mathrm{D}}_{\mathrm{e}}$ is the bond dissociation energy between the atoms, $\mathsf{\alpha }$ is a parameter defining the curvature of the potential well, $\mathrm{r}$ is the interatomic distance (nm), and ${\mathrm{r}}_0$ is the equilibrium bond length between the atoms.

In the case of bond deformation along the longitudinal axis, the strain component $\mathsf{\epsilon }$ is given by

$$\mathsf{\epsilon }=(\mathrm{r}-{\mathrm{r}}_0)/{\mathrm{r}}_0,$$

(21)

and Equation (21) takes the form:

$$\mathrm{U}\left(\mathsf{\epsilon }\right)={\mathrm{D}}_{\mathrm{e}}({\mathrm{e}}^{-2\mathsf{\alpha }{\mathrm{r}}_0\mathsf{\epsilon }}-{2\mathrm{e}}^{-\mathsf{\alpha }{\mathrm{r}}_0\mathsf{\epsilon }}),$$

(22)

The physical, unnormalized components of the elastic strain tensor are determined according to the well-known relation [20] ${\mathrm{C}}^{\mathrm{i}\mathrm{j}\mathrm{k}\mathrm{l}}$ = ${\displaystyle \frac{{\partial }^2\mathrm{U}\left({\mathsf{\epsilon }}_{\mathrm{i}\mathrm{j}}\right)}{\partial {\mathsf{\epsilon }}_{\mathrm{i}\mathrm{j}}\partial {\mathsf{\epsilon }}_{\mathrm{k}\mathrm{l}}}}{|}_{{\mathsf{\epsilon }}_{\mathrm{i}\mathrm{j}}=0,{\mathsf{\epsilon }}_{\mathrm{k}\mathrm{l}}=0}$, and the expression for calculating the elastic modulus is obtained after performing the normalization operation under the condition that $\mathsf{\epsilon }=0\mathrm{a}\mathrm{n}\mathrm{d}{\mathrm{k}}_{\mathrm{r}}={2D}_e{\left(\alpha \right)}^2$:

$$\mathrm{E}={\displaystyle \frac{1}{\mathrm{V}}}{\displaystyle \frac{{\partial }^2\mathrm{U}\left(\mathsf{\epsilon }\right)}{\partial \mathsf{\epsilon }\partial \mathsf{\epsilon }}}{|}_{\mathsf{\epsilon }=0}={\displaystyle \frac{2{\mathrm{D}}_{\mathrm{e}}{\left(\mathsf{\alpha }{\mathrm{r}}_0\right)}^2}{\mathrm{V}}}={\displaystyle \frac{8{\mathrm{D}}_{\mathrm{e}}{\left(\mathsf{\alpha }{\mathrm{r}}_0\right)}^2}{\mathsf{\pi }{\mathrm{d}}^2{\mathrm{r}}_0}}$$

(23)

The shear modulus $\mathrm{G}$ of the representative volume of the C–C bond is determined according to the established relation:

$$\mathrm{G}={\displaystyle \frac{\mathrm{E}}{2(1+\mathsf{\nu })}}={\displaystyle \frac{4{\mathrm{D}}_{\mathrm{e}}{\left(\mathsf{\alpha }{\mathrm{r}}_0\right)}^2}{\mathsf{\pi }{\mathrm{d}}^2{\mathrm{r}}_0(1+\mathsf{\nu })}}$$

(24)

The functional dependencies of the Young's modulus on the C–C bond diameter and the shear modulus on Poisson's ratio are shown in Figure 8 and Figure 9, respectively.

of the C–C bond.

4. Discussion

The obtained values of the transverse diameter of the C–C bond’s electron density localization region are crucial for refining the parameters of mechanical and multiscale models of carbon materials. In classical continuum models, which use the Timoshenko or Euler–Bernoulli beam approximations to describe the behavior of carbon nanostructures, the C–C bond is often treated as an element with equivalent mechanical properties, such as Young’s modulus, shear modulus, and geometric characteristics of its cross-section. However, this approach inherently assumes a simplified bond geometry and neglects the actual distribution of electron density. Studies conducted within quantum-mechanical frameworks, in particular using Density Functional Theory (DFT), show that the electron density of the σ C–C bond exhibits a spatial distribution resembling a Gaussian profile elongated along the interatomic axis. The isolines of this density display a smooth exponential decay from the bond center, as confirmed in [5]. Therefore, the actual cross-sectional area is determined not by atomic radii or tabulated parameters, but by the electron density distribution that governs the bond’s strength and stiffness. Consequently, defining an effective cross-sectional diameter corresponding to a specified fraction of the electron density is essential for the accurate interpretation of phenomenological elastic properties of the representative volume of diatomic systems.

Using the values obtained in this work,${\mathrm{d}}_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}=0.0844\mathrm{n}\mathrm{m}$, corresponding to 90–99% of the bond’s electron density, ensures consistency between molecular-mechanical and continuum models. The relevance of accurately determining the geometric parameters of the C–C bond stems from their impact on the calculation of bending stiffness, which is highly sensitive to the moment of inertia I ∼ d⁴; on torsional stiffness, determined by the polar moment J ∼ d⁴; and on the relationship between bending and torsional modes. Accounting for these factors is crucial when modeling the behavior of carbon nanotubes and graphene fragments. In particular, studies [29], [36], [40] have shown that errors in the chosen bond diameter of even 10-20% can lead to significant deviations in predicted stiffness values, especially for small-diameter nanotubes. Table 1 presents the values of the effective C–C bond thickness, Young’s modulus, and Poisson’s ratio used in the calculations of single-walled carbon nanotubes.

Figure 10 presents a comparison of the C–C bond diameter–Poisson’s ratio correspondence obtained in this work with results reported by other authors [29,33,38,41].

The oisson’s ratio ranges from 0.149 to 0.34 (Δ = 56%), while the C–C bond diameter varies from 0.066 nm to 0.0896 nm (Δ = 26%). The relatively high discrepancies in these parameters lead to an even greater spread in the values obtained for the Young’s modulus (Δ = 61%) (Table 1). This highlights the need for further studies aimed at refining models of diatomic covalent interactions.

In conclusion, it should be noted that the values of the effective bond diameter obtained in this work correspond to the actual electron density distribution in diatomic systems and enable:

• the elimination of discrepancies between continuum and atomistic models;
• the reconciliation of the geometric and physical properties of the C–C bond with molecular mechanics parameters, ensuring compliance with the thermodynamic stability conditions of deformable solid mechanics;
• the correct and physically justified application of the effective diameter concept in the engineering analysis and calculation of the mechanical properties of carbon nanostructures.

5. Conclusions

In this study, the fundamental problem of correlating the atomic bonding parameters of the C–C bond with its mechanical properties at the continuum model level is addressed. The use of traditional hypotheses for the Euler–Bernoulli beam theory proves inapplicable when directly utilizing molecular mechanics data, since neglecting shear deformation leads to nonphysical values for the Poisson's ratio (ν$\ge 0$.5).

To eliminate this inconsistency, an approach accounting for transverse shear based on the Timoshenko beam theory is proposed, providing an energetically consistent description of stretching, bending, and torsion of an individual C–C bond. By linking the molecular constants ${\mathrm{k}}_{\mathrm{r}}$, ${\mathrm{k}}_{\mathsf{\theta }}$, and ${\mathrm{k}}_{\mathsf{\tau }}$ to the structural parameters of an energetically equivalent representative-volume model, a sixth-order Equation for the effective bond diameter ${\mathrm{d}}_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}},$ dependent on the Poisson’s ratio, is obtained.

The solutions of this Equation show that the diameter $\mathrm{d}$ exhibits only a weak dependence on the Poisson’s ratio $\nu$ within its physically justified range of 0-0.5. Moreover, full compatibility between the atomistic and continuum descriptions of the two-atom C–C system is maintained across the entire range of these parameters. This range also corresponds to the conditions of thermodynamic stability, in which macroscopic system variables, such as temperature, pressure, volume, and entropy, undergo only small fluctuations around their equilibrium values when the system is isolated from the environment.

The proposed method therefore establishes a robust and physically consistent multiscale framework for modeling carbon nanostructures, such as carbon nanotubes and graphene, and provides a justified basis for determining effective geometric and elastic parameters of C–C bonds in continuum-based engineering analyses.

Author Contributions

Conceptualization, O.H. and S.A.; methodology, O.H. and S.A.; model validation, O.H.; formal analysis, O.H.; investigation, M.W. and S.P.; data curation, O.H.; writing—original draft preparation, O.H.; writing—review and editing, S.A.; supervision, O.H. and S.A. All authors have read and agreed to the published version of the manuscript.

Funding

This research study is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Project-ID: 554884501"Multiscale model of deformation and strength characteristics of dispersed carbon nanotube-reinforced polymers".

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.