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

Abstract

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.

Keywords

Halloysite nanotubes molecular dynamic optimized force field particle swarm optimization pulling out

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₂Si₂O₅(OH)₄·2H₂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.

Figure 1.

Schematic of nanosheet and nano scroll form of HNT.¹

Most research on Halloysite has concentrated on measured mechanical properties under different loading.² Lu et al.³ experimentally measured the elastic modulus. They found that the elastic modulus of the HNT layer is 130 ± 24 GPa. Lvov et al.¹ calculated Young’s modulus of a multi-wall alumino-silicate nanotube. They obtained Young’s modulus in the range of (230–340 GPa). Guimarã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.⁴ 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.⁵

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.⁶ PU has a variety of uses in a range of products.⁷ Furthermore, polyurethane has good corrosion resistance and toughness, which is effective functionality for the defense industries.⁸ Recently, new applications of polyurethane have been developed for impact protection of structures, including anti-explosion and ballistic loading.⁹ 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.

Figure 2.

Schematic of polyurethane.⁶

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–13

The total energy of a molecular structure is giving by;

\begin{matrix} U = \sum_{all bonds} \frac{1}{2} K_{b} {(b - b_{0})}^{2} + \sum_{all angles} \frac{1}{2} K_{θ} {(θ - θ_{0})}^{2} \\ + \sum_{all torsions} K_{φ} [1 - \cos (n φ)] + \sum_{χ} K_{χ} χ^{2} \\ + \sum_{i, j nonbonded} ε_{ij} [2 {(\frac{σ_{ij}}{r_{ij}})}^{9} - 3 {(\frac{σ_{ij}}{r_{ij}})}^{6}] \\ + \sum_{i, j nonbonded} \frac{q_{i} q_{j}}{4 π ε_{0} ε r_{ij}} \end{matrix}

(1)

The first four terms in equation (1) introduce the stretched bond energy (b), bond angles energy (θ), torsion angle energy (ϕ), and distortion energy respect to a reference plane (χ). 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¹⁶ studied the hardness of PU-HNT nanocomposites. They indicate that the scratch resistance has been improved by adding 5 wt% HNT. Gaaz et al¹⁷ 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.²⁰ Those ML models have been developed for small atomic configurations, which is not suitable for large atomic configurations.²¹ 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–24 Energy-fitting utilizes the atomic energy obtained from a high-accurate DFT simulation. Ercolessi and Adams,²⁵ introduced force-fitting for the first time, and it was employed in numerous simulations^26–29 to improve the efficiency of semi-empirical quantum methods.^30,31 Wu et al.³² optimized parameters of a potential field for long-range interaction through a ceteris paribus optimization method, and Daven et al.³³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³⁴), or global minimum explorations (e.g. simulated annealing³⁵), 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–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–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:

The optimized lattice parameter of simple unit cell HNT is found.

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.

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.

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).

We repeat step 4, this time with PCFF-IFF. Then, the interaction energy is calculated for different vertical positions of the PU.

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.

Figure 3.

(a, c) Nearest PU-HNT configuration and (b, d) farthest PU-HNT configuration.

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⁴² with the PBE gradient-corrected functionals⁴³ methods. The calculation was performed at zero temperature. Reciprocal space was represented by the Gamma k-points scheme with 5 × 5 × 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⁴⁴ was applied for describing the core electron interactions. The result is provided in Table 1. Then a supercell was generated from 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.

Table 1.

Computed and experimental lattice parameters of geometry optimization; parentheses are shown percentage deviations.

	DFT⁴⁵	DFT⁴⁶	Experimental⁴⁷	DFT
	B3LYP	B3LYP	Experimental⁴⁷	This work PBE
a (Å)	5.184 (0.3)	5.230 (1.2)	5.17	5.11314
b (Å)	7.396 (0.1)	7.480 (1.1)	7.40	7.30867
c (Å)	5.190 (0.8)	5.210 (1.2)	5.15–8.94	5.45623
α (°)	83.50 (−0.9)	84.10 (0.2)	84.23–91.69	84.5556
β (°)	99.80 (1.1)	98.80 (0.1)	98.70–104.61	99.9974
γ (°)	59.90 (0.0)	59.92 (89.2)	–	61.0147

Data are reported in Å and °.

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;

{E_{int eract}}^{DFT} = E_{tot - DFT} - (E_{PU} + E_{HNT})

(2)

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)”⁴² with the PBE gradient-corrected functionals⁴³ 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 k-points 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 = −4188.717 kcal/mol and energy for the flat surface of Halloysite is E_HNT = −50,656.40 kcal/mol. The result of interaction energy is presented in Figure 4.

Figure 4.

Interaction energy versus position of PU chain: (a) before and (b) after the fitting process.

The interaction energy calculations (MM approach)

For MM calculations, the same supercell was used (described in Figure 3). Then a minimizing process with 10⁻⁴ kcal/mol discrepancy was applied to the interaction between HNT/PU. The total energy between the PU and the HNT was calculated by the last two terms of equation (1):

{E_{int eract}}^{MM} = E_{tot - MM} - (E_{PU} + E_{HNT})

(3)

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_{interact}^{MM}$ indicates interaction energy between HNT and PU molecules.

In the classical molecular dynamic method, the non-bonded energy between the HNT and the PU is calculated with the 9-6 Lennard-Jones potential and Coulomb part as described in equation (4);

\begin{matrix} E_{Interact}^{MM} = \sum_{i, j nonbonded} ε_{ij} [2 {(\frac{σ_{ij}}{r_{ij}})}^{9} - 3 {(\frac{σ_{ij}}{r_{ij}})}^{6}] \\ + \sum_{i, j nonbonded} \frac{q_{i} q_{j}}{4 π ε_{0} ε r_{ij}} \end{matrix}

(4)

Where $E_{Interact}^{MM}$ is the non-bonding potential energy, r_ij is the gap between a pair of atoms, è identify the well depth of potential, σ 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₁₇N₃O₃H₁₈) 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 × 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):

f (Y) = \sum_{i} {({(E_{int eract}^{DFT})}_{i} - {(E_{int eract}^{MM})}_{i})}^{2}

(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.⁴⁸

Kennedy and Eberhart⁴⁸ 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.⁴⁸ During standard PSO, the new position of each particle is determined by a velocity function which is updated at each iteration.⁴⁹

The purpose of an optimization problem is to discover a variable described by a vector Y = [Y₁Y₂Y₃…Y_n] which minimizes or maximizes depending on the purposed optimization target of the function f(Y).⁴⁹ 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⁵⁰ implementation of the PSO algorithm was used to obtain Energy-fitting as described in equation (5). During iteration, we assigned initial parameters as the table from PCFF-IFF¹⁵ 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 $ε$ _i and $σ$ _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.

Results and discussion

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, cut-off 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 so-called “social learning Coefficient” and the “best Individual learning coefficient” of PSO. The number of PSO iteration was 150.

Table 2.

Improved pair coefficients for classical simulations.

Atom	Al	O	H
C
ε_ij (kcal·mol⁻¹)	0.033836	0.028527	0.001001
σ_ij (A°)	4.222966	3.892939	3.501295
N
ε_ij (kcal·mol⁻¹)	0.055334	0.038907	0.001428
σ_ij (A°)	4.419411	3.995249	3.77849
O
ε_ij (kcal·mol⁻¹)	0.064659	0.057637	0.003619
σ_ij (A°)	3.778643	3.642636	4.338017
H
ε_ij (kcal·mol⁻¹)	0.014066	0.015391	0.001482
σ_ij (A°)	4.268376	3.456217	2.469121

Non-bonded energy terms were calculated using E_ij = q_iq_j/r + ε_ij[2(σ_ij/r)⁹ − 3(σ_ij/r)⁶]. 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⁻¹). During the PSO iteration, the conjugate gradient minimization⁵¹ 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⁻⁴ kcal/mol. Similarity during the DFT, the conjugate gradient minimization⁵¹ 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⁵²:

Step 1: First, geometry optimization (GO) was conducted. GO was using the conjugate gradient⁵¹ method and it was continued until the subsequent difference energy steps are less than 10⁻⁴ 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.⁵³ 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⁵⁴ with the Lennard-Jones cutoff radius 18 Å. Process was applied to the nanocomposite under an NVE ensemble at 300 K temperature for 1.3 × 10⁵ steps. Also, the time step was 1.0 fs and a total of 5 × 10⁵ time steps were used. Figure 5 illustrates pull out simulation. The average interaction potential energy from the two systems is displayed in Figure 6.

Figure 5.

Pull-out of the HNT from the PU matrix: (a) fully embedded, (b) at 60 Å displacement of the HNT, (c) fully extracted HNT from PU, and (d) cross-cut of fully extracted HNT.

Figure 6.

Interaction potential energy for HNT/PU during the pull-out test.

Gou et al.⁵⁵ 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.⁵⁶ The updated potential is applied in this simulation. Following,⁵³ an equilibrium process was applied to the nanocomposite under an NPT ensemble at zero pressure and 500 K for 5 × 10⁶ steps. This step allows the grafted chains to move in the periodic cell. During this process, we adopted the PPPM summation method⁵⁴ 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⁻³.

After preparation, the PU/HNT structure is deformed by changing the cell length with constant speed (strain rate is 0.001 s⁻¹) 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⁵⁴ 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 × 10⁵. The resulting of the stress-strain curve is depicted in Figure 8.

Figure 7.

The molecular dynamics models of PU/HNT composite under the stretching test.

Figure 8.

Stress-strain diagram for MD and experimental PU/HNTs samples.

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₄Al₂O₉Si₂· 2H₂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.

Footnotes

Acknowledgements

We are indebted to our collaborators, friends, and colleagues for their many invaluable discussions, comments, and suggestions. In particular, we want to explicitly thank Rafaela Aguiar for preparing the experimental sample test. We are also grateful for the financial support by the NSERC Discovery Grant program, and Compute Canada for providing computing resources for this research.

Handling Editor: James Baldwin

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

ORCID iDs

Mohammadreza Heidari Pebdani

Ronald E Miller

Supplemental material

Supplemental material for this article is available online.

References

Lvov

Shchukin

Möhwald

, et al. Halloysite clay nanotubes for controlled release of protective agents. ACS Nano 2008; 2: 814–820.

Lecouvet

Horion

D’Haese

, et al. Elastic modulus of halloysite nanotubes. Nanotechnol 2013; 24: 105704.

Chen

, et al. Direct measurements of the Young’s modulus of a single halloysite nanotube using a transmission electron microscope with a bending stage. J Nanosci Nanotechnol 2011; 11: 7789–7793.

Guimarães

Enyashin

Seifert

, et al. Structural, electronic, and mechanical properties of single-walled halloysite nanotube models. J Phys Chem B 2010; 114: 11358–11363.

Presti

Pedone

Mancini

, et al. Insights into structural and dynamical features of water at halloysite interfaces probed by DFT and classical molecular dynamics simulations. Phys Chem Chem Phys 2016; 18: 2164–2174.

Njoroge

. Atomistic simulation of graphene-polyurethane nanocomposite for use in ballistic applications. Doctoral Dissertation, Texas A & M University, 2013.

Petrović

Ferguson

. Polyurethane elastomers. Prog Polym Sci 1991; 16: 695–836.

Roland

Twigg

, et al. High strain rate mechanical behavior of polyurea. Polymer 2007; 48: 574–578.

Chattopadhyay

Raju

. Structural engineering of polyurethane coatings for high performance applications. Prog Polym Sci 2007; 32: 352–418.

10.

Chmiela

Tkatchenko

Sauceda

, et al. Machine learning of accurate energy-conserving molecular force fields. Sci Adv 2017; 3: e1603015.

11.

Lubbers

Smith

Barros

. Hierarchical modeling of molecular energies using a deep neural network. J Chem Phys 2018; 148: 241715.

12.

Sun

. COMPASS: an ab initio force-field optimized for condensed-phase applications overview with details on Alkane and Benzene compounds. J Phys Chem B 1998; 102: 7338–7364.

13.

Crozier

. A brief survey of the LAMMPS particle simulation code: introduction case studies and future development. Sandia National Lab.(SNL-NM), Albuquerque, NM, 2011.

14.

Heinz

Lin

T-J

Mishra

, et al. Thermodynamically consistent force fields for the assembly of inorganic, organic, and biological nanostructures: the INTERFACE force field. Langmuir 2013; 29: 1754–1765.

15.

Heinz

Vaia

Koerner

, et al. Photoisomerization of azobenzene grafted to layered silicates: simulation and experimental challenges. Chem Mater 2008; 20: 6444–6456.

16.

Nikiforow

Pohl

, et al. Polyurethane coatings reinforced by halloysite nanotubes. Coatings 2013; 3: 16–25.

17.

Gaaz

Sulong

Kadhum

AAH

, et al. Absolute variation of the mechanical characteristics of halloysite reinforced polyurethane nanocomposites complemented by Taguchi and ANOVA approaches. Results Phys 2017; 7: 3287–3300.

18.

Muñoz-Muñoz

Guevara-Carrion

Llano-Restrepo

, et al. Lennard-Jones force field parameters for cyclic alkanes from cyclopropane to cyclohexane. Fluid Phase Equilib 2015; 404: 150–160.

19.

Guidez

Bertoni

, et al. Perspective: ab initioforce field methods derived from quantum mechanics. J Chem Phys 2018; 148: 090901.

20.

Behler

Parrinello

. Generalized neural-network representation of high-dimensional potential-energy surfaces. Phys Rev Lett 2007; 98: 146401.

21.

Yao

Herr

Toth

, et al. The TensorMol-0.1 model chemistry: a neural network augmented with long-range physics. Chem Sci 2018; 9: 2261–2269.

22.

Faber

Hutchison

Huang

, et al. Prediction errors of molecular machine learning models lower than hybrid DFT error. J Chem Theory Comput 2017; 13: 5255–5264.

23.

Hansen

Biegler

Ramakrishnan

, et al. Machine learning predictions of molecular properties: accurate many-body potentials and nonlocality in chemical space. J Phys Chem Lett 2015; 6: 2326–2331.

24.

Behler

. Erratum: “Perspective: Machine learning potentials for atomistic simulations” [J. Chem. Phys. 145, 170901 (2016)]. J Chem Phys 2016; 145: 219901.

25.

Ercolessi

Adams

. Interatomic potentials from first-principles calculations: the force-matching method. EPL 1994; 26: 583–588.

26.

Pinnick

Calderon

Rusnak

, et al. Achieving fast convergence of ab initio free energy perturbation calculations with the adaptive force-matching method. Theor Chem Acc 2012; 131: 1–11.

27.

Wei

Song

Wang

. A simple molecular mechanics potential for μm scale graphene simulations from the adaptive force matching method. J Chem Phys 2011; 134: 184704.

28.

Akin-Ojo

Wang

. The quest for the best nonpolarizable water model from the adaptive force matching method. J Comput Chem 2011; 32: 453–462.

29.

Izvekov

Parrinello

Burnham

, et al. Effective force fields for condensed phase systems from ab initio molecular dynamics simulation: a new method for force-matching. J Chem Phys 2004; 120: 10896–10913.

30.

Goldman

. Multi-center semi-empirical quantum models for carbon under extreme thermodynamic conditions. Chem Phys Lett 2015; 622: 128–136.

31.

Welborn

Chen

Wang

, et al. Why many semiempirical molecular orbital theories fail for liquid water and how to fix them. J Comput Chem 2015; 36: 934–939.

32.

Derrick

Nahin

, et al. Optimization of long range potential interaction parameters in ion mobility spectrometry. J Chem Phys 2018; 148: 074102.

33.

Daven

Tit

Morris

, et al. Structural optimization of Lennard-Jones clusters by a genetic algorithm. Chem Phys Lett 1996; 256: 195–200.

34.

Levenberg

. A method for the solution of certain non-linear problems in least squares. Q Appl Math 1944; 2: 164–168.

35.

van Laarhoven

PJM

Aarts

EHL

. Simulated annealing: theory and applications. Dordrecht, The Netherlands: Reidel, 1987.

36.

Odegard

Gates

Wise

, et al. Constitutive modeling of nanotube–reinforced polymer composites. Compos Sci Technol 2003; 63: 1671–1687.

37.

Schneider

Hamaekers

Chill

, et al. ATK-ForceField: a new generation molecular dynamics software package. Model Simul Mater Sci Eng 2017; 25: 085007.

38.

Frankland

SJV

Harik

Odegard

, et al. The stress–strain behavior of polymer–nanotube composites from molecular dynamics simulation. Compos Sci Technol 2003; 63: 1655–1661.

39.

Liu

Chen

. Evaluations of the effective material properties of carbon nanotube-based composites using a nanoscale representative volume element. Mech Mater 2003; 35: 69–81.

40.

Seidel

Lagoudas

. Micromechanical analysis of the effective elastic properties of carbon nanotube reinforced composites. Mech Mater 2006; 38: 884–907.

41.

Awasthi

Lagoudas

Hammerand

. Modeling of graphene–polymer interfacial mechanical behavior using molecular dynamics. Model Simul Mater Sci Eng 2009; 17: 15002.

42.

Kresse

Joubert

. From ultrasoft pseudopotentials to the projector augmented-wave method. Phys Rev B 1999; 59: 1758–1775.

43.

Perdew

Burke

Ernzerhof

. Generalized gradient approximation made simple. Phys Rev Lett 1996; 77: 3865–3868.

44.

Blöchl

. Projector augmented-wave method. Phys Rev B 1994; 50: 17953–17979.

45.

Ugliengo

Zicovich-Wilson

Tosoni

, et al. Role of dispersive interactions in layered materials: a periodic B3LYP and B3LYP-D* study of Mg(OH)2, Ca(OH)2 and kaolinite. J Mater Chem 2009; 19: 2564–2572.

46.

Tosoni

Doll

Ugliengo

. Hydrogen bond in layered materials: structural and vibrational properties of kaolinite by a periodic B3LYP approach. Chem Mater 2006; 18: 2135–2143.

47.

Neder

Burghammer

Grasl

, et al. Refinement of the kaolinite structure from single-crystal synchrotron data. Clays Clay Miner 1999; 47: 487–494.

48.

Kennedy

Eberhart

. Particle swarm optimization. In: Proceedings of ICNN’95-international conference on neural networks, 1995, vol. 4, pp.1942–1948. New York, NY: IEEE.

49.

Wang

Tan

Liu

. Particle swarm optimization algorithm: an overview. Soft Comput 2018; 22: 387–408.

50.

Symbolic Math Toolbox. Natick, Massachusetts, United States: Matlab. Mathworks Inc., 1993.

51.

Fletcher

Reeves

. Function minimization by conjugate gradients. Comput J 1964; 7: 149–154.

52.

Haghighatpanah

Bohlén

Bolton

. Molecular level computational studies of polyethylene and polyacrylonitrile composites containing single walled carbon nanotubes: effect of carboxylic acid functionalization on nanotube-polymer interfacial properties. Front Chem 2014; 2: 74.

53.

Peng

Ranganathan

Keblinski

, et al. Viscoelastic and dynamic properties of polymer grafted nanocomposites with high glass transition temperature graft chains. J Appl Phys 2019; 126: 195102.

54.

Hockney

Eastwood

. Computer simulation using particles. Boca Raton: CRC Press, 1988.

55.

Gou

Jiang

Minaie

, et al. Nanoscale modeling and simulation of interfacial bonding of single-walled nanotube reinforced composites. ASME Int Mech Eng Congress Expo 2003; 37149: 815–819.

56.

Plimpton

. Fast parallel algorithms for short-range molecular dynamics. J Comput Phys 1995; 117: 1–19.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB