Atomistic QM/MM simulations of the strength of covalent interfaces in carbon nanotubepolymer composites

We investigate the failure of carbon-nanotube/polymer composites by using a recently-developed hybrid quantum-mechanical/molecular-mechanical (QM/MM) approach to simulate nanotube pull-out from a cross-linked polyethene matrix. Our study focuses on the strength and failure modes of covalently-bonded nanotube–polymer interfaces based on amine, carbene and carboxyl functional groups and a [2+1] cycloaddition. We find that the choice of the functional group linking the polymer matrix to the nanotube determines the effective strength of the interface, which can be increased by up to 50% (up to the limit dictated by the strength of the polymer backbone itself) by choosing groups with higher interfacial binding energy. We rank the functional groups presented in this work based on the strength of the resulting interface and suggest broad guidelines for the rational design of nanotube functionalisation for nanotube–polymer composites.


Introduction
Functionalised carbon-nanotube/polymer composites (CNPC) are a promising structural material for demanding applications such as ballistic protection 1,2 and in the aerospace industry. The current generation of CNPCs, however, fall short of their theoretical limit 3,4 for mechanical properties, indicating that there is room for improving the efficiency of transferring external load to the interface between the carbon nanotubes (CNT) and the polymer matrix, and thereby more effectively realising the potential of CNT reinforcement. Key factors affecting the quality of load transfer in CNPCs include the dispersion 3,5-8 and alignment 3,6 of the fibres, as well as the strength of interfacial bonding. 7,[9][10][11][12][13][14] The latter is of critical importance 3,15 and is the focus of this work.
Experimental studies of CNPCs have shown that the quality of the interface can be significantly improved by the introduction of chemical cross-links between CNTs and the polymer matrix. CNT pull-out tests using atomic-force microscopy have shown that functionalisations that promote nanotube-matrix cross-links can improve the interfacial shear strength (ISS) from 40 MPa to 150 MPa 16 and similar trends have been observed by analysing CNPC Raman spectra 10,17,18 and in fracture tests. 19 The conclusions of these experiments are supported by computational investigations: Frankland et al. 20 have simulated CNT pull-out from crystalline polyethene (PE), finding that the ISS can be improved from 2.8 MPa to 110 MPa by introducing nanotube-polymer cross-links; Namilae et al. 21 have reported an increase of ISS from 12.5 MPa to 1000 MPa; Chowdhury et al. 22 report an ISS of 250-300 MPa for a crystalline CNT-PE interface and, while they do not provide a corresponding value for a nonbonded interface, their result can be contrasted with other CNT pull-out simulations which give ISS within the range of 1-150 MPa. [20][21][22][23][24][25][26][27][28][29][30][31][32][33] Whilst it has been shown that CNT functionalisation can have a significant impact on CNPC properties 3,15,34 and the effect of varying functionalisation has been studied extensively from the perspective of CNT dispersion in the polymer matrix, 7,8,35 less effort has been devoted to discriminating between various functionalisation strategies based on the strength of the resulting covalently-bonded interface, with most studies focusing on the comparison between bonded and nonbonded interfaces. 17,18,[20][21][22][23] Experimental studies concerned with development of high-performance functionalised CNTpolymer composites often rely on CNTs functionalised with carboxyl, [16][17][18][19][36][37][38][39] amine 10,38 or carbene 36,40,41 groups to create covalently bonded interfaces while better performing candidates are available.
It is important to note that the chemical composition of the interface is not the only factor affecting the strength of carbon nanotube-polymer connection and tuning combination of facets can lead to greatly improved composite performance. For CNPCs based on functionalized nanotubes, the ISS has been shown to increase with the increase in matrix density, 22 the crystallinity of the polymer 20 and functionalization density. 21,32 However, only weak dependence have been found between interfacial strength and CNT embedded length. 23 For composites with non-bonded interface, it has been found that the ISS is positively correlated with nanotube diameter 27,33 and only weakly coupled with CNT length. 23,27 Mechanical failure of the CNPC originates from local changes in chemical bonding that occur at the nanoscale but propagate via stress transfer through the polymer matrix at micrometre lengths. Accurately modelling the former requires a quantum-mechanical (QM) treatment of bond-breaking processes, but fully quantum-mechanical methods are too computationally expensive to apply on a routine basis to systems and processes at these length-and time-scales. This issue is often addressed through the use of classical molecular mechanical (MM) models based on parametrised empirical force-fields such as ReaxFF, 42 REBO 43 or AIREBO. 44 Such models can be calibrated to reproduce QM results, 45 but laborious refitting is required for each different chemistry and transferability to configurations not featured in the training set is not guaranteed. The limitations of such potentials for this particular problem were demonstrated in earlier work 46 where detachment of various functional groups from a CNT surface was simulated using ReaxFF and a fully QM approach. The results showed that bond-breaking processes were not described accurately with the reactive force-field, with significant qualitative differences as compared to the full QM simulations. Here, we use a concurrent hybrid QM/MM method 46 to tackle the multiscale problem of critical failure. In this approach, the system evolves according to a relatively simple MM model, but with regions where bond-breaking occurs identified on the fly and the atomic forces in these regions are computed using a QM approach. This method gave excellent agreement with the fully quantum calculations of the functional group detachment 46 at a fraction of computational cost. Similar QM/MM techniques have been used successfully to simulate crack propagation in silicon, 47,48 dislocation motion in metals 49 and molecular adsorption on a CNT surface. 50,51 In this work, we study how the chemical structure at the CNT-polymer interface determines its failure characteristics by simulating CNT pull-out from a cross-linked PE matrix. We focus on CNT-polymer interfaces composed of commonly used amine, carbene and carboxyl functional groups as well as a [2+1] cycloaddition. The systems studied differ only in the choice of the functional group linking the polymer matrix to the nanotube, allowing a direct comparison of the interfacial strength and failure characteristics of each functionalisation. Our findings allow us to make practical recommendations that will guide further experimental work. Additionally, the molecular resolution offered by our computational study enables us to analyse and explain the underlying interfacial failure mechanisms.
The remainder of this paper is structured as follows: our hybrid QM/MM approach is briefly summarised in Section 2 alongside a detailed description of the simulation setup; Section 3 is devoted to the analysis of the results of our nanotube fibre pull-out simulation; and we summarise our main conclusions in Section 4.

Methods
The single CNT pull-out simulations conducted in this study were carried out using the QM/MM hybrid method introduced and demonstrated by Gołębiowski et al. 46 In this approach, which builds on the 'Learn on the Fly' (LOTF) scheme originally introduced by De Vita and co-workers, 47,52 the focus is placed on ensuring forces for all atoms are sufficiently accurate to ensure accurate atomic trajectories and hence accurate measurement of derived observables. Here, this is achieved by combining forces computed with classical and QM simulations using the adaptive buffered force mixing variant of QM/MM, 53 rather than derived from a combined Hamiltonian as in the original LOTF method.
Classical force-fields are typically designed to accurately describe atomic systems in the vicinity of equilibrium configurations making them capable of describing the general morphology of the system. However, when the structure is significantly disturbed such as in the case of interfacial failure, the description can be incorrect. In the classical picture, local departure from the stable configuration is marked by elevated potential energies, therefore in our hybrid scheme particles with energies exceeding a pre-defined threshold are flagged for quantum-mechanical treatment; additional details about the flagging mechanism can be found in the ESI. † Marked atoms are used to form the inner QM clusters by selecting all particles within a set number of bond hops of each flagged particle and expanded by including a buffer region with a specific width (also defined in terms of bond hops) to form full clusters. The resulting full clusters are carved out of the system, isolated and have all dangling bonds passivated by hydrogenation to prepare structures for quantum mechanical calculations. After the QM computation, classical forces on the inner cluster atoms are replaced by the quantum forces computed for those particles to form a complete QM/MM description and the total sum of forces is adjusted to zero to conserve momentum.
We used the generalised Amber force-field (GAFF), 54 as implemented in LAMMPS, 55 as our MM model. The standard GAFF model was refitted for bond and angle terms in the CNT and the polyethene chain, in order to match the lattice constant and elastic response around the equilibrium obtained with the force-field and the QM method for a (10,0) CNT and a 25 monomer PE chain as described by Gołębiowski et al. 46 The QM method used is density-functional tight binding (DFTB), 56 as implemented in the PLATO 57 simulation package. DFTB calculations were performed with self-consistent charge equilibration with an electronic temperature of 0.2 eV and Pulay mixing, using the parameters by Krishnapriyan et al. 58 optimised for atomization energy and interatomic forces in C, H, N, and O system. In this work we do not use the electrostatic embedding approach for our quantum calculations but rather use a mechanical embedding technique and rely on a wide buffer to ensure an accurate description of forces in the core of the QM regions, following the discussion in Peguiron et al.; 59 see Gołębiowski et al. 46 for more details. Throughout the pull-out simulations, an NVT ensemble was enforced using the Langevin thermostat with the friction coefficient of 0.01 as implemented in the atomic simulation environment. 60 The simulations were carried out at a temperature of 100 K. The effect of changing the temperature from 100 K to 300 K was found to be negligible (see ESI † for details). DFTB calculations do not include a dispersion term and the forces on the atoms of the inner clusters were adjusted post hoc by including dispersion interactions with respect to all of the particles in the system as discussed by Zhechkov et al. 61 The correction was obtained from the Lennard-Jones term of the GAFF, as described in the ESI. †

Simulation setup
A CNT embedded in a cross-linked PE matrix was gradually pulled out of the matrix and the response of the system was analysed. The preparation of the simulation cell is described briefly below, with further details provided in the ESI. † The initial polymeric structure was a slab with a thickness of 50 Å and density equal to the target density of 0.3 g cm À3 . The slab was composed of fifty 48-monomer polyethene-like chains, each made out of six alternating units of propylene, acetylene and propylene organised in the following fashion: H-(-CH 2 -CH 2 -CH 2 -CHQCH-CH 2 -CH 2 -CH 2 -) 6 -H. The initial chain configuration was obtained using a Monte Carlo sampling procedure as implemented in the MedeA software package. 62 The next step involved creating a void inside the polymer cell by introducing seven Lennard-Jones (LJ) repulsive beads spaced equally between the centre of the slab's free surface and a point inside the polymeric structure located at the midpoints of the periodic axis and 35 Å away from the surface in the nonperiodic axis. The strength of LJ repulsion was increased quasistatically by changing the depth of the potential well from 4 meV to 200 meV and the equilibrium distance from 1 Å to 8 Å over the course of ten iterations, each including parameters change and geometry optimisation with the FIRE algorithm. 63 Finally the values were kept constant for 0.1 ns in an NVT simulation at 300 K to allow for local equilibration, after which the beads were removed from the matrix.
Once the process was completed, a 30 Å long, (10,0) carbon nanotube fragment with hydrogen passivated ends was placed inside the newly created void; its position was chosen so that it aligned with the repulsive beads resulting in one end of the CNT inside the slab and the other end aligned with the slab free surface. CNT functionalisation was then achieved by introducing three additional PE chains, identical to those used to construct the polymer bulk, which were attached to the CNT surface using amine, carbene or carboxyl functional groups or a [2+1] cycloaddition resulting in a functionalisation density of one group per 100 CNT carbon atoms. Schematic representations of the four covalently-bonded CNT-polymer interfaces is shown in Fig. 1. After functionalisation, the system was equilibrated using the procedure of Kucukpinar and Doruker, 64 which involves compression to 0.5 GPa over 0.3 ns and relaxation over an additional 0.3 ns, then heating the sample to 800 K over 0.5 ns and subsequent cooling to 300 K over the same period and a final compression-relaxation cycle identical to the first one. Finally, cross-linking of the polymer network was carried out by introducing atomic bonds between all carbon atoms bonded with three ions closer than 3.5 Å during a 0.5 ns simulation in an NVT ensemble at 800 K to generate a fully connected polymer network. All simulations in the abovedescribed process were performed with the GAFF, 54 as implemented in LAMMPS 55 with the timestep of 1 fs. We believe that a fully classical model is sufficiently accurate to identify generally correct morphology.
The procedure described above yields a cross-linked PE matrix reinforced with a 30 Å functionalised CNT that accounts for 12% of the total weight (see Fig. 2). The average density of the resulting cross-linked PE was found to be approximately 0.945 g cm À3 , which is consistent with experimental findings, 65 and the cross-link density of the polymer network, defined as a ratio of cross-linked carbon atoms to total carbon atoms was found to be approximately 10%.
CNT pull-out was simulated with the QM/MM method by extracting the CNT from the polymer slab at a constant velocity while analysing the response of the system. After equilibration for 10 ps, the terminating ring of carbon atoms aligned with the composite slab free surface was constrained to move outwards of the slab with a velocity which was gradually increased over the course of 1 ps from 0 Å fs À1 to 10 À3 Å fs À1 (equivalent to the strain rates that could be observed in ballistic tests 1,2 ). Such pull-out velocity is high when compared to experimental pullout tests but is comparable to velocities used in computational studies with classical potentials 21,23,24 that do not the include the accurate and computationally expensive quantum mechanical description of bond-breaking. The constrained CNT atoms are marked in green in Fig. 2. Following the velocity ramp-up, the pull-out simulation was carried out by keeping the pull-out velocity constant and recording the forces acting on the constrained ring of atoms.

Results and discussion
For each attachment chemistry shown in Fig. 1, the initial structure was prepared as discussed in Section 2. The simulation procedure including generating the initial atomic structure and simulating the CNT pull-out has been performed for an ensemble of ten independent configurations for each functionalisation chemistry.
For each simulation, the component parallel to the pull-out direction of the force acting on the constrained ring of atoms was recorded at each timestep as the pull-out force and plotted as a function of CNT displacement. The resulting curve was smoothed using a 100 fs time window with a top-hat smoothing filter to average out high-frequency vibrations associated with internal modes of the CNT. The smoothed curve displays local peaks that are associated with bond-breaking events; we use the value of the force at each of these peaks (which we refer to as the ''critical force'') as an indicator of the strength of the interface. The method for determining the exact location of each bond-breaking event is discussed in more detail in the ESI. † The top panel of Fig. 3 shows a typical example of the smoothed force-displacement curve, with peaks corresponding to failure marked with stars; there are three such events per simulation as each system is prepared with three functionalised attachments between the CNT and the polymer matrix. The bottom panel of the same figure shows the number of atoms flagged for QM calculation at each point of the simulation. It can be seen that the more computationally-demanding QM calculations are only performed when necessary, i.e., just before, during and just after each series of bond rearrangements and that atoms are promptly de-flagged following a failure event. It can also be seen that at any point during the simulation, the number of atoms for which QM forces are calculated is a small fraction of the total number of atoms in the system (typically less than 1% of the total number of atoms at any given time).
The critical forces recorded from each pull-out simulation have been aggregated and the distributions describing the  strength of the interface in each case forms the key result of this work: the histograms and aggregate statistics are shown in Fig. 4. We find that interfaces based on amine (NH), carbene (CH 2 ) and carboxyl (COOH) functional groups are characterised by similar strength; there is, however, an improvement in ISS of approximately 50% when the [2+1] cycloaddition (CH) is used to attach the polymer to the CNT. The first three functionalisations rely on a single covalent bond to graft PE chains to the CNT, whereas the cycloaddition is grafted with a pair of single bonds. The systems that we have investigated are identical except for the chemistry of the attachments, indicating that the functionalisation can indeed be chosen to enhance the ISS of the composite. In their study, Milowska et al. 66 have shown that the first group of attachments (NH, CH 2 and CH) have similar binding energies of approximately 1.3 eV while the [2+1] cycloaddition are characterised by a higher value of approximately 2.3 eV. This indicates that the binding energy of a functional group can be seen as a proxy for the critical force of the attachment and therefore used to gauge the strength of the resulting interface.
For each pull-out simulation, the area under the pull-out force-CNT displacement curve provides a measure of the interfacial fracture energy; the distribution of fracture energy across the ensemble of ten simulations for each functional group are portrayed in Fig. 5. The results show that interfaces based on amine, carbene and carboxyl functional groups have similar fracture energy; however, an improvement of approximately 70% is found when a [2+1] cycloaddition is used to anchor the polymer chains to the CNT. The difference between fracture energy values of the carboxyl-based interface and amine-or carbene-based ones was analysed with a two-sided t-test for the null hypothesis that two independent samples (with the COOH-based interface as the first sample and combined results for NH and CH 2 interfaces as the second sample) have identical expected values. The test concluded with a p-value of 0.037. In comparison, a similar test performed with amine, carbene and carboxyl interfaces in one sample and the [2+1] cycloaddition-based connection as the second resulted in a p-value of 5 Â 10 À7 , reinforcing the conclusion that our findings can be used to clearly divide the interfaces into two major classes with significantly different fracture energy performance.
The mechanism of interfacial failure has been further studied by inspecting the yield strength and atomistic structure of individual attachments in each of the composite systems. The analysis reveals that in systems where CH 2 , NH and COOH functional groups were used to form the interface, critical failure occurred by the failure of the functional group and detachment of polymer chains from the CNT surface demonstrating that the grafting point itself is the weakest element of the interface. On the other hand, in the case of [2+1] cycloadditions (CH), two equally common failure modes were observed: detachment of the functional group from the CNT and bond-breaking within the polymer matrix. The distribution of the strength of individual attachments, discriminated by their failure mode is shown in Fig. 6. The [2+1] cycloadditionbased interface and the polymer backbone fail at a similar load, demonstrating that their strengths are comparable and a further increase of the strength of the interface would not yield any effective improvement without fortifying the polymer matrix itself. Our results show that the functionalisation chemistry can be tuned to improve the ISS up to a limit defined by the strength of the polymer backbone.
Additional effort has been devoted to the analysis of the [2+1] cycloadditions (CH) based interface by exploring the effect of the attachment point geometry on the yield strength. Due to the structure of carbon nanotubes, [2+1] cycloaddition can be found in two configurations, attached to two CNT atoms located parallel to the CNT axis or at an angle of 451. Fig. 7 shows  the normalised distribution of the strength of individual attachments divided into these two configurations. It can be seen that the orientation of the attachments is not a major factor for the designing of high-strength CNT-polymer interfaces based on functional groups attached with two bonds.
Quantitative comparison of our calculated ISS with experimental measurements is very challenging for systems as complex as real CNT-polymer composites. The pull-out velocity commonly used in computational pull-out studies 21,23,24 is significantly higher than is experimentally achievable and limits the ability for the polymer to respond to the CNT motion. This is exacerbated by the differences in CNT diameters and lengths which vary by orders of magnitude between simulated and experimental systems. However, we have compared our results to the experimental work of Barber et al. 16 who performed pull-out experiments of carboxyl-functionalised CNTs. The experiments of Barber et al. 16 demonstrate that the pull-out force is proportional to the embedded surface area of the CNT and recorded the ISS (maximum pullout force over CNT area) of around 50 MPa. For comparison, our simulations give a value of around 1000 MPa, i.e., a factor of 20 larger. Putting aside the difference in chemistry, pull-out velocity and the size § between the experimental system and our model, a crucial piece of information is missing to enable us to relate these two values, namely the density of covalent attachments between the CNT and the polymer matrix in the experiment. In our model we have three attachments over the surface area of a 300 carbon atoms CNT, but the equivalent grafting ratio is not known for the experimental system in Barber et al., 16 and this makes it very difficult to make a reasonable attempt at a comparison. Other experimental studies concerned with processing of carbon nanotubes-polymer composites have measured grafting ratios between a factor of 1.5 to a factor of 30 smaller than in our study. 67 Such differences in grafting ratios are likely to be a significant determinant of the discrepancy between our simulated values of F max /A and the experimental measurements of Barber et al. 16

Conclusions
We investigated the critical failure of covalently bonded CNTpolymer interfaces by simulating pull-out of a CNT fibre from a cross-linked polyethene matrix. We focused on the strength and failure modes of composite interfaces composed of polymer chains grafted to the CNT surface by amine, carbene and carboxyl functional groups and via [2+1] cycloaddition. We used a hybrid QM/MM simulation method whereby the system evolves according to a relatively simple classical force-field, but regions that stray outside of the scope of accuracy of the force-field (e.g., where bond-breaking events are about to occur) are identified on the fly as the dynamical simulation progresses and the atomic forces in these regions are computed using a more accurate QM approach.
We found that interfaces based on CH 2 , NH and COOH functional groups exhibit similar strength, suggesting that other factors, such as CNT dispersion, should take precedence when discriminating between them. With the [2+1] cycloaddition, however, a significant improvement in interfacial shear strength and fracture energy was observed, with the resulting interface found to be 50% stronger and 70% tougher than the other cases studied. The [2+1] cycloaddition is characterised by a higher binding energy than the other functional groups, 66 indicating a  The total number of sideways attachments is 24 while the number of parallel ones is 6, therefore the counts are normalised for ease of comparison. The mean and standard deviation of the results are shown in the legend. The bins used in the histogram have a width of 1.5 eV Å À1 with the left-most edge of the first bin located at 0; the data points are placed at the midpoint of each bin. § The authors of Barber et al. 16 estimate their CNT diameters to be 80 AE 30 nm and their lengths to be 5-10 mm, whereas in our model we have a diameter that is two orders of magnitude smaller (0.8 nm) and a length that is over three orders of magnitude smaller (3 nm).
correlation between the functional group binding energy and the resulting ISS. Furthermore, in the case of the [2+1] cycloaddition, we observed two distinct failure modes in our simulations: detachment of the functional group from the CNT, and rupture of the polymer backbone inside the matrix. Both were found to occur at similar critical strengths, demonstrating that the grafting strength is on a par with the strength of the polymer backbone for the [2+1] cycloaddition-based interfaces studied here, effectively maximising the ISS for the composite. This is the first time this on-the-fly hybrid QM/MM approach has been used to simulate failure in CNT pull-out at time-and length-scales usually only accessible to purely classical methods, and expands the simulation toolkit available for studying these inherently multiscale systems.

Conflicts of interest
The authors declare no competing financial or non-financial interests.