Piezoelectric effects in boron nitride nanotubes predicted by the atomistic finite element method and molecular mechanics

We calculate the tensile and shear moduli of a series of boron nitride nanotubes and their piezoelectric response to applied loads. We compare in detail results from a simple molecular mechanics (MM) potential, the universal force field, with those from the atomistic finite element method (AFEM) using both Euler–Bernoulli and Timoshenko beam formulations. The MM energy minimisations are much more successful than those using the AFEM, and we analyse the failure of the latter approach both qualitatively and quantitatively.


Introduction
The interest in nanotubes has continued to grow ever since they were first shown to be a possible molecular structure for elemental carbon [1]. This discovery led to studies of other potential candidate materials with the potential to form nanotubes. One of the most promising was boron nitride as it was known to form both cubic and hexagonal polymorphs, as an almost exact III-V analogue of carbon. Boron nitride nanotubes (BNNTs) were predicted in 1994 [2] and shown to be semiconductors with a large band gap. The first successful synthesis of BNNTs was subsequently reported in 1995 by Chopra et al [3] using the plasma discharge method. In 2003 Mele and Král [4] used a Berry phase approach to predict piezoelectric properties of BNNTs caused by symmetry breaking of the structure under deformation. They showed that BNNTs would generate a coupled electric dipole dependent on their chirality and the loading type to which they were subjected. Armchair nanotubes subjected to torsion would generate a coupled electric dipole, while zigzag nanotubes would respond to elongation (as shown in figure 1). Further, their work showed how the strength of the coupling varied with the chiral angle for a given type of loading. The predicted values for the piezoelectric coupling tensor entries suggest that BNNTs are better piezoelectric materials than polymer piezoelectric materials [5]. Their piezoelectric properties coupled with their axial stiffness offers potential applications in composites [6] and nanoelectro-mechanical devices (NEMs) [7]. A further benefit of boron nitride materials is their large neutron absorption cross sectional area, adding to their value in space vehicle applications [6].
One technique for modelling nanostructures made of large numbers of atoms is the hybrid molecular mechanics (MM) finite element method which we shall refer to as the atomistic finite element method or AFEM. This method was first used to model carbon nanotubes by Li and Chou [8] in 2003 and has since been used to model graphene, hexagonal boron nitride, BNNTs and zirconia nanotubes by various authors [9][10][11][12]. The method has even been applied to deoxyribonucleic acid as a way to study the mechanics of DNA [13]. The method is based on the concept of taking energy equivalence between the strain energy of deformed beams and the deformation energy of molecular bonds as given by MM potential models. The chemical bonds between atoms are modelled as beams which can stretch, twist, bend and shear allowing the application of standard structural finite element techniques to determine the response of nanostructures to an applied force or deformation. Li and Chou demonstrated this method using both static [14] and dynamic simulations [15].
Much work using this method and related formulations has adopted different element types to model the molecular bonds. Meo and Rossi used non-linear spring elements and torsional spring elements to model bond extension and angular deformation properties of carbon-carbon bonds in graphene and carbon nanotubes [16]. Adhikari et al used the same method as Li and Chou to model DNA with the addition of elastic spar elements for the hydrogen bonds [13]. Boldrin et al [17] used the same principles as the original authors but used deep shear Timoshenko type beam elements to calculate the mechanical properties of boron nitride nanosheets. These beam elements differ from the slender Euler-Bernoulli elements used elsewhere by including shear deformation, a deformation mode that is relevant to beams with a low length-width ratio. It can be shown that this ratio is low for the effective beams employed in the AFEM. Boron nitride structures subject to an external electric field have been modelled by Zhang et al using the Li and Chou method [11]. However, their choice of atomic charges is based on a simplified ionic model that does not necessarily realistically represent the charge density distributions. Their assumption of charges of +3e and −3e, for the boron and nitrogen atoms, respectively, is intuitive but it is not consistent with results from population analysis of ab initio quantum mechanical calculations [18,19]. Nasdala et al developed a multi element approach (MDFEM) that uses different element types to represent the force terms in the underlying MM model [20]. To account for different types of forces, the elements consist of two node spring elements to carry bond elongation, three node coupled bond elements to carry angular deformation and four node elements to carry torsion. Although this method is more complicated than the method proposed by Li and Chou, a single bond will require many elements in order to capture the different deformation mechanisms, it only requires three translational degrees of freedom at each node instead of six. Another somewhat more simplified multi-element approach is that of Giannopoulos et al which uses one type of element to represent bond stretching and another element, which connects two atoms that share a bond with a common third atom, to carry the angular and torsional load [21].
The AFEM has been employed almost solely for the purpose of examining the stiffnesses of molecular structures. Moreover, the initial geometries used have been idealised nanotube geometries, where all atoms lie on the surface of a cylinder, as opposed to energy minimised structures. However, we seek to add additional functionality to calculate piezoelectric properties, such as those of BNNTs. This requires effective charges to be calculated for the atoms based on the molecular geometry and structure and a family of methods that provide a rapid method for effective charge calculations has been developed from the premise of electronegativity equalization [22][23][24][25]. The work presented here seeks to determine the piezoelectric tensor coefficients of BNNTs of a variety of chiralities using the AFEM and an MM energy minimisation approach. The dependence of the piezoelectric properties of BNNTs on the atomic displacements makes this a good test of the AFEM. The work compares the results of the AFEM with those of MM simulations of the same structures to highlight the differences in the outcomes from the different simulation methods.

Atomistic finite element method
The AFEM seeks to model molecular bonds using beam theory and the material and geometric properties of the beams have to be defined. The properties of the beams are deduced from empirical force models developed for MM simulations. Common empirical force models express the total energy of a molecule, U Total , as å å å Bonds Angles Torsions

Pairs vdW Coulomb
U ij is the bond stretching term which represents the energy associated with a change in the bond length between atoms i and j. Similarly, U ijk is the bond angle energy which depends on the angle between bonds ij and jk where atom j is common to both bonds. U ijkl is the torsion energy, the change in energy as a function of the angle between bonds ij and kl with atoms j and k connected by a third bond. U vdW and U Coulomb are van der Waals and electrostatic energy, respectively, and are considered non-bonded energy terms as they are calculated for pairs of atoms which need not be covalently bonded. These non-bonded terms have no counterpart in beam deformation and are neglected in the AFEM. These force models can use linear relationships between force and displacement or they can use non-linear relationships such as Morse potentials and truncated Fourier series. It is possible to approximate, using Taylor series, these non-linear potentials as linear at least close to equilibrium bond lengths, angles and dihedral angles thus allowing the determination of properties of the modelled structure due to small strains. The linearised force coefficients can then be equated to the structural properties of beams using and using either Here, q D D r , ij ijk and f D ijkl are the bond elongation, angular and dihedral deformations, respectively, and are determined from the difference between their deformed values and their relaxed values which were calculated by minimising the energy of the molecule. The force constants k ij , k ijk and k ijkl define the relationship between the different deformations and the change in energy. E A I L G , , , , and J are the material and geometric properties of beams, namely the Young's modulus, cross sectional area, second moment of area, length, shear modulus and angular moment of area, respectively. The length of the beam is the relaxed length of the molecular bond and the beam cross section is taken as a solid circle to provide symmetrical bending properties. Finally, Φ in equation (5) is the shear deformation constant and is given by is the shear area and κ is the shear correction factor, for which there are several formulations available, such as that by Cowper [26], where ν is the Poisson's ratio of the hypothetical material of the beam representing the atomic interactions.
The choice of beam model determines whether equations (4) or (5) is used (equation (4) for slender beams and (5) for deep shear beams) and in turn affects the need for determining the beam radius and Poisson's ratio. In the slender beam model, also referred to as an Euler-Bernoulli beam model, there is no need to determine a Poisson's ratio or radius as the MM derived constants can be used directly to construct the entries in the element stiffness matrix [8]. The radii of these beams can be calculated by combining equations (2) and (4) substituting appropriate functions of radius for the cross sectional area and second moment of area for a solid cylindrical tube, and then solving for the radius to give The radius of the beam produced using this beam model results in a significant radius-length ratio which, in standard beam theory, requires the use of a deep shear, Timoshenko type beam model. Euler-Bernoulli slender beams have been used in previous work with apparent success so whether one needs to be concerned with deep shear for modelling molecular bonds rather than using slender beams remains an open question.
The deep shear beam model requires calculation of the shear deformation constant and the shear correction factor which requires finding appropriate values for radius and Poisson's ratio. This is achieved by substituting the appropriate functions of radius for the area, second moment of area and polar moment of area into equations (2), (3), (5) and (6), and substituting a suitable equation for the shear correction factor into (6) before combining equations (2), (3) and (5). The result (after some protracted manipulation) provides a function relating radius and Poisson's ratio with the force constants and the beam/bond lengths appearing as coefficients. The resulting formula, assuming the shear correction factor as given in equation (7) is used, is n n This can be solved for coupled values of the radius and Poisson's ratio by minimising the difference between the value k ijk , as generated from the MM force field, and the right hand side of equation (9). A suitable pair of values can be chosen by applying an additional constraint of material isotropy, The isotropy condition is imposed for simplicity of the model rather than for any physical considerations. We can rewrite this by solving equations (2) and (3) for E and G, respectively, and then substituting into equation (10) The MM constants thus provide all the required information to model the molecular bonds using beams. These can then be constructed to form a frame structure using standard finite element techniques.

Calculating piezoelectric properties
In order to use AFEM and MM to calculate the piezoelectric properties of BNNTs the dipole density of the tubes and its rate of change due to deformation must be calculated. The dipole density can be found by summing the atomic dipole moments and dividing by the surface area of the tube. The areal dipole density is used, rather than the volumetric density, which would normally be used in calculations for bulk material properties, due to the difficulty in assigning a wall thickness to the hollow nanotube [5]. The use of per area rather than per volume densities also changes the nature of the piezoelectric tensor and the units of its entries from Cm -2 to Cm -1 . Nanotubes are considered to be a one-dimensional structure due to their high length-radius ratio. This low dimensionality allows the neglect of all but uniaxial tension and torsional deformation in determining their piezoelectric properties.
Atomic dipole moments can be calculated using atomic polarisabilities and the electric field at the atom sites using where p A is the atomic dipole moment, α is the atomic polarisability and ( ) E r A is the electric field at the atom position, r A . Atomic polarisabilities are available from the literature and are either calculated using quantum mechanics or are determined experimentally [27]. The polarisabilities of boron and nitrogen used in this work were´-3.03 10 24 cm 3 and´-1.10 10 24 cm 3 , respectively [28].
To calculate the electric field the position of the atom centres and the values of the effective charges on the atoms must be known. These effective charges, or partial atomic charges, represent the charge of the atoms if their electric field was due to a point charge at the atomic nuclei. The atom positions are calculated using MM or AFEM; they are the same as the location of the nodes in the finite element simulation. There are many methods available for calculating effective atomic charges but the most convenient are simple approaches based on electronegativities and atomic hardness or idempotential [23-25, 29, 30]. The basis of such methods is electronegativity equalisation, which involves the expansion of chemical potential, E, of an atom, A, in a molecule so where Q is the charge on the atom [23]. If the energy E A at Q = 0 is taken as the zero point of the energy scale, the atomic energy at = + Q 1 is equivalent to the ionisation energy of atom A; the change in the energy is due to the outermost electron being removed from the atom. In the same manner, the energy at = -Q 1 is equivalent to its electron affinity, the change in energy on the addition of an electron to the atom, and it is possible to show that where E i is the first ionisation potential and E ea is the electron affinity of the atom. Equation (14) is the definition of the Mulliken electronegativity [31], c A . This property describes approximately the ability of an atom to attract electrons. A further expression which can be derived from equation (13) is The idempotential, J AA , describes the Coulomb repulsion between two electrons in the f A orbital [24]. This allows the expression of the chemical potential of the molecule as the sum of the atomic potentials plus the Coulomb potential generated due to atoms with non-zero charge thus, å å å Here, Q Q J A B AB is the Coulomb interaction energy of atom A and atom B. The first sum on the left hand side of equation (16) provides the energy contribution of individual atoms and the second sum provides the contribution from the interaction between different atoms. Taking the derivative of equation (16) with respect to Q A it is possible, after some manipulation, to recast the problem as a linear set of equations that can be solved to yield partial atomic charges. This method assumes that unequal electronegativities will lead to the more electronegative atoms attracting electrons. This lowers the partial charge of atoms with greater electronegativity at the expense of the other atoms. The electric field created by the partial charges should counteract the electronegativity of the atoms and lead to a stable system of partial atomic charges. Once the atomic positions and charges are known the electric field at the individual atom sites is calculated using Coulomb's law.
The piezoelectric property of interest in this work is the piezoelectric coupling tensor for the BNNTs, Here P is the areal dipole density vector and S is the strain vector in a constant electric field. Considering the nanotubes as one-dimensional allows the neglect of all but the entry in P that is parallel to the tube axis. Likewise, the only entries for the strain vector are axial elongation and axial torsion. The non-zero entries of the piezoelectric coupling tensor all take the same value for BNNTs so it is only necessary to calculate a single derivative. In this work the effect of torsional and extensional strain on the dipole moment parallel to the tube axis was studied. The dipole density per unit area for the tubes is given by Tube where P is the dipole moment density, R Tube is the radius of the nanotube and L is the length of the nanotube. The shear strain, γ, due to the torsion in the nanotube was calculated using where θ is the torsion angle. The piezoelectric tensor entry was generated by simulating a number of tubes undergoing varying amounts of axial torsion and using finite differences to approximate the derivative.

Boron nitride nanotube simulations
To test the AFEM method a range of BNNTs were modelled using a bespoke Matlab script. The tubes varied in radius from 1.6-11.0 Å covering a number of chiralities in this range. The chirality of the nanotubes is determined by the vector that connects two equivalent lattice points of a hexagonal sheet. This vector describes the direction of rolling if the tube were to be produced by rolling a sheet to form the tube. The chirality can be described either by two integer values, that act as multipliers of the two basis vectors needed to span the two-dimensional hexagonal lattice, or by an angle formed between the vector and a predefined direction. In this work a 0°chiral angle is the same as an (n, n) chiral vector, a vector which points along an armchair edge of the hexagonal lattice. See figure 2 for illustration of this point. The length of tubes was varied as a function of tube radius in order to ensure a length-radius ratio of 20. The tubes were terminated with hydrogen atoms in order to provide a physically realistic structure, necessary for the charge equilibration methods. A range of electronegativity equalisation methods are available, such as the charge equilibration method (QEq) [24] and the partial atomic charges and hardness analysis (PACHA) formalism [25], which utilise the concepts presented in section 3. Atomic charges in this work were generated using the PACHA formalism as implemented in the General Lattice Utility Program (GULP) [32]. Initially the QEq method was used but its recursive procedure failed to converge for many of the tubes so PACHA was adopted instead. There was a small variation (<0.02 e) in the calculated average atomic charges depending on chirality and radius. The universal force field (UFF) [33] potential model was chosen to provide the force constants used in the AFEM model and the MM optimisations, thus allowing for direct comparison of results. UFF was chosen as it is relatively simple, has the option of linear or non-linear force expressions and provides the required parameters for boron, nitrogen and hydrogen atoms.
The first step in the simulations was optimisation of the nanotube geometry using MM energy minimisation in GULP. AFEM requires that the original atomic positions are those of the minimum energy for the molecule [34]. The original input geometry, generated by our script, assumed that all the atoms lie on the surface of a cylinder which is not an accurate model of a boron nitride nanotube [35]. Structure optimisation resulted in the atoms forming a buckled tube surface as can be seen in figure 3.
The depth of the buckling of the surface varies with the tube radii, as the radius increases the surface buckling lessens. The relationship between the depth of the buckled surface and the nanotube radius is shown in figure 4, including data generated using the density functional theory by Wirtz et al [36].
The force constants for the AFEM model were generated by making linear approximations of the force equations generated from the UFF potential about the equilibrium bond lengths, angles and torsion angles of the optimised geometry for each nanotube. Table 1 shows some examples of the force constants and beam properties that were used in the AFEM simulations. This harmonic potential approximation required all deformations to be kept small.
With the molecular geometry optimised, two different loading scenarios were simulated: torsion and extension. For the torsional simulations the hydrogen atoms and the first connected boron and nitrogen atoms at one end of the tube were rotated about the tube axis by progressively larger angles.
The angles of rotation were chosen by setting a limit of 0.2% on the shear strain of a continuous hollow tube with the same length and radius of the nanotube under consideration. The translational degrees of freedom perpendicular to the tube axis of the atoms at both ends of the tubes were then removed in order to provide the boundary conditions for both the MM and atomistic finite element simulations. This allowed the tubes to change in length while the relative rotation of the two ends was kept constant. For the extension simulations a similar approach was taken but with one end of the tubes subjected to an increasing displacement along the tube axis. The extension of the tubes was limited to 0.1% strain and only the axial coordinate was held constant for the boundary atoms at both ends of the tubes.
Finally, the tubes were optimised in GULP (using UFF) with the applied restrictions on the end atoms and using second derivative methods employing the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm [37]. For the largest tubes limited memory BFGS (lBFGS) was used for the first optimisation steps and the final steps employed the full BFGS algorithm. Simultaneously, the tube geometries were input into the AFEM code and a sparse stiffness matrix generated by assembling translated element stiffness matrices. The boundary conditions were applied, the stiffness matrix subdivided and the system of equations solved using Matlab's inbuilt left divide method. After the final atom positions were calculated in each method the dipole density was calculated and a finite difference scheme used to determine the rate of change of dipole density due to increasing strain in the tubes.

Results and discussion
The results of the torsional simulations are shown in figures 5 and 6. Figure 5 shows the relationship between the axial component of the polarisation vector and the shear strain in the nanotube. Figure 6 plots the piezoelectric response of a number of nanotubes of different chirality subjected to a torsional loading. The chiral angle varies between 0°, an armchair nanotube with chiral vector (n, n), and 30°, a zig-zag nanotube with chiral vector ( ) n, 0 . Although the overall trend in the results is similar to that in other work, neither method accurately follows those from more accurate quantum mechanical electronic structure methods [5]. The difference between the AFEM and the MM results are worth noting, given that the beam properties used in the AFEM models were based on the UFF force field that was used in the MM simulations. The most significant failure of AFEM was the inability to produce a zero value for the piezoelectric tensor coefficient for zig-zag type BNNTs under axial torsional loading. The results generated using MM, even using a relatively old potential model such as UFF, were surprisingly close to those generated using more accurate quantum mechanical electronic structure methods [5].
The results for the extension based simulations are shown in figures 7 and 8. The non-linear results in figure 7(b) indicate the AFEM method fails to predict the correct behaviour of the polarisation vector. The MM based results can be seen clearly in figure 8(b) and show the correct trend with Table 1. The linearised force constants and effective beam properties used for both the Euler-Bernoulli and Timoshenko beam models that were used for the (16,16) nanotube. The beam radius for the Euler-Bernoulli beam model is included for comparison with the Timoshenko beam model radius only. The beam radii and Poisson's ratio were found using the methods described in section 2.  Results from Sai and Mele [5] are included for comparison. Qualitatively, both the AFEM plots show similar responses compared to the MM as the chiral angle is varied, but quantitatively they show a smaller piezoelectric response and erroneous non-zero tensor coefficients for the 30°chiral angle tubes (i.e. zigzag tubes). respect to the chiral angle but the method appears to overestimate the value of the piezoelectric tensor coefficient. The work of Sai and Mele [5] predicts a lower (more negative) value for the piezoelectric coefficient at 30°for the extension case than at 0°for the torsion case and the MM results shown here also display this behaviour. The AFEM results fail to reproduce either the trend or the values expected based on either the MM results or results produced by other authors using higher fidelity models. Figure 8 To examine the differences between the AFEM and MM results a comparison of the final positions of the atoms, generated using AFEM and MM, was performed for each nanotube. The atomic positions were transformed from Cartesian to cylindrical polar coordinates and then the atomic displacements were calculated; for both AFEM and MM the change in position from the optimised geometry to the twisted are extended geometries. An example of the results from this analysis can be seen in figures 9 and 10. For the torsional case seen in figure 10, the AFEM results show a linear increase in angular displacement, which is expected of a continuum mechanics based model. The radial and axial displacements bare little similarity to those generated using the MM model. The MM method showed the tubes bulging in different positions along the tube from the AFEM and the bulges were coupled to angular displacement behaviour. The two methods gave very different results for the tension load case. The radial and angular displacements for the AFEM simulations were entirely different to those generated using the MM method.  To understand further the discrepancy between the two methods the energies relating to different molecular deformation modes were compared. This analysis was carried out on a (16,16) armchair BNNT and the results are in tables 2 and 3. These results clearly show that the Timoshenko formulation, although arguably correct for the aspect ratio of the beams used to model the molecular bonds, results in significantly lower total energies when compared with the Euler-Bernoulli formulation or MM. This result is expected as the more compliant nature of the Timoshenko beam formulation will result in lower input energies required to produce the same amount of deformation. The Euler-Bernoulli formulation performed reasonably well in the torsion load case but poorly in the tensile load case when compared with the MM minimisation results. This difference appeared to depend on the non-bonded energy terms ignored in the AFEM formulation.

Energy component Euler-Bernoulli Timoshenko
The energy calculations can be used to calculate values for the Young's modulus and the shear modulus of the   The AFEM method implements the bond stretch term in exactly the same manner as the MM potential. However, the bonded terms from the potential, angle bending, torsion, and out-of-plane bending, are present in the AFEM model only after considerable simplifications. These simplifications reduce the number of atoms involved in the bonded interaction terms by removing the three-body bond angle and fourbody bond torsion interactions. The beam models attempt to replace these terms with two-body bending and torsion terms which require the addition of rotational degrees of freedom. These degrees of freedom are not present in MM methods and the stiffness matrix coefficients for these terms rely on various assumptions which are questionable in a molecular context. This leads to situations where pure bond torsions that occur in MM can be modelled as beam bending in AFEM, as shown by Nasdala et al [41]. Complete neglect of the non-bonded interaction terms will always produce significant differences when comparing MM minimisations, using potentials with bonded and non-bonded terms, with AFEM regardless of how well the bonded terms are reproduced.
Overall, it is clear that the AFEM generates quantitatively different atomic displacements compared with MM methods. This results in the quantitatively and qualitatively incorrect behaviour of the piezoelectric tensor values, energies due to deformation and mechanical properties of BNNTs.

Conclusions
This paper has shown the limitations of using a simple implementation of an AFEM for investigating the piezoelectric behaviour of BNNTs. Further, it shows the failure of the method to accurately generate the displacement due to external loading for molecules with significant Coulombic interactions. We have applied methods for investigating the qualitative and quantitative behaviour of such models which allow direct comparison with higher fidelity models. The work raises the question of whether simple beam models, that use force-displacement constants that couple movement of only the bonded neighbour atoms, are suitable for molecular studies. The work also shows that neglecting non-bonded interaction terms from an empirical MM potential will generate spurious results.
The use of MM, coupled with dipole density calculations based on partial atomic charges calculated using electronegativity equalisation methods, shows some promise for calculating the piezoelectric properties of molecules.