Computational investigation of interfacial failure in composite materials is challenging because it is inherently multi-scale: the bond-breaking processes that occur at the covalently bonded interface and initiate failure involve quantum mechanical phenomena, yet the mechanisms by which external stresses are transferred through the matrix occur on length and time scales far in excess of anything that can be simulated quantum mechanically. In this work, we demonstrate and validate an adaptive quantum mechanics (QM)/molecular mechanics simulation method that can be used to address these issues and apply it to study critical failure at a covalently bonded carbon nanotube (CNT)-polymer interface. In this hybrid approach, the majority of the system is simulated with a classical forcefield, while areas of particular interest are identified on-the-fly and atomic forces in those regions are updated based on QM calculations. We demonstrate that the hybrid method results are in excellent agreement with fully QM benchmark simulations and offers qualitative insights missing from classical simulations. We use the hybrid approach to show how the chemical structure at the CNT-polymer interface determines its strength, and we propose candidate chemistries to guide further experimental work in this area.

The excellent mechanical properties of carbon fibre-polymer composites (CFPCs) make them remarkable and well-established structural materials. Despite their success, CFPCs may be improved by using alternative reinforcement such as carbon nanotubes (CNTs), which offer similar advantages to carbon fibres but their better mechanical properties translate to enhanced structural and functional characteristics of the composite.1–4 Current applications of carbon nanotube-polymer composites (CNPCs) involve mostly lab-based experimental systems;5–8 however, increased understanding and availability could make them an industrially feasible choice, replacing CFPCs in demanding applications such as ballistic protection8 or aerospace, where high-strength yet lightweight components are needed.

Composite reinforcement can be improved by increasing the quality of dispersion3,9–11 and alignment of the fibres,3,10 as well as the strength of the interface.11–17 Given the high surface area to the weight ratio of CNTs,3 the last factor is considered to be especially important3,11,12,18 when CNPCs are considered. Developing a fundamental understanding of failure mechanisms at the interface is essential to making better materials suited to a broader range of applications.

Experimental investigation of CNPC interfacial properties with Raman spectroscopy13,15,19 and CNT pull-out tests using atomic-force microscopy20–22 can be instrumental when evaluating the interfacial strength. However, measuring the properties of a single CNT-polymer interface, let alone a single attachment can be challenging due to the nanometre diameter of the nanotubes. As the mechanism of interfacial load transfer is believed to depend on the local atomic structure,13,23–25 the molecular resolution offered by computational studies can provide complementary insight into the underlying processes. However, studying this problem using computational simulations is challenging because it is inherently multi-scale: failure is a macroscopic phenomenon, yet it originates from local changes in chemical bonding that occur at the nanoscale. The bond-breaking processes at the polymer-CNT attachment point that initiate failure are quantum-mechanical in nature, yet the mechanisms by which stresses are transferred through the disordered polymer occur on length-scales far in excess of anything that can be simulated quantum-mechanically. Additionally, strain rates in this class of problem can span between 10−4 and 104 s−18,26 and even the highest ones can only be considered in a limited capacity due to time scale constraints of atomistic simulations.

In CNPCs, effects that would benefit from a quantum mechanics (QM) description, such as changes in bonding, are present only at points where interfacial failure is initiated while the behavior of the surrounding environment can be well-described by a less accurate, classical model. As such, using a single method is not a necessity and the multi-scale nature of the problem can be addressed using a hybrid approach that combines an accurate QM Hamiltonian with an inexpensive classical forcefield. Similar simulation techniques have been successfully applied in the past: QM/ molecular mechanics (MM) methods have been used to study crack propagation in a silicon crystal, resulting in an accurate treatment of atomic-scale effects and the prediction of the stability of crack propagation through various cleavage planes27–30 in agreement with the experimental results for the first time; atomistic QM/MM methods have also been used to study CNTs with the ONIOM approach,31 successfully predicting elastic properties of defective CNTs32 and interactions of CNTs with non-covalently bonded molecules.33,34 Investigation of CNT-polymer interfaces, however, has been mostly limited to classical molecular dynamics (MD)25,35–37 and atomistic-continuum hybrid approaches38–40 studying the properties of CNPC structures of 100 nm and larger. In this work, we investigate interfacial interactions on a smaller length scale, studying a simplified model of the CNT-polymer interface consisting of the fundamental building-block of the interface, namely, the attachment of a single CNT and polymer chain. In our method, we supplement the fully atomistic, classical approach with a quantum mechanical description of chemical bonding.

The complexity of polymer-based composites makes application of hybrid simulation techniques a challenging task as it is often unclear which elements of the system should be treated with which level of theory.41 Here, we demonstrate a QM/MM hybrid simulation technique with an automated, dynamic selection of QM regions based on atom-resolved potential energies. We show that the approach presented can be successfully applied to study the critical strength of a CNT-polymer attachment, representing a simplified model for the composite interface, and reproduces fully QM results at modest computational cost. The method is used to evaluate a variety of strategies for grafting polymer chains to the CNT surface and enables us to propose candidates chemistries that could be used to optimise the interfacial properties of high-performance CNPCs.

The remainder of the paper is structured as follows: in Sec. II, we describe our novel hybrid simulation method; in Sec. III, we show the results of a dynamical study of the mechanism of load transfer at the CNT-polymer interface; in Sec. IV, we demonstrate the locality of quantum mechanical effects in CNT/polymer systems; Sec. V describes the QM/MM study of fracture in defective CNT, while Sec. VI contains the conclusions of this paper.

The method demonstrated here is based on the “Learn on the fly” approach29,42 and aims to compute a purely local quantity—the force on each atom.43 Forces are calculated using the quantum method in the chosen region and a classical forcefield in the rest of the system instead of deriving all of them from a combined Hamiltonian. The procedure involving selecting the QM regions and using force mixing to find atomic forces is described in detail below and shown schematically in Fig. 1.

FIG. 1.

Schematic representation of QM/MM hybrid time step.

FIG. 1.

Schematic representation of QM/MM hybrid time step.

Close modal

Classical forcefields are designed to describe atomic interactions in the vicinity of an equilibrium configuration. However, when the structure is significantly disturbed (as would be the case during a fracture event), the description can be much less accurate. Within a classical picture, the total potential energy of an atom exceeding a certain threshold indicates the structure locally moving out of its stable configuration, at which point the classical treatment is no longer accurate and a quantum mechanical description is desired. In our scheme, atoms fulfilling this criterion are flagged and individual inner QM clusters are constructed by including all particles within a selected number of bond hops from each particle. Overlapping clusters are then joined to avoid unnecessary calculations although they could be treated separately to facilitate better parallelization as discussed by Csányi et al.43 and more recently by Caccin et al.44 It is important to note that the flagging mechanism is hysteretic: there is a lower energy threshold for an atom to remain within a QM region than to enter it. Selecting only out-of-equilibrium atoms has an additional benefit of ensuring that the expensive quantum mechanical calculations are only performed when needed. In most simulations, the majority of time steps do not involve any QM computations which allows for very efficient calculations. More detailed discussion of the flagging mechanism can be found in the supplementary material.

Once the inner QM clusters are found, the next step involves including the buffer region which is constructed by expanding the inner cluster through adding atoms within a fixed number of bond-hops. The full cluster constructed from inner and buffer regions is carved out of the structure, and any dangling bonds are hydrogenated to form a structure used to compute quantum mechanical properties. It is important to consider the size and the geometry of the buffer region used to ensure accurate force evaluation; it is advised to find the right balance between computational speed and accuracy. This discussion is picked up in Sec. IV.

After QM calculation, the quantum mechanical forces are used to replace the classical forces on the atoms in the inner region while QM forces on atoms in the buffer region are discarded. The boundary between inner and buffer regions is abrupt; however, ensuring smooth force mixing is not as significant as excluding low accuracy forces evaluated on the buffer region atoms; this has been discussed in more detail in Refs. 41 and 45. Discontinuities which may affect the action-reaction principle are treated by carefully re-fitting the fundamental forcefield to match the QM description around equilibrium. This approach to treating the boundary between two levels of theory is validated in Sec. V.

A typical simulation using the QM/MM routine consists of a multitude of iterations following the above description. It is summarized below and shown schematically in Fig. 1.

  1. Classical evaluation. Forces and energies on all atoms are calculated using a classical forcefield.

  2. QM flagging. Atoms that move out of an equilibrium configuration, which indicates uncommon interactions, are flagged based on their potential energies.

  3. Cluster carving. The simulation cell is created by carving out the selected area from the original geometry, placing it in a vacuum and terminating all the dangling bonds with hydrogen atoms.

  4. Cluster QM calculation. Prepared cluster is simulated using a quantum-mechanical method.

  5. Force mixing. Classical forces on flagged atoms in the simulation are replaced by the quantum mechanical forces.

  6. Time evolution. The system is evolved in time using the Verlet algorithm with the mixed quantum/classical forces; the end state of the system becomes the new initial state.

Within the following work, the Consistent Valence Forcefield (CVFF)46 as implemented in Large-scale Atomic/ Molecular Massively Parallel Simulator (LAMMPS)47 was used as a base classical potential. Two and three body forcefield terms for the polypropylene backbone and CNT carbon atoms were rescaled using the procedure described in Sec. V applied on the 10-monomer polypropylene chain and a 40 Å, end-capped (10, 0) CNT, respectively. In this study, density-functional tight binding (DFTB)48 implemented in the PLATO49 simulation package was used as the quantum mechanical method but could readily be extended to the large number of quantum mechanical codes interfaced to the Atomic Simulation Environment (ASE) framework.50 The implementation of the QM/MM coupling described here can be found in the Multi Cluster Force Mixing (MCFM) module in the MatSciPy library.51 DFTB calculations were performed self consistently with the parameters by Krishnapriyan et al.,52 using an electronic temperature of 0.015 Ry and Pulay mixing. No electrostatic embedding was used following the cluster embedding approach from Peguiron et al.;45 more details can be found in Sec. IV. It is important to note that in principle any combination of two atomistic simulation techniques can be coupled using this scheme as long as access to per-atom potential energy is provided by the method being used as an inexpensive force evaluator (MM element). The reactive forcefield simulations used for comparison were performed using the ReaxFF forcefield53 with parameterization from54 extending the CHO forcefield55 with nitrogen atom interactions. Throughout the simulation, an NVT ensemble was enforced using the Langevin thermostat as implemented in ASE.

The QM/MM technique was used to simulate a (10, 0) CNT with a polypropylene (PP) strand attached to the surface of the CNT as a representation of a CNT-polymeric interface. The attached polymer was composed of 25 PP monomers which accounts to more than two times the persistence length of the polymer which ensures an appropriate representation of an actual chain. A small velocity of 10−3 Å/fs in the outward radial direction with respect to CNT was applied to one of the chain atoms to simulate an external load leading to critical failure. Throughout the simulation, the force on the pulled out atom was recorded to measure the system’s response to applied strain. The simulation setup is shown in Fig. 2.

FIG. 2.

A system representing CNT-polypropylene composite composed of a (10, 0) CNT joined with a short PP strand using a bromide salt residue and a benzene ring. The gray atoms represent carbon, white hydrogen, red oxygen, and blue nitrogen. The arrow indicates a fixed velocity used to simulate external load. Additionally, the critical failure points for simulations with QM/MM approach and ReaxFF forcefield are shown.

FIG. 2.

A system representing CNT-polypropylene composite composed of a (10, 0) CNT joined with a short PP strand using a bromide salt residue and a benzene ring. The gray atoms represent carbon, white hydrogen, red oxygen, and blue nitrogen. The arrow indicates a fixed velocity used to simulate external load. Additionally, the critical failure points for simulations with QM/MM approach and ReaxFF forcefield are shown.

Close modal

The study was designed to compare CNT-polymer grafting mechanisms by investigating CNT-PP systems as the one described above, with various molecules used to join the CNT surface and polymer chain. The attachments shown in Fig. 3 are composed of the bromide salt residue with a benzene ring in panels (a) and (d), single carbon atom with a benzene ring in panels (b) and (e), as well as single nitrogen atom with a benzene ring in panels (c) and (f). Each of the attachment types is used to compose two systems: one where a single bond connects the molecule with the CNT surface [(a)–(c)] and one when two single bonds are used [(d)–(f)]. For each configuration, the force necessary to break the interfacial connection was found by analyzing the atomic force on the pull-out atom as a function of simulation time which has a linear relation with its displacement. The curves were initially smoothed by averaging over a 30 fs window to account for carbon-carbon bond vibrations. The critical force was measured by finding the maximum recorded value and taking an average over a 100 fs window to account for chain oscillations. For each grafting mechanism, the process was repeated for a set of six simulations at different temperatures T = 50, 100, 250, 200, and 300 K to get the mean and standard deviation.

FIG. 3.

CNT attachment point with different molecules used as link between the benzene ring and CNT surface.

FIG. 3.

CNT attachment point with different molecules used as link between the benzene ring and CNT surface.

Close modal

The results of the study using the QM/MM approach were compared with the outcomes of similar simulations with the ReaxFF forcefield.53 Both methods were evaluated against a fully QM simulation of a similar, slightly smaller system composed of a 20 Å CNT, 10 PP monomers, and the same connecting molecules, performed at T = 100 K; the system size was reduced for this study for computational feasibility. The critical force for systems with different attachments computed using the described methods is shown in Fig. 4. The results from the QM/MM simulations are in agreement with the fully QM results which indicates that the hybrid approach can effectively reproduce the results of the quantum mechanical model it is based on.

FIG. 4.

Critical force on the pulled-out particle for simulations of six similar systems using both QM/MM approach and ReaxFF for force evaluation. Each structure was composed of a (10, 0) CNT, an attachment molecule, a benzene ring, and a polypropylene strand in the form of CNT-attachment molecule-benzene-PP; different attachment point chemistries are shown in Fig. 3. Each point is located at coordinates dictated by the breaking force from a full QM simulation on the x-axis and breaking force from an approximated method on the y-axis. Each colour represents a single attachment chemistry, and shapes represent the simulation method, squares for ReaxFF and diamonds for QM/MM.

FIG. 4.

Critical force on the pulled-out particle for simulations of six similar systems using both QM/MM approach and ReaxFF for force evaluation. Each structure was composed of a (10, 0) CNT, an attachment molecule, a benzene ring, and a polypropylene strand in the form of CNT-attachment molecule-benzene-PP; different attachment point chemistries are shown in Fig. 3. Each point is located at coordinates dictated by the breaking force from a full QM simulation on the x-axis and breaking force from an approximated method on the y-axis. Each colour represents a single attachment chemistry, and shapes represent the simulation method, squares for ReaxFF and diamonds for QM/MM.

Close modal

We can use the QM/MM results to investigate the properties of interfaces composed of different attachments. In general, systems in which the doubly bonded attachments were used perform much better than their counterparts in all cases. The systems are otherwise identical, which indicates that altering the chemistry at the attachment point to increase the binding energy can lead to an increase in the CNT–polymer interface strength by a factor of two. In fact, in the doubly bonded systems, the critical failure is caused by the polymer chain breaking instead of detaching from the CNT which shows that carefully selecting the grafting strategy can lead to maximising the interfacial shear strength (ISS) to the limit dictated by the strength of the polymer used in the composite.

The key insights from those simulations are not captured when ReaxFF is used to drive the dynamics. The breaking force needed to detach the polymer chain from the CNT surface is systematically overestimated, and the distinction between high and low strength systems is not accurate. Attachment b is reported to break at relatively low load; however, all of the other systems are considered to be of similar, high strength. In general, ReaxFF with parameters from Ref. 54 offers a correct qualitative description of systems where the key elements are carbon-based; in cases (d), (e), and (f), the rupture occurs inside a carbon-based polymer chain, while in system b, only carbon atoms are present, and the ReaxFF simulations correctly predict relative strengths of those systems. On the other hand, in cases (a) and (c), where the attachment includes various atoms different than carbon, the electronic structure of the functional group is more intricate and a parameterized description of bonding is not accurate. As a result, ReaxFF predicts that for systems (a) and (c), the critical failure occurs at a much higher pull-out strength and is caused by a rupture inside the polypropylene chain, while the DFTB based hybrid approach shows that it is due to the polymer detaching from the CNT. This is shown in detail in Fig. 2.

While it should be possible to reproduce the QM results exactly with ReaxFF forcefield by carefully re-fitting the parameters, the results would need to be validated against fully QM simulations for each system in the study. This approach is time consuming and highly limiting in terms of exploring properties of unknown structures. The method demonstrated here is designed to be more transferable; the approach relies on the results of ab initio or semi-empirical computations and as such could be used to make predictions about previously unseen compositions and screen for high performance candidates as presented here. Additionally, the approach could be applied to complement the current state of the art studies of CNT-based structures.56,57

The methodology has not been used to study this class of materials before, and as such, a significant effort must be directed at validating and fine-tuning the approach. In this section, the locality of quantum mechanical calculations, as well as the selection of different QM clusters, will be investigated.

As mentioned in Sec. II, a hybrid scheme has areas described by two levels of theory and for this approach to be valid, the interactions described by the higher level theory must be localized to a small neighbourhood. Properties of atoms within the quantum region must be computed with high accuracy by only considering a finite neighbourhood, and if the region is large enough, it should be unaffected by the classical treatment of the remainder of the system. Both of those requirements are brought to completion when the strong locality condition is fulfilled. It can be defined as43 

nxinEtotxj0  as  xixjn,ij,
(1)

where Etol is the total energy and xi and xj are the positions of atoms i and j, respectively. The non-metallic systems obey strong locality, and the intrinsically long-range effects such as Coulombic or van der Waals interactions can be accurately treated with classical approximations outside the hybrid framework.43 

Here, the force locality has been tested for ten different configurations of a functionalized CNT structure shown in Fig. 5 as a representation of a CNT-polymer composite. For each structure, the forces acting on individual atoms were calculated for the equilibrated configuration, the central C1 carbon atom (see Fig. 5 for notation) was displaced by 0.01 Å, and the atomic forces were recalculated for the new structure. The relative change in forces acting on individual atoms as a function of their distance from the central C1 atom is given in Fig. 6. The change in forces decays exponentially with increasing distance until an accuracy comparable with simulation noise is achieved. This result indicates that the strong locality assumption is indeed correct for the functionalized CNT systems and CNPC structures. Strong force locality demonstrated here justifies the choice of buffered force QM/MM approach without the use of electrostatic embedding as using a buffer region with sufficient width allows forces to converge to reference values obtained from fully quantum mechanical calculations; this discussion is picked up in the supplementary material. The buffer region of 6 neighbour hops was chosen for practical calculation based on convergence tests presented in the supplementary material.

FIG. 5.

A system representing a CNT-polyethylene composite composed of a (6, 6) CNT functionalized with a short PE strand.

FIG. 5.

A system representing a CNT-polyethylene composite composed of a (6, 6) CNT functionalized with a short PE strand.

Close modal
FIG. 6.

A plot of the relative changes in atomic forces as a result of displacing the C1 carbon atom, given as a function of atom’s distance from the central C1 carbon. The system and the notation are shown in Fig. 5.

FIG. 6.

A plot of the relative changes in atomic forces as a result of displacing the C1 carbon atom, given as a function of atom’s distance from the central C1 carbon. The system and the notation are shown in Fig. 5.

Close modal

The hybrid method has been further validated by investigating the elastic response of a defective CNT. A 40 Å, end-capped (10, 0) CNT with a single atom vacancy was strained quasi-statically by iteratively straining the structure in 1% increments and optimising the geometry while constraining the positions of end-atoms. The structure is shown in Fig. 7. Throughout the simulation, the energy-based criterion was used to flag the quantum mechanical regions. In the beginning, all atoms were treated classically, and throughout the run, the QM inner region and the buffer were selected automatically at 7% strain, accurately flagging the defective part of the CNT once the system moved away from its equilibrium structure; this is shown in Fig. 8.

FIG. 7.

Figure representing the CNT system in equilibrium. Green region represents inner QM atoms, red the buffer region, while gray and white represent the fully classical carbons and hydrogens, respectively. Blue markings are to indicate the neighbours of the vacancy in the CNT.

FIG. 7.

Figure representing the CNT system in equilibrium. Green region represents inner QM atoms, red the buffer region, while gray and white represent the fully classical carbons and hydrogens, respectively. Blue markings are to indicate the neighbours of the vacancy in the CNT.

Close modal
FIG. 8.

A plot of average force on constrained, end atoms as a response to applied strain in the top panel and the number of QM atoms as a function of strain in the bottom panel.

FIG. 8.

A plot of average force on constrained, end atoms as a response to applied strain in the top panel and the number of QM atoms as a function of strain in the bottom panel.

Close modal

The study was performed three times using the CVFF, DFTB, and QM/MM approaches presented here. The two and three body terms in the forcefield CVFF were re-fitted to match the DFTB equilibrium length of the CNT as well as its elastic response in the range of 0%–3% strain. Even though the force-mixing approach is naive, a simple forcefield re-fitting ensures that stresses are accurately transferred through the classical-quantum barrier, demonstrated by the QM/MM scheme accurately reproducing the DFTB critical stress, as shown in Fig. 8.

The results presented in this work demonstrate that an accurate description of QM effects such as bond breaking is essential when simulating interfacial critical failure in CNT-polymer composites. Such effects govern the initial evolution of the failure process, and an accurate description of the surrounding electronic structure is required to simulate them correctly.

The need to include QM description in the large-scale MD simulation necessary to describe the load transfer through the polymer chains creates a difficult multiscale challenge. Here, it has been successfully addressed by employing a quantum/classical hybrid simulation technique with automated, dynamic flagging offering QM accuracy only in the crucial parts of the system, thus maintaining modest computational cost. The method was shown to produce similar results to a fully QM technique, and it is based on problems of critical failure in CNT and polymer-based systems.

The hybrid method demonstrated here is best suited to investigate problems where the crucial effects dictating the behavior of the system are local in nature; most critical failure effects, especially in non-homogenous systems, fall into that category. In this class of problems, the QM/MM approach could be used to take advantage of accuracy and predictive power of methods such as DFT at modest computational cost, allowing for studies at a large length and time scales. As a result, the method could be used as an improvement over classical forcefields in problems where transferability is important such as an investigation of various attachment strategies in the search for high-performing CNT-polymer interfaces.

Here, the QM/MM approach was used to study different CNT-polymer interfaces to identify the frailest element of the CNT-polymer interface and test potential high-performing replacements. The findings indicate that the chemistry at the attachment point significantly contributes to the effective ISS of the interface, and using a carefully chosen attachment can lead to increasing the ISS to the limit dictated by the strength of the polymer used in the composite. Computational investigations such as this one could be used as guidance for future developments of CNPCs with covalently bonded interfaces.

Please see supplementary material for a description of, and the results associated with, convergence test of the QM/MM atomic forces with respect to the size of the buffer region, details about the cluster QM calculations, QM flagging energy criterion as well as a discussion of force vs displacement curves and QM/MM force evaluation errors for simulations presented in Sec. III.

The authors would like to acknowledge Andrew Horsfield for his help with DFTB simulations.

This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie Grant Agreement No. 642890 through TheLink research and training project.

This work was supported through a studentship in the Centre for Doctoral Training on Theory and Simulation of Materials at Imperial College London funded by the EPSRC (Grant No. EP/L015579/1)

Calculations were performed on the Imperial College London High-Performance Computing Service.

We acknowledge the Thomas Young Centre under grant number TYC-101.

1.
J.
Njuguna
and
K.
Pielichowski
,
Adv. Eng. Mater.
5
,
769
(
2003
).
2.
Y.
Liu
and
S.
Kumar
,
ACS Appl. Mater. Interfaces
6
,
6069
(
2014
).
3.
J. N.
Coleman
,
U.
Khan
,
W. J.
Blau
, and
Y. K.
Gun’ko
,
Carbon
44
,
1624
(
2006
).
4.
M.
Griebel
and
J.
Hamaekers
,
Comput. Methods Appl. Mech. Eng.
193
,
1773
(
2004
).
5.
M.
Biercuk
,
S.
Ilani
,
C.
Marcus
, and
P. L.
McEuen
,
Top. Appl. Phys.
111
,
455
(
2008
).
6.
M. M.
Zamani
,
A.
Fereidoon
, and
A.
Sabet
,
Iran. Polym. J.
21
,
887
(
2012
).
7.
M.
Tehrani
,
A. Y.
Boroujeni
,
T. B.
Hartman
,
T. P.
Haugh
,
S. W.
Case
, and
M. S.
Al-Haik
,
Compos. Sci. Technol.
75
,
42
(
2013
).
8.
S.
Laurenzi
,
R.
Pastore
,
G.
Giannini
, and
M.
Marchetti
,
Compos. Struct.
99
,
62
(
2013
).
9.
M. L.
Manchado
,
L.
Valentini
,
J.
Biagiotti
, and
J.
Kenny
,
Carbon
43
,
1499
(
2005
).
10.
X. L.
Xie
,
Y. W.
Mai
, and
X. P.
Zhou
,
Mater. Sci. Eng. R: Rep.
49
,
89
(
2005
).
11.
P. C.
Ma
,
N. A.
Siddiqui
,
G.
Marom
, and
J. K.
Kim
,
Composites, Part A
41
,
1345
(
2010
).
12.
B.
Arash
,
Q.
Wang
, and
V. K.
Varadan
,
Sci. Rep.
4
,
6479
(
2014
).
13.
P. C.
Ma
,
Q. B.
Zheng
,
E.
Mader
, and
J. K.
Kim
,
Polymer
53
,
6081
(
2012
).
14.
M. R.
Loos
,
S. H.
Pezzin
,
S. C.
Amico
,
C. P.
Bergmann
, and
L. A. F.
Coelho
,
J. Mater. Sci.
43
,
6064
(
2008
).
15.
F.
Buffa
,
G. A.
Abraham
,
B. P.
Grady
, and
D.
Resasco
,
J. Polym. Sci., Part B: Polym. Phys.
45
,
490
(
2007
).
16.
T.
Chang
,
L.
Jensen
,
A.
Kisliuk
,
R.
Pipes
,
R.
Pyrz
, and
A.
Sokolov
,
Polymer
46
,
439
(
2005
).
17.
D.
Qian
,
E. C.
Dickey
,
R.
Andrews
, and
T.
Rantell
,
Appl. Phys. Lett.
76
,
2868
(
2000
).
18.
K.
Liao
and
S.
Li
,
Appl. Phys. Lett.
79
,
4225
(
2001
).
19.
D.
Roy
,
S.
Bhattacharyya
,
A.
Rachamim
,
A.
Plati
, and
M.-L. L.
Saboungi
,
J. Appl. Phys.
107
,
043501
(
2010
).
20.
T.
Tsuda
,
T.
Ogasawara
,
F.
Deng
, and
N.
Takeda
,
Compos. Sci. Technol.
71
,
1295
(
2011
).
21.
A. H.
Barber
,
S. R.
Cohen
, and
H. D.
Wagner
,
Appl. Phys. Lett.
82
,
4140
(
2003
).
22.
C. A.
Cooper
,
S. R.
Cohen
,
A. H.
Barber
, and
H. D.
Wagner
,
Appl. Phys. Lett.
81
,
3873
(
2002
).
23.
Q.
Xiong
and
S.
Meguid
,
Eur. Polym. J.
69
,
1
(
2015
).
24.
J.
Wernik
,
B.
Cornwell-Mott
, and
S.
Meguid
,
Int. J. Solids Struct.
49
,
1852
(
2012
).
25.
Y.
Li
and
G. D.
Seidel
,
Comput. Mater. Sci.
107
,
216
(
2015
).
26.
P.
Jindal
,
S.
Pande
,
P.
Sharma
,
V.
Mangla
,
A.
Chaudhury
,
D.
Patel
,
B. P.
Singh
,
R. B.
Mathur
, and
M.
Goyal
,
Composites, Part B
45
,
417
(
2013
).
27.
N.
Bernstein
and
D. W.
Hess
,
Phys. Rev. Lett.
91
,
025501
(
2003
).
28.
M. J.
Buehler
,
H.
Tang
,
A. C. T.
van Duin
, and
W. A.
Goddard
 III
, “
Threshold crack speed controls dynamical fracture of silicon single crystals
,”
Phys. Rev. Lett.
99
,
165502
(
2007
).
29.
G.
Csányi
,
T.
Albaret
,
M. C.
Payne
, and
A.
De Vita
,
Phys. Rev. Lett.
93
,
175503
(
2004
).
30.
J. R.
Kermode
,
T.
Albaret
,
D.
Sherman
,
N.
Bernstein
,
P.
Gumbsch
,
M. C.
Payne
,
G.
Csányi
, and
A.
De Vita
,
Nature
455
,
1224
(
2008
).
31.
L. W.
Chung
,
W. M.
Sameera
,
R.
Ramozzi
,
A. J.
Page
,
M.
Hatanaka
,
G. P.
Petrova
,
T. V.
Harris
,
X.
Li
,
Z.
Ke
,
F.
Liu
,
H. B.
Li
,
L.
Ding
, and
K.
Morokuma
, “
The ONIOM method and its applications
,”
Chem. Rev.
115
,
5678
(
2015
).
32.
R.
Khare
,
S. L.
Mielke
,
J. T.
Paci
,
S.
Zhang
,
R.
Ballarini
,
G. C.
Schatz
, and
T.
Belytschko
,
Phys. Rev. B
75
,
075412
(
2007
).
33.
P. N.
Samanta
and
K. K.
Das
,
Int. J. Quantum Chem.
116
,
411
(
2016
).
34.
K.
Hu
,
Q.
Zhuang
,
X.
Liu
, and
Z.
Han
,
Polym. Compos.
36
,
1454
(
2015
).
35.
S.
Namilae
,
U.
Chandra
,
A.
Srinivasan
, and
N.
Chandra
,
Comput. Model. Eng. Sci.
22
,
189
(
2007
).
36.
S. J. V.
Frankland
,
A.
Caglar
,
D. W.
Brenner
, and
M.
Griebel
,
J. Phys. Chem. B
106
,
3046
(
2002
).
37.
S.
Chowdhury
and
T.
Okabe
,
Composites, Part A
38
,
747
(
2007
).
38.
F.
Han
,
Y.
Azdoud
, and
G.
Lubineau
,
Comput. Mater. Sci.
81
,
641
(
2014
).
39.
V.
Unnikrishnan
and
J. N.
Reddy
,
Finite Elements Anal. Des.
49
,
13
(
2012
).
40.
C.
Li
and
T. W.
Chou
,
Compos. Sci. Technol.
66
,
2409
(
2006
).
41.
N.
Bernstein
,
J. R.
Kermode
, and
G.
Csányi
,
Rep. Prog. Phys.
72
,
026501
(
2009
).
42.
A.
De Vita
and
R.
Car
,
Mater. Res. Soc. Symp. Proc.
491
,
473
(
1998
).
43.
G.
Csányi
,
T.
Albaret
,
G.
Moras
,
M. C.
Payne
, and
A.
De Vita
,
J. Phys.: Condens. Matter
17
,
R691
(
2005
).
44.
M.
Caccin
,
Z.
Li
,
J. R.
Kermode
, and
A.
De Vita
,
Int. J. Quantum Chem.
115
,
1129
(
2015
).
45.
A.
Peguiron
,
L.
Colombi Ciacchi
,
A.
De Vita
,
J. R.
Kermode
, and
G.
Moras
,
J. Chem. Phys.
142
,
064116
(
2015
).
46.
P.
Dauber-Osguthorpe
,
V. A.
Roberts
,
D. J.
Osguthorpe
,
J.
Wolff
,
M.
Genest
, and
A. T.
Hagler
,
Proteins: Struct., Funct., Bioinf.
4
,
31
(
1988
).
47.
S.
Plimpton
,
J. Comput. Phys.
117
,
1
(
1995
); e-print arXiv:nag.2347 [10.1002].
48.
D.
Porezag
,
T.
Frauenheim
,
T.
Köhler
,
G.
Seifert
, and
R.
Kaschner
,
Phys. Rev. B
51
,
12947
(
1995
); e-print arXiv:1011.1669v3.
49.
S. D.
Kenny
and
A. P.
Horsfield
,
Comput. Phys. Commun.
180
,
2616
(
2009
).
50.
A. H.
Larsen
,
J. J.
Mortensen
,
J.
Blomqvist
,
I. E.
Castelli
,
R.
Christensen
,
M.
Dułak
,
J.
Friis
,
M. N.
Groves
,
B.
Hammer
,
C.
Hargus
,
E. D.
Hermes
,
P. C.
Jennings
,
P. B.
Jensen
,
J.
Kermode
,
J. R.
Kitchin
,
E. L.
Kolsbjerg
,
J.
Kubal
,
K.
Kaasbjerg
,
S.
Lysgaard
,
J. B.
Maronsson
,
T.
Maxson
,
T.
Olsen
,
L.
Pastewka
,
A.
Peterson
,
C.
Rostgaard
,
J.
Schiøtz
,
O.
Schütt
,
M.
Strange
,
K. S.
Thygesen
,
T.
Vegge
,
L.
Vilhelmsen
,
M.
Walter
,
Z.
Zeng
, and
K. W.
Jacobsen
,
J. Phys. Condens. Matter
29
,
273002
(
2017
).
51.
L.
Pastewka
and
J.
Kermode
, MatSciPy library: https://github.com/libAtoms/matscipy.
52.
A.
Krishnapriyan
,
P.
Yang
,
A. M. N.
Niklasson
, and
M. J.
Cawkwell
,
J. Chem. Theory Comput.
13
,
6191
(
2017
).
53.
A. C. T.
Van Duin
,
S.
Dasgupta
,
F.
Lorant
, and
W. A.
Goddard
,
J. Phys. Chem. A
105
,
9396
(
2001
).
54.
A. M.
Kamat
,
A. C. T.
Van Duin
, and
A.
Yakovlev
,
J. Phys. Chem. A
114
,
12561
(
2010
).
55.
K.
Chenoweth
,
A. C. T.
van Duin
, and
W. A.
Goddard
,
J. Phys. Chem. A
112
,
1040
(
2008
); e-print arXiv:0311084 [physics].
56.
E. M.
Byrne
,
M. A.
McCarthy
,
Z.
Xia
, and
W. A.
Curtin
,
Phys. Rev. Lett.
103
,
045502
(
2009
).
57.
F.
Pavia
and
W. A.
Curtin
,
Acta Mater.
59
,
6700
(
2011
).

Supplementary Material