Atomistic modeling of electromechanical properties of piezoelectric zinc oxide nanowires

Currently, numerous articles are devoted to examining the influence of geometry and charge distribution on the mechanical properties and structural stability of piezoelectric nanowires (NWs). The varied modeling techniques adopted in earlier molecular dynamics (MD) works dictated the outcome of the different efforts. In this article, comprehensive MD studies are conducted to determine the influence of varied interatomic potentials (partially charged rigid ion model, [PCRIM] ReaxFF, charged optimized many-body [COMB], and Buckingham), geometrical parameters (cross-section geometry, wire diameter, and length), and charge distribution (uniform full charges versus partially charged surface atoms) on the resulting mechanical properties and structural stability of zinc oxide (ZnO) NWs. Our optimized parameters for the Buckingham interatomic potential are in good agreement with the existing experimental results. Furthermore, we found that the incorrect selection of interatomic potentials could lead to excessive overestimate (61%) of the elastic modulus of the NW. While NW length was found to dictate the strain distribution along the wire, impacting its predicted properties, the cross-section shape did not play a major role. Assigning uniform charges for both the core and surface atoms of ZnO NWs leads to a drastic decrease in fracture properties.


Introduction
Owing to their superior mechanical strength, chemical and thermal stability, and piezoelectric properties, zinc oxide (ZnO) nanowires (NWs) have emerged as important nanostructure components that have an ample range of nanodevice applications, such as nanogenerators, nanotransistors, nanodiodes, nanosensors, and supercapacitors [1][2][3][4][5].The working mechanism of such piezoelectric nanodevices is largely influenced by the piezoelectric properties of ZnO NWs.Thus, designing a tunable mechanism based on ZnO NWs requires a reliable modeling technique that accounts for the surface effects of the actuating and sensing behaviors of a device.
The electromechanical properties of ZnO NWs are significantly influenced by defects [4,6], NW diameter [7], cross-sectional geometry [8][9][10], polarization orientation [11], and growth direction [11,12].This dependence was attributed to the surface effect and structural anisotropy at the nanoscale level [13] and to the associated reconfiguration of surface atoms that led to a stiffer crystal, resulting in a stronger NW than its bulk form [14,15].For example, Agrawal et al [16] performed a uniaxial tension test of ZnO NWs and found that their respective fracture stress and failure strain can be as high as 9.53 GPa and 6.2%, respectively, which are 500% larger than those reported for bulk ZnO.They also found that the elastic modulus decreased by 12.5% when the NW diameter increased from 20 nm to 80 nm, and NWs with larger diameters exhibited properties comparable to bulk material.
Atomistic simulations offer a reliable approach that complements experimental measurements.They are capable of treating all the parameters that govern the coupled behavior of piezoelectric NWs.The mechanical and piezoelectric properties of the ZnO NWs were determined using density functional theory calculations and molecular dynamics (MD) [8,[10][11][12][13][14][17][18][19][20][21][22].However, density functional theory investigations were mainly limited to NWs with very small diameters (D NW  5 nm) and short lengths (a single unit cell) due to the computational burden required to model larger systems [11].This order of magnitude difference in the dimensions of NWs from existing experimental works [1,15,16,22] limited our ability to understand the deformation mechanism of NWs used in engineering applications.On the other hand, MD can simulate much larger NWs, which allows us to tailor their properties and design novel nanoscale systems that utilize their multifunctionality [23].
The structural properties of ZnO NWs have been systematically investigated using MD simulations based on various interatomic potentials, such as Buckingham potential [14], partially charged rigid ion model (PCRIM) potential [24], among others.For example, MD investigations performed by Bandura et al [11] showed that the elastic modulus of NWs with a 20 nm thickness is around 161 GPa, which is close to the experimental findings.The predicted axial stiffness constant C 11 ( ) of the NWs using MD simulations was around 334 GPa [12], which is more than double the experimentally measured 140 GPa of bulk ZnO [25].He et al [17] concluded that the stress-strain curve of a ZnO NW subjected to a tensile load consists of the following four deformation stages: (i) initial elastic response of the Wurtzite structure, (ii) phase transformation into a body-centered tetragonal (BCT) structure, (iii) elongation of the BCT structure, and (iv) brittle fracture.Wang et al [21] reported 8% and 18% reductions in the elastic modulus and yield stress of ZnO NWs, respectively, when the temperature increased from 0.1 to 499 K.They also reported that the elastic properties of ZnO NWs were not a function of the strain rate or the length of the NW.Furthermore, several core-shell models have been proposed to describe the mechanical and piezoelectric properties of ZnO NWs [25,26].Fan et al [26] considered the NW crosssection to consist of three layers: the outer shell formed from the outmost three atomic layers, the inner shell, and core atoms.Additionally, some atomistic studies investigated the dependency of NW properties on its diameter and the conclusions varied widely [18,19,21,22,25].Zhang et al [8] attributed the size dependence of Young's moduli to the saturation of surface bonds.The free boundaries and the quantum confinement effect in the NWs also result in larger piezoelectric constants and band gaps than those measured for bulk ZnO [8].The dependency of the NW properties on its diameter is attributed to the combined effect of surface reconstruction and long-range ionic interactions.Furthermore, the phase transformation of NW crystals after reaching their maximum stress during tensile simulation has been reported in several MD studies [16,17,19,21,27].This Wurtzite-BCT transition has not been observed in any experimental investigation and was attributed to the inability to control the loading speed in experimental measurements or to the failure of the interatomic potential to capture the actual behavior of the NW at high strains [16].Moreover, ZnO NWs with circular, hexagonal, rectangular, square, and hexagonal crosssectional geometries have been investigated in previous MD studies to study the effect of NW shape on its performance and stability [10,11,18,20,21,24].
Despite the advantages of MD simulations over other analytical and atomistic modeling techniques, there are several limitations to the reported studies.For example, modeling a ZnO NW with a diameter of 50 nm and a length of 1 μm requires huge computational resources [21].Such NW consists of 164 million atoms, and very small strain rates are essential in guaranteeing a uniform deformation along its length.As a result, previous studies considered a wide range of NW lengths and assumed that the obtained properties were independent of the length parameter.Furthermore, it has also been concluded by several researchers [14,17,21] that the strain rate does not have a role in dictating the results, and strain rates up to 10 9 1 s −1 can be used to obtain accurate results.Moreover, Wang et al [21] reported that the elastic properties of ZnO NWs were independent of the strain rate and length of the NW.However, these general conclusions are based on NWs with smaller lengths (<10 nm) and will thus result in unrealistic behavior for longer NWs [14,28].
Additionally, modeling the surface and core atoms of ZnO NWs varied significantly in the literature, with some works [29,30] that assign (i) uniform charges for all atoms and others [31][32][33], (ii) partial charges for surface atoms, and (iii) full charges for core atoms.However, the charge distribution affects the structural stability of the NWs and leads to radically varying results.Finally, the results of the simulations performed using different interatomic potentials contradicted some cases due to the inherent limitations and assumptions of potentials.For example, the phase transformation in the ZnO crystals during stretching reported in several MD studies was not observed experimentally [16,17,19,21,27].Moreover, simple pairwise potentials such as, the Buckingham type potential, cannot accurately capture the breaking and formation of atomic bonds that occur during the fracture process of ZnO NWs [16].
The purpose of the current work is to identify and elucidate the reasons for the discrepancies and contradictions that exist in the literature concerning the mechanical behavior of ZnO NWs.The novelty of the current work lies in studying the impact of MD modeling parameters on the mechanical behavior of ZnO NWs.This is carried out by conducting comprehensive MD simulations considering different interatomic potentials, NW geometry, and shape, as well as the charge distribution upon the elastic and fracture behavior of ZnO NWs.First, we investigated the impact of the interatomic potentials on the electromechanical response of NWs using PCRIM, ReaxFF, charged optimized many-body (COMB), and Buckingham potentials.The effect of charge distribution on the properties of NWs and their structural stability was studied by considering the full and partial charges of the surface atoms.Furthermore, we examined the hypothesis offered by some researchers (e.g.[16] and [17]) regarding the phase transition during the elongation of the NWs.Special attention was given to the influence of NW geometry, interatomic potential type, temperature, and applied strain rate.The effect of NW length on the strain distribution along its length was studied by considering lengths ranging from 25 to 400 Å.The induced uniaxial stress at different sections of the NW was obtained, and the dependence of the elastic modulus was investigated.Finally, we determined the appropriate strain rate to be applied to the NW and found a relationship between the strain rate and NW length.

Molecular dynamics modeling of ZnO NWs
Ab initio methods are widely used for benchmark predictions of piezoelectric constants; however, for larger systems (>5 nm) these methods are incompatible due to the associated computational labor [11].Moreover, with ab initio, thermal effects in nanomaterials, such as pyroelectricity, cannot be studied.In this framework, classical MD simulations can be effectively utilized for large-scale atomic structures, accounting for temperature dependencies.MD simulations are preferred over nanoscale experiments for describing the morphology and understanding defect formation or the initiation of nanostructures in detail.In this section, we present the modeling procedure, atomic structure parameters, details of interatomic potentials, and all other parameters considered in this work.The atomic simulations of ZnO NWs with orientation [0001] were performed with LAMMPS [34] using the SHARCNET supercomputing platform.The following interatomic potentials were used in the LAMMPS code: Buckingham, PCRIM, COMB, and ReaxFF.The velocity Verlet algorithm was used for to integrate the equations of motion.The visual MD software and the open visualization tool OVITO were utilized to visualize the outcomes of the MD simulations.

Atomistic models of ZnO NWs
In this work, we conducted uniaxial tensile MD simulations of ZnO NWs, which have a growth direction along the [0001] axis.A wurtzite-structured crystal of space group P63Mc was used with a = 3.2417 Å, c/a = 1.6003Å, and u = 0.3819 lattice parameters, [35] to build the atomic structures of ZnO NWs.An in-house MATLAB code was developed to generate the atomic structures of ZnO NWs with hexagonal, triangle, and square cross-sectional areas, as shown in figures 1(b)-(d).The length and diameter effects on the elastic and fracture properties were investigated by modeling ZnO NWs with lengths and diameters ranging from ∼2.5 nm to 42 nm and 1 nm to 30 nm, respectively.NWs with smaller diameters were not considered in this work due to their structural instability and transformation into non-piezoelectric material with a d-BCT structure [20].
In the current work, MD simulations were performed to determine the elastic constant of ZnO NWs in the axial direction only, as they can be considered quasi-one-dimensional nanostructures.The size of the periodic box was 14.45 nm [20] to avoid the image effect.The popular velocity Verlet algorithm was used with a time step of 1 fs for numerical integration of the equations of motion.The geometric, loading, and MD parameters considered in the current study are summarized in table 1.

Tensile simulations of ZnO NWs
2.2.1.Equilibrium stage.Prior to the loading condition, the structure was first relaxed using the annealing process with a time step of 1 fs at room temperature (300 K).Utilizing the Maxwell-Boltzmann distribution, the initial velocities were randomly assigned to the atoms.The initial system of the ZnO NWs was equilibrated for 50 ps without applying any external force field to reach the equilibrium state.Next, under constant atmosphere pressure (NPT), a NW was relaxed for 1 ns to release the residual stresses and for another 1 ns under constant temperature and volume (NVT) to achieve a minimum energy state with given lattice constants and interatomic potentials [24].MD simulations were performed on WZ−ZnO with an NVT ensemble at room temperature (300 K) and zero pressure for a period of 0.5 ns.During the equilibration process, the velocities were rescaled for every 1 ps to reach a temperature of 300 K.The Nose-Hoover thermostat and barostat were applied to control the temperature and pressure, respectively.The equilibration process was continued for 20 ps with an NVT ensemble under a constant temperature of 300 K [24].
2.2.2.Defining the boundary condition.One of the ends of the ZnO NW was fixed to restrict its movement.For the tensile loading condition, the other end was subjected to constant strain loading.The middle layers of the NWs were not constrained, and the mechanical properties of NW can be predicted based on the deformation of such middle layers.The tensile loading scheme of the ZnO NW is shown in figure 2. We used periodic boundary conditions (PBCs) in the three mutually perpendicular directions, while a significant gap was maintained along the transverse direction to avoid interactions between neighboring NWs [9].

Tensile simulation.
Tensile simulations were performed to calculate the longitudinal elastic constant C 11 ( ) of ZnO NWs.In a stepwise manner, the right and left edges of the NW were separated into two rigid blocks.With every MD step, the NW underwent deformation of 1 Ǻ and then subsequently it was equilibrated for 500 ps, thus maintaining the deformation rate of 0.2 m s −1 .Note that we varied the deformation rate from 0.05 to 2 m s −1 and found no significant difference in the results.The high rate of loading may result in crystalline defects; therefore, a low strain rate of 0.001 ns −1 was applied.The new equilibrium state was used as an initial configuration for the next MD simulation.Repeating the above process, the NWs can be stretched continuously until the required strain is achieved.A strain rate of 0.001 ns −1 was retained to deform the NWs until the fracture.
2.2.4.Data capture and management.Due to the elongation of the ZnO NW, strains (ε) and corresponding stresses (σ) were developed.The atomic coordinates and forces at various strain levels were determined by averaging over 2 ps at the end of the relaxation step.The average uniaxial tensile stress on the nanowire cross-section is calculated as the total reaction axial force on all the fixed atoms divided by the NW's cross-sectional area.For the small stain rate, the deformation of ZnO NWs follows the Hooke's law, and σ-ε relationship can be obtained to predict the longitudinal elastic constant, C .
11 The strain (ε) in the NWs can be defined as where L and l are the respective lengths of the NW before and after deformation.The elastic modulus of NW can be calculated with the help of the slope of σ-ε curve.

Interatomic potential functions
The selection of interatomic potential is a crucial task of MD simulations, which represent the atom-atom interactions.In the current work, we investigated the dependence of the mechanical properties of ZnO NWs on the selection of interatomic potential.The reported results in the literature vary significantly due to the use of different interatomic potentials.Furthermore, every potential has its own limitations and strengths, including computational cost, capturing the transition of atomic structures under load, allowing interaction with the surrounding medium, and capturing bond-breaking and bond-formation processes.In particular, PCRIM, ReaxFF, COMB, and Buckingham interatomic potentials were considered in this work, and their descriptions are provided in the following subsections.

PCRIM potential.
The eight-parameter PCRIM potential was developed by Wang et al [24] in 2014.PCRIM potential can accurately reproduce the results of density functional theory and the experimental results of the mechanical properties and lattice parameters of bulk ZnO structures.Additionally, the PCRIM model was recently utilized in some studies [11,36] to calculate the elastic and fracture properties of ultrathin ZnO NWs, with diameters as low as 8.5 Å.The PCRIM potential was derived using ab initio energy surface fitting, empirical fitting, and lattice inversion.The PCRIM potential can be used to model interatomic interactions in ZnO.Note that the nature of this bond is largely ionic, and electrostatic interactions are dominant [24].The PCRIM potential consists of short-range two-body and Coulombic interactions.The short-range  ´--  interactions between Zn-Zn and Zn-O are described with the Born-Mayer potential as follows: The short-range non-coulombic interactions between O-O atoms are described with the pairwise Morse potential: where r ij is the bond length, while A , ij C , ij and ij r are the parameters that can be optimized for every type of pairwise interactions.The parameters used in the PCRIM potential for ZnO structures are summarized in table 2.
The long-range coulombic interactions can be calculated using the following equation: where q i and q j are the charges on the ions and 0 e denotes permittivity of free space.The long-range Coulombic interactions in the case of Zn-Zn were modeled using an effective charge of ±1.14 e on cations and anions.For short-range non-Coulombic interactions, a cut-off radius was taken as 12.0 Å, whereas the Ewald method was used for long-range Coulombic interactions [14,37].
2.3.2.ReaxFF potential.Raymand et al [38] developed a reactive force field (ReaxFF) to model ZnO structures and accurately determined their mechanical properties.The developed force field is capable of modeling the breakage and formation of atomic bonds that occur at high deformation in ZnO NWs.Hence, enabling the modeling of chemical reactions between the NWs and the surrounding medium, including water molecules, was considered in an extension of the ReaxFF force field.The total potential energy of the atomic system is given as follows: Coul In equation (4), the first five terms are bond-order dependent terms, which contribute to the local environment of every atom.E vdW and E Coul are the non-bond order dependent terms, and such terms are given by a taper function and shielded to avoid excessive repulsion between the pair of atoms at short distances.Details of each term can be found in [12].Overall, ReaxFF provides a good representation of the partial atomic charges of zinc and oxygen.Mulliken charges (B3LYP) and EEM charges (ReaxFF) for bulk ZnO are listed in table 3.

COMB potential.
The second-generation COMB potential is a variable charge-type potential that includes Coulombic interactions treated by the Wolf summation approach [39], correction terms, and the self-energy function.The total energy (E T ) of a system of atoms is given as follows: Coul and E ij polar are the self-energy of atom i, the bond-order potential between i and j atoms, the Coulomb interactions energy, the polarization terms, respectively, for organic systems used for COMB 3. E barr is a charge barrier function, E vdW is the van der Waals energy for COMB 3, and E corr is angular correction terms.In LAMMPS, this can be done by using the fix qeq/comb command in the input script while using the secondgeneration COMB potential.At each loading step, the  charge equilibration scheme allows the atoms to dynamically vary the charge such that atomic interactions are accurately computed.In the loading process, atoms are adjusted to relax their distance of separation between the atoms.
2.3.4.Buckingham potential.Buckingham [40] modified the potential developed by Lennard-Jones [41] for helium, neon and argon, and over the period of time, researchers modified Buckingham potential to account different issues such addition of coulombic term in it, improving its performance in short range interactions [42] etc.The modified Buckingham potential has been reported to accurately model the elastic, lattice, and dielectric constants of bulk ZnO [9,12,14,16,17,37].Therefore, we used the Binks parameters with the Buckingham potential with a cut-off set to 10 Å [17] for the short-range atomic interactions, t.The interaction between the two atoms in the Buckingham potential is as follows: where A, ρ, and C are potential parameters for short-range interactions.The values of these parameters are taken from [20] and summarized in table 4. The total energy of the Buckingham potential is given as follows: where E long is the long-range Coulombic interaction.For Coulomb potential with long-range interactions the Ewald and Wolf cut-off function parameters can be used.The Wolf Coulomb potential with long-range interactions can be calculated using the following relation [37]: where R c is the cut-off radius and a is the damping coefficient.In the present work, we considered, R c = 7 Å and a = 0.4 [37].

Cut-off function parameters
In the present work, we used two cut-off function parameters: Wolf and Ewald.These are important parameters that can be used with core-shell potentials when studying the piezoelectric properties of ZnO NWs.To calculate the surface properties of ZnO NWs, long-range Coulombic interactions play an important role.Hence, special attention is required when calculating long-ranged Coulombic interactions.For NWs with D < 10 nm, the Ewald summation method is a powerful tool for predicting long-range Coulombic force [14,37].The periodicity of the boundary conditions is required for the Ewald solver.The Ewald solver is computationally promising and effectively convergent, which makes the simulation stable and does not rely on the chosen integration time step [43].However, the main drawback of the Ewald summation solver is that it does not model the surface effect because it relies on lattice periodicity; hence, the Wolf summation method can be utilized to calculate the long-ranged Coulombic interactions for ZnO NWs.

Charges of surface atoms
It was observed that the Zn and O atoms of the polar (0001) surfaces of ZnO were unstable.Thus, the increasing size of NWs leads to a non-convergent electrostatic potential.To correct this non-convergent potential, surface reconstructions or charge reduction methods can be used, where the Zn and O atoms of polar surfaces are considered reduced charges [24,25].The partial charges on the surface Zn and O atoms were set to +1.5e * and −1.5e * , respectively, where e * = 0.58e is the effective electron charge, considering the fractional iconicity.We used the Buckingham-type potential to reproduce the atomic interactions in the ZnO NW, which predicts the experimentally observed surface reconstruction of the polar ZnO (0001) surface.

Results and discussions
This section is divided into four parts, and each section examines the effect of important parameters on the mechanical properties of ZnO NWs.The first is concerned with the effect of interatomic potentials, the second with the effect of an NW's geometry and shape, the third with the effect of surface charge distribution, and the fourth deals with the strain rate upon the mechanical properties of the NWs.Using the MD approach of the ZnO NW described in section 2, we first present a comparison of our predictions of Young's moduli of different NWs with those of the available MD and experimental data in the literature in table 5.It may be observed that the current predictions are in good agreement with the available MD results for the different diameters of NWs.The slight deviation in the values of Young's moduli is attributed to the use of different lengths of ZnO NW and all other selected parameters in the MD modeling and quantum calculations.Moreover, our results are in good agreement with the reported experimental results for a case of a 100 Å diameter NW.In the case of a lower diameter NW (50 Å), the current results underestimate the value of the elastic modulus by 19%.This discrepancy can be attributed to the inherent limitation in experimental techniques at the nanoscale level, as well as the difficulty in controlling the test of smaller diameter NWs and interpretating the results [16,37,46].Namely, errors and uncertainties in the calculation of NW diameter and cross-section may result in a bias toward higher elastic modulus values at smaller diameters resulting, in an overemphasized size effect [46].Note that very few experimental studies [11,37] reported the elastic moduli of smaller diameter NWs (<200 Å), and, in fact, few studies [11,46] defined the limit of NW diameter as 200 Å to use the MD or experimental technique.This is attributed to the fact that the minimum diameter NW of experimentally synthesized NWs was found to be ∼200 Å [11].

Effect of interatomic potentials
In this section, three different interatomic potentials, namely Buckingham, PCRIM, and COMB, were used to determine the mechanical properties of piezoelectric nanowires.The purpose of our selection was based on the fact that these potentials were used in the literature without providing bounds to their validity and accuracy.Figure 3 shows an equilibrated ZnO NW and an extended ZnO NW under uniaxial force modeled with the COMB interatomic potential.The charge distribution of the O and Zn atoms in NW with a diameter ∼30 Å (average charge on O is −0.579 e and Zn is 0.579 e) are shown in figure 4 for the three potentials under consideration.The figures show a clear distinction between the charge distributions associated with each of the three interatomic potentials considered.It should be noted that in the case of the Buckingham potential, all internal atoms possess a constant charge irrespective of their distances from the surface.
Unlike the Buckingham and PCRIM potentials, the COMB 3 potential systematically reassigns the value of the atomic charges based on the spatial locations and surroundings of the atoms.The difference in the way the atomic charges are treated in these interatomic potentials governs the behavior of the material.In fact, this is depicted in figure 5, which shows the relationship between the uniaxial tensile stress and the uniaxial tensile strain for ZnO NWs of a diameter of ∼30 Å and length of ∼50 Å.Three observations can be made based on the obtained results: (i) PCRIM potential predicts an upper bound, while Buckingham potential predicts a lower bound of the elastic properties, (ii) COMB 3 predicts unphysical strain hardening in the material behavior beyond 4.5% strain, and (iii) the material fractures rapidly following the limit of the elastic loading, indicative of the brittle nature of the NWs.Table 6 summarizes the predicted mechanical properties of NWs for the three interatomic potentials.In view of its popularity among the research community and our ability to facilitate the validation of some aspects of our work, we elected to use Buckingham interatomic potential in the remainder of the study, except for the case when we examined the effect of NW diameter on mechanical properties.

Effect of geometry
In this section, we examine the influence of the main geometrical parameters, such as length, diameter, and cross-sectional shape, on the resulting mechanical properties of NWs.  Figure 6 shows the stress-strain relations for NWs of the same diameter (30 Å) but with varying lengths.This shows that the elastic behavior of all lengths coincides with the fracture stress.This figure also shows that the fracture stress of the NW increases as the length of the NW decreases.Analyzing the stress strain curves shows decaying in the elastic modulus with increasing length, reaching an asymptotic value of 150 Å.This difference is mostly caused by the boundary effect, and care must be taken in future studies to avoid anomalous results.Figures 7(a) and (b) demonstrate the variation of fracture strain and fracture stress concerning the length of the NW with a constant diameter of 30 Å.As stated earlier, it can be seen that both the fracture strain and fracture stress of the NW degrade as the length of the NW increases.A steep fall in fracture strain and fracture stress is observed when the length of the NW changes from 20 to 100 Å; beyond this length, it decreases at a comparatively smaller rate.Figures S1 and S2 in the supplementary file illustrate the snapshots of ZnO NWs of diameter 70 Å and length of 25 Å and 420 Å with the undeformed equilibrated structure of the NWs, NW just before the fracture, and NW structure after failure.It can be seen that the fracture strain of the NW degrades with the increasing length of the NW.Additionally, the figure depicts, upon applying axial strain a notch formation before the fracture occurs.
We also examined the effect of the shape of the NW cross-section on its mechanical properties.Figure S3 in the supplementary file shows the three selected shapes: hexagon, triangle, and square.The mechanical properties corresponding to these geometries are shown in figure 8.This shows a very slight change in properties, effectively concluding that shape does not play a major role in influencing the mechanical properties of the wire.The elastic properties and percentages of surface atoms for the three shapes considered are summarized in table 7. Five NW diameters were selected (10-94 Å) to study their effect upon the atomistic mechanical properties of ZnO NWs. Figure 9 demonstrates that the elastic modulus of NW increases with an increase in diameter, and that the Buckingham potential prediction is supported by [20] for the surface passivation case.Figure 10 shows the variation in fracture stress and fracture strain with diameter.The fracture strength increases with the diameter of the NW, while the converse is true for the fracture strain.

Effect of surface atom charge
In this section, we discuss the effect of surface charges on the stability and properties of NWs.In this context, figure 11 shows a snapshot of equilibrated and fractured NW with Buckingham potential with Binks constants for hexagon ZnO NW of diameter 70 Å and length 55 Å using the Wolf and Ewald cut-off functions by considering the uniform and surface charges.The structural optimization of NWs results in the essential relaxation of external (facet) atomic layers while internal (core) atoms essentially remain in regular crosssection positions.Figure 12 illustrates the effect of charge  (uniform versus surface charges) and cut-off function type (Wolf versus Ewald) on the stress-strain relationship of ZnO NW of a hexagonal cross-section with a diameter of 70 Å and a length of 55 Å. Figure 11(d) shows that the surface charges improve the mechanical properties of the NW when we use the Wolf cut-off function.This is attributed to the fact that the Wolf summation is computationally more efficient for predicting long-range Coulombic interactions, and unlike Ewald summation, PBCs are not required in Wolf summation [39,47,48].Thus, Wolf's solver can simulate the surface properties of ZnO NWs more accurately.As a result, a link between theory and experiments can be achieved.On the other hand, the inability of Ewald summation to model the surface effect is attributed to its function of lattice periodicity; hence, the Wolf summation method is more appropriate for determining the long-range Coulombic interactions for ZnO NWs.It is worth noting that in this article we used the approach developed by Fennell and Gezelter [43], which is an extension of the damped and cutoff-neutralized Coulombic sum proposed by Wolf et al [39].
The figure also depicts that with the Wolf cut-off function, the elastic behavior of NW predicted by uniform and surface charges coincide up to the fracture stress.With the Ewald cut-off function, a discrepancy between the two results is observed.This is attributed to the fact that unlike the Ewald cut-off function, the Wolf cut-off function properly accounts for the surface effects and effectively predicts long-range Coulombic interactions.Lastly, the results obtained with Wolf surface charge and fractional iconicity underestimate the elastic behavior of NW.This is because, with fractional iconicity, the NW structure was unstable upon equilibrium, as shown in figure 11(c); hence, it fractured at a much lower axial stress.Table 8 summarizes the effects of the charge and cut-off function types on the mechanical properties of the ZnO NW.
For the NWs with smaller cross-sections (lateral dimensions <20 Å), the surface-to-volume ratio is higher; thus, the surface charges have a major influence at the core of the NW, as shown in figures 13(a) and (b).The figures show the atomic charges at the surface and core of the NW of diameter 16 Å upon fracture.As seen from the figures, the charges on the atoms adjacent to the surface atoms vary significantly, thus affecting the core of the structure.However, for the NWs with higher cross-sections (lateral dimensions >20 Å), the surface-to-volume ratio is smaller, and the bulk energy accounts for the major part of the total configurational energy, as shown in figures 14(a) and (b).The figures show the variation of atomic charges in the NW of diameter 94 Å upon fracture.It can be seen that due to the larger cross-section, the influence of the outer layer of the surface vanishes, and atomic charges at the major portion of the core remain constant.Thus, the material characteristics of such NWs approach to their bulk counterparts.

Effect of strain rate
Existing investigations and earlier work by the authors [28, 49; and the references therein] indicate that the selection of improper strain rates could lead to erroneous results.The application of smaller strain rates is necessary for modeling long NWs to ensure the transmission of the applied loads through its length and to suppress the phonon drag effect in nanostructures.Accordingly, we simulated NWs under varied strain rates to determine the appropriate value of strain rate.
Observing figure 15, one can easily seen that at a strain rate of 8 ´10 −5 ns −1 , the maximum elastic stress is 12.6 GPa, while it is 7.81 GPa at a strain rate of 8 ´10 −3 ns −1 , resulting in an error of ∼38%.Our extensive research indicates that for our current study a strain rate range of 8 ´10 −5 ns −1 to 8 10 −6 ns −1 would provide reliable predictions depending on the NW's length, as shown in [28,49].

Conclusions
We conducted comprehensive MD studies to elucidate and unravel many of the contradictions and misconceptions associated with modeling piezoelectric ZnO NWs.Specifically, the influence of interatomic potentials, geometry and shape of NW, charge distribution, and strain rate on the atomistic mechanical properties of ZnO were carefully examined and assessed.Our work reveals that the accuracy of the MD predictions is greatly governed by the selection of the appropriate interatomic potential, wire geometry, charge distribution, and strain rate.The following is a summary of our findings: 1.The PCRIM potential predicts an upper bound, while the Buckingham potential predicts a lower bound of the elastic properties of ZnO NW, whereas COMB 3 predicts unphysical strain hardening in material behavior beyond the 4.5% strain.2. Among all the potentials used herein, our optimized parameters for the Buckingham interatomic potential are in good agreement with the experimental results.properly accounts for the surface effects and effectively predicts long-range Coulombic interactions, unlike the Ewald summation method.5.The cross-section shape of the NW does not play a major role in influencing the mechanical properties of the wire.
For the NWs with smaller cross-sections (lateral dimensions <20 Å), the surface charges have a major influence at the core of the NW.However, for the NWs with higher cross-sections (lateral dimensions >20 Å), the bulk energy accounts for the major part of the total configurational energy.
We feel that these efforts will help the community in its selection of these parameters.Additionally, this study provides the relationship between these parameters and the mechanical properties of ZnO NWs.Discovery Grant #RGPIN-2018-03804 and Equipment grants #RTI-2018 -00703 and RTI-2020-00687.

Figure 2 .
Figure 2. MD model of a ZnO NW with tensile loading scheme.

Figure 3 .
Figure 3. (a) Equilibrated structure of ZnO NW, and (b) NW under axial force obtained with COMB 3 potential.

Figure 5 .
Figure 5. Demonstrable differences in constitutive behavior of ZnO NW predicted by different interatomic potentials.

Figure 6 .
Figure 6.Effect of ZnO NW lengths on the stress-strain relationship of the NWs.

Figure 7 .
Figure 7. Variation of (a) fracture strain and (b) fracture stress with the length of the ZnO NW.

Figure 8 .
Figure 8. Axial stress versus axial strain relationship of three NW shapes considered.

Figure 9 .
Figure 9.Effect of nanowire diameter on the elastic modulus using Buckingham interatomic potential.

Figure 10 .
Figure 10.Effect of nanowire diameter on the mechanical properties using Buckingham interatomic potential.

3 .
The elastic behavior of all lengths of NWs coincides with fracture stress, and the fracture stress of the NW increases as its length decreases.The elastic modulus of NW decreases as its length increases, which eventually reaches an asymptotic value at 150 Å. 4. Consideration of surface charges improves the mechanical properties of NWs.The Wolf cut-off function

Figure 11 .
Figure 11.Buckingham potential with Binks constants for hexagon ZnO NW of diameter 70 Å and length of 55 Å: (a) Wolf uniform charge, (b) surface charges Ewald, (c) Wolf surface charge fractional iconicity, and (d) surface charges Wolf.

Figure 12 .
Figure 12.Effect of charge (uniform versus surface charges) and cutoff function type (Wolf versus Ewald) on the stress-strain curves of a ZnO nanowire of hexagonal shape and diameter of 70 Å.

Figure 15 .
Figure15.Effect of strain rate on the stress-strain curves of a ZnO nanowire of hexagonal shape and diameter of 70 Å.

Figure 13 .
Figure 13.Atomic charges (e) at the left and right sides of the NW of diameter 16 Å upon fracture: (a) oxygen atoms and (b) zinc atoms.Atomic charges at the middle section of the NW of diameter 16 Å upon fracture: (c) oxygen atoms and (d) zinc atoms.

Figure 14 .
Figure 14.Atomic charges at the left and right sides of the NW of diameter 94 Å upon fracture: (a) oxygen atoms and (b) zinc atoms.Atomic charges at the middle section of the NW of diameter 94 Å upon fracture: (c) oxygen atoms and (d) zinc atoms (positive charges).

Table 1 .
Geometrical, loading, and MD parameters considered in the present study.

Table 5 .
Comparison of the elastic modulus of the ZnO NWs predicted by the buckingham interatomic potential with values obtained in the literature at 300 K.

Table 6 .
Elastic properties of ZnO NWs as determined by three interatomic potentials.

Table 7 .
Elastic properties and percentage surface atoms for the three shapes considered.Nanowire shape D circle (Å)

Table 8 .
The results of the effect of potential parameters on the mechanical properties of ZnO NWs.