Molecular dynamics simulation of pull-out Halloysite nanotube from polyurethane matrix

Molecular dynamics (MD) simulation has been applied to study of pull-out of Halloysite nanotubes (HNTs) from a polyurethane (PU) matrix. First, the Machine learning (ML) particle swarm optimization (PSO) method was used to obtain force field parameters for MD from data of density functional theory (DFT) calculations. The current study shows the possibility of using a PSO technique to modify the force field with DFT data with less than 5 kcal/mol discrepancy. Second, we considered the influence of atomic interface on pulling out of HNT from PU. Energy variation has been proposed as the cohesion strength between matrix and nanoparticle. In addition, the best Lennard Jones parameters in the MD simulation make good agreement with an experimental sample stress-strain response.


Introduction
Nanomaterials and Nano-composites have been studied in various applications. One innovative approach is to employ a polymer in combination with a strong nanomaterial. In this study, we will consider a composite where the polymer matrix is polyurethane (PU) and the nanoscale material is Halloysite nanotube (HNT).
Carbon nanotubes are the most commonly studied nanomaterial; however, toxicity and high cost make them less attractive in the industry, thus it has limited practical applications. Halloysite nanotubes are extensively available around the world at low cost, therefore, it has received more attention from researchers.
Halloysite is a clay tube with 600-1000 nm length and an inner radius of 30 nm, an outer radius of 100 nm. It is a mineral nanoparticle used in several nanomaterials, especially composites. Halloysite is a natural alumina-silicate (Al 2 Si 2 O 5 (OH) 4 Á2H 2 O) with a predominantly hollow structure. HNT consists of a two-layer structure, where the exterior face is a tetrahedral, siloxane structure, and the inner face has an octahedral composition with alumina assemblies as described in Figure 1.
Most research on Halloysite has concentrated on measured mechanical properties under different loading. 2 Lu et al. 3 experimentally measured the elastic modulus. They found that the elastic modulus of the HNT layer is 130 6 24 GPa. Lvov et al. 1 calculated Young's modulus of a multi-wall alumino-silicate nanotube. They obtained Young's modulus in the range of (230-340 GPa). Guimara˜es et al. calculated the young modulus of different arrangement of HNT nanotubes by applying DFT. They obtained the elastic modulus in the range 230-300 GPa for a zigzag arrangement and 300-340 GPa for the armchair arrangement. 4 Moreover, the elastic modulus of HNT declines with a growth radius. Also, other researchers found that HNTs are significantly flexible and can be bent 90°without breaking. 5 Polyurethane is a copolymer formed by a di-isocyanate (OCN-R-NCO) and a polyol or diol (OH-R#-OH), where R is aromatic or aliphatic and R# is a linear hydrocarbon chain. Also, PU is comprised of soft segments (SS) from the polyol and hard segments (HS) of the crystalline di-isocyanate. 6 PU has a variety of uses in a range of products. 7 Furthermore, polyurethane has good corrosion resistance and toughness, which is effective functionality for the defense industries. 8 Recently, new applications of polyurethane have been developed for impact protection of structures, including anti-explosion and ballistic loading. 9 The chemical arrangement of a polyurethane monomer applied in this research is presented in Figure 2. Later in this work, in the fitting process we are using one monomer (n = 1, m = 1) and for uniaxial stretching simulation the chain composed of two repeating monomers (n = 2), and five repeating polyol part (m = 5), and one di-isocyanate part.
Molecular Mechanics uses potentials that are functions of the intramolecular distance between pair atoms. These potentials are determined from experiments or DFT calculations. The potential energy function is derived from positions of the atoms that form a molecule. Some common forcefields are CVFF, Dreiding, UFF, and PCFF-IFF. [10][11][12][13] The total energy of a molecular structure is giving by; U = X all bonds The first four terms in equation (1) introduce the stretched bond energy (b), bond angles energy (u), torsion angle energy (f), and distortion energy respect to a reference plane (x). The final terms represent the non-bonded interactions as a sum of van der Waals and electrostatic. All terms are function of the distances (r ij ) between atoms. Each force field defines a specific  coefficient for each set of atom types. Because of computational accuracy, the polymer consistent forcefield (PCFF) with the interface extension (PCFF-IFF) has been used for this study. 14,15 Li et al 16 studied the hardness of PU-HNT nanocomposites. They indicate that the scratch resistance has been improved by adding 5 wt% HNT. Gaaz et al 17 reported that the optimal level of HNT is 1 wt% for improving tensile strength.
The purpose of machine learning (ML) is to improve efficiency and accuracy. 18,19 Lately, a variety of complex machine learning (ML) algorithms for small and large systems have been recommended for using the total energy data of DFT calculations. 20 Those ML models have been developed for small atomic configurations, which is not suitable for large atomic configurations. 21 In this paper, we introduce an alternative method to improve force fields in MD with the precise results of DFT calculations.
Lately, fitting to DFT results have become a common way to optimize force fields. 10,[22][23][24] Energy-fitting utilizes the atomic energy obtained from a highaccurate DFT simulation. Ercolessi and Adams, 25 introduced force-fitting for the first time, and it was employed in numerous simulations [26][27][28][29] to improve the efficiency of semi-empirical quantum methods. 30,31 Wu et al. 32 optimized parameters of a potential field for long-range interaction through a ceteris paribus optimization method, and Daven et al. 33 optimized Lennard-Jones parameters by a genetic algorithm.
There are several non-linear optimization techniques, such as iterative direct gradient minimization methods which correspond to a local minimum (e.g. Levenberg-Marquardt 34 ), or global minimum explorations (e.g. simulated annealing 35 ), which are computationally expensive for energy matching function. Therefore, we used a particle swarm optimization approach.
Most studies of combined loading response concentrate on the effective response of the composite. [36][37][38] Hence, they do not explicitly focus on the interfacial performance during failure. Subsequently, micromechanical simulations, which predict the effective behavior of nanocomposites are often based on an assumption that the interface is completely bonded or partially bonded. [39][40][41] There are two mechanisms to explain the strengthening of a composite by adding fibers. First is the energy absorption when the fiber pulls out, another is the crack-bridging by the fiber at the crack zone. Usually, fiber pull-out is an essential way for achieving favorable strength. As such the understanding of the mechanism of the pull-out in the HNT-polymer composites are essential. Therefore, at the end of this research, the pull-out simulation has been conducted using the molecular dynamics, to study the effects of the interaction between the HNT-PU nanocomposites.
Overall, the purpose of this research is to develop an accurate description of the interaction energy between HNT and PU chains. This new force field will be applied to molecular dynamics (MD) simulations of HNT/PU composites, which are too expensive to evaluate with DFT. For this purpose, in Section 2 we will introduce our computational technique in detail, then the results are provided in Section 3, followed by the MD pulling-out simulation in Section 4.

The methodology of the fitting process
In order to employ HNT in different applications, it is necessary to consider materials exposed to different loading conditions. The purpose of this section is to use computational techniques for fitting methods at the atomistic level for studying the dynamic behavior of the polymer matrix reinforced with a nanoparticle. Specifically, here we apply computational techniques to fitting the non-bonded energy between HNT and PU obtained from density function theory (DFT) with a form amenable to Molecular Mechanic calculations.
Here are steps for this computation method: 1. The optimized lattice parameter of simple unit cell HNT is found. 2. A slab surface is created according to the optimum unit cell of HNT and the computational cell was extended in the z-direction to create a vacuum space in which to embed a PU monomer. 3. For DFT calculations, the cut-off energy and minimal k-point mesh are determined through convergence studies to be 580 eV and 10 (auto mesh) respectively. These parameters are used for all trajectory calculations to minimize computational cost. 4. With DFT, we calculate the potential energy as a function of the vertical position of the PU monomer relative to the flat surface of the HNT (see Figure 3). 5. We repeat step 4, this time with PCFF-IFF.
Then, the interaction energy is calculated for different vertical positions of the PU. 6. In our approach, the particle swarm optimization (PSO) was applied to find the best fit between DFT and the non-bonded Lennard-Jones (LJ) parameters used in PCFF-IFF.

Preparation of supercell
We used DFT calculations to compute optimum lattice parameters. The system is relaxed with the cut-off energy of 17,300 kcal/mol. The cut-off energy is determined through convergence studies. DFT was carried out with VASP 42 with the PBE gradient-corrected functionals 43 methods. The calculation was performed at zero temperature. Reciprocal space was represented by the Gamma k-points scheme with 5 3 5 3 5 grid meshes. There are 34 atoms in the unit cells. The structures are relaxed without symmetry constraints with the cut-off energy of 750 eV. The projector augmented wave (PAW) pseudopotential 44 was applied for describing the core electron interactions. The result is provided in Table 1. Then a supercell was generated from 3 3 3 3 1 optimized unit cells and the cell extended cell in the z-direction to form a vacuum space to embed a PU chain ( Figure 3). We move the PU away with equal increments (0.35 Å ) to compute the energy of each configuration. The distance was calculated based on the center of mass of the HNT slab and the polymer chain.

The interaction energy calculations (DFT approach)
Using the generated model, we employed the density function (DFT) approach to calculate the non-bonding energy between HNT and PU. The interaction energy between the PU and the HNT was computed by the following equation; Where E tot-DFT refers to the total energy of the structure (HNT and PU), E HNT indicates the energy of the isolated slab HNT, E PU refers to the energy of the isolated polyurethane, and E inter DFT indicates interaction energy between HNT and PU molecules.
Density Functional was performed with the ''Vienna Ab-initio Simulation Package (VASP)'' 42 with the PBE gradient-corrected functionals 43 to compute the total energy of the above structure for each vertical 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 with kpoints grid meshes by an auto mesh with weighing 10. There are 350 atoms in the supercells and the cut-off energy was 13,380 kcal/mol. The cut-off energy is determined through convergence studies. The ions were relaxed with a conjugate gradient algorithm until the energy components on any ions were smaller than 0.25 kcal/mol. Additionally, DFT computation was performed to calculate the energy of a monomer PU and slab HNT separately with a similar set up. The energy for a single polyurethane chain is E PU = 24188.717 kcal/mol and energy for the flat surface of Halloysite is E HNT = 250,656.40 kcal/mol. The result of interaction energy is presented in Figure 4.

The interaction energy calculations (MM approach)
For MM calculations, the same supercell was used (described in Figure 3). Then a minimizing process with  PU and the HNT was calculated by the last two terms of equation (1): Where E tot-MM refers to the total energy of the structure (HNT and PU). E HNT indicates energy of HNT, E PU refers to energy of the polyurethane, and E MM interact indicates interaction energy between HNT and PU molecules.
In the classical molecular dynamic method, the nonbonded energy between the HNT and the PU is calculated with the 9-6 Lennard-Jones potential and Coulomb part as described in equation (4); Where E MM Interact is the non-bonding potential energy, r ij is the gap between a pair of atoms, e`identify the well depth of potential, s refers to equilibrium length, and q is the atomic charges.
In order to compute E Interact MM and eliminate local interaction, we build an isolated optimized polyurethane chain (C 17 N 3 O 3 H 18 ) cell. Then, conjugate gradient minimization was performed at zero temperature to calculate the energy of PU. This process was repeated for the isolated slab HNT as a supercell with a 3 3 3 3 1 repeated unit cell. After the equilibration process, energy for a single polyurethane chain is E PU = 4699.8 kcal/mol and energy for the flat surface of Halloysite is E HNT = 30,471.5 kcal/mol.

The fitting process (PSO approach)
The goal of this section is to minimize the discrepancy between equations (2) and (3) as defined in the following equation (5): The index i refers to unique set of atomic coordinates in the training set.
In order to solve the least-squares problem described in equation (5), we exploring optimum value of Lennard-Jones parameters of interaction pairs between HNT and PU, those optimum parameters minimize the difference between the DFT and molecular mechanic potential of all trajectories, for this purpose, the Particle Swarm Optimization (PSO) technique was used. 48 Kennedy and Eberhart 48 were the first to introduce the particle swarm optimization (PSO) algorithm. It is inspired by the swarm behavior of birds. The PSO algorithm is a metaheuristic algorithm that is appropriate to optimize nonlinear continuous functions. 48 During standard PSO, the new position of each particle is determined by a velocity function which is updated at each iteration. 49 The purpose of an optimization problem is to discover a variable described by a vector Y = [Y 1 Y 2 Y 3 .Y n ] which minimizes or maximizes depending on the purposed optimization target of the function f(Y). 49 The parameter vector Y is identified as the position vector of particle and n designates the number of parameters that defined as a problem. In our case, the function f(Y) is a fitness function or goal function, which is represented to minimize interaction energy as described in equation (5) and n represent to Lennard Jones coefficient of interaction atom types that have to be optimized. A MATLAB 50 implementation of the PSO algorithm was used to obtain Energyfitting as described in equation (5). During iteration, we assigned initial parameters as the table from PCFF-IFF 15 for the molecular mechanic potential, then interaction energy of MM arrangement was computed and compared with interaction energy of the same DFT arrangement. The coefficients e i and s i which are assigned to each Lennard Jones term were quickly obtained by linear least-squares fitting using the particle swarm optimization algorithm (PSO) method. PSO optimizes the coefficients while checking parameterization through the control of a bound (upper bound and lower bound) on each individual value in the set. This bound was set to be a 10% deviation from the original Lennard Jones coefficients of the PCFF-IFF force field. Table 1 presented previous experimental and DFT calculations results of lattice parameters of the unit cell of Halloysite with the current geometry optimization Discrepancies are due to the pseudopotential, cutoff energy, and k-point mesh density in the different algorithms. Table 2 describes the best-fit interaction energy between HNT/PU computed by DFT and MM approach. Figure 4 shows the fitting process results. In our case we chose a constant value of 2 for both the socalled ''social learning Coefficient'' and the ''best Individual learning coefficient'' of PSO. The number of PSO iteration was 150.

Results and discussion
Non-bonded energy terms were calculated using E ij = q i q j /r + e ij [2(s ij /r) 9 2 3(s ij /r) 6 ]. Where E ij is the non-bonded energy between atoms i and j, r is the distance between these atoms, q i is the charge on atom i and likewise for atom j. As shown in Figure 4, the modified parameters of Lennard-Jones are in good agreement between the MM calculation and corresponding DFT calculation with an accuracy of less than (5 kcal mol 21 ). During the PSO iteration, the conjugate gradient minimization 51 was performed at zero temperature to calculate the energy of systems in MM. This process was continued until the subsequent difference of energy between steps was less than 10 24 kcal/mol. Similarity during the DFT, the conjugate gradient minimization 51 was performed at zero temperature to calculate the energy of systems in each increment of movement. This process was continued until the subsequent difference in energy between steps was less than 0.25 kcal/mol.

Pull-Out Simulation
In nano-composites, the load transfer from the matrix to the nanoparticle occurs in the interface between matrix and nanoparticle. This depends on the cohesion strength between nanocomposite components. The interfacial interaction energy can determine the cohesion strength between matrix and nanoparticle. In the case of non-functionalized HNT-PU composites, the interfacial interaction energy is generated from van der Waals forces and Electrostatic forces.
MD was performed to simulate the pull-out of a single HNT from the PU matrix to determine the interface interaction energy for an HNT-PU nanocomposite. In this simulation, the PU was arranged in a random distribution around the HNT axis. The model applied in this section has a total of 102,095 atoms; The computational region was embedding inside a periodic structure. Periodic boundary conditions (PBC) were applied in all directions. The approach of Haghighatpanah et al. was adopted for the pull-out simulation, the method included the following steps 52 : Step 1: First, geometry optimization (GO) was conducted. GO was using the conjugate gradient 51 method and it was continued until the subsequent difference energy steps are less than 10 24 kcal/mol.
Step 2: The system obtained from Step 1 was equilibrated with an NpT ensemble at 1 atm pressure and 300 K for 100 ps. The HNT was constrained during the equilibration. The relaxation process was conducted to eliminate residual stress. 53 The equilibrium step was continued until the average total energy and density reach a constant value.
Step 3: After the equilibrium stage, the model was subjected to the pullout simulation. First, the PBC was extended along the HNT axis, then HNT was pulled out from the PU matrix along the HNT axial direction. The left end of the PU matrix was fixed as displayed in Figure 5(a). The HNT was pulled out of PU by a 0.005 Å displacement increment. During this process, we adopted the PPPM summation method 54 with the Lennard-Jones cutoff radius 18 Å . Process was applied to the nanocomposite under an NVE ensemble at 300 K temperature for 1.3 3 10 5 steps. Also, the time step was 1.0 fs and a total of 5 3 10 5 time steps were used. Figure 5 illustrates pull out simulation. The average interaction potential energy from the two systems is displayed in Figure 6.
Gou et al. 55 express the variational extraction energy as a function of the traveling distance of nanotube within the polymer. Figure 6 illustrates the variation of the interfacial energy of HNT-PU structure as a function of the extraction displacement of HNT. When the HNT is extracted from the polymer matrix, the interaction potential energy is increased which is directly related to the interfacial interaction. The interaction energy is increasing due to developing new interface areas. Figure 6 depicts the interaction potential energy of HNT-PU composite for the different increments of the HNT movement during the pulling out simulation. Once, the HNT is traveling out from the matrix, the potential energy reaches to the summit and with continued pulling out, the interaction energy of HNT with the matrix decreased in both simulations.

Model validation: Comparisons with the experimental result
We analyze the strain-stress of PU/HNT with 1 wt% of HNT by MD simulations and compare to an experimental sample. The amorphous models in this section were created with 1450 chains of PU. The system has 174,801 atoms. LAMMPS was used for all the MD simulations. 56 The updated potential is applied in this simulation. Following, 53 an equilibrium process was applied to the nanocomposite under an NPT ensemble at zero pressure and 500 K for 5 3 10 6 steps. This step allows the grafted chains to move in the periodic cell. During this process, we adopted the PPPM summation method 54 with the Lennard-Jones cutoff radius 18 Å . Also, the time step was 0.1 fs, and PBC was implemented along with all directions. The final density of the system was 0.93826 g/cm 23 .
After preparation, the PU/HNT structure is deformed by changing the cell length with constant speed (strain rate is 0.001 s 21 ) in the z-direction while the box length in the other two directions was fixed (as displayed in Figure 7). During this process, we adopted the PPPM summation method 54 with the Lennard-Jones cutoff radius 18 Å . Also, the time step was 1.0 fs, and the total number of time steps was 1.5 3 10 5 . The resulting of the stress-strain curve is depicted in Figure 8.   An experimental sample was prepared to conduct a uniaxial test. The experimental sample of PU was functionalized with an amino silane and the sample was mixed by loading HNT with PU at 0Á8 wt%. HNT was supplied by Sigma-Aldrich Canada Co. The chemical composition was H 4 Al 2 O 9 Si 2 Á 2H 2 O and the Molecular weight was 294.19 g/mol. Then a tensile test was conducted to measures stress versus relative elongation of the sample. Dusing the tensile test the specimens of HNT-PU was conducted with clamp pressure 400 psi, and the samples were stretched at 8.3 mm/s. Figure 8 depicts a comparison between MD simulations and experimental.
As displayed in the stress-strain result in Figure 8, there is the same agreement between MD results with the experimental result. However, the discrepancies between MD simulation and experimental are associated with accurate simulation of interatomic interaction, thus, this study was an attempt to eliminate them by PSO.

Conclusions
The development of an interatomic forcefield is a substantial subject for any molecular simulation to evaluate the dynamic response of molecular systems. The PCFF-IFF force field presents a set of interatomic potentials that enable study of polymeric materials. In this work, a new technique was developed for generating new interatomic coefficients for HNT-PU that can be simply obtainable from DFT-MM trajectories, by applying the particle swarm optimization (PSO) algorithm to fit MM calculation to the desirable DFT calculation. As a result, the optimal values of Lennard-Jones parameters were derived. This is a notable development for assembling standard force fields and DFT calculations. The new Lennard-Jones parameters present good agreement between MD stress-strain responses with the corresponding experimental result.
Finally, we studied the effect of changing the interatomic bonding between the HNT and PU composites on the strength of the composite. It is clear that weakening the HNT/PU bonding (to better match DFT results) leads to a reduction in composite strength, but we note that this effect is only a small contribution to the overall strength of the material.
The framework of this study, which considers an atomistic arrangement of HNT-PU interface to be described by a simple geometry of the HNT-PU, is extendible to various material arrangements especially design new composites reinforced by HNT.