Elastic Properties of Single-Walled Phosphide Nanotubes: Numerical Simulation Study

After a large-scale investigation into carbon nanotubes, significant research efforts have been devoted to discovering and synthesizing other nanotubes formed by chemical elements other than carbon. Among them, non-carbon nanotubes based on compounds of the elements of the 13th group of the periodic table and phosphorus. These inorganic nanotubes have proved to be more suitable candidates than carbon nanotubes for the construction of novel electronic and optical-electronic nano-devices. For this reason, until recently, mainly the structural and electrical properties of phosphide nanotubes were investigated, and studies to understand their mechanical behavior are infrequent. In the present work, the elastic properties of single-walled boron phosphide, aluminum phosphide, gallium phosphide and indium phosphide nanotubes were numerically evaluated using a nanoscale continuum modelling (also called molecular structural mechanics) approach. The force field constants required to assess the input parameters for numerical simulations were calculated for boron phosphide, aluminum phosphide, gallium phosphide and indium phosphide nanostructures using two different methods. The influence of input parameters on the elastic properties evaluated by numerical simulation was studied. A robust methodology to calculate the surface elastic moduli of phosphide nanotubes is proposed.


Introduction
Achievements in the fabrication of carbon nanotubes (CNTs) and their successful employing in numerous applications have given impetus to prediction and synthesis of new one-dimensional (1D) tubular structures with a honeycomb atomic arrangement. These nanostructures with graphene-like hexagonal lattice are based on compounds of chemical elements other than carbon. Among the examples of such compounds are those based on the elements of the 13th group of the periodic table, such as boron (B), aluminum (Al), gallium (Ga) and indium (In), which are able to establish with phosphorus (P) honeycomb diatomic arrangements, forming nanotubes (NTs) of boron phosphide (BP), aluminum phosphide (AlP), gallium phosphide (GaP) and indium phosphide (InP). These phosphorus-based diatomic NTs are, in general, broadband gap semiconductors with valuable electronics and optoelectronics properties and have promising applications as light-emitting devices operating in the visible range, such as light-emitting diodes (LEDs) [1][2][3][4], solar cells [5], building blocks for nano-integrated circuits [6,7] and parts in quantum electronic devices [8,9].
The synthesis of tubular phosphide nanostructures has been a challenge so far. Bakkers and Verheijen [8] synthesized for the first time crystalline InPNTs employing the vaporliquid-solid (VLS) laser ablation method and not using a template. Yin et al. [10] proposed a simple template-free thermal chemical process in a conventional furnace with controlled reaction temperature and gas flow for synthetizing single-crystalline InPNT. Palit et al. [11], in their recent study, grew single-crystalline InPNTs using the Metal-Organic Chemical Vapor Deposition (MOCVD) method in patterned template developed on germanium substrate by electron beam lithography (ELB) technique. Wu et al. [2] synthesized crystalline GaPNTs using a chemical reaction in a high-temperature tubular furnace and suggested that the growth of the GaP nanotube occurs according to the VLS mechanism [12].
The remaining phosphide nanotubes, BPNTs and AlPNTs, have not yet been synthesized, but these nanotubes were theoretically predicted. Mirzaei and Giahi [13], and Mirzaei and Meskinfam [14] used density functional theory (DFT) calculations to study optimized geometry and characterize the electronic structure of BPNTs. Theoretical results were also obtained regarding the functionalization of single-walled BPNTs and their suitability for drug carriage, using DFT calculations, in two works by Zahra Sayyad-Alangi et al. [15,16]. Lisenkov et al. [17] calculated the equilibrium geometry and energetic stability of AlPNTs, employing DFT, and Mirzaei and Mirzaei [18] examined the electronic structure of AlPNTs also using DFT calculations.
It is worth mentioning that the characterization of structural and electronic properties of GaPNTs and InPNTs was also carried out. Mirzaei and Mirzaei [19] used DFT calculations to optimize structure of single-walled GaPNTs and to compute the electric field gradient (EFG) tensors for the optimized nanotubes structure. Kamal et al. [20] studied the geometric and electronic structures of single-walled GaPNTs, based on ab initio calculations. Srivastava et al. [21] also employed the ab initio method to investigate the stability, electronic band structure and transport properties of single-walled zigzag GaPNTs. The band structure of the InPNTs was studied by Palit et al. [11], using a theoretical model based on the experimental observation of spectral emission lines. Erkoç [22] computed the optimized geometry and electronic properties of single-walled armchair and zigzag InPNTs, using semi-empirical calculations. In their work, Muhsen et al. [23] employed DFT to analyze the stability and electronic properties of single-walled zigzag InPNTs.
As it is of great importance to investigate the prospective employment of phosphide nanotubes in nano-electronics, their structure and electronic properties have been the focus of the research attention so far [18][19][20][21][22][23][24][25]. Consequently, the mechanical properties of the phosphide NTs have been less studied, although the understanding of the NTs' mechanical behavior can ensure the robustness and appropriate functioning of nano-devices involving nanotubes as components (building blocks).
Until now, the mechanical properties of phosphide NTs have been evaluated in the works of Kochaev [26], and Jiang and Guo [27], both studies exploring theoretical approaches to this end. Kochaev [26] evaluated the surface Young's modulus, that is, the product of Young's modulus by the wall thickness of the nanotubes, and the Poisson's ratio of AlPNTs and GaPNTs, using ab initio simulation within the atomistic approach. Jiang and Guo [27] applied the "stick-and-spring" model to obtain closed-form analytical solutions for the surface Young's modulus and Poisson's ratio and computed these values for BPNTs, GaPNTs and InPNTs. The model of Jiang and Guo [27] was developed under the nanoscale continuum modelling (NCM), also called molecular structural mechanics (MSM), approach, where the bonds between two atoms in the diatomic hexagonal nanostructure are regarded as beams or springs. For successful modelling of these bonds, first, an appropriate choice of the force field constants is necessary, which allows computing of the elastic properties of beams (springs). As for most non-carbon nanotubes, with the exception of boron nitride NTs, the force field constants for phosphide NTs practically do not appear in the literature. To the best of our knowledge, only Jiang and Guo [27] have suggested a method to calculate force constants for BP, GaP and InP nanotubes.
The present study is a systematic evaluation of the elastic properties of single-walled boron phosphide, aluminum phosphide, gallium phosphide and indium phosphide nanotubes (SWBPNTs, SWAlPNTs, SWGaPNTs and SWInPNTs, respectively) in a wide range of chiral indices and diameters, by finite element (FE) simulation. The force field constants for BP, AlP, GaP and InP nanostructures were determined using two different calculation approaches. The influence of the input parameters, selected for the FE modelling and computed basing on two sets of the force field constants, on the elastic properties of phosphide NTs was investigated. A comprehensive analysis was performed on the nanotube wall thickness, required for the calculation of the Young's and shear moduli. As a result, a robust methodology was proposed to assess the surface Young's and shear moduli.

Atomic Structure of Phosphide Nanotubes
The atomic structure of single-walled phosphide nanotubes is characterized by the chiral vector, C h , and the chiral angle, Θ, as follows: where n and m are the chiral indices, both having integers values; a 1 and a 2 are the unit vectors of the diatomic hexagonal lattice. This lattice consists of the A13 atom, which is one of the 13th group of the periodic table, such as boron (B), aluminum (Al), gallium (Ga) or indium (In), and phosphorus (P) atom, as shown in Figure 1 for AlP lattice. The length of the unit vector a is expressed by a = √ 3a A13−P , where a A13−P is the equilibrium bond length. As can be seen in Table 1, in which the bond lengths of the phosphide NTs available in the literature are presented, there is no agreement regarding the a A13−P values. approaches. The influence of the input parameters, selected for the FE modelling and computed basing on two sets of the force field constants, on the elastic properties of phosphide NTs was investigated. A comprehensive analysis was performed on the nanotube wall thickness, required for the calculation of the Young's and shear moduli. As a result, a robust methodology was proposed to assess the surface Young's and shear moduli.

Atomic Structure of Phosphide Nanotubes
The atomic structure of single-walled phosphide nanotubes is characterized by the chiral vector, Ch, and the chiral angle, Θ, as follows: where n and m are the chiral indices, both having integers values; a 1 and a 2 are the unit vectors of the diatomic hexagonal lattice. This lattice consists of the A13 atom, which is one of the 13th group of the periodic table, such as boron (B), aluminum (Al), gallium (Ga) or indium (In), and phosphorus (P) atom, as shown in Figure 1 for AlP lattice. The length of the unit vector a is expressed by a=√3a A13-P , where a A13-P is the equilibrium bond length. As can be seen in Table 1, in which the bond lengths of the phosphide NTs available in the literature are presented, there is no agreement regarding the a A13-P values.

BP
AlP GaP InP a A13−P , nm 0.183 [28] 0.193 [29] 0.234 [17] 0.240 [26] 0.220 [26] 0.225 [28] 0.229 [25] 0.236 [29] 0.246 [28] 0.256 [29] Nanomaterials 2022, 12, 2360 4 of 30 The phosphide NTs are formed by rolling up the respective diatomic hexagonal sheet into a cylinder of diameter, D n , which is considered the nanotube diameter and expressed by where n and m are the chiral indices and a A13−P is the equilibrium bond length of the diatomic structure of the phosphide. The chiral indices, n and m, together with the magnitudes of the chiral angle, which are within the range of 0 • to 30 • (see Figure 1), allow defining three main symmetry groups of single-walled phosphide nanotubes: for zigzag NTs Θ = 0 • and n = 0, for armchair NTs Θ = 30 • and n = m, and for chiral NTs 0 • < Θ < 30 • and n = m = 0. The zigzag and armchair configurations are called non-chiral NTs. Non-chiral (zigzag and armchair) and chiral SWBPNTs, SWAlPNTs, SWGaPNTs and SWInPNTs, displayed in ascending order of the chiral angle, Θ, are represented schematically in Figure 2.
where n and m are the chiral indices and a A13-P is the equilibrium bond length of the diatomic structure of the phosphide. The chiral indices, n and m, together with the magnitudes of the chiral angle, which are within the range of 0° to 30° (see Figure 1), allow defining three main symmetry groups of single-walled phosphide nanotubes: for zigzag NTs Θ = 0° and n = 0, for armchair NTs Θ = 30° and n = m, and for chiral NTs 0° < Θ < 30° and n ≠ m ≠ 0. The zigzag and armchair configurations are called non-chiral NTs. Non-chiral (zigzag and armchair) and chiral SWBPNTs, SWAlPNTs, SWGaPNTs and SWInPNTs, displayed in ascending order of the chiral angle, Θ, are represented schematically in Figure 2.

Force Field Constants
According to Mayo et al. [30] and Rappé et al. [31], the total potential energy of a molecular system is presented as the sum of the energy terms due to bonded and nonbonded interactions: The energy of bonded interactions, which are presented in Figure 3, is expressed through those associated with bond stretching, Ur, bond bending, Uθ, dihedral angle torsion, Uφ, and out-of-plane torsion, Uω, as follows: The energy of non-bonded interactions comprises the van der Waals, U vdW , electrostatic, U Q , and explicit hydrogen bonds, U Hbs , terms: U non-bond =U vdW +U Q +U Hbs .
In molecular systems, as in the case of phosphide nanotubes, U non-bond term can be omitted due to its smallness when compared with U bond term [27]; thus, U total =U bond . The four terms of Equation (5), which contribute to the total potential energy, U total , are given by the following expressions: where k r , k θ , k ϕ and k ψ are the bond stretching, bond bending, dihedral torsion and inversion force constants, respectively; r, θ, φ and ψ are bond length, bond angle, dihedral angle and inversion angle, respectively, and r 0 , θ 0 , ϕ 0 and ψ 0 correspond to the equilibrium position (see Figure 3).

Figure 2.
Configurations of (14, 0) zigzag, (12,6)  According to Mayo et al. [30] and Rappé et al. [31], the total potential energy of a molecular system is presented as the sum of the energy terms due to bonded and non-bonded interactions: The energy of bonded interactions, which are presented in Figure 3, is expressed through those associated with bond stretching, U r , bond bending, U θ , dihedral angle torsion, U φ , and out-of-plane torsion, U ω , as follows: The energy of non-bonded interactions comprises the van der Waals, U vdW , electrostatic, U Q , and explicit hydrogen bonds, U Hbs , terms: In molecular systems, as in the case of phosphide nanotubes, U non−bond term can be omitted due to its smallness when compared with U bond term [27]; thus, U total = U bond . The four terms of Equation (5), which contribute to the total potential energy, U total , are given by the following expressions: where k r , k θ , k φ and k ψ are the bond stretching, bond bending, dihedral torsion and inversion force constants, respectively; r, θ, φ and ψ are bond length, bond angle, dihedral angle and inversion angle, respectively, and r 0 , θ 0 , φ 0 and ψ 0 correspond to the equilibrium position (see Figure 3).  The dihedral angle torsion term describes the torsion interaction for two interatomic bonds linked by means of a common bond [30,31] (see Figure 3c). The out-of-plane torsion or inversion term illustrates how difficult it is to rotate three bonds (when one atom in a diatomic honeycomb lattice is bonded to three other atoms), keeping them in the same plane [30] (see Figure 3d). Under the assumption of small deformation and with the approximation that the variation of the dihedral angle, ϕ−ϕ 0 , is equal to that of the inversion angle, ψ−ψ 0 , the energies of the dihedral angle torsion and the out-of-plane torsion can be merged into a single equivalent term: where k τ is the torsional resistance force constant expressed by Thus, the total potential energy of the molecular system can be rewritten using three energy terms, U r , U θ and U τ , in the form of harmonic functions, as follows: The dihedral angle torsion term describes the torsion interaction for two interatomic bonds linked by means of a common bond [30,31] (see Figure 3c). The out-of-plane torsion or inversion term illustrates how difficult it is to rotate three bonds (when one atom in a diatomic honeycomb lattice is bonded to three other atoms), keeping them in the same plane [30] (see Figure 3d). Under the assumption of small deformation and with the approximation that the variation of the dihedral angle, (φ − φ 0 ), is equal to that of the inversion angle, (ψ − ψ 0 ), the energies of the dihedral angle torsion and the out-of-plane torsion can be merged into a single equivalent term: where k τ is the torsional resistance force constant expressed by Thus, the total potential energy of the molecular system can be rewritten using three energy terms, U r , U θ and U τ , in the form of harmonic functions, as follows: where ∆r, ∆θ and ∆φ are the bond stretching increment, bond angle bending variation and angle variation of the twist bond, respectively.
With regard to the calculation of the bond stretching, k r , and bond bending, k θ , force constants for the diatomic nanostructures, there are two established methods, one based on Universal Force Fields (UFF) [31], and the other using ab initio DFT computations in combination with the analytical expressions obtained using molecular mechanics (MM) [32,33]. Until now, these methods were mainly used for determination of the force field constants of boron nitride nanostructures (see, for example, [32][33][34][35]). Jiang and Guo [27], based on a UFF approach, proposed expressions for k r and k θ force constants for a wide class of non-carbon nanotubes, including BP, GaP and InPNTs.
According to Rappé et al. [31], the bond stretching constant, k r , in the UFF method is determined by the following expression: where Z * i and Z * j are the effective charges of the atoms of diatomic nanostructure; r ij and r ik are the lengths of the A13-P and P-A13 bonds, respectively, as presented in Figure 4, being r ij = r ik = a A13−P (see Figure 3a). where Δr, Δθ and Δϕ are the bond stretching increment, bond angle bending variation and angle variation of the twist bond, respectively. With regard to the calculation of the bond stretching, kr, and bond bending, kθ, force constants for the diatomic nanostructures, there are two established methods, one based on Universal Force Fields (UFF) [31], and the other using ab initio DFT computations in combination with the analytical expressions obtained using molecular mechanics (MM) [32,33]. Until now, these methods were mainly used for determination of the force field constants of boron nitride nanostructures (see, for example, [32][33][34][35]). Jiang and Guo [27], based on a UFF approach, proposed expressions for kr and kθ force constants for a wide class of non-carbon nanotubes, including BP, GaP and InPNTs.
According to Rappé et al. [31], the bond stretching constant, k r , in the UFF method is determined by the following expression: where Z i * and Z j * are the effective charges of the atoms of diatomic nanostructure; r ij and r ik are the lengths of the A13-P and P-A13 bonds, respectively, as presented in Figure 4, being r ij =r ik =a A13-P (see Figure 3a). The bond bending constant, k θ , in the UFF method is determined by the following expression [31]: where θ 0 is the angle between neighboring bonds in the diatomic nanostructure (see Figure 4) and r jk 2 =r ij 2 +r ik 2 −2r ij r ik cos θ 0 .
Rappé et al. [31], using UFF, showed that the bond bending constant, kθ, of the diatomic nanostructure depends on the three-body angles between the bond pairs A13-P-A13 and P-A13-P (see Figure 3b), which results in two different values for the bond bending constant. Moreover, in the same work, the relationship between two values for the bond bending constant, kθ1 and kθ2, and the effective charges of the atoms A13 and P (Z 1,2 * ) was proposed: The bond bending constant, k θ , in the UFF method is determined by the following expression [31]: where θ 0 is the angle between neighboring bonds in the diatomic nanostructure (see Figure 4) and r 2 jk = r 2 ij +r 2 ik −2r ij r ik cos θ 0 . Rappé et al. [31], using UFF, showed that the bond bending constant, k θ , of the diatomic nanostructure depends on the three-body angles between the bond pairs A13-P-A13 and P-A13-P (see Figure 3b), which results in two different values for the bond bending constant. Moreover, in the same work, the relationship between two values for the bond bending constant, k θ1 and k θ2 , and the effective charges of the atoms A13 and P (Z * 1,2 ) was proposed: The alternative method for obtaining the bond stretching and bond bending force constants is based on molecular mechanics (MM) and employs results from DFT calculations. The expressions, which result from MM analytical models for two-dimensional honeycomb Nanomaterials 2022, 12, 2360 8 of 30 diatomic nanostructures, allow relating the surface Young's modulus, E s , and the Poisson's ratio, ν, with the force field constants, k r , k θ1 and k θ2 , as follows [32,33]: The values of E s and ν can be obtained experimentally or from ab initio DFT computations (see, for example, [28,36]).
Thus, solving the system of Equations (14) and taking into account Equation (13), the following expressions can be derived to determine the bond stretching, k r , and bond bending, k θ1 and k θ2 , constants of the diatomic nanostructure: where E s and ν are the surface Young's modulus and Poisson's ratio of the diatomic sheet, respectively; Z * 1 and Z * 2 are the effective charges of the atoms; and r ij = a A13−P is the bond length. Literature data, which are required to calculate the bond stretching, k r , and bond bending, k θ1 and k θ2 , force constants for phosphide NTs, are summarized in Table 2. Table 2. Effective charges of atoms [31], and bond length, surface Young's modulus and Poisson's ratio, obtained from DFT computations [28], for phosphide nanostructures. To the best of our knowledge, the methods for calculating the torsional force constant, k τ , for phosphide nanostructures have not been appropriately explored and the values of k τ have not been reported so far. Jiang and Guo [27] proposed an expression to calculate the inversion force constant, k ψ , adopting the UFF method. Among the generic molecular force fields, DREIDING force field [30] allows the determination of the force constants, based only on the hybridization of atoms, regardless of the atoms involved. For diatomic phosphide nanostructures, the DREIDING force field provides the dihedral torsion force constant, k φ = 25 kcal/mol, and the inversion force constant, k ψ = 40 (kcal/mol)/rad 2 . Consequently, the torsional resistance force constant, k τ , can already be calculated using Equation (9). Table 3 presents the bond stretching, k r , bond bending, k θ1 and k θ2 , force constants obtained using both calculation methods, UFF (case 1) and DFT+MM (case 2), and the torsional resistance force constant, k τ , taken from DREIDING, for phosphide nanostructures. As far as we know, there are no values of E s and ν, obtained experimentally or with the help of DFT calculations, which would allow the calculation of bond stretching and bond bending force constants for AlP nanostructures, using Equations (15) and (16). Thus, the study of numerical simulation of the elastic properties of the AlPNTs is limited by case 1 of the set of force constants.

Equivalent Properties of Elastic Beams
In the present study, the NCM/MSM approach was used to evaluate the elastic properties of phosphide NTs. This approach makes use of the connection between the potential energies of bond interactions (see Equation (10) and Figure 3) and the strain energies associated with axial, bending and torsional elastic deformations of equivalent beam elements. Figure 5 shows the beam element undergoing pure tension (a), pure bending (b) and pure torsion (c). and bond bending force constants for AlP nanostructures, using Equations (15) and (16). Thus, the study of numerical simulation of the elastic properties of the AlPNTs is limited by case 1 of the set of force constants.

Equivalent Properties of Elastic Beams
In the present study, the NCM/MSM approach was used to evaluate the elastic properties of phosphide NTs. This approach makes use of the connection between the potential energies of bond interactions (see Equation (10) and Figure 3) and the strain energies associated with axial, bending and torsional elastic deformations of equivalent beam elements. Figure 5 shows the beam element undergoing pure tension (a), pure bending (b) and pure torsion (c). The energies related to stretching, UA, bending, UB, and torsion, UT, of the beam with length, l, are expressed as follows: The energies related to stretching, U A , bending, U B , and torsion, U T , of the beam with length, l, are expressed as follows: where F A is the axial force, M B is the bending moment, and T is the torsion moment; E b and G b are the beam Young's and shear moduli, respectively; A b is the beam cross-section area, I b is the beam moment of inertia and J b is the beam polar moment of inertia; E b A b is the beam tensile rigidity, E b I b is the beam bending rigidity and G b J b is the beam torsional rigidity; ∆l is the beam axial stretching displacement, ω is the rotational angle at the ends of the beam and ∆ϑ is the relative rotation between the ends of the beam.
Equations (10) and (17)- (19) allow one to establish the equality between the potential energies related to bond interactions and strain energies associated with elastic deformation of the beam elements, i.e., U r = U A , U θ = U B and U τ = U T . In other words, the tensile, E b A b , bending, E b I b , and torsional, G b J b , rigidities of the beam with length l are expressed with the help of bond stretching, k r , bond bending, k θ , and torsional resistance, k τ , force constants [37]: Equations (20)- (22) together with the assumption that the beam length, l, is equal to the bond length, a A13−P , support the use of the NCM/MSM approach to model the mechanical response of phosphide NTs.
Assuming that the beam has a circular cross-section with diameter d, its cross-section area, A b , moment of inertia, I b , and polar moment of inertia, J b , are determined by the following expressions: The diameter, d, the Young's modulus, E b , and the shear modulus, G b , of the beam can be calculated comparing Equations (20)-(22) with Equations (23)- (25) and taking into account the two values of the bond bending constant, k θ1 and k θ2 , as follows: The Poisson's ratio of the beam can be assessed by the relationship derived from molecular mechanics [32,33,38] as follows: Thus, knowing the values of the force field constants, k r , k θ1 and k θ2 , k τ , for phosphide nanostructures (Table 3), the geometrical and elastic properties of the beam elements can be deduced, as shown in Table 4. Table 5 shows the geometrical characteristics of SWBPNTs, SWAlPNTs, SWGaPNTs and SWInPNTs of three main configurations, armchair (Θ = 30 • ), zigzag (Θ = 0 • ) and chiral (family of Θ = 19.1 • ), used in the FE analysis. The chiral indices of NTs were chosen in order to obtain structures with comparable diameters. The length of the NTs was chosen to be about 30 times greater than the NT diameter to ensure that the mechanical response of nanotubes is independent of the NT length [39].

Geometrical Characteristics of Phosphide Nanotubes and FE Analysis
The input parameters (two sets of these parameters for each NTs, except for SWAlPNTs) used for FE simulation of phosphide nanotubes are those presented in Table 4.  The meshes of SWBPNTs, SWAlPNTs, SWGaPNTs and SWInPNTs used in FE analyses were constructed with the help of the Nanotube Modeler© software (version 1.8) developed by JCrystalSoft (www.jcrystal.com accessed on 21 June 2022). The in-house application InterfaceNanotubes.NM [39] was used for conversion of the PDB (Program Database) files, obtained from the Nanotube Modeler© software, into the format usable in the ABAQUS ® code. Afterwards, the FE code ABAQUS ® was used to study the mechanical response of the phosphide NTs under conventional tensile, bending and torsion tests. For this, in order to carry out the respective numerical tests, the axial tensile force, F a , the transverse force, F t , and the torsional moment, T, were applied to one end of the NT, when the other end was fixed. In the torsion test, the nodes under loading were prevented from moving in the radial direction. The results taken from the FE analysis of the tensile, bending and torsion test are, respectively, the axial displacement, u a , the transverse displacement, u t , and the twist angle, ϕ. This makes it possible to determine the tensile, EA, bending, EI, and torsional, GJ, rigidities of the phosphide NT with a length L n by the following expressions:

Rigidities of SWBPNTs, SWAlPNTs, SWGaPNTs and SWInPNTs
The tensile, EA, bending, EI, and torsional, GJ, rigidities of SWBPNTs, SWAlPNTs, SWGaPNTs and SWInPNTs obtained by Equations (30)-(32) for the two cases (except SWAlPNTs) of the numerical simulation input values presented in Table 4, are plotted as a function of the nanotube diameter, D n , in Figures 6-8, respectively. For case 1 and case 2 of numerical simulation input parameters from Table 4, the rigidity values for chiral or non-chiral (zigzag and armchair) NTs follow the same trend with increasing D n . The EA, EI and GJ rigidities obtained for case 1 (UFF) are higher than those for case 2 (DFT + MM).
As for the cases of the single-walled carbon nanotubes (SWCNTs) [40,41] and singlewalled boron nitride nanotubes (SWBNNTs) [39], and for the phosphide NTs under study, the values of the tensile rigidity, EA, can be represented by a linear function of nanotube diameter, D n (see Figure 9a,b), and the values of bending, EI, and torsional, GJ, rigidities can be represented by a linear function of D 3 n (see Figure 9c-f). Similar to what was found in the authors' previous work for the SWBNNTs [39], the straight lines in Figure 9a-f can be expressed as follows: where α A13P , β A13P and γ A13P are the fitting parameters. The values of these parameters, obtained from the graphs in Figure 9a-f for single-walled phosphide NTs, are shown in Table 6. The accuracy of the evaluation of the EA, EI and GJ rigidities values analytically estimated with Equations (33)- (35), respectively, was verified through the comparison with the values of the EA, EI and GJ rigidities calculated by Equations (30)-(32), respectively, using the data taken from FE analysis. The mean differences between the EA, EI and GJ values evaluated analytically and those obtained by the FE analysis are shown in Table 7. It can be seen from Table 7 that Equations (33)-(35) allow an accurate calculation of the three rigidities of SWBPNTs, SWAlPNTs, SWGaPNTs and SWInPNTs. The mean difference does not exceed 0.85%, which is the greatest value observed for the bending rigidity.
In order to better comprehend the results of the tensile, EA, bending, EI, and torsional, GJ, rigidities presented in Figure 9, the evolutions of α A13P , β A13P and γ A13P fitting parameters with the bond length value, a A13P (see Table 2), considering the cases of SWBPNTs, SWAlPNTs, SWGaPNTs and SWInPNTs, are shown in Figure 10, for the two cases of input parameters. For case 1, all three fitting parameters decrease from SWBPNTs to SWAlPNTs, i.e., as the value of a A13−P increases, and then the values of α A13P , β A13P and γ A13P remain nearly unchanged when moving to SWInPNTs. For case 2, the fitting parameters α A13P , β A13P and γ A13P decrease with increasing bond length from SWBPNTs to SWInPNTs.
It should be noted that for case 1 of SWBPNTs, the ratio β A13P /γ A13P is about 1, which means that the EI and GJ rigidities are nearly equal. The value of this ratio, β A13P /γ A13P = 1.1, for SWBPNTs (case 2), SWAlPNTs, SWGaPNTs and SWInPNTs (case 1), indicates a certain difference between the values of bending, EI, and torsional, GJ, rigidities. This difference becomes higher for case 2 of SWInPNTs, for which β A13P /γ A13P = 1.2.   As for the cases of the single-walled carbon nanotubes (SWCNTs) [40,41] and singlewalled boron nitride nanotubes (SWBNNTs) [39], and for the phosphide NTs under study, the values of the tensile rigidity, EA, can be represented by a linear function of nanotube diameter, D n (see Figure 9a,b), and the values of bending, EI, and torsional, GJ, rigidities can be represented by a linear function of D n 3 (see Figure 9c-f).
(a) (b) Figure 10. Evolutions of (a) α A13P , and (b) β A13P and γ A13P fitting parameters with the bond length, a A13-P , for the two cases of input values in the numerical simulation of SWBPNTs, SWAlPNTs, SWGaPNTs and SWInPNTs.
It should be noted that for case 1 of SWBPNTs, the ratio β A13P γ A13P is about 1, which means that the EI and GJ rigidities are nearly equal. The value of this ratio, β A13P γ A13P = 1.1, for SWBPNTs (case 2), SWAlPNTs, SWGaPNTs and SWInPNTs (case 1), indicates a certain difference between the values of bending, EI, and torsional, GJ, rigidities. This difference becomes higher for case 2 of SWInPNTs, for which β A13P γ A13P = 1.2.

Effect of Nanotube Wall Thickness on the Calculation of Elastic Moduli
As previously deduced for SWCNTs [40,41] and SWBNNTs [39], the Young's and shear moduli of NTs structures can be calculated using the following expressions, respectively: Figure 10. Evolutions of (a) α A13P , and (b) β A13P and γ A13P fitting parameters with the bond length, a A13−P , for the two cases of input values in the numerical simulation of SWBPNTs, SWAlPNTs, SWGaPNTs and SWInPNTs.

Effect of Nanotube Wall Thickness on the Calculation of Elastic Moduli
As previously deduced for SWCNTs [40,41] and SWBNNTs [39], the Young's and shear moduli of NTs structures can be calculated using the following expressions, respectively: where EA, EI and GJ are the tensile, bending and torsional rigidities, respectively, and t n is the nanotube wall thickness. The EA, EI and GJ rigidities can be calculated from the results of the FE analysis by Equations (30)-(32), respectively, or evaluated analytically using Equations (33)-(35), respectively. To date, regarding the NT wall thickness, no t n value, observed experimentally or calculated by theoretical approaches, has been reported for phosphide nanotubes. Furthermore, there are no reliable data indicating that the wall thickness value of phosphide NTs can be adopted equal to 0.34 nm (the graphite interlayer spacing). For this reason, the Young's, E, and shear, G, moduli determined by Equations (36) and (37), respectively, are plotted as a function of the inverse of the nanotube wall thickness, 1/t n , (in the range 0.1 ≤ t n ≤ 1.5 nm) for phosphide NTs selected from Table 5, as shown in Figure 11 (for E values) and Figure 12 (for G values). All three NT configurations, zigzag, chiral and armchair, and the two cases of input parameters are considered.   The evolutions of both the Young's, E, and the shear, G, moduli as a function of the inverse nanotube wall thickness, 1/t n , follow a quasi-linear trend for NTs with diameter D n ≳ 1.7 nm. For phosphide NTs with diameters D n ≲ 1.7 nm, the deviation from the quasi-linear trend occurs for high values of NTs wall thickness, when t n is in the range between the nanotube diameter and half the nanotube diameter, D n ≲ t n ≲ D n 2 ⁄ . Thus, the smaller the value of D n , the more noticeable is the deviation from the quasi-linearity of the evolution of the Young's and shear moduli with 1/t n . The SWBPNTs, SWAlPNTs, SWGaPNTs and SWInPNTs with diameters D n ≲ 1.7 nm and wall thickness in the range D n ≲ t n ≲ D n 2 ⁄ behave as solid cylinders instead of hollow tubes, which influences the results of elastic moduli.

Surface Young's Modulus of Phosphide Nanotubes
The lack of knowledge about the reliable wall thickness value has caused the studies, concerning evaluation of the elastic properties of the non-carbon nanotubes (N-CNTs), to mainly focus on determining the N-CNTs' surface elastic moduli [42]. Regarding the elastic properties of the phosphide NTs, Kochaev [26] and Jiang and Guo [27], are the only authors who calculated the surface Young´s modulus of BP [27], AlP [26], GaP [26,27] and InPNTs [27]. Thus, to calculate the Young´s, E, and shear, G, moduli without the necessity of knowing the NT wall thickness and to facilitate the comparison with the literature results, a methodology for calculation of the surface Young's, E s , and shear, G s , moduli, based on the linear evolutions of E and G moduli as a function of the inverse of the NT thickness, 1/t n , is suggested.
As long as the wall thickness of phosphide nanotube is greater than half of its diameter, the Young's and shear moduli become quasi-linear functions of the inverse wall thickness (see Figures 11 and 12). The linear parts of the evolutions of E and G moduli The evolutions of both the Young's, E, and the shear, G, moduli as a function of the inverse nanotube wall thickness, 1/t n , follow a quasi-linear trend for NTs with diameter D n 1.7 nm. For phosphide NTs with diameters D n 1.7 nm, the deviation from the quasi-linear trend occurs for high values of NTs wall thickness, when t n is in the range between the nanotube diameter and half the nanotube diameter, D n t n D n /2. Thus, the smaller the value of D n , the more noticeable is the deviation from the quasi-linearity of the evolution of the Young's and shear moduli with 1/t n . The SWBPNTs, SWAlPNTs, SWGaPNTs and SWInPNTs with diameters D n 1.7 nm and wall thickness in the range D n t n D n /2 behave as solid cylinders instead of hollow tubes, which influences the results of elastic moduli.

Surface Young's Modulus of Phosphide Nanotubes
The lack of knowledge about the reliable wall thickness value has caused the studies, concerning evaluation of the elastic properties of the non-carbon nanotubes (N-CNTs), to mainly focus on determining the N-CNTs' surface elastic moduli [42]. Regarding the elastic properties of the phosphide NTs, Kochaev [26] and Jiang and Guo [27], are the only authors who calculated the surface Young's modulus of BP [27], AlP [26], GaP [26,27] and InPNTs [27]. Thus, to calculate the Young's, E, and shear, G, moduli without the necessity of knowing the NT wall thickness and to facilitate the comparison with the literature results, a methodology for calculation of the surface Young's, E s , and shear, G s , moduli, based on the linear evolutions of E and G moduli as a function of the inverse of the NT thickness, 1/t n , is suggested.
As long as the wall thickness of phosphide nanotube is greater than half of its diameter, the Young's and shear moduli become quasi-linear functions of the inverse wall thickness (see Figures 11 and 12). The linear parts of the evolutions of E and G moduli with 1/t n can be described by expressions E = α 1 (1/t n ) and G = α 2 (1/t n ), respectively, where α 1 and α 2 are the slopes of the corresponding straight lines. Taking into account that the surface Young's, E s , modulus is the Young's modulus, E, multiplied by the NT wall thickness, E s = Et n , and, likewise, the surface shear modulus, G s , is the shear modulus, G, multiplied by the NT wall thickness, G s = Gt n , it can be written as Thus, Equations (38) and (39) are the basis of the methodology to calculate the surface Young's, E s , and shear, G s moduli through the slope of the linear portion of the evolution of the corresponding elastic modulus (E or G) as a function of the inverse of the nanotube wall thickness.
In this context, the evolutions of the E and G moduli as a function of 1/t n , as shown in the examples in Figures 11 and 12, were plotted for all SWBPNTs, SWAlPNTs, SWGaPNTs and SWInPNTs in Table 5. Only the linear portions of the evolutions of the elastic moduli, with an R-squared value of at least 0.9997, which approximately corresponds to the NT wall thickness range t n D n /3, were considered for analysis. Then, to assess the surface Young's, E s , and shear, G s , moduli, the slopes of straight lines were determined (see equations (38) and (39), respectively). Figure 13 shows the evolutions of E s with the diameter of the nanotubes, D n , for SWBPNTs, SWAlPNTs, SWGaPNTs and SWInPNTs. The surface Young's modulus is nearly stable over the entire diameter range of the NTs studied, regardless of the case of the input parameters, the NTs chirality and the first element (B, Al, Ga, In) of phosphide compound forming the nanotube. As can be seen in Figure 13a, for case 1 of input parameters, the SWBPNTs have the highest value of E s of about 0.218 TPa·nm, and the value of E s decreases approximately by half for SWAlPNTs (E s = 0.094 TPa·nm), SWGaPNTs (E s = 0.107 TPa·nm) and SWInPNTs (E s = 0.096 TPa·nm). For case 2 of input parameters, the surface Young's modulus decreases from SWBPNTs to SWInPNTs, the E s value being approximately 2.2 and 3.2 times smaller for SWGaPNTs (E s = 0.078 TPa·nm) and SWInPNTs (E s = 0.055 TPa·nm), respectively, when compared with SWBPNTs (E s = 0.176 TPa·nm) (see Figure 13b). The values of E s calculated for case 1 are higher than those evaluated for case 2, whatever the NTs, as can be seen in Figure 13c for the SWBPNTs and SWInPNTs. Figure 13d compares the surface Young's modulus, E s , results obtained for case 2 of (n, n) and (n, 0) SWInPNTs with those available in the literature.
Among the literature results, the values of the surface Young's modulus reported by Jiang and Guo [27] for (n, n) and (n, 0) SWInPNTs, which are in better agreement with the E s values calculated in the present study, were chosen for comparison. A mean difference of about 7.9% was observed for nanotubes with diameter D n 1.0 nm, when comparing the E s results by Jiang and Guo [27] with those obtained for case 2. Comprehensive comparison with the literature results appears to be difficult due to the scarcity of studies on the evaluation of the phosphide NTs' surface Young's modulus, but discrepancies were noticed in the reported E s values and trends. For example, in an ab initio simulation study, Kochaev [26] observed, for SWAlPNTs and SWGaPNTs, a considerable increase of the surface Young's modulus, for (n, n) nanotubes when D n 1.2 nm and (n, 0) nanotubes when D n 0.9 nm, and then E s reaches the maximum, followed by a sharp decrease. An alternative trend reported in the literature for SWBPNTs, SWGaNTs and SWInPNTs [27] consists of a slight increase of the E s value when the nanotube diameter increases, and then E s becomes practically constant for NTs with diameters D n 0.7 nm (see Figure 13d for the case of SWInPNTs). Finally, the comparison with the available results was carried out as far as possible, as documented in Table 8.  In contrast to the reasonable agreement of the E s values for SWInPNTs assessed by Jiang and Guo [27] and those obtained for case 2 of the input parameters (see Figure 13c), less agreement is observed when the E s results reported in the same work [27] for SWBPNTs and SWGaPNTs are compared with the current results calculated for case 2. For SWBPNTs and SWGaPNTs with diameter D n ≳ 1.0 nm, the difference between the  In contrast to the reasonable agreement of the E s values for SWInPNTs assessed by Jiang and Guo [27] and those obtained for case 2 of the input parameters (see Figure 13c), less agreement is observed when the E s results reported in the same work [27] for SWBPNTs and SWGaPNTs are compared with the current results calculated for case 2. For SWBPNTs and SWGaPNTs with diameter D n 1.0 nm, the difference between the E s values by Jiang and Guo [27] and the current ones reaches ≈33% and ≈23%, respectively (see Table 8). These dissimilarities can be attributed to different calculation methods for the surface Young's modulus and force field constants. Jiang and Guo [27] assessed E s employing closed-form analytical solutions based on the "stick-and-spring" model. To calculate the bond stretching, k r , and bond bending, k θ , force constants for phosphide NTs, Jiang and Guo [27] modified the UFF method and used bond length values not equal to the present study. Moreover, they introduced a negative inversion force constant, k ψ , without taking into account the dihedral torsion force constant, k φ . With regard to the results reported by Kochaev [26], the best agreement is observed when comparing with the current E s values for SWGaPNTs. The difference between the values of E s for case 1 and the maximum E s values evaluated by Kochaev [26] reaches ≈34% and ≈23% for (n, n) and (n, 0) GaP nanotubes, respectively (see Table 8).
To clarify the results shown in Figure 13, the surface Young's modulus of the phosphide NTs is plotted as a function of the bond length, a A13−P , for cases 1 and 2 of input parameters, in Figure 14a. The values of E s evaluated by Jiang and Guo [27] are also shown in Figure 14a. The evolution of the ratio of the surface Young's moduli obtained for cases 1 and 2, E s(UFF) /E s(DFT) , with a A13−P is presented in Figure 14b.  [27] and the current ones reaches ≈33% and ≈23%, respectively (see Table 8). These dissimilarities can be attributed to different calculation methods for the surface Young's modulus and force field constants. Jiang and Guo [27] assessed E s employing closed-form analytical solutions based on the "stick-and-spring" model. To calculate the bond stretching, k r , and bond bending, k θ , force constants for phosphide NTs, Jiang and Guo [27] modified the UFF method and used bond length values not equal to the present study. Moreover, they introduced a negative inversion force constant, k ψ , without taking into account the dihedral torsion force constant, k ϕ . With regard to the results reported by Kochaev [26], the best agreement is observed when comparing with the current E s values for SWGaPNTs. The difference between the values of E s for case 1 and the maximum E s values evaluated by Kochaev [26] reaches ≈34% and ≈23% for (n, n) and (n, 0) GaP nanotubes, respectively (see Table 8).
To clarify the results shown in Figure 13, the surface Young´s modulus of the phosphide NTs is plotted as a function of the bond length, a A13-P , for cases 1 and 2 of input parameters, in Figure 14a. The values of E s evaluated by Jiang and Guo [27] are also shown in Figure 14a. The evolution of the ratio of the surface Young´s moduli obtained for cases 1 and 2, E s(UFF) E s(DFT) ⁄ , with a A13-P is presented in Figure 14b. For phosphide nanotubes, E s decreases with increasing bond length, a A13-P , i.e., from SWBPNTs to SWInPNTs, except for case 1, for which the values of E s calculated for SWAlPNTs and SWInPNTs are nearly equal. Furthermore, the decreasing trend of the surface Young´s modulus with a A13-P can be easily established from the results of Jiang and Guo [27] (see Figure 14a).
Regarding the E s results calculated for case 1 (UFF) and case 2 (DFT+MM) of the input parameters for numerical simulation, the ratio of E s(UFF) E s(DFT) ⁄ increases with the bond length of NTs, achieving the highest difference between the values of E s for the SWInPNTs. For these NTs, E s calculated for case 1 is 1.8 times larger than that for case 2. Figure 15 shows the evolutions of the surface shear modulus, G s , calculated based on Figure 12 and Equation (39), as a function of the nanotube diameter, D n , for SWBPNTs, SWAlPNTs, SWGaPNTs and SWInPNTs. For phosphide nanotubes, E s decreases with increasing bond length, a A13−P , i.e., from SWBPNTs to SWInPNTs, except for case 1, for which the values of E s calculated for SWAlP-NTs and SWInPNTs are nearly equal. Furthermore, the decreasing trend of the surface Young's modulus with a A13−P can be easily established from the results of Jiang and Guo [27] (see Figure 14a).

Surface Shear Modulus of Phosphide Nanotubes
Regarding the E s results calculated for case 1 (UFF) and case 2 (DFT+MM) of the input parameters for numerical simulation, the ratio of E s(UFF) /E s(DFT) increases with the bond length of NTs, achieving the highest difference between the values of E s for the SWInPNTs. For these NTs, E s calculated for case 1 is 1.8 times larger than that for case 2. Figure 15 shows the evolutions of the surface shear modulus, G s , calculated based on Figure 12 and Equation (39), as a function of the nanotube diameter, D n , for SWBPNTs, SWAlPNTs, SWGaPNTs and SWInPNTs. For small NT diameters, D n , the surface shear modulus, G s , decreases for zigzag (n, 0) NTs ( Figure 15a) and increases for armchair (n, n) NTs (Figure 15b), being stable for high values of D n , regardless of the NTs symmetry group. The surface shear modulus for chiral (n, m) NTs is nearly constant over the entire range of NT diameters and is equal to the converged average value of G s , obtained for (n, n) and (n, 0) NTs (Figure 15a,b). For case 1, the SWBPNTs have the highest convergent average value of the surface shear modulus (G s = 0.108 TPa⋅nm) For small NT diameters, D n , the surface shear modulus, G s , decreases for zigzag (n, 0) NTs ( Figure 15a) and increases for armchair (n, n) NTs (Figure 15b), being stable for high values of D n , regardless of the NTs symmetry group. The surface shear modulus for chiral (n, m) NTs is nearly constant over the entire range of NT diameters and is equal to the converged average value of G s , obtained for (n, n) and (n, 0) NTs (Figure 15a,b).

Surface Shear Modulus of Phosphide Nanotubes
For case 1, the SWBPNTs have the highest convergent average value of the surface shear modulus (G s = 0.108 TPa·nm) and G s decreases approximately 2.4 times for SWGaPNTs (G s = 0.048 TPa·nm), SWAlPNTs (G s = 0.041 TPa·nm) and SWInPNTs (G s = 0.044 TPa·nm). For case 2, the converged average value of G s decreases from SWBPNTs (G s = 0.083 TPa·nm) to SWGaPNTs (G s = 0.034 TPa·nm) and also for SWInPNTs, which have the lowest G s value of 0.023 TPa·nm. The value of G s , of the three symmetry groups of phosphide NTs, for small diameters nanotubes depends on the chiral angle and decreases from zigzag (n, 0) NTs with Θ = 0 • to chiral (n, m) NTs with Θ = 19.1 • , and then to armchair (n, n) NTs with Θ = 30 • (Figure 15c-f). It can be noted in Figure 15c-f that the greater the value of the bond length, a A13−P ,of the phosphide nanotube (see Table 2), the greater the value of the NT diameter, D st n , for which the shear modulus becomes nearly constant, regardless of the NTs symmetry. The values of D st n are the same regardless of the input parameters case, although the values of G s calculated for case 1 are always higher than those for case 2.
In order to understand better the results shown in Figure 15, the convergent average value of the surface shear modulus, G s , the NT diameter, D st n , for which the value of G s becomes stable, and the ratio between the surface shear moduli calculated for cases 1 and 2, G s(UFF) /G s(DFT) , were plotted in Figure 16 as a function of the NT bond length, a A13−P .  ) and also for SWInPNTs, which have the lowest G s value of 0.023 TPa⋅nm. The value of G s , of the three symmetry groups of phosphide NTs, for small diameters nanotubes depends on the chiral angle and decreases from zigzag (n, 0) NTs with Θ = 0° to chiral (n, m) NTs with Θ = 19.1°, and then to armchair (n, n) NTs with Θ = 30° (Figure 15c-f). It can be noted in Figure 15c-f that the greater the value of the bond length, a A13-P , of the phosphide nanotube (see Table 2), the greater the value of the NT diameter, D n st , for which the shear modulus becomes nearly constant, regardless of the NTs symmetry. The values of D n st are the same regardless of the input parameters case, although the values of G s calculated for case 1 are always higher than those for case 2.
In order to understand better the results shown in Figure 15, the convergent average value of the surface shear modulus, G s , the NT diameter, D n st , for which the value of G s becomes stable, and the ratio between the surface shear moduli calculated for cases 1 and 2, G s(UFF) G s(DFT) ⁄ , were plotted in Figure 16 as a function of the NT bond length, a A13-P . As can be seen from Figure 16a, the decrease in the value of the surface shear modulus calculated for case 1 occurs when the bond length increases from 0.183 nm (SWBPNTs) to 0.225 nm (SWGaPNTs). Nearby are the SWAlPNTs and SWInPNTs nanotubes, for which the values of G s are close to each other and to those evaluated for SWGaPNTs. For As can be seen from Figure 16a, the decrease in the value of the surface shear modulus calculated for case 1 occurs when the bond length increases from 0.183 nm (SWBPNTs) to 0.225 nm (SWGaPNTs). Nearby are the SWAlPNTs and SWInPNTs nanotubes, for which the values of G s are close to each other and to those evaluated for SWGaPNTs. For case 2, the G s value decreases with increasing of a A13−P . On the contrary, the diameter, D st n , increases with increasing the bond length (see Figure 16b). With regard to the ratio G s(UFF) /G s(DFT) , the values of G s obtained for case 1 (UFF) is 1.3, 1.4 and 1.9 times bigger than those calculated for case 2 (DFT + MM) for SWBPNTs, SWGaPNTs and SWInPNTs, respectively (see Figure 16c).
It is worth noting that, similarly to the surface Young's modulus, the shear modulus of phosphide NTs is also sensitive to the input parameters for numerical simulation, which in turn depend on the bond length and the force field constants.
Since the phosphide NTs have a potential application in the NT-based devices and hybrid nanostructures, where they can be combined with carbon or other non-carbon NTs, the surface Young's and shear moduli of the SWBPNTs, SWAlPNTs, SWGaPNTs and SWInPNTs were compared with respective surface elastic moduli of the SWCNTs and SWBNNTs in Table 9. E s and G s of carbon and boron-nitride NTs were calculated from the results of the authors' previous work [39], using the expressions E s = Et n and G s = Gt n , respectively, and taking into account that the value of the nanotube wall thickness t n = 0.34 nm, for both classes of NTs. It can be concluded from Table 9 that SWBPNTs and, particularly, SWAlPNTs, SWGaP-NTs and SWInPNTs have weak mechanical properties when compared with SWCNTs and SWBNNTs. Thus, when designing and constructing NT-based devices and hybrid nanostructures, it is desirable to combine the phosphide NTs, with low mechanical strength, with CNTs or N-CNTs, with high mechanical strength.

Poisson's Ratio of SWBPNTs, SWAlPNTs, SWGaPNTs and SWInPNTs
The Poisson's ratio of the phosphide NTs can be calculated, assuming the isotropy condition and considering Equations (36) and (37), by the following expression: where EI and GJ are bending and torsional rigidities, respectively. The EI and GJ rigidities can be calculated either from the results of FE analysis by Equations (31) and (32), respectively, or analytically using Equations (34) and (35), respectively. Combining Equation (40) with relationships (34) and (35), the Poisson's ratio can be expressed by an equation independent of the NT diameter, through the fitting parameters β A13P and γ A13P , as follows: Figure 17 shows the evolution of the Poisson's ratio, ν, calculated by Equation (40), with the NT diameter, D n , for SWBPNTs, SWAlPNTs, SWGaPNTs and SWInPNTs in Table 5. The two cases of the analyzed input parameters are considered. The values of ν calculated by Equation (41), which is independent of D n , are also presented in Figure 17. The Poisson's ratio evaluated for case 2 is higher than that for case 1. For (n, n), (n, 0) and (n, m) phosphide NTs with high diameters, the Poisson's ratio converges to the constant value calculated by Equation (41). The greater the value of the bond length, a A13−P , the greater the value of the nanotube diameter, D st n , for which ν becomes stable (see Figure 17). The value of D st n does not depend on the case of input parameters used for numerical simulations. greater the value of the nanotube diameter, D n st , for which ν becomes stable (see Figure   17). The value of D n st does not depend on the case of input parameters used for numerical simulations. For phosphide NT diameter lower than D n st , the Poisson's ratio decreases for (n, n) armchair and (n, m) chiral NTs, but for (n, 0) zigzag NTs, the ν value increases. In addition, the (n, 0) phosphide NTs with small diameters D n < 1.0 nm have a negative Poisson's ratio (demonstrate an auxetic behavior). It can also be seen from Figure 17 that for NTs with a diameter smaller than the diameter D n st , the Poisson's ratio noticeably depends on the chiral angle and is greater for (n, n) NTs (Θ = 30°) and smaller for (n, 0) NTs (Θ = 0°). This effect is particularly evident for the small phosphide NT diameters, D n ≤ 1.0 nm. A similar dependence of the value of ν on the NT chiral angle was reported for the SWBNNTs with diameter below 1.5 nm [39]. This result was explained by the fact that the ratio between bending and torsion rigidities, EI/GJ, did not have a constant value for D n < 1.5 nm.
To analyze the influence of the input parameters on the Poisson's ratio results, the evolutions of the Poisson´s ratio, ν, as a function of the NT diameter, D n , for cases 1 and 2 of (n, 0) and (n, n) SWBPNTs, SWAlPNTs, SWGaPNTs and SWInPNTs are plotted in Figure 18. The available literature results for (n, n) SWBPNTs, SWGaPNTs and SWInPNTs are also shown ( Figure 18e). As can be seen in Figure 18a,b, for case 1 of the phosphide NTs, the convergent average value of ν increases from SWBPNTs (ν = 0.01) to SWAlPNTs (ν = 0.14), and then decreases for SWInPNTs (ν = 0.09). For case 2, the Poisson´s ratio grows For phosphide NT diameter lower than D st n , the Poisson's ratio decreases for (n, n) armchair and (n, m) chiral NTs, but for (n, 0) zigzag NTs, the ν value increases. In addition, the (n, 0) phosphide NTs with small diameters D n < 1.0 nm have a negative Poisson's ratio (demonstrate an auxetic behavior). It can also be seen from Figure 17 that for NTs with a diameter smaller than the diameter D st n , the Poisson's ratio noticeably depends on the chiral angle and is greater for (n, n) NTs (Θ = 30 • ) and smaller for (n, 0) NTs (Θ = 0 • ). This effect is particularly evident for the small phosphide NT diameters, D n ≤ 1.0 nm. A similar dependence of the value of ν on the NT chiral angle was reported for the SWBNNTs with diameter below 1.5 nm [39]. This result was explained by the fact that the ratio between bending and torsion rigidities, EI/GJ, did not have a constant value for D n < 1.5 nm.
To analyze the influence of the input parameters on the Poisson's ratio results, the evolutions of the Poisson's ratio, ν, as a function of the NT diameter, D n , for cases 1 and 2 of (n, 0) and (n, n) SWBPNTs, SWAlPNTs, SWGaPNTs and SWInPNTs are plotted in Figure 18. The available literature results for (n, n) SWBPNTs, SWGaPNTs and SWInPNTs are also shown ( Figure 18e). As can be seen in Figure 18a,b, for case 1 of the phosphide NTs, the convergent average value of ν increases from SWBPNTs (ν = 0.01) to SWAlPNTs (ν = 0.14), and then decreases for SWInPNTs (ν = 0.09). For case 2, the Poisson's ratio grows from 0.06 for SWBPNTs to 0.20 for SWInPNTs (see Figure 18c,d). The values of ν obtained for cases 1 and 2 of SWGaPNTs are close and equal to 0.12 and 0.13, respectively.  Although the Poisson's ratio values evaluated by Jiang and Guo [27] for SWBPNTs, SWGaPNTs and SWInPNTs are 83%, 69% and 56% higher, respectively, than those calculated in the present study for case 2, it is possible to compare the trends of the evolutions of ν as a function of NT diameter, D n , as shown in Figure 18d. Jiang and Guo [27] found that for (n, n) and (n, 0) phosphide NTs, the Poisson's ratio decreases with increasing D n , and then the ν value converges to an approximately constant value. This trend is in agreement with the current trend of evolution of ν with NT diameter for (n, n) phosphide NTs, although, the decreasing rate is slower for the ν evolution reported by Jiang and Guo [27].
To clarify the results shown in Figures 17 and 18, the convergent average value of the Poisson's ratio, ν, and the NT diameter, D st n , for which the value of ν becomes stable, are plotted as a function of the NT bond length, a A13−P in Figure 19. The values of ν evaluated by Jiang and Guo [27] are presented in both plots for comparison purposes.
Although the Poisson's ratio values evaluated by Jiang and Guo [27] for SWBPNTs, SWGaPNTs and SWInPNTs are 83%, 69% and 56% higher, respectively, than those calculated in the present study for case 2, it is possible to compare the trends of the evolutions of ν as a function of NT diameter, D n , as shown in Figure 18d. Jiang and Guo [27] found that for (n, n) and (n, 0) phosphide NTs, the Poisson's ratio decreases with increasing D n , and then the ν value converges to an approximately constant value. This trend is in agreement with the current trend of evolution of ν with NT diameter for (n, n) phosphide NTs, although, the decreasing rate is slower for the ν evolution reported by Jiang and Guo [27].
To clarify the results shown in Figures 17 and 18, the convergent average value of the Poisson´s ratio, ν, and the NT diameter, D n st , for which the value of ν becomes stable, are plotted as a function of the NT bond length, a A13-P in Figure 19. The values of ν evaluated by Jiang and Guo [27] are presented in both plots for comparison purposes. The Poisson´s ratio, ν, obtained for case 1 of phosphide NTs increases with the increasing of the a A13-P value up to 0.234 (SWAlPNTs), and then ν drops after a A13-P = 0.246 nm (SWInPNTs). For case 2, the value of ν increases with increasing bond length, which is in a good agreement with the results of Jiang and Guo [27] (see Figure 19a). Moreover, the same good agreement is found between trends of D n st values, obtained in the present work and reported by Jiang and Guo [27]. In both studies, it is shown that the greater the bond length, a A13-P , the greater the NT diameter, D n st , for which ν becomes stable (see Figure 19b). It is worth mentioning that the D n st values obtained in the present study and reported in the work of Jiang and Guo [27] are close. The current Poisson´s ratio results for (n, n) and (n, 0) phosphide NTs are summarized in Table 10 together with the values of ν available in the literature.  The Poisson's ratio, ν, obtained for case 1 of phosphide NTs increases with the increasing of the a A13−P value up to 0.234 (SWAlPNTs), and then ν drops after a A13−P = 0.246 nm (SWInPNTs). For case 2, the value of ν increases with increasing bond length, which is in a good agreement with the results of Jiang and Guo [27] (see Figure 19a). Moreover, the same good agreement is found between trends of D st n values, obtained in the present work and reported by Jiang and Guo [27]. In both studies, it is shown that the greater the bond length, a A13−P , the greater the NT diameter, D st n , for which ν becomes stable (see Figure 19b). It is worth mentioning that the D st n values obtained in the present study and reported in the work of Jiang and Guo [27] are close.
The current Poisson's ratio results for (n, n) and (n, 0) phosphide NTs are summarized in Table 10 together with the values of ν available in the literature. It can be seen from Table 10 that the Poisson's ratio results are scarce so far and the existing ν values show large scattering, regardless of the modelling and calculation approaches used for this end. Analyzing current Poisson's ratio results and those reported by Jiang and Guo [27], it can be concluded that, similar to the surface elastic moduli, the ν values are sensitive to the values of bond length and force field constants.

Conclusions
The elastic properties, comprising the three rigidities, tensile, bending and torsional, the surface Young's and shear moduli and the Poisson's ratio of SWBPNTs, SWAlPNTs, SWGaPNTs and SWInPNTs were evaluated resorting to a numerical simulation study, based on the NCM/MSM approach. The main conclusions of the present study are indicated below.
The force field constants, required for calculating the input parameters for numerical simulation, were computed for BP, AlP, GaP and InP nanostructures, using two approaches, which resulted in two input sets.
Equations describing the relationship between each of the three rigidities and the nanotube diameter were obtained for SWBPNTs, SWAlPNTs, SWGaPNTs and SWInPNTs, and the fitting parameters for the Equations (30)-(32) were calculated for two sets of input parameters for numerical simulation. This allowed to expand, to phosphide nanotubes, the previously established method, allowing the calculation of tensile, bending and torsional rigidities without resorting to numerical simulation.
A robust methodology was proposed to assess the surface Young's and shear moduli of the SWBPNTs, SWAlPNTs, SWGaPNTs and SWInPNTs, based on the Young's and shear moduli evolutions as a function of the inverse nanotube thickness. It is expected that this methodology will be useful to evaluate the surface elastic moduli of the N-CNTs, for which the exact value of the nanotube wall thickness is unknown.
The tensile, bending and torsional rigidities, the surface Young's and shear moduli, and the Poisson's ratio of SWBPNTs, SWAlPNTs, SWGaPNTs and SWInPNTs are sensitive to the bond length and force field constants of the diatomic phosphide nanostructures.
The results obtained provide a substantial contribution to a benchmark with regard to the determination of the elastic properties of the phosphide nanotubes by numerical and analytical methods. Funding: This research is sponsored by FEDER funds through the program COMPETE-Programa Operacional Factores de Competitividade-and by national funds through FCT, Fundação para a Ciência e a Tecnologia, under the project UIDB/00285/2020 and LA/P/0112/2020.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author after obtaining permission of an authorized person.