Molecular Insight into Structural and Mechanical Properties of Halloysite Structure

In this study, we simulated the rolling mechanism of Halloysite by molecular dynamics (MD) under different conditions. We have illustrated that the transformation from slab Halloysite to scroll shape depends on the initial geometry, dimension and proper selection of the interatomic potential. Also, the molecular simulation was conducted to determine the mechanical properties of Halloysite under different conditions. The results show that the Elastic modulus of the armchair Nano scroll was higher than the zigzag with similar dimensions and that Young's modulus of both arrangements decreases with increased radius. Moreover, with an increasing radius (>20 Å), Young's modulus of a Halloysite nano-scroll approaches that of the Halloysite slab configuration. Finally, the tensile strain of a Halloysite nanosheet was 0.08±0.04. The result of this study is a great help for understanding Halloysite, which can be used for designing nanocomposites.

Halloysite nanotube (HNT) has become an interesting natural nanoparticle for researchers due to inexpensively and availability in several regions. Also, it is a reliable and compatible particle that can enhance polymer functionality. HNT is basically a mineral alumina-silicate (Al2Si2O5(OH)4.2H2O) composed of a two-layer structure, where the exterior structure is composed of siloxane, and the inner structure consists of alumina.
Halloysite forms tubes with a length of 600-1000 nm and an internal radius of 12-30 nm, an external radius of 30-100 nm as shown in 1,2 . The inner cylindrical hollow of Halloysite can store and release other molecules, for example, drugs, proteins, and anti-corrosion and antioxidants material. Thus, Halloysite is a good candidate for purposes such as radius is in the range of 230-300 GPA, and 300-340 GPa for the armchair arrangement with 21-46 Å radius 1 . Other researchers also determined that the HNTs are significantly flexible and can be bent 90 degrees without breaking 8 while a few investigations have been conducted to analyze the influence of geometry on stiffness 9 , strength 1,7,10,11 . Molecular simulation is a promising method to investigate properties of the chemical structure at the atomic level, and it has been successfully used for several molecular structures 12,13 and interlayer formation 14,15 . We apply molecular simulation to study the influence of geometry and size of the HNT on the mechanical properties.
Molecular Mechanics is basically a method that uses potentials that are functions of intramolecular bonds lengths, bond angles, and bond dihedral angles of the molecules.
These potentials are determined from experiments or DFT calculations. Some common potential energy functions used are CVFF, Dreiding, UFF, and PCFF [16][17][18][19] . Generally, the choice of a particular potential depends on the properties of the structure.
The first four terms in equation 1 are, respectively to stretch bonds energy (b), bend angles energy(θ) torsion angles ( ) and distort planar atoms out of the reference formed plane (χ). The final terms represent the nonbonded interactions as a sum of a van der Waals (vdW) term and electrostatic (Coulomb). All functions depend on the distances (rij) between atom pairs. Each potential provides those coefficients and is designed for various purposes. Despite the compatibility of different potentials to simulate bonds with experimental observation, accurate simulation of non-covalent bonds is more challenging in classical potentials. Therefore, each forcefield introduces a different approach to simulate interatomic interactions which make the different results of molecular simulation, therefore, the polymer consistent forcefield (PCFF) with an interface extension called (PCFF-IFF) has been used for molecular dynamic in this study 20,21 and Dreiding and COMPASS were used in the molecular mechanics 18,19 to compare the results with previous research.
The previous computer simulation of rolling Halloysite 2 was concentrated on applying the different schemes to take into account long-ranged Coulomb interaction with keeping the periodical boundary conditions which leads to calculating 13 Å distance for inter-layer that is 7.5 Å 2,7 for the Halloysite nanotubes. The most notable difference between our simulation and rest is that the internal radius and interlayer distance, which is 8 Å with employing the PCFF-IFF potential and study curvature effects on the physical properties. Also, despite the increasing applicability and studies on the HNT, there is still no appropriate research on the microstructural details of Halloysite under deformation. The question is that the defects would change the mechanical properties of Halloysite.
Therefore, in this study also, the effects of defect size on the mechanical properties of Halloysite were considered by using MD simulations.
This study is prepared as follows. In Section 2, we will investigate the rolling mechanisms of Halloysite nanotubes with Molecular Mechanics, then in Section 3, atomic simulation with DFT, MM and MD methods will be performed to calculate the mechanical properties of Halloysite. The results are provided in Section 4 and compared to the previous results. The conclusions and suggestions of this research are presented in Section 5.

Computational Method of rolling mechanism and discussion
HNT can be experimentally generated by mixing kaolin clay in water for a long time with a detergent. 2   properties with using the simplified HNT structure, for example, the study on the singlewalled Halloysite nanotube 1,3,30,31 , Halloysite nano scroll with overlapping arms 32 and two slab surface arrangement 8 . In those studies, the researchers mainly used a pre-defined structure and assigned structural optimization. The main issue in these cases is the limited size of the model, which restricts the ability to describe a realistic Halloysite structure.
Here, we evaluated the structure of the HNT to study the rolling mechanism. All the MD simulations are carried out by using the "large scale atomic/molecular massively parallel simulator (LAMMPS) 33 " packages. In this section, because of computational accuracy, we used the polymer consistent forcefield (PCFF) with an interface extension (PCFF-IFF) 20,21 .
Also, we used the particle-particle, particle-mesh (PPPM) 34 algorithm with a relative force accuracy of 10 -5 and real space cut-off of 18 Å.
In the previous literature, the formation of HNT is associate with rolling a Halloysite sheet around the c-axis, and ac-axis 2 . In our simulation, we applied a supercell including 24480 atoms. It was generated by replicating 30 and 24 unit cells along with a and c crystal axis, respectively. Then, we removed the periodic boundary conditions, as seen in figure   2. The initial unit cell of Halloysite was adapted from ref. 29 . We applying the conjugate gradients algorithm 35   The single Halloysite nanotube can be generated from the smaller initial layer as described in figure 5. A supercell was created of 18 and 15 replicated unit cells along a and c crystallin axis, respectively. Then, we applying the Molecular mechanic to optimize the structure by conjugate gradients algorithm 33 . We can observe that the horizontal geometry of the HNT layer is unstable and rolling starts at corners, similar to the experimental result demonstrated roll structure in figure 6 2 . One simple explanation for rolling is that the bond length of Al-dioctahedral and Si-tetrahedral layers are different.
The final configuration has good agreement with atomic force microscopy (AFM) imaging and transmission scanning of Halloysite nanotubes in figure  observations. However, the inner radius of the experimental dimensions is larger than our simulation.

Mechanical Properties of the HNT
The understanding of the influence of dimensions on the mechanical properties of nanoparticles is the main subject to classify their applicability and functionality. Also, experimental analyses in the nanoscale are considerably challenging due to the complexities of the preparation process for conventional mechanical tests. During the last three decades, several techniques have been proposed to investigate the behavior of nanomaterials by using atomic force microscopy 36,37 . Various functions are available for the estimation of elastic moduli, which we can categorize as static and dynamic ways.

Static approach
The static approach ignores the fluctuations which happened due to thermal movement.
This technique is usually applied to a geometry optimized structure. Analytical calculation on the response in an optimized structure by using the cell Hessian matrix, or by employing a small and finite strain. Young's moduli in the static approach are given by eq 2 38 ; Where V0 is the initial volume, U is the total energy and ε is the strain. In this method, a small finite strain is applied along crystal axes then the elastic constants are calculated. The six stress components are determined from a summation over all atoms in where i,j = x, y, and z, Mk refers to the atom mass, νki identifies the atomic velocity, V refers to the volume of the structure, rki refers to the atom location, and fkj refers to the atomic force. With applying Hooke's low and Voigt notation, the stiffness tensor Cijkl, can be obtained as the following equation: where σij refers to the stress tensor, εij refers to the strain tensor, and Cijkl (i, j, k, l = x, y, z) refers to the stiffness tensor. Because of symmetry, the stiffness tensor reduced to 21 components which can be obtained from eq 5 37,38 ; where U refers to the potential energy. With applying Voigt notation, the Young moduli, Ei, determined by the stiffness tensor, Siijj, where Sij =Cij -1 , therefore Ei = 1/Sii 13,37 .
In our case, DFT and MM simulations were conducted to calculate the mechanical properties of HNTs with a defined periodic optimized structure in a static approach. For this purpose, we are using DFT analyses of the local structure of a single repeating unit of

Dynamic approach
The dynamic technique applies to an arrangement taken from the molecular dynamic simulation. Dynamic approach examining the elastic moduli through applying a constant strain in a specific direction and calculate the consequent average stress at, for example, constant temperature. In our case, uniaxial tensile with the constant strain rate was applied for all models at a constant temperature and the Elastic modulus was calculated from corresponding stress-strain curves. 1.6×10 6 s −1 respectively at 300k temperature. Also, as illustrated in Figure 7, three uniaxial tension tests were conducted along crystal axes. The models were deformed at 300k temperature with a strain rate equal to 10 9 s −1 with a canonical NVE ensemble with 0.1 fs time step and total time 1.6×10 5 s -1 . We used the PPPM method with 10 -5 force accuracy 34 .
The results are presented in figure 10.

Figure 7: Schematic of loading direction (a) OA axis (b) OB axis (c) OC axis (d) The simulation cell of HNT is arranged into four layers, silicon (yellow) aluminum (pink). Oxygen and hydrogen atoms are displayed
Moreover, molecular dynamics simulation was conduct to reveal the effects of different defective Halloysite under tension on their mechanical properties. Four types of defects were constructed on the Halloysite sheet of mode 1. The defect sizes follow the order Si3 < Al3 < Si6 < Al6, with the removal of 3 and 6 of Si atoms and Al atoms, respectively. Also, the defect was located in the middle of each layer as shown in figure 8.

Model 2:
In addition, the MD method is used to study the influence of zigzag (ZZ) and armchair(ac) arrangement on the mechanical properties of Halloysite nano scroll. To this purpose, the pre-rolled models were generated from 18×15×1optimized unit cells of HNT including 9180 atoms and rolled with a 12 Å radius with 75 Å ×75 Å ×150 Å crystal dimension, as displayed in figure 8. 44 Then geometry optimization by using conjugate gradient 45 methods was conducted and it was continued until the subsequent difference energy steps are less than 10 −4 kcal/mol. Then the structure was equilibrated with an NpT ensemble at 1 atm pressure and 300 K for 50 ps. Then all models HNTs were deformed along the tube axis at 300k temperature with a strain rate equal to10 9 s −1 . All simulations were incorporated with a canonical NVE ensemble with 0.1 fs time step. We used the PPPM method with 10 -5 force accuracy 34 . In this section, to avoid the effect of the cutoff radius, Lennard Joan's cutoff radius is chosen as 18 Å 46 . LAMMPS is used for all the MD simulations 33 . PCFF-IFF has been used for both arrangements 20,21 . The results are presented in figure 15.

Result and discussion of atomic simulation
Elastic stiffness tensors derived from the static approach are presented in tables 1 and 2.  to simulation non-bonded interaction in different potential and charges atoms which is constants in the MM but it could be flexible in the relaxation process in DFT calculations, k-point and cut-off radius in the computing process.
To determine young's modulus of MD simulation of model 1 and 2, an available approach is applying a fitted straight line to the stress-strain curve if the deformation is small enough. In our case, the strain that is less than 0.04 can be described as the linear  In this classical force field, bonds can be stretched, twisted and interacts with the other atoms or molecules in the neighborhood, but cannot be broken. Thus, the failure states are more challenging through applying PCFF-IFF. Therefore, the methodology of Hua Yang and et al 7 was selected to determine the failure process of the system. According to their definition, the failure state was defined as the cutoff distances for bond breakage.
Similarly, we employed a radial distribution function (RDF) to specify the criteria to define the break of individual bonds. When the first smallest amount of the radial distribution function curves arrives mean the possibility of detecting atom pairs with the detachment length is approximately zero. Hence, we used this assumption that the bond is broken if the length between atom pairs is larger than this amount. This allows Halloysite bonds to be broken in model 1 of MD simulation without using quantum mechanics. Fig. 11 displays the radial distribution function and table 3 is listed the cutoff distances for bond breakage.

Surface energy
The After the constructed model, DFT was carried out with VASP 48 with the PBE gradientcorrected functionals 49 to calculate the total energy of the mentioned system for each increment. The projector augmented wave (PAW) pseudopotential was applied for describing the core electrons interactions. The calculation was performed at zero temperature. Reciprocal space was represented Gamma k-points scheme with 5x5x5 grid meshes. There are 68 atoms in the supercells and the cut-off energy was 600 eV. therefore, and the adhesion energy (γ(r)) is calculated as; Where γ is cohesive energy, Etot is the total energy of the interacting to unit HNT at the relative position, EHNT is the energy of the isolated units HNT and A is the surface area.
Cohesive energy obtained from the DFT calculation was 0.0477 eV/Å 2 , which is in agreement with the MD result (0.029 eV/Å 2 ) reported by Hua Yang and et al 7

and Junhua
Zhao and et al 50 .
Moreover, of the result of model 1 in Fig. 9 shown us, Elastic modulus of armchair arrangements of model 2 was completely higher than that of zigzag arrangements. The zigzag arrangement showed a relatively Young's modulus in is 219.90 GPa. The elastic modulus of the armchair structure is 328.46 GPa. To study the influence of temperature on mechanical properties, one ac and one zz arrangement of HNT with a 12 Å radius was chosen (model 2). Their Elastic modulus was calculated with an increase in temperature from 300 K to 700 K, under tensile loading with a similar strain rate equal to 10 9 S −1 , the result is provided in Fig. 20.

Figure 20: Young's modulus of HNT ac and zz under uniaxial tensile in different temperature
We can observe a decline in Elastic modulus for both models. There is an obvious decrease in Elastic modulus for ac respect to zz. By increasing temperature to 700 k, Elastic modulus decrement is around about 5% and 21% for zz and ac respectively.
At this step, we calculated the Elastic modulus of the ac and zz structure of model 2 while the radius slightly increased from 12 Å to 25 Å and all models have 132 A length.
During the calculation, boundary condition was removed then energy minimization was conducted. Then models were deformed along the tube axis (c crystal axis) at 300k temperature under three different strain rates (10 10 , 10 9 and 10 8 S -1 ). LAMMPS is used for all the MD simulations 33 . The results are provided in Fig. 21. As observed in Fig. 21 and 300-340 GPa for Ar structure 1 . As observed in Fig. 21, Young's modulus for the larger strain rates loading is higher than lower strain rates loading and this is more noticeable at the smaller radius as figure 21, we can observe an obvious reduction of Young's modulus for the larger radius (>20 Å), that explained with increasing radius, the elastic modulus of Halloysite Nano scroll is approaching the Halloysite flat sheet.

Conclusion
In this research, we studied molecular mechanic simulations to illustrate the transformation of a slab surface Halloysite to a single nano scroll structure. It was observed that the initial size of flat surface alumosilicate plays a crucial role to generate a Halloysite nano scroll.
In this Study, we also chose the molecular dynamics simulation to determine the Elastic modulus of Halloysite nanotube with zigzag and armchair arrangements. The elastic modulus of both arrangements is decreased with an increase in radius. Although the Elastic modulus of armchair arrangement was higher than the zigzag arrangement, the zigzag has more stable with increasing radius. The elastic modulus of Halloysite Nano scroll is approaching the Halloysite slab configuration which was calculated by DFT and MM methods. Moreover, the deformation along three directions was determined, and the weakest direction was along the OB crystal axis. Besides, temperature climbing made a significant decrease in elastic modulus in both arrangements. The results of this research propose a great help for understanding Halloysite which can be used for designing nanocomposites.

Acknowledgment
We gratefully acknowledge the funding support received for this project from the NSERC Discovery Grant program for supporting, Computing resources were provided by Compute Canada.

Supplementary Material
Movie M1 and M2 show axial stretch of the armchair and zigzag. M3 illustrates fracture simulation.

Consent for publication N/A
Funding: The Natural Sciences and Engineering Research Council of Canada (NSERC) Grant program for supporting this research.

Conflicts of interest/Competing interests the authors declare that they have no conflicts of interest.
Availability of data and material Data available on request from the authors.
Code availability LAMMPS main inputs available on request from the authors.